      SUBROUTINE PCG2AL(ISUM,LENX,LCV,LCSS,LCP,LCCD,LCHCHG,LCLHCH,
     1       LCRCHG,LCLRCH,MXITER,ITER1,NCOL,NROW,NLAY,IN,IOUT,NPCOND)
C
C-----VERSION 0002 01MAY1989 PCG2AL
C
C     ******************************************************************
C     ALLOCATE STORAGE IN THE X ARRAY FOR PCG ARRAYS
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
C     ------------------------------------------------------------------
C
C-------PRINT A MESSAGE IDENTIFYING PCG PACKAGE
      WRITE(IOUT,1)
    1 FORMAT(1H0,'PCG2 -- CONJUGATE GRADIENT SOLUTION PACKAGE'
     1,', VERSION 2, 5/1/88')
C
C-------READ AND PRINT MXITER,ITER1 AND NPCOND
      READ(IN,2) MXITER,ITER1,NPCOND
    2 FORMAT(3I10)
      WRITE(IOUT,3) MXITER,ITER1,NPCOND
    3 FORMAT(' MAXIMUM OF',I4,' CALLS OF SOLUTION ROUTINE'/
     1       ,' MAXIMUM OF',I4,' INTERNAL ITERATIONS PER '
     2       ,'CALL TO SOLUTION ROUTINE'/
     3       ,' MATRIX PRECONDITIONING TYPE :',I5)
C
C-------ALLOCATE SPACE FOR THE PCG ARRAYS
      ISOLD=ISUM
      NRC=NROW*NCOL
      ISIZ=NRC*NLAY
      LCV=ISUM
      ISUM=ISUM+ISIZ
      LCSS=ISUM
      ISUM=ISUM+ISIZ
      LCP=ISUM
      ISUM=ISUM+ISIZ
      LCCD=ISUM
      IF(NPCOND.NE.2) ISUM=ISUM+ISIZ
      LCHCHG=ISUM
      ISUM=ISUM+MXITER*ITER1
      LCLHCH=ISUM
      ISUM=ISUM+3*MXITER*ITER1
      LCRCHG=ISUM
      ISUM=ISUM+MXITER*ITER1
      LCLRCH=ISUM
      ISUM=ISUM+3*MXITER*ITER1
C
C-------CALCULATE AND PRINT THE SPACE USED IN THE X ARRAY
      ICG=ISUM-ISOLD
      WRITE(IOUT,4) ICG
    4 FORMAT(1X,I7,' ELEMENTS IN X ARRAY ARE USED BY PCG')
      ISUM1=ISUM-1
      WRITE(IOUT,5) ISUM1,LENX
    5 FORMAT(1X,I7,' ELEMENTS OF X ARRAY USED OUT OF',I7)
      IF(ISUM1.GT.LENX) WRITE(IOUT,6)
    6 FORMAT(1X,'   ***X ARRAY MUST BE DIMENSIONED LARGER***')
C
C-------RETURN
      RETURN
      END
      SUBROUTINE PCG2RP(MXITER,ITER1,HCLOSE,RCLOSE,NPCOND,NBPOL,
     1                RELAX,IPRPCG,IN,IOUT,MUTPCG,NITER)
C
C-----VERSION 0002 01MAY1989 PCG2RP
C     01SEPT1990 IPCGCD OMITTED; NITER ADDED
C
C     ******************************************************************
C     READ DATA FOR PCG
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
C     ------------------------------------------------------------------
C
C-------READ HCLOSE,RCLOSE,RELAX,NBPOL,IPRPCG,MUTPCG
      READ(IN,1) HCLOSE,RCLOSE,RELAX,NBPOL,IPRPCG,MUTPCG
    1 FORMAT(3F10.0,4I10)
C
C-------PRINT MXITER,ITER1,NPCOND,HCLOSE,RCLOSE,RELAX,NBPOL,IPRPCG,
C-------MUTPCG
      WRITE(IOUT,100)
  100 FORMAT(1H0,///57X,'SOLUTION BY THE CONJUGATE-GRADIENT METHOD'
     1/57X,43('-'))
      WRITE(IOUT,115) MXITER
  115 FORMAT(1H0,38X,'MAXIMUM NUMBER OF CALLS TO PCG ROUTINE =',I9)
      WRITE(IOUT,120) ITER1
  120 FORMAT(1H ,42X,'MAXIMUM ITERATIONS PER CALL TO PCG =',I9)
      WRITE(IOUT,122) NPCOND
  122 FORMAT(1H ,49X,'MATRIX PRECONDITIONING TYPE =',I9)
      IF(NPCOND.EQ.2) WRITE(IOUT,123)
  123 FORMAT(1H ,58X,'THE MATRIX WILL BE SCALED')
      WRITE(IOUT,124) RELAX,NBPOL
  124 FORMAT(1H ,26X,'RELAXATION FACTOR (ONLY USED WITH',
     1' PRECOND. TYPE 1) =',E15.5,/,
     2 1H ,19X,'PARAMETER OF POLYMOMIAL PRECOND.'
     3 ,' = 2 (2) OR IS CALCULATED :',I9)
      WRITE(IOUT,125) HCLOSE
  125 FORMAT(1H ,43X,'HEAD CHANGE CRITERION FOR CLOSURE =',E15.5)
      WRITE(IOUT,127) RCLOSE
  127 FORMAT(1H ,39X,'RESIDUAL CHANGE CRITERION FOR CLOSURE =',E15.5)
      IF(IPRPCG.LE.0)IPRPCG=999
      WRITE(IOUT,130) IPRPCG,MUTPCG
  130 FORMAT(1H ,30X,'PCG HEAD AND RESIDUAL CHANGE PRINTOUT INTERVAL ='
     1,I9,/,1H ,23X,'PRINTING FROM SOLVER IS LIMITED(1) OR ',
     2'SUPPRESSED (>1) =',I9)
      NITER=0
C
      RETURN
      END
      SUBROUTINE PCG2AP(HNEW,IBOUND,CR,CC,CV,HCOF,RHS,V,SS,P,CD,
     1      HCHG,LHCH,RCHG,LRCH,KITER,NITER,HCLOSE,RCLOSE,ICNVG,
     2      KSTP,KPER,IPRPCG,MXITER,ITER1,NPCOND,NBPOL,NSTP,NCOL,NROW,

cigwmc Version 3.0  Variables IU,IP,SN,SP, and SR added per Mary Hill
cigwmc instructions

     3      NLAY,NODES,RELAX,IOUT,MUTPCG,iu,ip,sn,sp,sr)

cigwmc End of change

C-----VERSION 0002 01MAY1989 PCG2AP
C     01JULY1990 COMMENT STATEMENTS ADDED AND MODIFIED
C     01SEPT1990 IPCGCD OMITTED; STATEMENT 590 ADDED
C     27SEPT1990 STATEMENT IN DO 155 LOOP CHANGED
cigwmc Version 3.0 Lower case COMMENTS added by IGWMC per
cigwmc Mary Hill instructions.
c     01sept1992 Added statements related to sensitivity calculations
c                at the end of the 115 loop.  Changed the 510 FORMAT
c                statement and the preceding IF statement.

cigwmc End of change
C
C     ******************************************************************
C     SOLUTION BY THE CONJUGATE GRADIENT METHOD -
C                                          UP TO ITER1 ITERATIONS
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      PARAMETER (DZERO=0.D0,DONE=1.D0)
      DOUBLE PRECISION HNEW,HHCOF,RRHS
      DOUBLE PRECISION Z,B,D,E,F,H,S,ALPHA
      DOUBLE PRECISION ZHNEW,BHNEW,DHNEW,FHNEW,HHNEW,SHNEW
      DOUBLE PRECISION SRNEW,SROLD,SSCR,SSCC,SSCV,VCC,VCR,VCV
      DOUBLE PRECISION CDCC,CDCR,CDCV
      DOUBLE PRECISION PN,VN,HCHGN,RCHGN,PAP
      DOUBLE PRECISION FCC,FCR,FCV,FV
C
      DIMENSION HNEW(NODES), IBOUND(NODES), CR(NODES), CC(NODES),
     1  CV(NODES), HCOF(NODES), RHS(NODES),
     2  V(NODES), SS(NODES), P(NODES), CD(NODES), HCHG(MXITER*ITER1),
     3  LHCH(3,MXITER*ITER1), RCHG(MXITER*ITER1), LRCH(3,MXITER*ITER1)
C     ------------------------------------------------------------------
C
C-------ASSIGN VARIABLE EQUAL TO THE NUMBER OF CELLS IN ONE LAYER
      NRC=NROW*NCOL
C-------INITIALIZE VARIABLES USED TO CALCULATE ITERATION PARAMETERS
      SRNEW=DZERO
      BPOLY=0.
      IF(NPCOND.NE.1) RELAX=1.
      NORM=0
      IF(NPCOND.EQ.2) NORM=1
C-------INITIALIZE VARIABLE USED TO TEST FOR NEGATIVE CHOLESKY DIAGONAL
      CD1=0.
C------CLEAR PCG WORK ARRAYS.
      DO 100 N=1,NODES
      SS(N)=0.
      P(N)=0.
  100 V(N)=0.
C------FOR NPCOND=1, INITIALIZE CHOLESKY DIAGONAL
      IF(NPCOND.EQ.1) THEN
        DO 105 N=1,NODES
  105   CD(N)=0.
      ENDIF
C
C------CALCULATE THE RESIDUAL. IF NORM=1, CALCULATE THE DIAGONALS OF
C------THE A MATRIX,AND STORE THEM IN HCOF.
      DO 115 K=1,NLAY
      DO 115 I=1,NROW
      DO 115 J=1,NCOL
C
C-------CALCULATE 1 DIMENSIONAL SUBSCRIPT OF CURRENT CELL AND
C-------SKIP CALCULATIONS IF CELL IS INACTIVE
      N=J+(I-1)*NCOL+(K-1)*NRC
      IF(IBOUND(N).EQ.0) THEN
        CC(N)=0.
        CR(N)=0.
        IF(N.LE.(NODES-NRC)) CV(N)=0.
        IF(N.GE.2) CR(N-1)=0.
        IF(N.GE.NCOL+1) CC(N-NCOL)=0.
        IF(N.LE.(NODES-NRC).AND.N.GE.NRC+1) CV(N-NRC)=0.
        HCOF(N)=0.
        RHS(N)=0.
        GO TO 115
      ENDIF
C
C-------CALCULATE 1 DIMENSIONAL SUBSCRIPTS FOR LOCATING THE 6
C-------SURROUNDING CELLS
      NRN=N+NCOL
      NRL=N-NCOL
      NCN=N+1
      NCL=N-1
      NLN=N+NRC
      NLL=N-NRC
C
C-------CALCULATE 1 DIMENSIONAL SUBSCRIPTS FOR CONDUCTANCE TO THE 6
C-------SURROUNDING CELLS.
      NCF=N
      NCD=N-1
      NRB=N-NCOL
      NRH=N
      NLS=N
      NLZ=N-NRC
C
C-----GET CONDUCTANCES TO NEIGHBORING CELLS
C-------NEIGHBOR IS 1 ROW BACK
      B=DZERO
      BHNEW=DZERO
      IF(I.NE.1) THEN
        B=CC(NRB)
        BHNEW=B*(HNEW(NRL)-HNEW(N))
      ENDIF
C
C-------NEIGHBOR IS 1 ROW AHEAD
      H=DZERO
      HHNEW=DZERO
      IF(I.NE.NROW) THEN
        H=CC(NRH)
        HHNEW=H*(HNEW(NRN)-HNEW(N))
      ENDIF
C
C-------NEIGHBOR IS 1 COLUMN BACK
      D=DZERO
      DHNEW=DZERO
      IF(J.NE.1) THEN
        D=CR(NCD)
        DHNEW=D*(HNEW(NCL)-HNEW(N))
      ENDIF
C
C-------NEIGHBOR IS 1 COLUMN AHEAD
      F=DZERO
      FHNEW=DZERO
      IF(J.NE.NCOL) THEN
        F=CR(NCF)
        FHNEW=F*(HNEW(NCN)-HNEW(N))
      ENDIF
C
C-------NEIGHBOR IS 1 LAYER BEHIND
      Z=DZERO
      ZHNEW=DZERO
      IF(K.NE.1) THEN
        Z=CV(NLZ)
        ZHNEW=Z*(HNEW(NLL)-HNEW(N))
      ENDIF
C
C-------NEIGHBOR IS 1 LAYER AHEAD
      S=DZERO
      SHNEW=DZERO
      IF(K.NE.NLAY) THEN
        S=CV(NLS)
        SHNEW=S*(HNEW(NLN)-HNEW(N))
      ENDIF
C
      IF(I.EQ.NROW) CC(N)=0.
      IF(J.EQ.NCOL) CR(N)=0.
C
C-------CALCULATE THE RESIDUAL AND STORE IT IN RHS.  TO SCALE A,
C-------CALCULATE THE DIAGONAL OF THE A MATRIX, AND STORE IT IN HCOF.
      E=-Z-B-D-F-H-S
      RRHS=RHS(N)
      HHCOF=HNEW(N)*HCOF(N)
      RHS(N)=RRHS-ZHNEW-BHNEW-DHNEW-HHCOF-FHNEW-HHNEW-SHNEW
      IF(NORM.EQ.1) HCOF(N)=HCOF(N)+E
      IF(IBOUND(N).LT.0.) RHS(N)=0.
  115 CONTINUE
C
C-------SCALE CC,CR,CV,RHS AND HNEW IF NORM=1.
      IF(NORM.EQ.1) THEN
        DO 120 K=1,NLAY
        DO 120 I=1,NROW
        DO 120 J=1,NCOL
          N=J+(I-1)*NCOL+(K-1)*NRC
          IF(IBOUND(N).EQ.0) GO TO 120
          HHCOF=SQRT(-HCOF(N))
          IF(N.LE.(NODES-NCOL).AND.CC(N).GT.0.)
     #      CC(N)=CC(N)/(HHCOF*(SQRT(-HCOF(N+NCOL))))
          IF(CR(N).GT.0.) CR(N)=CR(N)/(HHCOF*(SQRT(-HCOF(N+1))))
          IF(N.LE.(NODES-NRC).AND.CV(N).GT.0.)
     #      CV(N)=CV(N)/(HHCOF*(SQRT(-HCOF(N+NRC))))
          HNEW(N)=HNEW(N)*HHCOF
          RHS(N)=RHS(N)/HHCOF
  120   CONTINUE
      ENDIF
C
C-------CALCULATE PARAMETER B OF THE POLYNOMIAL PRECONDITIONING METHOD
      IF(NPCOND.NE.2) GO TO 152
      IF(NBPOL.EQ.2) THEN
        BPOLY=2
        GO TO 151
      ENDIF
      DO 150 K=1,NLAY
      DO 150 I=1,NROW
      DO 150 J=1,NCOL
C
      N=J+(I-1)*NCOL+(K-1)*NRC
      IF(IBOUND(N).LE.0)GO TO 150
C
      NCF=N
      NCD=N-1
      NRB=N-NCOL
      NRH=N
      NLS=N
      NLZ=N-NRC
C
      B=DZERO
      IF(I.NE.1) B=CC(NRB)
      H=DZERO
      IF(I.NE.NROW) H=CC(NRH)
      D=DZERO
      IF(J.NE.1) D=CR(NCD)
      F=DZERO
      IF(J.NE.NCOL) F=CR(NCF)
      Z=DZERO
      IF(K.NE.1) Z=CV(NLZ)
      S=DZERO
      IF(K.NE.NLAY) S=CV(NLS)
C
C-------NOTE : ABS. VAL. OF THE DIAG. OF THE SCALED A MATRIX IS 1.
      HHCOF=HCOF(N)
      IF(NORM.EQ.1) HHCOF=DONE
      T=DABS(Z)+DABS(B)+DABS(D)+ABS(HHCOF)+DABS(F)+DABS(H)+DABS(S)
      IF(T.GT.BPOLY) BPOLY=T
  150 CONTINUE
  151 CONTINUE
C
C-------CALCULATE ITERATION PARAMETERS FOR POLYNOMIAL PRECONDITIONING
C-------METHOD FOR A NEGATIVE DEFINITE MATRIX.
        C0=(15./32.)*(BPOLY**3)
        C1=(27./16.)*(BPOLY**2)
        C2=(9./4.)*BPOLY
  152 CONTINUE
C
C-------START INTERNAL ITERATIONS
      IITER=0
      IF(KITER.EQ.1) NITER=0
      ICNVG=0
      IICNVG=0
153   CONTINUE
      IITER=IITER+1
      NITER=NITER+1
C
C-------INITIALIZE VARIABLES THAT TRACK MAXIMUM HEAD CHANGE AND RESIDUAL
C-------VALUE DURING EACH ITERATIONS
      BIGH=0.
      BIGR=0.
C
C
C-------CHECK NPCOND FOR PRECONDITIONING TYPE AND EXECUTE PROPER CODE
      IF(NPCOND.EQ.2) GO TO 165
C
C-------CHOLESKY PRECONDITIONING
C
C-------STEP THROUGH CELLS TO CALCULATE THE DIAGONAL OF THE CHOLESKY
C-------MATRIX (FIRST INTERNAL ITERATION ONLY) AND THE INTERMEDIATE
C-------SOLUTION.  STORE THEM IN CD AND V, RESPECTIVELY.
      DO 155 K=1,NLAY
      DO 155 I=1,NROW
      DO 155 J=1,NCOL
C
      N=J+(I-1)*NCOL+(K-1)*NRC
      IF(IBOUND(N).LE.0)GO TO 155
C
C-------CALCULATE V
      H=DZERO
      VCC=DZERO
      IC=N-NCOL
      IF(I.NE.1) THEN
        H=CC(IC)
        IF(CD(IC).NE.0.) VCC=H*V(IC)/CD(IC)
      ENDIF
C
      F=DZERO
      VCR=DZERO
      IR=N-1
      IF(J.NE.1) THEN
        F=CR(IR)
        IF(CD(IR).NE.0.) VCR=F*V(IR)/CD(IR)
      ENDIF
C
      S=DZERO
      VCV=DZERO
      IL=N-NRC
      IF(K.NE.1) THEN
        S=CV(IL)
        IF(CD(IL).NE.0.) VCV=S*V(IL)/CD(IL)
      ENDIF
      V(N)=RHS(N)-VCR-VCC-VCV
C
C-------CALCULATE CD - FIRST INTERNAL ITERATION ONLY
      IF(IITER.EQ.1) THEN
          CDCR=DZERO
          CDCC=DZERO
          CDCV=DZERO
          FCC=DZERO
          FCR=DZERO
          FCV=DZERO
          IF(IR.GT.0.AND.CD(IR).NE.0.) CDCR=(F**2)/CD(IR)
          IF(IC.GT.0.AND.CD(IC).NE.0.) CDCC=(H**2)/CD(IC)
          IF(IL.GT.0.AND.CD(IL).NE.0.) CDCV=(S**2)/CD(IL)
          IF(NPCOND.EQ.1) THEN
              IF(IR.GT.0) THEN
                  FV=CV(IR)
cigwmc Version 3.0 comment statement modified per Mary Hill instructions.

C                 MODIFIED FROM HILL(1990) 9/27/90: 2 REPLACES 1 AND

cigwmc End of change

C                 2 REPLACES 1
                  IF(K.EQ.NLAY.AND.((J+I).GT.2)) FV=DZERO
                  IF(CD(IR).NE.0.) FCR=(F/CD(IR))*(CC(IR)+FV)
              ENDIF
              IF(IC.GT.0) THEN
                  FV=CV(IC)
                  IF(K.EQ.NLAY.AND.(I.GT.1)) FV=DZERO
                  IF(CD(IC).NE.0.) FCC=(H/CD(IC))*(CR(IC)+FV)
              ENDIF
              IF(IL.GT.0) THEN
                  IF(CD(IL).NE.0.) FCV=(S/CD(IL))*(CR(IL)+CC(IL))
              ENDIF
          ENDIF
          IF(NORM.EQ.0) THEN
             B=DZERO
             H=DZERO
             D=DZERO
             F=DZERO
             Z=DZERO
             S=DZERO
             IF(I.NE.1) B=CC(IC)
             IF(I.NE.NROW) H=CC(N)
             IF(J.NE.1) D=CR(IR)
             IF(J.NE.NCOL) F=CR(N)
             IF(K.NE.1) Z=CV(IL)
             IF(K.NE.NLAY) S=CV(N)
             HHCOF=HCOF(N)-Z-B-D-F-H-S
          ENDIF
          IF(NORM.EQ.1) HHCOF=-DONE
          CD(N)=HHCOF-CDCR-CDCC-CDCV-RELAX*(FCR+FCC+FCV)
         IF(CD1.EQ.0..AND.CD(N).NE.0.) CD1=CD(N)

cigwmc Version 3.0 LT changed to LE per Mary Hill instructions.
cigwmc   IF(CD(N)*CD1.LT.0.) THEN

c--------.LT. changed to .LE. 01SEPT1991
         IF(CD(N)*CD1.Le.0.) THEN

cigwmc End of change

           WRITE(IOUT,510)

cigwmc Version 3.0 510 comments and format statement changed per
cigmwc Mary Hill instructions.

c--------Changed 510 FORMAT 01SEPT1991 and 20MAR1992

  510      FORMAT(//,'  CHOLESKY DIAGONAL LESS THAN OR EQUAL TO ZERO',
     1        ' -- EXECUTION TERMINATED (MATRIX NOT DIAGONALLY DOMINAN',
     2        'T -- IN SOME SITUATIONS SETTING RELAX < 1.0 RESOLVES ',
     3         'THE PROBLEM)')

cigwmc End of change

           STOP
         ENDIF
      ENDIF
C
  155 CONTINUE
C
C-------STEP THROUGH EACH CELL AND SOLVE FOR S OF THE CONJUGATE
C-------GRADIENT ALGORITHM BY BACK SUBSTITUTION. STORE RESULT IN SS.
      DO 160 KK=NLAY,1,-1
      DO 160 II=NROW,1,-1
      DO 160 JJ=NCOL,1,-1
C
      N=JJ+(II-1)*NCOL+(KK-1)*NRC
      IF(IBOUND(N).LE.0)GO TO 160
C
      NC=N+1
      NR=N+NCOL
      NL=N+NRC

cigwmc Version 3.0 Original comments in this location deleted per
cigwmc Mary Hill instructions.

cigwmc End of change

      SSCR=DZERO
      SSCC=DZERO
      SSCV=DZERO
      IF(JJ.NE.NCOL) SSCR=CR(N)*SS(NC)/CD(N)
      IF(II.NE.NROW) SSCC=CC(N)*SS(NR)/CD(N)
      IF(KK.NE.NLAY) SSCV=CV(N)*SS(NL)/CD(N)
      VN=V(N)/CD(N)
      SS(N)=VN-SSCR-SSCC-SSCV
  160 CONTINUE
C-------SKIP OVER OTHER PRECONDITIONING TYPES
      GO TO 199
  165 CONTINUE
C
C-------POLYNOMIAL PRECONDITIONING
      DO 170 N=1,NODES
      V(N)=RHS(N)
  170 CONTINUE
      CALL SPCG2E(IBOUND,RHS,HCOF,CR,CC,CV,V,SS,C2,NORM,NCOL,NROW,
     1            NLAY,NODES)
      CALL SPCG2E(IBOUND,RHS,HCOF,CR,CC,CV,SS,V,C1,NORM,NCOL,NROW,
     1            NLAY,NODES)
      CALL SPCG2E(IBOUND,RHS,HCOF,CR,CC,CV,V,SS,C0,NORM,NCOL,NROW,
     1            NLAY,NODES)
  199 CONTINUE
C
C-------CALCULATE P OF THE CONJUGATE GRADIENT ALGORITHM
      SROLD=SRNEW
      SRNEW=DZERO
      DO 200 N=1,NODES
        IF(IBOUND(N).LE.0)GO TO 200
        SRNEW=SRNEW+SS(N)*RHS(N)
  200 CONTINUE
C
      IF(IITER.EQ.1) THEN
         DO 205 N=1,NODES
  205    P(N)=SS(N)
      ELSE
         DO 210 N=1,NODES
  210    P(N)=SS(N)+(SRNEW/SROLD)*P(N)
      ENDIF
C
C-------CALCULATE ALPHA OF THE CONJUGATE GRADIENT ROUTINE.
C-------FOR THE DENOMINATOR OF ALPHA, MULTIPLY THE MATRIX A BY THE
C-------VECTOR P, AND STORE IN V; THEN MULTIPLY P BY V.  STORE IN PAP.
      PAP=DZERO
      DO 290 K=1,NLAY
      DO 290 I=1,NROW
      DO 290 J=1,NCOL
C
      N=J+(I-1)*NCOL+(K-1)*NRC
      V(N)=0.
      IF(IBOUND(N).LE.0)GO TO 290
C
      NRN=N+NCOL
      NRL=N-NCOL
      NCN=N+1
      NCL=N-1
      NLN=N+NRC
      NLL=N-NRC
C
      NCF=N
      NCD=NCL
      NRB=NRL
      NRH=N
      NLS=N
      NLZ=NLL
C
      B=DZERO
      IF(I.NE.1) B=CC(NRB)
      H=DZERO
      IF(I.NE.NROW)H=CC(NRH)
      D=DZERO
      IF(J.NE.1) D=CR(NCD)
      F=DZERO
      IF(J.NE.NCOL) F=CR(NCF)
      Z=DZERO
      IF(K.NE.1) Z=CV(NLZ)
      S=DZERO
      IF(K.NE.NLAY) S=CV(NLS)
C
      IF(NORM.EQ.0) PN=P(N)
      IF(NORM.EQ.1) PN=DZERO
      BHNEW=DZERO
      HHNEW=DZERO
      DHNEW=DZERO
      FHNEW=DZERO
      ZHNEW=DZERO
      SHNEW=DZERO
      IF(NRL.GT.0) BHNEW=B*(P(NRL)-PN)
      IF(NRN.LE.NODES) HHNEW=H*(P(NRN)-PN)
      IF(NCL.GT.0) DHNEW=D*(P(NCL)-PN)
      IF(NCN.LE.NODES) FHNEW=F*(P(NCN)-PN)
      IF(NLL.GT.0) ZHNEW=Z*(P(NLL)-PN)
      IF(NLN.LE.NODES) SHNEW=S*(P(NLN)-PN)
C
C-------CALCULATE THE PRODUCT OF MATRIX A AND VECTOR P AND STORE
C-------RESULT IN V.
      PN=HCOF(N)*P(N)
      IF(NORM.EQ.1) PN=-P(N)
      VN=ZHNEW+BHNEW+DHNEW+PN+FHNEW+HHNEW+SHNEW
      V(N)=VN
      PAP=PAP+P(N)*VN
  290 CONTINUE
C
C-------CALCULATE ALPHA
      ALPHA=SRNEW/PAP
C
C-------CALCULATE NEW HEADS AND RESIDUALS, AND SAVE THE LARGEST
C-------CHANGE IN HEAD AND THE LARGEST VALUE OF THE RESIDUAL.
      DO 300 K=1,NLAY
      DO 300 I=1,NROW
      DO 300 J=1,NCOL
C
      N=J+(I-1)*NCOL+(K-1)*NRC
      IF(IBOUND(N).LE.0) GO TO 300
C
C-------HEAD
      HCHGN=ALPHA*P(N)
      IF(DABS(HCHGN).GT.ABS(BIGH)) THEN
          BIGH=HCHGN
          IH=I
          JH=J
          KH=K
          NH=N
      ENDIF
      HNEW(N)=HNEW(N)+HCHGN
C
C--------RESIDUAL (V IS THE PRODUCT OF MATRIX A AND VECTOR P)
      RCHGN=-ALPHA*V(N)
      RHS(N)=RHS(N)+RCHGN
      IF(ABS(RHS(N)).GT.ABS(BIGR)) THEN
          BIGR=RHS(N)
          IR=I
          JR=J
          KR=K
          NR=N
      ENDIF
  300 CONTINUE
C
C-------UNSCALE LARGEST CHANGE IN HEAD AND LARGEST RESIDUAL, AND
C-------CHECK THE CONVERGENCE CRITERION
       IF(NORM.EQ.1) THEN
          BIGH=BIGH/SQRT(-HCOF(NH))
          BIGR=BIGR*SQRT(-HCOF(NR))
       ENDIF
       IF(MXITER.EQ.1) THEN
          IF(ABS(BIGH).LE.HCLOSE.AND.ABS(BIGR).LE.RCLOSE) ICNVG=1
       ELSE
          IF(IITER.EQ.1.AND.
     1       ABS(BIGH).LE.HCLOSE.AND.ABS(BIGR).LE.RCLOSE) ICNVG=1
       ENDIF
       IF(ABS(BIGH).LE.HCLOSE.AND.ABS(BIGR).LE.RCLOSE) IICNVG=1
C
C-------STORE THE LARGEST UNSCALED HEAD CHANGE AND RESIDUAL VALUE
C-------(THIS ITERATION) AND THEIR LOCATIONS.
      II=NITER
      HCHG(II)=BIGH
      LHCH(1,II)=KH
      LHCH(2,II)=IH
      LHCH(3,II)=JH
C
      RCHG(II)=BIGR
      LRCH(1,II)=KR
      LRCH(2,II)=IR
      LRCH(3,II)=JR
C
C-------GO TO NEXT INTERNAL ITERATION IF CONVERGENCE HAS NOT BEEN
C-------REACHED AND IITER IS LESS THAN ITER1
      IF(MXITER.EQ.1) THEN
        IF(ICNVG.EQ.0.AND.IITER.LT.ITER1) GO TO 153
      ELSE
        IF(IICNVG.EQ.0.AND.IITER.LT.ITER1) GO TO 153
      ENDIF
C
C-------UNSCALE CR,CC,CV AND HNEW
      IF(NORM.EQ.1) THEN
      DO 310 N=1,NODES
      IF(IBOUND(N).EQ.0) GO TO 310
      HHCOF=SQRT(-HCOF(N))
      IF(N.LE.(NODES-NCOL).AND.CC(N).GT.0.)
     #   CC(N)=CC(N)*(HHCOF*(SQRT(-HCOF(N+NCOL))))
      IF(N.LE.(NODES-1).AND.CR(N).GT.0.)
     #   CR(N)=CR(N)*(HHCOF*(SQRT(-HCOF(N+1))))
      IF(N.LE.(NODES-NRC).AND.CV(N).GT.0.)
     #   CV(N)=CV(N)*(HHCOF*(SQRT(-HCOF(N+NRC))))
      HNEW(N)=HNEW(N)/HHCOF
  310 CONTINUE
      ENDIF
C
C-------IF END OF TIME STEP, PRINT # OF ITERATIONS THIS STEP
      IF(ICNVG.EQ.0 .AND. KITER.NE.MXITER) GO TO 600
      IF(MUTPCG.GT.1) GO TO 590
      IF(KSTP.EQ.1) WRITE(IOUT,500)
  500 FORMAT(1H0)
      WRITE(IOUT,501) KITER,KSTP,KPER,NITER
  501 FORMAT(1X,I5,' CALLS TO PCG ROUTINE FOR TIME STEP',I4,
     1' IN STRESS PERIOD',I3,/1X,I5,' TOTAL ITERATIONS')
      IF(MUTPCG.EQ.1) GO TO 590
C
C-------PRINT HEAD CHANGE EACH ITERATION IF PRINTOUT INTERVAL IS REACHED
      IF(ICNVG.EQ.0 .OR. KSTP.EQ.NSTP .OR. MOD(KSTP,IPRPCG).EQ.0)

cigwmc Version 3.0 Call to SPCG2P modified per Mary Hill instructions

cigwmc1      CALL SPCG2P(HCHG,LHCH,RCHG,LRCH,IITER,KITER,KSTP,KPER,
     1      CALL SPCG2P(HCHG,LHCH,RCHG,LRCH,

cigwmc End of change

     2                  ITER1,NITER,MXITER,IOUT,NPCOND,BPOLY)
590   NITER=0
C
C-------RETURN
600   RETURN
C
      END
cigwmc Version 3.0 SUBROUTINE statement modified per Mary Hill instructions.

cigwmcSUBROUTINE SPCG2P(HCHG,LHCH,RCHG,LRCH,IITER,KITER,KSTP,KPER,ITER1,
      SUBROUTINE SPCG2P(HCHG,LHCH,RCHG,LRCH,ITER1,

cigmwc End of change

     1                  NITER,MXITER,IOUT,NPCOND,BPOLY)
C
C
C-----VERSION 0001 01MAY1988 SPCG2P
C     ******************************************************************
C     PRINT MAXIMUM HEAD CHANGE AND RESIDUAL VALUE FOR EACH ITERATION
C                           DURING A TIME STEP
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
C
      DIMENSION HCHG(MXITER*ITER1), LHCH(3,MXITER*ITER1)
      DIMENSION RCHG(MXITER*ITER1), LRCH(3,MXITER*ITER1)
C     ------------------------------------------------------------------
C
      IF(NPCOND.EQ.2) WRITE(IOUT,2) BPOLY
    2 FORMAT(1H0,'B OF THE POLYNOMIAL PRECONDITIONING METHOD: ',E12.4)
      WRITE(IOUT,5)
    5 FORMAT(1H0,'MAXIMUM HEAD CHANGE FOR EACH ITERATION:'/
     1    1H0,4('    HEAD CHANGE  LAYER,ROW,COL')/1X,120('-'))
      WRITE (IOUT,10) (HCHG(J),(LHCH(I,J),I=1,3),J=1,NITER)
      WRITE(IOUT,11)
      WRITE(IOUT,15)
   15 FORMAT(1H0,'MAXIMUM RESIDUAL FOR EACH ITERATION:'/
     1    1H0,4('     RESIDUAL    LAYER,ROW,COL')/1X,120('-'))
      WRITE (IOUT,10) (RCHG(J),(LRCH(I,J),I=1,3),J=1,NITER)
      WRITE(IOUT,11)
C
      RETURN
C
   10 FORMAT((1X,4(4X,G12.4,' (',I3,',',I3,',',I3,')')))
   11 FORMAT(1H0)
C
      END
      SUBROUTINE SPCG2E(IBOUND,RHS,HCOF,CR,CC,CV,VIN,VOUT,C,NORM,NCOL,
     1                  NROW,NLAY,NODES)
C
C
C-----VERSION 0001 01MAY1989 SPCG2E
C     ******************************************************************
C           MATRIX MULTIPLICATIONS FOR POLYNOMIAL PRECONDITIONING
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
C
      DOUBLE PRECISION VN,CRHS,Z,B,D,F,H,S,ZV,BV,DV,FV,HV,SV,DZERO
      DIMENSION IBOUND(NODES),CR(NODES),CC(NODES),CV(NODES),
     1  RHS(NODES),VIN(NODES),VOUT(NODES),HCOF(NODES)
C     ------------------------------------------------------------------
C
      DZERO=0.
      NRC=NROW*NCOL
      DO 290 K=1,NLAY
      DO 290 I=1,NROW
      DO 290 J=1,NCOL
C
      N=J+(I-1)*NCOL+(K-1)*NRC
      VOUT(N)=0.
      IF(IBOUND(N).LE.0)GO TO 290
C
      NRN=N+NCOL
      NRL=N-NCOL
      NCN=N+1
      NCL=N-1
      NLN=N+NRC
      NLL=N-NRC
C
      NCF=N
      NCD=NCL
      NRB=NRL
      NRH=N
      NLS=N
      NLZ=NLL
C
      B=DZERO
      BV=DZERO
      IF(I.NE.1.AND.IBOUND(NRL).GE.0) THEN
        B=CC(NRB)
        BV=B*VIN(NRL)
      ENDIF
      H=DZERO
      HV=DZERO
      IF(I.NE.NROW.AND.IBOUND(NRN).GE.0) THEN
        H=CC(NRH)
        HV=H*VIN(NRN)
      ENDIF
      D=DZERO
      DV=DZERO
      IF(J.NE.1.AND.IBOUND(NCL).GE.0) THEN
        D=CR(NCD)
        DV=D*VIN(NCL)
      ENDIF
      F=DZERO
      FV=DZERO
      IF(J.NE.NCOL.AND.IBOUND(NCN).GE.0) THEN
        F=CR(NCF)
        FV=F*VIN(NCN)
      ENDIF
      Z=DZERO
      ZV=DZERO
      IF(K.NE.1.AND.IBOUND(NLL).GE.0) THEN
        Z=CV(NLZ)
        ZV=Z*VIN(NLL)
      ENDIF
      S=DZERO
      SV=DZERO
      IF(K.NE.NLAY.AND.IBOUND(NLN).GE.0) THEN
        S=CV(NLS)
        SV=S*VIN(NLN)
      ENDIF
C
C-------CALCULATE THE PRODUCT OF MATRIX A AND VECTOR VIN AND STORE
C------ RESULT IN VOUT
      VN=HCOF(N)*VIN(N)
      IF(NORM.EQ.1) VN=-VIN(N)
      CRHS=C*RHS(N)
      VOUT(N)=CRHS+ZV+BV+DV+VN+FV+HV+SV
  290 CONTINUE
      RETURN
      END
