SUBROUTINE delves_surface(jval, ichanl)
!
!========================================================================================
! Written by: Jeff Crawford  : adapted from OneDim codes written by G. Parker
!
! This routine obtains the asymptotic Delves' theta eigenfunctions using the two-body
! Jacobi mass-scaled s potential. These functions are defined along a constant rho 
! surface and for a fixed rotational angular momentum j.
! The eigenfunctions are solutions to equation 160 in the APH paper (JCP 87, 3888 (1987))
! Currently, this only finds the eigenfunctions for a single channel. Simple to change
!   to an arbitrary channel by making an option to pot to look at a different diatom
!
! Variables:
!           td_dim = dimension of coordinate grid
Use Numeric_Kinds_Module
USE CommonInfo_Module
USE PiFactrs_Module
USE APH_Module
USE SurfaceDelves_Module
USE QuantumNumber_Module
USE Print_Module

USE Masses_Module
IMPLICIT NONE

!=========================================================================================
!   I N P U T

INTEGER,INTENT(IN) :: jval, ichanl

!=========================================================================================
!   P A R A M E T E R S

CHARACTER(LEN=1),PARAMETER:: Jobz='V', UpLo='L', Range='V'
INTEGER, PARAMETER:: ndaf=2, mval=32, ldab=42
REAL(dp),PARAMETER:: AbsTol=1.d-12, tol=1.d-14

!=========================================================================================
!   I N T E R N A L S

CHARACTER(LEN=2) l2,  label2
CHARACTER(LEN=100) filename
INTEGER i, ipt, ivib  , kmax, otest
INTEGER neig,lwork,liwork,info,nu,istat
REAL(dp) e, coord, mag
REAL(dp) s, v, v1, v2
REAL(dp) re, v0re, v1re, v2re, rc, v0rc, v1rc, v2rc

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

INTEGER,ALLOCATABLE:: iwork(:), ifail(:)
REAL(dp),ALLOCATABLE :: psi(:,:), psisq(:), eig(:), work(:)
REAL(dp),ALLOCATABLE :: psi_temp(:,:), eig_temp(:)
REAL(dp),ALLOCATABLE :: kinetic(:,:), toeplitz(:), dafx(:,:), hermit(:)

!=========================================================================================
!   N A M E L I S T S

! Combine this with jacobi surface (essentially the same)
! Use the same wrapper as jacobi_surface to get these functions

!=========================================================================================

 l2=label2(jval)

 PRINT*,'------------------------------------------------------------------------------'
 PRINT*,'Begin Routine delves_surface for j = '//l2

!=========================================================================================
! Allocate arrays and initialize variables

 ALLOCATE( eig(td_dim), psi(td_dim,td_dim), psisq(td_dim) )
 ALLOCATE( ifail(td_dim) )
 ALLOCATE( kinetic(td_dim,td_dim) )
 ALLOCATE( toeplitz(td_dim) )
 ALLOCATE( dafx(td_dim,0:ndaf) )
 ALLOCATE( hermit(0:mval+ndaf+1) )

 liwork=MAX(1,5*td_dim)
 lwork=8*td_dim

 ALLOCATE( work(lwork), iwork(liwork) )

 psi=0.d0
 eig=0.d0
 ifail=0
 info=0
 neig=0
 mag=0.d0
 work=0.d0
 iwork=0

!=========================================================================================
! Create coordinate and diatomic potential grids
! Uses two-body potential

  IF (firstcalld) THEN

    ALLOCATE( kineticd(td_dim,td_dim)  )

    ! Delves' theta grid
    delta_td=(pi/2.d0)/(td_dim-1)
    DO ipt=1,td_dim
      td_vals(ipt) = delta_td*(ipt-1)
      ss_vals(ipt) = rho_infty*Sin(td_vals(ipt))
    ENDDO

    ss_max=ss_vals(td_dim)

    !Asymptotic value of potential
    CALL pot(ss_vals(td_dim),vasyd,v1,v2,jval,usys)

    CALL get_daf_s(delta_td,mval,dafx,hermit,ndaf,td_dim-1,toeplitz,usys2*rho_infty*rho_infty,kmax)

    CALL kin_radial_s(td_dim,toeplitz,kinetic,kmax)

    kineticd = kinetic

    firstcalld=.false.

  ELSE
    
    kinetic=kineticd

  ENDIF

! Jacobi (mass-scaled) small s grid and potential
 DO ipt=1,td_dim
   s = ss_vals(ipt)
   CALL pot(s, v, v1, v2, jval, usys, ichanl)
   td_pot(ipt)=v
   td_vp(ipt)=v1
 ENDDO

!=========================================================================================
! Output potential data

 OPEN(UNIT=output_unit0,FILE=TRIM(outdir)//'potential/delves_v_'//l2//'.txt',IOSTAT=istat,STATUS='replace',ACTION='write')
   IF (istat.ne.0) STOP 'OPEN failed - surface_delves'
 OPEN(UNIT=output_unit1,FILE=TRIM(outdir)//'potential/delves_v_s_'//l2//'.txt',IOSTAT=istat,STATUS='replace',ACTION='write')
   IF (istat.ne.0) STOP 'OPEN failed - surface_delves'
   
 DO ipt=1,td_dim
   WRITE(output_unit0,*) td_vals(ipt), td_pot(ipt)
   WRITE(output_unit1,*) ss_vals(ipt), td_pot(ipt)
 ENDDO

 CLOSE(output_unit0,IOSTAT=istat,STATUS='keep')
 IF (istat.ne.0) STOP 'CLOSE failed - surface_delves'
 CLOSE(output_unit1,IOSTAT=istat,STATUS='keep')
 IF (istat.ne.0) STOP 'CLOSE failed - surface_delves'

!=========================================================================================
! Find critical points of potential

 CALL critical(ss_vals, td_vp, td_dim, jval, vasyd, re, v0re, v1re, v2re, rc, v0rc, v1rc, v2rc, usys, ichanl)

 PRINT*,''
 PRINT*,'DELVES CRITIAL POINTS for j = '//l2
 PRINT*,''
 WRITE(*,'(1X,2(a,ES13.6,3X))') 're   = ',re,  'rc   = ',rc
 WRITE(*,'(1X,3(a,ES13.6,3X))') 'v0re = ',v0re,'v1re = ',v1re,'v2re = ',v2re
 WRITE(*,'(1X,3(a,ES13.6,3X))') 'v0rc = ',v0rc,'v1rc = ',v1rc,'v2rc = ',v2rc
 PRINT*,''

!=========================================================================================
! Add potential to kinetic energy operator to obtain full 1D Hamiltonian: 
! Store Hamiltonian in kinetic

  DO i=1,td_dim
    kinetic(i,i)=kinetic(i,i)+td_pot(i)
  ENDDO

!=========================================================================================
! Diagonalize with dsyevx

  CALL dsyevx(Jobz, Range, UpLo, td_dim, kinetic, td_dim, v0re, v0rc, 1, 1, AbsTol, neig, eig, psi, td_dim, work, lwork, iwork, ifail, info)

!=========================================================================================
! Error message check for dsyevx

  DO i=1,td_dim
    IF (ifail(i).ne.0.or.info.ne.0) THEN 
      PRINT*,'ERROR in routine dsyevx'  
      PRINT*,'ifail =',ifail(i),'  info=',info   
      STOP 'delves_surface'
    ENDIF
  ENDDO

  DEALLOCATE( ifail, iwork, work )

!=========================================================================================
! Divide by sqrt(delta_td) to provide proper amplitude (for proper normalization)

  psi=psi/Sqrt(delta_td)

!=========================================================================================
! Save number of bound states for current j value

 ndvib(jval)=neig

!=========================================================================================
! Output eigenfunctions and energies to file:

! Energies

   OPEN(UNIT=output_unit2,FILE=TRIM(outdir)//'surface_functions/delves/energy_'//l2//'.txt',IOSTAT=istat,STATUS='replace',ACTION='write')
   IF (istat.ne.0) STOP 'OPEN failed - surface_delves'
   OPEN(UNIT=output_unit3,FILE=TRIM(outdir)//'surface_functions/delves/energy_long_'//l2//'.txt',IOSTAT=istat,STATUS='replace',ACTION='write')
   IF (istat.ne.0) STOP 'OPEN failed - surface_delves'
   DO nu=1,neig
     e=eig(nu)
     WRITE(output_unit2,'(i4,e14.7)') nu, e
     WRITE(output_unit3,'(a,i4,a,3(es23.15,a))')   "nu=", nu, '  Energy=',e, "(Ha)  ", e*27.2117d0,"(eV)  ",e*219474.6314d0,"(cm-1)"
   ENDDO
   CLOSE(UNIT=output_unit2,IOSTAT=istat,STATUS='keep')
   IF (istat.ne.0) STOP 'CLOSE failed - surface_delves'
   CLOSE(UNIT=output_unit3,IOSTAT=istat,STATUS='keep')
   IF (istat.ne.0) STOP 'CLOSE failed - surface_delves'

! Wave Functions

     OPEN(UNIT=output_unit4,FILE=TRIM(outdir)//'surface_functions/delves/function_'//l2//'.txt',IOSTAT=istat,STATUS='replace',ACTION='write')
       IF (istat.ne.0) STOP 'OPEN failed - surface_delves'
     OPEN(UNIT=output_unit5,FILE=TRIM(outdir)//'surface_functions/delves/function_s_'//l2//'.txt',IOSTAT=istat,STATUS='replace',ACTION='write')
       IF (istat.ne.0) STOP 'OPEN failed - surface_delves'
     DO ipt=1,td_dim
       WRITE(output_unit4,'(1x,e14.7,1000(1X,e14.7))') td_vals(ipt), ( psi(ipt,ivib), ivib=1,neig )
       WRITE(output_unit5,'(1x,e14.7,1000(1X,e14.7))') ss_vals(ipt), ( psi(ipt,ivib), ivib=1,neig )
     ENDDO
     CLOSE(output_unit4,IOSTAT=istat,STATUS='keep')
     IF (istat.ne.0) STOP 'CLOSE failed - surface_delves'
     CLOSE(output_unit5,IOSTAT=istat,STATUS='keep')
     IF (istat.ne.0) STOP 'CLOSE failed - surface_delves'

! Binary to be read in during analysis

  ALLOCATE( eig_temp(neig) , psi_temp(td_dim,neig) )
  DO ipt=1,neig
    eig_temp(ipt)=eig(ipt)
    psi_temp(:,ipt)=psi(:,ipt)
  ENDDO

   filename='delves_surface_'//l2//'.bin'
   OPEN(UNIT=bin_unit0,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')
   WRITE(bin_unit0) jval ! Use this to make sure the correct data is read in later
   WRITE(bin_unit0) neig
   WRITE(bin_unit0) td_dim
   WRITE(bin_unit0) eig_temp
   WRITE(bin_unit0) psi_temp
   CLOSE(bin_unit0,STATUS='keep')

  DEALLOCATE( eig_temp , psi_temp )

!=========================================================================================
!  Check normalization of eigenfunctions

  vibloop: DO i=1,neig
    mag=0.d0
    DO ipt=1,td_dim
      psisq(ipt)=psi(ipt,i)*psi(ipt,i)
    ENDDO

    CALL SimpExtInt(psisq,td_dim,delta_td,mag)

    IF (ABS(mag-1.d0).gt.tol) THEN ! Check that it is normalized
      PRINT*,'NORMALIZATION ERROR: DELVES SURFACE FUNCTION:' 
      PRINT*,'JVAL =', jval ,'     IVIB = ',i
      PRINT*,'1.0-NORM = ',1.d0-mag
    ENDIF
  ENDDO vibloop

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

  DEALLOCATE( psi, psisq, eig)
  DEALLOCATE(kinetic)  
  DEALLOCATE(toeplitz)  
  DEALLOCATE(dafx)  
  DEALLOCATE(hermit)

  IF (jval.eq.jmax) DEALLOCATE(kineticd)

  PRINT*,'End Routine delves_surface for j = '//l2
  PRINT*,' '
 
END SUBROUTINE delves_surface