SUBROUTINE combine_abc(jsfdim, j_vector, njbf, whichir)
!
!=========================================================================================
!
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
IMPLICIT NONE

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

INTEGER,INTENT(IN) :: jsfdim, whichir

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

INTEGER :: j_vector(jsfdim)

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

CHARACTER(LEN=100):: filename
LOGICAL  :: readabc
INTEGER  :: ipt, j_a, j_b, j_c, nchir, ndim_sector, rcount
INTEGER  :: njbf(nchanls) 

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

COMPLEX(dp),ALLOCATABLE :: wf_a(:), wf_b(:), wf_c(:), wf_ir(:)

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

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

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

readabc=.true.
rcount=0
nchir=nchi_ir(whichir)
ndim_sector=ntheta*nchir


ALLOCATE( wf_a(ndim_sector), wf_b(ndim_sector), wf_c(ndim_sector) )
ALLOCATE( wf_ir(ndim_sector) )

!=========================================================================================
! Open files to obtain and store complete symmetrized jacobi surface functions
!-----------------------------------------------------------------------------------------

 filename='jacobi_sym_c1.bin'
 OPEN(UNIT=bin_unit1,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')

 filename='jacobi_sym_c2.bin'
 OPEN(UNIT=bin_unit2,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')

 filename='jacobi_sym_c3.bin'
 OPEN(UNIT=bin_unit3,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')

 filename='jacobi_sym.bin'
 OPEN(UNIT=bin_unit4,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')

 IF (nirrep.eq.8) THEN !(IF 1) nirrep for C6v

   DO ipt=1,njsf_ir !(DO 1)
   !DO ipt=1,njbf(2)
     IF (whichir.le.4) THEN !(IF 2) A1, A2, B1, B2------------------------------------------
       READ(bin_unit1) j_a
       READ(bin_unit2) j_b
       READ(bin_unit3) j_c
    
       READ(bin_unit1) wf_a
       READ(bin_unit2) wf_b
       READ(bin_unit3) wf_c

       IF (j_a.ne.j_b.or.j_a.ne.j_c) THEN
         PRINT*,'Error: incorrect mixing of j states'
         STOP 'combine_abc'
       ENDIF

       wf_ir=(wf_a + wf_b + wf_c)/Sqrt(3.d0)

       WRITE(bin_unit4) j_b
       WRITE(bin_unit4) wf_ir

     ELSE                   !(IF 2) E1, E2--------------------------------------------------
 
       IF (readabc.and.rcount.lt.njsf_a) THEN
       !IF (readabc.and.rcount.lt.njbf(1)) THEN
         READ(bin_unit1) j_a
         READ(bin_unit2) j_b
         READ(bin_unit3) j_c
       ELSE
         READ(bin_unit2) j_b
         READ(bin_unit3) j_c
       ENDIF

       IF (j_a.eq.j_b.and.j_a.eq.j_c) THEN ! (IF 2) when j_a = j_b = j_c

         readabc=.true.
         rcount=rcount+1
         READ(bin_unit1) wf_a
         READ(bin_unit2) wf_b
         READ(bin_unit3) wf_c

         wf_ir=(2.d0*wf_a - wf_b - wf_c)/Sqrt(6.d0)

         WRITE(bin_unit4) j_b
         WRITE(bin_unit4) wf_ir 

       ELSE                               ! (IF 2) when j_a does not exist

         readabc=.false.
         READ(bin_unit2) wf_b
         READ(bin_unit3) wf_c

         wf_ir=(wf_b - wf_c)/Sqrt(2.d0)

         WRITE(bin_unit4) j_b
         WRITE(bin_unit4) wf_ir

       ENDIF                             ! (IF 2)

     ENDIF                  !(IF 1)---------------------------------------------------------

     j_vector(ipt)=j_b

   ENDDO           ! (DO 1)

 ELSEIF (nirrep.eq.4) THEN ! (IF 1) C2v

      DO ipt=1,njsf_a
      !DO ipt=1,njbf(1)
        READ(bin_unit1) j_a
        READ(bin_unit1) wf_a
        WRITE(bin_unit4) j_a
        WRITE(bin_unit4) wf_a
        j_vector(ipt)=j_a
      ENDDO

      DO ipt=1,njsf_ir
      !DO ipt=1,njbf(2)
        READ(bin_unit2) j_b
        READ(bin_unit2) wf_b
        READ(bin_unit3) j_c
        READ(bin_unit3) wf_c

        IF (j_b.ne.j_c) THEN
         PRINT*,'Error: incorrect mixing of j states'
         STOP 'combine_abc'
       ENDIF

       IF (mod(j_b,2).eq.0) THEN

         wf_ir=(wf_b + wf_c)/Sqrt(2.d0)

       ELSE

         wf_ir=(wf_b - wf_c)/Sqrt(2.d0)

       ENDIF

       WRITE(bin_unit4) j_b
       WRITE(bin_unit4) wf_ir
       j_vector(ipt+njsf_a)=j_b
       !j_vector(ipt+njbf(1))=j_b
      ENDDO

 ELSE ! (IF 1) C2

      DO ipt=1,njbf(1)
        READ(bin_unit1) j_a
        READ(bin_unit1) wf_a
        WRITE(bin_unit4) j_a
        WRITE(bin_unit4) wf_a
        j_vector(ipt)=j_a
      ENDDO      
      DO ipt=1,njbf(2)
        READ(bin_unit2) j_b
        READ(bin_unit2) wf_b
        WRITE(bin_unit4) j_b
        WRITE(bin_unit4) wf_b
        j_vector(ipt+njbf(1))=j_b
      ENDDO
      DO ipt=1,njbf(3)
        READ(bin_unit3) j_c
        READ(bin_unit3) wf_c
        WRITE(bin_unit4) j_c
        WRITE(bin_unit4) wf_c
        j_vector(ipt+njbf(1)+njbf(2))=j_c
      ENDDO

 ENDIF ! (IF 1)

DEALLOCATE( wf_a, wf_b, wf_c, wf_ir )

CLOSE(bin_unit1)
CLOSE(bin_unit2)
CLOSE(bin_unit3)
CLOSE(bin_unit4)

END SUBROUTINE combine_abc