C     PROGRAM TO COMPUTE PRESSURE BROADENING  CROSS SECTIONS,
C       SIG(JA,JB;JA1,JB1;K)   <==>  SIG(VA,JA,VB,JB;VA,JA1,VB,JB1;K)
C       FROM S-MATRICES SAVED ON TWO (2) TAPES BY SCATTERING PROG.
C       EACH TAPE HAS S-MATRICES FROM A DIFFERENT VIBRAITONAL LEVEL
C       MOLSCAT, VERSION 11 (UNFORMATTED SAVE TAPE) COMPATIBLE.
C       UPDATED FOR VERSION 14 FORMAT JUL 94
C
C    ------------------------------------------------------------------
C
C     NAMELIST &CNTROL PARAMETERS  ---
C       NLPRBR (EQUIV IFLS), LINE, and LTYPE ARE MOLSCAT &INPUT VARIABLES
C         TO DESCRIBE THE CROSS SECTIONS (SEE MOLSCAT DOCUMENTATION) 
C       IIN(I),I=1,2 ARE INPUT TAPES (DEFAULT IS 10, 11)
C      parameters to control printing:
C       JSTEPZ - IF > 0 CONTRIBUTION FROM EVERY JSTEPZ'TH JTOT ONLY.
C       JMOD   - PRINTS ACCUMULATED PRBR SIGS ONLY FOR MOD(JTOT,JMOD)=JMOD
C       JPRINT - SUPPRESSES PRINTING OF PRBR SIGS FOR JTOT<JPRINT
C      parameters to control 'shifting' vibration from one of the tapes:
C       LSHIFT   - SHIFT INPUT 'LINE' VALUES FROM LSHIFT-TH TAPE (default 2)
C       IVIB(2)  - VIBRATIONAL QUANT. NOS. TO ADD TO JLEV
C       VSHIFT(2)- VIB. ENERGY SHIFT TO ADD TO ELEVEL/ENERGY (default 1000.)
C       IIN(2)   - INPUT DATA SET NOS FOR TWO VIB CALCS (DEFAULT IT1,IT2)
C       LFMT - FORMATTED/UNFORMATTED INPUT TAPE. FALSE (DEFAULT) UNFORM
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                 (REQUIRES MODIFIED SUBROUTINE PRBR)
C     25 JUN 1992 ADDS ABILITY TO READ UNFORMATTED (V11 AND LATER)
C                 SAVE TAPES.  CONTROLLED BY LFMT.
C     1 SEPT 1994 UPDATED TO VERSION 14 
C    ------------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL LFMT
C
C     ALLOW FOR UP TO 200 POSSIBLE ENERGIES (MXNRG),
C               UP TO 400 CHANNELS PER S-MATRIX (MXCH)
C               UP TO 1000 'LEVELS' (MXLEV) DETERMINED BY CMBASE DIMS.
C               UP TO 1000  VALUES IN JLEV  (MXLVQN)
C               UP TO 500 'LINES' (MXLN)
      PARAMETER(MXNRG=200,MXCH=400,MXLVQN=1000,MXLN=500,MXLEV=1000)
C
      CHARACTER*1 CA(2)
      LOGICAL EOF(2)
      INTEGER  PRNTLV,ITYPE,NLEV,NQN,NNRG,NOPEN,JLEV,J,L
      DIMENSION LABEL1(20),LABEL2(20)
      DIMENSION JLEV(MXLVQN),J(MXCH),L(MXCH),NBASIS(MXCH),
     1          ENERGY(MXNRG),WVEC(MXCH),
     2          SREAL(MXCH*MXCH),SIMAG(MXCH*MXCH),
     3          IC(MXCH),IL(MXCH),IC1(MXCH),IL1(MXCH)
C     VARIABLES TO HOLD INPUT FROM TWO TAPES PRIOR TO MERGING.
      DIMENSION ELVL1(MXLEV),ELVL2(MXLEV),
     1          JLEV1(MXLVQN),JLEV2(MXLVQN),ENERG1(MXNRG),ENERG2(MXNRG)
C
C     PRBR VARIABLES
      DIMENSION LINE(4*MXLN),LTYPE(MXLN)
      EQUIVALENCE(NLPRBR,IFLS)
C
C     VIBRATIONAL STATE VARIABLES
      DIMENSION IVIB(2),VSHIFT(2),IIN(2)
C
C  INFORMATION ORIGINALLY PASSED AS ENTRY PRBRBS, NOW IN COMMON
C     VERSION 14 DIMENSIONS FOR /CMBASE/
      LOGICAL MPLMIN
      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     DA READ/WRITE SIMULATOR NEEDS NIPR.  SHOULD BE SET TO NIPR=2
      COMMON /INTPAC/NIPR
C
      NAMELIST/CNTROL/LINE,LTYPE,ILSU,PRNTLV,IFLS,NLPRBR,IFEGEN,
     1         JSTEPZ,JMOD,JPRINT,IVIB,VSHIFT,IIN,LSHIFT,LFMT
C
C
      DATA CA/'A','B'/
C
C     DEFAULT INPUT UNITS ...
      DATA IT1/10/, IT2/11/
C
C *** TOLERANCE FOR ENERGY (AND OTHER REAL VALUE) CHECK
      DATA TOL/1.E-6/
C
C     FORMATS FOR SAVE TAPE . . .
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 
C     CALL XUFLOW(0) -- not required on rs/6000
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     SET PARAMETERS AND DEFAULT VALUES; THEN READ &CNTROL
      NIPR=2
      LFMT=.FALSE.
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.
C
      JSTEPZ=0
      JMOD=1
      JPRINT=0
      PRNTLV=1
      DO 1101 I=1,MXLN
 1101 LTYPE(I)=-1
C     DEFAULT IS THAT 2ND TAPE IS VIBRATIONALLY EXCITED; VSHIFT IS ARBITRARY
      IIN(1)=IT1
      IIN(2)=IT2
      LSHIFT=2
      IVIB(1)=0
      IVIB(2)=1
      VSHIFT(1)=0.
      VSHIFT(2)=1000.
      READ(5,CNTROL,END=9000)
C     FORCE IFEGEN TO ZERO SINCE ENERGIES CANNOT BE ADDED TO TAPE.
      IFEGEN=0
C
      WRITE(6,604)  IIN,IVIB,VSHIFT
  604 FORMAT('  PROGRAM TO CALCULATE VIBRATION-ROTATION PRESSURE',
     1       ' BROADENING CROSS SECTIONS'/'  INPUT FROM TWO SAVE TAPES',
     2       ', ONE FOR EACH OF TWO UNCOUPLED VIBRATIONAL LEVELS'/
     3       '  THIS IS MOLSCAT VERSION 14 (SEPT 94) COMPATIBLE'//
     4       25X,'   LEVEL-1   LEVEL-2'/
     5       ' S-MATRIX FROM INPUT UNIT',I9,I10/
     6       ' VIB QUANTUM NUMBER      ',I9,I10/
     7       ' VIB ENERGY SHIFT (1/CM) ',2F10.2)
C
C      * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C      READ BEGINNING OF TAPE 1                                        *
C
      IT=IIN(1)
      IF (LFMT) THEN
         READ(IT,100) LABEL1,ITYPE1,NLEV1,NQN1,URED1,IP1
      ELSE
         READ(IT) LABEL1,ITYPE1,NLEV1,NQN1,URED1,IP1
      ENDIF
      WRITE(6,600) IT,LABEL1,ITYPE1,NLEV1,NQN1,URED1,IP1
  600 FORMAT(//'  INPUT 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     THERE ARE LIMITS TO BACKWARD COMPATIBILITY ...
      IF (IP1.LE.3) THEN
        WRITE(6,*) ' *** THIS VERSION OF PRBRVIB REQUIRES IPROGM.GT.3'
        STOP
      ENDIF
      NSQ=NLEV1*NQN1
      IF (NSQ.LE.MXLVQN) GO TO 1001
      WRITE(6,619) NSQ,MXLVQN
  619 FORMAT(/' ERROR. NLEV*NQN EXCEEDS MXLVQN',2I6)
      STOP
 1001 IF (LFMT) THEN
         READ(IT,101) (JLEV1(I),I=1,NSQ)
      ELSE
         READ(IT) (JLEV1(I),I=1,NSQ)
      ENDIF
      DO 1002 I=1,NLEV1
 1002 WRITE(6,601) I,(JLEV1(NLEV1*(K-1)+I),K=1,NQN1)
  601 FORMAT( ' FOR LEVEL',I4,'  QUANTUM NOS. ARE',10I4)
      IF (LFMT) THEN
         READ(IT,102) NLVL1,(ELVL1(I),I=1,NLVL1)
      ELSE
         READ(IT) NLVL1,(ELVL1(I),I=1,NLVL1)
      ENDIF
 1003 IF (NLVL1.LE.MXLEV)  GO TO 1005
      WRITE(6,699) NLVL1,MXLEV
  699 FORMAT(/' * * * ERROR. NLEVEL. GT. MXLEV',2I5)
      STOP
 1005 WRITE(6,603) NLVL1,(ELVL1(I),I=1,NLVL1)
  603 FORMAT( /  I5,'   ROTATIONAL LEVELS'/(' ',1P,5D14.5))
      IF (LFMT) THEN
         READ(IT,102) NNRG1,(ENERG1(I),I=1,NNRG1)
      ELSE
         READ(IT) NNRG1,(ENERG1(I),I=1,NNRG1)
      ENDIF
      IF (NNRG1.LE.MXNRG) GO TO 1004
      WRITE(6,618) NNRG1,MXNRG
  618 FORMAT(/' ERROR. NNRG EXCEEDS MXNRG',2I6)
      STOP
 1004 WRITE(6,602) NNRG1,(ENERG1(I),I=1,NNRG1)
  602 FORMAT( /  I5,'  ENERGIES ARE'/(' ',1P,5D14.5))
C
C      * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C      READ BEGINNING OF TAPE 2                                        *
C
      IT=IIN(2)
      IF (LFMT) THEN
         READ(IT,100) LABEL2,ITYPE2,NLEV2,NQN2,URED2,IP2
      ELSE
         READ(IT) LABEL2,ITYPE2,NLEV2,NQN2,URED2,IP2
      ENDIF
      WRITE(6,600) IT,LABEL2,ITYPE2,NLEV2,NQN2,URED2,IP2
      IF (IP2.LE.3) THEN
        WRITE(6,*) ' *** THIS VERSION OF PRBRVIB REQUIRES IPROGM.GT.3'
        STOP
      ENDIF
      NSQ=NLEV2*NQN2
      IF (NSQ.LE.MXLVQN) GO TO 1011
      WRITE(6,619) NSQ,MXLVQN
      STOP
 1011 IF (IP1.LT.14.AND.IP2.LT.14) THEN
        IP=MIN(IP1,IP2)
      ELSEIF (IP1.GE.14.AND.IP2.GE.14) THEN 
        IP=MAX(IP1,IP2)
      ELSE
        WRITE(6,*) ' *** CANNOT PROCESS PRE- AND POST-VERSION 14',
     1             ' TAPES IN SAME JOB'
        STOP
      ENDIF    
      IF (LFMT) THEN
         READ(IT,101) (JLEV2(I),I=1,NSQ)
      ELSE
         READ(IT) (JLEV2(I),I=1,NSQ)
      ENDIF
      DO 1012 I=1,NLEV2
 1012 WRITE(6,601) I,(JLEV2(NLEV2*(K-1)+I),K=1,NQN2)
      IF (LFMT) THEN
         READ(IT,102) NLVL2,(ELVL2(I),I=1,NLVL2)
      ELSE
         READ(IT) NLVL2,(ELVL2(I),I=1,NLVL2)
      ENDIF
 1013 IF (NLVL2.LE.MXLEV)  GO TO 1015
      WRITE(6,699) NLVL2,MXLEV
      STOP
 1015 WRITE(6,603) NLVL2,(ELVL2(I),I=1,NLVL2)
      IF (LFMT) THEN
         READ(IT,102) NNRG2,(ENERG2(I),I=1,NNRG2)
      ELSE
         READ(IT) NNRG2,(ENERG2(I),I=1,NNRG2)
      ENDIF
      IF (NNRG2.LE.MXNRG) GO TO 1014
      WRITE(6,618) NNRG1,MXNRG
      STOP
 1014 WRITE(6,602) NNRG2,(ENERG2(I),I=1,NNRG2)
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     MERGE TWO SETS OF VALUES: JLEVEL/ELEVEL/JLEV/ENERGY
C
      WRITE(6,620)
  620 FORMAT(/ 15(' ****')/21X,'MERGED (AND SHIFTED) DESCRIPTORS'
     1        /15(' ****') )
C     N.B. CODE ONLY SUPPORTS ITYPE=10*N+1 ==> NQN=2
C          CHANGES THIS TO ITYPE=10*N+2 ==> NQN=3
      IF (ITYPE1.EQ.ITYPE2 .AND. NQN1.EQ.NQN2 .AND. NQN1.EQ.2 .AND.
     1    MOD(ITYPE1,10).EQ.1. AND. ABS(URED1-URED2).LE.TOL) GO TO 1200
      WRITE(6,621)
  621 FORMAT(//' ***** TERMINAL ERROR: ITYPE, NQN, OR URED MISMATCHED',
     1         ' OR ILLEGAL.')
      STOP
 1200 NLEV=NLEV1+NLEV2
      NQN=3
      ITYPE=ITYPE1+1
      NSQ=NLEV*NQN
      IF (NSQ.LE.MXLVQN) GO TO 1201
      WRITE(6,619) NSQ,MXLVQN
 1201 DO 1202 I=1,NLEV1
      JLEV(I)=JLEV1(I)
      JLEV(NLEV+I)=IVIB(1)
 1202 JLEV((2*NLEV)+I)=JLEV1(NLEV1+I)
      I2=0
      DO 1203 I=NLEV1+1,NLEV
      I2=I2+1
      JLEV(I)=JLEV2(I2)
      JLEV(NLEV+I)=IVIB(2)
 1203 JLEV((2*NLEV)+I)=JLEV2(NLEV2+I2)+NLEV1
      DO 1204 I=1,NLEV
 1204 WRITE(6,601) I,(JLEV(NLEV*(K-1)+I),K=1,NQN)
      NLEVEL=NLVL1+NLVL2
      IF (NLEVEL.LE.MXLEV)  GO TO 1205
      WRITE(6,699) NLEVEL,MXLEV
      STOP
 1205 DO 1206 I=1,NLVL1
 1206 ELEVEL(I)=ELVL1(I)+VSHIFT(1)
      DO 1207 I=1,NLVL2
 1207 ELEVEL(NLVL1+I)=ELVL2(I)+VSHIFT(2)
      WRITE(6,603) NLEVEL,(ELEVEL(I),I=1,NLEVEL)
      NNRG=NNRG1+NNRG2
      IF (NNRG.LE.MXNRG) GO TO 1211
      WRITE(6,618) NNRG,MXNRG
      STOP
 1211 DO 1212 I=1,NNRG1
 1212 ENERGY(I)=ENERG1(I)+VSHIFT(1)
      DO 1213 I=1,NNRG2
 1213 ENERGY(NNRG1+I)=ENERG2(I)+VSHIFT(2)
      WRITE(6,602) NNRG,(ENERGY(I),I=1,NNRG)
C                                                                     *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     SET UP FOR PRESSURE BROAD.; CALL PRBRIN TO PROCESS &CNTROL DATA  *
C                                                                     *
      ILSU=11
      DO 1220 I=1,MXCH
 1220 NBASIS(I)=I
C     IF LSHIFT=0, GO DIRECTLY TO PRBRIN; OTHERWISE MODIFY LINE()
      IF (LSHIFT.LE.0) GO TO 1230
      IF (LSHIFT.EQ.1 .OR. LSHIFT.EQ.2) THEN
        WRITE(6,622) LSHIFT,CA(LSHIFT),NLEV1
  622   FORMAT(//' INPUT LSHIFT =',I4/' INPUT LINE() VALUES FOR J',A1/
     1      6X,' WILL BE SHIFTED BY',I3,'  TO REQUEST 2ND VIB. LEVEL')
        IF (NLPRBR.LE.0) THEN
          WRITE(6,623)
  623     FORMAT(/' ***** NO LINES REQUESTED BY NLPRBR OR IFLS')
          STOP
        ENDIF
        NL=0
        DO 1221 I=1,NLPRBR
        WRITE(6,624) (LINE(NL+K),K=1,4)
  624   FORMAT(/'  LINE() =',4I4,' ======>')
        LINE(NL+LSHIFT)=LINE(NL+LSHIFT)+NLEV1
        LINE(NL+LSHIFT+2)=LINE(NL+LSHIFT+2)+NLEV1
        WRITE(6,625) (LINE(NL+K),K=1,4)
  625   FORMAT(34X,4I4)
 1221   NL=NL+4
      ELSE
        WRITE(6,626) LSHIFT
  626   FORMAT(//' ***** TERMINAL ERROR.  ILLEGAL LSHIFT =',I6)
        STOP
      ENDIF
C ***
 1230 CALL PARMPR(JSTEPZ,JMOD,JPRINT)
C ***
      CALL PRBRIN(IFLS,LINE,LTYPE,MXLN,ILSU,NNRG,ENERGY,NNRG,IFEGEN,
     1            JLEV,PRNTLV)
C
C      * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C                                                                      *
C     START INPUT OF JTOT'S (ALTERNATE TAPES) AND ENERGIES             *
C                                                                      *
      JSTEPX=-1
      JTOLD=-1
      INFROM=1
      EOF(1)=.FALSE.
      EOF(2)=.FALSE.
      ISVD=0
C
 2000 IT=IIN(INFROM)
      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
      IF (JTOT.LT.0)  THEN
C       PROCESS CHKPT/RESTART MARKERS (FROM EARLY VERISONS OF MOLSCAT)
        WRITE(6,698)  JTOT
  698   FORMAT(/' * * * NOTE. ROLLOUT MARKER ENCOUNTERED:',I6)
        GO TO 2000
      ENDIF
      ECHK=ECHK+VSHIFT(INFROM)
      IF (INFROM.EQ.2) INRG=INRG+NNRG1
C     COUPLED STATES CASES NEED MVALUE; SET ACCORDING FOR MPLMIN=.T.
      MVALUE=M-1
      IF (PRNTLV.GE.3)  WRITE(6,662) INFROM,JTOT,M,INRG,ECHK,IEXCH,WT
  662 FORMAT(' TAPE-',I1,' JTOT=',I4,'.',I2,'  ENERGY(',I3,')=',F10.3,
     1       ' IEX=',I2,' WT=',F7.3)
C
C     SEE IF THIS IS A NEW (OR THE FIRST) JTOT.
      IF (JTOT.EQ.JTOLD.OR.JTOLD.EQ.-1) GO TO 2011
C     NEW JTOT:  DETERMINE/CHECK JSTEP
      JSTEP=JTOT-JTOLD
      IF (JSTEP.EQ.JSTEPX) GO TO 2911
      IF (JSTEPX.EQ.-1) GO TO 2912
      WRITE(6,691) JSTEP,JSTEPX
  691 FORMAT('0 ***** ERROR.  JSTEP INCONSISTENCY.  CURRENT, OLD',2I4)
C     JSTEP=0
 2912 JSTEPX=JSTEP
 2911 IF (INFROM.EQ.2) GO TO 2913
C
C     NEW JTOT ON TAPE-1.
      IF (EOF(2)) GO TO 2913
C     UNLESS WE'VE HIT EOF ON TAPE-2, SAVE PARMS FROM 1 & READ TAPE-2
C     ISVD=0 => WE HAVE NOT YET READ A JTOT RECORD
      IF (ISVD.EQ.0) THEN
        JTSV=JTOT
        INSV=INRG
        ECHSV=ECHK
        IEXSV=IEXCH
        WTSV=WT
        MSV=M
        ISVD=INFROM
        INFROM=2
        GO TO 2000
C     ISVD=2 => WE HAVE ALREADY READ A JTOT RECORD
      ELSEIF (ISVD.EQ.2) THEN
C     IF (ISVD.EQ.2) THEN
        CALL IXCH(JTSV,JTOT)
        CALL IXCH(INSV,INRG)
        CALL IXCH(IEXSV,IEXCH)
        CALL IXCH(MSV,M)
        CALL DXCH(ECHSV,ECHK)
        CALL DXCH(WTSV,WT)
        ISVD=INFROM
        INFROM=2
        IT=IIN(INFROM)
        GO TO 2011
      ELSE
        WRITE(6,628)  INFROM, JTOT
        STOP
      ENDIF
C
C     WE HAVE NOW FINISHED A JTOT FROM BOTH TAPES.
C       PRINT CROSS SECTINS (ON JMOD AND JPRINT CRITERIA)
 2913 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 2920
      IF (JTOLD-JSTEPZ*(JTOLD/JSTEPZ).EQ.0) GO TO 2920
      WRITE(6,627) JSTEPZ
  627 FORMAT('  **** **** **** JTOT SKIPPED BECAUSE OF JSTEPZ =',I4)
      GO TO 2012
 2920 IF (JTOLD-JMOD*(JTOLD/JMOD).NE.0.AND.JTOLD.LE.JPRINT) GO TO 2012
      IF (PRNTLV.GE.2)  CALL PRBOUT(JPB(JSTEP,JSTEPZ))
C
C     SEE IF WE GOT HERE BECAUSE OF EOF(2) ==> INFROM=1
 2012 IF (INFROM.EQ.1) GO TO 2011
C     IF THERE IS NO MORE ON TAPE-1, CONTINUE READING TAPE-2
      IF (EOF(1)) GO TO 2011
C     OTHERWISE, WE SHOULD SAVE TAPE-2 VALUES, RESTORE TAPE-1 VALUES
C       AND CONTINUE INPUT FROM TAPE-1
      IF (ISVD.NE.1) THEN
        WRITE(6,628)  INFROM, JTOT
  628   FORMAT(/' ***** ERROR.  NO SAVED JTOT RECORD FOR TAPE',I4,
     1         ', JTOT =',I4)
        STOP
      ENDIF
      CALL IXCH(JTSV,JTOT)
      CALL IXCH(INSV,INRG)
      CALL IXCH(IEXSV,IEXCH)
      CALL IXCH(MSV,M)
      CALL DXCH(ECHSV,ECHK)
      CALL DXCH(WTSV,WT)
      ISVD=INFROM
      INFROM=1
      IT=IIN(INFROM)
C
 2011 JTOLD=JTOT
C     BELOW CHECKS THAT ENERGY(INRG) CORRESPONDS WITH HEADER RECORD. . .
      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
C     CONTINUE WITH INPUT OF CHANNEL BASIS AND S-MATRICES
 2002 IF (LFMT) THEN
         READ(IT,104,END=9009) NOPEN,(J(I),L(I),WVEC(I),I=1,NOPEN)
      ELSE
        IF (IP.LT.14) THEN
           READ(IT,END=9009) NOPEN,(J(I),L(I),WVEC(I),I=1,NOPEN)
        ELSE
           READ(IT,END=9009) (J(I),L(I),WVEC(I),I=1,NOPEN)
        ENDIF
      ENDIF
      IF (NOPEN.LE.MXCH)  GO TO 2003
      WRITE(6,696)  MXCH,INFROM,JTOT,INRG
  696 FORMAT(/' ***** ERROR. NO. OPEN CHANNELS EXCEEDS MXCH =',I4,
     1  ' FOR TAPE,JTOT,INRG =',3I5)
 2003 NSQ=NOPEN*NOPEN
      IF (INFROM.EQ.1) GO TO 2005
C     FOR 2ND TAPE, SHIFT LEVEL VALUES, J(), BY NLEV1
      DO 2004 I=1,NOPEN
 2004 J(I)=J(I)+NLEV1
C
 2005 IF (LFMT) THEN
         READ(IT,105) (SREAL(I),I=1,NSQ)
         READ(IT,105) (SIMAG(I),I=1,NSQ)
      ELSE
         CALL SREAD(IT,NOPEN,SREAL,IEND)
         CALL SREAD(IT,NOPEN,SIMAG,IEND)
         IF (IEND.NE.0) THEN
            WRITE(6,644) IT
  644       FORMAT(' *****  UNEXPECTED EOF FOUND BY SREAD ON FILE'
     1             ,I4)
            GO TO 9000
         ENDIF
      ENDIF
C
C     SHOULD WE SKIP CALLING PRBR FROM THESE OWING TO JSTEPZ ?
      IF (JSTEPZ.LE.1) GO TO 4000
      IF (JTOT-JSTEPZ*(JTOT/JSTEPZ).EQ.0) GO TO 4000
      IF (PRNTLV.GE.2) 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
 4000 CONTINUE
      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 ANOTHER JTOT,INRG,MVAL SET ...
      GO TO 2000
C
C    UNEXPECTED EOF READING NOPEN, ... RECORD
 9009 WRITE(6,645) IT
  645 FORMAT(' ***** UNEXPECTED EOF READING LEV,L,WV  ON TAPE',I4)
C
C     END OF FILE ON UNIT IIN(INFROM)
 9000 WRITE(6,695)  INFROM,IIN(INFROM)
  695 FORMAT(/' *****  END OF FILE TAPE',I2,', ON UNIT(',I3,').')
      EOF(INFROM)=.TRUE.
C     SEE IF WE'VE REACHED EOF ON OTHER TAPE ALSO.
      IF (EOF(3-INFROM)) GO TO 2006
      IF (ISVD.EQ.3-INFROM) THEN
        CALL IXCH(JTSV,JTOT)
        CALL IXCH(INSV,INRG)
        CALL IXCH(IEXSV,IEXCH)
        CALL IXCH(MSV,M)
        CALL DXCH(ECHSV,ECHK)
        CALL DXCH(WTSV,WT)
        ISVD=INFROM
        INFROM=3-INFROM
        IT=IIN(INFROM)
        GO TO 2011
      ELSE
        WRITE(6,692) INFROM,ISVD
  692   FORMAT(/' ***** ERROR.  INFROM,ISVD =',2I6)
      ENDIF
        INFROM=3-INFROM
        GO TO 2000
C
 2006 WRITE(6,617) JTOT
      CALL PRBOUT(JPB(JSTEP,JSTEPZ))
      CALL DACLOS
C
      STOP
      END
      SUBROUTINE CHKSTR
C     DUMMY VERSION FOR PRBRVIB.V14
C     PRBR GETS TEMP STORAGE VIA /MEMORY/ FOR PRBR3; NOT IMPLEMENTED HERE
C     THIS ROUTINE WILL TRAP AND STOP; ITYPE=3 ONLY
      WRITE(6,*) ' *** CHKSTR CALLED.  NOT USED BY PRBRVIB'
      WRITE(6,*) ' ***                 CANNOT HANDLE ITYPE=3'
      STOP
      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 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
      SUBROUTINE IXCH(I,J)
      IX=I
      I=J
      J=IX
      RETURN
      END
      SUBROUTINE DXCH(A,B)
      DOUBLE PRECISION A,B,X
      X=A
      A=B
      B=X
      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
      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)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XJ6J(200)
      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
      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)
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
      SUBROUTINE DASIZE(ILSU,MXREC)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE MXUSED
      PARAMETER (NREC=250000)
      DIMENSION IX(6,NREC)
      DIMENSION IR2(2),IS2(2)
      EQUIVALENCE (R,IR1,IR2(1)),(S,IS1,IS2(1))
      COMMON /DA/IX
      COMMON /ASSVAR/IDA
      COMMON /INTPAC/NIPR
      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
      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