Subroutine evolve( tterm, phi_ir, ndim_ir, whichir,nchir)
Use Numeric_Kinds_Module
Use CommonInfo_Module
Use Time_Module
Use User_Module
Use APH_Module
Use Jacobi_Module
Use SymGroup_Module
Use SurfaceAph_Module
Use NIP_Module
Use Phi_Module
Implicit None
save
!=========================================================================================
! Edited by: Jeff Crawford
!
! This routine uses the Chebychev method to evolve the
!      time-dependent wavepacket.
!
! See C. Leforestier, et al J. Comp. Phys. 94, 59 (1991) for 
!      explanation of method
!
! Variables:
!     bessel_arg: argument of bessel functions
!     bessel_func: Array containing bessel function values
!     ncheby: number of terms in Chebychev expansion
!    

Integer,Intent(In) :: tterm, whichir, ndim_ir,nchir

LOGICAL, PARAMETER :: debug=.true.

Character(len=100) :: filename
Character(len=1)   :: ir1, label1
Complex(dp) :: phi_ir(ndim_ir)
Complex(dp),Allocatable :: phi(:)

Integer ir_pt, ispt, irho, itheta, ichi, indexchi

Integer i,aphindex

!=====
!TEMP
integer  :: chishift
real(dp) :: x, y, probamp, norm, chi_sum, theta_sum, rho_sum

!=========================================================================================
! Extract relevant irreducible representation part of wave function
! Do on first time step for each iredducible representation

IF (tterm.eq.1) THEN
 
  ALLOCATE(phi(ndim))
!==============================================================================
! Read in symetrized wave function
  filename = 'iwf_sym.bin'
  OPEN(UNIT=bin_unit0,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')
  READ(bin_unit0) phi
  CLOSE(bin_unit0,STATUS='keep')

  !nchir=nchi_ir(whichir)
  indexchi=indexchi_array(whichir)
  chishift=chishift_ir(whichir)
  ir1=label1(whichir)

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

        ispt=aphindex(ntheta,nchi,irho,itheta,ichi+indexchi,.false.)
        ir_pt=aphindex(ntheta,nchir,irho,itheta,ichi,.false.)

        phi_ir(ir_pt) = phi(ispt)

      ENDDO
    ENDDO
  ENDDO

  DEALLOCATE(phi)

!=============================================================================
! Get ir nip

  ALLOCATE(vnip_ir(ndim_ir))
   DO irho=1,nrho
     DO itheta=1,ntheta
       DO ichi=1,nchir

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

         vnip_ir(ir_pt) = Exp(-vnip(irho)*deltat)

       ENDDO
     ENDDO
   ENDDO

!==============================================================================
! Check normalization

  rho_sum=0.d0
  DO irho=1,nrho
   theta_sum=0.d0
    DO itheta=1,ntheta
      chi_sum=0.d0
      DO ichi=1,nchir
        ir_pt=aphindex(ntheta,nchir,irho,itheta,ichi,.false.)
        chi_sum=chi_sum+(Conjg(phi_ir(ir_pt))*phi_ir(ir_pt))*deltachi
      ENDDO
      theta_sum=theta_sum+chi_sum!*weights(itheta)
    ENDDO
    rho_sum=rho_sum+theta_sum*deltarho
  ENDDO
  PRINT*,'---------------------------------------------------------------------'
  PRINT*,'Initial wave packet normalization = ', rho_sum

!==============================================================================
! Allocate arrays for propagation

  !ALLOCATE( phi_real(ndim_ir), phi_imag(ndim_ir) )
  !ALLOCATE( phi_n(ndim_ir), phi_nm1(ndim_ir), phi_nm2(ndim_ir) )
  ALLOCATE( phi_real_sum(ndim_ir), phi_imag_sum(ndim_ir))
  ALLOCATE( phi_real_n(ndim_ir), phi_real_nm1(ndim_ir), phi_real_nm2(ndim_ir) )
  ALLOCATE( phi_imag_n(ndim_ir), phi_imag_nm1(ndim_ir), phi_imag_nm2(ndim_ir) )

ENDIF

!==============================================================================
!  Calculate chebychev polynomials for time evolution.
!  Propagate the wavepacket.

  !call findtime(.true.,'zz')
  !DO irho=1,10
  Call ztcheby( ndim_ir, phi_ir, whichir,nchir)
  !enddo
  !call findtime(.false.,'zz')

  !if (tterm.eq.2) print*, phi_ir
  !phi_real_n=REAL(phi_ir)
  !phi_imag_n=AIMAG(phi_ir)
  !call findtime(.true.,'zz')
  !DO irho=1,10
  !Call zztcheby( ndim_ir, whichir,nchir)
  !enddo
  !call findtime(.false.,'zz')
  !phi_ir=Cmplx(phi_real_n,phi_imag_n,dp)

!==============================================================================
! Deallocate Arrays

if (tterm.eq.tstep) then
  !DEALLOCATE( phi_real, phi_imag )
  !DEALLOCATE( phi_n, phi_nm1, phi_nm2 )
  DEALLOCATE( phi_real_sum, phi_imag_sum )
  DEALLOCATE( phi_real_n, phi_real_nm1, phi_real_nm2 )
  DEALLOCATE( phi_imag_n, phi_imag_nm1, phi_imag_nm2 )
  DEALLOCATE(vnip_ir)
endif

Return
End Subroutine evolve