SUBROUTINE h_apply(phi,nlength,whichir)
USE Numeric_Kinds_Module
USE FileUnits_Module
USE APH_Module
USE Aph_Ham_Module
USE Cheby_Module
USE SymGroup_Module 
USE SurfaceAPH_Module
USE CommonInfo_Module
USE Convrsns_Module
USE masses_module
USE QuantumNumber_Module
IMPLICIT NONE
!===============================================================================
! This routine was written by J. Crawford and G. A. Parker
!
!   One useful paper is by David K. Hoffman, Mark Arnold
!   Donald J. Kouri. J. Phys. Chem 96, 6539-6545 1992.
!===============================================================================
!   I N P U T / O U T P U T

INTEGER,INTENT(IN) :: nlength, whichir
REAL (dp),INTENT(INOUT) :: phi(nlength)

!===============================================================================
!   I N T E R N A L S
INTEGER :: i,irho,itheta,ichi
INTEGER :: ifirst, checklength, check_ir, nchir, indexchi, chishift, ndim_tc
!REAL (dp) :: shift, norm

!===============================================================================
!    A L L O C A T A B L E

REAL (dp), ALLOCATABLE:: wavefunction (:)
REAL (dp), ALLOCATABLE:: t_chi_p(:,:), phi3(:,:,:), y3(:,:,:), y3store(:,:,:)
!
data ifirst/0/
!
save

!==============================================================================
! Initialize
!------------------------------------------------------------------------------
 IF (ifirst.eq.0) THEN
   checklength=0
   check_ir=0
 ENDIF
!==============================================================================
! Assign correct chi dimensioning based upon the appropriate irreducible
! representation.
!------------------------------------------------------------------------------
!define these early and pass them?
 nchir=nchi_ir(whichir)
 indexchi=indexchi_array(whichir)
 chishift=chishift_ir(whichir)
!==============================================================================
! Allocate irreducible representation dependent arrays
!------------------------------------------------------------------------------
 IF (whichir.ne.check_ir) THEN
   IF (Allocated(y3)) THEN 
     DEALLOCATE (wavefunction, t_chi_p, phi3, y3)
     DEALLOCATE (t_chi_temp, num_eig, chi_eig, chi_evec, compr, compr_mat)
     DEALLOCATE(t_chi_clip, tempvcut, tempeig)
     DEALLOCATE(y3store)      
   ENDIF
     ALLOCATE( wavefunction (nlength))
     ALLOCATE( t_chi_p(nchir,nchir),phi3(nchir,nrho,ntheta) )
     ALLOCATE( y3(nchir,nrho,ntheta), chi_evec(nchir,nchir) )
     ALLOCATE( t_chi_temp(nchir,nchir), chi_eig(nchir) )
     ALLOCATE( num_eig(nrho*ntheta), compr(nchir,nchir) )
     ALLOCATE( compr_mat(nchir,nchir,nchir), tempeig(nchir,nchir) )
     ALLOCATE( t_chi_clip(nchir,nchir), tempvcut(nchir,nchir) )
     ALLOCATE( y3store(nchir,nrho,ntheta))
 ELSE ! Check that nlength has maintained it's value
   IF (nlength.ne.checklength) THEN
     STOP'ERROR: nlength has changed before whichir'
   ENDIF
 ENDIF
!==============================================================================
! Apply the Hamiltonian to the wavefunction
!------------------------------------------------------------------------------

!==============================================================================
! Place wavefunction into a 3d array phi3

!IF(whichir.ne.check_ir) call findtime(.true.,'full')

 ndim_tc=ntheta*nchir

 DO itheta=1,ntheta
   DO irho=1,nrho
     DO ichi=1,nchir
       phi3(ichi,irho,itheta)=phi(ndim_tc*(irho-1)+nchir*(itheta-1)+ichi)
     ENDDO
   ENDDO
 ENDDO

!==============================================================================
! Apply T_rho:
!==============================================================================
 !IF(whichir.ne.check_ir) call findtime(.true.,'trho')
 CALL trho_apply(nchir,ndim_tc,phi3,y3,ifirst)
 !IF(whichir.ne.check_ir) call findtime(.false.,'trho')
 y3store=y3

!if (indexchi.ne.0) then
!print*,y3
!Stop
!endif
!==============================================================================
! Apply T_theta:
!------------------------------------------------------------------------------
 !IF(whichir.ne.check_ir)call findtime(.true.,'ttheta')
 CALL ttheta_apply(nchir,ndim_tc,phi3,y3,ifirst)
!IF(whichir.ne.check_ir)call findtime(.false.,'ttheta')
 y3store=y3store+y3
!==============================================================================
! Apply T_chi
!------------------------------------------------------------------------------
 !IF(whichir.ne.check_ir) call findtime(.true.,'tchi')
 CALL tchi_apply(nchir,ndim_tc,phi3,y3,t_chi_p,whichir,check_ir,ifirst,indexchi)
 !IF(whichir.ne.check_ir) call findtime(.false.,'tchi')

 y3store=y3store+y3

 IF (whichir.ne.check_ir) THEN
 IF (ALLOCATED(v_pot2)) DEALLOCATE(v_pot2)
 ALLOCATE(v_pot2(nchir,nrho,ntheta))
 DO irho=1,nrho
   DO itheta=1,ntheta
     DO ichi=1,nchir
       v_pot2(ichi,irho,itheta)=v_pot(irho,itheta,ichi+chishift)+((15.d0/(4.d0*usys2))*rhom2(irho))
     ENDDO
   ENDDO
 ENDDO
 ENDIF
!==============================================================================
!  Apply  V and the Chebychev "normalization"
!------------------------------------------------------------------------------
!IF(whichir.ne.check_ir)call findtime(.true.,'v')
! CALL v_apply(nchir,nlength,ndim_tc,phi,phi3,y3,wavefunction,chishift, &
!              indexchi)
!IF(whichir.ne.check_ir)call findtime(.false.,'v')
!IF(whichir.ne.check_ir)call findtime(.false.,'full')

! y3store=y3store+y3

 DO irho=1,nrho
   DO itheta=1,ntheta
     DO ichi=1,nchir
       wavefunction(ndim_tc*(irho-1)+nchir*(itheta-1)+ichi)=  &
                                         y3store(ichi,irho,itheta)
     ENDDO
   ENDDO
 ENDDO

 CALL v_apply(nchir,nlength,ndim_tc,phi,phi3,y3,wavefunction,chishift, &
              indexchi)


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

!==============================================================================
! Assign the new wave function to the output array
!------------------------------------------------------------------------------
 phi=wavefunction
 check_ir=whichir
 checklength=nlength
!==============================================================================
! Set ifirst = 1 to let the routine know it has been called at least once
!------------------------------------------------------------------------------

 ifirst=1

RETURN
END SUBROUTINE h_apply