SUBROUTINE smatrixABC(ieng, njvib_max, maxopen, 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, jpt, ichanl, lastj, vshift, vibcount, chshift, jbfdim
INTEGER  :: jbf(nchanls), maxopen(req_chanls)
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_a(:), s_b(:), s_c(:)
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. (A or B)
   !--------------------------------------------------------------------------------------
   filename='smatrix_'//ir1//'.bin'
   OPEN(UNIT=smat_unit2,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted')

 ENDIF

 READ(smat_unit2) jbf(1)
 READ(smat_unit2) jbf(2)
 READ(smat_unit2) jbf(3)
 READ(smat_unit2) jbfdim

 ALLOCATE(j_vector(jbfdim), atf(jbfdim))

 READ(smat_unit2) j_vector
 READ(smat_unit2) atf

 sumj=0.d0
 sumtot=0.d0

ALLOCATE( s_a(maxopen(1)), s_b(maxopen(2)), s_c(maxopen(3)) )

 s_a = 0.d0
 s_b = 0.d0
 s_c = 0.d0
 chshift=0

DO ichanl=1,nchanls

  vibcount=1
  vshift=0
  lastj=0

  DO ipt=1,jbf(ichanl)

    jpt=ipt+chshift

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

    IF (j_vector(jpt).ne.lastj) THEN
      vshift=vshift+njvib_max(lastj,ichanl)
      vibcount=1
    ENDIF
    IF (ichanl.eq.1) s_a(vibcount+vshift)=ss
    IF (ichanl.eq.2) s_b(vibcount+vshift)=ss
    IF (ichanl.eq.3) s_c(vibcount+vshift)=ss 
    lastj=j_vector(jpt)
    vibcount=vibcount+1
    sumtot=sumtot+ss
  ENDDO

  chshift=chshift+jbf(ichanl)

ENDDO

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

  IF (firstcall) THEN
   filename='s_a2.txt'
   OPEN(UNIT=smatrix_unit0,FILE=TRIM(outdir)//TRIM(filename),STATUS='unknown')
   filename='s_b2.txt'
   OPEN(UNIT=smatrix_unit1,FILE=TRIM(outdir)//TRIM(filename),STATUS='unknown')
   filename='s_c2.txt'
   OPEN(UNIT=smatrix_unit2,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_a(ipt),ipt=1,maxopen(1))
  WRITE(smatrix_unit1,'(1X,e14.7,1000(1X,e14.7))') eval, (s_b(ipt),ipt=1,maxopen(2))
  WRITE(smatrix_unit2,'(1X,e14.7,1000(1X,e14.7))') eval, (s_c(ipt),ipt=1,maxopen(3))


  IF (ieng.eq.edim) THEN
    CLOSE(smatrix_unit0)
    CLOSE(smatrix_unit1)
    CLOSE(smatrix_unit2)
  ENDIF

   DEALLOCATE(s_a)
   DEALLOCATE(s_b)
   DEALLOCATE(s_c)
   DEALLOCATE(j_vector)
   DEALLOCATE(atf)

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

firstcall=.false.

END SUBROUTINE smatrixABC