C     PROGRAM TO COMPUTE PRESSURE BROADENING  CROSS SECTIONS,
C       SIG(JA,JB;JA1,JB1;K) (CF SHAFER & GORDON, JCP 58, 5422 (1973))
C       FROM S-MATRICES SAVED ON TAPE BY SCATTERING PROG.
C     COMPATBILE WITH MOLSCAT SAVE TAPES THROUGH VERSION 14. 
C     N.B. S-MATRICES AT MATCHING PAIRS OF ENERGIES FOR REQUESTED
C       JA,JB MUST BE PRESENT ON THE SAVE TAPE.
C     THIS PROGRAM IS PARTICULARLY USEFUL FOR CALCULATING ADDITIONAL
C       CROSS SECTIONS (E.G. LINE COUPLING) AFTER AN INITIAL MOLSCAT
C       RUN, OR FOR CALCULATING VALUES FROM A TAPE ASSEMBLED FROM
C       MULTIPLE RUNS AT DIFFERENT JTOT RANGES.
C     THIS PROGRAM CALLS THE REGULAR MOLSCAT PRESSURE BROADENING
C       ROUTINES, BUT USES S-MATRICES FROM A SAVE TAPE RATHER THAN
C       CALCULTAING THEM.
C     IF CALCULATING MANY CROSS SECTIONS YOU MIGHT WANT TO INCREASE
C       PARAMETER NREC IN SUBROUTINE DASIZE TO ALLOW FOR MORE
C       INTERNAL WORKING STORAGE.
C    ------------------------------------------------------------------
C
C     NAMELIST &INPUT PARAMETERS  ---
C       NLPRBR (EQUIV TO IFLS), LINE, LTYPE, ARE STANDARD MOLSCAT
C          &INPUT VALUES -- SEE MOLSCAT DOCUMENTATION
C       ILSU IS NO LONGER USED.
C       IFEGEN FORCED TO ZERO HERE AS ENERGIES CANNOT BE ADDED
C       JMOD   - PRINTS ACCUM PRBR SIG ONLY FOR MOD(JTOT,JMOD)=0
C       JPRINT - SUPPRESSES PRINTING OF PRBR SIGS FOR JTOT<JPRINT
C       JSTEPZ - IF > 0 CONTRIBUTION FROM EVERY JSTEPZ'TH JTOT ONLY.
C    ------------------------------------------------------------------
C
C     PROGRAM HISTORY -
C     12 FEB 1985 MODIFIED TO SKIP SOME INPUT JTOT ON JSTEPZ PARM
C     01 MAR 1989 JMOD,JPRINT ADDED TO INPUT: CONTROL PRINTING
C                 DA FILE (11) SIMULATED BY IN-CORE STORAGE
C     25 JUN 1992 ADDS ABILITY TO READ UNFORMATTED (V11 AND LATER)
C                 SAVE TAPES.  CONTROLLED BY LFMT.
C     02 SEP 1994 UPDATED TO VERSION 14; N.B. DUMMY CHKSTR USED
C                 WHICH WILL PRECLUDE ITYPE=3 CALCS.
C    ------------------------------------------------------------------
C
C     KNOWN BUG:
C       IMPLEMENTAION LIMIT ON ITYPE=3; JLEVEL IS NOT SET PROPERLY.
C    ------------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER(MXNRG=100,MXCH=500,MXLVQN=700,MXLN=300,MXLEV=1000)
C       MXNRG IS MAX. NO. OF ENERGIES
C       MXCH IS MAX. NO OF CHANNESL PER S-MATRIX
C       MXLEV IS DIM OF ELEVEL - MUST BE CONSISTENT W/ CMBASE
C       MXLVQN IS MAX. NLEV*NQN (DIMENSIN OF JLEV)
      INTEGER LABEL(20),ITYPE,NLEV,NQN,NNRG,NOPEN,
     1  JLEV(MXLVQN),J(MXCH),L(MXCH)
      DIMENSION ENERGY(MXNRG),WVEC(MXCH),
     1          SREAL(MXCH*MXCH),SIMAG(MXCH*MXCH)
C     PRBR VARIABLES
      INTEGER PRNTLV,LINE(4*MXLN),LTYPE(MXLN)
      EQUIVALENCE(NLPRBR,IFLS)
      DIMENSION IC(MXCH),IL(MXCH),NBASIS(MXCH)
      DIMENSION IC1(MXCH),IL1(MXCH)
      LOGICAL LFMT
C
C  INFORMATION ORIGINALLY PASSED AS ENTRY PRBRBS, NOW IN COMMON
      LOGICAL MPLMIN
C     VERSION 14 DIMENSIONS FOR /CMBASE/
      COMMON /CMBASE/ ROTI(12),ELEVEL(1000),EMAX,WX(2),SPNUC,
     1   NLEVEL,JLEVEL(4000),MISC(30)
      EQUIVALENCE (IDENT,MISC(28))
      COMMON /PRBASE/ ITYPE,NQN,NLEV,MVALUE,IEXCH,MPLMIN
C
C   DYNAMIC STORAGE COMMON BLOCK; HOLDS NIPR NEEDED IN DASIZE
      COMMON /MEMORY/ MX,IXNEXT,NIPR,IVLFL,X(1)
C
      NAMELIST/INPUT/LINE,LTYPE,ILSU,PRNTLV,IFLS,NLPRBR,IFEGEN
     1              ,JSTEPZ,JMOD,JPRINT,LFMT
C
C     INPUT UNIT (IT)
      DATA IT/10/
C     TOLERANCE FOR ENERGY CHECK
      DATA TOL/1.E-6/
C     FORMATTED/UNFORMATTED TAPE (DEFAULT UNFORMATTED)
      DATA LFMT/.FALSE./
C
C     FORMATS FOR SAVE TAPE (IN CASE LFMT=.TRUE.)
C
  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 VIA APPROPRIATE MACHINE DEPENDENT METHOD
C     CALL XUFLOW(0)                 
C
      NIPR=2
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *  READ BEGINNING OF TAPE.                                      *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
      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) IT,LABEL,ITYPE,NLEV,NQN,URED,IP
  600 FORMAT('1 PROGRAM TO CALCULATE PRESSURE BROADENING FROM ',
     1       'SAVE TAPE ON UNIT',I4/'  LABEL =',20A4/
     2       '  ITYPE =',I3,5X,'NLEV, NQN =',2I4,5X,'URED =',F8.4,
     3       5X,'IPROGM =',I3)
C
C     JLEV CONTAINS THE QUANTUM NUMBERS FOR EACH 'LEVEL'
C       ((JLEV(I,NQ),NQ=1,NQN),I=1,NLEV)
C
      NSQ=NLEV*NQN
      IF (NSQ.LE.MXLVQN) GO TO 1991
      WRITE(6,619) NSQ,MXLVQN
  619 FORMAT(/'  *** ERROR. NLEV*NQN EXCEEDS MXLVQN',2I6)
      STOP
 1991 IF (LFMT) THEN
         READ(IT,101) (JLEV(I),I=1,NSQ)
      ELSE
         READ(IT) (JLEV(I),I=1,NSQ)
      ENDIF
C
      DO 1002 I=1,NLEV
 1002 WRITE(6,601) I,(JLEV(NLEV*(K-1)+I),K=1,NQN)
  601 FORMAT(/'  FOR LEVEL',I4,'  QUANTUM NOS. ARE',10I4)
C
      IF (IP.LT.3) THEN
        WRITE(6,*) ' *** NO IMPLEMENTATION FOR IPROGM.LT.3'
        STOP
      ENDIF
      IF (LFMT) THEN
         READ(IT,102) NLEVEL,(ELEVEL(I),I=1,NLEVEL)
      ELSE
         READ(IT) NLEVEL,(ELEVEL(I),I=1,NLEVEL)
      ENDIF
      IF (NLEVEL.LE.MXLEV)  GO TO 1000
      WRITE(6,699)
  699 FORMAT(/'  *** ERROR. NLEVEL. GT. MXLEV')
      STOP
 1000 WRITE(6,603) NLEVEL,(ELEVEL(I),I=1,NLEVEL)
  603 FORMAT(' ',I4,'   ROTATIONAL LEVELS'/(' ',8E16.8))
C
C     NNRG IS NUMBER OF ENERGIES FOR WHICH SCATTERING CALC WAS DONE.
C     ENERGY(I),I=1,NNRG ARE THE ENERGIES (IN (1/CM))
      IF (LFMT) THEN
         READ(IT,102) NNRG,(ENERGY(I),I=1,NNRG)
      ELSE
         READ(IT) NNRG,(ENERGY(I),I=1,NNRG)
      ENDIF
      IF (NNRG.LE.MXNRG) GO TO 1881
      WRITE(6,618) NNRG,MXNRG
  618 FORMAT(/'  *** ERROR. NNRG EXCEEDS MXNRG',2I6)
      STOP
 1881 WRITE(6,602) NNRG,(ENERGY(I),I=1,NNRG)
  602 FORMAT(/,1X,I4,'  ENERGIES ARE'/(' ',8D16.8))
C
C      * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C                                                                      *
C     SET UP FOR PRESSURE BROADENING; READ &INPUT AND CALL PRBRIN      *
C                                                                      *
C      * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
      ILSU=11
      JSTEPZ=0
      JMOD=1
      JPRINT=999
      PRNTLV=3
      DO 1101 I=1,MXLN
 1101 LTYPE(I)=-1
      DO 1102 I=1,MXCH
 1102 NBASIS(I)=I
C     MPLMIN IS TRUE IN CURRENT VERSIONS OF MOLSCAT
      MPLMIN=.TRUE.
C     FORCE RM=1., SINCE WVEC ON TAPE HAS BEEN CONVERTED TO 1/ANG.
      RM=1.
      READ(5,INPUT,END=9000)
C     FORCE IFEGEN TO ZERO SINCE ENERGIES CANNOT BE ADDED TO TAPE.
      IFEGEN=0
C 
      CALL PARMPR(JSTEPZ,JMOD,JPRINT)
      CALL PRBRIN(IFLS,LINE,LTYPE,MXLN,ILSU,NNRG,ENERGY,NNRG,IFEGEN,
     1            JLEV,PRNTLV)
C
C      * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     START LOOP OVER JTOT / ENERGY VALUES . . .                       *
C     JTOT,INRG,ENERGY(INRG) HEAD EACH BLOCK OF S-MATRICES.            *
C     N.B.  SAME JTOT,INRG MAY BE REPEATED FOR NON-INTERACTING PARITIES*
C      * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     JTOLD PROCESSING USED TO GET JSTEP.
      JSTEPX=-1
      JTOLD=-1
 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 ***
C     FOR COUPLED STATES CASES, SET MVALUE; IN MOLSCAT, DONE IN BASE.
C       CODE HERE ASSUMES MPLMIN=.TRUE. AS IN CCP6 VERSION 9
      MVALUE=M-1
C ***
C     IF THIS IS A NEW JTOT, OUTPUT PRBR CROSS SECTIONS
      IF (JTOT.EQ.JTOLD) GO TO 2011
      IF (JTOLD.EQ.-1) GO TO 2011
      JSTEP=JTOT-JTOLD
      IF (JSTEP.EQ.JSTEPX) GO TO 2911
      IF (JSTEPX.EQ.-1) GO TO 2912
      WRITE(6,691) JSTEP,JSTEPX
  691 FORMAT(/'  ***** ERROR.  JSTEP INCONSISTENCY.  CURRENT, OLD',2I4)
      JSTEP=0
 2912 JSTEPX=JSTEP
 2911 WRITE(6,617) JTOLD
  617 FORMAT(//'  **** **** **** ACCUMULATED THROUGH JTOT =',I4)
C     SKIP OUTPUT IF JTOT (I.E., JTOLD) WAS SKIPPED FOR JSTEPZ
      IF (JSTEPZ.LE.1) GO TO 4001
      IF (JTOLD-JSTEPZ*(JTOLD/JSTEPZ).EQ.0) GO TO 4001
      WRITE(6,627) JSTEPZ
  627 FORMAT(/'  **** **** **** JTOT SKIPPED BECAUSE OF JSTEPZ =',I4)
      GO TO 2011
 4001 CONTINUE
C     (3/1/89)  CUT DOWN ON OUTPUT 
      IF (JTOLD-JMOD*(JTOLD/JMOD).NE.0.AND.JTOLD.LE.JPRINT) GO TO 2011
      IF (PRNTLV.GE.3)  CALL PRBOUT(JPB(JSTEP,JSTEPZ))
 2011 JTOLD=JTOT
      IF (PRNTLV.GE.4)  WRITE(6,662) JTOT,ECHK,IEXCH,WT,M
  662 FORMAT(' JTOT,ENERGY,IEXCH,WT,M-VAL=',I4,E16.8,I4,E16.8,I4)
C     N.B. JTOT .LT. 0  INDICATES A 'CHECKPOINT' RECORD. NO S-MATRICES.
      IF (JTOT.GE.0)  GO TO 2001
      JTOT=-JTOT
      IF (PRNTLV.GE.4) WRITE(6,698)  JTOT
  698 FORMAT(' *** NOTE. ROLLOUT MARKER ENCOUNTERED.  NROLL =',I6)
      GO TO 2000
C     BELOW CHECKS THAT ENERGY(INRG) CORRESPONDS WITH HEADER RECORD. . .
 2001 IF (ABS((ECHK-ENERGY(INRG))/ECHK).LE.TOL) GO TO 2002
      WRITE(6,697)  INRG,ECHK
  697 FORMAT('  *** WARNING.  FOR ',I4,'-TH ENERGY, ECHECK =',D16.8)
C
 2002 IF (LFMT) THEN
         READ(IT,104,END=9999) NOPEN,(J(I),L(I),WVEC(I),I=1,NOPEN)
      ELSE
        IF (IP.LT.14) THEN
           READ(IT,END=9999) NOPEN,(J(I),L(I),WVEC(I),I=1,NOPEN)
        ELSE
           READ(IT,END=9999) (J(I),L(I),WVEC(I),I=1,NOPEN)
        ENDIF
      ENDIF
      IF (NOPEN.LE.MXCH)  GO TO 2003
  696 FORMAT(/'  *** ERROR. NO. OF OPEN CHANNELS EXCEEDS',I6,
     1  '   JTOT,INRG =',2I6)
      WRITE(6,696)  MXCH,JTOT,INRG
      STOP
 2003 NSQ=NOPEN*NOPEN
C
C     GET S-MATRICES
      IF (LFMT) THEN
         READ(IT,105,END=9999) (SREAL(I),I=1,NSQ)
         READ(IT,105,END=9999) (SIMAG(I),I=1,NSQ)
      ELSE
         CALL SREAD(IT,NOPEN,SREAL,IEND)
         IF (IEND.NE.0) GO TO 9999
         CALL SREAD(IT,NOPEN,SIMAG,IEND)
         IF (IEND.NE.0) GO TO 9999
      ENDIF
C     ALLOW FOR SKIPPING THIS JTOT ON JSTEPZ
      IF (JSTEPZ.LE.1) GO TO 4000
      IF (JTOT-JSTEPZ*(JTOT/JSTEPZ).EQ.0) GO TO 4000
      IF (PRNTLV.GE.4) WRITE(6,693) JTOT,INRG,M,JSTEPZ
  693 FORMAT(' JTOT,INRG,M =',3I4,'  SKIPPED DUE TO JSTEPZ =',I4)
      GO TO 2000
C
C      * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     PROCESS THIS SET OF S-MATRICES FOR PRESSURE BROADENING CROSS SECT*
C      * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
 4000 CALL PRBR(JTOT,M,NOPEN,INRG,RM,NBASIS,J,L,WVEC,SREAL,SIMAG,IC,IL,
     1          IC1,IL1,
     2          JLEV,MXPAR,WT,PRNTLV,ILSU)
C
C     GO BACK FOR MORE JTOT, INRG SETS . . .
      GO TO 2000
C
 9999 WRITE(6,*) ' *** UNEXPECTED EOF ENCOUNTERED ON UNIT (10)'
      GO TO 9001
C     END OF FILE ON (10)
 9000 WRITE(6,695)
  695 FORMAT(/'  *** NOTE.  END OF FILE REACHED ON UNIT (10).')
      WRITE(6,617) JTOLD
 9001 CALL PRBOUT(JPB(JSTEP,JSTEPZ))
      CALL DACLOS
C
      STOP
      END
      SUBROUTINE PARMPR(JSTEPZ,JMOD,JPRINT)
      IF (JSTEPZ.GT.1) WRITE(6,600) JSTEPZ
  600 FORMAT('0 ***'/'  *** WILL SKIP SOME JTOT BASED ON JSTEPZ =',I4/
     1       '  ***')
      IF(JPRINT.LE.1) RETURN
      WRITE(6,601) JPRINT,JMOD
  601 FORMAT('0 ***'/'  *** FOR JTOT .LT.',I4,'  ONLY MOD(JTOT,',I2,
     1       ')=0 WILL BE PRINTED.'/'  ***')
      RETURN
      END
      FUNCTION JPB(JXTEP,JXTEPZ)
      JSTEP=ABS(JXTEP)
      JSTEPZ=MAX(ABS(JXTEPZ),1)
      JBIG=MAX(JSTEP,JSTEPZ)
      JSMALL=MIN(JSTEP,JSTEPZ)
      IF (JBIG.EQ.JSMALL) GO TO 2000
      IF (JBIG-JSMALL*(JBIG/JSMALL).EQ.0) GO TO 3000
C     HERE THERE IS NO COMMON DENOMINATOR, USE PRODUCT
      JPB=JBIG*JSMALL
      GO TO 9000
C     HERE JBIG IS EVEN MULTIPLE OF JSMALL, USE BIGGER
 3000 JPB=JBIG
      GO TO 9000
C     HERE BOTH ARE SAME, USE EITHER
 2000 JPB=JBIG
 9000 RETURN
      END
      SUBROUTINE CHKSTR(NUSED)
      WRITE(6,*) ' *** CHKSTR.  ERROR'
      WRITE(6,*) ' ***          PRBR_SAVE IMPLEMENTATION DOES NOT USE',
     1           ' /MEMORY/'
      WRITE(6,*) ' ***          TERMINATING.'
      STOP
      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 PRBR(JTOT,M,N,INRG,RM,
     1          NBASIS,LEV,L,WVEC,SREAL,SIMAG,IC,IL,IC1,IL1,
     2          JLEV,MXPAR,WGHT,PRINT,ILSU)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE
C
C ***   AUG 76 NEW COUPLED STATES TREATMENT (KOURI ET AL.)
C ***   JUL 86 (CCP6 VERSION 9)  MOD 26 AUG 86 TO GET MXREC PROPERLY.
C ***                            AND MOD 21 OCT 86 : EDIFMX
C ***   OCT 86 VERSION FOR 'OFF-DIAGONAL' CROSS SECTIONS
C ***   JAN 87 CHANGES TO GET MPLMIN HANDLING CORRECT FOR ITYPE=25,26
C ***          AND ADD JSTEP TO ENTRY PRBOUT (REQUIRES CHANGE IN DRIVER)
C ***   MAR 87 CORRECTIONS FOR ITYPE=26
C ***   DEC 88 INCLUDE ITYPE=7 AND Q=0
C ***   MAR 89 HAS 'IN-CORE' D.A. SIMULATION
C ***          (NEED SUBROUTINE DASIZE/ENTRIES DARD1,DARD2,DAWR1,DAWR2)
C ***   JUL 92 REMOVES ALL REFERENCES TO LCSOLD (OLD, INCORRECT,
C ***          FORMULATION FOR COUPLED STATES:  SEE, E.G.,
C ***               GREEN, ET AL. JCP, 66, 1409 (1977))
C ***          CALLS TO ENTRIES (IN PRBR3) ALSO HAVE BEEN TRAPPED THERE.
C ***   JUN 93 FIXES BUG IN PRBR3 AND USES /MEMORY/ TO ELIMINATE LIMITS.
C ***   AUG 94 V14: ENTRY PRBCNT ADDED AND COMMON CMBASE CHANGED
C
C     CALCULATES SIGMA(JA1,JB1;JA,JB;K)
C       WHERE A/B INDICATE INITIAL/FINAL SPECTRAL LINES,
C             A1/B1 ARE AFTER COLLISION,  AND K IS TENSOR ORDER
C       SEE, E.G., SHAFER AND GORDON, JCP 58, 5422 (1973).
C
C       SUPPOSED TO BE UPWARD COMPATIBLE IF LDIAG=.TRUE.:
C       LDIAG=.TRUE.  TAKES *OLD* INPUT LINE=LEVA,LEVB, LEVA,LEVB, ... ,
C                     AND SETS LEVA1=LEVA, LEVB1=LEVB FOR ALL LINES.
C       LDIAG=.FALSE. INPUT IS LINE=LEVA,LEVB,LEVA1,LEVB1,
C                                   LEVA,LEVB,LEVA1,LEVB1, ...
C       N.B. LDIAG FORCED TO TRUE FOR ITYPE=3 CALCULATIONS.
C
C     ENTRY PRBRIN ACCEPTS &INPUT DATA AND SETS UP PRES. BROAD. CALC.
C     ENTRY PRBOUT PRINTS OUT ACCUMULATED SIGR, SIGI.
C     ENTRY PRBCNT FINDS WHETHER AN S-MATRIX WILL BE USED FOR PB CALC
C
C     PRBR SPECIFICATIONS --------------------------------------
C
      DIMENSION NBASIS(1),LEV(1),L(1),IC(1),IL(1),IC1(1),IL1(1),
     1   JLEV(NLEV,NQN)
      DIMENSION WVEC(1),SREAL(1),SIMAG(1)
C
C     JTOT IS TOTAL ANGULAR MOMENTUM
C     M = 0 FOR LAST PARITY STEP AT THIS JTOT.
C     N IS NUMBER OF OPEN CHANNELS, DETERMINES DIMENSION OF VECTORS.
C     INRG IS INDEX FOR ENERGY VALUES
C     RM IS SCALING FACTOR FOR RADIAL WAVEFUNCTION.
C     NBASIS (I) POINTS TO LEV,L VALUES FOR ITH OPEN CHANNEL.
C     LEV IS VECTOR OF BASIS SET LEVELS
C     L IS VECTOR OF BASIS ORBITAL ANGULAR MOMENTA.
C     WVEC IS VECTOR OF WAVEVECTORS
C     SREAL(N,N) IS REAL PART OF S MATRIX.
C     SIMAG(N,N) IS IMAGINARY PART OF S MATRIX.
C
      LOGICAL ITYPE3,EPM,LCSNEW,MPLMIN,LCSSYM
      INTEGER JT(2)
C
C     PRBRIN SPECIFICATIONS ------------------------------------
C
      INTEGER NLPRBR,MXLN,LINE(MXLN),ILSU,NNRG,PRINT,MXNRG,IFEGEN
      INTEGER T(MXLN)
      DIMENSION ENERGY(NNRG)
C
C     NLPRBR =0 FOR NO LINE SHAPE CALC.
C          =N (GT.0) GIVES NO. OF LINES FOR WHICH TO COMPUTE L.S.
C     LINE(4*I-3),...  ,I=1,NLPRBR IS LEVEL DATA FOR LINES.
C     ILSU (NOW REDUNDANT) WAS DIRECT ACCESS FILE FOR WORKING STORAGE
C     ENERGY(NNRG) ARE ENERGIES AT WHICH S MATRIX IS CALCULATED.
C     MXNRG IS MAXIMUM DIMENSION OF ENERGY ARRAY
C     IFEGEN .GT. 0  REQUESTS GENERATION OF ADDITIONAL ENERGY VALUES.
C     PRINT IS INTEGER PRINT CONTROL.
C
      LOGICAL NOCALC,PF,NDEBUG
      LOGICAL LDIAG,EXISTS,LDIAGX
      CHARACTER*8 PTP(3)
C     STORAGE DIMENSIONED FOR NO. OF LINES = MAXLN.
      DIMENSION LN(400,9)
      DIMENSION EREL(400),SIGR(400),SIGI(400)
      DIMENSION P(2),PRTY(4)
C
C  INFORMATION ORIGINALLY PASSED AS ENTRY PRBRBS, NOW IN COMMON
C
      COMMON /CMBASE/ ROTI(12),ELEVEL(1000),EMAX,WT(2),SPNUC,
     1   NLEVEL,JLEVEL(4000),MISC(26),JHALF,IDENT,MXJL,MXEL
      COMMON /PRBASE/ ITYPE,NQN,NLEV,MVALUE,IEXCH,MPLMIN
      COMMON /ASSVAR/ IDA
C
C     NLEV AND NLEVEL ARE NO. OF LEVELS IN BASIS SET.
C     JLEV AND JLEVEL ARE QUANTUM NUMBERS FOR THESE LEVELS.
C     ELEVEL ARE ENERGIES OF THESE LEVELS.
C
C   DYNAMIC STORAGE COMMON BLOCK ...
      COMMON /MEMORY/ MX,IXNEXT,NIPR,IVLFL,X(1)
C
C     --- DATA INITIALIZATIONS ---
C
      DATA PTP/' Q = 0  ',' DIPOLE ',' RAMAN  '/
      DATA P/1.D0,-1.D0/, PRTY/1.D0,-1.D0,-1.D0,1.D0/
C *** BELOW REPLACES JMH'S CRITERION OF 1.D-10 FOR ENERGY DIFFERENCE
C *** SMALLER VALUE MAY BE NEEDED FOR RESONANCE CALCULATIONS.
      DATA EDIFMX/5.D-6/
C     FOR COMPATBILITY WITH OLD INPUT, SET LDIAG=.TRUE.
      DATA LDIAGX/.FALSE./
C     IF NDEBUG .EQ. .FALSE. CHECK FOR 'IMPOSSIBLE' NUMBERS OF MATCHES.
      DATA NDEBUG/.FALSE./
C     DIMENSION LIMITATION ...
      DATA MAXLN/400/
C     FOR CHECKING OVER-WRITE OF "DA FILE"
      DATA JCHKSV/-1/
C
C     STATEMENT FUNCTION (LOGICAL)
      EXISTS(I) = I.GT.0 .AND. I.LE.NLEVEL
C
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
      IF (NOCALC) RETURN
      IF (JCHKSV.EQ.-1) JCHKSV=JTOT
      DO 3000 IA=1,2
C     IA=1 CHECKS FOR USE OF THIS JTOT,INRG WITH J(ALPHA).
C     IA=2 FOR J(BETA).
      IB=3-IA
C     FIND LINES, I, WHICH USE THIS INRG, JTOT S MATRIX.
      IKEEP=0
      DO 3100 I=1,NLINE
      IF (LN(I,IA+3).NE.INRG) GO TO 3100
      K=LN(I,3)
      JDIFMX=K
      IF (LCSNEW) JDIFMX=0
      JDM=MAX(JDM,JDIFMX)
      IF (ITYPE3) GO TO 3211
C     FOR ITYPE=1,2,5 GET J-VALUE FROM 1ST COL OF JLEV.
      JA=JLEV(LN(I,1),1)
      JB=JLEV(LN(I,2),1)
      JA1=JLEV(LN(I,8),1)
      JB1=JLEV(LN(I,9),1)
C     PARITY FACTOR FOR CS WITH MPLMIN; THIS IS NORMALLY +1.
      F3PJ=PARITY(JA+JA1+JB+JB1)
C     FIND BASIS FNS. CORRESPONDING TO JA/JA1 (JB/JB1) AND GET L VALUES.
C        ROWS=>JA1,IC1,IL1    COLS=>JA,IC,IL
C     FOR DIAG CASE (JA=JA1), IC/IC1 AND IL/IL1 HAVE SAME VALUES.
      NLVAL=0
      NLVAL1=0
      DO 3200 II=1,N
      JJ=NBASIS(II)
      IF (LEV(JJ).NE.LN(I,IA)) GO TO 3201
      NLVAL=NLVAL+1
      IC(NLVAL)=II
      IL(NLVAL)=L(JJ)
 3201 IF (LEV(JJ).NE.LN(I,IA+7)) GO TO 3200
      NLVAL1=NLVAL1+1
      IC1(NLVAL1)=II
      IL1(NLVAL1)=L(JJ)
 3200 CONTINUE
      GO TO 3212
C
C     FOR ITYPE=3 GET J-VALUE FROM JLEVEL.  RECALL J1,J2 PACKED IN ORDER
C
 3211 JA=JLEVEL(2*LN(I,1)-1)
      JB=JLEVEL(2*LN(I,2)-1)
C     BELOW MAY BE NEEDED FOR COMPATIBILITY IN OFF-DIAG CODE
      JA1=JA
      JB1=JB
C     ALLOCATE TEMPORARY STORAGE FOR SR,SI,TR,JBAR,ISTB,NBLK,LVAL
      NSQ=N*N
      NINT=(N+NIPR-1)/NIPR
      IT1=IXNEXT
      IT2=IT1+NSQ
      IT3=IT2+NSQ
      IT4=IT3+NSQ
      IT5=IT4+NINT
      IT6=IT5+NINT
      IT7=IT6+NINT
      IXNEXT=IT7+NINT
C     WRITE(6,*) ' IT1-IT7,IXNEXT',IT1,IT2,IT3,IT4,IT5,IT6,IT7,IXNEXT
      NUSED=0
      CALL CHKSTR(NUSED)
      CALL PRBR3(N,SREAL,SIMAG,JTOT,NLEV,NQN,JLEV,NBASIS,LEV,L,NPACK,
     1   LN(I,IA),NLVAL,IC,IL,
     2   X(IT1),X(IT2),X(IT3),X(IT4),X(IT5),X(IT6),X(IT7))
C     N.B. ONLY SR,SI NEED TO BE KEPT, IXNEXT COULD BE REDUCED HERE ...
C     IXNEXT=IT3
C>>SG 1 JUN 93: NLVAL1,IL1,IC1 MUST BE ALSO BE SET (CF. DIAGONAL CASE).
      NLVAL1=NLVAL
      DO 3311 II=1,NLVAL
      IL1(II)=IL(II)
 3311 IC1(II)=IC(II)
C<<SG 1 JUN 93
C
 3212 IF (NLVAL.LE.0 .OR. NLVAL1.LE.0) GO TO 3100
C     GET WAVE NUMBER POINTED TO BY IPWVEC (SHOULD BE SAME FOR ALL IC'S)
      IPWVEC=NBASIS(IC(1))
      FK=WVEC(IPWVEC)
      FK=FK/RM
      FK=PI/(FK*FK)
C
C     CHECK FOR AVAILABILITY OF J(BETA) MATRICES WITH WHICH TO CALC.
      IPTR=LN(I,IB+5)
 3400 IF (IPTR.LE.0) GO TO 5000
C
C     PROCESS CURRENT S MATRIX WITH PREVIOUSLY STORED MATRICES.
 4000 IDA=IPTR
C     READ(ILSU,REC=IDA) IPTIND,IIN,IAIN,JTIN,IPTBK,NNTRY
      CALL DARD1(IPTIND,IIN,IAIN,JTIN,IPTBK,NNTRY)
      IDA=IDA+1
      JT(IA)=JTIN
      JT(IB)=JTOT
      IF (LCSNEW) GO TO 4981
      MVALIN=1
      IF (IPTIND.EQ.-1) GO TO 4001
      WRITE(6,402) IPTR,IPTIND
  402 FORMAT(/'  * * * ERROR.  FOR IPTR =',I6,', ILLEGAL IPTIND =',I6)
 4003 IPTR=IPTBK
      GO TO 3400
C     'LCSNEW' CODE CAN USE JLEV(I,1) AS J-VALUE FOR ITYPE=1,2,5,6
 4981 MVALIN = IPTIND+1000
      IF (IABS(MVALIN).LE.JLEV(LN(I,IB),1)) GO TO 4001
      WRITE(6,405) IPTIND,JA,JB,IA
  405 FORMAT('  * * * ERROR. MVALUE LARGER THAN J-VALUE IN CS(NEW)',
     1  4I12)
      GO TO 4003
 4001 IF (IIN.EQ.I .AND. IAIN.EQ.IB) GO TO 4002
      WRITE(6,403) IIN,I,IAIN,IA
  403 FORMAT(/'  * * * ERROR.  IIN, I, IAIN, IA =',4I12)
      GO TO 4003
 4002 JDIF=JTOT-JTIN
      IF (JDIF.GE.0) GO TO 4004
      WRITE(6,404) JTOT,JTIN
  404 FORMAT(/'  * * * ERROR.  JTOT.LT.JTIN =',2I6)
      GO TO 4003
 4004 IF (JDIF.GT.JDIFMX) GO TO 5000
C
 4005 ONE=1.D0
C     N.B. FOR ITYPE=1,2,5,6,  JLEV(LEV,NQN)=LEV
      IF (.NOT.LCSNEW) GO TO 4006
C     SET UP FOR NEW COUPLED STATES CODE
      XJA=DBLE(JLEV(LN(I,IA),1))
      XJB=DBLE(JLEV(LN(I,IB),1))
      XJA1=DBLE(JLEV(LN(I,IA+7),1))
      XJB1=DBLE(JLEV(LN(I,IB+7),1))
      XN=DBLE(LN(I,3))
      XMA=DBLE(MVALUE)
      XMB=DBLE(MVALIN)
C     FACTORS BELOW FOR ITYPE=25,26 MPLMIN=.TRUE.
      TJ2P=1.D0
      F3P=F3PJ
      IF (.NOT.LCSSYM) GO TO 4006
      IXA=LN(I,IA)
      IXA1=LN(I,IA+7)
      IXB=LN(I,IB)
      IXB1=LN(I,IB+7)
      IF (ITYPE.EQ.26) GO TO 4016
      TJ2P=PARITY(JLEV(IXB,2)+JLEV(IXB,3)+JLEV(IXB1,2)+JLEV(IXB1,3))
      F3P=F3P*TJ2P*
     1    PARITY(JLEV(IXA,2)+JLEV(IXA,3)+JLEV(IXA1,2)+JLEV(IXA1,3))
      GO TO 4006
 4016 TJ2P=PRTY(JLEV(IXB,3)+1)*PRTY(JLEV(IXB1,3)+1)
      F3P=F3P*TJ2P*PRTY(JLEV(IXA,3)+1)*PRTY(JLEV(IXA1,3)+1)
C     FOR OLD CS CODE, ADD 1 TO SREAL ONLY FOR M=MP=0.
 4006 DO 3600 ID=1,NNTRY
C     READ(ILSU,REC=IDA) L1,L2,SR,SI
      CALL DARD2(L1,L2,SR,SI)
      IDA=IDA+1
      NMATCH=0
C
C     II LOOPS FOR L1 => INITIAL (UNPRIMED) I.E. JA:IL,IC (COLS)
C     JJ LOOPS FOR L2 => FINAL (PRIMED) I.E. JA1:IL1,IC1  (ROWS)
      DO 3700 II=1,NLVAL
      IF (L1.NE.IL(II)) GO TO 3700
      DO 3800 JJ=1,NLVAL1
      IF (L2.NE.IL1(JJ)) GO TO 3800
      NMATCH=NMATCH+1
      IX=(IC(II)-1)*N+IC1(JJ)
C *** N.B. ABOVE HAS REVERSED ROWS/COLS FROM OLD CODE,
C ***      BUT NOTE SYMMETRY OF S-MATRIX
      SRNOW=SREAL(IX)
      SINOW=SIMAG(IX)
C     THERE SHOULD ONLY BE 1 MATCH.  IF NDEBUG.EQ..FALSE. CHECK THIS.
      IF (NDEBUG) GO TO 4099
 3800 CONTINUE
 3700 CONTINUE
C
 3999 IF (NMATCH-1) 3600,4099,4098
C
 4098 WRITE(6,409) NMATCH,IPTIND,IIN,IAIN,JTIN,IPTBK,NNTRY,L1,L2
  409 FORMAT(/' * * * WARNING. NMATCH.GT.1 =',I3,'  FOR IPTIND,IIN,',
     1  'IAIN,JTIN,IPTBK,NNTRY,L1,L2 =',8I6)
      GO TO 3600
C
 4099 CONTINUE
      IF (EPM) GO TO 3630
      IF (LCSNEW) GO TO 3640
C     BELOW IS CLOSE COUPLING CODE . . .
      DJA1=JA1+JA1+1
      DJA=JA+JA+1
      F1=DBLE((2*JTOT+1)*(2*JTIN+1)) * SQRT(DJA1/DJA)
C     FOR ITYPE=3 2 INDICES PACKED INTO L1,L2 - REMOVE 2ND.
      L1N=L1-NPACK*(L1/NPACK)
      L2N=L2-NPACK*(L2/NPACK)
C *** JA,JA1 ADDED TO PARITY FACTOR OCT 86 FOR OFF-DIAGONAL CASE.
C ***  C.F. GREEN, C.P.L. 47, 119 (1977); COREY, MCCOURT, AND LIU,
C ***       J. PHYS. CHEM. 88, 2031 (1984).
      PARIT=PARITY(JA-JA1+L1N-L2N)
      F2=SIXJ(JA,K,L1N,JT(1),JB,JT(2))
C     N.B. DELTA FUNCTION CONTAINS ALL QUANTUM NUMBERS.
C          FOR SOME CASES CODE BELOW RECALCULATES SIXJ
      IF (L1.EQ.L2 .AND.
     1    LN(I,1).EQ.LN(I,8) .AND. LN(I,2).EQ.LN(I,9)) GO TO 3610
      F3=SIXJ(JA1,K,L2N,JT(1),JB1,JT(2))
      FR=0.D0
      GO TO 3620
C
 3610 F3=F2
      FR=ONE
      GO TO 3620
C
C     BELOW FOR EFFECTIVE POTENTIAL METHOD
 3630 F1=2*JTOT+1
      F2=1.D0
      F3=1.D0
      FR=1.D0
      PARIT=1.D0
      GO TO 3620
C
C     BELOW FOR COUPLED STATES (LCSNEW) CASE
 3640 DJA1=JA1+JA1+1
      DJA=JA+JA+1
      F1=DBLE(2*JTOT+1) * SQRT(DJA1/DJA)
      FR=1.D0
      IF (LN(I,1).NE.LN(I,8) .OR. LN(I,2).NE.LN(I,9)) FR=0.D0
      PARIT=PARITY(JB+JB1)
C     *** ACCOUNT FOR PLUS/MINUS MVALUE COUNTING. . .
C     *** AFTER OCT 76, THIS IS CONTROLLED BY MPLMIN PASSED FROM BASE.
      TJ1=THRJ(XJA,XJB,XN,XMA,-XMB,XMB-XMA)
     1   *THRJ(XJA1,XJB1,XN,XMA,-XMB,XMB-XMA)
      F2=TJ1
      F3=1.D0
      IF (.NOT.MPLMIN) GO TO 3620
      TJ2=THRJ(XJA,XJB,XN,XMA,XMB,-XMB-XMA)
     1   *THRJ(XJA1,XJB1,XN,XMA,XMB,-XMB-XMA)
      F2=F2 + TJ2 * TJ2P
      F3=1.D0 + F3P
      IF (MVALUE.EQ.0) F3=F3/2.D0
      IF (MVALIN.EQ.0) F3=F3/2.D0
      GO TO 3620
C
 3620 FCT=F1*F2*F3*FK*PARIT
      FR=(FR-(SR*SRNOW+SI*SINOW))*FCT
      FI=P(IA)*(SI*SRNOW-SR*SINOW)*FCT
      SIGR(I)=SIGR(I)+FR
      SIGI(I)=SIGI(I)+FI
      IF (I.NE.IKEEP .AND. PRINT.GE.4) WRITE(6,602)
  602 FORMAT(1X)
      IKEEP=I
      IF (PRINT.GE.4) WRITE(6,408) I,JT(IA),JT(IB),FR,SIGR(I),FI,SIGI(I)
  408 FORMAT('  LINE',I3,' BETWEEN JTOT =',I3,' AND',I3,E14.6,
     1  ' ADDED TO SIGR =',E14.6,8X,E14.6,' ADDED TO SIGI =',E14.6)
 3600 CONTINUE
C     THIS COMPLETES CALC. FOR ALL SAVED S MATRIX ELEMENTS IN THIS SET.
C     FOLLOW POINTERS BACKWARDS TO FIND OTHER S MATRICES.
      IPTR=IPTBK
      GO TO 3400
C
C     NO MORE STORED S MATRICES TO PROCESS.
C     STORE CURRENT MATRIX VALUES FOR FUTURE USE.
C
 5000 NNTRY=NLVAL*NLVAL1
      IF (NEXTDA+NNTRY.LT.MXREC) GO TO 4200
      WRITE(6,401) MXREC
  401 FORMAT(/' *** PRBR. WARNING: TEMPORARY STORAGE OF',I8,
     1  ' WORDS EXCEEDED.'/
     2  ' ***       WRAP-AROUND FEATURE WILL BE USED')
      IF (JTOT-JCHKSV.LE.JDM+1) THEN
        WRITE(6,491)
  491   FORMAT(' *** PRBR. POSSIBLE ERROR')
        WRITE(6,492) JCHKSV,JTOT,JDM
  492   FORMAT(11X,'INITIAL, FINAL JTOT WERE',2I5,'.',
     1    ' NEED  DELTA-JTOT =',I3)
C       ONE MIGHT WANT TO TERMINATE CALCULATION HERE.
      ELSE
        WRITE(6,493)
  493   FORMAT(' *** PRBR. PROBABLY OK.')
        WRITE(6,492) JCHKSV,JTOT,JDM
      ENDIF
      JCHKSV=JTOT
      NEXTDA=1
C
 4200 IPTIND=-1
      IF (LCSNEW) IPTIND=MVALUE-1000
      IPTBK=LN(I,IA+5)
C
C     IDA IS ASSOCIATED VARIABLE, SET TO NEXT AVAILABLE RECORD AFTER
C     READ/WRITE.
      IDA=NEXTDA
C     WRITE(ILSU,REC=IDA) IPTIND,I,IA,JTOT,IPTBK,NNTRY
      CALL DAWR1(IPTIND,I,IA,JTOT,IPTBK,NNTRY)
      IDA=IDA+1
      DO 4100 II=1,NLVAL
      DO 4100 JJ=1,NLVAL1
C     IX IS INDEX TO SREAL, SIMAG. CHANGED FOR OFF-DIAG CODE TO
C        CORRESPOND WITH CHOICE IN DO 3700 LOOP ABOVE
      IX=(IC(II)-1)*N+IC1(JJ)
      L1=IL(II)
      L2=IL1(JJ)
C     WRITE(ILSU,REC=IDA) L1,L2,SREAL(IX),SIMAG(IX)
      CALL DAWR2(L1,L2,SREAL(IX),SIMAG(IX))
 4100 IDA=IDA+1
C
      LN(I,IA+5)=NEXTDA
      NEXTDA=IDA
C     THIS FINISHES PROCESSING LINE=I FOR IA.
C     FOR ITYPE3 RESTORE SREAL,SIMAG AND RELEASE TEMP STORAGE
      IF (ITYPE3) THEN
C       WRITE(6,*) ' IT1,IT2,IXNEXT',IT1,IT2,IXNEXT
        CALL PRBR3R(N,SREAL,SIMAG,X(IT1),X(IT2))
C       RELEASE TEMPORARY STORAGE
        IXNEXT=IT1
      ENDIF
C
 3100 CONTINUE
C     THIS ENDS LOOP OVER I = 1,NLINE
 3000 CONTINUE
C     THIS ENDS LOOP OVER IA = ALPHA, BETA.
      RETURN
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C
      ENTRY PRBRIN(NLPRBR,LINE,T,MXLN,ILSU,NNRG,ENERGY,MXNRG,IFEGEN,
     1   JLEV,PRINT)
C
      LDIAG=LDIAGX
      NEXTDA=1
      JDM=0
C     NPACK IS USED TO PACK TWO INDICES (J1,J2) INTO ONE INTEGER VAR.
      NPACK=2**16
      PI=ACOS(-1.D0)
      IF (NLPRBR.GT.0) GO TO 1000
      NOCALC=.TRUE.
      RETURN
 1000 WRITE(6,111)
  111 FORMAT(//'  REVIEW OF REQUESTED PRESSURE BROADENING CALCULATION.')
      WRITE(6,610)
  610 FORMAT(/'  ****** THIS IS OFF-DIAGONAL VERSION (DEC 88) ******')
      IF (LDIAG) WRITE(6,609)
  609 FORMAT(9X,'LDIAG = TRUE OPTION TO HANDLE OLD INPUT')
C     COUPLED STATES (ITYPE=25,26 ADDED OCT 76) / 'OLD' CS CODE REMOVED
      LCSNEW=ITYPE.EQ.21 .OR. ITYPE.EQ.22
      LCSSYM=ITYPE.EQ.25 .OR. ITYPE.EQ.26
      LCSNEW=LCSNEW .OR. LCSSYM
      EPM=ITYPE.EQ.11 .OR. ITYPE.EQ.12 .OR. ITYPE.EQ.15 .OR. ITYPE.EQ.16
      ITYPE3=ITYPE.EQ.3
C     TRAP TYPES INCONSISTENT W/LDIAG=.FALSE.
      IF (LDIAG. OR. .NOT.(ITYPE3.OR.EPM)) GO TO 1992
      IF (ITYPE3) THEN
        LDIAG=.TRUE.
        WRITE(6,*) '  *** FOR ITYPE=3 LDIAG HAS BEEN FORCED TO TRUE.'
        WRITE(6,*) '  *** CHECK THAT &INPUT LINE DATA IS COMPATIBLE.'
      ELSE
        WRITE(6,197) ITYPE
  197   FORMAT(/'  *** LDIAG=.TRUE. INCOMPATIBLE WITH ITYPE =',I4,';',
     1    ' REQUEST CANCELED.')
        NOCALC=.TRUE.
        RETURN
      ENDIF
 1992 IF (ITYPE.EQ.1 .OR. ITYPE.EQ.2 .OR. ITYPE.EQ.5) GO TO 1990
      IF (ITYPE.EQ.6) GO TO 1990
      IF (ITYPE.EQ.7) GO TO 1990
      IF (ITYPE.EQ.31) GO TO 1990
      IF (ITYPE3 .AND. IDENT.EQ.0) GO TO 1991
      IF (LCSNEW) GO TO 1990
      IF (EPM) GO TO  1990
      WRITE(6,199) ITYPE
  199 FORMAT(/'  * * * WARNING.  REQUESTED PRESSURE BROADENING ',
     1  'CALCULATION  NOT SUPPORTED FOR ITYPE =',I4,'.'/
     2  18X,'REQUEST CANCELLED.')
      NOCALC=.TRUE.
      RETURN
 1991 WRITE(6,198)
  198 FORMAT(/'  * * * WARNING.  LIMITED IMPLEMENTATION FOR PRESSURE ',
     1  'BROADENING WITH ITYPE = 3 .AND. IDENT = 0')
 1990 NOCALC=.FALSE.
C     NPL IS NUMBER OF INDICES PER LINE
      NPL=4
      IF (LDIAG) NPL=2
      WRITE(6,100) NLPRBR
  100 FORMAT(/'  PRESSURE-BROADENING LINE-SHAPE CALCULATION REQUESTED'
     1  ,' FOR',I4,' LINES.')
C     N.B. DIMENSION ON LINE PASSED FROM DRIVER IS 2*MXLN
      IF ((NPL*NLPRBR)/2.LE.MXLN) GO TO 1100
      I=(2*MXLN)/NPL
      WRITE(6,101) NLPRBR,I
  101 FORMAT(/'  * * * WARNING. ',I5,' LINES REQUESTED, REDUCED TO MAX',
     1  ' OF',I3)
      NLPRBR=I
 1100 WRITE(6,110)
  110 FORMAT(/5X,' LINE   LEV(A) LEV(B)  LEV(A1) LEV(B1)')
      DO 1103 I=1,NLPRBR
      IST=NPL*(I-1)
 1103 WRITE(6,112) I,(LINE(IST+II),II=1,NPL)
  112 FORMAT(5X,I4,I8,I6,I9,I7)
      NLINE=0
      NEWE=NNRG
      DO 1200 I=1,NLPRBR
      PF=.FALSE.
      L1=LINE((NPL*(I-1))+1)
      L2=LINE((NPL*(I-1))+2)
      IF (LDIAG) GO TO 1104
      L3=LINE((NPL*(I-1))+3)
      L4=LINE((NPL*(I-1))+4)
      GO TO 1105
 1104 L3=L1
      L4=L2
 1105 IF (EXISTS(L1).AND.EXISTS(L2).AND.EXISTS(L3).AND.EXISTS(L4))
     1    GO TO 1150
      WRITE(6,108) L1,L2,L3,L4
  108 FORMAT(/'  * * * WARNING.  REQUESTED LINE FOR LEVELS',4I4,
     1  '  CANNOT BE PROCESSED - OUTSIDE NLEV RANGE.')
      GO TO 1200
 1150 IF (ITYPE3) GO TO 1151
C     FOR ITYPE=1,2,5 TAKE J-VALUE FROM JLEV(LEV,1)
      JA=JLEV(L1,1)
      JB=JLEV(L2,1)
      JA1=JLEV(L3,1)
      JB1=JLEV(L4,1)
      GO TO 1152
C     FOR ITYPE=3 TAKE J-VALUES FROM JLEVEL.  RECALL J1,J2 PACKED.
 1151 JA=JLEVEL(2*L1-1)
      JB=JLEVEL(2*L2-1)
      JA1=JA
      JB1=JB
 1152 K=T(I)
      K1=K
      IF (K.GE.0) GO TO 1106
C     IF 'LTYPE' NOT INPUT, CALCULATE FROM J-VALUES
      K=IABS(JA-JB)
      K1=IABS(JA1-JB1)
C     OTHER THAN DIPOLE AND RAMAN (Q=0,2) TRANSITIONS NOT IMPLEMENTED.
C           ALSO, MUST HAVE SAME 'K' FOR JA/JB AND JA1/JB1
 1106 IF ((K.EQ.0.OR.K.EQ.1.OR.K.EQ.2).AND.(K.EQ.K1)) GO TO 1300
      WRITE(6,102) JA,JB,JA1,JB1,I
  102 FORMAT(/'  * * * *  ONLY DIPOLE AND RAMAN TRANSITIONS IMPLEMENTED'
     1  ,'.  CANNOT PROCESS JA,JB,JA1,JB1=',4I4/
     2  10X,'SPECIFY LTYPE(',I3,' ) = 0, 1 OR 2 IN &INPUT.')
      GO TO 1200
 1300 NEPL=0
      IF (IFEGEN.LE.0) GO TO 1301
C     GENERATE ENERGY VALUES IF NECESSARY TO LINE AT REQUESTED E'S.
      DO 1401 II=1,NNRG
      EA=ENERGY(II)+ELEVEL(L1)
      DO 1501 JJ=1,NEWE
      ED=(ENERGY(JJ)-EA)/ENERGY(JJ)
      IF (ABS(ED).GT.EDIFMX) GO TO 1501
      JL=JJ
      GO TO 1502
 1501 CONTINUE
      NEWE=NEWE+1
      IF (NEWE.LT.MXNRG) GO TO 1201
 1202 WRITE(6,611) I,II
  611 FORMAT(/'  * * * WARNING.  CANNOT ADD ENERGY VALUE FOR LINE',I4,
     1  '  RELATIVE ENRGY NO.',I4)
      NEWE=NEWE-1
      GO TO 1401
 1201 ENERGY(NEWE)=EA
      JL=NEWE
 1502 EB=ENERGY(II)+ELEVEL(L2)
      DO 1101 JJ=1,NEWE
      ED=(ENERGY(JJ)-EB)/ENERGY(JJ)
      IF (ABS(ED).GT.EDIFMX) GO TO 1101
      JK=JJ
      GO TO 1102
 1101 CONTINUE
      NEWE=NEWE+1
      IF (NEWE.GE.MXNRG) GO TO 1202
      ENERGY(NEWE)=EB
      JK=NEWE
 1102 NEPL=NEPL+1
      IF (NLINE.LT.MAXLN) GO TO 1108
      WRITE(6,109) L1,L2,L3,L4,JL,JK
      GO TO 1200
 1108 NLINE=NLINE+1
      LN(NLINE,1)=L1
      LN(NLINE,2)=L2
      LN(NLINE,3)=K
      LN(NLINE,4)=JL
      LN(NLINE,5)=JK
      LN(NLINE,6)=0
      LN(NLINE,7)=0
      LN(NLINE,8)=L3
      LN(NLINE,9)=L4
      SIGR(NLINE)=0.D0
      SIGI(NLINE)=0.D0
      EREL(NLINE)=ENERGY(II)
      IF (PF) GO TO 1601
      WRITE(6,103) L1,JA,L2,JB,L3,JA1,L4,JB1,PTP(K+1)
      PF=.TRUE.
 1601 WRITE(6,104) EREL(NLINE),JL,JK
 1401 CONTINUE
      GO TO 1403
C
C     CHECK FOR AVAILABLE PAIRS OF ENERGIES WITH SAME INITIAL REL KE
C
 1301 DO 1400 II=1,NNRG
      EA=ENERGY(II)-ELEVEL(L1)
C *** 2 SEPT. 86 TO PREVENT CHOOSING NEGATIVE ENERGIES
      IF (EA.LE.0.D0) GO TO 1400
      DO 1500 JJ=1,NNRG
      EB=ENERGY(JJ)-ELEVEL(L2)
      ED=(EA-EB)/MAX(EA,EB,1.D0)
      IF (ABS(ED).GT.EDIFMX) GO TO 1500
      NEPL=NEPL+1
      IF (NLINE.LT.MAXLN) GO TO 1109
      WRITE(6,109) L1,L2,L3,L4,II,JJ
  109 FORMAT(/'  * * * WARNING.  NO MORE TABLE SPACE FOR LEVELS =',4I3,
     1  ', AND ENERGIES=',2I4,';   WILL BE SKIPPED.')
      GO TO 1200
 1109 NLINE=NLINE+1
      LN(NLINE,1)=L1
      LN(NLINE,2)=L2
      LN(NLINE,3)=K
      LN(NLINE,4)=II
      LN(NLINE,5)=JJ
      LN(NLINE,6)=0
      LN(NLINE,7)=0
      LN(NLINE,8)=L3
      LN(NLINE,9)=L4
      SIGR(NLINE)=0.D0
      SIGI(NLINE)=0.D0
      EREL(NLINE)=(EA+EB)/2.D0
      IF (PF) GO TO 1600
      WRITE(6,103) L1,JA,L2,JB,L3,JA1,L4,JB1,PTP(K+1)
  103 FORMAT(/'  LEVELS',I3,' (JA =',I3,' ), ',I3,' (JB =',I3,' )',
     1  ' **TO** LEVELS',I3,' (JA1 =',I3,' ),',I3,' (JB1 =',I3,' )',
     2  ' WILL BE PROCESSED FOR',A8,'RADIATION.')
      PF=.TRUE.
 1600 WRITE(6,104) EREL(NLINE),LN(NLINE,4),LN(NLINE,5)
  104 FORMAT(12X,'AT RELATIVE K.E. =',F18.9,' (1/CM)  WITH',I4,
     1 '-TH AND',I4,'-TH ENERGY VALUES RESPECTIVELY.')
 1500 CONTINUE
 1400 CONTINUE
 1403 IF (NEPL.GT.0) GO TO 1200
      WRITE(6,105) L1,L2
  105 FORMAT(/'  * * * WARNING. NO RELEVANT PAIRS OF ENERGY VALUES ',
     1  'FOR REQUESTED LEVELS =',2I4)
 1200 CONTINUE
C
      IF (IFEGEN.LE.0) GO TO 1701
C
C     REMOVE ENERGIES THAT ARE NOT NEEDED FOR PRESSURE-BROADENING
      NNRG=NEWE
      II=1
 1609 IF (II.GT.NNRG) GO TO 1615
      DO 1610 I=1,NLINE
      IF (LN(I,4).EQ.II .OR. LN(I,5).EQ.II) THEN
C       WE ARE USING THIS ENERGY. CHECK THE NEXT
        II=II+1
        GO TO 1609
      ENDIF
 1610 CONTINUE
C     REACH HERE IF IT IS NOT USED.  COMPRESS LIST & REVISE INDEXES
      NNRG=NNRG-1
      DO 1611 JJ=II,NNRG
      ENERGY(JJ)=ENERGY(JJ+1)
 1611 CONTINUE
      DO 1612 I=1,NLINE
      IF (LN(I,4).GT.II) LN(I,4)=LN(I,4)-1
 1612 IF (LN(I,5).GT.II) LN(I,5)=LN(I,5)-1
      GO TO 1609
C
C     SORT ENERGIES INTO DESCENDING ORDER, RENUMBER LN(,4) AND LN(,5)
C
 1615 DO 1809 JJ=1,NNRG
      II=0
      ETOP=0.D0
      DO 1805 I=JJ,NNRG
      IF (ETOP.GE.ENERGY(I)) GO TO 1805
      ETOP=ENERGY(I)
      II=I
 1805 CONTINUE
      ENERGY(II)=ENERGY(JJ)
      ENERGY(JJ)=ETOP
      DO 1808 I=1,NLINE
      KK=LN(I,4)
      IF (KK.EQ.JJ) LN(I,4)=II
      IF (KK.EQ.II) LN(I,4)=JJ
      KK=LN(I,5)
      IF (KK.EQ.JJ) LN(I,5)=II
 1808 IF (KK.EQ.II) LN(I,5)=JJ
 1809 CONTINUE
      WRITE(6,612) NNRG,(I,ENERGY(I),I=1,NNRG)
  612 FORMAT(//'  MODIFIED ENERGY LIST NOW CONTAINS',I4,'  VALUES,'
     1  /(15X,'ENERGY(',I3,') =',F18.9))
      WRITE(6,921)
  921 FORMAT(/'  LINE-SHAPE TABLES HAVE BEEN MODIFIED ACCORDINGLY.')
 1701 IF (NLINE.GT.0) GO TO 1700
      WRITE(6,106)
  106 FORMAT(/'  * * * WARNING. NONE OF REQUESTED LINES CAN BE ',
     1  'PROCESSED.  REQUEST CANCELLED.')
      NOCALC=.TRUE.
      RETURN
C     CODE FOR IN-CORE ROUTINES REPLACES OLDER CODE BELOW
 1700 CALL DAOPEN
      CALL DASIZE(ILSU,MXREC)
      RETURN
C1700 IF (ILSU.LE.0) GO TO 1900
C
C    OLD CODE FOR OPENING PRESSURE BROADENING SCRATCH FILE.
C    THIS IS NOW REPLACED BY IN-CORE STORAGE (JUN 92).
C    THE APPROPRIATE RECORD LENGTH IS MACHINE-DEPENDENT, AND
C    MUST BE SUFFICIENT TO STORE 6 INTEGERS. FOR EXAMPLE:
C
C       IBM:   LREC = 24 (BYTES)
C       VAX:   LREC =  6 (LONGWORDS)
C       CRAY:  LREC = 48 (BYTES)
C
C     LREC=24
C     OPEN(ILSU,STATUS='SCRATCH',ACCESS='DIRECT',FORM='UNFORMATTED',
C    1   RECL=LREC,ERR=1900)
C
C1900 WRITE(6,107) ILSU
C 107 FORMAT(/'  * * * WARNING.  UNABLE TO OPEN DIRECT ACCESS FILE ON ',
C    1   'CHANNEL',I3,'.  REQUESTED LINE-SHAPE CALCULATION CANCELLED.')
C     NOCALC=.TRUE.
C     RETURN
C     IF (ILSU.LE.100) WRITE(6,200) ILSU
C 200 FORMAT(/'  LINE SHAPE CALCULATION USES DIRECT ACCESS FILE ON ',
C    1   'CHANNEL',I4,' FOR TEMPORARY STORAGE')
C     RETURN
C
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
      ENTRY PRBCNT(INRG,LEV,N,IUSE)
C
C     SEE WHETHER THIS S-MATRIX WILL CONTRIBUTE TO ANY OF THE LINES
C     REQUESTED: IUSE=1 MEANS THAT IT WILL.
C
C     THIS CODE MAY NOT WORK FOR DIATOM-DIATOM: BE CAUTIOUS FOR NOW
      IF (ITYPE3) GO TO 5920
C
      DO 5910 I=1,NLINE
      DO 5910 IA=1,2
      IF (LN(I,IA+3).NE.INRG) GO TO 5910
      DO 5908 II=1,N
      L1=LEV(II)
      IF (L1.EQ.LN(I,IA) .OR. L1.EQ.LN(I,IA+7)) GO TO 5920
 5908 CONTINUE
 5910 CONTINUE
      IUSE=0
      RETURN
C
 5920 IUSE=1
      RETURN 
C
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
      ENTRY PRBOUT(JSTEP)
      IF (NOCALC) RETURN
      WRITE(6,671)
  671 FORMAT(/'  ACCUMULATED PRESSURE-BROADENING CROSS SECTIONS',
     1  ' (IN ANG**2)')
      FACTJ=1.D0
      IF (JSTEP.GT.1) THEN
        FACTJ=DBLE(JSTEP)/DBLE(JHALF)
        WRITE(6,672) FACTJ
  672   FORMAT(/'  *** NOTE *** CROSS SECTIONS MULTIPLIED BY',F5.1,
     1    ' TO ACCOUNT FOR JSTEP')
      ENDIF
      WRITE(6,673)
  673 FORMAT(/'  LINE LEV(A) LEV(B) LEV(A1) LEV(B1)  TYPE',8X,
     1  'EREL(1/CM)',11X,'RE(S)',12X,'IM(S)')
C
      DO 6996 I=1,NLINE
      SIGRJS=SIGR(I)*FACTJ
      SIGIJS=SIGI(I)*FACTJ
 6996 WRITE(6,601) I,LN(I,1),LN(I,2),
     1   LN(I,8),LN(I,9),PTP(LN(I,3)+1),EREL(I),SIGRJS,SIGIJS
  601 FORMAT(2I5,2I7,I8,3X,A8,F18.9,2D17.6)
      RETURN
      END
      SUBROUTINE PRBR3(N,SREAL,SIMAG,JTOT,NLEV,NQN,JLEV,NBASIS,LEV,L,
     &      NPACK,LINE,NLVAL,IC,IL,SR,SI,TR,JBAR,ISTB,NBLK,LVAL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION SREAL(N,N),SIMAG(N,N)
      DIMENSION SR(N,N),SI(N,N),TR(N,N),JBAR(N),ISTB(N),NBLK(N),LVAL(N)
      INTEGER JLEV(NLEV,NQN),NBASIS(N),LEV(N),L(N),IC(N),IL(N)
C
C     PRBR3 IS AN INTERFACE ROUTINE WHICH TRANSFORMS ITYPE=3 S-MATRICES
C     INTO FORMAT NECESSARY FOR PRESSURE BROADENING CALC..
C     THIS NECESSITATES TRANSFORM FROM J1,J2(J12)L,JTOT TO
C     J1(J2,L)JBAR,JTOT REPRESENTATION.  IC, IL MADE COMPATIBLE.
C
C     SREAL, SIMAG ARE STORED TEMPORARILY IN SR, SI WHICH MUST BE
C     RESTORED (VIA ENTRY PRBR3R) BEFORE RETURNING CONTROL FROM PRBR.
C     4 JUN 93: REMOVE LIMITS ON INTERNALLY DEFINED WORKING STORAGE
C
C     ENTRIES PRBRCS, PBCSSM, AND PBCSDA, WHICH HAD BEEN USED IN
C       'OLD' CS TRANSORMATION (SHOWN TO BE INCORRECT) ARE NOW
C       TRAPPED WITH ERROR MSG (7/92).  THEY SHOULD NOT BE CALLED.
C
      DATA ZERO/0.D0/
C
C     STATEMENT FUNCTION DEFINITION.
      Z(I)=I+I+1
C
C     SAVE SREAL, SIMAG IN SR, SI AND CLEAR SREAL, SIMAG.
 1000 DO 1100 I1=1,N
      DO 1100 I2=1,N
      SR(I1,I2)=SREAL(I1,I2)
      SREAL(I1,I2)=0.D0
      SI(I1,I2)=SIMAG(I1,I2)
 1100 SIMAG(I1,I2)=0.D0
C
C     PICK OUT ROWS OF S WHICH CORRESPOND TO 'LINE' AND GROUP INTO
C       BLOCKS WITH SAME L-VALUE.
C       NB IS NO. OF BLOCKS, IC IS TABLE OF ALL
C       (NLVAL) ROWS, IL HAS J12 VALUES, ISTB HAS (START-1) IN IC OF
C       EACH BLOCK, NBLK HAS NO. IN EACH BLOCK.
      NB=0
      NLVAL=0
      DO 2000 II=1,N
      JJ=NBASIS(II)
      IF (JLEV(LEV(JJ),NQN).NE.LINE)  GO TO 2000
C     CHECK TO SEE THAT THIS ISN'T ALREADY IN IC. . .
      IM=NLVAL
 2100 IF (IM.LE.0)  GO TO 2200
      IF (II.EQ.IC(IM))  GO TO 2000
      IM=IM-1
      GO TO 2100
C     IF NOT, START A NEW BLOCK. . .
 2200 NB=NB+1
      ISTB(NB)=NLVAL
      NBLK(NB)=1
      LVAL(NB)=L(JJ)
      NLVAL=NLVAL+1
      IC(NLVAL)=II
      IL(NLVAL)=JLEV(LEV(JJ),3)
C     PROCESS REMAINING TO FIND OTHER WITH SAME L-VALUE. . .
      II1=II+1
      IF (II1.GT.N)  GO TO 2000
      DO 2500 II2=II1,N
      JJ=NBASIS(II2)
      IF (JLEV(LEV(JJ),NQN).NE.LINE)  GO TO 2500
      IF (L(JJ).NE.LVAL(NB))  GO TO 2500
      NBLK(NB)=NBLK(NB)+1
      NLVAL=NLVAL+1
      IC(NLVAL)=II2
      IL(NLVAL)=JLEV(LEV(JJ),3)
 2500 CONTINUE
 2000 CONTINUE
      IF (NLVAL.LE.0)  RETURN
C
C     DO BOOKEEPING FOR TRANSFORMATION AND GET JBAR VALUES.
C     GET JA, J2 FROM 1ST MEMBER - SHOULD BE SAME FOR ALL.
      JJ=LEV(NBASIS(IC(1)))
      JA=JLEV(JJ,1)
      J2=JLEV(JJ,2)
      JLOW=IABS(JTOT-JA)
      JTOP=JTOT+JA
      DO 3000 II=1,NB
      LLL=LVAL(II)
      NJB=0
      JMIN=IABS(J2-LLL)
      JMAX=J2+LLL
      DO 3100 II2=JLOW,JTOP
      IF (II2.LT.JMIN)  GO TO 3100
      IF (II2.GT.JMAX)  GO TO 3100
C     ALL TRIANGLE INEQUALITIES SATISFIED.  INCLUDE JBAR
      NJB=NJB+1
      JBAR(ISTB(II)+NJB)=II2
 3100 CONTINUE
C     NO.(JBAR) SHOULD EQUAL NO.(J12)  FOR EVERY BLOCK.
      IF (NJB.EQ.NBLK(II) )  GO TO 3000
      WRITE(6,601) NJB,NBLK(II),II
  601 FORMAT('0 * * * ERROR.  NO.(JBAR) .NE. NO.(J12)',2I6,
     1  '  FOR BLOCK =',I4)
      STOP
 3000 CONTINUE
C
C     SET UP TRANFORM MATRIX TR(JBAR,J12)
      DO 3500 II1=1,NLVAL
      DO 3500 II2=1,NLVAL
 3500 TR(II1,II2)=0.D0
      DO 3600 II=1,NB
      LLL=LVAL(II)
      NTOP=NBLK(II)
      DO 3600 II1=1,NTOP
      JB=JBAR(ISTB(II)+II1)
      DO 3600 II2=1,NTOP
      J12=IL(ISTB(II)+II2)
 3600 TR(ISTB(II)+II1,ISTB(II)+II2)=SQRT(Z(JB)*Z(J12))
     1             *  SIXJ(JA,J2,JTOT,LLL,J12,JB)
C
C     TRANSFORM SR,SI TO SREAL, SIMAG
C     LOOP OVER L
      DO 3700 II1=1,NB
      NTOP1=NBLK(II1)
      IST1=ISTB(II1)
C     LOOP OVER L-PRIME
      DO 3700 II2=1,NB
      NTOP2=NBLK(II2)
      IST2=ISTB(II2)
C     LOOP OVER JBAR
      DO 3700 JB=1,NTOP1
      IRNEW=IC(IST1+JB)
C     LOOP OVER JBAR-PRIME
      DO 3700 JBP=1,NTOP2
      ICNEW=IC(IST2+JBP)
C     LOOP OVER J12
      DO 3700 J12=1,NTOP1
      IROLD=IC(IST1+J12)
C     LOOP OVER J12-PRIME
      DO 3700 J12P=1,NTOP2
      ICOLD=IC(IST2+J12P)
      FACTOR = TR(IST1+JB,IST1+J12) * TR(IST2+JBP,IST2+J12P)
      SREAL(IRNEW,ICNEW)=SREAL(IRNEW,ICNEW)+FACTOR*SR(IROLD,ICOLD)
      SIMAG(IRNEW,ICNEW)=SIMAG(IRNEW,ICNEW)+FACTOR*SI(IROLD,ICOLD)
 3700 CONTINUE
C
C     PACK L,JBAR INDICES INTO IL()
      II2=0
      DO 4000 II=1,NB
      NTOP=NBLK(II)
      LLL=LVAL(II)
      DO 4000 II1=1,NTOP
      II2=II2+1
 4000 IL(II2)=LLL*NPACK+JBAR(II2)
C     SREAL, SIMAG, IC, AND IL ARE NOW IN FORMAT EXPECTED BY PRBR
      RETURN
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C
      ENTRY PRBR3R(N,SREAL,SIMAG,SR,SI)
C     RESTORE SIMAG,SREAL IN CASE THEY ARE NEEDED FOR FURTHER PROCESSING
      DO 5000 II1=1,N
      DO 5000 II2=1,N
      SREAL(II1,II2)=SR(II1,II2)
 5000 SIMAG(II1,II2)=SI(II1,II2)
      RETURN
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     ENTRIES BELOW PROVIDED INTERFACE TO OLD (INCORRECT) CS CODE
C                  ----- DUMMIED 7/92 -----
C       ENTRY PRBRCS TRANSFORMS BODY-FIXED TO SPACE-FIXED S-MATRICES.
C       ENTRY PBCSSM RETURNS MATRIX ELEMENTS STORED INTERNALLY HERE.
C       ENTRY PBCSDA OUTPUTS STORED S-MATRICES ON DISK.
C
      ENTRY PRBRCS
C     ENTRY PRBRCS(N,SREAL,SIMAG,JTOT,MVALUE,IEXCH,WGHT,IPWVEC,
C    &     NLEV,NQN,JLEV,NBASIS,LEV,L,LINE,NLVAL,ITYPE)
      ENTRY PBCSSM
C     ENTRY PBCSSM(L1,L2,NDEBUG,NMATCH,SROUT,SIOUT)
      ENTRY PBCSDA
C     ENTRY PBCSDA(ILSU)
      WRITE(6,7300)
 7300 FORMAT(' *** ERROR.  PBBRCS, PBCSSM, OR PBCSDA CALLED.'/
     1       ' ***         THIS (INCORRECT) CS TRANSFORM NO LONGER',
     2       ' SUPPORTED.')
      STOP
      END
      SUBROUTINE DASIZE(ILSU,MXREC)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE MXUSED,IX
      PARAMETER (NREC=300000)
      DIMENSION IX(6,NREC)
      DIMENSION IR2(2),IS2(2)
      EQUIVALENCE (R,IR1,IR2(1)),(S,IS1,IS2(1))
      COMMON /ASSVAR/IDA
C
C   DYNAMIC STORAGE COMMON BLOCK ...
C     NEEDED FOR NIPR; PREVIOUSLY PASSED IN COMMON /INTPAC/
      COMMON /MEMORY/ MX,IXNEXT,NIPR,IVLFL,X(1)
C
      DATA MAX/NREC/
C
      MXREC=MAX
      ILSU=999
      WRITE(6,601) MAX
  601 FORMAT( ' *** *** NUMBER OF SIMULATED RECORDS =',I7)
      RETURN
C
      ENTRY DAOPEN
      MXUSED=0
      WRITE(6,600)
  600 FORMAT(/' *** *** IN-CORE DA SIMULATION ROUTINE HAS CONTROL.',
     1       /' *** *** DA FILE WILL NOT BE USED.')
C
      IF(NIPR.EQ.1 .OR. NIPR.EQ.2) GOTO 1000
      WRITE(6,602) NIPR
  602 FORMAT(' *** ERROR IN DASIZE/DAOPEN: NIPR =',I3,' INVALID')
      STOP
 1000 RETURN
C
      ENTRY DARD1(I1,I2,I3,I4,I5,I6)
      I1=IX(1,IDA)
      I2=IX(2,IDA)
      I3=IX(3,IDA)
      I4=IX(4,IDA)
      I5=IX(5,IDA)
      I6=IX(6,IDA)
      RETURN
C
      ENTRY DAWR1(I1,I2,I3,I4,I5,I6)
      MXUSED=MAX0(MXUSED,IDA)
      IX(1,IDA)=I1
      IX(2,IDA)=I2
      IX(3,IDA)=I3
      IX(4,IDA)=I4
      IX(5,IDA)=I5
      IX(6,IDA)=I6
      RETURN
C
      ENTRY DARD2(I1,I2,X1,X2)
      I1=IX(1,IDA)
      I2=IX(2,IDA)
      IF(NIPR.EQ.1) THEN
        IR1=IX(3,IDA)
        IS1=IX(4,IDA)
      ELSE
        IR2(1)=IX(3,IDA)
        IR2(2)=IX(4,IDA)
        IS2(1)=IX(5,IDA)
        IS2(2)=IX(6,IDA)
      ENDIF
      X1=R
      X2=S
      RETURN
C
      ENTRY DAWR2(I1,I2,X1,X2)
      MXUSED=MAX0(MXUSED,IDA)
      IX(1,IDA)=I1
      IX(2,IDA)=I2
      R=X1
      S=X2
      IF(NIPR.EQ.1) THEN
        IX(3,IDA)=IR1
        IX(4,IDA)=IS1
      ELSE
        IX(3,IDA)=IR2(1)
        IX(4,IDA)=IR2(2)
        IX(5,IDA)=IS2(1)
        IX(6,IDA)=IS2(2)
      ENDIF
      RETURN
C
      ENTRY DACLOS
      WRITE(6,610) MXUSED,MAX
  610 FORMAT(/' *** IN-CORE DA SIMULATOR USED',I10,' OF THE',I10,
     1        ' ALLOCATED RECORDS')
      RETURN
      END
      FUNCTION THREEJ (J1,J2,J3)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     COMPUTATION OF SPECIAL WIGNER 3J COEFFICIENT WITH
C     VANISHING PROJECTIONS.  SEE EDMONDS, P. 50.
C
C ***
C *** THIS VERSION EVALUATES BINOM AND PARITY IN-LINE
C ***   SHOULD IMPROVE EFFICIENCY, ESPECIALLY ON CRAY;
C ***   ALSO GIVES IMPROVEMENT ON AMDAHL  (SG: 20 DEC 92)
C ***
C     STATEMENT FUNCTION FOR DELTA ASSOCIATED W/ RACAH AND SIXJ SYMBOLS
C     DELTA(I,J,K)= SQRT(1.D0/ ( BINOM(I+J+K+1,I+J-K) *
C    1                 BINOM(K+K+1,I-J+K) * DBLE(K+J-I+1) )  )
C
      I1=J1+J2+J3
      IF (I1-2*(I1/2).NE.0)    GO TO 8
    1 I2=J1-J2+J3
      IF (I2) 8,2,2
    2 I3=J1+J2-J3
      IF (I3) 8,3,3
    3 I4=-J1+J2+J3
      IF (I4) 8,4,4
    4 I5=I1/2
      I6=I2/2
      SIGN=1.D0
      IF (I5-2*(I5/2).NE.0) SIGN=-SIGN
C   7 THREEJ=SIGN*DELTA(J1,J2,J3)*BINOM(I5,J1)*BINOM(J1,I6)
C     B1,B2 ARE BINOM ASSOCIATED W/ DELTA
      N=J1+J2+J3+1
      M=J1+J2-J3
      NM = N-M
      MNM = MIN(NM,M)
      IF(MNM.LE.0) THEN
        B1=1.D0
      ELSE
        FN = N+1
        F = 0.D0
        B = 1.D0
        DO 101 I = 1,MNM
        F = F+1.D0
        C = (FN-F)*B
  101   B = C/F
        B1 = B
      ENDIF
      N=J3+J3+1
      M=J1-J2+J3
      NM = N-M
      MNM = MIN(NM,M)
      IF(MNM.LE.0) THEN
        B2=1.D0
      ELSE
        FN = N+1
        F = 0.D0
        B = 1.D0
        DO 102 I = 1,MNM
        F = F+1.D0
        C = (FN-F)*B
  102   B = C/F
        B2 = B
      ENDIF
      DELTA=SQRT(1.D0/(B1*B2*(J3+J2-J1+1)))
C     B3=BINOM(I5,J1),  B4=BINOM(J1,I6)
      N=I5
      M=J1
      NM = N-M
      MNM = MIN(NM,M)
      IF(MNM.LE.0) THEN
        B3=1.D0
      ELSE
        FN = N+1
        F = 0.D0
        B = 1.D0
        DO 103 I = 1,MNM
        F = F+1.D0
        C = (FN-F)*B
  103   B = C/F
        B3 = B
      ENDIF
      N=J1
      M=I6
      NM = N-M
      MNM = MIN(NM,M)
      IF(MNM.LE.0) THEN
        B4=1.D0
      ELSE
        FN = N+1
        F = 0.D0
        B = 1.D0
        DO 104 I = 1,MNM
        F = F+1.D0
        C = (FN-F)*B
  104   B = C/F
        B4 = B
      ENDIF
      THREEJ=SIGN*DELTA*B3*B4
      RETURN
    8 THREEJ=0.D0
      RETURN
      END
      FUNCTION PARITY(I)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARITY=1.D0
      IF((I/2)*2-I.NE.0) PARITY=-1.D0
      RETURN
      END
      FUNCTION SIXJ(J1,J2,J5,J4,J3,J6)
C     CALCULATES 6-J SYMBOL:   _(J1 J2 J3 )_
C                               (J4 J5 J6 )
C     INTERFACE TO J6J ROUTINE.
C     MODFIED BY S. GREEN 20 AUG 93; PASS DIMENSION OF XJ6J FOR CHECKING
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(MXDIM=200)
      DIMENSION XJ6J(MXDIM)
      IVAL=MXDIM
      CALL J6J(         DBLE(J2),DBLE(J3),
     1         DBLE(J4),DBLE(J5),DBLE(J6),
     3         IVAL,XJ1MIN,XJ6J)
      IND=1+J1-INT(XJ1MIN+0.1D0)
      SIXJ=0.D0
      IF(IND.GE.1 .AND. IND.LE.IVAL) SIXJ=XJ6J(IND)
      RETURN
      END
      SUBROUTINE J6J(J2,J3,L1,L2,L3,IVAL,J1MIN,D6J)
      IMPLICIT DOUBLE PRECISION (A-H,J-Z)
      DIMENSION D6J(2)
      DATA ZERO/0.D0/,TENTH/0.1D0/,HALF/0.5D0/,ONE/1.D0/,TWO/2.D0/,
     $ CONST/1.0D-12/
      E(J1S)=SQRT((J1S-MJ23S)*(J23P1S-J1S)*(J1S-ML23S)*(L23P1S-J1S))
      F(J1,JJP1)=(TWO*J1+ONE)*(JJP1*(FACT-JJP1-TWO*LLP1)+FACT2)
C
C  THIS ROUTINE CALCULATES THE 6-J COEFFICIENTS FOR ALL PERMISSIBLE
C  VALUES OF J1 FOR FIXED VALUES OF J2, J3, L1, L2, AND L3 USING THE
C  RECURSIVE ALGORITHM OF K. SCHULTEN AND R. G. GORDON, J. MATH. PHYS.
C  VOL. 16, P. 1961, (1975).
C  PROGRAMMED BY D. E. FITZ, 10/22/79
C  MODIFIED BY S. GREEN 20 AUG 93 TO TEST DIMENSION ON D6J
C
      MXDIM=IVAL
      JJP2=J2*(J2+ONE)
      JJP3=J3*(J3+ONE)
      LLP1=L1*(L1+ONE)
      LLP2=L2*(L2+ONE)
      LLP3=L3*(L3+ONE)
      MJ23S=(J2-J3)**2
      ML23S=(L2-L3)**2
      J23P1S=(J2+J3+ONE)**2
      L23P1S=(L2+L3+ONE)**2
      FACT2=(LLP2-LLP3)*(JJP2-JJP3)
      FACT=JJP2+JJP3+LLP2+LLP3
      J1MIN=MAX(ABS(J2-J3),ABS(L2-L3))
      J1MAX=MIN(J2+J3,L2+L3)
      IVAL=INT(J1MAX-J1MIN+ONE+TENTH)
      IF (IVAL.GT.MXDIM) THEN
        WRITE(6,*) ' SIXJ/J6J INADEQUATE STORAGE.  NEEDS',IVAL,
     1             '  GIVEN',MXDIM
        STOP
      ENDIF
C
C  TEST FOR OTHER TRIANGULAR INEQUALITES.
C
      IL1=INT(TWO*L1+TENTH)
      IL2=INT(TWO*L2+TENTH)
      IL3=INT(TWO*L3+TENTH)
      IJ2=INT(TWO*J2+TENTH)
      IJ3=INT(TWO*J3+TENTH)
      IF((IJ2.LE.IL1+IL3.AND.IJ2.GE.IABS(IL1-IL3)).AND.
     $ (IJ3.LE.IL1+IL2.AND.IJ3.GE.IABS(IL1-IL2))) GO TO 11
      DO 12 I=1,IVAL
   12 D6J(I)=ZERO
      RETURN
C
   11 INMID=(IVAL+3)/2
      SGNV=J2+J3+L2+L3
      SGN=ONE
      ISIGN=INT(SGNV+TENTH)
      IF(MOD(ISIGN,2).NE.0) SGN=-ONE
      D6J(1)=HALF
C
C  UPWARD RECURSION.
C
      IF(IVAL.EQ.1) GO TO 40
      JJP1=J1MIN*(J1MIN+ONE)
      F1=F(J1MIN,JJP1)
      J1=J1MIN+ONE
      J1S=J1*J1
      E2=E(J1S)
      IF(J1MIN.LT.TENTH) GO TO 15
      D6J(2)=-F1*D6J(1)/(E2*J1MIN)
      GO TO 16
   15 D6J(2)=-HALF*(LLP2+JJP2-LLP1)*D6J(1)/SQRT(JJP2*LLP2)
   16 SCALE=D6J(2)
      IF(IVAL.EQ.2) GO TO 40
      DO 21 IJ2=3,INMID
      JJP1=J1*(J1+ONE)
      F1=F(J1,JJP1)
      J1=J1+ONE
      E1=E2
      J1S=J1*J1
      E2=E(J1S)
   21 D6J(IJ2)=-(F1*D6J(IJ2-1)+J1*E1*D6J(IJ2-2))/(E2*(J1-ONE))
      SCALE=D6J(INMID)
      IEXC=5
      IF(ABS(SCALE).GT.CONST) GO TO 18
      INMID=INMID-1
      SCALE=D6J(INMID)
      IEXC=3
      GO TO 30
   18 IF(IVAL.EQ.3) GO TO 40
C
C  DOWNWARD RECURSION.
C
   30 D6J(IVAL)=HALF
      J1=J1MAX
      J1S=J1*J1
      JJP1=J1*(J1+ONE)
      F1=F(J1,JJP1)
      E1=E(J1S)
      D6J(IVAL-1)=-F1*D6J(IVAL)/(E1*(J1+ONE))
      IEND=IVAL-INMID
      IF(IVAL.LE.IEXC) GO TO 31
      DO 32 IJ2=2,IEND
      J1=J1-ONE
      E2=E1
      J1S=J1*J1
      JJP1=J1*(J1+ONE)
      E1=E(J1S)
      F1=F(J1,JJP1)
   32 D6J(IVAL-IJ2)=-(J1*E2*D6J(IVAL-IJ2+2)+F1*D6J(IVAL-IJ2+1))/
     $ (E1*(J1+ONE))
C
C  MATCH UPWARD AND DOWNWARD RECURSIVE RESULTS BY SCALING.
C
   31 SCALE=SCALE/D6J(INMID)
      DO 33 IJ2=INMID,IVAL
   33 D6J(IJ2)=SCALE*D6J(IJ2)
C
C  NORMALIZE RESULTS AND SET PHASE.
C
   40 SUM=ZERO
      DO 41 IJ2=1,IVAL
      J1=J1MIN+DBLE(IJ2-1)
   41 SUM=SUM+(TWO*J1+ONE)*D6J(IJ2)**2
      RNORM=ONE/SQRT(SUM*(TWO*L1+ONE))
      IF((SGN*D6J(IVAL)).LT.ZERO) RNORM=-RNORM
      DO 42 IJ2=1,IVAL
   42 D6J(IJ2)=D6J(IJ2)*RNORM
      RETURN
      END
      FUNCTION THRJ(F1,F2,F3,G1,G2,G3)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE MUNG
      DIMENSION X(202),Y(202)
      DATA MUNG/0/,MXIX/202/
      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