C     PROGRAM TO COMPUTE  DIFFERENTIAL (AND TOTAL) CROSS SECTIONS
C       FROM S-MATRICES SAVED ON TAPE BY MOLSCAT SCATTERING PROG.
C     COMBINED CLOSE COUPLING / COUPLED STATES VERSION
C       JAN 95 (SG) COMPATIBLE WITH MOLSCAT VERSION 14
C     LIMITED TO ITYP=1,2,7  (PAINLESS UPGRADE FOR ITYP=5,6 ?)
C     CALCULATES ONLY A SINGLE SAVED ENERGY, ENERGY(ICALC)
C
C     NAMELIST (&DIFF) INPUT VARIABLES:
C     --------------------------------
C       JIMIN,JIMAX,JSTEP - CALC FOR JI FROM JIMIN TO JIMAX BY JSTEP
C       JFMAX - GET CROSS SECTIONS FOR JI->JF JF=JI,JFMAX,JSTEP
C       ANGMIN,ANGMX,NANG - DO CALC AT NANG ANGLES ON (ANGMIN,ANGMAX)
C       JTOTU - WILL LIMIT CALCULATION TO JTOT.LT.JTOTU
C       ICALC - DO CALCULATION FOR ENERGY(ICALC)
C       IT    - INPUT TAPE UNIT (DEFAULT 10)
C       PRNTLV - .EQ.2 ECHO INPUT USED; .GE.5 ECHO ALL INPUT
C                .GE.3 PRINT SIG(I->F); .GE.15 PRINT S-MATRICES
C       LSIN   - DEFAULT IS .FALSE. TO PRINT DSIG/DTH
C                SET TO     .TRUE. TO PRINT LN (SIN(TH)*DSIG/DTH)
C       LFMT   - LOGICAL VARIABLE FOR FORMATTED SAVE TAPE
C                DEFAULT IS .FALSE. (UNFORMATTED)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     MOLSCAT WORKING STORAGE; USED FOR SIGMA'S AND DCS
      PARAMETER (MXDIM=300000)
      DIMENSION X(MXDIM)
      COMMON /MEMORY/ MX,IXNEXT,NIPR,IVLFL,X
C
C     COMMON FOR NAMELIST VARIABLES
      LOGICAL LSIN
      COMMON/DFF/ANGMIN,ANGMAX,NANG,ICALC,JIMIN,JIMAX,JSTEP,JFMAX,LSIN
C
C     STORAGE FOR SAVE TAPE VALUES (SREAL/SIMAG ARE IN X())
      PARAMETER (MXCH=500,MXJLEV=500)
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)
C     N.B. WE ACCUMULATE SIGMA FOR ONLY ONE OF THE ENERGIES ...
      DIMENSION ENERGY(MXNRG),MINJT(1),MAXJT(1)
      DIMENSION ELEVEL(MXEL)
      INTEGER ITYPE,NLEV,NQN,NNRG,NOPEN,PRNTLV
      CHARACTER*4 LABEL(20)
      LOGICAL LFMT
C
      DATA JSTEPT/-1/
      DATA PI/3.14159 26535 89793 D0/
C
      NAMELIST /DIFF/JIMIN,JIMAX,JFMAX,JSTEP,NANG,ANGMIN,ANGMX,
     1               IT,ICALC,PRNTLV,JTOTU,LFMT,LSIN
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     INITIALIZE VALUES IN /MEMORY/
      MX=MXDIM
      IXNEXT=1
      NIPR=2
      CALL CHKSTR(-1)
C
C     DEFAULT VALUES FOR NAMELIST (&DIFF) INPUT
      LFMT=.FALSE.
      LSIN=.FALSE.
      JTOTU=999999
      IT=10
      ICALC=1
      JIMIN=0
      JIMAX=0
      JFMAX=0
      JSTEP=1
      NANG=181
      ANGMIN=0.D0
      ANGMAX=180.D0
      PRNTLV=2
C     OTHER INITIALIZATIONS  ...
      DO 10 I=1,MXCH
   10 NBASIS(I)=I
      DO 11 I=1,MXEL
   11 ELEVEL(I)=0.D0
      RM=1.D0
C
C     GET (POSSIBLE) NAMELIST INPUT DATA
      READ(5,DIFF,END=8000)
 8000 WRITE(6,608) IT,ICALC,JIMIN,JIMAX,JFMAX,JSTEP,
     1             NANG,ANGMIN,ANGMAX,PRNTLV
  608 FORMAT('  *** PROGRAM TO COMPUTE STATE-TO-STATE DIFFERENTIAL',
     1       ' CROSS SECTIONS'/6X,'FROM S-MATRICES SAVED ON TAPE UNIT',
     2       I3/6X,'D(SIGMA)/D(OMEGA) FOR ENERGY NO.',I4/
     3       /'  *** DIFFERENTIAL SCATTERING AMPLITUDES REQUESTED FOR'/
     4        6X,'JIMIN =',I5,5X,'JIMAX =',I5,5X,'JFMAX =',I5,
     5        5X,'JSTEP =',I5/6X,I5,' ANGLES BETWEEEN',F8.2,' AND',
     6        F8.2,' DEGREES'//6X,'PRINT CONTROL, PRNTLV =',I3)
      IF (JFMAX.LT.JIMAX) THEN
        WRITE(6,615) JFMAX,JIMAX
  615   FORMAT(/'  JFMAX =',I4,'  INCREASED TO EQUAL JIMAX =',I4)
        JFMAX=JIMAX
      ENDIF
      IF (LSIN) THEN
        WRITE(6,*)
        WRITE(6,*) ' *** LSIN REQUESTS LN (SIN(TH) * DSIG/DTH) OUTPUT'
      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 W/ REQUESTED TAPE FORMAT'
        STOP
      ENDIF
      ITYP=ITYPE-10*(ITYPE/10)
      IF (.NOT.(ITYP.EQ.1.OR.ITYP.EQ.2.OR.
     1                       ITYP.EQ.7)) THEN
        WRITE(6,*) ' *** CODE NOT IMPLEMENTED FOR ITYPE',ITYPE
        STOP
      ENDIF
      IF (ITYP.EQ.2.OR.ITYP.EQ.7) WRITE(6,613)
  613 FORMAT(/'  ***** LIMITED IMPLEMENTATION FOR ITYP = 2,7',
     1       /'       CROSS SECTIONS IN V=0 ONLY.')
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
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 ON SAVE TAPE IS',I4/
     1       ('  ',1P,5D14.4))
      WRITE(6,607) ICALC,ENERGY(ICALC)
  607 FORMAT(/'  ***'/'  *** DIFFERENTIAL CROSS SECTIONS ONLY FOR',
     1       ' ENERGY(',I3,') =',F10.4,' (1/CM)'/'  ***')
C
C     INITIALIZE OUTINX (WHICH INITIALIZES DFFCC)
      IXSIG=IXNEXT
      CALL OUTINX(ENERGY(ICALC),1,NLEV,NQN,JLEV,X(IXSIG),ITYPE,
     1            MINJT,MAXJT)
C
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     START READ SAVE TAPE LOOP OVER JTOT / ENERGY VALUES
C
      JSTEPX=-1
      JTOLD=-1
      IXSM=IXNEXT-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,*)  ' *** ROLLOUT MARKER ENCOUNTERED',JTOT
        GO TO 2000
      ENDIF
      IF (JTOTU.LT.999999.AND.JTOT.GT.JTOTU) THEN
        WRITE(6,*) ' *** STOPPING INPUT ON JTOT.GT.JTOTU',JTOT,JTOTU
        GO TO 9000
      ENDIF
C     ECHO INPUT
      IF (PRNTLV.GE.5) 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
C *** SHOULD PROBABLY LOOK AT WT FOR REASONABLENESS
C
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('  *** 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(/'  *** 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
      IXNEXT=IXNEXT+2*NSQ
      NUSED=0
      CALL CHKSTR(NUSED)
      IF (LFMT) THEN
C       WE SHOULD ARRANGE TO "SKIP" IF INRG.NE.ICALC
        READ(IT,105,END=9090) (X(IXSM+I),I=1,NSQ)
        READ(IT,105,END=9090) (X(IXSM+NSQ+I),I=1,NSQ)
      ELSE
         CALL SREAD(IT,NOPEN,X(IXSM+1),IEND)
         IF (IEND.GT.0) GO TO 9090
         CALL SREAD(IT,NOPEN,X(IXSM+NSQ+1),IEND)
         IF (IEND.GT.0) GO TO 9090
      ENDIF
C
C     CALL OUTPUT, BUT ONLY FOR INRG=ICALC
      IF (INRG.EQ.ICALC) THEN
        IF (PRNTLV.EQ.2) WRITE(6,663) IT,JTOT,M,INRG,ECHK,IEXCH,WT,NOPEN
        CALL OUTPUX(JTOT,NBASIS,J,L,WVEC,X(IXSM+1),X(IXSM+NSQ+1),
     1        NOPEN,M,WT,1,RM,PRNTLV,
     2        ENERGY(INRG),X(IXSIG),JLEV,MINJT,MAXJT,NLEV,NQN)
      ENDIF
C     GO BACK FOR MORE JTOT, INRG SETS . . .
      IXNEXT=IXSM+1
      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.
 9123 CALL OUTPCX(X(IXSIG),ENERGY(ICALC),1,MINJT,MAXJT,
     1            PRNTLV,LABEL,JSTEP)
      STOP
      END
      SUBROUTINE OUTPUX(JTOT,NBASIS,J,L,WVEC,SREAL,SIMAG,
     1        NOPEN,M,WT,INRG,RM,PRINT,
     2        ENERGY,SIG,JLEV,MINJT,MAXJT,NLEV,NQN)
C     THIS IS A MODIFICATION OF (VERSION 12) OUTPUT ROUTINE FOR
C      (DCS_SAVE.F) DIFFERENTIAL SCTATTERING UTILITY IN VERSION 14.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE IXDIFF,NSTOR,NLEVEL,LCS
C
C     ENTRY OUTPUX - PROCESSES S MATRICES TO X-SECTIONS OUTPUTS THEM.
C     ENTRY OUTINX - INITIALIZATION ENTRY FROM MAIN PROGRAM
C     ENTRY OUTPCX - PRINTS FINAL CROSS SECTIONS.
C
      DIMENSION NBASIS(1),J(1),L(1),WVEC(1),SREAL(1),SIMAG(1)
      INTEGER MAXJT(1),MINJT(1)
      INTEGER PRINT
      INTEGER NLEV,NQN,JLEV(1)
      DIMENSION SIG(1),ENERGY(1)
      LOGICAL LCS
      CHARACTER*4 LABEL(20)
      CHARACTER*1 BLANK
C
C     NEED IXNEXT FROM /MEMORY/
      COMMON/MEMORY/MX,IXNEXT,NIPR,IVLFL,X(1)
C
      DATA BLANK/' '/
      DATA EPS/1.D-10/
      DATA PI/3.14159 26535 89793 D0/
C
C     THIS PROGRAM FOR DFF_SAVE REQUIRES THAT INRG.EQ.1
      IF (INRG.NE.1) THEN
        WRITE(6,*) ' *** OUTPUX. ILLEGAL INRG VALUE',INRG
        STOP
      ENDIF
C     BOOKEEPING FOR MINJT AND MAXJT
      IF (MINJT(INRG).LT.0)  MINJT(INRG)=JTOT
      IF (JTOT.GT.MAXJT(INRG)) MAXJT(INRG)=JTOT
C
C     PRINT OUT OPEN-CHANNEL BASIS FUNCTIONS AND WAVEVECTORS IN 1/A.
      IF (PRINT.GE.15)  WRITE(6,601)
  601 FORMAT('0 CHANNEL NO.  TARGET LEVEL  ORBITAL L   WVEC(1/ANG.)')
      IF (ABS(RM-1.D0).GT.1.D-8) THEN
        WRITE(6,*) ' *** OUTPUX. RM.NE.1.D0',RM
        STOP
      ENDIF
C
C     PROCESS S-MATRIX. ACCUMULATE X-SECTIONS IN TSIG. PRINT.
C     J(NBASIS(I)) IS LEVEL NUMBER OF ITH BASIS FN. IN ASYMPTOTIC AREA.
C  CLEAR TSIG.
      DO 1400 I=1,NSTOR
 1400 SIG(I)=0.D0
      IF (PRINT.GE.15) WRITE(6,606)
  606 FORMAT('0 ROW  COL',10X,'S**2',15X,'PHASE/2PI',12X,'RE (S)',14X,
     1      'IM (S)'  )
C
      IJ=0
      NTOP=(NQN-1)*NLEV
C     CALCULATE GLOBAL MULTIPLICATIVE FACTOR FOR X-SECTIONS.
      XJ=DBLE(2*JTOT+1)*PI
      PI2=2.D0*PI
      IF (WT.GT.0.D0) XJ=XJ*WT
      DO 2000 ICOL=1,NOPEN
      LEVC=J(NBASIS(ICOL))
      LCOL=ABS(JLEV(NTOP+LEVC))
      DO 2000 IROW=1,NOPEN
      DD=WVEC(NBASIS(IROW))
      DD=DD*DD
      LEVR=J(NBASIS(IROW))
C     IJ IS INDEX OF SREAL,SIMAG(IROW,ICOL)
      IJ=IJ+1
      SMAG=( SREAL(IJ)*SREAL(IJ)+SIMAG(IJ)*SIMAG(IJ) )
      IF (PRINT.LT.15 .OR. SMAG.LE.EPS)  GO TO 2300
      PHASE=ATAN2(SIMAG(IJ),SREAL(IJ)) / PI2
 2100 WRITE(6,607) IROW,ICOL,SMAG,PHASE,SREAL(IJ),SIMAG(IJ)
  607 FORMAT(2I5,4E20.6)
 2300 IF (IROW.NE.ICOL)  GO TO 2400
C     FOR IROW = ICOL,  CALCULATE  T = 1 - S.
      SMAG=1.D0-SREAL(IJ)
      SMAG=SMAG*SMAG + SIMAG(IJ)*SIMAG(IJ)
 2400 CONTINUE
      LROW=ABS(JLEV(NTOP+LEVR))
      IF (LROW.GT.NLEVEL .OR. LCOL.GT.NLEVEL) GO TO 2000
C     II IS INDEX OF SIG(ICOL,IROW).  N.B. JLEV(LEV,NQN) HAS POINTER
C          TO 'SERIAL' NUMBER OF 'LEVEL'.
      II=(LROW-1)*NLEVEL+LCOL
C     ACCOUNT FOR K(J,J), DEGEN. LATTER IN SIG(NSTOR+II).
      SIG(II) = SIG(II) + SMAG*XJ/(SIG(NSTOR+II)*DD)
 2000 CONTINUE
C
C     ACCUMULATE X-SECTIONS.  SET IJ TO START OF INRG-TH ENERGY IN SIG
 4100 IJ=(INRG+1)*NSTOR
      II=0
      DO 3000 JI=1,NLEVEL
      DO 3000 I=1,NLEVEL
      II=II+1
      IJ=IJ+1
      SIG(IJ)=SIG(IJ)+SIG(II)
 3000 CONTINUE
C
C     PROCESS S-MATRIX TO DIFFERENTIAL SCATTERING CROSS SECTIONS
      IF (LCS) THEN
        CALL DFFCS(JTOT,M,WT,NBASIS,J,L,WVEC,NOPEN,SREAL,SIMAG,
     1            JLEV,NLEV,NQN,
     2            INRG,ENERGY,SIG(IXDIFF),SIG(IXDIFF))
      ELSE
        CALL DFFCC(JTOT,M,WT,NBASIS,J,L,WVEC,NOPEN,SREAL,SIMAG,
     1            JLEV,NLEV,NQN,
     2            INRG,ENERGY,SIG(IXDIFF),SIG(IXDIFF))
      ENDIF
C
      RETURN
C
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     INITIALIZATION ENTRY.
C
      ENTRY OUTINX(ENERGY,NNRG,NLEV,NQN,JLEV,SIG,ITYP,MINJT,MAXJT)
C
      IF (NNRG.NE.1) THEN
        WRITE(6,*) ' *** OUTINX. ILLEGAL NNRG VALUE',INRG
        STOP
      ENDIF
      ITYPE=ITYP
      LCS=ITYP.GE.21.AND.ITYP.LE.29
C     FIND NLEVEL AS MAX. INDEX FROM JLEV(I,NQN) . . .
      IJ=(NQN-1)*NLEV
      NLEVEL=0
      DO 6196 II=1,NLEV
      IJ=IJ+1
 6196 NLEVEL=MAX0(NLEVEL,ABS(JLEV(IJ)))
      IF (NLEVEL.LE.0)  THEN
        WRITE(6,603)
  603   FORMAT('0 * * * ERROR.  NO LEVELS IN SIG. MATRICES.')
        STOP
      ENDIF
      NSTOR=NLEVEL*NLEVEL
C
C     CALCULATE REQUIRED STORAGE AND INITIALIZE IT.  NSTOR LOCATIONS
C     FOR TSIG; NSTOR FOR DEGEN; NNRG*NSTOR FOR SIG.
      IJ=NSTOR*(NNRG+2)
C     IXNEXT INCREMENTED TO REFLECT ADDITIONAL STORAGE TAKEN UP BY SIG.
      IXNEXT=IXNEXT+IJ
C     N.B NUSED MUST BE INITIALIZED NON-NEGATIVE BEFORE CHKSTR
      NUSED=0
      CALL CHKSTR(NUSED)
C     SET POINTER FOR DIFFERENTIAL XSECT (DFFCC) STORAGE AFTER SIG()
      IXDIFF=IJ+1
C     INITIALIZE SIG TO ZERO
      II=2*NSTOR+1
      DO 6200 I=II,IJ
 6200 SIG(I)=0.D0
C
C     DEGENERACY INFORMATION IN SIG(NSTOR+(LROW-1)*NLEVEL+LCOL)
      CALL DEGENG(ITYPE,JLEV,NLEV,NQN,NLEVEL,SIG(NSTOR+1))
C
C     SET UP 'BOOKKEEPING' VARIABLES.
 6321 DO 6100 I=1,NNRG
      MINJT(I)=-1
 6100 MAXJT(I)=0
C
C     INITIALIZE DFFCC/DFFCS -- N.B STORAGE FOR BOTH NTB AND DFX START
C                          AT SIG(IXDIFF).
      IF (LCS) THEN
        CALL DFFCSI(ENERGY,NNRG,JLEV,NLEV,NQN,ITYP,
     1            SIG(IXDIFF),SIG(IXDIFF))
      ELSE
        CALL DFFCCI(ENERGY,NNRG,JLEV,NLEV,NQN,ITYP,
     1            SIG(IXDIFF),SIG(IXDIFF))
      ENDIF
      RETURN
C
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
      ENTRY OUTPCX(SIG,ENERGY,NNRG,MINJT,MAXJT,PRINT,LABEL,JSTEP)
      IF(PRINT.LT.3) GO TO 9000
      WRITE(6,9678) LABEL
 9678 FORMAT(//20A4)
      IF(JSTEP.GT.1) WRITE(6,9604) JSTEP
 9604 FORMAT(/'  ***  CROSS SECTIONS MULTIPLIED BY JSTEP =',I3)
      WRITE(6,9679)
 9679 FORMAT(/12X,'ENERGY  JTOTL  JSTEP  JTOTU',7X,
     1       'F    I',12X,'SIG(F,I)'/)
 9001 IJ=2*NSTOR
      XJSTEP=JSTEP
      DO 9900 K=1,NNRG
      MN=MINJT(K)
      MXJ=MAXJT(K)
      EK=ENERGY(K)
      DO 9900 I=1,NLEVEL
      DO 9900 II=1,NLEVEL
      IJ=IJ+1
      SIJ=SIG(IJ)*XJSTEP
      IF (SIJ.GE.EPS)
     1WRITE(6,103) BLANK,EK,MN,JSTEP,MXJ,II,I,SIJ
  103 FORMAT(A1,F19.6,I5,2I7, 5X,2I5,E20.6)
 9900 CONTINUE
C
 9000 IF (LCS) THEN
        CALL DFFCSO(SIG(IXDIFF),SIG(IXDIFF),ENERGY,NNRG,JSTEP)
      ELSE
        CALL DFFCCO(SIG(IXDIFF),SIG(IXDIFF),ENERGY,NNRG,JSTEP)
      ENDIF
      RETURN
      END
      SUBROUTINE DEGENG(ITYPE,JLEV,NLEV,NQN,NLEVEL,DG)
C     FOR ITYP=1,2,7 (USES FACT THAT JLEV(I) HAS J-VALUE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION JLEV(NLEV,NQN),DG(NLEVEL,NLEVEL)
      DO 1000 I1=1,NLEVEL
      DO 1000 I2=1,NLEVEL
 1000 DG(I1,I2)=-9.9
      DO 2000 I1=1,NLEV
      IX=ABS(JLEV(I1,NQN))
      IF (IX.LT.0.OR.IX.GT.NLEVEL) THEN
        WRITE(6,*) ' *** DEGENG. INDEX.GT.NLEVEL',IX,NLEVEL
        STOP
      ENDIF
      JI=JLEV(I1,1)
      XF=JI+JI+1
      DO 2100 I2=1,NLEVEL
 2100 DG(I2,IX)=XF
 2000 CONTINUE
      DO 3000 I1=1,NLEVEL
      DO 3000 I2=1,NLEVEL
      IF (DG(I1,I2).GT.0.D0) GO TO 3000
      WRITE(6,*) ' *** DEGENG. UNINITIALIZED DG',I1,I2
      STOP
 3000 CONTINUE
      RETURN
      END
      SUBROUTINE DFFCC(JTOT,MV,WT,NB,J,L,WV,NOP,SREAL,SIMAG,
     1                  JLEV,NLEV,NQN,
     2                  INRG,ENERGY,NTB,DFX)
C
C     DIFFERENTIAL SCATTERING CODE.  JAN 95 VERSION FOR SAVE TAPE.
C     LIMITED IMPLEMENTATION FOR ITYPE=1,2,7 (CLOSE COUPLING) ONLY
C       (NB ITYPE=5,6 SHOULD BE PAINLESS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      SAVE NDF,IXCOS,IXPL,SQRTHF
C
C     NCOL IS NUMBER OF COLS FOR OUTPUT OF DSIG/DTHETA; LEGAL = 4 OR 6
      DIMENSION NB(1),J(1),L(1),WV(1),SREAL(1),SIMAG(1)
      DIMENSION ENERGY(1),JLEV(NLEV,NQN)
      DIMENSION NTB(5,1),DFX(1)
      LOGICAL FIRST,LAST,ACCUM
C
C     NEED IXNEXT,NIPR FROM /MEMORY/
      COMMON/MEMORY/MX,IXNEXT,NIPR,IVLFL,X(1)
C
C     COMMON FOR NAMELIST VARIABLES
      LOGICAL LSIN
      COMMON/DFF/ANGMIN,ANGMAX,NANG,ICALC,JIMIN,JIMAX,JSTEP,JFMAX,LSIN
C
      DATA PI/3.14159 26535 89793 D0/
C
C     STATEMENT FUNCTIONS ...
      ANG(I)=ANGMIN+(ANGMAX-ANGMIN)*(I-1)/(NANG-1)
C
C     VERSION ONLY CALCULATES DIFFERENTIAL X-SECTION FOR SINGLE E (1)
      IF (INRG.NE.1) THEN
        WRITE(6,*) ' *** DFFCC CALLED WITH INRG.NE.1'
        STOP
      ENDIF
C
C     ARTHURS AND DALGARNO, EQ. 13 APPEARS TO GIVE 2*PI TIMES THE
C       DIFFERENTIAL CROSS SECTIONS WE WANT (INTEGRATED OVER PHI?)
C       HENCE, INSTEAD OF FACTOR SQRT(PI), USE SQRT(1/2).
      XJTOT=JTOT
      XJFACT=(XJTOT+XJTOT+1.D0)*SQRTHF
C
C     LOOP OVER NTB ENTRIES
      DO 2200 I=1,NDF
      JI=JLEV(NTB(1,I),1)
      XJI=JI
      JF=JLEV(NTB(2,I),1)
      XJF=JF
      MI=NTB(3,I)
      XMI=MI
      MF=NTB(4,I)
      XMF=MF
      MDIF=MI-MF
C     MAY NEED TO CHANGE SIGN OF PLM IF MDIF.LT.0
      SGN=1.D0
      IF (MDIF.LT.0.AND.MDIF-2*(MDIF/2).NE.0) SGN=-SGN
      XMDIF=XMI-XMF
      IXR=NTB(5,I)
      IXI=IXR+NANG
C     NOW LOOP OVER S-MATRIX ELEMENTS TO SEE IF THEY CONTRIBUTE
      DO 2100 IR=1,NOP
      IXS=IR
      IF (J(IR).NE.NTB(2,I)) GO TO 2100
      LF=L(IR)
      IF (LF.LT.ABS(MDIF)) GO TO 2100
      XLF=LF
C     GET 3-J SYMBOL
      TJ2=THRJ(XJF,XLF,XJTOT,XMF,XMDIF,-XMI)
C     GET ASSOCIATED LEGENDRE POLYNOMIALS
      DO 2105 IA=1,NANG
 2105 DFX(IXPL+IA)=SGN*PLM(LF,MDIF,DFX(IXCOS+IA))
      DO 2101 IC=1,NOP
      IF (J(IC).NE.NTB(1,I)) GO TO 2101
      LI=L(IC)
      XLI=LI
      XLFACT=SQRT(XLI+XLI+1.D0)
      WVEC=WV(NB(IC))
C     GET 3-J
      TJ1=THRJ(XJI,XLI,XJTOT,XMI,0.D0,-XMI)
C     GET DELTA FOR T=DELTA-S
      IF (IR.EQ.IC) THEN
        DELTA=1.D0
      ELSE
        DELTA=0.D0
      ENDIF
      LDIF=LI-LF
 2199 IF (LDIF.LT.0) THEN
        LDIF=LDIF+400
        GO TO 2199
      ENDIF
      LDIF=LDIF-4*(LDIF/4)
      IF (LDIF.EQ.0) THEN
        XTR=DELTA-SREAL(IXS)
        XTI=-SIMAG(IXS)
      ELSEIF (LDIF.EQ.1) THEN
        XTR=SIMAG(IXS)
        XTI=DELTA-SREAL(IXS)
      ELSEIF (LDIF.EQ.2) THEN
        XTR=-DELTA+SREAL(IXS)
        XTI=SIMAG(IXS)
      ELSEIF (LDIF.EQ.3) THEN
        XTR=-SIMAG(IXS)
        XTI=-DELTA+SREAL(IXS)
      ELSE
        WRITE(6,*) ' *** DFFCC.  ILLEGAL LOGICAL PATH.'
        STOP
      ENDIF
      FACT=XJFACT*XLFACT*TJ1*TJ2/WVEC
      DO 2150 IA=1,NANG
      DFX(IXR+IA)=DFX(IXR+IA)+XTR*DFX(IXPL+IA)*FACT
 2150 DFX(IXI+IA)=DFX(IXI+IA)+XTI*DFX(IXPL+IA)*FACT
 2101 IXS=IXS+NOP
 2100 CONTINUE
 2200 CONTINUE
      RETURN
 
C     -------------------------------------------------------
      ENTRY DFFCCI(ENERGY,NNRG,JLEV,NLEV,NQN,ITYP,NTB,DFX)
      WRITE(6,600)
  600 FORMAT(/' DIFFERENTIAL SCATTERING CROSS SECTION (JAN 95).  ',
     1        'INITIALIZATION.')
C
C     SET-UP TABLE FOR DIFFERENTIAL SCAT AMPLITUDES; NDF COUNTS THESE
C     N.B. WE BEGIN FILLING THIS TABLE WITHOUT CHECKING SPACE IN STORAGE
      NDF=0
      DO 1000 JI=JIMIN,JIMAX,JSTEP
        DO 1001 I=1,NLEV
        IF (JLEV(I,1).NE.JI) GO TO 1001
        IF ((ITYP.EQ.2.OR.ITYP.EQ.7).AND.JLEV(I,2).NE.0) GO TO 1001
        LEVI=I
C       N.B. LEVI=JLEV(I,NQN); SHOULD EQUAL I FOR SUPPORTED ITYPE
        GO TO 1002
 1001   CONTINUE
      WRITE(6,604) JI
  604 FORMAT(' *** WARNING.  REQUESTED LEVEL NOT IN BASIS.  JI =',I4)
      GO TO 1000
 1002 DO 1009 JF=JI,JFMAX,JSTEP
        DO 1003 I=1,NLEV
        IF (JLEV(I,1).NE.JF) GO TO 1003
        IF ((ITYP.EQ.2.OR.ITYP.EQ.7).AND.JLEV(I,2).NE.0) GO TO 1003
        LEVF=I
        GO TO 1004
 1003   CONTINUE
      WRITE(6,605) JF
  605 FORMAT(' *** WARNING.  REQUESTED LEVEL NOT IN BASIS.  JF =',I4)
      GO TO 1009
 1004 DO 1008 MVAL=-JI,JI
      DO 1008 MVAL1=0,JF
      NDF=NDF+1
      NTB(1,NDF)=LEVI
      NTB(2,NDF)=LEVF
      NTB(3,NDF)=MVAL
 1008 NTB(4,NDF)=MVAL1
 1009 CONTINUE
 1000 CONTINUE
      IF (NDF.LE.0) THEN
        WRITE(6,606)
  606   FORMAT('  DFFCC.  ***** ERROR.  NO LEVELS -- TERMINATING')
        STOP
      ENDIF
      NSTBL=(5*NDF+NIPR-1)/NIPR
      IXCOS=NSTBL
      IXPL=IXCOS+NANG
      IX=IXPL+NANG
      IXMN=IX+1
C     REAL AMPLITUDE AT ANGLE I IS IN  DFX(NTB(I,5)+I)
C        IMAGINARY AMPLITUDE    IS IN  DFX(NTB(I,5)+NANG+I)
      DO 1007 I=1,NDF
      NTB(5,I)=IX
 1007 IX=IX+2*NANG
      WRITE(6,607) NDF,IX
  607 FORMAT(' DFFCC.  TABLES SET FOR',I5,
     1       ' DIFFERENTIAL SCATTERING AMPLITUDES.'/
     2       10X,'REQUIRED STORAGE IS',I10//)
      IXMX=IX
C     INCREMENT IXNEXT TO REFLECT STORAGE USED
      IXNEXT=IXNEXT+IX
      NUSED=0
      CALL CHKSTR(NUSED)
C
C     SET-UP COS(THETA-I) IN ALLOCATED STORAGE, AND ZERO AMPLITUDE STOR
C       WE DO NOT SAVE ANGLES; RECALC THEM FROM ANGMIN,ANGMAX,NANG
      DO 1010 I=1,NANG
 1010 DFX(IXCOS+I)=COS(PI*ANG(I)/180.D0)
      DO 1020 I=IXMN,IXMX
 1020 DFX(I)=0.D0
      SQRTHF=SQRT(5.D-1)
      RETURN
C
C     -------------------------------------------------------
      ENTRY DFFCCO(NTB,DFX,ENERGY,NNRG,JTSTEP)
      IF (NDF.LE.0) RETURN
      IF (JTSTEP.NE.1) WRITE(6,610) JTSTEP
  610 FORMAT(' DFFCC. ***** WARNING'/8X,'DIFFERENTIAL CROSS SECTIONS',
     1       ' MAY ERR FOR JTOT STEP =',I4)
C     INITIALIZE LEVELS TO ZERO, AS REAL LEVELS ARE POSITIVE VALUES
      LEVI=0
      LEVF=0
      FIRST=.TRUE.
      LAST=.FALSE.
      I=1
 3000 IF (I.GT.NDF) THEN
        LAST=.TRUE.
        GO TO 3500
      ENDIF
      NEWI=NTB(1,I)
      NEWF=NTB(2,I)
      IF (NEWI.NE.LEVI .OR. NEWF.NE.LEVF) GO TO 3500
C
C     ACCUMULATE DIFFERENTIAL X-SECTIONS FOR THIS LEVI,LEVF
C       PUT DSIG/DTHETA INTO 'PLM STGORAGE' AS IT IS NO LONGER USED
 3300 IXR=NTB(5,I)
      IXI=NTB(5,I)+NANG
C     WEIGHT BY 1/(2*JI+1) AND, FOR NON-ZERO MF-VALUE, BY 2
C       JI FROM NTB() AS JLEV() IS NOT AVAILABLE HERE.
      WGHT=NTB(1,I)-1
      WGHT=1.D0/(2.D0*WGHT+1.D0)
      IF (NTB(4,I).NE.0) WGHT=2.D0*WGHT
      WRITE(6,625) NTB(1,I)-1,NTB(2,I)-1,NTB(3,I),NTB(4,I),WGHT
  625 FORMAT(' DFFCCO.  ACCUMULATING JI, JF =',2I3,', MI, MF =',
     1       2I3,', WEIGHTED BY',F9.5)
      DO 3310 K=1,NANG
      ADD=DFX(IXR+K)*DFX(IXR+K)+DFX(IXI+K)*DFX(IXI+K)
 3310 DFX(IXPL+K)=DFX(IXPL+K)+WGHT*ADD
C
      I=I+1
      GO TO 3000
C
C     NEW LEVI,LEVF.  IF NOT FIRST, OUTPUT BEFORE ACCUMULATING NEW ONE
 3500 IF (FIRST) GO TO 3600
C
C     PRINT DIFFERENTIAL CROSS SECTIONS AND TRAPEZOID RULE INTEGRAL
      CALL DCSPRT(DFX(IXPL+1),DFX(IXCOS+1))
C
      IF (LAST) GO TO 9000
C     ZERO STORAGE FOR NEXT SET
 3600 DO 3610 K=1,NANG
 3610 DFX(IXPL+K)=0.D0
      LEVI=NEWI
      LEVF=NEWF
      FIRST=.FALSE.
      WRITE(6,*)
      GO TO 3300
 9000 RETURN
      END
      SUBROUTINE DFFCS(JTOT,MV,WT,NB,J,L,WV,NOP,SREAL,SIMAG,
     1                  JLEV,NLEV,NQN,
     2                  INRG,ENERGY,NTB,DFX)
C
C     DIFFERENTIAL SCATTERING CODE.  27 OCT 94 VERSION FOR SAVE TAPE.
C     LIMITED IMPLEMENTATION FOR ITYPE=21,22,27 (COUPLED STATES) ONLY
C       (NB ITYPE=25,26 SHOULD BE PAINLESS)
C     MERGED WITH CC CODE JAN 95
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      SAVE NDF,IXCOS,IXPL
C
      DIMENSION NB(1),J(1),L(1),WV(1),SREAL(1),SIMAG(1)
      DIMENSION ENERGY(1),JLEV(NLEV,NQN)
      DIMENSION NTB(5,1),DFX(1)
      LOGICAL FIRST,LAST,ACCUM,NEQ,LFMT
C
C     NEED IXNEXT,NIPR FROM /MEMORY/
      COMMON/MEMORY/MX,IXNEXT,NIPR,IVLFL,X(1)
C
C     COMMON FOR NAMELIST VARIABLES
      COMMON/DFF/ANGMIN,ANGMAX,NANG,ICALC,JIMIN,JIMAX,JSTEP,JFMAX
C
      DATA PI/3.14159 26535 89793 D0/
C
C     STATEMENT FUNCTIONS ...
      ANG(I)=ANGMIN+(ANGMAX-ANGMIN)*(I-1)/(NANG-1)
      NEQ(A,B)=ABS(A-B).GT.1.D-8
C
C     VERSION ONLY CALCULATES DIFFERENTIAL X-SECTION FOR SINGLE E (1)
      IF (INRG.NE.1) THEN
        WRITE(6,*) ' *** DFFCS CALLED WITH INRG.NE.1'
        STOP
      ENDIF
C
C     MV IS 'SYMMETRY COUNTER'; ASSUME THAT MVALUE=MV-1
C       WHICH SHOULD BE CONSISTENT WITH MOLSCAT CODE
      MVALUE=MV-1
C     DO SOME CHECKING FOR CONSISTENCY OF MV AND WT
      IF (MV.EQ.1 .AND. NEQ(WT,1.D0)) WRITE(6,699) MV,WT
  699 FORMAT(' DFFCS ***** INCONSISTENT MV, WT',I4,F10.3)
      IF (MV.GT.1 .AND. NEQ(WT,2.D0)) WRITE(6,699) MV,WT
      IF (MV.LT.0) WRITE(6,699) MV,WT
C
C     GENERATE LEGENDRE FUNCTION AT COSINES OF ANGLES
C       ASSUME HERE THAT JTOT = L(NB(I)), I=1,NOPEN;  COULD BE CHECKED.
C     XLEG(I,Z)=SQRT(2.D0/(2*I+1))*PLM(I,0,Z), COMBINE W/(2*JTOT+1) FACT
      XLFACT=SQRT(2.D0*(JTOT+JTOT+1))
      DO 2001 I=1,NANG
 2001 DFX(IXPL+I)=XLFACT*PLM(JTOT,0,DFX(IXCOS+I))
C
C     LOOP OVER S-MATRIX AND SEE IF ELEMENTS CONTRIBUTE TO ACCUMULATED
C       DIFFERENTIAL CROSS SECTIONS.  N.B LEVEL IS ACTUALLY
C       JLEV(NB(ICOL),NQN), BUT WE USE  JLEV(I,NQN)=I  FOR ITYPE=21,22
      IJ=0
      DO 2100 ICOL=1,NOP
      LEVC=J(NB(ICOL))
      DO 2100 IROW=1,NOP
      LEVR=J(NB(IROW))
      IJ=IJ+1
      IF (IROW.EQ.ICOL) THEN
        DELTA=1.D0
      ELSE
        DELTA=0.D0
      ENDIF
      WVEC=.5D0/WV(NB(IROW))
      NUSE=0
        DO 2101 I=1,NDF
        IF (NTB(2,I).NE.MVALUE) GO TO 2101
        IF (NTB(1,I).NE.LEVR.OR.NTB(3,I).NE.LEVC) GO TO 2101
C
C       MATCH: ADD TO APPROPRIATE DFX().  CHECK FOR ONLY 1 MATCH.
C       USE EQ (8) FROM KOURI & MCGUIRE, CHEM PHYS LETT 29, 414 (1974)
C         THE AMPLITUDE HAS FACTOR 1/(2*K) IN FRONT OF T-MATRIX
C         (AND ALSO SQRT(-1), WHICH CAN BE IGNORED).
C       DIFFERS BY FACTOR OF 2 FROM MCGUIRE & KOURI,JCP 60, 2488 (1974)
C       FORMER GIVES PROPER INTEGRATED VALUE AS 2*PI*
C                          INTEGRAL(-1,1) OF D(COS(TH)) SIG(COS(TH))
        NUSE=NUSE+1
        NTB(5,I)=NTB(5,I)+1
        IXR=NTB(4,I)
        IXI=NTB(4,I)+NANG
        ADDR=WVEC*(DELTA-SREAL(IJ))
        ADDI=-WVEC*SIMAG(IJ)
        DO 2200 K=1,NANG
          DFX(IXR+K)=DFX(IXR+K)+ADDR*DFX(IXPL+K)
 2200     DFX(IXI+K)=DFX(IXI+K)+ADDI*DFX(IXPL+K)
 2101   CONTINUE
      IF (NUSE.GT.1) WRITE(6,698) IROW,ICOL
  698 FORMAT(' DFFCS.  ERROR.  MORE THAN ONE MATCH FOR S-MATRIX',2I4)
 2100 CONTINUE
      RETURN
 
C     -------------------------------------------------------
      ENTRY DFFCSI(ENERGY,NNRG,JLEV,NLEV,NQN,ITYP,NTB,DFX)
      WRITE(6,600)
  600 FORMAT(/' DIFFERENTIAL SCATTERING CROSS SECTION (APR 93).  ',
     1        'INITIALIZATION.')
C
C     SET-UP TABLE FOR DIFFERENTIAL SCAT AMPLITUDES; NDF COUNTS THESE
C     N.B. WE BEGIN FILLING THIS TABLE WITHOUT CHECKING SPACE IN STORAGE
      NDF=0
      DO 1000 JI=JIMIN,JIMAX,JSTEP
        DO 1001 I=1,NLEV
        IF (JLEV(I,1).NE.JI) GO TO 1001
        IF ((ITYP.EQ.22.OR.ITYP.EQ.27).AND.JLEV(I,2).NE.0) GO TO 1001
        LEVI=I
C       N.B. LEVI=JLEV(I,NQN); SHOULD EQUAL I FOR SUPPORTED ITYPE
        GO TO 1002
 1001   CONTINUE
      WRITE(6,604) JI
  604 FORMAT(' *** WARNING.  REQUESTED LEVEL NOT IN BASIS.  JI =',I4)
      GO TO 1000
 1002 DO 1009 JF=JI,JFMAX,JSTEP
        DO 1003 I=1,NLEV
        IF (JLEV(I,1).NE.JF) GO TO 1003
        IF ((ITYP.EQ.22.OR.ITYP.EQ.27).AND.JLEV(I,2).NE.0) GO TO 1003
        LEVF=I
        GO TO 1004
 1003   CONTINUE
      WRITE(6,605) JF
  605 FORMAT(' *** WARNING.  REQUESTED LEVEL NOT IN BASIS.  JF =',I4)
      GO TO 1009
 1004 DO 1008 MVAL=0,JI
      NDF=NDF+1
      NTB(1,NDF)=LEVI
      NTB(2,NDF)=MVAL
 1008 NTB(3,NDF)=LEVF
 1009 CONTINUE
 1000 CONTINUE
      IF (NDF.LE.0) THEN
        WRITE(6,606)
  606   FORMAT('  DFFCS.  ***** ERROR.  NO LEVELS -- TERMINATING')
        STOP
      ENDIF
      NSTBL=(5*NDF+NIPR-1)/NIPR
      IXCOS=NSTBL
      IXPL=IXCOS+NANG
      IX=IXPL+NANG
      IXMN=IX+1
C     REAL AMPLITUDE AT ANGLE I IS IN  DFX(NTB(I,4)+I)
C        IMAGINARY AMPLITUDE    IS IN  DFX(NTB(I,4)+NANG+I)
      DO 1007 I=1,NDF
      NTB(4,I)=IX
      NTB(5,I)=0
 1007 IX=IX+2*NANG
      WRITE(6,607) NDF,IX
  607 FORMAT(' DFFCS.  TABLES SET FOR',I5,
     1       ' DIFFERENTIAL SCATTERING AMPLITUDES.'/
     2       10X,'REQUIRED STORAGE IS',I10)
C     ---------------  DEBUGGING  -----------------------
C     WRITE(6,609) (I,(NTB(K,I),K=1,4),I=1,NDF)
C 609 FORMAT(' INDEX LEV-I  MVAL LEV-F        ADDR'/(4I6,I12))
C     ---------------  DEBUGGING  -----------------------
      IXMX=IX
C     INCREMENT IXNEXT TO REFLECT STORAGE USED
      IXNEXT=IXNEXT+IX
      NUSED=0
      CALL CHKSTR(NUSED)
C
C     SET-UP COS(THETA-I) IN ALLOCATED STORAGE, AND ZERO AMPLITUDE STOR
C       WE DO NOT SAVE ANGLES; RECALC THEM FROM ANGMIN,ANGMAX,NANG
      DO 1010 I=1,NANG
 1010 DFX(IXCOS+I)=COS(PI*ANG(I)/180.D0)
      DO 1020 I=IXMN,IXMX
 1020 DFX(I)=0.D0
      RETURN
C
      ENTRY DFFCSO(NTB,DFX,ENERGY,NNRG,JTSTEP)
      IF (NDF.LE.0) RETURN
      IF (JTSTEP.NE.1) WRITE(6,610) JTSTEP
  610 FORMAT(' DFFCS. ***** WARNING'/9X,'DIFFERENTIAL CROSS SECTIONS',
     1       ' MAY ERR FOR JTOT STEP =',I4)
C     INITIALIZE LEVELS TO ZERO, AS REAL LEVELS ARE POSITIVE VALUES
      LEVI=0
      LEVF=0
      FIRST=.TRUE.
      LAST=.FALSE.
      I=1
 3000 IF (I.GT.NDF) THEN
        LAST=.TRUE.
        GO TO 3500
      ENDIF
      NEWI=NTB(1,I)
      NEWF=NTB(3,I)
      IF (NEWI.NE.LEVI .OR. NEWF.NE.LEVF) GO TO 3500
C
C     ACCUMULATE DIFFERENTIAL X-SECTIONS FOR THIS LEVI,LEVF
 3300 IXR=NTB(4,I)
      IXI=NTB(4,I)+NANG
      ACCUM=NTB(5,I).GT.0
C     WEIGHT BY 1/(2*JI+1) AND, FOR NON-ZERO M-VALUE, BY 2
C       JI FROM NTB() AS JLEV() IS NOT AVAILABLE HERE.
      WGHT=NTB(1,I)-1
      WGHT=1.D0/(2.D0*WGHT+1.D0)
      IF (NTB(2,I).NE.0) WGHT=2.D0*WGHT
      WRITE(6,625) NTB(1,I)-1,NTB(3,I)-1,NTB(2,I),WGHT
  625 FORMAT(' DFFCSO.  ACCUMULATING JI, JF =',2I3,', M =',
     1       I3,', WEIGHTED BY',F9.5)
      ADDMAX=0.D0
      DO 3310 K=1,NANG
      ADD=DFX(IXR+K)*DFX(IXR+K)+DFX(IXI+K)*DFX(IXI+K)
      ADDMAX=MAX(ADDMAX,ADD)
      IF (ACCUM) DFX(IXPL+K)=DFX(IXPL+K)+WGHT*ADD
 3310 CONTINUE
      IF (.NOT.ACCUM) THEN
        WRITE(6,627)
  627   FORMAT(' ****  WARNING.  REQUESTED M-VALUE NOT AVAILABLE ****')
C       BELOW IS A RELIC FROM DEBUGGING ---------
        IF (NEQ(ADDMAX,0.D0)) WRITE(6,628) ADDMAX
  628   FORMAT(' ****  BUT MEMORY HAS BEEN CORRUPTED.  MAX(ADD) =',
     1          1PD12.4)
      ENDIF
C
      I=I+1
      GO TO 3000
C
C     NEW LEVI,LEVF.  IF NOT FIRST, OUTPUT BEFORE ACCUMULATING NEW ONE
 3500 IF (FIRST) GO TO 3600
C
C     PRINT DIFFERENTIAL CROSS SECTIONS AND TRAPEZOID RULE INTEGRAL
      CALL DCSPRT(DFX(IXPL+1),DFX(IXCOS+1))
C
      IF (LAST) GO TO 9000
C     ZERO STORAGE FOR NEXT SET
 3600 DO 3610 K=1,NANG
 3610 DFX(IXPL+K)=0.D0
      LEVI=NEWI
      LEVF=NEWF
      FIRST=.FALSE.
      GO TO 3300
 9000 RETURN
      END
      SUBROUTINE DCSPRT(DFX,COSANG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION DFX(NANG),COSANG(NANG)
      PARAMETER (NCOL=4)
      DIMENSION IPR(NCOL)
      LOGICAL COLPRT
C
C     COMMON FOR NAMELIST VARIABLES
      LOGICAL LSIN,LFMT
      COMMON/DFF/ANGMIN,ANGMAX,NANG,ICALC,JIMIN,JIMAX,JSTEP,JFMAX,LSIN
C
      DATA PI/3.14159 26535 89793 D0/
C     LFMT=.TRUE. PRINTS ANGLES DOWN EACH OUTPUT COLUMN
      DATA LFMT/.TRUE./
C
C     STATEMENT FUNCTIONS ...
      ANG(I)=ANGMIN+(ANGMAX-ANGMIN)*(I-1)/(NANG-1)
C
C     CHECK BY TRAPEZOIDAL RULE INTEGRAL; NOTE COS IS STORED (+1,-1)
      XINT=0.D0
      DO 3510 K=2,NANG
 3510 XINT=XINT+0.5D0*(DFX(K)+DFX(K-1))
     1               *(COSANG(K)-COSANG(K-1))
      XINT=-2.D0*PI*XINT
      IF (LSIN) THEN
        DO 3411 K=1,NANG
        SINDCS=SQRT(1.D0-COSANG(K)*COSANG(K))*DFX(K)
        IF (SINDCS.GT.0.D0) THEN
          DFX(IXPL+K)=LOG(SINDCS)
        ELSE
          DFX(IXPL+K)=0.D0
        ENDIF
 3411 CONTINUE
      WRITE(6,*)
      WRITE(6,*) ' ***** OUTPUT IS LOG(SIN(THETA) * DSIGMA/DTHETA)'
      ENDIF
C
C       PRINT HEADER FOR 4 OR 6 COLUMN OUTPUT
        IF (NCOL.EQ.4) THEN
          WRITE(6,631)
  631     FORMAT(/4('    ANG  DSIG/DTHETA'))
        ELSE
          WRITE(6,632)
  632     FORMAT(/6('    ANG  DSIG/DTHETA'))
        ENDIF
C
      IF (LFMT) THEN
C       PRINT DOWN COLUMNS IN ANGLE ORDER
        NLINE=(NANG+NCOL-1)/NCOL
        DO 3520 IL=1,NLINE
        IF (COLPRT(NANG,NCOL,IL,IPR,NPR))
     1      WRITE(6,630) (ANG(IPR(K)),DFX(IPR(K)),K=1,NPR)
  630     FORMAT(6(0P,F8.1,2X,1P,D10.3))
 3520   CONTINUE
      ELSE
C       LFMT=.FALSE. PRINT ACROSS ROWS IN ANGLE ORDER
        IF (NCOL.EQ.4) THEN
          WRITE(6,635) (ANG(K),DFX(K),K=1,NANG)
  635     FORMAT(4(0P,F8.1,2X,1P,D10.3))
        ELSE
          WRITE(6,636) (ANG(K),DFX(K),K=1,NANG)
  636     FORMAT(6(0P,F8.1,2X,1P,D10.3))
        ENDIF
      ENDIF
      WRITE(6,643) XINT
  643 FORMAT(10X,'INTEGRATED (TRAPEZOIDAL RULE) =',1P,D10.3/)
C
      RETURN
      END
      FUNCTION COLPRT(NVAL,NCOL,LINE,IX,NX)
C
C     TO FACILITATE PRINTING OF NVAL VALUES DOWN NCOL COLUMNS.
C     FOR INPUT LINE, SEE IF SOME OF NVAL-VALUES GET PRINTED HERE.
C       IF SO, SET COLPRT=.TRUE., NX = THE NUMBER OF VALUES (.LE.NCOL),
C                  AND (IX(I),I=1,NX) = INDICES OF VALUES ON THIS LINE.
C       IF NOT, (I.E., LINE BEYOND NVAL/NCOL+1) SET COLPRT=.FALSE.
C
C      SAMPLE USE:
C      ==========
C      DO 1000 I=1,LARGE_NUMBER
C      IF (COLPRT(NVAL,NCOL,I,IX,NX)) THEN
C        WRITE(6,*) (VALUES(IX(J)),J=1,NX)
C      ELSE
C        GO TO 2000
C     ENDIF
C 1000 CONTINUE
C 2000 .....
C
      LOGICAL COLPRT
      INTEGER NVAL,NCOL,LINE,IX(NCOL),NX
      INTEGER I,NLINE,IXTRA
C
      IF (NCOL.LT.1) THEN
        WRITE(6,*) '  COLPRT.  ERROR. NCOL =',NCOL
        COLPRT=.FALSE.
        RETURN
      ENDIF
C
      NLINE=NVAL/NCOL
      IXTRA=NVAL-NCOL*(NVAL/NCOL)
      IF (LINE.LE.NLINE) THEN
        IX(1)=LINE
        NX=1
        COLPRT=.TRUE.
        IF (NCOL.LT.2) GO TO 2000
        DO 1000 I=2,NCOL
        IX(I)=IX(I-1)+NLINE
        IF (IXTRA.GT.0) THEN
          IX(I)=IX(I)+1
          IXTRA=IXTRA-1
        ENDIF
        IF (IX(I).GT.NVAL) GO TO 2000
 1000   NX=I
 2000   RETURN
      ELSEIF (LINE.EQ.NLINE+1.AND.IXTRA.GT.0) THEN
        NX=IXTRA
        DO 3000 I=1,NX
 3000   IX(I)=LINE+(I-1)*(NLINE+1)
      ELSE
        COLPRT=.FALSE.
        RETURN
      ENDIF
      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
      SUBROUTINE CHKSTR(NUSED)
C     NEW ROUTINE FOR DYNAMIC STORAGE HANDLING (SG:1/21/93)
C     REVISED 2 FEB 94 TO REFLECT USAGE OF UPPER X() IN HIH2O
C
C     NUSED.LT.0 ON ENTRY RESETS HIH2O (HIGH-WATER MARK)
C                    ALSO RESETS MX=MAXMAX
C     HIH2O DOES NOT REFLECT USE OF THE TOP OF X()
C     NOTE: IT IS *NOT* NECESSARY TO PASS NUSED DOWN TO ALL ROUTINES
C           WHICH CALL CHKSTR.  HOWEVER, A NON-NEGATIVE VALUE MUST
C           BE SUPPLIED TO PREVENT RESETTING HIH2O, MX.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER HIH2O
      SAVE HIH2O,MAXMAX
C
C   DYNAMIC STORAGE COMMON BLOCK ...
      COMMON /MEMORY/ MX,IXNEXT,NIPR,IVLFL,X(1)
C
      DATA MAXMAX/0/
C
C     NEGATIVE INPUT NUSED IS SIGNAL TO RESET HIH2O, MX
      IF (NUSED.LT.0) THEN
        HIH2O=0
        MX=MAX0(MX,MAXMAX)
      ENDIF
C
C     STORAGE REQUIREMENT FOR CURRENT CALL IS ITOP=IXNEXT-1
      ITOP=IXNEXT-1
      IF (ITOP.GT.MX) THEN
        WRITE(6,600) ITOP,MX,MAXMAX
  600   FORMAT(/' CHKSTR.  CANNOT PROVIDE REQUESTED STORAGE.'/
     1          10X,'CURRENT REQUEST =',I12/
     2          10X,'CURRENT  LIMIT  =',I12/
     3          10X,'ORIGINAL LIMIT  =',I12)
        STOP
      ELSE
        MAXMAX=MAX(MAXMAX,MX)
        NEED=ITOP+(MAXMAX-MX)
        HIH2O=MAX(HIH2O,NEED)
        NUSED=HIH2O
      ENDIF
      RETURN
      END
      FUNCTION PLM(L,MM,X)
C
C     COMPUTES NORMALIZED ASSOC. LEGENDRE POLYNOMIALS BY RECURSION.
C     THE VALUES RETURNED ARE NORMALIZED FOR INTEGRATION OVER X
C     (I.E. INTEGRATION OVER COS THETA BUT NOT PHI).
C     NOTE THAT THE NORMALIZATION GIVES
C           PLM(L,0,1)=SQRT(L+0.5)
C           PLM(L,0,X)=SQRT(L+0.5) P(L,X)
C     FOR M.NE.0, THE VALUE RETURNED DIFFERS FROM THE USUAL
C           DEFINITION OF THE ASSOCIATED LEGENDRE POLYNOMIAL
C           (E.G. EDMONDS PAGES 23-24)
C           BY A FACTOR OF (-1)**M*SQRT(L+0.5)*SQRT((L-M)!/(L+M)!)
C     THUS THE SPHERICAL HARMONICS ARE
C          CLM = PLM * EXP(I*M*PHI) / SQRT(L+0.5)
C          YLM = PLM * EXP(I*M*PHI) / SQRT(2*PI)
C
C     PROGRAM OF R. NERF MODIFED BY S. GREEN.
C       MOD FEB. 82 BY S.G. ACCORDING TO R.T PACK'S SUGGESTION
C       FOR IMPROVED ACCURACY LARGE ABS(X)
C     MODIFIED AUG 88 BY S.G. TO KEEP RAT FROM BLOWING UP AT LARGE L
C       BY INCLUDING FACTOR (-Z)**L IN PRECEDING LOOP
C     ERROR IN STMT. NO. 211 CORRECTED MAY 93 BY SG
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      M=IABS(MM)
      IF ( M.GT.L .OR. L.LT.0)  GO TO 9999
      XL=L
      XM=M
      P1=1.D0
      P2=X
      IF(L.EQ.0) GO TO 20
C *** USE ALTERNATE RECURSION FOR LARGE M OR LARGE ABS(X)
      IF (M.EQ.0)  GO TO 10
      XTEST=0.49999D0*(1.D0+1.D0/XM)
      MTEST=L/3
      IF (M.GT.MTEST .OR. ABS(X).GT.XTEST)  GO TO 210
C ***
   10 DO 100 I=1,L
      XI=I
      P3=((2.D0*XI+1.D0)*X*P2-XI*P1)/(XI+1.D0)
      P1=P2
  100 P2=P3
      IF (M.EQ.0) GO TO 20
C AT END OF LOOP P1=P(L,0,X)
      IF (ABS(X).GT.1.D0)  GO TO 9999
      Z=SQRT(1.D0-X*X)
      IF (Z.LE.1.D-10) GO TO 999
  201 P2=(XL+1.D0)*(P2-X*P1)/Z
      DO 200 I=1,M
      XI=I
      P3=-2.D0*X/Z*P2*XI-(XL+XI)*(XL-XI+1.D0)*P1
      P1=P2
  200 P2=P3
      GO TO 20
C ***
C     BELOW RECURS DOWN IN M, AS SUGGESTED BY R. T PACK
  210 IF (ABS(X).GT.1.D0)  GO TO 9999
      Z=SQRT(1.D0-X*X)
      IF (Z.LE.1.D-10)  GO TO 999
C     CALCULATE RATIO OF FACTORIALS FOR PLL
      RAT=1.D0
      XI=0.D0
      DO 211 I=1,L
      XI=XI+1.D0
  211 RAT=-0.5D0*Z*(XL+XI)*RAT
C     CALCULATE PLL   ----  N.B. ABOVE INCL (-Z)**L COMPUTATION
C     P1=RAT*(-Z)**L
      P1=RAT
      IF (M.EQ.L)  GO TO 20
      P2=P1
      P1=-X*P2/Z
      IF (M.EQ.L-1)  GO TO 20
      LM1=L-M-1
C     RECUR DOWNWARD IN M
      DO 213 I=1,LM1
      MU=L-I-1
      XMU=MU
      P3=P2
      P2=P1
      P1=2.D0*(XMU+1.D0)*X*P2/Z+P3
  213 P1=-P1/((XL-XMU)*(XL+XMU+1.D0))
      GO TO 20
C ***
C     NORMALIZATION . . .
   20 XNORM=(2.D0*XL+1.D0)/2.D0
      IF (M.LE.0)  GO TO 1000
      XLM=XL+1.D0
      XLP=XL
      DO 1100 I=1,M
      XLM=XLM-1.D0
      XLP=XLP+1.D0
 1100 XNORM=XNORM/(XLM*XLP)
 1000 PLM=P1*SQRT(XNORM)
      RETURN
 9999 WRITE(6,699)  L,MM,X
  699 FORMAT('0 * * * ERROR.  ARGUMENT OUT OF RANGE FOR PLM(',2I6,D16.8,
     1  ' ).')
C     IF Z=0, THEN X=1 AND PLM(1.)=0 FOR M.GT.0.
C       IN THAT CASE, WE HAVE BRANCHED TO 999
  999 PLM=0.D0
      RETURN
      END
      FUNCTION THRJ(F1,F2,F3,G1,G2,G3)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE MUNG,X,Y
      DIMENSION X(302),Y(302)
      DATA MUNG/0/,MXIX/302/
      IF (MUNG.EQ.21) GO TO 69
      MUNG = 21
      X(1) = 0.D0
      DO 100 I = 1, 201
      A = I
      X(I+1) = LOG(A) +X(I)
      Y(I+1) = LOG(A)
  100 CONTINUE
   69 IF(F1-ABS(G1)) 1,13,13
   13 IF(F2-ABS(G2))1,14,14
   14 IF(F3-ABS(G3))1,15,15
   15 SUM=F1+F2+F3
      NSUM=SUM+.001D0
      IF(SUM-NSUM)2,2,1
    1 THRJ=0.D0
      RETURN
    2 IF(ABS(G1+G2+G3)-1.D-08)3,3,1
    3 IF(F1+F2-F3)1,4,4
    4 IF(F1+F3-F2)1,5,5
    5 IF(F2+F3-F1)1,6,6
    6 J1=2.D0*F3+2.001D0
      J2=F1+F2-F3+1.001D0
      J3=F1-F2+F3+1.001D0
      J4=-F1+F2+F3+1.001D0
      J5=F1+F2+F3+2.001D0
      J6=F1+G1+1.001D0
      J7=F1-G1+1.001D0
      J8=F2+G2+1.001D0
      J9=F2-G2+1.001D0
      J10=F3+G3+1.001D0
      J11=F3-G3+1.001D0
      IF(J5.GT.MXIX) THEN
        WRITE(6,601) J5,MXIX
  601   FORMAT(' *** DIMENSION ERROR IN THRJ - INDEX.GT.MXIX',2I5)
        STOP
      ENDIF
      R=0.5D0*(Y(J1)+X(J2)+X(J3)+X(J4)-X(J5)
     1+X(J6)+X(J7)+X(J8)+X(J9)+X(J10)+X(J11))
      SUM=0.D0
      F=-1
      KZ=-1
    7 KZ=KZ+1
      F=-F
      J1=KZ+1
      J2=F1+F2-F3-KZ+1.001D0
      IF(J2)20,20,8
    8 J3=F1-G1-KZ+1.001D0
      IF(J3)20,20,9
    9 J4=F2+G2-KZ+1.001D0
      IF(J4)20,20,10
   10 J5=F3-F2+G1+KZ+1.001D0
      IF(J5)7,7,11
   11 J6=F3-F1-G2+KZ+1.001D0
      IF(J6)7,7,12
   12 JMAX=MAX(J1,J2,J3,J4,J5,J6)
      IF(JMAX.GT.MXIX) THEN
        WRITE(6,601) JMAX,MXIX
        STOP
      ENDIF
      S=-(X(J1)+X(J2)+X(J3)+X(J4)+X(J5)+X(J6))
      SUM=SUM+F*EXP(R+S)
      GO TO 7
   20 INT=ABS(F1-F2-G3)+0.0001D0
      VAL=((-1.D0)**INT)*SUM/SQRT(2.D0*F3+1.D0)
      IF(ABS(VAL).LE.1.D-6) VAL=0.D0
      THRJ=VAL
      RETURN
      END