SUBROUTINE za_matrix(ieng, jsfdim, a_mat, whichir)

!=========================================================================================
! Written by: Jeff Crawford
!
! This routine obtains the A matrix used to find the S matrix by solving the matrix 
!    equation AS=F.
! The A matrix contains the overlap integrals of the Delves'theta with the mass-scaled 
!    Jacobi S and s eigenfunctions.
! There will be a matrix for each energy grid point and each matrix will contain elements
!   A_{fn}, where f is the final state quantum number of the Delves' functions and n is 
!   the quantum number of the Jacobi functions

! The total number of vib fcns is determined by the total energy

USE Numeric_Kinds_Module
USE CommonInfo_Module
USE EnergyGrid_Module
USE Jacobi_Module
USE APH_Module
USE Print_Module
USE SymGroup_Module
USE QuantumNumber_Module
USE SurfaceAPH_Module
USE SurfaceJacobi_Module
IMPLICIT NONE

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

INTEGER,INTENT(IN) :: ieng, jsfdim, whichir

!=========================================================================================
!   O U T P U T

COMPLEX(dp) :: a_mat(totnaphsf,jsfdim)

!=========================================================================================
!   I N T E R N A L S
CHARACTER(LEN=2)   :: l2
CHARACTER(LEN=100) :: filename

INTEGER :: lambdaval, nchir, ndim_sector, ijrot, isf, jpt
INTEGER :: lcheck, neig_aph, ndim_sec

REAL(dp)    :: e_total, probamp
COMPLEX(dp) :: a_element

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

REAL(dp), ALLOCATABLE    :: sfunction(:,:), senergy(:)
COMPLEX(dp), ALLOCATABLE :: wf_j(:)

!=========================================================================================
!   F U N C T I O N S
CHARACTER(LEN=2) :: label2

!=========================================================================================
! Initialize
!-----------------------------------------------------------------------------------------
  a_mat=0.d0
  a_element=0.d0
  e_total=e_vals(ieng)

  nchir=nchi_ir(whichir)
  ndim_sector=ntheta*nchir

 ALLOCATE(wf_j(ndim_sector))

!=========================================================================================
! READ APH surface functions and jacobi vibrational functions
!-----------------------------------------------------------------------------------------

DO lambdaval=0,jtotal ! (DO 1) aph lambda loop
    
  l2=label2(lambdaval)
  !=======================================================================================
  ! Obtain APH surface functions and energies (one lambda at a time)
  !---------------------------------------------------------------------------------------
  filename='aph_surface_'//l2//'.bin'
  OPEN(UNIT=bin_unit0,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted')
  READ(bin_unit0) lcheck
  IF (lcheck.ne.lambdaval) STOP 'wrong lambdaval readin : a_matrix'
  READ(bin_unit0) neig_aph
  READ(bin_unit0) ndim_sec
  IF (ndim_sec.ne.ndim_sector) STOP 'wrong ndim_sec for surface functions in project_aph'
  ALLOCATE( sfunction(ndim_sector,neig_aph), senergy(neig_aph))
  READ(bin_unit0) senergy
  READ(bin_unit0) sfunction
  CLOSE(bin_unit0)

  DO isf=1,naphsf(lambdaval) ! (DO 2) aph sf loop

    !===================================================================================
    ! Open file containing symmetrized jacobi surface functions for irrep whichir
    !-----------------------------------------------------------------------------------
    filename='jacobi_sym.bin'    
    OPEN(UNIT=bin_unit7,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')      

    DO jpt=1,jsfdim
  
      READ(bin_unit7) ijrot
      READ(bin_unit7) wf_j

      CALL overlap_aph(wf_j, sfunction(:,isf), ndim_sector, a_element, whichir, nchir)
           
      a_mat(isf,jpt)=a_element

    ENDDO
    
    CLOSE(bin_unit7)
 !----------------------------------------------------------------------------------------

  ENDDO !  (DO 2) aph sf loop
ENDDO ! (DO 1) aph lambda loop


DEALLOCATE(wf_j,sfunction,senergy)

END SUBROUTINE za_matrix