SUBROUTINE szsmatrix(ieng, njvib_max, whichir)

!=========================================================================================
! 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, whichir, njvib_max(0:jmax,req_chanls)

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

CHARACTER(LEN=1)   :: ir1
CHARACTER(LEN=100) :: filename
INTEGER  :: ipt, lastj, shift, vibcount, jsfdim, jsf_a, jsf_ir
REAL(dp) :: sumj, sumtot, eval, ss

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

INTEGER,ALLOCATABLE     :: j_vector(:)
REAL(dp),ALLOCATABLE    :: s_react(:), s_nonreact(:), s_react2(:)
COMPLEX(dp),ALLOCATABLE :: atf(:)

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

CHARACTER(LEN=1) :: label1

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

 IF (firstcall) THEN

   ir1=label1(whichir)

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

 ENDIF

 READ(smat_unit2) jsf_a
 READ(smat_unit2) jsf_ir
 READ(smat_unit2) jsfdim

 ALLOCATE(j_vector(jsfdim), atf(jsfdim))

 READ(smat_unit2) j_vector
 READ(smat_unit2) atf

 sumj=0.d0
 sumtot=0.d0

ALLOCATE(s_react(react_dim),s_nonreact(nonreact_dim))

 s_react=0.d0
 s_nonreact=0.d0

vibcount=1
shift=0
lastj=Mod(j_in,2)

DO ipt=1,jsf_a

  ss=(Real(atf(ipt))**2)+(Aimag(atf(ipt))**2)

  IF (j_vector(ipt).ne.lastj) THEN
    shift=shift+njvib_max(lastj,1)
    vibcount=1
  ENDIF
  s_nonreact(vibcount+shift)=ss
  lastj=j_vector(ipt)
  vibcount=vibcount+1
  sumtot=sumtot+ss
ENDDO

IF (nirrep.eq.4) THEN

  vibcount=1
  shift=0
  lastj=0

  DO ipt=1,jsf_ir

    ss=(Real(atf(ipt+jsf_a)/Sqrt(2.d0))**2)+(Aimag(atf(ipt+jsf_a)/Sqrt(2.d0))**2)

    IF (j_vector(ipt+jsf_a).ne.lastj) THEN
      shift=shift+njvib_max(lastj,2)
      vibcount=1
    ENDIF
  !  print*,ipt,ss
    s_react(vibcount+shift)=2.d0*ss
    lastj=j_vector(ipt+jsf_a)
    vibcount=vibcount+1
    sumtot=sumtot+2.d0*ss
  ENDDO

ELSEIF (nirrep.eq.2) THEN

  ALLOCATE(s_react2(react_dim))

  vibcount=1
  shift=0
  lastj=0

  DO ipt=1,jsf_ir

    ss=(Real(atf(ipt+jsf_a))**2)+(Aimag(atf(ipt+jsf_a))**2)

    IF (j_vector(ipt+jsf_a).ne.lastj) THEN
      shift=shift+njvib_max(lastj,2)
      vibcount=1
    ENDIF
    s_react(vibcount+shift)=ss
    lastj=j_vector(ipt+jsf_a)
    vibcount=vibcount+1
  ENDDO

  vibcount=1
  shift=0
  lastj=0

  DO ipt=1,jsf_ir

    ss=(Real(atf(ipt+jsf_a+jsf_ir))**2)+(Aimag(atf(ipt+jsf_a+jsf_ir))**2)

    IF (j_vector(ipt+jsf_a+jsf_ir).ne.lastj) THEN
      shift=shift+njvib_max(lastj,3)
      vibcount=1
    ENDIF
    s_react2(vibcount+shift)=ss
    lastj=j_vector(ipt+jsf_a+jsf_ir)
    vibcount=vibcount+1
  ENDDO

ENDIF

!=======================================================================================
! 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='s_sumtot.txt'
   !OPEN(UNIT=smatrix_unit2,FILE=TRIM(outdir)//TRIM(filename),STATUS='unknown')
   IF (nirrep.eq.2) THEN
     filename='s_react2.txt'
     OPEN(UNIT=smatrix_unit3,FILE=TRIM(outdir)//TRIM(filename),STATUS='unknown')
   ENDIF
  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,2(e14.7))') eval, sumtot
  IF (nirrep.eq.2) THEN
    WRITE(smatrix_unit3,'(1X,e14.7,1000(1X,e14.7))') eval, (s_react2(ipt),ipt=1,react_dim)
  ENDIF


  IF (ieng.eq.edim) THEN
    CLOSE(smatrix_unit0)
    CLOSE(smatrix_unit1)
   ! CLOSE(smatrix_unit2)
    IF (nirrep.eq.2) THEN
      CLOSE(smatrix_unit3)
    ENDIF
  ENDIF

   DEALLOCATE(s_react)
   DEALLOCATE(s_nonreact)
   DEALLOCATE(j_vector)
   DEALLOCATE(atf)
   IF (nirrep.eq.2) DEALLOCATE(s_react2)

   If (ieng.eq.edim) CLOSE(smat_unit2)

firstcall=.false.

END SUBROUTINE szsmatrix