C     PROGRAM TO COMPUTE TOTAL, DEGENERACY AVERAGED CROSS SECTIONS,
C       SIG(J2,J1),  FROM S-MATRICES SAVED ON TAPE BY SCATTERING PROG.
C     14 OCT 94 (SG) MOLSCAT VERSION 14 COMPATIBLE
C
C     NAMELIST (SAVE) VARIABLES:
C       JTOTL - SKIP JTOT.LT.JTOTL
C       JTOTU - SKIP JTOT.GT.JTOTU
C       JSTEP - SKIP JTOT NOT EQUAL TO JTOTL(JSTEP)JTOTU
C       JTOUT - IF .LT. 999999 (DEFAULT) NO ECHO FOR JTOT.LT.JTOUT
C       IT    - INPUT TAPE UNIT (DEFAULT 10)
C       MXSIG - IF .GT. 0 WILL LIMIT THE NUMBER OF ACCUM SIG
C       LFMT  - FORMATTED/UNFORMATTED DATA SET
C
C *** ***** ***** ***** ***** ***** ***** ***** ***** ***** ***** ***
C *** *****  LFMT CONTROLS FORMATTED (VERSION 10 AND EARLIER)   * ***
C *** *****                UNFORMATTED (VERSION 11 AND LATER)   * ***
C *** ***** ***** ***** ***** ***** ***** ***** ***** ***** ***** ***
C
C     *** NOTE.  WORKS FOR ITYPE=1,2,5 (AND EP AND CS EQUIVALENTS)
C     *** CHANGES 4/20/92 FOR ITYPE=3
C
C     SIG(JF,JI) OUTPUT ON CARDS IN FORMAT FOR STATISTICAL EQUILIBRIUM
C       PROGRAM.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXCH=300,MXJLEV=500,MXS=MXCH*MXCH,MXSG=50000)
C     DIMS OF ELEVEL (IN COMMON/CMBASE/) AND ENERGY FIXED BY MOLSCAT
      PARAMETER (MXEL=1000,MXNRG=100)
      DIMENSION NBASIS(MXCH),J(MXCH),L(MXCH),WVEC(MXCH),JLEV(MXJLEV)
      DIMENSION ENERGY(MXNRG),MINJT(MXNRG),MAXJT(MXNRG),ISTSIG(MXNRG)
      DIMENSION SREAL(MXS),SIMAG(MXS),SIG(MXSG)
      DIMENSION ELEVEL(MXEL)
      INTEGER ITYPE,NLEV,NQN,NNRG,NOPEN
      CHARACTER*4 LABEL(20)
      CHARACTER*1 COUT,CC,CS,BL,ST,XST
      LOGICAL LOGI,LOGJ,LFMT,LFORM,LPRINT
C
C     LFORM AND LPRINT CONTROL CROSS SECTION OUTPUT
      DATA LFORM/.TRUE./, LPRINT/.FALSE./
C
      NAMELIST /SAVE/ JTOTL,JTOTU,JSTEP,JTOUT,IT,MXSIG,LFMT
C
C     DEFAULT VALUES FOR NAMELIST VARIABLES
      DATA IT/10/, MXSIG/ 0/
      DATA JTOTU/999999/,JTOTL/0/,JSTEP/1/,JTOUT/999999/
      DATA LFMT/.FALSE./
C
C     CRITERION FOR PRINTING SIGMA
      DATA EPS/1.D-8/
      DATA JSTEPT/-1/
      DATA COUT/'X'/, CC/'C'/, CS/'S'/, BL/' '/,ST/'*'/
      DATA PI/3.14159 26535 89793 D0/
C
C     FORMATS FOR SAVE TAPE (FOR VERSION 11 AND EARLIER)
  100 FORMAT(20A4/3I4,F8.4,I4)
  101 FORMAT(20I4)
  102 FORMAT(I4/(5E16.8))
  103 FORMAT(2I4,E16.8,I4,E16.8,I4)
  104 FORMAT(I4/(2I4,E16.8) )
  105 FORMAT(5E16.8)
C
C     SUPPRESS UNDERFLOWS
C     CALL XUFLOW(0)
C
C     INITIALIZATION . . .
      DO 10  I=1,MXNRG
      MINJT(I)=-1
   10 MAXJT(I)=-1
      DO 11 I=1,MXEL
   11 ENERGY(I)=0.D0
C
C     GET (POSSIBLE) NAMELIST INPUT DATA
      READ(5,SAVE,END=8000)
 8000 WRITE(6,608) IT,JTOTL,JTOTU,JSTEP,JTOUT
  608 FORMAT('  *** PROGRAM TO COMPUTE STATE-TO-STATE CROSS SECTIONS'/
     1       '      FROM S-MATRICES SAVED ON TAPE UNIT',I3/
     2       '      NAMELIST PARAMETERS ARE JTOTL =',I8/
     3       30X,'JTOTU =',I8/30X,'JSTEP =',I8/30X,'JTOUT =',I8)
      IF (MXSIG.GT.0) WRITE(6,643) MXSIG
  643 FORMAT('  *** NOTE POSSIBLE LIMIT BY MXSIG =',I6)
      IF (JTOUT.LT.999999) THEN
        WRITE(6,688) JTOUT
  688   FORMAT('  ***  INPUT ECHO SUPPRESSED FOR JTOT.LE.',I7)
      ENDIF
C
C     GET AND PROCESS HEADER RECORD FROM SAVE TAPE
      IF (LFMT) THEN
        READ(IT,100) LABEL,ITYPE,NLEV,NQN,URED,IP
      ELSE
        READ(IT) LABEL,ITYPE,NLEV,NQN,URED,IP
      ENDIF
      WRITE(6,600) LABEL,ITYPE,NLEV,NQN,URED,IP
  600 FORMAT(/' SAVE TAPE LABEL =',20A4/
     2       '  ITYPE =',I3,5X,'NLEV, NQN =',2I4,5X,'URED =',F8.4,
     3       5X,'IPROGM =',I3)
      IF (LFMT.AND.IP.GE.11) THEN
        WRITE(6,*) ' *** IPROGM INCONSISTENT WITH FORMATTED TAPE'
        STOP
      ENDIF
      COUT=BL
      IF (ITYPE.LT.10) COUT=CC
      IF (ITYPE.GE.20 .AND. ITYPE.LT.30) COUT=CS
      ITYPE=ITYPE-10*(ITYPE/10)
      IF (.NOT.(ITYPE.EQ.1.OR.ITYPE.EQ.2.OR.ITYPE.EQ.5.OR.ITYPE.EQ.6
     1    .OR.ITYPE.EQ.7.OR.ITYPE.EQ.3)) THEN
        WRITE(6,*) ' *** CODE NOT IMPLEMENTED FOR ITYPE',ITYPE
        STOP
      ENDIF
C
C     GET JLEV INFORMATION FROM TAPE
      NSQ=NLEV*NQN
      IF (NSQ.GT.MXJLEV) THEN
        WRITE(6,*) ' *** NLEV*NQN.GT.MXJLEV',MXJLEV
        STOP
      ENDIF
      IF (LFMT) THEN
        READ(IT,101) (JLEV(I),I=1,NSQ)
      ELSE
        READ(IT) (JLEV(I),I=1,NSQ)
      ENDIF
      IDENT=0
      IF (ITYPE.EQ.3) CALL IDCHK(IDENT,NLEV,JLEV)
C
C     GET NLEVEL/ELEVEL INFORMATION
      IF (LFMT) THEN
        IF (IP.GE.3) READ(IT,102) NLEVEL,(ELEVEL(I),I=1,NLEVEL)
      ELSE
         READ(IT) NLEVEL,(ELEVEL(I),I=1,NLEVEL)
      ENDIF
      DO 1002 I=1,NLEV
 1002 WRITE(6,601) I,ELEVEL(ABS(JLEV(NLEV*(NQN-1)+I))),
     1             (JLEV(NLEV*(K-1)+I),K=1,NQN)
  601 FORMAT(' LEVEL',I4,'  ENERGY',F10.3,'  QUANTUM NOS. ',10I4)
C
C     GET COLLISION ENERGIES FROM SAVE TAPE
      IF (LFMT) THEN
        READ(IT,102) NNRG,(ENERGY(I),I=1,NNRG)
      ELSE
        READ(IT) NNRG,(ENERGY(I),I=1,NNRG)
      ENDIF
      IF (NNRG.GT.MXNRG) THEN
        WRITE(6,*) ' *** NNRG.GT.MXNRG',NNRG,MXNRG
        STOP
      ENDIF
      WRITE(6,602) NNRG,(ENERGY(I),I=1,NNRG)
  602 FORMAT(/' NUMBER OF COLLISION ENERGIES IS',I4/(' ',1P,5D14.4))
C
C**** GET TOP 'LEVEL NO.' FROM JLEV(I,NQN)
      ILOW=NLEV*(NQN-1)+1
      ITOP=NLEV*NQN
      NTOP=0
      DO 1101 I=ILOW,ITOP
 1101 NTOP=MAX(NTOP,ABS(JLEV(I)))
      WRITE(6,681) NTOP,NLEV
  681 FORMAT(/'  *** NO. OF LEVELS FOR SIGMA(,) IS',I4,'  NLEV =',I4)
      IF (MXSIG.GT.0) NTOP=MIN(NTOP,MXSIG)
      NSQ=NTOP*NTOP
      ISTSIG(1)=0
      DO 1001 I=2,NNRG
 1001 ISTSIG(I)=ISTSIG(I-1)+NSQ
C     CLEAR SIG()
      ISTOP=NTOP*NTOP*NNRG
      IF (ISTOP.GT.MXSG) THEN
        WRITE(6,*) ' *** STORAGE FOR SIGMA() EXCEEDED',ISTOP,MXSG
        STOP
      ENDIF
      DO 1100 I=1,ISTOP
 1100 SIG(I)=0.
C
C     START READ SAVE TAPE LOOP OVER JTOT / ENERGY VALUES . . .
C
      JSTEPX=-1
      JTOLD=-1
C
C     FIRST, GET JTOT,INRG,M RECORD -----
 2000 IF (LFMT) THEN
        READ(IT,103,END=9000) JTOT,INRG,ECHK,IEXCH,WT,M
      ELSE
        IF (IP.LT.14) THEN
          READ(IT,END=9000) JTOT,INRG,ECHK,IEXCH,WT,M
        ELSE
          READ(IT,END=9000) JTOT,INRG,ECHK,IEXCH,WT,M,NOPEN
        ENDIF
      ENDIF
C     N.B. JTOT .LT. 0  INDICATES A OBSOLETE 'CHECKPOINT' RECORD.
      IF (JTOT.LT.0)  THEN
        WRITE(6,698)  JTOT
  698   FORMAT('0 * * * NOTE. ROLLOUT MARKER ENCOUNTERED.  NROLL =',I6)
        GO TO 2000
      ENDIF
C     ECHO INPUT
      IF (JTOUT.GE.999999.OR.(JTOUT.LT.999999.AND.JTOT.GT.JTOUT)) THEN
        IF (IP.LT.14) THEN
          WRITE(6,662) IT,JTOT,M,INRG,ECHK,IEXCH,WT
  662     FORMAT(' UNIT(',I2,'): JTOT ',I4,'.',I3,'  E(',I3,')=',
     1         F10.3,' IEX,WT=',I2,F6.2)
        ELSE
          WRITE(6,663) IT,JTOT,M,INRG,ECHK,IEXCH,WT,NOPEN
  663     FORMAT(' UNIT(',I2,'): JTOT ',I4,'.',I3,'  E(',I3,')=',
     1         F10.3,' IEX,WT,NOP=',I2,F7.2,I5)
        ENDIF
      ENDIF
C     CALCULATE JSTEP VALUE FROM TAPE (JSTEPT)
      IF (JTOT.EQ.JTOLD .OR. JTOLD.EQ.-1) GO TO 2011
      JSTEPT=JTOT-JTOLD
      IF (JSTEPT.EQ.JSTEPX) GO TO 2011
      IF (JSTEPX.EQ.-1) GO TO 2912
      WRITE(6,649) JSTEPT,JSTEPX
  649 FORMAT('0 *** ERROR JSTEPT,JSTEPX =',2I6)
      STOP
 2912 JSTEPX=JSTEPT
 2011 JTOLD=JTOT
C     CHECK THAT ENERGY(INRG) CORRESPONDS WITH HEADER RECORD. . .
      IF (ABS(ECHK-ENERGY(INRG)).LT.1.D-6) GO TO 2002
      WRITE(6,697)  INRG,ECHK
  697 FORMAT('0 * * * WARNING.  FOR ',I4,'-TH ENERGY, ECHECK =',D16.8)
C
C     GET J,L,WVEC VALUES -----
 2002 IF (LFMT) THEN
        READ(IT,104,END=9090) NOPEN,(J(I),L(I),WVEC(I),I=1,NOPEN)
        IF (NOPEN.GT.MXCH)  THEN
          WRITE(6,696)  MXCH,JTOT,INRG
  696     FORMAT('0 *** ERROR. NO. OF OPEN CHANNELS EXCEEDS',I6,
     1           '   JTOT,INRG =',2I6)
          STOP
        ENDIF
      ELSE
        IF (IP.GE.14) THEN
          IF (NOPEN.GT.MXCH) THEN
            WRITE(6,696)  MXCH,JTOT,INRG
            STOP
          ENDIF
          READ(IT,END=9090) (J(I),L(I),WVEC(I),I=1,NOPEN)
        ELSE
          READ(IT,END=9090) NOPEN,(J(I),L(I),WVEC(I),I=1,NOPEN)
          IF (NOPEN.GT.MXCH)  THEN
            WRITE(6,696)  MXCH,JTOT,INRG
            STOP
          ENDIF
        ENDIF
      ENDIF
C
C     FINALLY GET S-MATRICES ---
      NSQ=NOPEN*NOPEN
      IF (LFMT) THEN
      READ(IT,105,END=9090) (SREAL(I),I=1,NSQ)
      READ(IT,105,END=9090) (SIMAG(I),I=1,NSQ)
      ELSE
         CALL SREAD(IT,NOPEN,SREAL,IEND)
         IF (IEND.GT.0) GO TO 9090
         CALL SREAD(IT,NOPEN,SIMAG,IEND)
         IF (IEND.GT.0) GO TO 9090
      ENDIF
C
C     ALLOW FOR SKIPPING  S-MATRICES ON JTOTL,JTOTU,JSTEP INPUT.
      JDIF=JTOT-JTOTL
      IF (JTOT.LT.JTOTL .OR. JTOT.GT.JTOTU .OR.
     1    JDIF-JSTEP*(JDIF/JSTEP).NE.0) THEN
         WRITE(6,682) JTOT
  682    FORMAT(' --- SKIPPING JTOT',I6,'   NOT IN JTOTL (JSTEP) JTOTU')
         GO TO 2000
      ENDIF
C
C     ACCUMULATE TO CROSS SECTION; BOOKKEEPING FOR MINJT, MAXJT
      IF (MINJT(INRG).LT.0)  MINJT(INRG)=JTOT
      MAXJT(INRG)=MAX0(MAXJT(INRG),JTOT)
      IJ=0
C     LET II MEASURE INITIAL AND JJ MEASURE FINAL QUANT. STATES
      XJT=(2*JTOT+1)*PI
      IF (IP.LT.14) THEN
        IF (IEXCH.NE.0) XJT=XJT*WT
      ELSE
        IF (WT.NE.0.D0) XJT=XJT*WT
      ENDIF
      DMX=0.D0
      OMX=0.D0
      DO 3000 II=1,NOPEN
C     LEVI/LEVJ ARE 'LEVEL NO.'; IXI/IXJ 'SIG INDX' FROM JLEV(LEV,NQN)
C     BEGINNING IN V14, NEGATIVE SIG INDX FOR "INCOMPLETE" CS VALUES
      LEVI=J(II)
      IXI=JLEV(NLEV*(NQN-1)+LEVI)
      LOGI=IXI.LT.0
      IXI=ABS(IXI)
      DO 3000 JJ=1,NOPEN
      LEVJ=J(JJ)
      IXJ=JLEV(NLEV*(NQN-1)+LEVJ)
      LOGJ=IXJ.LT.0
      IXJ=ABS(IXJ)
      FACT=XJT/(WVEC(JJ)*WVEC(JJ))
      IJ=IJ+1
      IF (II.EQ.JJ)  SREAL(IJ)=1.D0-SREAL(IJ)
C     INDEX OF SIG(JJ,II) IS ISTSIG(INRG)+NTOP*(IXJ-1)+IXI
      IX=ISTSIG(INRG)+NTOP*(IXJ-1)+IXI
      ADD = (SREAL(IJ)*SREAL(IJ)+SIMAG(IJ)*SIMAG(IJ))*FACT
C     GET DEGENERACY DENOMINATOR
      IF (ITYPE.EQ.3) THEN
C       BELOW IS ITYPE=3 W/ TAKAYANAGI IDENTICAL PARTICLE COUNTING
        JI1=JLEV(LEVJ)
        JI2=JLEV(LEVJ+NLEV)
        DEGEN=DBLE((2*JI1+1)*(2*JI2+1))
        IF (IDENT.NE.0) THEN
          JF1=JLEV(LEVI)
          JF2=JLEV(LEVI+NLEV)
          IF (JI1.EQ.JI2) DEGEN=DEGEN/2.D0
          IF (JF1.EQ.JF2) DEGEN=DEGEN/2.D0
        ENDIF
      ELSE
C       FOLLOWING APPLIES TO LINEARS AND TO SYM/ASYM TOPS
        DEGEN=DBLE(2*JLEV(LEVJ)+1)
      ENDIF
      ADD=ADD/DEGEN
      FACT=1.D0
      IF (LOGI.AND.LOGJ) FACT=-FACT
      SIG(IX)=SIG(IX)+ADD*FACT
      IF (II.EQ.JJ) OMX=MAX(OMX,ADD)
      IF (II.NE.JJ) DMX=MAX(DMX,ADD)
 3000 CONTINUE
      IF (JTOUT.GE.999999.OR.(JTOUT.LT.999999.AND.JTOT.GT.JTOUT))
     1               WRITE(6,653) OMX,DMX
  653 FORMAT(10X,'MAX CHANGE IN DIAG, OFF-DIAG SIGMA',1P,2D12.4)
C     GO BACK FOR MORE JTOT, INRG SETS . . .
      GO TO 2000
C
C     END OF FILE ON (10)
 9000 WRITE(6,695) IT
  695 FORMAT('  *** NOTE.  NORMAL END OF FILE REACHED ON UNIT (',I3,')')
      GO TO 9123
 9090 WRITE(6,690) IT
  690 FORMAT('  ***** ***** ABNORMAL EOF REACHED ON UNIT (',I3,')')
C
C     PRINT CROSS SECTIONS.  CODE REMAINS FOR TOTALS OUT OF LEVEL
 9123 JSTOUT=JSTEP*JSTEPT
      XJSTEP=JSTOUT
      WRITE(6,123) JSTEPT,JSTEP,XJSTEP
  123 FORMAT(/'  *****'/'  JSTEP ON TAPE =',I4/
     1        '  JSTEP REQUESTED IN &SAVE INPUT =',I4/
     2        '  CROSS SECTIONS MULTIPLIED BY',F9.2/'  *****')
      IF (LFORM) WRITE(6,702)
  702 FORMAT(/12X,'ENERGY  JTOTL  JSTEP  JTOTU',9X,
     1       'F    I',10X,'SIG(F,I)'/)
      DO 4000 I=1,NNRG
      DO 4100 II=1,NTOP
      EREL=ENERGY(I)-ELEVEL(II)
      SUMII=0.
      DO 4101 JJ=1,NTOP
      IST=ISTSIG(I)+(II-1)*NTOP+JJ
      SIG(IST)=SIG(IST)*XJSTEP
C     ***OFF-DIAGONAL**** TOTAL OUT OF LEVEL II
      IF (II.NE.JJ.AND.SIG(IST).GT.0.D0) SUMII=SUMII+SIG(IST)
      XST=BL
      IF (SIG(IST).LT.0.D0) XST=ST
      IF (ABS(SIG(IST)).LT.EPS) GO TO 4101
      IF (LFORM) THEN
C       FORMAT BELOW MATCHES MOLSCAT VERSION 14 OUTPUT ROUTINE
        WRITE(6,703)BL,ENERGY(I),MINJT(I),JSTOUT,MAXJT(I),
     1              JJ,II,ABS(SIG(IST)),XST
  703   FORMAT(A1,F19.6,I5,2I7,5X,2I5,1P,D20.6,1X,A1)
      ELSE
        WRITE(6,700) ENERGY(I),MINJT(I),JSTOUT,MAXJT(I),JJ,II,SIG(IST)
     2                 ,COUT,NLEV
  700   FORMAT(' ',F19.6,I5,'(',I3,')',I5,5X,2I5,E20.6,5X,A1,I2)
      ENDIF
 4101 CONTINUE
      IF (LPRINT) WRITE(6,692) EREL,II,SUMII
  692 FORMAT('  *** REL E =',F10.2,'  SUM OUT OF LEVEL',I4,' =',F8.3)
 4100 CONTINUE
 4000 CONTINUE
      STOP
      END
      SUBROUTINE IDCHK(IDENT,NLEV,JLEV)
      DIMENSION JLEV(NLEV)
      LOGICAL J1GEJ2,J1LEJ2,J1AND2
      WRITE(6,600)
  600 FORMAT('0 ***'/'  *** FOR MOD(ITYPE,10) = 3 ATTEMP TO DETERMINE'
     1      ,' IDENT FROM JLEV DATA'/'  ***')
      J1GEJ2=.TRUE.
      J1LEJ2=.TRUE.
      DO 1000 I=1,NLEV
      J1=JLEV(I)
      J2=JLEV(NLEV+I)
      J1GEJ2=J1GEJ2.AND.J1.GE.J2
 1000 J1LEJ2=J1LEJ2.AND.J1.LE.J2
      IF (J1LEJ2.AND.J1GEJ2) THEN
      WRITE(6,601)
  601 FORMAT('0 J1=J2 FOR ALL BASIS FUNCTIONS ',
     1       ' -- ASSUME NON-IDENTICAL')
      IDENT=0
      GO TO 9000
      ENDIF
      IF ((.NOT.J1LEJ2).AND.(.NOT.J1GEJ2)) THEN
      WRITE(6,602)
  602 FORMAT('0 BASIS NOT ORDERED J1.GT.J2 OR J1.LE.J2',
     1       ' -- ASSUME NON-IDENTICAL')
      IDENT=0
      GO TO 9000
      ENDIF
C     IF WE REACH BELOW, EITHER J1.GE.J2 **OR** J1.LE.J2 FOR ALL BASIS
      WRITE(6,603)
  603 FORMAT('0 BASIS IS ORDERED EITHER J1.GT.J2 OR J1.LE.J2',
     1       ' -- ASSUME IDENTICAL PARTICLES')
      IDENT=1
 9000 WRITE(6,604) IDENT
  604 FORMAT('0 *** IDENT HAS BEEN SET TO',I3)
      RETURN
      END
      SUBROUTINE SWRITE(IU,N,S)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION S(N,N)
C
      WRITE(IU) ((S(I,J),J=1,I),I=1,N)
      RETURN
C
      ENTRY SREAD(IU,N,S,IEND)
      IEND=0
      READ(IU,END=9999) ((S(I,J),J=1,I),I=1,N)
      DO 1000 I=1,N
      DO 1000 J=1,I-1
 1000 S(J,I)=S(I,J)
      RETURN
 9999 IEND=1
      RETURN
      END