C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  1 JUN 05 - TJP - DSPDFV,DIRFCK: FIX, INCLUDE UHF,ROHF,GVB PARTS
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL, ROOT
C  7 SEP 04 - MWS - DIRFCKV: CHANGE INTEGRAL UNPACKING LOOP
C  7 SEP 04 - TJP - DIRFCKV: CHANGE DIMENSIONS IN DIRFCKV BY ONE
C  4 NOV 03 - TJP - ADD DIRFCKV,GENRAL,FORMSV,XYZINTV,RT123V
C 26 MAR 03 - KRG - TFTRI,MTARBR: USE DSYMM INSTEAD OF DGEMM
C 26 MAR 02 - CMC - CRAY COMPILER BUG REQUIRES DIRECTIVES IN JKXYZV
C 12 NOV 98 - MWS - REMOVE *ALL, *AMD, *ARD, *CEL, *CVX, *VAX SUPPORT
C 13 APR 98 - AR,AR,AJL - FUJITSU USES STANDARD DGEMM MATRIX OPERATIONS
C 14 AUG 95 - JF  - CHANGES FOR NEC SX3 VECTORIZATION
C 26 JUL 95 - MWS - REMOVE DREAL
C 11 NOV 94 - MWS - REMOVE SOME FTNCHEK WARNINGS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C 21 JUL 94 - MWS - ADD DSPDFV,JDXYZV,JKGNMV,JKXYZV; DELETE SPDD,TRNRK
C  1 JUN 94 - MWS - DELETE TRNSIM
C 11 MAR 94 - FM  - HSTARU: CRAY VECTOR LENGTH FROM 64 TO 128
C  5 APR 94 - MWS - ADD A VARIABLE TO MACHSW, INCLUDE MRTRBR
C 15 JUN 92 - MWS - DELETE SUBROUTINE MTATBS
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C 10 JAN 92 - TLW - CHANGE REWINDS TO CALL SEQREW
C  8 JAN 92 - TLW - MAKE WRITES PARALLEL; ADD COMMON PAR
C  8 JAN 92 - CWS - VAX SUPPORT: USE DGEMM IF DXML PRESENT, OTHERWISE
C                   USE MXMULT, INCLUDE VECTORIZED HSTARU ROUTINE
C 28 AUG 91 - RHN - ADDED FUJITSU (*FUJ) SUPPORT
C  9 AUG 91 - TLW - REDIMENSION FORMS FOR F AND G FUCTIONS
C 14 JUL 91 - MWS - GLDIAG: JACOBI NEEDS STORAGE FOR ALL VECTORS
C 29 JAN 91 - TLW - GLDIAG: CHANGED COMMON DIAGSW TO MACHSW AND
C                   ADDED ICORFL
C 13 NOV 90 - FJ  - FIX AMDAHL VERSION IN TFSQB,TFTRI
C  2 NOV 90 - MWS - CHANGE /INDEX/ TO /INTDEX/
C 23 MAR 90 - MWS - GLDIAG CALLS JACDIA IF KDIAG=3
C  5 DEC 89 - MWS - IBM COMPILER BUG REQUIRES DIRECTIVES IN SPDD.
C 21 OCT 89 - DAVID MULLALLY AND MWS - FORMS,SPDD,JDDSPD: VECTORIZE
C                   ON IJ WHEN LONGER THAN KL, REMOVE DIRECTIVES.
C  8 AUG 89 - MWS - DELETE VECTORIZED SLVLEQ ROUTINE
C 20 JUL 89 - TLW - ADDED *ARD LINES, AND ROUTINE MXMULT
C  5 JUL 89 - FJ  - ADDED AMDAHL (*AMD) SUPPORT (SSL II AND VOCL)
C 22 MAY 89 - MWS - ADD FPS MODEL 500 (*CEL) SUPPORT (VAST-2 LIBRARY)
C 21 FEB 89 - STE - CRAY DIRECTIVES IN FORMS AND SPDD
C 17 FEB 89 - MWS - CHANGE GLDIAG TO USE IBM'S ESSL VERSION 3.
C 19 JAN 89 - MWS - ADD VECTORIZED JDDSPD
C 11 OCT 88 - MWS - PASS WORKING STORAGE TO TRNRK AND TRNSIM
C 12 SEP 88 - JAM - ADD CONVEX (*CVX) AND ALLIANT (*ALL) SUPPORT
C  8 JUN 88 - MWS - ADD VECTORIZED TRNSIM AND TRNRK ROUTINES
C  1 JUN 88 - MWS - ADD VECTORIZED SPDD ROUTINE
C 22 MAY 88 - MWS - ADD VECTORIZED FORMS ROUTINE
C 11 MAY 88 - MWS AND RICHARD WALSH - ADD CRAY SUPPORT
C  7 APR 88 - MWS - NEW MODULE CREATED TO SUPPORT VECTOR MACHINES.
C*MODULE VECTOR
C  ---------------------------------------------------------------
C  THIS CODE SHOULD NOT BE COMPILED ON SCALAR MACHINES.
C  ---------------------------------------------------------------
C  THE SUBROUTINES IN THIS MODULE HAVE THE **SAME** NAMES AS OTHER
C  SUBROUTINES IN OTHER GAMESS MODULES, SHOWN ON THE *DECK CARDS.
C  THESE VECTORIZED ROUTINES ARE INTENDED TO REPLACE THEIR SCALAR
C  COUTERPARTS ON VECTOR MACHINES.  THIS IS DONE BY LOADING THESE
C  ROUTINES **BEFORE** THE OBJECT CODE FOR THEIR SCALAR VERSIONS.
C  ---------------------------------------------------------------
C  SUPPORTED MACHINES ARE
C      *CRY - CRAY RESEARCH
C      *FUJ - FUJITSU VECTOR PROCESSORS
C      *IBM - IBM 3090 SERIES WITH VECTOR FACILITY
C      *NEC - NEC SX-3 SUPER-UX
C  ---------------------------------------------------------------
C  SCALAR MACHINES ABORT IN CASE THIS IS ACCIDENTALLY COMPILED.
C  ---------------------------------------------------------------
C  THERE USED TO BE DIRECTIVES LIKE THIS
C         CDIR$ IVDEP
C         C$DIR NO_RECURRENCE
C         CVD$  NODEPCHK
C         *VOCL LOOP,NOVREC
C  IN FRONT OF THE SCALAR LOOPS IN FORMS,SPDD,JDDSPD.
C  IN CASE I WANT THEM AGAIN, WE KEEP THESE COMMENTS HERE.
C  CDIR$  DIRECTIVES ARE FOR CRAY MACHINES
C  C$DIR  DIRECTIVES ARE FOR CONVEX
C  CVD$   DIRECTIVES ARE FOR THE VAST PREPROCESSOR
C         (ALLIANT, FPS 500)
C  *VOCL  DIRECTIVES ARE FOR AMDAHL/FUJITSU
C  CDEC$  DIRECTIVES ARE FOR VAX (ONLY INHIBITOR DIRECTIVES USED)
C  ---------------------------------------------------------------
C*MODULE VECTOR  *DECK DSPDFV       (GRD2C)
      SUBROUTINE DSPDFV(NG,NR,IJKLG,GIJKL,FIJKL,SIJKL,XYZ,
     1 FIXYZ,FJXYZ,FKXYZ,FLXYZ,SIIXYZ,SJJXYZ,SKKXYZ,SLLXYZ,
     2 SIJXYZ,SIKXYZ,SILXYZ,SJKXYZ,SJLXYZ,SKLXYZ,
     3 XY,XZ,YZ,X,Y,Z,SJ,SK,SL,DIJSI,DIJSJ,DKLSK,DKLSL,DAB)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL IS,JS,KS,LS
      LOGICAL IJS,IJKS,IJKLS
      LOGICAL IIEQJJ,KKEQLL,IJEQKL,IJGTKL,IJLTKL
      LOGICAL SPI,SPJ,SPK,SPL,SPIJ,SPKL,SPIJKL
      LOGICAL SKIPI,SKIPJ,SKIPK,SKIPL
      LOGICAL OUT,DBG
      COMMON /DSHLNO/ LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJ,IJD,KLD,IJ,KL
      COMMON /DSHLT / RTOL,DTOL,VTOL1,VTOL2,VTOLS,OUT,DBG
      COMMON/DERSKP/IIAT,JJAT,KKAT,LLAT,SKIPI,SKIPJ,SKIPK,SKIPL
      COMMON/DERPAR/IDER,JDER,KDER,LDER,NDER
      COMMON/DERSHL/FD(12),SD(78)
      COMMON/SHLEQU/IIEQJJ,KKEQLL,IJEQKL,IJGTKL,IJLTKL
      COMMON/SHLTYP/SPI,SPJ,SPK,SPL,SPIJ,SPKL,SPIJKL
      DIMENSION GIJKL(1),FIJKL(12,1),SIJKL(78,1),IJKLG(4,1)
      DIMENSION DAB(1)
      DIMENSION    XYZ(NG*NR,1)
      DIMENSION  FIXYZ(NG*NR,1), FJXYZ(NG*NR,1)
      DIMENSION  FKXYZ(NG*NR,1), FLXYZ(NG*NR,1)
      DIMENSION SIIXYZ(NG*NR,1),SJJXYZ(NG*NR,1)
      DIMENSION SKKXYZ(NG*NR,1),SLLXYZ(NG*NR,1)
      DIMENSION SIJXYZ(NG*NR,1),SIKXYZ(NG*NR,1)
      DIMENSION SILXYZ(NG*NR,1),SJKXYZ(NG*NR,1)
      DIMENSION SJLXYZ(NG*NR,1),SKLXYZ(NG*NR,1)
      DIMENSION XY(NG*NR),XZ(NG*NR),YZ(NG*NR)
      DIMENSION  X(NG*NR), Y(NG*NR), Z(NG*NR)
      DIMENSION SJ(NG*NR),SK(NG*NR),SL(NG*NR)
      DIMENSION DIJSI(NG*NR),DIJSJ(NG*NR),DKLSK(NG*NR),DKLSL(NG*NR)
C
      DIMENSION FDUM((MAXI+1)*(MAXJ+1)*(MAXK+1)*(MAXL+1)*2)
      DIMENSION FDTMP((MAXI+1)*(MAXJ+1)*(MAXK+1)*(MAXL+1)*2,12)
      DIMENSION TMPSL((MAXI+1)*(MAXJ+1)*(MAXK+1)*(MAXL+1)*2,NG*NR)
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- FIND IJKLNMAX -----
C     LET'S COMPUTE IJKLNMAX AHEAD OF TIME SO WE CAN REMOVE THE
C     VECTOR DEPENDENCE AND COLLAPSE QUAD-LOOPS INTO SINGLE LOOPS.
C     (ALL STORES IN THIS ROUTINE ARE DEPENDENT ONLY ON IJKLN)
C
      IJKLN=0
      DO I=MINI,MAXI
        JMAX=MAXJ
        IF(IIEQJJ) JMAX=I
        DO J=MINJ,JMAX
          KMAX=MAXK
          IF(IJEQKL) KMAX=I
          DO K=MINK,KMAX
            LMAX=MAXL
            IF(KKEQLL           ) LMAX=K
            IF(IJEQKL.AND.K.EQ.I) LMAX=J
            DO L=MINL,LMAX
              IJKLN=IJKLN+1
            END DO
          END DO
        END DO
      END DO
      IJKLNMAX=IJKLN
C
      IF(SPIJKL) GO TO 1000
C
C     ----- NO SHARED EXPONENTS ; SUM UP ( IX * IY * IZ ) -----
C
      IF(NDER.LT.1) GO TO 700
C
C     ----- GRADIENT : COMPUTE VALUES -----
C
      IJKLN=0
      DUMFD1 =ZERO
      DUMFD2 =ZERO
      DUMFD3 =ZERO
      DUMFD4 =ZERO
      DUMFD5 =ZERO
      DUMFD6 =ZERO
      DUMFD7 =ZERO
      DUMFD8 =ZERO
      DUMFD9 =ZERO
      DUMFD10 =ZERO
      DUMFD11 =ZERO
      DUMFD12 =ZERO
      DO 590 IJKLN=1,IJKLNMAX
        NN=IJKLG(1,IJKLN)
        NX=IJKLG(2,IJKLN)
        NY=IJKLG(3,IJKLN)
        NZ=IJKLG(4,IJKLN)
        DUMFXI =ZERO
        DUMFYI =ZERO
        DUMFZI =ZERO
        DUMFXJ =ZERO
        DUMFYJ =ZERO
        DUMFZJ =ZERO
        DUMFXK =ZERO
        DUMFYK =ZERO
        DUMFZK =ZERO
        DUMFXL =ZERO
        DUMFYL =ZERO
        DUMFZL =ZERO
        DO 540 IGR=1,NG*NR
          DUMXY=XYZ(IGR,NX)*XYZ(IGR,NY)
          DUMXZ=XYZ(IGR,NX)*XYZ(IGR,NZ)
          DUMYZ=XYZ(IGR,NY)*XYZ(IGR,NZ)
          IF(SKIPI) GO TO 500
            DUMFXI =DUMFXI + FIXYZ(IGR,NX)*DUMYZ
            DUMFYI =DUMFYI + FIXYZ(IGR,NY)*DUMXZ
            DUMFZI =DUMFZI + FIXYZ(IGR,NZ)*DUMXY
  500     IF(SKIPJ) GO TO 510
            DUMFXJ =DUMFXJ + FJXYZ(IGR,NX)*DUMYZ
            DUMFYJ =DUMFYJ + FJXYZ(IGR,NY)*DUMXZ
            DUMFZJ =DUMFZJ + FJXYZ(IGR,NZ)*DUMXY
  510     IF(SKIPK) GO TO 520
            DUMFXK =DUMFXK + FKXYZ(IGR,NX)*DUMYZ
            DUMFYK =DUMFYK + FKXYZ(IGR,NY)*DUMXZ
            DUMFZK =DUMFZK + FKXYZ(IGR,NZ)*DUMXY
  520     IF(SKIPL) GO TO 530
            DUMFXL =DUMFXL + FLXYZ(IGR,NX)*DUMYZ
            DUMFYL =DUMFYL + FLXYZ(IGR,NY)*DUMXZ
            DUMFZL =DUMFZL + FLXYZ(IGR,NZ)*DUMXY
  530     CONTINUE
  540   CONTINUE
        IF(SKIPI) GO TO 550
          DUMFD1=DUMFD1+DAB(NN)*DUMFXI
          DUMFD2=DUMFD2+DAB(NN)*DUMFYI
          DUMFD3=DUMFD3+DAB(NN)*DUMFZI
  550   IF(SKIPJ) GO TO 560
          DUMFD4=DUMFD4+DAB(NN)*DUMFXJ
          DUMFD5=DUMFD5+DAB(NN)*DUMFYJ
          DUMFD6=DUMFD6+DAB(NN)*DUMFZJ
  560   IF(SKIPK) GO TO 570
          DUMFD7=DUMFD7+DAB(NN)*DUMFXK
          DUMFD8=DUMFD8+DAB(NN)*DUMFYK
          DUMFD9=DUMFD9+DAB(NN)*DUMFZK
  570   IF(SKIPL) GO TO 580
          DUMFD10=DUMFD10+DAB(NN)*DUMFXL
          DUMFD11=DUMFD11+DAB(NN)*DUMFYL
          DUMFD12=DUMFD12+DAB(NN)*DUMFZL
  580   CONTINUE
  590 CONTINUE
C
C     ----- GRADIENT:STORE VALUES -----
C       MOVE THE STORES OUT OF THE LOOP TO PREVENT
C       LIMITATIONS BY CRITICAL VALUES
C
      IF(SKIPI) GO TO 600
        FD( 1)=FD(1)+DUMFD1
        FD( 2)=FD(2)+DUMFD2
        FD( 3)=FD(3)+DUMFD3
  600 IF(SKIPJ) GO TO 620
        FD( 4)=FD(4)+DUMFD4
        FD( 5)=FD(5)+DUMFD5
        FD( 6)=FD(6)+DUMFD6
  620 IF(SKIPK) GO TO 640
        FD( 7)=FD(7)+DUMFD7
        FD( 8)=FD(8)+DUMFD8
        FD( 9)=FD(9)+DUMFD9
  640 IF(SKIPL) GO TO 660
        FD(10)=FD(10)+DUMFD10
        FD(11)=FD(11)+DUMFD11
        FD(12)=FD(12)+DUMFD12
  660 CONTINUE
C
C     ----- ZEROTH AND FIRST DERIVATIVE INTEGRALS -----
C             CALCULATE VALUES
C
  700 IF(.NOT.OUT) RETURN
C
      DO 940 IJKLN=1,IJKLNMAX
        NN=IJKLG(1,IJKLN)
        NX=IJKLG(2,IJKLN)
        NY=IJKLG(3,IJKLN)
        NZ=IJKLG(4,IJKLN)
        DUM=ZERO
        DUMFXI =ZERO
        DUMFYI =ZERO
        DUMFZI =ZERO
        DUMFXJ =ZERO
        DUMFYJ =ZERO
        DUMFZJ =ZERO
        DUMFXK =ZERO
        DUMFYK =ZERO
        DUMFZK =ZERO
        DUMFXL =ZERO
        DUMFYL =ZERO
        DUMFZL =ZERO
        DO 890 IGR=1,NG*NR
          DUM=DUM+XYZ(IGR,NX)*XYZ(IGR,NY)*XYZ(IGR,NZ)
          DUMXY=XYZ(IGR,NX)*XYZ(IGR,NY)
          DUMXZ=XYZ(IGR,NX)*XYZ(IGR,NZ)
          DUMYZ=XYZ(IGR,NY)*XYZ(IGR,NZ)
          IF(SKIPI) GO TO 850
            DUMFXI =DUMFXI + FIXYZ(IGR,NX)*DUMYZ
            DUMFYI =DUMFYI + FIXYZ(IGR,NY)*DUMXZ
            DUMFZI =DUMFZI + FIXYZ(IGR,NZ)*DUMXY
  850     IF(SKIPJ) GO TO 860
            DUMFXJ =DUMFXJ + FJXYZ(IGR,NX)*DUMYZ
            DUMFYJ =DUMFYJ + FJXYZ(IGR,NY)*DUMXZ
            DUMFZJ =DUMFZJ + FJXYZ(IGR,NZ)*DUMXY
  860     IF(SKIPK) GO TO 870
            DUMFXK =DUMFXK + FKXYZ(IGR,NX)*DUMYZ
            DUMFYK =DUMFYK + FKXYZ(IGR,NY)*DUMXZ
            DUMFZK =DUMFZK + FKXYZ(IGR,NZ)*DUMXY
  870     IF(SKIPL) GO TO 880
            DUMFXL =DUMFXL + FLXYZ(IGR,NX)*DUMYZ
            DUMFYL =DUMFYL + FLXYZ(IGR,NY)*DUMXZ
            DUMFZL =DUMFZL + FLXYZ(IGR,NZ)*DUMXY
  880     CONTINUE
  890   CONTINUE
        FDUM(IJKLN)=DUM
        IF(SKIPI) GO TO 900
          FDTMP(IJKLN,1)=DUMFXI
          FDTMP(IJKLN,2)=DUMFYI
          FDTMP(IJKLN,3)=DUMFZI
  900   IF(SKIPJ) GO TO 910
          FDTMP(IJKLN,4)=DUMFXJ
          FDTMP(IJKLN,5)=DUMFYJ
          FDTMP(IJKLN,6)=DUMFZJ
  910   IF(SKIPK) GO TO 920
          FDTMP(IJKLN,7)=DUMFXK
          FDTMP(IJKLN,8)=DUMFYK
          FDTMP(IJKLN,9)=DUMFZK
  920   IF(SKIPL) GO TO 930
          FDTMP(IJKLN,10)=DUMFXL
          FDTMP(IJKLN,11)=DUMFYL
          FDTMP(IJKLN,12)=DUMFZL
  930   CONTINUE
  940 CONTINUE
C
C     ----- ZEROTH AND FIRST DERIVATIVE INTEGRALS -----
C     ----- STORE VALUES -----
C
      DO IJKLN=1,IJKLNMAX
        NN=IJKLG(1,IJKLN)
        GIJKL(NN)=GIJKL(NN)+FDUM(IJKLN)
        IF(SKIPI) GO TO 950
          FIJKL( 1,NN)=FIJKL( 1,NN)+FDTMP(IJKLN,1)
          FIJKL( 2,NN)=FIJKL( 2,NN)+FDTMP(IJKLN,2)
          FIJKL( 3,NN)=FIJKL( 3,NN)+FDTMP(IJKLN,3)
  950   IF(SKIPJ) GO TO 960
          FIJKL( 4,NN)=FIJKL( 4,NN)+FDTMP(IJKLN,4)
          FIJKL( 5,NN)=FIJKL( 5,NN)+FDTMP(IJKLN,5)
          FIJKL( 6,NN)=FIJKL( 6,NN)+FDTMP(IJKLN,6)
  960   IF(SKIPK) GO TO 970
          FIJKL( 7,NN)=FIJKL( 7,NN)+FDTMP(IJKLN,7)
          FIJKL( 8,NN)=FIJKL( 8,NN)+FDTMP(IJKLN,8)
          FIJKL( 9,NN)=FIJKL( 9,NN)+FDTMP(IJKLN,9)
  970   IF(SKIPL) GO TO 980
          FIJKL(10,NN)=FIJKL(10,NN)+FDTMP(IJKLN,10)
          FIJKL(11,NN)=FIJKL(11,NN)+FDTMP(IJKLN,11)
          FIJKL(12,NN)=FIJKL(12,NN)+FDTMP(IJKLN,12)
  980   CONTINUE
      END DO
      RETURN
C
 1000 CONTINUE
C
C     ----- SHARED EXPONENTS ; FORM ( IX * IY * IZ ) -----
C
C     ----- COMPUTE SL(IGR) AHEAD OF TIME -----
C
      IJKLN=0
      DO I=MINI,MAXI
        IS=SPI.AND.I.EQ.1
        JMAX=MAXJ
        IF(IIEQJJ) JMAX=I
        KMAX=MAXK
        IF(IJEQKL) KMAX=I
        DO J=MINJ,JMAX
          JS=SPJ.AND.J.EQ.1
          IJS=IS.OR.JS
          IF(IJS) THEN
            DO IGR=1,NG*NR
              IF(JS.AND.IS) THEN
                SJ(IGR)=DIJSJ(IGR)*DIJSI(IGR)
              ELSE IF(JS) THEN
                SJ(IGR)=DIJSJ(IGR)
              ELSE IF(IS) THEN
                SJ(IGR)=DIJSI(IGR)
              ENDIF
            END DO
          ENDIF
          DO K=MINK,KMAX
            KS=SPK.AND.K.EQ.1
            IJKS=IJS.OR.KS
            IF(IJKS) THEN
              DO IGR=1,NG*NR
                IF(KS.AND.IJS) THEN
                  SK(IGR)=DKLSK(IGR)*SJ(IGR)
                ELSE IF(KS) THEN
                  SK(IGR)=DKLSK(IGR)
                ELSE IF(IJS) THEN
                  SK(IGR)=SJ(IGR)
                ENDIF
              END DO
            ENDIF
            LMAX=MAXL
            IF(KKEQLL           ) LMAX=K
            IF(IJEQKL.AND.K.EQ.I) LMAX=J
            DO L=MINL,LMAX
              IJKLN=IJKLN+1
              LS=SPL.AND.L.EQ.1
              IJKLS=IJKS.OR.LS
              IF (.NOT.IJKLS) THEN
                DO IGR=1,NG*NR
                  TMPSL(IJKLN,IGR)=1
                END DO
              ELSE
                DO IGR=1,NG*NR
                  IF(LS.AND.IJKS) THEN
                    SL(IGR)=DKLSL(IGR)*SK(IGR)
                  ELSE IF(LS) THEN
                    SL(IGR)=DKLSL(IGR)
                  ELSE IF(IJKS) THEN
                    SL(IGR)=SK(IGR)
                  ENDIF
                  TMPSL(IJKLN,IGR)=SL(IGR)
                END DO
              ENDIF
            END DO  ! END OF L LOOP
          END DO    ! END OF K LOOP
        END DO      ! END OF J LOOP
      END DO        ! END OF I LOOP
C
      IF(NDER.LT.1) GO TO 3000
C
C
C     ----- GRADIENT : COMPUTE VALUES -----
C
      DUMFD1 =ZERO
      DUMFD2 =ZERO
      DUMFD3 =ZERO
      DUMFD4 =ZERO
      DUMFD5 =ZERO
      DUMFD6 =ZERO
      DUMFD7 =ZERO
      DUMFD8 =ZERO
      DUMFD9 =ZERO
      DUMFD10 =ZERO
      DUMFD11 =ZERO
      DUMFD12 =ZERO
      DO IGR=1,NG*NR
        DUMFXI =ZERO
        DUMFYI =ZERO
        DUMFZI =ZERO
        DUMFXJ =ZERO
        DUMFYJ =ZERO
        DUMFZJ =ZERO
        DUMFXK =ZERO
        DUMFYK =ZERO
        DUMFZK =ZERO
        DUMFXL =ZERO
        DUMFYL =ZERO
        DUMFZL =ZERO
        DO IJKLN=1,IJKLNMAX
          NN=IJKLG(1,IJKLN)
          NX=IJKLG(2,IJKLN)
          NY=IJKLG(3,IJKLN)
          NZ=IJKLG(4,IJKLN)
          DABNN = DAB(NN)
          DUMXY=DABNN*XYZ(IGR,NX)*XYZ(IGR,NY)*TMPSL(IJKLN,IGR)
          DUMXZ=DABNN*XYZ(IGR,NX)*XYZ(IGR,NZ)*TMPSL(IJKLN,IGR)
          DUMYZ=DABNN*XYZ(IGR,NY)*XYZ(IGR,NZ)*TMPSL(IJKLN,IGR)
          IF(.NOT. SKIPI) THEN
            DUMFXI =DUMFXI + FIXYZ(IGR,NX)*DUMYZ
            DUMFYI =DUMFYI + FIXYZ(IGR,NY)*DUMXZ
            DUMFZI =DUMFZI + FIXYZ(IGR,NZ)*DUMXY
          END IF
          IF(.NOT. SKIPJ) THEN
            DUMFXJ =DUMFXJ + FJXYZ(IGR,NX)*DUMYZ
            DUMFYJ =DUMFYJ + FJXYZ(IGR,NY)*DUMXZ
            DUMFZJ =DUMFZJ + FJXYZ(IGR,NZ)*DUMXY
          END IF
          IF(.NOT. SKIPK) THEN
            DUMFXK =DUMFXK + FKXYZ(IGR,NX)*DUMYZ
            DUMFYK =DUMFYK + FKXYZ(IGR,NY)*DUMXZ
            DUMFZK =DUMFZK + FKXYZ(IGR,NZ)*DUMXY
          END IF
          IF(.NOT. SKIPL) THEN
            DUMFXL =DUMFXL + FLXYZ(IGR,NX)*DUMYZ
            DUMFYL =DUMFYL + FLXYZ(IGR,NY)*DUMXZ
            DUMFZL =DUMFZL + FLXYZ(IGR,NZ)*DUMXY
          END IF
        END DO
        IF(.NOT. SKIPI) THEN
          DUMFD1=DUMFD1 + DUMFXI
          DUMFD2=DUMFD2 + DUMFYI
          DUMFD3=DUMFD3 + DUMFZI
        END IF
        IF(.NOT. SKIPJ) THEN
          DUMFD4=DUMFD4 + DUMFXJ
          DUMFD5=DUMFD5 + DUMFYJ
          DUMFD6=DUMFD6 + DUMFZJ
        END IF
        IF(.NOT. SKIPK) THEN
          DUMFD7=DUMFD7 + DUMFXK
          DUMFD8=DUMFD8 + DUMFYK
          DUMFD9=DUMFD9 + DUMFZK
        END IF
        IF(.NOT. SKIPL) THEN
          DUMFD10=DUMFD10 + DUMFXL
          DUMFD11=DUMFD11 + DUMFYL
          DUMFD12=DUMFD12 + DUMFZL
        END IF
      END DO
C
C     ----- GRADIENT:STORE VALUES -----
C
      IF(.NOT. SKIPI) THEN
        FD( 1) = FD( 1) + DUMFD1
        FD( 2) = FD( 2) + DUMFD2
        FD( 3) = FD( 3) + DUMFD3
      END IF
      IF(.NOT. SKIPJ) THEN
        FD( 4) = FD( 4) + DUMFD4
        FD( 5) = FD( 5) + DUMFD5
        FD( 6) = FD( 6) + DUMFD6
      END IF
      IF(.NOT. SKIPK) THEN
        FD( 7) = FD( 7) + DUMFD7
        FD( 8) = FD( 8) + DUMFD8
        FD( 9) = FD( 9) + DUMFD9
      END IF
      IF(.NOT. SKIPL) THEN
        FD(10) = FD(10) + DUMFD10
        FD(11) = FD(11) + DUMFD11
        FD(12) = FD(12) + DUMFD12
      END IF
C
 3000 CONTINUE
C
C     ----- ZEROTH AND FIRST DERIVATIVE INTEGRALS -----
C     ----- COMPUTE VALUES -----
C
      IF(.NOT.OUT) RETURN
C
      DO IJKLN=1,IJKLNMAX
        NX=IJKLG(2,IJKLN)
        NY=IJKLG(3,IJKLN)
        NZ=IJKLG(4,IJKLN)
        DUM    =ZERO
        DUMFXI =ZERO
        DUMFYI =ZERO
        DUMFZI =ZERO
        DUMFXJ =ZERO
        DUMFYJ =ZERO
        DUMFZJ =ZERO
        DUMFXK =ZERO
        DUMFYK =ZERO
        DUMFZK =ZERO
        DUMFXL =ZERO
        DUMFYL =ZERO
        DUMFZL =ZERO
        DO IGR=1,NG*NR
          DUMXY=XYZ(IGR,NX)*XYZ(IGR,NY)*TMPSL(IJKLN,IGR)
          DUMXZ=XYZ(IGR,NX)*XYZ(IGR,NZ)*TMPSL(IJKLN,IGR)
          DUMYZ=XYZ(IGR,NY)*XYZ(IGR,NZ)*TMPSL(IJKLN,IGR)
          DUM=DUM+XYZ(IGR,NX)*DUMYZ
          IF(.NOT. SKIPI) THEN
            DUMFXI =DUMFXI + FIXYZ(IGR,NX)*DUMYZ
            DUMFYI =DUMFYI + FIXYZ(IGR,NY)*DUMXZ
            DUMFZI =DUMFZI + FIXYZ(IGR,NZ)*DUMXY
          END IF
          IF(.NOT. SKIPJ) THEN
            DUMFXJ =DUMFXJ + FJXYZ(IGR,NX)*DUMYZ
            DUMFYJ =DUMFYJ + FJXYZ(IGR,NY)*DUMXZ
            DUMFZJ =DUMFZJ + FJXYZ(IGR,NZ)*DUMXY
          END IF
          IF(.NOT. SKIPK) THEN
            DUMFXK =DUMFXK + FKXYZ(IGR,NX)*DUMYZ
            DUMFYK =DUMFYK + FKXYZ(IGR,NY)*DUMXZ
            DUMFZK =DUMFZK + FKXYZ(IGR,NZ)*DUMXY
          END IF
          IF(.NOT. SKIPL) THEN
            DUMFXL =DUMFXL + FLXYZ(IGR,NX)*DUMYZ
            DUMFYL =DUMFYL + FLXYZ(IGR,NY)*DUMXZ
            DUMFZL =DUMFZL + FLXYZ(IGR,NZ)*DUMXY
          END IF
        END DO
        FDUM(IJKLN)=DUM
        IF(.NOT. SKIPI) THEN
          FDTMP(IJKLN,1)=DUMFXI
          FDTMP(IJKLN,2)=DUMFYI
          FDTMP(IJKLN,3)=DUMFZI
        END IF
        IF(.NOT. SKIPJ) THEN
          FDTMP(IJKLN,4)=DUMFXJ
          FDTMP(IJKLN,5)=DUMFYJ
          FDTMP(IJKLN,6)=DUMFZJ
        END IF
        IF(.NOT. SKIPK) THEN
          FDTMP(IJKLN,7)=DUMFXK
          FDTMP(IJKLN,8)=DUMFYK
          FDTMP(IJKLN,9)=DUMFZK
        END IF
        IF(.NOT. SKIPL) THEN
          FDTMP(IJKLN,10)=DUMFXL
          FDTMP(IJKLN,11)=DUMFYL
          FDTMP(IJKLN,12)=DUMFZL
        END IF
      END DO
C
C     ----- ZEROTH AND FIRST DERIVATIVE INTEGRALS -----
C     ----- STORE VALUES -----
C
      DO IJKLN=1,IJKLNMAX
        NN=IJKLG(1,IJKLN)
        GIJKL(NN)=GIJKL(NN)+FDUM(IJKLN)
        IF(.NOT. SKIPI) THEN
          FIJKL( 1,NN)=FIJKL( 1,NN)+FDTMP(IJKLN,1)
          FIJKL( 2,NN)=FIJKL( 2,NN)+FDTMP(IJKLN,2)
          FIJKL( 3,NN)=FIJKL( 3,NN)+FDTMP(IJKLN,3)
        END IF
        IF(.NOT. SKIPJ) THEN
          FIJKL( 4,NN)=FIJKL( 4,NN)+FDTMP(IJKLN,4)
          FIJKL( 5,NN)=FIJKL( 5,NN)+FDTMP(IJKLN,5)
          FIJKL( 6,NN)=FIJKL( 6,NN)+FDTMP(IJKLN,6)
        END IF
        IF(.NOT. SKIPK) THEN
          FIJKL( 7,NN)=FIJKL( 7,NN)+FDTMP(IJKLN,7)
          FIJKL( 8,NN)=FIJKL( 8,NN)+FDTMP(IJKLN,8)
          FIJKL( 9,NN)=FIJKL( 9,NN)+FDTMP(IJKLN,9)
        END IF
        IF(.NOT. SKIPL) THEN
          FIJKL(10,NN)=FIJKL(10,NN)+FDTMP(IJKLN,10)
          FIJKL(11,NN)=FIJKL(11,NN)+FDTMP(IJKLN,11)
          FIJKL(12,NN)=FIJKL(12,NN)+FDTMP(IJKLN,12)
        END IF
      END DO
C
      RETURN
      END
C*MODULE VECTOR  *DECK FORMS        (INT2A)
*VOCL TOTAL,SCALAR
      SUBROUTINE FORMS(GHONDO)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION GHONDO(*)
      DIMENSION GG(225)
C
      COMMON /DENS  / DKL(225),DIJ(225)
      COMMON /INTDEX/ IJGT(225),IJX(225),IJY(225),IJZ(225),IK(225),
     *                KLGT(225),KLX(225),KLY(225),KLZ(225)
      COMMON /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJ,IJ,KL,IJKL
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /XYZ   / XIN(625,9),YIN(625,9),ZIN(625,9)
C
C     SAMPLE VECTOR LENGTHS FOR VARIOUS IJ OR KL SHELL TYPES ARE
C     DS = 6*1=6  PP=3*3=9, DD=6*6=36, DP=6*3=18, AND SO ON.
C     THE VALUE FOR IBM WAS OPTIMIZED, OTHERS ARE EDUCATED (?) GUESSES.
C
*CRY  PARAMETER (NVECT=7)
*IBM  PARAMETER (NVECT=7)
*FUJ  PARAMETER (NVECT=7)
*NEC  PARAMETER (NVECT=7)
C
C     ----- FORM INTEGRALS OVER FUNCTIONS -----
C     DIMENSIONING XIN(81,5), AND ROLLING UP THE COMPUTATION
C     OF GHONDO IN A LOOP OF LENGTH NROOTS ADDS 33 SECONDS TO
C     A 240 SECOND INTEGRAL COMPUTATION JOB.  LEAVE IT UNROLLED.
C
C     THE COMPUTATION OF GG(K) IS VECTORIZABLE, WHILE THE
C     COMPUTATION OF GHONDO(N)=GHONDO(N)+... IS DATA DEPENDENT.
C
      GO TO (100,200,300,400,500,600,700,800,900), NROOTS
C
C     SCALAR CODE FOR NROOTS=1
C
  100 CONTINUE
CDEC$ NOVECTOR
      DO 140 I = 1,IJ
         D1 = DIJ(I)
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         N1 = IJGT(I)
         KMAX = IK(I)
*VDIR NOVECTOR
         DO 130 K = 1,KMAX
            MX = NX+KLX(K)
            MY = NY+KLY(K)
            MZ = NZ+KLZ(K)
            N = N1+KLGT(K)
            GHONDO(N) = GHONDO(N)
     *            + (XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1))*D1*DKL(K)
  130    CONTINUE
  140 CONTINUE
      RETURN
C
C     SCALAR CODE FOR NROOTS=2
C
  200 CONTINUE
CDEC$ NOVECTOR
      DO 240 I = 1,IJ
         D1 = DIJ(I)
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         N1 = IJGT(I)
*VDIR NOVECTOR
         KMAX = IK(I)
         DO 230 K = 1,KMAX
            MX = NX+KLX(K)
            MY = NY+KLY(K)
            MZ = NZ+KLZ(K)
            N = N1+KLGT(K)
            GHONDO(N) = GHONDO(N)
     *             +(XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1)
     *             + XIN(MX,2)*YIN(MY,2)*ZIN(MZ,2))*D1*DKL(K)
  230    CONTINUE
  240 CONTINUE
      RETURN
C
C     CODE FOR NROOTS=3.  BOTH VECTOR AND SCALAR CODE
C
  300 CONTINUE
      KMAX = IK(IJ)
      IF(IK(1).EQ.KMAX  .AND.  IJ.GE.KMAX
     *                  .AND.  IJ.GE.NVECT) GO TO 360
CDEC$ NOVECTOR
      DO 340 I = 1,IJ
         D1 = DIJ(I)
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         N1 = IJGT(I)
         KMAX = IK(I)
         IF(KMAX.GE.NVECT) THEN
*VOCL LOOP,VECTOR
            DO 310 K = 1,KMAX
               MX = NX+KLX(K)
               MY = NY+KLY(K)
               MZ = NZ+KLZ(K)
               GG(K) = (XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1)
     *                + XIN(MX,2)*YIN(MY,2)*ZIN(MZ,2)
     *                + XIN(MX,3)*YIN(MY,3)*ZIN(MZ,3))*D1*DKL(K)
*VDIR NOVECTOR
  310       CONTINUE
            DO 320 K = 1,KMAX
               N = N1+KLGT(K)
               GHONDO(N) = GHONDO(N) + GG(K)
  320       CONTINUE
*VDIR NOVECTOR
         ELSE
            DO 330 K = 1,KMAX
               MX = NX+KLX(K)
               MY = NY+KLY(K)
               MZ = NZ+KLZ(K)
               N = N1+KLGT(K)
               GHONDO(N) = GHONDO(N)
     *                +(XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1)
     *                + XIN(MX,2)*YIN(MY,2)*ZIN(MZ,2)
     *                + XIN(MX,3)*YIN(MY,3)*ZIN(MZ,3))*D1*DKL(K)
  330       CONTINUE
         END IF
  340 CONTINUE
      RETURN
C
  360 CONTINUE
CDEC$ NOVECTOR
      DO 390 K = 1,KMAX
         D1 = DKL(K)
         NX = KLX(K)
         NY = KLY(K)
         NZ = KLZ(K)
         N1 = KLGT(K)
*VOCL LOOP,VECTOR
         DO 370 I = 1,IJ
            MX = NX+IJX(I)
            MY = NY+IJY(I)
            MZ = NZ+IJZ(I)
            GG(I) = (XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1)
     *             + XIN(MX,2)*YIN(MY,2)*ZIN(MZ,2)
     *             + XIN(MX,3)*YIN(MY,3)*ZIN(MZ,3))*D1*DIJ(I)
  370    CONTINUE
*VDIR NOVECTOR
         DO 380 I = 1,IJ
            N = N1+IJGT(I)
            GHONDO(N) = GHONDO(N) + GG(I)
  380    CONTINUE
  390 CONTINUE
      RETURN
C
C     VECTOR CODE FOR NROOTS=4
C
  400 CONTINUE
      KMAX = IK(IJ)
      IF(IK(1).EQ.KMAX  .AND.  IJ.GT.KMAX) GO TO 460
C
CDEC$ NOVECTOR
      DO 440 I = 1,IJ
         D1 = DIJ(I)
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         N1 = IJGT(I)
         KMAX = IK(I)
*VOCL LOOP,VECTOR
         DO 410 K = 1,KMAX
            MX = NX+KLX(K)
            MY = NY+KLY(K)
            MZ = NZ+KLZ(K)
            GG(K) = (XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1)
     *             + XIN(MX,2)*YIN(MY,2)*ZIN(MZ,2)
     *             + XIN(MX,3)*YIN(MY,3)*ZIN(MZ,3)
     *             + XIN(MX,4)*YIN(MY,4)*ZIN(MZ,4))*D1*DKL(K)
  410    CONTINUE
*VDIR NOVECTOR
         DO 420 K = 1,KMAX
            N = N1 + KLGT(K)
            GHONDO(N) = GHONDO(N) + GG(K)
  420    CONTINUE
  440 CONTINUE
      RETURN
C
  460 CONTINUE
CDEC$ NOVECTOR
      DO 490 K = 1,KMAX
         D1 = DKL(K)
         NX = KLX(K)
         NY = KLY(K)
         NZ = KLZ(K)
         N1 = KLGT(K)
*VOCL LOOP,VECTOR
         DO 470 I = 1,IJ
            MX = NX+IJX(I)
            MY = NY+IJY(I)
            MZ = NZ+IJZ(I)
            GG(I) = (XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1)
     *             + XIN(MX,2)*YIN(MY,2)*ZIN(MZ,2)
     *             + XIN(MX,3)*YIN(MY,3)*ZIN(MZ,3)
     *             + XIN(MX,4)*YIN(MY,4)*ZIN(MZ,4))*D1*DIJ(I)
  470    CONTINUE
*VDIR NOVECTOR
         DO 480 I = 1,IJ
            N = N1+IJGT(I)
            GHONDO(N) = GHONDO(N) + GG(I)
  480    CONTINUE
  490 CONTINUE
      RETURN
C
C     VECTOR CODE FOR NROOTS=5
C
  500 CONTINUE
      KMAX = IK(IJ)
      IF(IK(1).EQ.KMAX  .AND.  IJ.GT.KMAX) GO TO 560
C
CDEC$ NOVECTOR
      DO 540 I = 1,IJ
         D1 = DIJ(I)
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         N1 = IJGT(I)
         KMAX = IK(I)
*VOCL LOOP,VECTOR
         DO 510 K = 1,KMAX
            MX = NX+KLX(K)
            MY = NY+KLY(K)
            MZ = NZ+KLZ(K)
            GG(K) = (XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1)
     *             + XIN(MX,2)*YIN(MY,2)*ZIN(MZ,2)
     *             + XIN(MX,3)*YIN(MY,3)*ZIN(MZ,3)
     *             + XIN(MX,4)*YIN(MY,4)*ZIN(MZ,4)
     *             + XIN(MX,5)*YIN(MY,5)*ZIN(MZ,5))*D1*DKL(K)
  510    CONTINUE
*VDIR NOVECTOR
         DO 520 K = 1,KMAX
            N = N1+KLGT(K)
            GHONDO(N) = GHONDO(N) + GG(K)
  520    CONTINUE
  540 CONTINUE
      RETURN
C
  560 CONTINUE
CDEC$ NOVECTOR
      DO 590 K = 1,KMAX
         D1 = DKL(K)
         NX = KLX(K)
         NY = KLY(K)
         NZ = KLZ(K)
         N1 = KLGT(K)
*VOCL LOOP,VECTOR
         DO 570 I = 1,IJ
            MX = NX+IJX(I)
            MY = NY+IJY(I)
            MZ = NZ+IJZ(I)
            GG(I) = (XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1)
     *             + XIN(MX,2)*YIN(MY,2)*ZIN(MZ,2)
     *             + XIN(MX,3)*YIN(MY,3)*ZIN(MZ,3)
     *             + XIN(MX,4)*YIN(MY,4)*ZIN(MZ,4)
     *             + XIN(MX,5)*YIN(MY,5)*ZIN(MZ,5))*D1*DIJ(I)
  570    CONTINUE
*VDIR NOVECTOR
         DO 580 I = 1,IJ
            N = N1+IJGT(I)
            GHONDO(N) = GHONDO(N) + GG(I)
  580    CONTINUE
  590 CONTINUE
      RETURN
C
C     VECTOR CODE FOR NROOTS=6
C
  600 CONTINUE
      KMAX = IK(IJ)
      IF(IK(1).EQ.KMAX  .AND.  IJ.GT.KMAX) GO TO 660
C
CDEC$ NOVECTOR
      DO 640 I = 1,IJ
         D1 = DIJ(I)
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         N1 = IJGT(I)
         KMAX = IK(I)
*VOCL LOOP,VECTOR
         DO 610 K = 1,KMAX
            MX = NX+KLX(K)
            MY = NY+KLY(K)
            MZ = NZ+KLZ(K)
            GG(K) = (XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1)
     *             + XIN(MX,2)*YIN(MY,2)*ZIN(MZ,2)
     *             + XIN(MX,3)*YIN(MY,3)*ZIN(MZ,3)
     *             + XIN(MX,4)*YIN(MY,4)*ZIN(MZ,4)
     *             + XIN(MX,5)*YIN(MY,5)*ZIN(MZ,5)
     *             + XIN(MX,6)*YIN(MY,6)*ZIN(MZ,6))*D1*DKL(K)
  610    CONTINUE
*VDIR NOVECTOR
         DO 620 K = 1,KMAX
            N = N1+KLGT(K)
            GHONDO(N) = GHONDO(N) + GG(K)
  620    CONTINUE
  640 CONTINUE
      RETURN
C
  660 CONTINUE
CDEC$ NOVECTOR
      DO 690 K = 1,KMAX
         D1 = DKL(K)
         NX = KLX(K)
         NY = KLY(K)
         NZ = KLZ(K)
         N1 = KLGT(K)
*VOCL LOOP,VECTOR
         DO 670 I = 1,IJ
            MX = NX+IJX(I)
            MY = NY+IJY(I)
            MZ = NZ+IJZ(I)
            GG(I) = (XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1)
     *             + XIN(MX,2)*YIN(MY,2)*ZIN(MZ,2)
     *             + XIN(MX,3)*YIN(MY,3)*ZIN(MZ,3)
     *             + XIN(MX,4)*YIN(MY,4)*ZIN(MZ,4)
     *             + XIN(MX,5)*YIN(MY,5)*ZIN(MZ,5)
     *             + XIN(MX,6)*YIN(MY,6)*ZIN(MZ,6))*D1*DIJ(I)
  670    CONTINUE
*VDIR NOVECTOR
         DO 680 I = 1,IJ
            N = N1+IJGT(I)
            GHONDO(N) = GHONDO(N) + GG(I)
  680    CONTINUE
  690 CONTINUE
      RETURN
C
C     VECTOR CODE FOR NROOTS=7
C
  700 CONTINUE
      KMAX = IK(IJ)
      IF(IK(1).EQ.KMAX  .AND.  IJ.GT.KMAX) GO TO 760
C
CDEC$ NOVECTOR
      DO 740 I = 1,IJ
         D1 = DIJ(I)
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         N1 = IJGT(I)
         KMAX = IK(I)
*VOCL LOOP,VECTOR
         DO 710 K = 1,KMAX
            MX = NX+KLX(K)
            MY = NY+KLY(K)
            MZ = NZ+KLZ(K)
            GG(K) = (XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1)
     *             + XIN(MX,2)*YIN(MY,2)*ZIN(MZ,2)
     *             + XIN(MX,3)*YIN(MY,3)*ZIN(MZ,3)
     *             + XIN(MX,4)*YIN(MY,4)*ZIN(MZ,4)
     *             + XIN(MX,5)*YIN(MY,5)*ZIN(MZ,5)
     *             + XIN(MX,6)*YIN(MY,6)*ZIN(MZ,6)
     *             + XIN(MX,7)*YIN(MY,7)*ZIN(MZ,7))*D1*DKL(K)
  710    CONTINUE
*VDIR NOVECTOR
         DO 720 K = 1,KMAX
            N = N1+KLGT(K)
            GHONDO(N) = GHONDO(N) + GG(K)
  720    CONTINUE
  740 CONTINUE
      RETURN
C
  760 CONTINUE
CDEC$ NOVECTOR
      DO 790 K = 1,KMAX
         D1 = DKL(K)
         NX = KLX(K)
         NY = KLY(K)
         NZ = KLZ(K)
         N1 = KLGT(K)
*VOCL LOOP,VECTOR
         DO 770 I = 1,IJ
            MX = NX+IJX(I)
            MY = NY+IJY(I)
            MZ = NZ+IJZ(I)
            GG(I) = (XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1)
     *             + XIN(MX,2)*YIN(MY,2)*ZIN(MZ,2)
     *             + XIN(MX,3)*YIN(MY,3)*ZIN(MZ,3)
     *             + XIN(MX,4)*YIN(MY,4)*ZIN(MZ,4)
     *             + XIN(MX,5)*YIN(MY,5)*ZIN(MZ,5)
     *             + XIN(MX,6)*YIN(MY,6)*ZIN(MZ,6)
     *             + XIN(MX,7)*YIN(MY,7)*ZIN(MZ,7))*D1*DIJ(I)
  770    CONTINUE
*VDIR NOVECTOR
         DO 780 I = 1,IJ
            N = N1+IJGT(I)
            GHONDO(N) = GHONDO(N) + GG(I)
  780    CONTINUE
  790 CONTINUE
      RETURN
C
C     VECTOR CODE FOR NROOTS=8
C
  800 CONTINUE
      KMAX = IK(IJ)
      IF(IK(1).EQ.KMAX  .AND.  IJ.GT.KMAX) GO TO 860
C
CDEC$ NOVECTOR
      DO 840 I = 1,IJ
         D1 = DIJ(I)
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         N1 = IJGT(I)
         KMAX = IK(I)
*VOCL LOOP,VECTOR
         DO 810 K = 1,KMAX
            MX = NX+KLX(K)
            MY = NY+KLY(K)
            MZ = NZ+KLZ(K)
            GG(K) = (XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1)
     *             + XIN(MX,2)*YIN(MY,2)*ZIN(MZ,2)
     *             + XIN(MX,3)*YIN(MY,3)*ZIN(MZ,3)
     *             + XIN(MX,4)*YIN(MY,4)*ZIN(MZ,4)
     *             + XIN(MX,5)*YIN(MY,5)*ZIN(MZ,5)
     *             + XIN(MX,6)*YIN(MY,6)*ZIN(MZ,6)
     *             + XIN(MX,7)*YIN(MY,7)*ZIN(MZ,7)
     *             + XIN(MX,8)*YIN(MY,8)*ZIN(MZ,8))*D1*DKL(K)
  810    CONTINUE
*VDIR NOVECTOR
         DO 820 K = 1,KMAX
            N = N1+KLGT(K)
            GHONDO(N) = GHONDO(N) + GG(K)
  820    CONTINUE
  840 CONTINUE
      RETURN
C
  860 CONTINUE
CDEC$ NOVECTOR
      DO 890 K = 1,KMAX
         D1 = DKL(K)
         NX = KLX(K)
         NY = KLY(K)
         NZ = KLZ(K)
         N1 = KLGT(K)
*VOCL LOOP,VECTOR
         DO 870 I = 1,IJ
            MX = NX+IJX(I)
            MY = NY+IJY(I)
            MZ = NZ+IJZ(I)
            GG(I) = (XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1)
     *             + XIN(MX,2)*YIN(MY,2)*ZIN(MZ,2)
     *             + XIN(MX,3)*YIN(MY,3)*ZIN(MZ,3)
     *             + XIN(MX,4)*YIN(MY,4)*ZIN(MZ,4)
     *             + XIN(MX,5)*YIN(MY,5)*ZIN(MZ,5)
     *             + XIN(MX,6)*YIN(MY,6)*ZIN(MZ,6)
     *             + XIN(MX,7)*YIN(MY,7)*ZIN(MZ,7)
     *             + XIN(MX,8)*YIN(MY,8)*ZIN(MZ,8))*D1*DIJ(I)
  870    CONTINUE
*VDIR NOVECTOR
         DO 880 I = 1,IJ
            N = N1+IJGT(I)
            GHONDO(N) = GHONDO(N) + GG(I)
  880    CONTINUE
  890 CONTINUE
      RETURN
C
C     VECTOR CODE FOR NROOTS=9
C
  900 CONTINUE
      KMAX = IK(IJ)
      IF(IK(1).EQ.KMAX  .AND.  IJ.GT.KMAX) GO TO 960
C
CDEC$ NOVECTOR
      DO 940 I = 1,IJ
         D1 = DIJ(I)
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         N1 = IJGT(I)
         KMAX = IK(I)
*VOCL LOOP,VECTOR
         DO 910 K = 1,KMAX
            MX = NX+KLX(K)
            MY = NY+KLY(K)
            MZ = NZ+KLZ(K)
            GG(K) = (XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1)
     *             + XIN(MX,2)*YIN(MY,2)*ZIN(MZ,2)
     *             + XIN(MX,3)*YIN(MY,3)*ZIN(MZ,3)
     *             + XIN(MX,4)*YIN(MY,4)*ZIN(MZ,4)
     *             + XIN(MX,5)*YIN(MY,5)*ZIN(MZ,5)
     *             + XIN(MX,6)*YIN(MY,6)*ZIN(MZ,6)
     *             + XIN(MX,7)*YIN(MY,7)*ZIN(MZ,7)
     *             + XIN(MX,8)*YIN(MY,8)*ZIN(MZ,8)
     *             + XIN(MX,9)*YIN(MY,9)*ZIN(MZ,9))*D1*DKL(K)
  910    CONTINUE
*VDIR NOVECTOR
         DO 920 K = 1,KMAX
            N = N1+KLGT(K)
            GHONDO(N) = GHONDO(N) + GG(K)
  920    CONTINUE
  940 CONTINUE
      RETURN
C
  960 CONTINUE
CDEC$ NOVECTOR
      DO 990 K = 1,KMAX
         D1 = DKL(K)
         NX = KLX(K)
         NY = KLY(K)
         NZ = KLZ(K)
         N1 = KLGT(K)
*VOCL LOOP,VECTOR
         DO 970 I = 1,IJ
            MX = NX+IJX(I)
            MY = NY+IJY(I)
            MZ = NZ+IJZ(I)
            GG(I) = (XIN(MX,1)*YIN(MY,1)*ZIN(MZ,1)
     *             + XIN(MX,2)*YIN(MY,2)*ZIN(MZ,2)
     *             + XIN(MX,3)*YIN(MY,3)*ZIN(MZ,3)
     *             + XIN(MX,4)*YIN(MY,4)*ZIN(MZ,4)
     *             + XIN(MX,5)*YIN(MY,5)*ZIN(MZ,5)
     *             + XIN(MX,6)*YIN(MY,6)*ZIN(MZ,6)
     *             + XIN(MX,7)*YIN(MY,7)*ZIN(MZ,7)
     *             + XIN(MX,8)*YIN(MY,8)*ZIN(MZ,8)
     *             + XIN(MX,9)*YIN(MY,9)*ZIN(MZ,9))*D1*DIJ(I)
  970    CONTINUE
*VDIR NOVECTOR
         DO 980 I = 1,IJ
            N = N1+IJGT(I)
            GHONDO(N) = GHONDO(N) + GG(I)
  980    CONTINUE
  990 CONTINUE
      RETURN
      END
C*MODULE VECTOR  *DECK GLDIAG       (EIGEN)
      SUBROUTINE GLDIAG(LDVECT,NVECT,NDIM,H,WRK,EIG,VECTOR,IERR,IWRK)
C
C     ----- GENERAL ROUTINE TO DIAGONALIZE A MATRIX -----
C     IF KDIAG = 0, USE A ROUTINE FROM THE VECTOR LIBRARY,
C                   IF AVAILABLE, OR EVVRSP OTHERWISE.
C              = 1, USE EVVRSP
C              = 2, USE GIVEIS
C              = 3, USE JACOBI
C     THIS ROUTINE DOES NOT REPLACE EIGEN, ONLY THE SUBROUTINE
C     'GLDIAG' CONTAINED WITHIN EIGEN.SRC
C
C           LDVECT = ROW DIMENSION OF VECTOR
C           NDIM   = DIMENSION (ORDER) OF MATRIX TO BE SOLVED
C           NVECT  = NUMBER OF VECTORS DESIRED
C           H      = MATRIX TO BE DIAGONALIZED
C           WRK    = N*8 W.P. REAL WORDS OF SCRATCH SPACE
C           EIG    = EIGENVECTORS (OUTPUT)
C           VECTOR = EIGENVECTORS (OUTPUT)
C           IERR   = ERROR FLAG (OUTPUT)
C           IWRK   = N INTEGER WORDS OF SCRATCH SPACE
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION VECTOR(LDVECT,*),WRK(NDIM,8),EIG(NDIM),H(*),IWRK(NDIM)
C
      LOGICAL USEEVV,GOPARR,DSKWRK,MASWRK
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /FMCOM / XX(1)
      COMMON /MACHSW/ KDIAG,ICORFL,IXDR
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      IERR = 0
      USEEVV=.FALSE.
C
C         -------- USE VECTORIZED LIBRARY ROUTINE -----
C
      IF(KDIAG.LE.0) THEN
C
C        THE FOLLOWING IBM CODE USES RELEASE 3 OF THE ESSL. IF YOU
C           NEED AN ESSL RELEASE 1 OR 2 VERSION SEND E-MAIL TO MWS.
C        IBM'S ESSL DSPEV IS SLIGHTLY FASTER THAN EVVRSP FOR SMALLER
C           MATRICES, BUT RATHER SLOWER FOR BIGGER ONES.
C        IBM'S ESSL DSPSV OUTPERFORMS EVVRSP FOR SINGLE ROOT EXTRACTION.
C
*IBM     IF(NVECT.EQ.0) THEN
*IBM        IOPT=20
*IBM        CALL DSPEV(IOPT,H,EIG,VECTOR,LDVECT,NDIM,WRK,2*NDIM)
*IBM     END IF
C
*IBM     IF(NVECT.EQ.NDIM) THEN
*IBM        IF(NDIM.LT.300) THEN
*IBM           IOPT=21
*IBM           CALL DSPEV(IOPT,H,EIG,VECTOR,LDVECT,NDIM,WRK,2*NDIM)
*IBM        ELSE
*IBM           USEEVV=.TRUE.
*IBM        END IF
*IBM     END IF
C
*IBM     IF(NVECT.GT.0  .AND.  NVECT.LT.NDIM) THEN
*IBM        IOPT=21
*IBM        CALL VALFM(LOADFM)
*IBM        LWRK = LOADFM + 1
*IBM        LAST = LWRK   + 9*NDIM
*IBM        NEED = LAST-LOADFM
*IBM        CALL GETFM(NEED)
*IBM        CALL DSPSV(IOPT,H,EIG,VECTOR,LDVECT,NDIM,NVECT,
*IBM *                 XX(LWRK),9*NDIM)
*IBM        CALL RETFM(NEED)
*IBM     END IF
C
*CRY     USEEVV=.TRUE.
*FUJ     USEEVV=.TRUE.
*NEC     USEEVV=.TRUE.
      END IF
C
C         ------ USE STEVE ELBERT'S EVVRSP ROUTINES -----
C
      IF(USEEVV .OR. KDIAG.EQ.1  .OR.  KDIAG.GT.3) THEN
         LENH = (NDIM*NDIM+NDIM)/2
         KORDER =0
         CALL EVVRSP(IW,NDIM,NVECT,LENH,LDVECT,H,WRK,IWRK,EIG,VECTOR
     *              ,KORDER,IERR)
      END IF
C
C         ------ USE MODIFIED EISPAK ROUTINE -----
C
      IF(KDIAG.EQ.2)
     *   CALL GIVEIS(NDIM,NVECT,LDVECT,H,WRK,IWRK,EIG,VECTOR,IERR)
C
C         ------ USE A JACOBI ROTATION ROUTINE -----
C
      IF(KDIAG.EQ.3) THEN
         IF(NVECT.EQ.NDIM) THEN
            CALL JACDG(H,VECTOR,EIG,IWRK,WRK,LDVECT,NDIM)
         ELSE
            IF (MASWRK) WRITE(IW,9000) NDIM,NVECT,LDVECT
            CALL ABRT
         END IF
      END IF
      RETURN
C
 9000 FORMAT(1X,'IN -GLDIAG-, NDIM,NVECT,LDVECT=',3I8/
     *       1X,'THE JACOBI CODE CANNOT COPE WITH NDIM.NE.NVECT!'/
     *       1X,'SO THIS RUN DOES NOT PERMIT KDIAG=3.')
      END
C*MODULE VECTOR  *DECK JDXYZV       (GRD2C)
      SUBROUTINE JDXYZV(GI,GIJ,GIJK,GIJKL,HIJK,HIJKL,FIJK,FIJKL,
     1 NG,NIMAX,NJMAX,NKMAX,NLMAX,NI,NJ,NK,NL,AAI,AAJ,AAK,AAL,
     2 FI,FJ,FK,FL,SII,SJJ,SKK,SLL,SIJ,SIK,SIL,SJK,SJL,SKL)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL SKIPI,SKIPJ,SKIPK,SKIPL
      LOGICAL IS,JS,KS,LS
      COMMON/DERSKP/IIAT,JJAT,KKAT,LLAT,SKIPI,SKIPJ,SKIPK,SKIPL
C
C     ----- VERSION 1 -----
C
      DIMENSION GIJKL(NG,NLMAX,NKMAX*NJMAX*NIMAX)
      DIMENSION  GIJK(NG,NLMAX,NKMAX,NJMAX*NIMAX)
      DIMENSION   GIJ(NG,NLMAX*NKMAX,NJMAX,NIMAX)
      DIMENSION    GI(NG,NLMAX*NKMAX*NJMAX,NIMAX)
      DIMENSION HIJKL(NG,NLMAX,NKMAX*NJMAX,NIMAX)
      DIMENSION  HIJK(NG,NLMAX,NKMAX,NJMAX,NIMAX)
      DIMENSION FIJKL(NG,NLMAX,NKMAX,NJMAX,NIMAX)
      DIMENSION  FIJK(NG,NLMAX,NKMAX,NJMAX*NIMAX)
      DIMENSION   AAI(NG)
      DIMENSION   AAJ(NG)
      DIMENSION   AAK(NG)
      DIMENSION   AAL(NG)
      DIMENSION    FI(NG,NLMAX*NKMAX*NJMAX,NIMAX)
      DIMENSION    FJ(NG,NLMAX*NKMAX,NJMAX,NIMAX)
      DIMENSION    FK(NG,NLMAX,NKMAX,NJMAX*NIMAX)
      DIMENSION    FL(NG,NLMAX,NKMAX*NJMAX*NIMAX)
      DIMENSION   SII(NG,NLMAX*NKMAX*NJMAX,NIMAX)
      DIMENSION   SJJ(NG,NLMAX*NKMAX,NJMAX,NIMAX)
      DIMENSION   SKK(NG,NLMAX,NKMAX,NJMAX*NIMAX)
      DIMENSION   SLL(NG,NLMAX,NKMAX*NJMAX*NIMAX)
      DIMENSION   SIJ(NG,NLMAX*NKMAX,NJMAX,NIMAX)
      DIMENSION   SIK(NG,NLMAX,NKMAX,NJMAX,NIMAX)
      DIMENSION   SIL(NG,NLMAX,NKMAX*NJMAX,NIMAX)
      DIMENSION   SJK(NG,NLMAX,NKMAX,NJMAX,NIMAX)
      DIMENSION   SJL(NG,NLMAX,NKMAX,NJMAX,NIMAX)
      DIMENSION   SKL(NG,NLMAX,NKMAX,NJMAX*NIMAX)
C
      IS=NI.EQ.1
      JS=NJ.EQ.1
      KS=NK.EQ.1
      LS=NL.EQ.1
C
C     ----- FIRST DERIVATIVES ONLY -----
C
      IF(SKIPI) GO TO 1030
C
C     ----- -FI- ONLY -----
C
      DO 1010 IG =1,NG
        DO 1010 LKJ=1,NLMAX*NKMAX*NJMAX
           FI(IG,LKJ,1)=  GI(IG,LKJ,2)*AAI(IG)
 1010 CONTINUE
      IF(IS) GO TO 1030
      DO 1020 IG =1,NG
        DO 1020 I  =2,NI
          DO 1020 LKJ=1,NLMAX*NKMAX*NJMAX
             FI(IG,LKJ,I)= GI(IG,LKJ,I+1)*AAI(IG)
     1              -GI(IG,LKJ,I-1)*(I-1)
 1020 CONTINUE
C
 1030 IF(SKIPJ) GO TO 1130
C
C     ----- -FJ- ONLY -----
C
      DO 1110 IG=1,NG
        DO 1110 I =1,NIMAX
          DO 1110 LK=1,NLMAX*NKMAX
             FJ(IG,LK,1,I)=  GIJ(IG,LK,2,I)*AAJ(IG)
 1110 CONTINUE
      IF(JS) GO TO 1130
      DO 1120 IG=1,NG
        DO 1120 I =1,NIMAX
          DO 1120 J =2,NJ
            DO 1120 LK=1,NLMAX*NKMAX
               FJ(IG,LK,J,I)= GIJ(IG,LK,J+1,I)*AAJ(IG)
     1               -GIJ(IG,LK,J-1,I)*(J-1)
 1120 CONTINUE
C
 1130 IF(SKIPK) GO TO 1230
C
C     ----- -FK- ONLY -----
C
      DO 1210 IG=1,NG
        DO 1210 JI=1,NJMAX*NIMAX
          DO 1210 L =1,NLMAX
             FK(IG,L,1,JI)=  GIJK(IG,L,2,JI)*AAK(IG)
 1210 CONTINUE
      IF(KS) GO TO 1230
      DO 1220 IG=1,NG
        DO 1220 JI=1,NJMAX*NIMAX
          DO 1220 K =2,NK
            DO 1220 L =1,NLMAX
               FK(IG,L,K,JI)= GIJK(IG,L,K+1,JI)*AAK(IG)
     1               -GIJK(IG,L,K-1,JI)*(K-1)
 1220 CONTINUE
C
 1230 IF(SKIPL) GO TO 1330
C
C     ----- -FL- AND -SLL- -----
C
      DO 1310 IG =1,NG
        DO 1310 KJI=1,NKMAX*NJMAX*NIMAX
           FL(IG,1,KJI)=  GIJKL(IG,2,KJI)*AAL(IG)
 1310 CONTINUE
      IF(LS) GO TO 1330
      DO 1320 IG =1,NG
        DO 1320 KJI=1,NKMAX*NJMAX*NIMAX
          DO 1320 L  =2,NL
             FL(IG,L,KJI)= GIJKL(IG,L+1,KJI)*AAL(IG)
     1              -GIJKL(IG,L-1,KJI)*(L-1)
 1320 CONTINUE
C
 1330 CONTINUE
      RETURN
      END
*OPTION* -NV
C     THE ABOVE LINE TURNS OFF VECTORIZATION OF THIS ROUTINE ON NEC-SX3
C*MODULE VECTOR  *DECK JDDSPD       (HSS2B)
*VOCL TOTAL,SCALAR
      SUBROUTINE JDDSPD(IJ,KL,DIJ,DKL,
     *                IJGT,KLGT,IJX,KLX,IJY,KLY,IJZ,KLZ,GOUT,
     *                NFDIJ,NFDKL,NIJKL)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION GG(169)
      DIMENSION DIJ(NFDIJ),DKL(NFDKL),IJGT(NFDIJ),KLGT(NFDKL),
     *          IJX(NFDIJ),KLX(NFDKL),IJY(NFDIJ),KLY(NFDKL),
     *          IJZ(NFDIJ),KLZ(NFDKL)
      DIMENSION GOUT(NIJKL)
C
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /DDXYZ / XINT(2250),YINT(2250),ZINT(2250)
C
C     ----- FORM INTEGRALS OVER FUNCTIONS -----
C     INTEGRAL TYPES ARE [IJ/KL] [IJ/K'L] [IJ/K''L], WITH OR
C     WITHOUT VARIOUS 'S ON THE I AND J (NO MORE THAN 2 PRIMES).
C     TYPE [IJ/K'L'] DOES NOT OCCUR.  SAMPLE VECTOR LENGTHS KL ARE
C        SS = 1    S'S=PS    = 3    S''S=(S+D)S  =  7
C        SP = 3    S'P=PP    = 9    S''P=(S+D)P  = 21
C        DD =36    D'D=(P+F)D=78    D''D=(S+D+G)D=132
C     WORST CASE IS [D'D'/KL] GIVING IJ = 169 FROM (P+F)**2
C
      GO TO (100,200,300,400,500,600),NROOTS
C
C     SCALAR CODE FOR NROOTS=1, VECTOR CODE FOR ALL OTHER CASES
C
  100 CONTINUE
CDEC$ NOVECTOR
      DO 160 I = 1,IJ
         D1 = DIJ(I)
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         N1 = IJGT(I)
         DO 120 K = 1,KL
            MX = NX+KLX(K)
            MY = NY+KLY(K)
            MZ = NZ+KLZ(K)
            N = N1+KLGT(K)
            GOUT(N) = (XINT(MX)*YINT(MY)*ZINT(MZ))*D1*DKL(K) + GOUT(N)
  120    CONTINUE
  160 CONTINUE
      RETURN
C
C     CODE FOR NROOTS=2
C
  200 CONTINUE
CDEC$ NOVECTOR
      DO 260 I = 1,IJ
         D1 = DIJ(I)
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         N1 = IJGT(I)
*VOCL LOOP,VECTOR
         DO 220 K = 1,KL
            MX = NX+KLX(K)
            MY = NY+KLY(K)
            MZ = NZ+KLZ(K)
            GG(K) = (XINT(MX  )*YINT(MY  )*ZINT(MZ  )
     *              +XINT(MX+1)*YINT(MY+1)*ZINT(MZ+1))*D1*DKL(K)
  220    CONTINUE
         DO 240 K = 1,KL
            N = N1+KLGT(K)
            GOUT(N) = GOUT(N) + GG(K)
  240    CONTINUE
  260 CONTINUE
      RETURN
C
C     CODE FOR NROOTS=3
C
  300 CONTINUE
      IF(IJ.GT.KL) GO TO 350
C
CDEC$ NOVECTOR
      DO 330 I = 1,IJ
         D1 = DIJ(I)
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         N1 = IJGT(I)
*VOCL LOOP,VECTOR
         DO 310 K = 1,KL
            MX = NX+KLX(K)
            MY = NY+KLY(K)
            MZ = NZ+KLZ(K)
            GG(K) = (XINT(MX  )*YINT(MY  )*ZINT(MZ  )
     *              +XINT(MX+1)*YINT(MY+1)*ZINT(MZ+1)
     *              +XINT(MX+2)*YINT(MY+2)*ZINT(MZ+2))*D1*DKL(K)
  310    CONTINUE
         DO 320 K = 1,KL
            N = N1+KLGT(K)
            GOUT(N) = GOUT(N) + GG(K)
  320    CONTINUE
  330 CONTINUE
      RETURN
C
  350 CONTINUE
CDEC$ NOVECTOR
      DO 380 K = 1,KL
         D1 = DKL(K)
         NX = KLX(K)
         NY = KLY(K)
         NZ = KLZ(K)
         N1 = KLGT(K)
*VOCL LOOP,VECTOR
         DO 360 I = 1,IJ
            MX = NX+IJX(I)
            MY = NY+IJY(I)
            MZ = NZ+IJZ(I)
            GG(I) = (XINT(MX  )*YINT(MY  )*ZINT(MZ  )
     *              +XINT(MX+1)*YINT(MY+1)*ZINT(MZ+1)
     *              +XINT(MX+2)*YINT(MY+2)*ZINT(MZ+2))*D1*DIJ(I)
  360    CONTINUE
         DO 370 I = 1,IJ
            N = N1+IJGT(I)
            GOUT(N) = GOUT(N) + GG(I)
  370    CONTINUE
  380 CONTINUE
      RETURN
C
C     CODE FOR NROOTS=4
C
  400 CONTINUE
      IF(IJ.GT.KL) GO TO 450
C
CDEC$ NOVECTOR
      DO 430 I = 1,IJ
         D1 = DIJ(I)
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         N1 = IJGT(I)
*VOCL LOOP,VECTOR
         DO 410 K = 1,KL
            MX = NX+KLX(K)
            MY = NY+KLY(K)
            MZ = NZ+KLZ(K)
            GG(K) = (XINT(MX  )*YINT(MY  )*ZINT(MZ  )
     1              +XINT(MX+1)*YINT(MY+1)*ZINT(MZ+1)
     2              +XINT(MX+2)*YINT(MY+2)*ZINT(MZ+2)
     3              +XINT(MX+3)*YINT(MY+3)*ZINT(MZ+3))*D1*DKL(K)
  410    CONTINUE
         DO 420 K = 1,KL
            N = N1+KLGT(K)
            GOUT(N) = GOUT(N) + GG(K)
  420    CONTINUE
  430 CONTINUE
      RETURN
C
  450 CONTINUE
CDEC$ NOVECTOR
      DO 480 K = 1,KL
         D1 = DKL(K)
         NX = KLX(K)
         NY = KLY(K)
         NZ = KLZ(K)
         N1 = KLGT(K)
*VOCL LOOP,VECTOR
         DO 460 I = 1,IJ
            MX = NX+IJX(I)
            MY = NY+IJY(I)
            MZ = NZ+IJZ(I)
            GG(I) = (XINT(MX  )*YINT(MY  )*ZINT(MZ  )
     *              +XINT(MX+1)*YINT(MY+1)*ZINT(MZ+1)
     *              +XINT(MX+2)*YINT(MY+2)*ZINT(MZ+2)
     *              +XINT(MX+3)*YINT(MY+3)*ZINT(MZ+3))*D1*DIJ(I)
  460    CONTINUE
         DO 470 I = 1,IJ
            N = N1+IJGT(I)
            GOUT(N) = GOUT(N) + GG(I)
  470    CONTINUE
  480 CONTINUE
      RETURN
C
C     CODE FOR NROOTS=5
C
  500 CONTINUE
      IF(IJ.GT.KL) GO TO 550
C
CDEC$ NOVECTOR
      DO 530 I = 1,IJ
         D1 = DIJ(I)
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         N1 = IJGT(I)
*VOCL LOOP,VECTOR
         DO 510 K = 1,KL
            MX = NX+KLX(K)
            MY = NY+KLY(K)
            MZ = NZ+KLZ(K)
            GG(K) = (XINT(MX  )*YINT(MY  )*ZINT(MZ  )
     *              +XINT(MX+1)*YINT(MY+1)*ZINT(MZ+1)
     *              +XINT(MX+2)*YINT(MY+2)*ZINT(MZ+2)
     *              +XINT(MX+3)*YINT(MY+3)*ZINT(MZ+3)
     *              +XINT(MX+4)*YINT(MY+4)*ZINT(MZ+4))*D1*DKL(K)
  510    CONTINUE
         DO 520 K = 1,KL
            N = N1+KLGT(K)
            GOUT(N) = GOUT(N) + GG(K)
  520    CONTINUE
  530 CONTINUE
      RETURN
C
  550 CONTINUE
CDEC$ NOVECTOR
      DO 580 K = 1,KL
         D1 = DKL(K)
         NX = KLX(K)
         NY = KLY(K)
         NZ = KLZ(K)
         N1 = KLGT(K)
*VOCL LOOP,VECTOR
         DO 560 I = 1,IJ
            MX = NX+IJX(I)
            MY = NY+IJY(I)
            MZ = NZ+IJZ(I)
            GG(I) = (XINT(MX  )*YINT(MY  )*ZINT(MZ  )
     *              +XINT(MX+1)*YINT(MY+1)*ZINT(MZ+1)
     *              +XINT(MX+2)*YINT(MY+2)*ZINT(MZ+2)
     *              +XINT(MX+3)*YINT(MY+3)*ZINT(MZ+3)
     *              +XINT(MX+4)*YINT(MY+4)*ZINT(MZ+4))*D1*DIJ(I)
  560    CONTINUE
         DO 570 I = 1,IJ
            N = N1+IJGT(I)
            GOUT(N) = GOUT(N) + GG(I)
  570    CONTINUE
  580 CONTINUE
      RETURN
C
C     CODE FOR NROOTS=6
C
  600 CONTINUE
      IF(IJ.GT.KL) GO TO 650
C
CDEC$ NOVECTOR
      DO 630 I=1,IJ
         D1 = DIJ(I)
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         N1 = IJGT(I)
*VOCL LOOP,VECTOR
         DO 610 K = 1,KL
            MX = NX+KLX(K)
            MY = NY+KLY(K)
            MZ = NZ+KLZ(K)
            GG(K) = (XINT(MX  )*YINT(MY  )*ZINT(MZ  )
     *              +XINT(MX+1)*YINT(MY+1)*ZINT(MZ+1)
     *              +XINT(MX+2)*YINT(MY+2)*ZINT(MZ+2)
     *              +XINT(MX+3)*YINT(MY+3)*ZINT(MZ+3)
     *              +XINT(MX+4)*YINT(MY+4)*ZINT(MZ+4)
     *              +XINT(MX+5)*YINT(MY+5)*ZINT(MZ+5))*D1*DKL(K)
  610    CONTINUE
         DO 620 K = 1,KL
            N = N1+KLGT(K)
            GOUT(N) = GOUT(N) + GG(K)
  620    CONTINUE
  630 CONTINUE
      RETURN
C
  650 CONTINUE
CDEC$ NOVECTOR
      DO 680 K = 1,KL
         D1 = DKL(K)
         NX = KLX(K)
         NY = KLY(K)
         NZ = KLZ(K)
         N1 = KLGT(K)
*VOCL LOOP,VECTOR
         DO 660 I = 1,IJ
            MX = NX+IJX(I)
            MY = NY+IJY(I)
            MZ = NZ+IJZ(I)
            GG(I) = (XINT(MX  )*YINT(MY  )*ZINT(MZ  )
     *              +XINT(MX+1)*YINT(MY+1)*ZINT(MZ+1)
     *              +XINT(MX+2)*YINT(MY+2)*ZINT(MZ+2)
     *              +XINT(MX+3)*YINT(MY+3)*ZINT(MZ+3)
     *              +XINT(MX+4)*YINT(MY+4)*ZINT(MZ+4)
     *              +XINT(MX+5)*YINT(MY+5)*ZINT(MZ+5))*D1*DIJ(I)
  660    CONTINUE
         DO 670 I = 1,IJ
            N = N1+IJGT(I)
            GOUT(N) = GOUT(N) + GG(I)
  670    CONTINUE
  680 CONTINUE
      RETURN
      END
C*MODULE VECTOR  *DECK JKGNMV       (GRD2C)
      SUBROUTINE JKGNMV(GNM,NG,NMAX,MMAX,B00,B01,B10,C00,D00,F00)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL NMAXS,NMAXP,MMAXS,MMAXP
      COMMON /SHLGNM/NMAXS,NMAXP,MMAXS,MMAXP
      DIMENSION GNM(NG,NMAX,MMAX)
      DIMENSION C00(NG),D00(NG),F00(NG)
      DIMENSION B00(NG,1),B01(NG,1),B10(NG,1)
C
C     ----- G(0,0) -----
C
      DO 10 IG=1,NG
   10   GNM(IG,1,1)=F00(IG)
      IF(NMAXS.AND.MMAXS) RETURN
      IF(NMAXS) GO TO 30
C
C     ----- G(1,0) = C00 * G(0,0) -----
C
      DO 20 IG=1,NG
   20   GNM(IG,2,1)=C00(IG)*GNM(IG,1,1)
C
   30 IF(MMAXS) GO TO 60
C
C     ----- G(0,1) = D00 * G(0,0) -----
C
      DO 40 IG=1,NG
   40   GNM(IG,1,2)=D00(IG)*GNM(IG,1,1)
      IF(NMAXS) GO TO 60
C
C     ----- G(1,1) = B00 * G(0,0) + D00 * G(1,0) -----
C
      DO 50 IG=1,NG
   50   GNM(IG,2,2)=B00(IG,1)*GNM(IG,1,1)+D00(IG)*GNM(IG,2,1)
C
   60 MAX=MAX0(NMAX-1,MMAX-1)
      DO 70 IG=1,NG
        DO 70 M=2,MAX
   70     B00(IG,M)=B00(IG,M-1)+B00(IG,1)
C
      IF(NMAXP) GO TO 120
C
C     ----- G(N+1,0) = N * B10 * G(N-1,0) + C00 * G(N,0) -----
C
      DO 80 IG=1,NG
        DO 80 N=2,NMAX-1
   80     B10(IG,N)=B10(IG,N-1)+B10(IG,1)
      DO 90 IG=1,NG
        DO 90 N=2,NMAX-1
   90     GNM(IG,N+1,1)=B10(IG,N-1)*GNM(IG,N-1,1)+C00(IG)*GNM(IG,N,1)
      IF(MMAXS) GO TO 120
C
C     ----- G(N,1) = N * B00 * G(N-1,0) + D00 * G(N,0) -----
C
      DO 100 IG=1,NG
        DO 100 N=2,NMAX-1
  100     GNM(IG,N+1,2)=B00(IG,N)*GNM(IG,N,1)+D00(IG)*GNM(IG,N+1,1)
C
  120 IF(MMAXP) GO TO 170
C
C     ----- G(0,M+1) = M * B01 * G(0,M-1) + D00 * G(O,M) -----
C
      DO 130 IG=1,NG
        DO 130 M=2,MMAX-1
  130     B01(IG,M)=B01(IG,M-1)+B01(IG,1)
      DO 140 IG=1,NG
        DO 140 M=2,MMAX-1
  140     GNM(IG,1,M+1)=B01(IG,M-1)*GNM(IG,1,M-1)+D00(IG)*GNM(IG,1,M)
      IF(NMAXS) GO TO 170
C
C     ----- G(1,M) = M * B00 * G(0,M-1) + C00 * G(0,M) -----
C
      DO 150 IG=1,NG
        DO 150 M=2,MMAX-1
  150     GNM(IG,2,M+1)=B00(IG,M)*GNM(IG,1,M)+C00(IG)*GNM(IG,1,M+1)
C
  170 IF(NMAXP.OR.MMAXP) RETURN
C
C     ----- G(N+1,M) = N * B10 * G(N-1,M  ) -----
C                    +     C00 * G(N  ,M  )
C                    + M * B00 * G(N  ,M-1)
C
      DO 180 IG=1,NG
        DO 180 N=2,NMAX-1
          DO 180 M=2,MMAX-1
            GNM(IG,N+1,M+1)=B10(IG,N-1)*GNM(IG,N-1,M+1)+
     1                C00(IG    )*GNM(IG,N  ,M+1)+
     2                B00(IG,M  )*GNM(IG,N  ,M  )
  180 CONTINUE
C
      RETURN
      END
C*MODULE VECTOR  *DECK JKXYZV       (GRD2C)
      SUBROUTINE JKXYZV(GIJKL,HIJKL,GNKL,HNKL,FNKL,GNM,HNM,
     1 NG,NMAX,MMAX,NIMAX,NJMAX,NKMAX,NLMAX,DIJ,DKL,EXPNDI,EXPNDK)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL EXPNDI,EXPNDK
      DIMENSION GIJKL(NG*NLMAX*NKMAX,NJMAX,NIMAX)
      DIMENSION HIJKL(NG*NLMAX*NKMAX*NJMAX,NIMAX)
      DIMENSION  GNKL(NG,NLMAX,NKMAX,NMAX)
      DIMENSION  HNKL(NG*NLMAX*NKMAX,NMAX)
      DIMENSION  FNKL(NG*NLMAX*NKMAX*NMAX)
      DIMENSION   GNM(NG,NMAX,MMAX)
      DIMENSION   DIJ(NG)
      DIMENSION   DKL(NG)
C
C     ----- G(N,K,L) -----
C
      IF (.NOT.EXPNDK) THEN
C
CDIR$ NEXTSCALAR
        DO 30 NK=1,NKMAX
          DO 10 IG=1,NG
            DO 10 NL=1,NLMAX
              DO 10  N=1,NMAX
   10           GNKL(IG,NL,NK,N)=GNM(IG,N,NL)
          IF(NK.EQ.NKMAX) GO TO 30
          MAX=MMAX-NK
          DO 20 IG=1,NG
            DO 20  M=1,MAX
              DO 20  N=1,NMAX
   20           GNM(IG,N,M)=DKL(IG)*GNM(IG,N,M)+GNM(IG,N,M+1)
   30   CONTINUE
C
      ELSE
C
CDIR$ NEXTSCALAR
        DO 70 NL=1,NLMAX
          DO 50 IG=1,NG
            DO 50 NK=1,NKMAX
              DO 50  N=1,NMAX
   50           GNKL(IG,NL,NK,N)=GNM(IG,N,NK)
          IF(NL.EQ.NLMAX) GO TO 70
          MAX=MMAX-NL
          DO 60 IG=1,NG
            DO 60  N=1,NMAX
              DO 60  M=1,MAX
   60           GNM(IG,N,M)=DKL(IG)*GNM(IG,N,M)+GNM(IG,N,M+1)
   70   CONTINUE
C
      ENDIF
C
C     ----- G(I,J,K,L) -----
C
      IF (.NOT.EXPNDI) THEN
C
CDIR$ NEXTSCALAR
        DO 130 NI=1,NIMAX
          DO 110 IGLKJ=1,NG*NLMAX*NKMAX*NJMAX
  110       HIJKL(IGLKJ,NI)=FNKL(IGLKJ)
          IF(NI.EQ.NIMAX) GO TO 130
          MAX=NMAX-NI
          DO 120 IG=1,NG
            DO 120 NL=1,NLMAX
              DO 120 NK=1,NKMAX
                DO 120  N=1,MAX
  120             GNKL(IG,NL,NK,N)=DIJ(IG)*GNKL(IG,NL,NK,N)+
     1                             GNKL(IG,NL,NK,N+1)
  130   CONTINUE
C
      ELSE
C
CDIR$ NEXTSCALAR
        DO 170 NJ=1,NJMAX
          DO 150 IGLK=1,NG*NLMAX*NKMAX
            DO 150 NI=1,NIMAX
  150         GIJKL(IGLK,NJ,NI)=HNKL(IGLK,NI)
          IF(NJ.EQ.NJMAX) GO TO 170
          MAX=NMAX-NJ
          DO 160 IG=1,NG
            DO 160 NL=1,NLMAX
              DO 160 NK=1,NKMAX
                DO 160  N=1,MAX
  160             GNKL(IG,NL,NK,N)=DIJ(IG)*GNKL(IG,NL,NK,N)+
     1                             GNKL(IG,NL,NK,N+1)
  170   CONTINUE
      ENDIF
C
      RETURN
      END
C*MODULE VECTOR  *DECK MTARBR       (MTHLIB)
      SUBROUTINE MTARBR(A,NA,B,MB,C,LDC,INCA)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(*),B(NA,MB),C(LDC,MB)
C
*CRY  PARAMETER (ZERO=0.0E+00, ONE=1.0E+00)
*FUJ  PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
*IBM  PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
*NEC  PARAMETER (ZERO=0.0E+00, ONE=1.0E+00)
C
      COMMON /FMCOM / XX(1)
C
C     MATRIX MULTIPLY C = A*B.  A IS SYMMETRIC, B AND C ARE RECTANGULAR
C
C PARAMETERS
C     A     - THE INPUT REAL SYMMETRIC MATRIX OF ORDER NA
C             STORED IN SYMMETRIC STORAGE MODE.
C     B     - THE INPUT REAL NA BY MB RECTANGULAR MATRIX
C     NA    - THE ORDER OF MATRIX A, AND THE LEADING DIMENSION OF B
C     MB    - THE COLUMN DIMENSION OF MATRICES B AND C
C     C     - THE OUTPUT PRODUCT NA BY MB MATRIX
C     LDC   - THE LEADING DIMENSION DIMENSION OF MATRIX C
C     INCA  - ADDRESS DIFFERENCE OF TWO ADJACENT ELEMENTS OF A
C
C     EXPAND A TO SQUARE STORAGE, THEN MULTIPLY.
C
      CALL VALFM(LOADFM)
      LOC  = LOADFM + 1
      LAST = LOC    + NA*NA
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      CALL CPYTSQ(A,XX(LOC),NA,INCA)
C
*CRY  CALL SSYMM('L','U',NA,MB,   ONE,XX(LOC),NA,B,NA,ZERO,C,LDC)
C
*FUJ  CALL DSYMM('L','U',NA,MB,   ONE,XX(LOC),NA,B,NA,ZERO,C,LDC)
C
*IBM  CALL DSYMM('L','U',NA,MB,   ONE,XX(LOC),NA,B,NA,ZERO,C,LDC)
C
*NEC  CALL DSYMM('L','U',NA,MB,   ONE,XX(LOC),NA,B,NA,ZERO,C,LDC)
C
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE VECTOR  *DECK TFSQB        (MTHLIB)
      SUBROUTINE TFSQB(V,Q,T,M,N,NDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION V(NDIM,M),Q(NDIM,M),T(N)
C
      COMMON /FMCOM / XX(1)
C
*CRY  PARAMETER (ZERO=0.0E+00, ONE=1.0E+00)
*FUJ  PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
*IBM  PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
*NEC  PARAMETER (ZERO=0.0E+00, ONE=1.0E+00)
C
C     ----- BACK TRANSFORM THE SQUARE MATRIX V WITH Q -----
C                          V = Q*V
C     ON ENTRY, Q IS N X M, AND V IS M X M (M.LE.N)
C     ON EXIT, V IS N X M.   Q IS UNALTERED.
C
C     THE VECTORIZED VERSION GROWS A TEMPORARY MATRIX C,
C     COMPUTES C=Q*V, AND THEN MOVES C INTO V.
C
      CALL VALFM(LOADFM)
      LOC  = LOADFM + 1
      LAST = LOC    + N*M
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
*CRY  CALL SGEMM('N','N',N,M,M,ONE,Q,NDIM,V,NDIM,ZERO,XX(LOC),N)
C
*FUJ  CALL DGEMM('N','N',N,M,M,ONE,Q,NDIM,V,NDIM,ZERO,XX(LOC),N)
C
*IBM  CALL DGEMM('N','N',N,M,M,ONE,Q,NDIM,V,NDIM,ZERO,XX(LOC),N)
C
*NEC  CALL SGEMM('N','N',N,M,M,ONE,Q,NDIM,V,NDIM,ZERO,XX(LOC),N)
C
      CALL CPYSS(XX(LOC),N,V,NDIM,N,M)
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE VECTOR  *DECK TFTRI        (MTHLIB)
      SUBROUTINE TFTRI(H,F,Q,T,M,N,NDIM)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION H(*),F(*),Q(NDIM,M),T(N)
C
*CRY  PARAMETER (ZERO=0.0E+00, ONE=1.0E+00)
*FUJ  PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
*IBM  PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
*NEC  PARAMETER (ZERO=0.0E+00, ONE=1.0E+00)
C
      COMMON /FMCOM / XX(1)
C
C     ----- TRANSFORM THE TRIANGULAR MATRIX F USING VECTORS Q -----
C                      H = Q-DAGGER * F * Q
C     THE ORDER OF THE TRIANGULAR MATRICES H AND F ARE M AND N.
C
C     THE VECTORIZED VERSION GROWS 3 TEMPORARY MATRICES, CONVERTS
C     F TO SQUARE STORAGE, COMPUTES C=F*Q, THEN H=Q-DAGGER*C,
C     THEN CONVERTS H FROM SQUARE TO TRIANGULAR STORAGE MODE.
C
      CALL VALFM(LOADFM)
      LOCF = LOADFM + 1
      LOCH = LOCF   + N*N
      LOCC = LOCH   + M*M
      LAST = LOCC   + N*M
      NEED = LAST   - LOCF
      CALL GETFM(NEED)
C
      CALL CPYTSQ(F,XX(LOCF),N,1)
C
*CRY  CALL SSYMM('L','U',N,M,  ONE,XX(LOCF),N,Q,NDIM,ZERO,XX(LOCC),N)
*CRY  CALL SGEMM('T','N',M,M,N,ONE,Q,NDIM,XX(LOCC),N,ZERO,XX(LOCH),M)
C
*FUJ  CALL DSYMM('L','U',N,M,  ONE,XX(LOCF),N,Q,NDIM,ZERO,XX(LOCC),N)
*FUJ  CALL DGEMM('T','N',M,M,N,ONE,Q,NDIM,XX(LOCC),N,ZERO,XX(LOCH),M)
C
*IBM  CALL DSYMM('L','U',N,M,  ONE,XX(LOCF),N,Q,NDIM,ZERO,XX(LOCC),N)
*IBM  CALL DGEMM('T','N',M,M,N,ONE,Q,NDIM,XX(LOCC),N,ZERO,XX(LOCH),M)
C
*NEC  CALL SSYMM('L','U',N,M,  ONE,XX(LOCF),N,Q,NDIM,ZERO,XX(LOCC),N)
*NEC  CALL SGEMM('T','N',M,M,N,ONE,Q,NDIM,XX(LOCC),N,ZERO,XX(LOCH),M)
C
      CALL CPYSQT(XX(LOCH),H,M,1)
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE VECTOR  *DECK CPYSS
      SUBROUTINE CPYSS(A,LDA,B,LDB,N,M)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,M),B(LDB,M)
      PARAMETER (SMALL=1.0D-10, ZERO=0.0D+00)
C
C     COPY SQUARE MATRIX A (N X M) TO B, CHECKING FOR SMALL NOS.
C
      DO 200 J=1,M
         DO 100 I=1,N
            VAL = A(I,J)
            IF(ABS(VAL).LT.SMALL) VAL=ZERO
            B(I,J) = VAL
  100    CONTINUE
  200 CONTINUE
      RETURN
      END
C*MODULE VECTOR  *DECK DIRFCK      (SCFLIB)
      SUBROUTINE DIRFCK(SCFTYP,IA,DA,FA,DB,FB,GHONDO,L2,NINT,NXYZ)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL IANDJ,KANDL,SAME,OUT,UROHF
C
      DIMENSION IA(*),DA(*),FA(*),DB(*),FB(*),GHONDO(*)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXAO=8192)
      PARAMETER (ZERO=0.0D+00)
C
      COMMON /DFTPAR/ DFTTYP(20),EXENA,EXENB,EXENC,IDFT34,NAUXFUN,
     *                                                    NAUXSHL
      COMMON /ERIOUT/ ISH,JSH,KSH,LSH,LSTRI,LSTRJ,LSTRK,LSTRL
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MISC  / IANDJ,KANDL,SAME
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJX,IJ,KL,IJKL
      COMMON /SHLT  / TOL,CUTOFF,ICOUNT,OUT
C
      DATA HALF /0.5D+00/
      DATA RHF,UHF,ROHF/8HRHF     ,8HUHF     ,8HROHF     /
      DATA GVB/8HGVB     /
C
      DIMENSION VALV((KMAX(ISH)+1)*(KMAX(JSH)+1)*
     *               (KMAX(KSH)+1)*(KMAX(LSH)+1)*24)
      DIMENSION LOCV((KMAX(ISH)+1)*(KMAX(JSH)+1)*
     *               (KMAX(KSH)+1)*(KMAX(LSH)+1)*24,4)
      DIMENSION NIJKLV((KMAX(ISH)+1)*(KMAX(JSH)+1)*
     *                 (KMAX(KSH)+1)*(KMAX(LSH)+1)*24,6)
      DIMENSION VIJKLV((KMAX(ISH)+1)*(KMAX(JSH)+1)*
     *                 (KMAX(KSH)+1)*(KMAX(LSH)+1)*24,6)
C
      SAME  = ISH .EQ. KSH .AND. JSH .EQ. LSH
      IANDJ = ISH .EQ. JSH
      KANDL = KSH .EQ. LSH
C
      MINI = KMIN(ISH)
      MINJ = KMIN(JSH)
      MINK = KMIN(KSH)
      MINL = KMIN(LSH)
      MAXI = KMAX(ISH)
      MAXJ = KMAX(JSH)
      MAXK = KMAX(KSH)
      MAXL = KMAX(LSH)
      LOCI = KLOC(ISH)-MINI
      LOCJ = KLOC(JSH)-MINJ
      LOCK = KLOC(KSH)-MINK
      LOCL = KLOC(LSH)-MINL
C
C     ----- FORM FOCK OPERATOR DIRECTLY FROM INTEGRALS -----
C     THIS ROUTINE WAS PIECED TOGETHER FROM QOUT AND HSTAR
C     BY FRANK JENSEN AT ODENSE UNIVERSITY IN MARCH 1990.
C     SCF FUNCTIONS BESIDES RHF ADDED BY MWS IN AUGUST 1991.
C
C     NOTE THAT OFF-DIAGONAL ELEMENTS WILL NEED TO BE HALVED LATER.
C
      HFSCAL=DFTTYP(3)
      UROHF = SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF
      CUTINT = CUTOFF
C
      IF(SCFTYP.EQ.GVB) THEN
         NSHL = NHAM
         IF(NCO.GT.0) NSHL=NSHL-1
      END IF
C
C     UNPACK AHEAD OF TIME TO ALLOW VECTORIZATION OF BULK OF WORK.
C     THIS WAS CHANGED BY MIKE SCHMIDT TO EXACTLY RESEMBLE THE
C     CODE IN THE REAL DIRFCK ROUTINE (SEE SCFLIB) WHEN THE POPLE
C     INDEXING WAS REMOVED.
C
      IJN = 0
      JMAX = MAXJ
      NINTORIG = NINT + 1
      IF (SAME) THEN
        DO I = MINI,MAXI
          I_ERIC = (I-MINI)*LSTRI + 1
          IF (IANDJ) JMAX = I
          DO 340 J = MINJ,JMAX
            IJ_ERIC = (J-MINJ)*LSTRJ + I_ERIC
            LMAX = MAXL
            IJN = IJN+1
            KLN=0
            DO K =  MINK,MAXK
               IJK_ERIC = (K-MINK)*LSTRK + IJ_ERIC
               IF (KANDL) LMAX = K
               LMAX2=LMAX
               IF( (KLN+1).GT.IJN)THEN
                 GO TO 340
               ELSE
                 LMAX2 = MIN( MINL+(IJN-(KLN+1)), LMAX )
               ENDIF
               DO L = MINL,LMAX2
                  KLN = KLN+1
                  IJKL_ERIC = (L-MINL)*LSTRL + IJK_ERIC
                  VAL = GHONDO( IJKL_ERIC )
                  IF(ABS(VAL).GE.CUTINT) THEN
                     NINT=NINT+1
                     NINTCURR=NINT-NINTORIG+1
                     LOCV(NINTCURR,1) = I
                     LOCV(NINTCURR,2) = J
                     LOCV(NINTCURR,3) = K
                     LOCV(NINTCURR,4) = L
                     VALV(NINTCURR) = VAL
                  END IF
               END DO
            END DO
  340     CONTINUE
        END DO
C
C       FOR SAME = .FALSE., THE IF(KLN.GT.IJN)  WILL NEVER OCCUR
C       SO WE CAN SAFELY STRIP IJN AND KLN OUT OF THE LOOP NEST
C
      ELSE   ! .NOT. SAME
        DO I = MINI,MAXI
          I_ERIC = (I-MINI)*LSTRI + 1
          IF (IANDJ) JMAX = I
          DO J = MINJ,JMAX
            IJ_ERIC = (J-MINJ)*LSTRJ + I_ERIC
            LMAX = MAXL
            DO K =  MINK,MAXK
               IJK_ERIC = (K-MINK)*LSTRK + IJ_ERIC
               IF (KANDL) LMAX = K
               DO L = MINL,LMAX
                  IJKL_ERIC = (L-MINL)*LSTRL + IJK_ERIC
                  VAL = GHONDO( IJKL_ERIC )
                  IF(ABS(VAL).GE.CUTINT) THEN
                     NINT=NINT+1
                     NINTCURR=NINT-NINTORIG+1
                     LOCV(NINTCURR,1) = I
                     LOCV(NINTCURR,2) = J
                     LOCV(NINTCURR,3) = K
                     LOCV(NINTCURR,4) = L
                     VALV(NINTCURR) = VAL
                  END IF
               END DO
            END DO
          END DO
        END DO
      END IF
C
C     NOW WE KNOW THE ACTUAL NUMBER OF LOOPS (NINTMAX).
C
      NINTMAX=NINT-NINTORIG+1
C
C     NOW USE VALUES TO FIND NIJ,NJK, ETC.
C     THIS LOOP NOW VECTORIZES QUITE NICELY.
C
      DO M=1,NINTMAX
        II = LOCI + LOCV(M,1)
        JJ = LOCJ + LOCV(M,2)
        KK = LOCK + LOCV(M,3)
        LL = LOCL + LOCV(M,4)
        VAL = VALV(M)
C
        IF (II .LT. JJ) THEN
          N = II
          II = JJ
          JJ = N
        END IF
        IF (KK .LT. LL) THEN
          N = KK
          KK = LL
          LL = N
        END IF
        IIMINKK = II-KK
        IF ( IIMINKK .LT. 0 .OR.
     *          (IIMINKK .EQ. 0 .AND. JJ .LT. LL) ) THEN
          N = II
          II = KK
          KK = N
          N = JJ
          JJ = LL
          LL = N
        END IF
C
        IF(II.EQ.JJ) VAL = VAL*HALF
        IF(KK.EQ.LL) VAL = VAL*HALF
        IF(II.EQ.KK  .AND.  JJ.EQ.LL) VAL = VAL*HALF
C
        NIJ = IA(II)+JJ
        NKL = IA(KK)+LL
        NIK = IA(II)+KK
        NIL = IA(II)+LL
        IF(JJ.LT.KK) THEN
          NJK = IA(KK)+JJ
        ELSE
          NJK = IA(JJ)+KK
        END IF
        IF(JJ.LT.LL) THEN
          NJL = IA(LL)+JJ
        ELSE
          NJL = IA(JJ)+LL
        END IF
        NIJKLV(M,1) = NIJ
        NIJKLV(M,2) = NKL
        NIJKLV(M,3) = NIK
        NIJKLV(M,4) = NIL
        NIJKLV(M,5) = NJK
        NIJKLV(M,6) = NJL
        VALV(M) = VAL
        IF(SCFTYP.EQ.RHF) THEN
          IF(NXYZ.EQ.1) THEN
              VALHFS = VAL*HFSCAL
              VAL4 = VAL*4
              VIJKLV(M,1) = VAL4*DA(NKL)
              VIJKLV(M,2) = VAL4*DA(NIJ)
              VIJKLV(M,3) = ZERO - VALHFS*DA(NJL)
              VIJKLV(M,4) = ZERO - VALHFS*DA(NJK)
              VIJKLV(M,5) = ZERO - VALHFS*DA(NIL)
              VIJKLV(M,6) = ZERO - VALHFS*DA(NIK)
          END IF
        END IF
      END DO
C
C     WE NOW HAVE EVERYTHING READY, PUT INTEGRALS INTO FOCK MATRIX.
C
C       HFSCAL DISTINGUISHES CODE FOR HF OR DFT
C                  IF(HFSCAL .EQ. 1.0D+00) THEN
C           WE CAN REMOVE THIS TEST BECAUSE MULTIPLYING BY 1.0
C             DOES NOT CHANGE THE ANSWER FOR HF CASE.
C
C                    RHF
C
      IF(SCFTYP.EQ.RHF) THEN
C     NXYZ DISTINGUISHES CODE FOR RHF OR RHF RESPONSE EQUATIONS
        IF(NXYZ.EQ.1) THEN
C         REPEATED INDICES ON BOTH DIMENSIONS PREVENTS THIS FROM
C         VECTORIZING CLEANLY, THUS WE SPLIT IT OFF FROM THE MAIN
C         COMPUTATION LOOP ABOVE.
C         THE M LOOP RUNS TERRIBLE IN VECTOR.  *ALWAYS* RUN IT SCALAR.
          DO N=1,6
!DIR$ NEXTSCALAR
            DO M=1,NINTMAX
              NIJ = NIJKLV(M,N)
              FA(NIJ) = FA(NIJ) + VIJKLV(M,N)
            END DO
          END DO
        ELSE   ! NXYZ .NE. 1
          NXYZM1 = NXYZ-1
          DO M=1,NINTMAX
            NIJ = (NIJKLV(M,1)-1)*NXYZ+1
            NKL = (NIJKLV(M,2)-1)*NXYZ+1
            NIK = (NIJKLV(M,3)-1)*NXYZ+1
            NIL = (NIJKLV(M,4)-1)*NXYZ+1
            NJK = (NIJKLV(M,5)-1)*NXYZ+1
            NJL = (NIJKLV(M,6)-1)*NXYZ+1
            VALHFS = VALV(M)*HFSCAL
            VAL4 = VALV(M)*4
            NIJN = NIJ+NXYZM1
            NKLN = NKL+NXYZM1
            NIKN = NIK+NXYZM1
            NILN = NIL+NXYZM1
            NJKN = NJK+NXYZM1
            NJLN = NJL+NXYZM1
            FA(NIJ:NIJN) = FA(NIJ:NIJN)+VAL4*DA(NKL:NKLN)
            FA(NKL:NKLN) = FA(NKL:NKLN)+VAL4*DA(NIJ:NIJN)
            FA(NIK:NIKN) = FA(NIK:NIKN)-VALHFS*DA(NJL:NJLN)
            FA(NIL:NILN) = FA(NIL:NILN)-VALHFS*DA(NJK:NJKN)
            FA(NJK:NJKN) = FA(NJK:NJKN)-VALHFS*DA(NIL:NILN)
            FA(NJL:NJLN) = FA(NJL:NJLN)-VALHFS*DA(NIK:NIKN)
          END DO
        END IF
C
C                    UROHF
C
C       HFSCAL DISTINGUISHES CODE FOR HF OR DFT
C                  IF(HFSCAL .EQ. 1.0D+00) THEN
C           WE CAN REMOVE THIS TEST BECAUSE MULTIPLYING BY 1.0
C             DOES NOT CHANGE THE ANSWER FOR HF CASE.
C
      ELSE IF(UROHF) THEN
        DO M=1,NINTMAX
          NIJ = NIJKLV(M,1)
          NKL = NIJKLV(M,2)
          NIK = NIJKLV(M,3)
          NIL = NIJKLV(M,4)
          NJK = NIJKLV(M,5)
          NJL = NIJKLV(M,6)
          VAL = VALV(M)
          VAL2 = VAL+VAL
          VAL2HFS = VAL2*HFSCAL
          VAL4 = VAL*4
          DUM = VAL4*(DA(NKL)+DB(NKL))
          FA(NIJ) = FA(NIJ)+DUM
          FB(NIJ) = FB(NIJ)+DUM
          DUM = VAL4*(DA(NIJ)+DB(NIJ))
          FA(NKL) = FA(NKL)+DUM
          FB(NKL) = FB(NKL)+DUM
          FA(NIK) = FA(NIK)-VAL2HFS*DA(NJL)
          FB(NIK) = FB(NIK)-VAL2HFS*DB(NJL)
          FA(NIL) = FA(NIL)-VAL2HFS*DA(NJK)
          FB(NIL) = FB(NIL)-VAL2HFS*DB(NJK)
          FA(NJK) = FA(NJK)-VAL2HFS*DA(NIL)
          FB(NJK) = FB(NJK)-VAL2HFS*DB(NIL)
          FA(NJL) = FA(NJL)-VAL2HFS*DA(NIK)
          FB(NJL) = FB(NJL)-VAL2HFS*DB(NIK)
        END DO
C
C                    GVB
C
      ELSE IF(SCFTYP.EQ.GVB) THEN
        IF(HFSCAL .EQ. 1.0D+00) THEN
          DO M=1,NINTMAX
            NIJ = NIJKLV(M,1)
            NKL = NIJKLV(M,2)
            NIK = NIJKLV(M,3)
            NIL = NIJKLV(M,4)
            NJK = NIJKLV(M,5)
            NJL = NIJKLV(M,6)
            VAL = VALV(M)
            VAL4 = VAL*4
            IF(NCO.GT.0) THEN
               FA(NIJ) = FA(NIJ)+VAL4*DA(NKL)
               FA(NKL) = FA(NKL)+VAL4*DA(NIJ)
               FA(NIK) = FA(NIK)-VAL *DA(NJL)
               FA(NIL) = FA(NIL)-VAL *DA(NJK)
               FA(NJK) = FA(NJK)-VAL *DA(NIL)
               FA(NJL) = FA(NJL)-VAL *DA(NIK)
               IOFF1 = L2
               IOFF2 = L2
            ELSE
               IOFF1 = 0
               IOFF2 = 0
            END IF
            DO IFO = 1,NSHL
               FA(NIJ+IOFF1)    = FA(NIJ+IOFF1)
     *                     + VAL4*DA(NKL+IOFF2)
               FA(NKL+IOFF1)    = FA(NKL+IOFF1)
     *                     + VAL4*DA(NIJ+IOFF2)
               FA(NIK+IOFF1+L2) = FA(NIK+IOFF1+L2)
     *                     + VAL *DA(NJL+IOFF2)
               FA(NIL+IOFF1+L2) = FA(NIL+IOFF1+L2)
     *                     + VAL *DA(NJK+IOFF2)
               FA(NJK+IOFF1+L2) = FA(NJK+IOFF1+L2)
     *                     + VAL *DA(NIL+IOFF2)
               FA(NJL+IOFF1+L2) = FA(NJL+IOFF1+L2)
     *                     + VAL *DA(NIK+IOFF2)
               IOFF1 = IOFF1+L2+L2
               IOFF2 = IOFF2+L2
            END DO
          END DO
        ELSE
           WRITE(IW,*)'DIRECT GVB DFT NOT IMPLEMENTED'
           CALL ABRT()
C          I COULDN'T FIGURE OUT THE DO 220 LOOP  :-)
        END IF
C
      END IF
      RETURN
      END
C
C*MODULE VECTOR  *DECK GENRAL       (INT2A)
      SUBROUTINE GENRAL(GHONDO,DDIJ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION GHONDO(*),DDIJ(*)
C
      LOGICAL IANDJ,KANDL,SAME,OUT,NORM,DOUBLE,KMAXCHK
C
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      COMMON /DENS  / DKL(225),DIJ(225)
      COMMON /IJGNRL/ AA(MXG2),R(MXG2),X1(MXG2),Y1(MXG2),Z1(MXG2),
     *                IJD(225)
      COMMON /INTDEX/ IJGT(225),IJX(225),IJY(225),IJZ(225),IK(225),
     *                KLGT(225),KLX(225),KLY(225),KLZ(225)
      COMMON /MISC  / IANDJ,KANDL,SAME
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /SETINT/ IN(9),KN(9),NI,NJ,NK,NL,NMAX,MMAX,
     +                BP01,B00,B10,XCP00,XC00,YCP00,YC00,ZCP00,ZC00,F00,
     +                DXIJ,DYIJ,DZIJ,DXKL,DYKL,DZKL
      COMMON /SHLINF/ AG(MXGSH),CSA(MXGSH),CPA(MXGSH),CDA(MXGSH),
     *                CFA(MXGSH),CGA(MXGSH),
     *                BG(MXGSH),CSB(MXGSH),CPB(MXGSH),CDB(MXGSH),
     *                CFB(MXGSH),CGB(MXGSH),
     *                CG(MXGSH),CSC(MXGSH),CPC(MXGSH),CDC(MXGSH),
     *                CFC(MXGSH),CGC(MXGSH),
     *                DG(MXGSH),CSD(MXGSH),CPD(MXGSH),CDD(MXGSH),
     *                CFD(MXGSH),CGD(MXGSH),
     *                XI,YI,ZI,XJ,YJ,ZJ,RRI,XK,YK,ZK,XL,YL,ZL,RRK,
     *                NGA,NGB,NGC,NGD
      COMMON /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     +                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     +                NIJ,IJ,KL,IJKL
      COMMON /SHLT  / TOL,CUTOFF,ICOUNT,OUT
C
      DIMENSION IN1(9)
      DIMENSION MXV(225*(NIJ+1)*(IJ+1))
      DIMENSION MYV(225*(NIJ+1)*(IJ+1))
      DIMENSION MZV(225*(NIJ+1)*(IJ+1))
      DIMENSION NINDV(225*(NIJ+1)*(IJ+1))
      DIMENSION IJDV(225*(NIJ+1)*(IJ+1))
      DIMENSION DKLV(225*(NIJ+1)*(IJ+1))
      DIMENSION NIJV(NIJ+1),NNV(NIJ+1),UV(NIJ+1,9),WV(NIJ+1,9)
      DIMENSION XINTV(NIJ+1,625,9),YINTV(NIJ+1,625,9),ZINTV(NIJ+1,625,9)
C
      PARAMETER (SQRT3=1.73205080756888D+00, SQRT5=2.23606797749979D+00,
     *           SQRT7=2.64575131106459D+00, PI252=34.986836655250D+00,
     *           ZERO=0.0D+00, HALF=0.5D+00, ONE=1.0D+00)
C
C     GENERAL INTEGRAL ROUTINE FOR SPD FUNCTIONS
C
      FACTOR = PI252*QQ4
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
      NI = LIT-1
      NJ = LJT-1
      NK = LKT-1
      NL = LLT-1
      DXIJ = XI-XJ
      DYIJ = YI-YJ
      DZIJ = ZI-ZJ
      DXKL = XK-XL
      DYKL = YK-YL
      DZKL = ZK-ZL
      NMAX = NI+NJ
      MMAX = NK+NL
      MAX = NMAX+1
      DO 100 I = 1,MAX
         N = I-1
         IF (N .LE. NI) IN1(I) = 125*N+1
         IF (N .GT. NI) IN1(I) = 125*NI+25*(N-NI)+1
  100 CONTINUE
      MAX = MMAX+1
      DO 120 K = 1,MAX
         N = K-1
         IF (N .LE. NK) KN(K) = 5*N
         IF (N .GT. NK) KN(K) = 5*NK+N-NK
  120 CONTINUE
C
C     STORE THE INDICES NEEDED FOR XINT/YINT/ZINT HERE...  THEY
C     ARE INVARIANT TO KG OR LG LOOPS, AND WERE ORIGINALLY
C     COMPUTED IN FORMS.  MOVING THEM HERE ELIMINATES A GREAT
C     DEAL OF UNNECESSARY COMPUTATION.
C     THE K LOOP VECTORIZES.
C
      KMAX = IK(IJ)
      NCNT = 0
      IF(NROOTS .GT. 2 .AND.
     *   IK(1).EQ.KMAX  .AND.  IJ.GT.KMAX) THEN
        KMAXCHK = .TRUE.
      ELSE
        KMAXCHK = .FALSE.
      END IF
C
      DO I = 1,IJ
        IF (.NOT.KMAXCHK) KMAX = IK(I)
        DO K = 1,KMAX
          NCNT = NCNT + 1
          MXV(NCNT)   = IJX(I)+KLX(K)
          MYV(NCNT)   = IJY(I)+KLY(K)
          MZV(NCNT)   = IJZ(I)+KLZ(K)
          IJDV(NCNT)  = IJD(I)
          NINDV(NCNT) = IJGT(I)+KLGT(K)
        END DO
      END DO
C
C     ----- K PRIMITIVE
C
      LGMAX = NGD
      DO 790 KG = 1,NGC
        AK = CG(KG)
        BRRK = AK*RRK
        AKXK = AK*XK
        AKYK = AK*YK
        AKZK = AK*ZK
C
C        ----- L PRIMITIVE
C
        IF (KANDL) LGMAX = KG
        DO 690 LG = 1,LGMAX
          AL = DG(LG)
          B = AK+AL
          BINV = ONE/B
          BBRRK = AL*BRRK*BINV
          IF (BBRRK .LE. TOL) THEN
            XB = (AKXK+AL*XL)*BINV
            YB = (AKYK+AL*YL)*BINV
            ZB = (AKZK+AL*ZL)*BINV
            BXBK = B*(XB-XK)
            BYBK = B*(YB-YK)
            BZBK = B*(ZB-ZK)
            BXBI = B*(XB-XI)
            BYBI = B*(YB-YI)
            BZBI = B*(ZB-ZI)
C
C           ----- DENSITY FACTOR
C
            DOUBLE=KANDL.AND.KG.NE.LG
            N = 0
            MAX = MAXL
            DUM1 = ZERO
            DUM2 = ZERO
            FBINV = FACTOR*BINV
C
C           THIS LARGE SERIES OF IF-TESTS WAS CONVERTED INTO
C           IF-ELSES FOR COMPILER PERFORMANCE AND EASIER
C           INTERPRETATION.
C
            DO 370 K = MINK,MAXK
               IF ( K .EQ. 1 ) THEN
                 DUM1 = CSC(KG)*FBINV
               ELSE IF ( K .EQ. 2 ) THEN
                 DUM1 = CPC(KG)*FBINV
               ELSE IF ( K .EQ. 5 ) THEN
                 DUM1 = CDC(KG)*FBINV
               ELSE IF ( K .EQ. 11 ) THEN
                 DUM1 = CFC(KG)*FBINV
               ELSE IF ( K .EQ. 21 ) THEN
                 DUM1 = CGC(KG)*FBINV
               ELSE IF (NORM) THEN
                 IF ( K .EQ. 8 ) THEN
                   DUM1 = DUM1*SQRT3
                 ELSE IF ( K .EQ. 14 ) THEN
                   DUM1 = DUM1*SQRT5
                 ELSE IF ( K .EQ. 20 ) THEN
                   DUM1 = DUM1*SQRT3
                 ELSE IF ( K .EQ. 24 ) THEN
                   DUM1 = DUM1*SQRT7
                 ELSE IF ( K .EQ. 30 ) THEN
                   DUM1 = DUM1*SQRT5/SQRT3
                 ELSE IF ( K .EQ. 33 ) THEN
                   DUM1 = DUM1*SQRT3
                 END IF
               END IF
               IF (KANDL) MAX = K
               DO 360 L = MINL,MAX
                  IF ( L .EQ. 1 ) THEN
                    DUM2 = DUM1*CSD(LG)
                    IF ( DOUBLE ) THEN
                      IF (K .LE. 1) THEN
                        DUM2 = DUM2+DUM2
                      ELSE
                        DUM2 = DUM2+CSC(KG)*CPD(LG)*FBINV
                      END IF
                    END IF
                  ELSE IF ( L .EQ. 2 ) THEN
                    DUM2 = DUM1*CPD(LG)
                    IF (DOUBLE) DUM2 = DUM2+DUM2
                  ELSE IF ( L .EQ. 5 ) THEN
                    DUM2 = DUM1*CDD(LG)
                    IF (DOUBLE) DUM2 = DUM2+DUM2
                  ELSE IF ( L .EQ. 11 ) THEN
                    DUM2 = DUM1*CFD(LG)
                    IF (DOUBLE) DUM2 = DUM2+DUM2
                  ELSE IF ( L .EQ. 21 ) THEN
                    DUM2 = DUM1*CGD(LG)
                    IF (DOUBLE) DUM2 = DUM2+DUM2
                  ELSE IF ( NORM ) THEN
                    IF ( L .EQ. 8 ) THEN
                      DUM2 = DUM2*SQRT3
                    ELSE IF ( L .EQ. 14 ) THEN
                      DUM2 = DUM2*SQRT5
                    ELSE IF ( L .EQ. 20 ) THEN
                      DUM2 = DUM2*SQRT3
                    ELSE IF ( L .EQ. 24 ) THEN
                      DUM2 = DUM2*SQRT7
                    ELSE IF ( L .EQ. 30 ) THEN
                      DUM2 = DUM2*SQRT5/SQRT3
                    ELSE IF ( L .EQ. 33 ) THEN
                      DUM2 = DUM2*SQRT3
                    END IF
                  END IF
                  N = N+1
                  DKL(N) = DUM2
  360          CONTINUE
  370       CONTINUE
          END IF
C
C     CONVERT DKL TO DKLV FOR USE IN FORMSV.  K LOOP VECTORIZES.
C
          KMAX = IK(IJ)
          NCNT = 0
          DO 390 I = 1,IJ
            IF (.NOT.KMAXCHK) KMAX = IK(I)
            DO 380 K = 1,KMAX
              NCNT = NCNT + 1
              DKLV(NCNT) = DKL(K)
  380       CONTINUE
  390     CONTINUE
C
C         ----- PAIR OF I,J PRIMITIVES
C         ONCE AGAIN, FIGURE OUT HOW MANY LOOPS TO DO FIRST.
C         THIS VECTORIZES.
C
          NN = 0
          NIJTOTAL = 0
          DO 400 N = 1,NIJ
            DUM = BBRRK+R(N)
            IF (DUM .LE. TOL) THEN
              NIJTOTAL = NIJTOTAL + 1
              NIJV(NIJTOTAL) = N
              NNV(NIJTOTAL) = NN
            END IF
  400     NN = NN+16
C
C         CALL SPECIAL VECTOR VERSION OF RT123.  BY SPLITTING THE
C         NIJ LOOP WE CAN ADD AN EXTRA DIMENSION TO U AND W AND
C         VECTORIZE OVER THIS LOOP.
C
          IF (NROOTS .LE. 3) THEN
            DO 450 NNEW = 1,NIJTOTAL
              N = NIJV(NNEW)
              NN = NNV(NNEW)
              DUM = BBRRK+R(N)
              A = AA(N)
              AB = A*B
              AANDB = A+B
              EXPE = EXP(-DUM)/SQRT(AANDB)
              RHO = AB/AANDB
              XA = X1(N)
              YA = Y1(N)
              ZA = Z1(N)
              XX = RHO*((XA-XB)*(XA-XB) + (YA-YB)*(YA-YB)
     *                                + (ZA-ZB)*(ZA-ZB))
C
C              ----- ROOTS AND WEIGHTS FOR QUADRATURE
C
              CALL RT123V(XX,UV,WV,NROOTS,NIJ,NNEW)
  450       CONTINUE
C
C         ELSE CALL THE REGULAR VERSIONS OF ROOT4/5/6.
C         CALLED MUCH LESS OFTEN... EFFORT TO VECTORIZE FAR EXCEEDS
C          THE COMPUTATIONAL GAIN.  LEFT IT ALONE.
C
          ELSE
            DO 490 NNEW = 1,NIJTOTAL
              N = NIJV(NNEW)
              NN = NNV(NNEW)
              DUM = BBRRK+R(N)
              A = AA(N)
              AB = A*B
              AANDB = A+B
              EXPE = EXP(-DUM)/SQRT(AANDB)
              RHO = AB/AANDB
              XA = X1(N)
              YA = Y1(N)
              ZA = Z1(N)
              XX = RHO*((XA-XB)*(XA-XB) + (YA-YB)*(YA-YB)
     *                                + (ZA-ZB)*(ZA-ZB))
C
C              ----- ROOTS AND WEIGHTS FOR QUADRATURE
C
              IF (NROOTS .EQ. 4) THEN
                CALL ROOT4
              ELSE IF (NROOTS .EQ. 5) THEN
                CALL ROOT5
              ELSE
                CALL ROOT6
              END IF
              DO 480 M = 1, NROOTS
                UV(NNEW,M) = U(M)
                WV(NNEW,M) = W(M)
  480         CONTINUE
  490       CONTINUE
          END IF
C
          DO 590 NNEW = 1,NIJTOTAL
            N = NIJV(NNEW)
            NN = NNV(NNEW)
            DUM = BBRRK+R(N)
            DO 500 I = 1,IJ
               DIJ(I) = DDIJ(IJD(I)+NN)
  500       CONTINUE
            A = AA(N)
            AB = A*B
            AANDB = A+B
            EXPE = EXP(-DUM)/SQRT(AANDB)
            RHO = AB/AANDB
            XA = X1(N)
            YA = Y1(N)
            ZA = Z1(N)
            XX = RHO*((XA-XB)*(XA-XB) + (YA-YB)*(YA-YB)
     *                                + (ZA-ZB)*(ZA-ZB))
            AXAK = A*(XA-XK)
            AYAK = A*(YA-YK)
            AZAK = A*(ZA-ZK)
            AXAI = A*(XA-XI)
            AYAI = A*(YA-YI)
            AZAI = A*(ZA-ZI)
            C1X = BXBK+AXAK
            C2X = A*BXBK
            C3X = BXBI+AXAI
            C4X = B*AXAI
            C1Y = BYBK+AYAK
            C2Y = A*BYBK
            C3Y = BYBI+AYAI
            C4Y = B*AYAI
            C1Z = BZBK+AZAK
            C2Z = A*BZBK
            C3Z = BZBI+AZAI
            C4Z = B*AZAI
            MM = 0
            MAX = NMAX+1
C
C              COMPUTE TWO-ELECTRON INTEGRALS FOR EACH ROOT
C
            DO 520 M = 1,NROOTS
               U2 = UV(NNEW,M)*RHO
               F00 = EXPE*WV(NNEW,M)
               DO 510 I = 1,MAX
                  IN(I) = IN1(I)+MM
  510          CONTINUE
               DUMINV = ONE/(AB+U2*AANDB)
               DM2INV = HALF*DUMINV
               BP01 = (A+U2)*DM2INV
               B00 = U2*DM2INV
               B10 = (B+U2)*DM2INV
               XCP00 = (U2*C1X+C2X)*DUMINV
               XC00 = (U2*C3X+C4X)*DUMINV
               YCP00 = (U2*C1Y+C2Y)*DUMINV
               YC00 = (U2*C3Y+C4Y)*DUMINV
               ZCP00 = (U2*C1Z+C2Z)*DUMINV
               ZC00 = (U2*C3Z+C4Z)*DUMINV
               CALL XYZINTV(NIJ,NNEW,XINTV,YINTV,ZINTV)
               MM = MM+625
  520       CONTINUE
  590     CONTINUE
C
C              ----- FORM (I,J//K,L) INTEGRALS OVER FUNCTIONS
C
CDIR$ NOINLINE
          CALL FORMSV(GHONDO,DDIJ,XINTV,YINTV,ZINTV,
     *           NIJTOTAL,NNV,NIJ,225*NIJ*IJ,NROOTS,
     *           MXV,MYV,MZV,IJDV,NINDV,DKLV,NCNT )
  690   CONTINUE
  790 CONTINUE
C
      RETURN
      END
C*MODULE VECTOR  *DECK FORMSV   CALLED ONLY FROM WITHIN THIS FILE
      SUBROUTINE FORMSV(GHONDO,DDIJ,XINTV,YINTV,ZINTV,
     *  NIJTOTAL,NNV,NIJ,IJDIM,NROOTS,
     *  MXV,MYV,MZV,IJDV,NINDV,DKLV,NCNT )
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00)
C
      DIMENSION GHONDO(*),DDIJ(*)
      DIMENSION MXV(IJDIM),MYV(IJDIM),MZV(IJDIM),NNV(NIJ)
      DIMENSION IJDV(IJDIM),DKLV(IJDIM),GGV(IJDIM),NINDV(IJDIM)
      DIMENSION XINTV(NIJ+1,625,9),YINTV(NIJ+1,625,9),ZINTV(NIJ+1,625,9)
C
C     ----- FORM INTEGRALS OVER FUNCTIONS
C
      GGV(1:NCNT) = ZERO
C
C     THE DO N=1,NIJ LOOP HAS BEEN PULLED IN FROM GENRAL, BUT
C       ONLY AFTER WE DETERMINE THE HOW MANY LOOPS WE ACTUALLY
C       NEED TO DO (NIJTOTAL), WHICH IS FIGURED OUT IN GENRAL.
C     WE BUILD MX,MY,MZ AT THE BEGINNING OF GENRAL NOW, SINCE
C       THEIR VALUES ARE INVARIANT TO LOOPS THROUGH GENRAL; AT
C       THE SAME TIME WE DETERMINE NCNT, IJDV, AND NINDV.
C     DKLV IS BUILT FROM DKL IN GENRAL.
C     AS A RESULT OF ALL THIS, WE CAN NOW VECTORIZE CLEANLY OVER
C       NCNT.
C
      DO 700 NNEW = 1,NIJTOTAL
        NN = NNV(NNEW)
        DO 500 NC = 1,NCNT
          MX = MXV(NC)
          MY = MYV(NC)
          MZ = MZV(NC)
          DUMGG = ZERO
          DO L = 1, NROOTS
            DUMGG = DUMGG + XINTV(NNEW,MX,L)
     *              *YINTV(NNEW,MY,L)*ZINTV(NNEW,MZ,L)
          END DO
          GGV(NC) = GGV(NC) + DUMGG*DDIJ(IJDV(NC)+NN)*DKLV(NC)
  500   CONTINUE
  700 CONTINUE
C
C     STORE VALUES INTO GHONDO
C     REPEATED INDICES PREVENTS THIS FROM VECTORIZING, THUS
C      WE SPLIT IT OFF FROM THE MAIN COMPUTATION LOOP ABOVE.
C     NO NEED TO LOOP OVER NIJTOTAL HERE, THE UPDATES TO GHONDO
C      ARE NC-DEPENDENT.  THE GGV ARRAY WILL CONTAIN THE APPROPRIATE
C      SUMMATION FOR EACH VALUE OF NC.
C
      DO 900 NC = 1,NCNT
        NIND = NINDV(NC)
        GHONDO(NIND) = GHONDO(NIND) + GGV(NC)
  900 CONTINUE
C
      RETURN
      END
C*MODULE VECTOR  *DECK XYZINTV  CALLED ONLY FROM WITHIN THIS FILE
      SUBROUTINE XYZINTV(NIJ,NNEW,XINTV,YINTV,ZINTV)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL N0,N1,M0,M1,FIRST1,FIRST2,FIRST3,FIRST4
C
      COMMON /SETINT/ I(9),K(9),NIMAX,NJMAX,NKMAX,NLMAX,NMAX,MMAX
     +               ,BP01,B00,B10,XCP00,XC00,YCP00,YC00,ZCP00,ZC00,F00
     +               ,DXIJ,DYIJ,DZIJ,DXKL,DYKL,DZKL
C
      DIMENSION XINTV(NIJ+1,5625),YINTV(NIJ+1,5625),ZINTV(NIJ+1,5625)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
      N0 = NMAX .EQ. 0
      N1 = NMAX .LE. 1
      M0 = MMAX .EQ. 0
      M1 = MMAX .LE. 1
C
C     ----- I(0,0) -----
C
      I1 = I(1)
      XINTV(NNEW,I1) = ONE
      YINTV(NNEW,I1) = ONE
      ZINTV(NNEW,I1) = F00
      IF (N0 .AND. M0) RETURN
      I2 = I(2)
      K2 = K(2)
      CP10 = B00
      IF (.NOT. N0) THEN
C
C     ----- I(1,0) -----
C
        XINTV(NNEW,I2) = XC00
        YINTV(NNEW,I2) = YC00
        ZINTV(NNEW,I2) = ZC00*F00
        IF (M0) GO TO 120
      END IF
C
C     ----- I(0,1) -----
C
      I3 = I1+K2
      XINTV(NNEW,I3) = XCP00
      YINTV(NNEW,I3) = YCP00
      ZINTV(NNEW,I3) = ZCP00*F00
      IF (.NOT. N0) THEN
C
C     ----- I(1,1) -----
C
        I3 = I2+K2
        XINTV(NNEW,I3) = XCP00*XINTV(NNEW,I2)+CP10
        YINTV(NNEW,I3) = YCP00*YINTV(NNEW,I2)+CP10
        ZINTV(NNEW,I3) = ZCP00*ZINTV(NNEW,I2)+CP10*F00
      END IF
  120 IF (.NOT. N1) THEN
        C10 = ZERO
        I3 = I1
        I4 = I2
        DO 160 N = 2,NMAX
          C10 = C10+B10
C
C     ----- I(N,0) -----
C
          I5 = I(N+1)
          XINTV(NNEW,I5) = C10*XINTV(NNEW,I3)+XC00*XINTV(NNEW,I4)
          YINTV(NNEW,I5) = C10*YINTV(NNEW,I3)+YC00*YINTV(NNEW,I4)
          ZINTV(NNEW,I5) = C10*ZINTV(NNEW,I3)+ZC00*ZINTV(NNEW,I4)
          IF ( .NOT. M0) THEN
            CP10 = CP10+B00
C
C     ----- I(N,1) -----
C
            I3 = I5+K2
            XINTV(NNEW,I3) = XCP00*XINTV(NNEW,I5)+CP10*XINTV(NNEW,I4)
            YINTV(NNEW,I3) = YCP00*YINTV(NNEW,I5)+CP10*YINTV(NNEW,I4)
            ZINTV(NNEW,I3) = ZCP00*ZINTV(NNEW,I5)+CP10*ZINTV(NNEW,I4)
          END IF
          I3 = I4
  160   I4 = I5
      END IF
      IF ( .NOT. M1) THEN
        CP01 = ZERO
        C01 = B00
        I3 = I1
        I4 = I1+K2
        DO 220 M = 2,MMAX
          CP01 = CP01+BP01
C
C     ----- I(0,M) -----
C
          I5 = I1+K(M+1)
          XINTV(NNEW,I5) = CP01*XINTV(NNEW,I3)+XCP00*XINTV(NNEW,I4)
          YINTV(NNEW,I5) = CP01*YINTV(NNEW,I3)+YCP00*YINTV(NNEW,I4)
          ZINTV(NNEW,I5) = CP01*ZINTV(NNEW,I3)+ZCP00*ZINTV(NNEW,I4)
          IF (.NOT. N0) THEN
            C01 = C01+B00
C
C     ----- I(1,M) -----
C
            I3 = I2+K(M+1)
            XINTV(NNEW,I3) = XC00*XINTV(NNEW,I5)+C01*XINTV(NNEW,I4)
            YINTV(NNEW,I3) = YC00*YINTV(NNEW,I5)+C01*YINTV(NNEW,I4)
            ZINTV(NNEW,I3) = ZC00*ZINTV(NNEW,I5)+C01*ZINTV(NNEW,I4)
          END IF
          I3 = I4
  220   I4 = I5
      END IF
      IF (.NOT. N1 .AND. .NOT. M1) THEN
C
C     ----- I(N,M) -----
C
        C01 = B00
        K3 = K2
        DO 280 M = 2,MMAX
          K4 = K(M+1)
          C01 = C01+B00
          I3 = I1
          I4 = I2
          C10 = B10
          DO 260 N = 2,NMAX
            I5 = I(N+1)
            XINTV(NNEW,I5+K4) = C10*XINTV(NNEW,I3+K4)
     *           +XC00*XINTV(NNEW,I4+K4)+C01*XINTV(NNEW,I4+K3)
            YINTV(NNEW,I5+K4) = C10*YINTV(NNEW,I3+K4)
     *           +YC00*YINTV(NNEW,I4+K4)+C01*YINTV(NNEW,I4+K3)
            ZINTV(NNEW,I5+K4) = C10*ZINTV(NNEW,I3+K4)
     *           +ZC00*ZINTV(NNEW,I4+K4)+C01*ZINTV(NNEW,I4+K3)
            C10 = C10+B10
            I3 = I4
  260     I4 = I5
  280   K3 = K4
      END IF
      IF (NJMAX .GT. 0) THEN
C
C     ----- I(NI,NJ,M) -----
C
        M = 0
        I5 = I(NMAX+1)
        FIRST1 = .TRUE.
        DO 430 WHILE (FIRST1 .OR. M .LE. MMAX)
          MIN = NIMAX
          KM = K(M+1)
          FIRST2 = .TRUE.
          DO 360 WHILE (FIRST2 .OR. MIN .LT. NMAX)
            N = NMAX
            I3 = I5+KM
            FIRST3 = .TRUE.
            DO 340 WHILE (FIRST3 .OR. N .GT. MIN)
              I4 = I(N)+KM
              XINTV(NNEW,I3) = XINTV(NNEW,I3)+DXIJ*XINTV(NNEW,I4)
              YINTV(NNEW,I3) = YINTV(NNEW,I3)+DYIJ*YINTV(NNEW,I4)
              ZINTV(NNEW,I3) = ZINTV(NNEW,I3)+DZIJ*ZINTV(NNEW,I4)
              I3 = I4
              N = N-1
              FIRST3 = .FALSE.
  340       END DO
            MIN = MIN+1
            FIRST2 = .FALSE.
  360     END DO
          IF (NIMAX .GT. 0) THEN
            I3 = 25+KM+I1
            DO 400 NJ = 1,NJMAX
              I4 = I3
              DO 380 NI = 1,NIMAX
                XINTV(NNEW,I4)=XINTV(NNEW,I4+100)+DXIJ*XINTV(NNEW,I4-25)
                YINTV(NNEW,I4)=YINTV(NNEW,I4+100)+DYIJ*YINTV(NNEW,I4-25)
                ZINTV(NNEW,I4)=ZINTV(NNEW,I4+100)+DZIJ*ZINTV(NNEW,I4-25)
  380         I4 = I4+125
  400       I3 = I3+25
          END IF
          M = M+1
          FIRST1 = .FALSE.
  430   END DO
      END IF
      IF (NLMAX .GT. 0) THEN
C
C     ----- I(NI,NJ,NK,NL) -----
C
        I5 = K(MMAX+1)
        IA = I1
        NI = 0
        FIRST4 = .TRUE.
        DO 580 WHILE (FIRST4 .OR. NI .LE. NIMAX)
          NJ = 0
          IB = IA
          FIRST1 = .TRUE.
          DO 570 WHILE (FIRST1 .OR. NJ .LE. NJMAX)
            MIN = NKMAX
            FIRST2 = .TRUE.
            DO 530 WHILE (FIRST2 .OR. MIN .LT. MMAX)
              M = MMAX
              I3 = IB+I5
              FIRST3 = .TRUE.
              DO 520 WHILE (FIRST3 .OR. M .GT. MIN)
                I4 = IB+K(M)
                XINTV(NNEW,I3) = XINTV(NNEW,I3)+DXKL*XINTV(NNEW,I4)
                YINTV(NNEW,I3) = YINTV(NNEW,I3)+DYKL*YINTV(NNEW,I4)
                ZINTV(NNEW,I3) = ZINTV(NNEW,I3)+DZKL*ZINTV(NNEW,I4)
                I3 = I4
                M = M-1
                FIRST3 = .FALSE.
  520         END DO
              MIN = MIN+1
              FIRST2 = .FALSE.
  530       END DO
            IF (NKMAX .GT. 0) THEN
              I3 = IB+1
              DO 560 NL = 1,NLMAX
                I4 = I3
                DO 540 NK = 1,NKMAX
                  XINTV(NNEW,I4)= XINTV(NNEW,I4+4)+DXKL*XINTV(NNEW,I4-1)
                  YINTV(NNEW,I4)= YINTV(NNEW,I4+4)+DYKL*YINTV(NNEW,I4-1)
                  ZINTV(NNEW,I4)= ZINTV(NNEW,I4+4)+DZKL*ZINTV(NNEW,I4-1)
                  I4 = I4+5
  540           END DO
              I3 = I3+1
  560         END DO
            END IF
            NJ = NJ+1
            IB = IB+25
            FIRST1 = .FALSE.
  570     END DO
          NI = NI+1
          IA = IA+125
          FIRST4 = .FALSE.
  580   END DO
      END IF
C
      RETURN
      END
C*MODULE RYSPOL  *DECK RT123V   CALLED ONLY FROM WITHIN THIS FILE
C
      SUBROUTINE RT123V(X,UV,WV,NROOTS,NIJ,NN)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C             *****   VERSION FEBRUARY 13,1975   *****
      DIMENSION UV(NIJ+1,9),WV(NIJ+1,9)
      DATA R12,PIE4/2.75255128608411D-01, 7.85398163397448D-01/
      DATA R22,W22/ 2.72474487139158D+00, 9.17517095361369D-02/
      DATA R13/     1.90163509193487D-01/
      DATA R23,W23/ 1.78449274854325D+00, 1.77231492083829D-01/
      DATA R33,W33/ 5.52534374226326D+00, 5.11156880411248D-03/
C
CDIR$ INLINEALWAYS RT123V
C
      IF (X .GT. 5.0D+00) GO TO 400
      IF (X .GT. 1.0D+00) GO TO 280
      IF (X .GT. 3.0D-07) GO TO 180
C
C     X IS APPROXIMATELY ZERO.         NROOTS=1,2, OR 3
C
      IF (NROOTS == 1) THEN
        UV(NN,1) = 0.5D+00 -X/5.0D+00
        WV(NN,1) = 1.0D+00 -X/3.0D+00
      ELSE IF (NROOTS == 2) THEN
        UV(NN,1) = 1.30693606237085D-01 -2.90430236082028D-02 *X
        UV(NN,2) = 2.86930639376291D+00 -6.37623643058102D-01 *X
        WV(NN,1) = 6.52145154862545D-01 -1.22713621927067D-01 *X
        WV(NN,2) = 3.47854845137453D-01 -2.10619711404725D-01 *X
      ELSE  ! NROOTS == 3
        UV(NN,1) = 6.03769246832797D-02 -9.28875764357368D-03 *X
        UV(NN,2) = 7.76823355931043D-01 -1.19511285527878D-01 *X
        UV(NN,3) = 6.66279971938567D+00 -1.02504611068957D+00 *X
        WV(NN,1) = 4.67913934572691D-01 -5.64876917232519D-02 *X
        WV(NN,2) = 3.60761573048137D-01 -1.49077186455208D-01 *X
        WV(NN,3) = 1.71324492379169D-01 -1.27768455150979D-01 *X
      END IF
      RETURN
C
C     X = 0.0 TO 1.0                   NROOTS=1,2, OR 3
C
  180 IF (NROOTS < 3) THEN
        F1 = ((((((((-8.36313918003957D-08*X+
     +    1.21222603512827D-06 )*X-
     +    1.15662609053481D-05 )*X+9.25197374512647D-05 )*X-
     +    6.40994113129432D-04 )*X+3.78787044215009D-03 )*X-
     +    1.85185172458485D-02 )*X+7.14285713298222D-02 )*X-
     +    1.99999999997023D-01 )*X+3.33333333333318D-01
        WV(NN,1) = (X+X)*F1+EXP(-X)
        IF (NROOTS < 2) THEN
          UV(NN,1) = F1/(WV(NN,1)-F1)
        ELSE
          UV(NN,1) = (((((((-2.35234358048491D-09*X+
     +      2.49173650389842D-08)*X-
     +      4.558315364581D-08)*X-2.447252174587D-06)*X+
     +      4.743292959463D-05)*X-5.33184749432408D-04 )*X+
     +      4.44654947116579D-03 )*X-2.90430236084697D-02 )*X+
     +      1.30693606237085D-01
          UV(NN,2) = (((((((-2.47404902329170D-08*X+
     +      2.36809910635906D-07)*X+
     +      1.835367736310D-06)*X-2.066168802076D-05)*X-
     +      1.345693393936D-04)*X-5.88154362858038D-05 )*X+
     +      5.32735082098139D-02 )*X-6.37623643056745D-01 )*X+
     +      2.86930639376289D+00
          WV(NN,2) = ((F1-WV(NN,1))*UV(NN,1)+F1)*(1.0D+00+UV(NN,2))/
     +            (UV(NN,2)-UV(NN,1))
          WV(NN,1) = WV(NN,1)-WV(NN,2)
        END IF
      ELSE
        UV(NN,1) = ((((((-5.10186691538870D-10*X+
     +    2.40134415703450D-08)*X-
     +    5.01081057744427D-07 )*X+7.58291285499256D-06 )*X-
     +    9.55085533670919D-05 )*X+1.02893039315878D-03 )*X-
     +    9.28875764374337D-03 )*X+6.03769246832810D-02
        UV(NN,2) = ((((((-1.29646524960555D-08*X+
     +    7.74602292865683D-08)*X+
     +    1.56022811158727D-06 )*X-1.58051990661661D-05 )*X-
     +    3.30447806384059D-04 )*X+9.74266885190267D-03 )*X-
     +    1.19511285526388D-01 )*X+7.76823355931033D-01
        UV(NN,3) = ((((((-9.28536484109606D-09*X-
     +    3.02786290067014D-07)*X-
     +    2.50734477064200D-06 )*X-7.32728109752881D-06 )*X+
     +    2.44217481700129D-04 )*X+4.94758452357327D-02 )*X-
     +    1.02504611065774D+00 )*X+6.66279971938553D+00
        F2 = ((((((((-7.60911486098850D-08*X+
     +    1.09552870123182D-06 )*X-
     +    1.03463270693454D-05 )*X+8.16324851790106D-05 )*X-
     +    5.55526624875562D-04 )*X+3.20512054753924D-03 )*X-
     +    1.51515139838540D-02 )*X+5.55555554649585D-02 )*X-
     +    1.42857142854412D-01 )*X+1.99999999999986D-01
        E = EXP(-X)
        F1 = ((X+X)*F2+E)/3.0D+00
        WV(NN,1) = (X+X)*F1+E
        T1 = UV(NN,1)/(UV(NN,1)+1.0D+00)
        T2 = UV(NN,2)/(UV(NN,2)+1.0D+00)
        T3 = UV(NN,3)/(UV(NN,3)+1.0D+00)
        A2 = F2-T1*F1
        A1 = F1-T1*WV(NN,1)
        WV(NN,3) = (A2-T2*A1)/((T3-T2)*(T3-T1))
        WV(NN,2) = (T3*A1-A2)/((T3-T2)*(T2-T1))
        WV(NN,1) = WV(NN,1)-WV(NN,2)-WV(NN,3)
      END IF
      RETURN
  280 IF (X .GT. 3.0D+00) GO TO 340
C
C     X = 1.0 TO 3.0                   NROOTS=1,2, OR 3
C
      Y = X-2.0D+00
      IF (NROOTS < 3) THEN
        F1 = ((((((((((-1.61702782425558D-10*Y+1.96215250865776D-09 )*Y-
     +     2.14234468198419D-08 )*Y+2.17216556336318D-07 )*Y-
     +     1.98850171329371D-06 )*Y+1.62429321438911D-05 )*Y-
     +     1.16740298039895D-04 )*Y+7.24888732052332D-04 )*Y-
     +     3.79490003707156D-03 )*Y+1.61723488664661D-02 )*Y-
     +     5.29428148329736D-02 )*Y+1.15702180856167D-01
        WV(NN,1) = (X+X)*F1+EXP(-X)
        IF (NROOTS < 2) THEN
          UV(NN,1) = F1/(WV(NN,1)-F1)
          RETURN
        ELSE
          UV(NN,1) = (((((((((-6.36859636616415D-12*Y+
     +     8.47417064776270D-11)*Y-
     +     5.152207846962D-10)*Y-3.846389873308D-10)*Y+
     +     8.472253388380D-08)*Y-1.85306035634293D-06 )*Y+
     +     2.47191693238413D-05 )*Y-2.49018321709815D-04 )*Y+
     +     2.19173220020161D-03 )*Y-1.63329339286794D-02 )*Y+
     +     8.68085688285261D-02
          UV(NN,2) = ((((((((( 1.45331350488343D-10*Y+
     +     2.07111465297976D-09)*Y-
     +     1.878920917404D-08)*Y-1.725838516261D-07)*Y+
     +     2.247389642339D-06)*Y+9.76783813082564D-06 )*Y-
     +     1.93160765581969D-04 )*Y-1.58064140671893D-03 )*Y+
     +     4.85928174507904D-02 )*Y-4.30761584997596D-01 )*Y+
     +     1.80400974537950D+00
          WV(NN,2) = ((F1-WV(NN,1))*UV(NN,1)+F1)*(1.0D+00+UV(NN,2))/
     +            (UV(NN,2)-UV(NN,1))
          WV(NN,1) = WV(NN,1)-WV(NN,2)
          RETURN
        END IF
      ELSE
        UV(NN,1) = (((((((( 1.44687969563318D-12*Y+
     +     4.85300143926755D-12)*Y-
     +     6.55098264095516D-10 )*Y+1.56592951656828D-08 )*Y-
     +     2.60122498274734D-07 )*Y+3.86118485517386D-06 )*Y-
     +     5.13430986707889D-05 )*Y+6.03194524398109D-04 )*Y-
     +     6.11219349825090D-03 )*Y+4.52578254679079D-02
        UV(NN,2) = ((((((( 6.95964248788138D-10*Y-
     +     5.35281831445517D-09)*Y-
     +     6.745205954533D-08)*Y+1.502366784525D-06)*Y+
     +     9.923326947376D-07)*Y-3.89147469249594D-04 )*Y+
     +     7.51549330892401D-03 )*Y-8.48778120363400D-02 )*Y+
     +     5.73928229597613D-01
        UV(NN,3) = ((((((((-2.81496588401439D-10*Y+
     +     3.61058041895031D-09)*Y+
     +     4.53631789436255D-08 )*Y-1.40971837780847D-07 )*Y-
     +     6.05865557561067D-06 )*Y-5.15964042227127D-05 )*Y+
     +     3.34761560498171D-05 )*Y+5.04871005319119D-02 )*Y-
     +     8.24708946991557D-01 )*Y+4.81234667357205D+00
        F2 = ((((((((((-1.48044231072140D-10*Y+1.78157031325097D-09 )*Y-
     +     1.92514145088973D-08 )*Y+1.92804632038796D-07 )*Y-
     +     1.73806555021045D-06 )*Y+1.39195169625425D-05 )*Y-
     +     9.74574633246452D-05 )*Y+5.83701488646511D-04 )*Y-
     +     2.89955494844975D-03 )*Y+1.13847001113810D-02 )*Y-
     +     3.23446977320647D-02 )*Y+5.29428148329709D-02
        E = EXP(-X)
        F1 = ((X+X)*F2+E)/3.0D+00
        WV(NN,1) = (X+X)*F1+E
        T1 = UV(NN,1)/(UV(NN,1)+1.0D+00)
        T2 = UV(NN,2)/(UV(NN,2)+1.0D+00)
        T3 = UV(NN,3)/(UV(NN,3)+1.0D+00)
        A2 = F2-T1*F1
        A1 = F1-T1*WV(NN,1)
        WV(NN,3) = (A2-T2*A1)/((T3-T2)*(T3-T1))
        WV(NN,2) = (T3*A1-A2)/((T3-T2)*(T2-T1))
        WV(NN,1) = WV(NN,1)-WV(NN,2)-WV(NN,3)
      END IF
      RETURN
C
C     X = 3.0 TO 5.0                   NROOTS =1,2, OR 3
C
  340 Y = X-4.0D+00
      IF (NROOTS .EQ. 3) GO TO 380
      F1 = ((((((((((-2.62453564772299D-11*Y+3.24031041623823D-10 )*Y-
     +     3.614965656163D-09)*Y+3.760256799971D-08)*Y-
     +     3.553558319675D-07)*Y+3.022556449731D-06)*Y-
     +     2.290098979647D-05)*Y+1.526537461148D-04)*Y-
     +     8.81947375894379D-04 )*Y+4.33207949514611D-03 )*Y-
     +     1.75257821619926D-02 )*Y+5.28406320615584D-02
      WV(NN,1) = (X+X)*F1+EXP(-X)
      IF (NROOTS .EQ. 2) GO TO 360
      UV(NN,1) = F1/(WV(NN,1)-F1)
      RETURN
  360 UV(NN,1) = ((((((((-4.11560117487296D-12*Y+
     +     7.10910223886747D-11)*Y-
     +     1.73508862390291D-09 )*Y+5.93066856324744D-08 )*Y-
     +     9.76085576741771D-07 )*Y+1.08484384385679D-05 )*Y-
     +     1.12608004981982D-04 )*Y+1.16210907653515D-03 )*Y-
     +     9.89572595720351D-03 )*Y+6.12589701086408D-02
      UV(NN,2) = (((((((((-1.80555625241001D-10*Y+
     +     5.44072475994123D-10)*Y+
     +     1.603498045240D-08)*Y-1.497986283037D-07)*Y-
     +     7.017002532106D-07)*Y+1.85882653064034D-05 )*Y-
     +     2.04685420150802D-05 )*Y-2.49327728643089D-03 )*Y+
     +     3.56550690684281D-02 )*Y-2.60417417692375D-01 )*Y+
     +     1.12155283108289D+00
      WV(NN,2) = ((F1-WV(NN,1))*UV(NN,1)+F1)*(1.0D+00+UV(NN,2))/
     +        (UV(NN,2)-UV(NN,1))
      WV(NN,1) = WV(NN,1)-WV(NN,2)
      RETURN
  380 UV(NN,1) = ((((((( 1.44265709189601D-11*Y-
     +     4.66622033006074D-10)*Y+
     +     7.649155832025D-09)*Y-1.229940017368D-07)*Y+
     +     2.026002142457D-06)*Y-2.87048671521677D-05 )*Y+
     +     3.70326938096287D-04 )*Y-4.21006346373634D-03 )*Y+
     +     3.50898470729044D-02
      UV(NN,2) = ((((((((-2.65526039155651D-11*Y+
     +     1.97549041402552D-10)*Y+
     +     2.15971131403034D-09 )*Y-7.95045680685193D-08 )*Y+
     +     5.15021914287057D-07 )*Y+1.11788717230514D-05 )*Y-
     +     3.33739312603632D-04 )*Y+5.30601428208358D-03 )*Y-
     +     5.93483267268959D-02 )*Y+4.31180523260239D-01
      UV(NN,3) = ((((((((-3.92833750584041D-10*Y-
     +     4.16423229782280D-09)*Y+
     +     4.42413039572867D-08 )*Y+6.40574545989551D-07 )*Y-
     +     3.05512456576552D-06 )*Y-1.05296443527943D-04 )*Y-
     +     6.14120969315617D-04 )*Y+4.89665802767005D-02 )*Y-
     +     6.24498381002855D-01 )*Y+3.36412312243724D+00
      F2 = ((((((((((-2.36788772599074D-11*Y+2.89147476459092D-10 )*Y-
     +     3.18111322308846D-09 )*Y+3.25336816562485D-08 )*Y-
     +     3.00873821471489D-07 )*Y+2.48749160874431D-06 )*Y-
     +     1.81353179793672D-05 )*Y+1.14504948737066D-04 )*Y-
     +     6.10614987696677D-04 )*Y+2.64584212770942D-03 )*Y-
     +     8.66415899015349D-03 )*Y+1.75257821619922D-02
        E = EXP(-X)
        F1 = ((X+X)*F2+E)/3.0D+00
        WV(NN,1) = (X+X)*F1+E
        T1 = UV(NN,1)/(UV(NN,1)+1.0D+00)
        T2 = UV(NN,2)/(UV(NN,2)+1.0D+00)
        T3 = UV(NN,3)/(UV(NN,3)+1.0D+00)
        A2 = F2-T1*F1
        A1 = F1-T1*WV(NN,1)
        WV(NN,3) = (A2-T2*A1)/((T3-T2)*(T3-T1))
        WV(NN,2) = (T3*A1-A2)/((T3-T2)*(T2-T1))
        WV(NN,1) = WV(NN,1)-WV(NN,2)-WV(NN,3)
      RETURN
  400 IF (X .GT. 15.0D+00) GO TO 560
      E = EXP(-X)
      IF (X .GT. 10.0D+00) GO TO 480
C
C     X = 5.0 TO 10.0                  NROOTS =1,2, OR 3
C
      WV(NN,1) = (((((( 4.6897511375022D-01/X-6.9955602298985D-01)/X +
     +     5.3689283271887D-01)/X-3.2883030418398D-01)/X +
     +     2.4645596956002D-01)/X-4.9984072848436D-01)/X -
     +     3.1501078774085D-06)*E + SQRT(PIE4/X)
      F1 = (WV(NN,1)-E)/(X+X)
      IF (NROOTS-2) 420,440,460
  420 UV(NN,1) = F1/(WV(NN,1)-F1)
      RETURN
  440 Y = X-7.5D+00
      UV(NN,1) =(((((((((((((-1.43632730148572D-16*Y+
     +     2.38198922570405D-16)*
     +     Y+1.358319618800D-14)*Y-7.064522786879D-14)*Y-
     +     7.719300212748D-13)*Y+7.802544789997D-12)*Y+
     +     6.628721099436D-11)*Y-1.775564159743D-09)*Y+
     +     1.713828823990D-08)*Y-1.497500187053D-07)*Y+
     +     2.283485114279D-06)*Y-3.76953869614706D-05 )*Y+
     +     4.74791204651451D-04 )*Y-4.60448960876139D-03 )*Y+
     +     3.72458587837249D-02
      UV(NN,2) =(((((((((((( 2.48791622798900D-14*Y-
     +     1.36113510175724D-13)*Y-
     +     2.224334349799D-12)*Y+4.190559455515D-11)*Y-
     +     2.222722579924D-10)*Y-2.624183464275D-09)*Y+
     +     6.128153450169D-08)*Y-4.383376014528D-07)*Y-
     +     2.49952200232910D-06 )*Y+1.03236647888320D-04 )*Y-
     +     1.44614664924989D-03 )*Y+1.35094294917224D-02 )*Y-
     +     9.53478510453887D-02 )*Y+5.44765245686790D-01
      WV(NN,2) = ((F1-WV(NN,1))*UV(NN,1)+F1)*(1.0D+00+UV(NN,2))/
     +        (UV(NN,2)-UV(NN,1))
      WV(NN,1) = WV(NN,1)-WV(NN,2)
      RETURN
  460 F2 = (F1+F1+F1-E)/(X+X)
      Y = X-7.5D+00
      UV(NN,1) =((((((((((( 5.74429401360115D-16*Y+
     +     7.11884203790984D-16)*Y-
     +     6.736701449826D-14)*Y-6.264613873998D-13)*Y+
     +     1.315418927040D-11)*Y-4.23879635610964D-11 )*Y+
     +     1.39032379769474D-09 )*Y-4.65449552856856D-08 )*Y+
     +     7.34609900170759D-07 )*Y-1.08656008854077D-05 )*Y+
     +     1.77930381549953D-04 )*Y-2.39864911618015D-03 )*Y+
     +     2.39112249488821D-02
      UV(NN,2) =((((((((((( 1.13464096209120D-14*Y+
     +     6.99375313934242D-15)*Y-
     +     8.595618132088D-13)*Y-5.293620408757D-12)*Y-
     +     2.492175211635D-11)*Y+2.73681574882729D-09 )*Y-
     +     1.06656985608482D-08 )*Y-4.40252529648056D-07 )*Y+
     +     9.68100917793911D-06 )*Y-1.68211091755327D-04 )*Y+
     +     2.69443611274173D-03 )*Y-3.23845035189063D-02 )*Y+
     +     2.75969447451882D-01
      UV(NN,3) =(((((((((((( 6.66339416996191D-15*Y+
     +     1.84955640200794D-13)*Y-
     +     1.985141104444D-12)*Y-2.309293727603D-11)*Y+
     +     3.917984522103D-10)*Y+1.663165279876D-09)*Y-
     +     6.205591993923D-08)*Y+8.769581622041D-09)*Y+
     +     8.97224398620038D-06 )*Y-3.14232666170796D-05 )*Y-
     +     1.83917335649633D-03 )*Y+3.51246831672571D-02 )*Y-
     +     3.22335051270860D-01 )*Y+1.73582831755430D+00
        T1 = UV(NN,1)/(UV(NN,1)+1.0D+00)
        T2 = UV(NN,2)/(UV(NN,2)+1.0D+00)
        T3 = UV(NN,3)/(UV(NN,3)+1.0D+00)
        A2 = F2-T1*F1
        A1 = F1-T1*WV(NN,1)
        WV(NN,3) = (A2-T2*A1)/((T3-T2)*(T3-T1))
        WV(NN,2) = (T3*A1-A2)/((T3-T2)*(T2-T1))
        WV(NN,1) = WV(NN,1)-WV(NN,2)-WV(NN,3)
      RETURN
C
C     X = 10.0 TO 15.0                 NROOTS=1,2, OR 3
C
  480 WV(NN,1) = (((-1.8784686463512D-01/X+2.2991849164985D-01)/X -
     +     4.9893752514047D-01)/X-2.1916512131607D-05)*E + SQRT(PIE4/X)
      F1 = (WV(NN,1)-E)/(X+X)
      IF (NROOTS-2) 500,520,540
  500 UV(NN,1) = F1/(WV(NN,1)-F1)
      RETURN
  520 UV(NN,1) = ((((-1.01041157064226D-05*X+1.19483054115173D-03)*X -
     +     6.73760231824074D-02)*X+1.25705571069895D+00)*X + (((-
     +     8.57609422987199D+03/X+5.91005939591842D+03)/X -
     +     1.70807677109425D+03)/X+2.64536689959503D+02)/X -
     +     2.38570496490846D+01)*E + R12/(X-R12)
      UV(NN,2) = ((( 3.39024225137123D-04*X-9.34976436343509D-02)*X -
     +     4.22216483306320D+00)*X + (((-2.08457050986847D+03/X -
     +     1.04999071905664D+03)/X+3.39891508992661D+02)/X -
     +     1.56184800325063D+02)/X+8.00839033297501D+00)*E + R22/(X-R22)
      WV(NN,2) = ((F1-WV(NN,1))*UV(NN,1)+F1)*(1.0D+00+UV(NN,2))/
     +     (UV(NN,2)-UV(NN,1))
      WV(NN,1) = WV(NN,1)-WV(NN,2)
      RETURN
  540 F2 = (F1+F1+F1-E)/(X+X)
      Y = X-12.5D+00
      UV(NN,1) = ((((((((((( 4.42133001283090D-16*Y-
     +     2.77189767070441D-15)*Y-
     +     4.084026087887D-14)*Y+5.379885121517D-13)*Y+
     +     1.882093066702D-12)*Y-8.67286219861085D-11 )*Y+
     +     7.11372337079797D-10 )*Y-3.55578027040563D-09 )*Y+
     +     1.29454702851936D-07 )*Y-4.14222202791434D-06 )*Y+
     +     8.04427643593792D-05 )*Y-1.18587782909876D-03 )*Y+
     +     1.53435577063174D-02
      UV(NN,2) = ((((((((((( 6.85146742119357D-15*Y-
     +     1.08257654410279D-14)*Y-
     +     8.579165965128D-13)*Y+6.642452485783D-12)*Y+
     +     4.798806828724D-11)*Y-1.13413908163831D-09 )*Y+
     +     7.08558457182751D-09 )*Y-5.59678576054633D-08 )*Y+
     +     2.51020389884249D-06 )*Y-6.63678914608681D-05 )*Y+
     +     1.11888323089714D-03 )*Y-1.45361636398178D-02 )*Y+
     +     1.65077877454402D-01
      UV(NN,3) =(((((((((((( 3.20622388697743D-15*Y-
     +     2.73458804864628D-14)*Y-
     +     3.157134329361D-13)*Y+8.654129268056D-12)*Y-
     +     5.625235879301D-11)*Y-7.718080513708D-10)*Y+
     +     2.064664199164D-08)*Y-1.567725007761D-07)*Y-
     +     1.57938204115055D-06 )*Y+6.27436306915967D-05 )*Y-
     +     1.01308723606946D-03 )*Y+1.13901881430697D-02 )*Y-
     +     1.01449652899450D-01 )*Y+7.77203937334739D-01
        T1 = UV(NN,1)/(UV(NN,1)+1.0D+00)
        T2 = UV(NN,2)/(UV(NN,2)+1.0D+00)
        T3 = UV(NN,3)/(UV(NN,3)+1.0D+00)
        A2 = F2-T1*F1
        A1 = F1-T1*WV(NN,1)
        WV(NN,3) = (A2-T2*A1)/((T3-T2)*(T3-T1))
        WV(NN,2) = (T3*A1-A2)/((T3-T2)*(T2-T1))
        WV(NN,1) = WV(NN,1)-WV(NN,2)-WV(NN,3)
      RETURN
  560 IF (X .GT. 33.0D+00) GO TO 660
C
C     X = 15.0 TO 33.0                 NROOTS=1,2, OR 3
C
      E = EXP(-X)
      WV(NN,1) = (( 1.9623264149430D-01/X-4.9695241464490D-01)/X -
     +     6.0156581186481D-05)*E + SQRT(PIE4/X)
      F1 = (WV(NN,1)-E)/(X+X)
      IF (NROOTS-2) 580,600,620
  580 UV(NN,1) = F1/(WV(NN,1)-F1)
      RETURN
  600 UV(NN,1) = ((((-1.14906395546354D-06*X+1.76003409708332D-04)*X -
     +     1.71984023644904D-02)*X-1.37292644149838D-01)*X + (-
     +     4.75742064274859D+01/X+9.21005186542857D+00)/X -
     +     2.31080873898939D-02)*E + R12/(X-R12)
      UV(NN,2) = ((( 3.64921633404158D-04*X-9.71850973831558D-02)*X -
     +     4.02886174850252D+00)*X + (-1.35831002139173D+02/X -
     +     8.66891724287962D+01)/X+2.98011277766958D+00)*E + R22/(X-R22)
      WV(NN,2) = ((F1-WV(NN,1))*UV(NN,1)+F1)*(1.0D+00+UV(NN,2))/
     +     (UV(NN,2)-UV(NN,1))
      WV(NN,1) = WV(NN,1)-WV(NN,2)
      RETURN
  620 F2 = (F1+F1+F1-E)/(X+X)
      IF (X .GT. 20.0D+00) GO TO 640
      UV(NN,1) = ((((((-2.43270989903742D-06*X+3.57901398988359D-04)*X -
     +     2.34112415981143D-02)*X+7.81425144913975D-01)*X -
     +     1.73209218219175D+01)*X+2.43517435690398D+02)*X + (-
     +     1.97611541576986D+04/X+9.82441363463929D+03)/X -
     +     2.07970687843258D+03)*E + R13/(X-R13)
      UV(NN,2) = (((((-2.62627010965435D-04*X+3.49187925428138D-02)*X -
     +     3.09337618731880D+00)*X+1.07037141010778D+02)*X -
     +     2.36659637247087D+03)*X + ((-2.91669113681020D+06/X +
     +     1.41129505262758D+06)/X-2.91532335433779D+05)/X +
     +     3.35202872835409D+04)*E + R23/(X-R23)
      UV(NN,3) = ((((( 9.31856404738601D-05*X-2.87029400759565D-02)*X -
     +     7.83503697918455D-01)*X-1.84338896480695D+01)*X +
     +     4.04996712650414D+02)*X + (-1.89829509315154D+05/X +
     +     5.11498390849158D+04)/X-6.88145821789955D+03)*E + R33/(X-R33)
        T1 = UV(NN,1)/(UV(NN,1)+1.0D+00)
        T2 = UV(NN,2)/(UV(NN,2)+1.0D+00)
        T3 = UV(NN,3)/(UV(NN,3)+1.0D+00)
        A2 = F2-T1*F1
        A1 = F1-T1*WV(NN,1)
        WV(NN,3) = (A2-T2*A1)/((T3-T2)*(T3-T1))
        WV(NN,2) = (T3*A1-A2)/((T3-T2)*(T2-T1))
        WV(NN,1) = WV(NN,1)-WV(NN,2)-WV(NN,3)
      RETURN
  640 UV(NN,1) = ((((-4.97561537069643D-04*X-5.00929599665316D-02)*X +
     +     1.31099142238996D+00)*X-1.88336409225481D+01)*X -
     +     6.60344754467191D+02 /X+1.64931462413877D+02)*E + R13/(X-R13)
      UV(NN,2) = ((((-4.48218898474906D-03*X-5.17373211334924D-01)*X +
     +     1.13691058739678D+01)*X-1.65426392885291D+02)*X -
     +     6.30909125686731D+03 /X+1.52231757709236D+03)*E + R23/(X-R23)
      UV(NN,3) = ((((-1.38368602394293D-02*X-1.77293428863008D+00)*X +
     +     1.73639054044562D+01)*X-3.57615122086961D+02)*X -
     +     1.45734701095912D+04 /X+2.69831813951849D+03)*E + R33/(X-R33)
        T1 = UV(NN,1)/(UV(NN,1)+1.0D+00)
        T2 = UV(NN,2)/(UV(NN,2)+1.0D+00)
        T3 = UV(NN,3)/(UV(NN,3)+1.0D+00)
        A2 = F2-T1*F1
        A1 = F1-T1*WV(NN,1)
        WV(NN,3) = (A2-T2*A1)/((T3-T2)*(T3-T1))
        WV(NN,2) = (T3*A1-A2)/((T3-T2)*(T2-T1))
        WV(NN,1) = WV(NN,1)-WV(NN,2)-WV(NN,3)
      RETURN
C
C     X = 33.0 TO INFINITY             NROOTS=1,2, OR 3
C
  660 WV(NN,1) = SQRT(PIE4/X)
      IF (NROOTS-2) 680,700,720
  680 UV(NN,1) = 0.5D+00/(X-0.5D+00)
      RETURN
  700 IF (X .GT. 40.0D+00) GO TO 740
      E = EXP(-X)
      UV(NN,1) = (-8.78947307498880D-01*X+1.09243702330261D+01)*E +
     +     R12/(X-R12)
      UV(NN,2) = (-9.28903924275977D+00*X+8.10642367843811D+01)*E +
     +     R22/(X-R22)
      WV(NN,2) = ( 4.46857389308400D+00*X-7.79250653461045D+01)*E +
     +     W22*WV(NN,1)
      WV(NN,1) = WV(NN,1)-WV(NN,2)
      RETURN
  720 IF (X .GT. 47.0D+00) GO TO 760
      E = EXP(-X)
      UV(NN,1) = ((-7.39058467995275D+00*X+3.21318352526305D+02)*X -
     +     3.99433696473658D+03)*E + R13/(X-R13)
      UV(NN,2) = ((-7.38726243906513D+01*X+3.13569966333873D+03)*X -
     +     3.86862867311321D+04)*E + R23/(X-R23)
      UV(NN,3) = ((-2.63750565461336D+02*X+1.04412168692352D+04)*X -
     +     1.28094577915394D+05)*E + R33/(X-R33)
      WV(NN,3) = ((( 1.52258947224714D-01*X-8.30661900042651D+00)*X +
     +  1.92977367967984D+02)*X-1.67787926005344D+03)*E + W33*WV(NN,1)
      WV(NN,2) = (( 6.15072615497811D+01*X-2.91980647450269D+03)*X +
     +  3.80794303087338D+04)*E + W23*WV(NN,1)
      WV(NN,1) = WV(NN,1)-WV(NN,2)-WV(NN,3)
      RETURN
  740 UV(NN,1) = R12/(X-R12)
      UV(NN,2) = R22/(X-R22)
      WV(NN,2) = W22*WV(NN,1)
      WV(NN,1) = WV(NN,1)-WV(NN,2)
      RETURN
  760 UV(NN,1) = R13/(X-R13)
      UV(NN,2) = R23/(X-R23)
      UV(NN,3) = R33/(X-R33)
      WV(NN,2) = W23*WV(NN,1)
      WV(NN,3) = W33*WV(NN,1)
      WV(NN,1) = WV(NN,1)-WV(NN,2)-WV(NN,3)
      RETURN
      END
