!---------------------------------------------------------------------------
SUBROUTINE v_apply(nchir,nlength,ndim_tc,phi,phi3,y3,wavefunction,chishift,&
                   indexchi)
USE Numeric_Kinds_Module
USE CommonInfo_Module
USE APH_Module
USE Aph_Ham_Module
USE Cheby_Module
USE masses_module
IMPLICIT NONE

!==============================================================================
! I N P U T / O U T P U T
INTEGER,INTENT(IN) :: nchir, nlength, ndim_tc, chishift,indexchi

REAL(dp),INTENT(IN)    :: phi(nlength)
REAL(dp),INTENT(INOUT) :: wavefunction(nlength)
REAL(dp),INTENT(INOUT) :: phi3(nchir,nrho,ntheta), y3(nchir,nrho,ntheta)

!==============================================================================
! I N T E R N A L S
INTEGER :: i, irho, itheta, ichi, aphindex

!==============================================================================
!  Apply  V and the Chebychev "normalization"
!------------------------------------------------------------------------------
! DO i=1,nlength
!   irho=(i-1)/ndim_tc+1
!   itheta=(i-(irho-1)*ndim_tc-1)/nchir+1
!   ichi=i-(irho-1)*ndim_tc-(itheta-1)*nchir+chishift
!   IF (v_pot(irho,itheta,ichi).gt.vcut)THEN
!     wavefunction(i)=0.0d0
!   ELSE
!     wavefunction(i)=wavefunction(i)+phi(i)*(v_pot(irho,itheta,ichi)+ &
!                       ((15.d0/(4.d0*usys2))*rhom2 (irho)))
!     wavefunction(i)=(2.0d0*wavefunction(i)-shift*phi(i))/norm
!   ENDIF
! ENDDO

 DO irho=1,nrho
   DO itheta=1,ntheta
     DO ichi=1,nchir

        i=aphindex(ntheta,nchir,irho,itheta,ichi,.false.)

        wavefunction(i)=wavefunction(i)+phi(i)*v_pot2(ichi,irho,itheta)
        if (v_pot(irho,itheta,ichi).gt.vcut) wavefunction(i)=0.d0

     ENDDO
   ENDDO
 ENDDO

!  y3=phi3*v_pot2
 
!  wavefunction=2.0d0*wavefunction-shift*phi/norm 

! DO i=1,nlength
!    wavefunction(i)=(2.0d0*wavefunction(i)-shift*phi(i))/norm
! ENDDO

! wavefunction=(2.d0*wavefunction)-((shift/norm)*phi)

RETURN
END SUBROUTINE v_apply