SUBROUTINE zget_sym_jacobi(ieng, njvib_open, njbf, whichir)
!
!=========================================================================================
! This routine obtains the symmetrized channel eigenfunctions. The symmetrized functions 
!   contain all the irreducible representations that comprise each channel function.
!
Use Numeric_Kinds_Module
USE CommonInfo_Module
USE SurfaceJacobi_Module
USE SurfaceDelves_Module
USE SurfaceAph_Module
USE QuantumNumber_Module
USE Print_Module
USE SymGroup_Module
USE APH_Module
USE EnergyGrid_Module
USE Masses_Module
USE Jacobi_Module
IMPLICIT NONE

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

INTEGER,INTENT(IN) :: ieng, whichir, njvib_open(0:jmax,req_chanls)

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

CHARACTER(LEN=1)   :: c1, creq
CHARACTER(LEN=2)   :: j2
CHARACTER(LEN=3)   :: v3
CHARACTER(LEN=100) :: filename
LOGICAL  :: bdelete 
INTEGER  :: indexchi, ichanl, ivib, ijrot, ispt, itc, jbegin, jend, jstep, chanlsum
INTEGER  :: irpt,nchir,ndim_sector,irho,itheta,ichi, tcdim, neig, neig_open, prevsum
INTEGER  :: njbf(nchanls)
REAL(dp) :: eval, k_out, chi_sum, theta_sum, fac

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

REAL(dp),ALLOCATABLE    :: jac_energy(:), symmat(:,:), symmat_ir(:,:),wfread(:)
COMPLEX(dp),ALLOCATABLE :: wftran(:), wfvibrot(:), resultvec(:), splitvec(:),wf_ir(:)

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

CHARACTER(LEN=1) :: label1
CHARACTER(LEN=2) :: label2
CHARACTER(LEN=3) :: label3
INTEGER :: aphindex, aphindex_tc

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

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

indexchi=indexchi_array(whichir)
nchir=nchi_ir(whichir)
ndim_sector=ntheta*nchir

tcdim=ntheta*nchi
eval=e_vals(ieng)

ALLOCATE( wfvibrot(tcdim), wftran(tcdim),wfread(tcdim))
ALLOCATE( symmat(nchi,nchi), symmat_ir(nchir,nchi))
ALLOCATE( resultvec(nchir), splitvec(nchi) )
ALLOCATE( wf_ir(ndim_sector))
   
!=========================================================================================
! Read in symmetrization matrix and pick off relevant portion of symmat
!-----------------------------------------------------------------------------------------

filename='symmat.bin'
OPEN(unit=bin_unit0,file=TRIM(BinOutdir)//TRIM(filename),form='unformatted')
READ(bin_unit0) symmat
CLOSE(bin_unit0)

DO ichi=1,nchir
  symmat_ir(ichi,:)=symmat(:,ichi+indexchi)
ENDDO

DEALLOCATE(symmat)

!=========================================================================================
! Obtain symmetrized jacobi surface functions
!-----------------------------------------------------------------------------------------

totnjsf=0
prevsum=0

DO ichanl=1,nchanls ! (1 DO) chanl loop 

  CALL jacobi_grid(ichanl) ! Switch to appropriate Jacobi channel coordinates

  chanlsum=0
  c1=label1(ichanl)

  !=======================================================================================
  ! Open file to store complete symmetrized jacobi surface functions
  !---------------------------------------------------------------------------------------

  filename='jacobi_sym_c'//c1//'.bin'
  OPEN(UNIT=bin_unit7,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')

  CALL jchoose( ichanl, whichir, jbegin, jend, jstep )

  DO ijrot = jbegin,jend,jstep ! (2 DO) rot loop

    j2=label2(ijrot)

    IF (ichanl.gt.req_chanls) THEN ! (1 IF) 
      neig_open=njvib_open(ijrot,req_chanls)
      creq=label1(req_chanls)
    ELSE
      neig_open=njvib_open(ijrot,ichanl)
      creq=label1(ichanl)
    ENDIF                          ! (1 IF)

    totnjsf=totnjsf+neig_open

    IF (neig_open.gt.0) THEN ! (2 IF)
      
      !=====================================================================================
      ! Obtain Jacobi rovibrational eigenenergies from file
      !-------------------------------------------------------------------------------------
      filename = 'jacobi_energy_c'//creq//'_'//j2//'.bin'
      OPEN(UNIT=bin_unit0,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')
      READ(bin_unit0) neig
      IF (ALLOCATED(jac_energy)) DEALLOCATE(jac_energy)
      ALLOCATE(jac_energy(neig))
      READ(bin_unit0) jac_energy
      CLOSE(bin_unit0,STATUS='keep')
      !=====================================================================================
      ! Open file containing Jacobi rovibrational eigenfunctions
      !-------------------------------------------------------------------------------------
      filename='jacobi_totsurface_c'//c1//'_'//j2//'.bin'
      OPEN(UNIT=bin_unit8,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')

      DO ivib=1,neig_open ! (3 DO) vib loop
        
        chanlsum = chanlsum + 1

        v3=label3(ivib)
        k_out=Sqrt(usys2*(eval-jac_energy(ivib)))

        CALL translational(k_out, tcdim, wftran, ijrot, ivib)
        READ(bin_unit8) wfread ! read in smalls * captheta function

        wfvibrot=Cmplx(wfread,0.d0,dp)
        wfvibrot=wfvibrot*wftran

        irho=ind_rho_infty
        DO itheta=1,ntheta
          DO ichi=1,nchi
            itc=aphindex_tc(nchi,itheta,ichi)
            ispt=aphindex(ntheta,nchi,irho,itheta,ichi,.false.)
            wfvibrot(itc)=(1.d0/Sqrt(2.d0))*(1.d0/(smalls(ispt)*larges(ispt)))*wfvibrot(itc)*Sqrt(weights(itheta))
          ENDDO
        ENDDO

       !====================================================================================
       ! Symmetrize each channel surface function
       !------------------------------------------------------------------------------------

        DO itheta=1,ntheta
          DO ichi=1,nchi
            itc=aphindex_tc(nchi,itheta,ichi)
            splitvec(ichi)=wfvibrot(itc)
          ENDDO
          resultvec=Matmul(symmat_ir,splitvec)
          DO ichi=1,nchir
            irpt=aphindex_tc(nchir,itheta,ichi)
            wf_ir(irpt)=resultvec(ichi)
          ENDDO
        ENDDO

        WRITE(bin_unit7) ijrot
        WRITE(bin_unit7) wf_ir

      ENDDO ! (3 DO) vib loop         

      CLOSE(bin_unit8)
    
    ENDIF ! (2 IF) neig_open

  ENDDO ! (2 DO) rot loop

  CLOSE(bin_unit7)
 
  njbf(ichanl)=chanlsum

  IF (chanlsum.gt.prevsum) prevsum=chanlsum
  IF (ichanl.eq.1) njsf_a=chanlsum

ENDDO ! (1 DO) chanl loop

 njsf_ir=prevsum

DEALLOCATE(wfvibrot,wftran,wfread)
DEALLOCATE(resultvec,splitvec)
DEALLOCATE(symmat_ir,wf_ir)


!STOP 'zget'

END SUBROUTINE zget_sym_jacobi