SUBROUTINE jacobi_basis(jval, ichanl, sdim )
!
!=========================================================================================
! This routine obtains asymptotic diatomic eigenfunctions for all arrangement channels  
!   as a function of Jacobi mass-scaled s.
!
! Variables:
!   eig      = array containing diatomic eigenvalues
!   ichanl   = arrangement channel
!   jval     = diatomic rotational quantum number
!   kineticj = array contianing the Jacobi small s kinetic energy matrix
!   neig     = number of eigenvalues
!   psi      = array containing diatmic basis functions
!   rc       = location of the potential maximum in Jacobi small s coords
!   re       = location of the potential minimum in Jacobi small s coords
!   sdim     = dimension of Jacobi small s grid
!   smax     = Jacobi small s grid maximum
!   smin     = Jacobi small s grid minimum
!   svals    = array containing Jacobi small s grid
!   v0re     = minimum value of the diatomic potential energy
!   v0rc     = maximum value of the diatomic potential energy
!   vib_pot  = array containing the diatomic potential energy (including centrifugal part)
!
!=========================================================================================
Use Numeric_Kinds_Module
USE CommonInfo_Module
USE Convrsns_Module
USE SurfaceJacobi_Module
USE QuantumNumber_Module
USE Print_Module
USE Masses_Module
IMPLICIT NONE

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

INTEGER,INTENT(IN) :: jval

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

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

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

CHARACTER(LEN=1)  :: c1
CHARACTER(LEN=2)  :: l2
CHARACTER(LEN=100):: filename
INTEGER  :: ivib, i, ipt, istat, kmax, ichanl, neig, lwork, liwork, info, nu, sdim
REAL(dp) :: e, mag, s, v, v1, v2, re, v0re, v1re, v2re, rc, v0rc, v1rc, v2rc, smin, smax

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

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

!=========================================================================================
!   F U N C T I O N S

CHARACTER(LEN=1) :: label1
CHARACTER(LEN=2) :: label2
save
!=========================================================================================
! Obtain file labels

 l2=label2(jval)
 c1=label1(ichanl)

!=========================================================================================
! Allocate Arrays and Initialize Variables

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

  liwork=MAX(1,5*(sdim)) 
  lwork=8*sdim           
  
  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

!=========================================================================================
! Find grid maximum (finite value considered to be infinity).
! Calculate coordinate grid

  IF (firstcallj) THEN

    IF (ichanl.eq.1) ALLOCATE(vib_eig_in(sdim))

    ALLOCATE( svals(sdim), vib_pot(sdim), vp(sdim) )

    kinetic=0.d0
    ALLOCATE( kineticj(sdim,sdim) )

    CALL findvasy(vasyj,smax,ichanl)

    smin=0.d0
    deltas=(smax-smin)/(sdim-1)

    DO ipt=1,sdim
      svals(ipt)=smin+deltas*(ipt-1)
    ENDDO

    CALL get_daf_s(deltas,mval,dafx,hermit,ndaf,nd_daf,toeplitz,usys2,kmax)

    CALL kin_radial_s(sdim,toeplitz,kinetic,nd_daf,kmax)

    kineticj=kinetic

    firstcallj=.false.

  ELSE

    kinetic=kineticj

  ENDIF

!=========================================================================================
! Calculate diatomic potential grid

    DO ipt=1,sdim
      s=svals(ipt)
      CALL pot(s, v, v1, v2, jval, usys, ichanl)
      If (v.ne.v) v=4.6d0
      vib_pot(ipt)=v 
      vp(ipt)=v1
    ENDDO

!=========================================================================================
! DEBUG: Output potential data

 IF (debug) THEN
   filename = 'potential/diatomic_v_c'//c1//'_'//l2//'.txt'
   OPEN(UNIT=output_unit0,FILE=TRIM(outdir)//TRIM(filename),IOSTAT=istat,STATUS='replace',ACTION='write')
     IF (istat.ne.0) STOP 'OPEN failed - surface_functions'

   DO ipt=1,sdim
     WRITE(output_unit0,*) svals(ipt),vib_pot(ipt)
   ENDDO

   CLOSE(UNIT=output_unit0,IOSTAT=istat,STATUS='keep')
     IF (istat.ne.0) STOP 'CLOSE failed - surface_functions'
 ENDIF
!========================================================================================
! Find critical points of potential (and output?)

 CALL critical(svals, vp, sdim, jval, vasyj, re, v0re, v1re, v2re, rc, v0rc, v1rc, v2rc, usys,ichanl)

! PRINT*,''
! PRINT*,'JACOBI CRITICAL POINTS'
! PRINT*,''
! WRITE(*,'(1X,3(a,ES13.6,3X))') 're   = ', re ,  'rc   = ',rc  ,'smax = ',smax
! 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*,''
! 21 FORMAT()

  DO i=1,sdim 
    kinetic(i,i)=kinetic(i,i)+vib_pot(i)
  ENDDO

!=========================================================================================
! Diagonalize with dsyevx
  
  CALL dsyevx(Jobz, Range, UpLo, sdim, kinetic, sdim, v0re, vasyj, 1, 1, AbsTol, neig, eig, psi, sdim, work, lwork, iwork, ifail, info)
 
!========================================================================================= 
! Error message check for dsyevx

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

  DEALLOCATE(ifail,iwork,work)

!=========================================================================================
! divide by sqrt(deltas) to provide proper amplitude (for integration)
 
  psi=psi/Sqrt(deltas)

!=========================================================================================
! Save number of bound states for current j value
! Place initial state info into arrays to be used later (if applicable)
! nvib_in+1 allows for nvib_in to start at zero

    IF (ichanl.eq.1.and.jval.eq.j_in) THEN
      vib_eig_in=psi(:,nvib_in+1)  !Initial rovibrational eigenfunction
      vib_energy_in=eig(nvib_in+1) !Initial eigenenergy
    ENDIF

!=========================================================================================
! Output eigenfunctions and energies to file:
IF (debug) THEN
! Energies
  filename = 'surface_functions/jacobi/energy_c'//c1//'_'//l2//'.txt'
  OPEN(UNIT=output_unit2,FILE=TRIM(outdir)//TRIM(filename),IOSTAT=istat,STATUS='replace',ACTION='write')
    IF (istat.ne.0) STOP 'OPEN failed - surface_functions'
  filename = 'surface_functions/jacobi/energy_long_c'//c1//'_'//l2//'.txt'
  OPEN(UNIT=output_unit3,FILE=TRIM(outdir)//TRIM(filename),IOSTAT=istat,STATUS='replace',ACTION='write')
    IF (istat.ne.0) STOP 'OPEN failed - surface_functions'
  DO nu=1,neig
    e=eig(nu)
    WRITE(output_unit2,'(i4,e23.15)') nu-1, e
    WRITE(output_unit3,'(a,i4,a,3(es23.15,a))')   "nu=", nu-1, '  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_functions'
  CLOSE(UNIT=output_unit3,IOSTAT=istat,STATUS='keep')
  IF (istat.ne.0) STOP 'CLOSE failed - surface_functions'

! Wave Functions
    filename = 'surface_functions/jacobi/function_c'//c1//'_'//l2//'.txt'
    OPEN(UNIT=output_unit4,FILE=TRIM(outdir)//TRIM(filename),IOSTAT=istat,STATUS='replace',ACTION='write')
    IF (istat.ne.0) STOP 'OPEN failed - jacobi_surface'
    DO ipt=1,sdim
      WRITE(output_unit4,'(1x,e14.7,1000(1X,e14.7))') svals(ipt), ( psi(ipt,ivib), ivib=1,neig )
    ENDDO
    CLOSE(UNIT=output_unit4,IOSTAT=istat,STATUS='keep')
    IF (istat.ne.0) STOP 'CLOSE failed - jacobi_surface'
ENDIF
!=========================================================================================
! Create binaries to be read in during analysis

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

   filename = 'jacobi_surface_c'//c1//'_'//l2//'.bin'
   OPEN(UNIT=bin_unit0,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')
   WRITE(bin_unit0) ichanl
   WRITE(bin_unit0) jval ! Use this to make sure the correct data is read in later
   WRITE(bin_unit0) neig
   WRITE(bin_unit0) sdim
   WRITE(bin_unit0) svals
   WRITE(bin_unit0) smin
   WRITE(bin_unit0) smax
   WRITE(bin_unit0) psi_temp
   CLOSE(bin_unit0,STATUS='keep')

   filename = 'jacobi_energy_c'//c1//'_'//l2//'.bin'
   OPEN(UNIT=bin_unit0,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')
   WRITE(bin_unit0) neig
   WRITE(bin_unit0) eig_temp
   CLOSE(bin_unit0,STATUS='keep')

  DEALLOCATE( eig_temp , psi_temp )

!=========================================================================================
!  Check normalization of initial state function and all analysis state functions

  vibloop: DO i=1,neig
    mag=0.d0
    DO ipt=1,sdim
      mag = mag + psi(ipt,i)*psi(ipt,i)*deltas
    ENDDO

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

!=========================================================================================
! Deallocate arrays

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

  IF (jval.eq.jmax) DEALLOCATE(svals, vib_pot, vp, kineticj)

END SUBROUTINE jacobi_basis