SUBROUTINE zsmatrix(ieng, jsfdim, j_vector, atf, njvib_max)

!=========================================================================================
! Written by: Jeff Crawford
!

USE Numeric_Kinds_Module
USE CommonInfo_Module
USE Convrsns_Module
USE Jacobi_Module
USE APH_Module
USE Delves_Module
USE Overlap_Module
USE EnergyGrid_Module
USE QuantumNumber_Module
USE SurfaceJacobi_Module
USE SurfaceAPH_Module
USE SymGroup_Module
USE Masses_Module
use SMat_Module
IMPLICIT NONE
SAVE

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

INTEGER,INTENT(IN)     :: ieng, jsfdim, j_vector(jsfdim), njvib_max(0:jmax,req_chanls)
COMPLEX(dp),INTENT(IN) :: atf(jsfdim)

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

CHARACTER(LEN=1)   :: ir1
CHARACTER(LEN=100) :: filename
INTEGER     :: ipt, jpt, ichanl, pir, irpt, lastj, shift, vibcount, jsfdim_p
REAL(dp)    :: sum, mag, eval
COMPLEX(dp) :: ss

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

INTEGER,ALLOCATABLE     :: j_vector_p(:)
REAL(dp),ALLOCATABLE    :: s_react(:), sr_real(:), sr_imag(:)
REAL(dp),ALLOCATABLE    :: s_nonreact(:), snr_real(:), snr_imag(:)
COMPLEX(dp),ALLOCATABLE :: atf_p(:)

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

CHARACTER(LEN=1) :: label1

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

 IF (firstcall) THEN
   DO irpt=1,4
     IF (to_propagate(irpt)) pir=irpt
   ENDDO

   ir1=label1(pir)

   !======================================================================================
   ! Open file containg s matrix from first propagated irred. rep. (A1 A2 B1 or B2)
   !--------------------------------------------------------------------------------------
   filename='smatrix_'//ir1//'.bin'
   OPEN(UNIT=smat_unit1,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted') 

 ENDIF

 READ(smat_unit1) jsfdim_p

 ALLOCATE(j_vector_p(jsfdim_p), atf_p(jsfdim_p))

 READ(smat_unit1) j_vector_p
 READ(smat_unit1) atf_p

 sum=0.d0

 ALLOCATE(s_react(react_dim),s_nonreact(nonreact_dim))
 ALLOCATE(sr_real(react_dim),sr_imag(react_dim))
 ALLOCATE(snr_real(nonreact_dim),snr_imag(nonreact_dim))
 s_react=0.d0
 s_nonreact=0.d0
 sr_real=0.d0
 sr_imag=0.d0
 snr_real=0.d0
 snr_imag=0.d0

   DO ichanl=1,nchanls-1 ! (DO 1: chanls loop)
     ipt=1
   
     IF (ichanl.eq.1) THEN
       lastj=Mod(j_in,2)
     ELSE
       lastj=0
     ENDIF

     shift=0
     vibcount=1

     DO jpt=1,jsfdim

         !============================================================================
         ! Channel A
         !----------------------------------------------------------------------------
         IF (ichanl.eq.1) THEN

           IF (j_vector(jpt).eq.j_vector_p(ipt)) THEN

             ss=( atf_p(ipt)+(2.d0*atf(jpt)/Sqrt(2.d0)) )/Sqrt(3.d0)
             mag=Real(ss)**2+Aimag(ss)**2
             ! place into smatrix array
             IF (j_vector_p(ipt).ne.lastj) THEN 
               shift=shift+njvib_max(lastj,1)
               vibcount=1
             ENDIF
             s_nonreact(vibcount+shift)=mag
             snr_real(vibcount+shift)=REAL(ss)
             snr_imag(vibcount+shift)=AIMAG(ss)
             lastj=j_vector_p(ipt)
             vibcount=vibcount+1
             !
             sum=sum+mag
             IF (ipt.lt.jsfdim_p) THEN
               ipt=ipt+1
             ELSE
               ipt=jsfdim_p
             ENDIF
           ENDIF
         !============================================================================
         ! Channel B
         !----------------------------------------------------------------------------
         ELSEIF (ichanl.eq.2) THEN

           IF (j_vector(jpt).eq.j_vector_p(ipt)) THEN
             ss=(atf_p(ipt)/Sqrt(3.d0))-(atf(jpt)/Sqrt(6.d0))
             mag=Real(ss)**2+Aimag(ss)**2
             sum=sum+mag
             IF (ipt.lt.jsfdim_p) THEN
               ipt=ipt+1
             ELSE
               ipt=jsfdim_p
             ENDIF
           ELSE
             ss=(atf(jpt)/Sqrt(2.d0))
             mag=Real(ss)**2+Aimag(ss)**2
             sum=sum+mag
           ENDIF
           IF (j_vector(jpt).ne.lastj) THEN
               shift=shift+njvib_max(lastj,1)
               vibcount=1
           ENDIF
             s_react(vibcount+shift)=2.d0*mag
             sr_real(vibcount+shift)=REAL(ss)
             sr_imag(vibcount+shift)=AIMAG(ss)
             lastj=j_vector(jpt)
             vibcount=vibcount+1

         !============================================================================
         ! Channel C
         !----------------------------------------------------------------------------
         ELSEIF (ichanl.eq.3) THEN

           IF (j_vector(jpt).eq.j_vector_p(ipt)) THEN
             ss=(atf_p(ipt)/Sqrt(3.d0))-(atf(jpt)/Sqrt(6.d0))
             sum=sum+Real(ss)**2+Aimag(ss)**2
             IF (ipt.lt.jsfdim_p) THEN
               ipt=ipt+1
             ELSE
               ipt=jsfdim_p
             ENDIF
           ELSE
             ss=-(atf(jpt)/Sqrt(2.d0))
             mag=Real(ss)**2+Aimag(ss)**2
             sum=sum+mag
           ENDIF

         ENDIF ! channel conditional
   
     ENDDO        
   ENDDO     !(DO 1): end ichanl loop

!=======================================================================================
! Output to File
!---------------------------------------------------------------------------------------

  IF (firstcall) THEN
   filename='s_nonreact.txt'
   OPEN(UNIT=smatrix_unit0,FILE=TRIM(outdir)//TRIM(filename),STATUS='unknown')
   filename='s_react.txt'
   OPEN(UNIT=smatrix_unit1,FILE=TRIM(outdir)//TRIM(filename),STATUS='unknown')
   filename='snr_real.txt'
   OPEN(UNIT=smatrix_unit2,FILE=TRIM(outdir)//TRIM(filename),STATUS='unknown')
   filename='snr_imag.txt'
   OPEN(UNIT=smatrix_unit3,FILE=TRIM(outdir)//TRIM(filename),STATUS='unknown')
   filename='sr_real.txt'
   OPEN(UNIT=smatrix_unit4,FILE=TRIM(outdir)//TRIM(filename),STATUS='unknown')
   filename='sr_imag.txt'
   OPEN(UNIT=smatrix_unit5,FILE=TRIM(outdir)//TRIM(filename),STATUS='unknown')
!   filename='smatrix_sum.txt'
!   OPEN(UNIT=smatrix_unit6,FILE=TRIM(outdir)//TRIM(filename),STATUS='unknown')
  ENDIF


  eval=e_vals(ieng)*autoev
  WRITE(smatrix_unit0,'(1X,e14.7,1000(1X,e14.7))') eval, (s_nonreact(ipt),ipt=1,nonreact_dim)
  WRITE(smatrix_unit1,'(1X,e14.7,1000(1X,e14.7))') eval, (s_react(ipt),ipt=1,react_dim)
  WRITE(smatrix_unit2,'(1X,e14.7,1000(1X,e14.7))') eval, (snr_real(ipt),ipt=1,nonreact_dim)
  WRITE(smatrix_unit3,'(1X,e14.7,1000(1X,e14.7))') eval, (snr_imag(ipt),ipt=1,nonreact_dim)
  WRITE(smatrix_unit4,'(1X,e14.7,1000(1X,e14.7))') eval, (sr_real(ipt),ipt=1,react_dim)
  WRITE(smatrix_unit5,'(1X,e14.7,1000(1X,e14.7))') eval, (sr_imag(ipt),ipt=1,react_dim)
 ! WRITE(smatrix_unit6,'(1X,e14.7,e14.7)') eval, sum
  IF (ieng.eq.edim) THEN
    CLOSE(smatrix_unit0)
    CLOSE(smatrix_unit1)
    CLOSE(smatrix_unit2)
    CLOSE(smatrix_unit3)
    CLOSE(smatrix_unit4)
    CLOSE(smatrix_unit5)
  !  CLOSE(smatrix_unit6)
  ENDIF

   DEALLOCATE(j_vector_p)
   DEALLOCATE(atf_p)
   DEALLOCATE(s_react,sr_real,sr_imag)
   DEALLOCATE(s_nonreact,snr_real,snr_imag)

  ! If (ieng.eq.edim) CLOSE(smat_unit1)

firstcall=.false.

END SUBROUTINE zsmatrix