C FILE SIP1.F77                                                                 
      SUBROUTINE SIP1AL(ISUM,LENX,LCEL,LCFL,LCGL,LCV,LCHDCG,LCLRCH,             
     1       LCW,MXITER,NPARM,NCOL,NROW,NLAY,IN,IOUT)                           
C                                                                               
C-----VERSION 1745 24APR1985 SIP1AL                                             
C                                                                               
C     ******************************************************************        
C     ALLOCATE STORAGE IN THE X ARRAY FOR SIP ARRAYS                            
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
C     ------------------------------------------------------------------        
C                                                                               
C1------PRINT A MESSAGE IDENTIFYING SIP PACKAGE                                 
      WRITE(IOUT,1)IN                                                           
    1 FORMAT(1H0,'SIP1 -- STRONGLY IMPLICIT PROCEDURE SOLUTION PACKAGE'         
     1,', VERSION 1, 04/24/85',' INPUT READ FROM UNIT',I3)                      
C                                                                               
C2------READ AND PRINT MXITER AND NPARM                                         
      READ(IN,2) MXITER,NPARM                                                   
    2 FORMAT(2I10)                                                              
      WRITE(IOUT,3) MXITER,NPARM                                                
    3 FORMAT(1X,'MAXIMUM OF',I4,' ITERATIONS ALLOWED FOR CLOSURE'/              
     1       1X,I2,' ITERATION PARAMETERS')                                     
C                                                                               
C3------ALLOCATE SPACE FOR THE SIP ARRAYS                                       
      ISOLD=ISUM                                                                
      NRC=NROW*NCOL                                                             
      ISIZ=NRC*NLAY                                                             
      LCEL=ISUM                                                                 
      ISUM=ISUM+ISIZ                                                            
      LCFL=ISUM                                                                 
      ISUM=ISUM+ISIZ                                                            
      LCGL=ISUM                                                                 
      ISUM=ISUM+ISIZ                                                            
      LCV=ISUM                                                                  
      ISUM=ISUM+ISIZ                                                            
      LCHDCG=ISUM                                                               
      ISUM=ISUM+MXITER                                                          
      LCLRCH=ISUM                                                               
      ISUM=ISUM+3*MXITER                                                        
      LCW=ISUM                                                                  
      ISUM=ISUM+NPARM                                                           
C                                                                               
C4------CALCULATE AND PRINT THE SPACE USED IN THE X ARRAY                       
      ISP=ISUM-ISOLD                                                            
      WRITE(IOUT,4) ISP                                                         
    4 FORMAT(1X,I6,' ELEMENTS IN X ARRAY ARE USED BY SIP')                      
      ISUM1=ISUM-1                                                              
      WRITE(IOUT,5) ISUM1,LENX                                                  
    5 FORMAT(1X,I6,' 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                                                                               
C5------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE SIP1RP(NPARM,MXITER,ACCL,HCLOSE,W,IN,IPCALC,IPRSIP,            
     1                IOUT)                                                     
C                                                                               
C-----VERSION 0925 16DEC1982 SIP1RP                                             
C                                                                               
C     ******************************************************************        
C     READ DATA FOR SIP                                                         
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
      DIMENSION W(NPARM)                                                        
C     ------------------------------------------------------------------        
C                                                                               
C1------READ ACCL,HCLOSE,WSEED,IPCALC,IPRSIP                                    
      READ(IN,1) ACCL,HCLOSE,IPCALC,WSEED,IPRSIP                                
    1 FORMAT(2F10.0,I10,F10.0,I10)                                              
      IF(ACCL.EQ.0.) ACCL=1.                                                    
C                                                                               
C2------PRINT DATA VALUES JUST READ                                             
      WRITE(IOUT,100)                                                           
  100 FORMAT(1H0,///57X,'SOLUTION BY THE STRONGLY IMPLICIT PROCEDURE'           
     1/57X,43('-'))                                                             
      WRITE(IOUT,115) MXITER                                                    
  115 FORMAT(1H0,47X,'MAXIMUM ITERATIONS ALLOWED FOR CLOSURE =',I9)             
      WRITE(IOUT,120) ACCL                                                      
  120 FORMAT(1H ,63X,'ACCELERATION PARAMETER =',G15.5)                          
      WRITE(IOUT,125) HCLOSE                                                    
  125 FORMAT(1H ,52X,'HEAD CHANGE CRITERION FOR CLOSURE =',E15.5)               
      IF(IPRSIP.LE.0)IPRSIP=999                                                 
      WRITE(IOUT,130) IPRSIP                                                    
  130 FORMAT(1H ,52X,'SIP HEAD CHANGE PRINTOUT INTERVAL =',I9)                  
C                                                                               
C3------CHECK IF SPECIFIED VALUE OF WSEED SHOULD BE USED OR IF                  
C3------SEED SHOULD BE CALCULATED                                               
      IF(IPCALC.EQ.0) GO TO 150                                                 
C                                                                               
C3A-----CALCULATE SEED & ITERATION PARAMETERS PRIOR TO 1ST ITERATION            
      WRITE(IOUT,140)                                                           
  140 FORMAT(1H0,52X,'CALCULATE ITERATION PARAMETERS FROM MODEL',               
     1' CALCULATED WSEED')                                                      
      GO TO 1000                                                                
C                                                                               
C3B-----USE SPECIFIED VALUE OF WSEED                                            
C3B-----CALCULATE AND PRINT ITERATION PARAMETERS                                
  150 P1=-1.                                                                    
      P2=NPARM-1                                                                
      DO 160 I=1,NPARM                                                          
      P1=P1+1.                                                                  
  160 W(I)=1.-WSEED**(P1/P2)                                                    
      WRITE(IOUT,161) NPARM,WSEED,(W(J),J=1,NPARM)                              
  161 FORMAT(1H0,/,I5,' ITERATION PARAMETERS CALCULATED FROM',                  
     1     ' SPECIFIED WSEED =',F11.8,' :'//(10X,6E15.7))                       
C                                                                               
C4------RETURN                                                                  
 1000 RETURN                                                                    
      END                                                                       
      SUBROUTINE SIP1AP(HNEW,IBOUND,CR,CC,CV,HCOF,RHS,EL,FL,GL,V,               
     1      W,HDCG,LRCH,NPARM,KITER,HCLOSE,ACCL,ICNVG,KSTP,KPER,                
     2      IPCALC,IPRSIP,MXITER,NSTP,NCOL,NROW,NLAY,NODES,IOUT)                
C-----VERSION 1534 26OCT1984 SIP1AP                                             
C                                                                               
C     ******************************************************************        
C     SOLUTION BY THE STRONGLY IMPLICIT PROCEDURE -- 1 ITERATION                
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
      DOUBLE PRECISION HNEW,DITPAR,AC,HHCOF,RRHS,XI,DZERO,DONE,RES              
      DOUBLE PRECISION Z,B,D,E,F,H,S,AP,TP,CP,GP,UP,RP                          
      DOUBLE PRECISION ZHNEW,BHNEW,DHNEW,FHNEW,HHNEW,SHNEW                      
      DOUBLE PRECISION AL,BL,CL,DL,ELNCL,FLNCL,GLNCL                            
      DOUBLE PRECISION ELNRL,FLNRL,GLNRL,ELNLL,FLNLL,GLNLL                      
      DOUBLE PRECISION VNRL,VNCL,VNLL,ELXI,FLXI,GLXI,VN,HCFHNW                  
C                                                                               
      DIMENSION HNEW(NODES), IBOUND(NODES), CR(NODES), CC(NODES),               
     1  CV(NODES), HCOF(NODES), RHS(NODES), EL(NODES), FL(NODES),               
     2  GL(NODES), V(NODES), W(NPARM), HDCG(MXITER), LRCH(3,MXITER)             
C     ------------------------------------------------------------------        
C                                                                               
C1------CALCULATE ITERATION PARAMETERS IF FLAG IS SET.  THEN                    
C1------CLEAR THE FLAG SO THAT CALCULATION IS DONE ONLY ONCE.                   
      IF(IPCALC.NE.0)                                                           
     1     CALL SSIP1I(CR,CC,CV,IBOUND,NPARM,W,NCOL,NROW,NLAY,IOUT)             
      IPCALC=0                                                                  
C                                                                               
C2------ASSIGN VALUES TO FIELDS THAT ARE CONSTANT DURING AN ITERATION           
      DZERO=0.                                                                  
      DONE=1.                                                                   
      AC=ACCL                                                                   
      NRC=NROW*NCOL                                                             
      NTH=MOD(KITER-1,NPARM)+1                                                  
      DITPAR=W(NTH)                                                             
C                                                                               
C3------INITIALIZE VARIABLE THAT TRACKS MAXIMUM HEAD CHANGE DURING              
C3------THE ITERATION                                                           
      BIGG=0.                                                                   
C                                                                               
C4------CLEAR SIP WORK ARRAYS.                                                  
      DO 100 I=1,NODES                                                          
      EL(I)=0.                                                                  
      FL(I)=0.                                                                  
      GL(I)=0.                                                                  
  100 V(I)=0.                                                                   
C                                                                               
C5------SET NORMAL/REVERSE EQUATION ORDERING FLAG (1 OR -1) AND                 
C5------CALCULATE INDEXES DEPENDENT ON ORDERING                                 
      IDIR=1                                                                    
      IF(MOD(KITER,2).EQ.0)IDIR=-1                                              
      IDNRC=IDIR*NRC                                                            
      IDNCOL=IDIR*NCOL                                                          
C                                                                               
C6------STEP THROUGH CELLS CALCULATING INTERMEDIATE VECTOR V                    
C6------USING FORWARD SUBSTITUTION                                              
      DO 150 K=1,NLAY                                                           
      DO 150 I=1,NROW                                                           
      DO 150 J=1,NCOL                                                           
C                                                                               
C6A-----SET UP CURRENT CELL LOCATION INDEXES.  THESE ARE DEPENDENT              
C6A-----ON THE DIRECTION OF EQUATION ORDERING.                                  
      IF(IDIR.LE.0)GO TO 120                                                    
      II=I                                                                      
      JJ=J                                                                      
      KK=K                                                                      
      GO TO 122                                                                 
  120 II=NROW-I+1                                                               
      JJ=J                                                                      
      KK=NLAY-K+1                                                               
C                                                                               
C6B-----CALCULATE 1 DIMENSIONAL SUBSCRIPT OF CURRENT CELL AND                   
C6B-----SKIP CALCULATIONS IF CELL IS NOFLOW OR CONSTANT HEAD                    
  122 N=JJ+(II-1)*NCOL+(KK-1)*NRC                                               
      IF(IBOUND(N).LE.0)GO TO 150                                               
C                                                                               
C6C-----CALCULATE 1 DIMENSIONAL SUBSCRIPTS FOR LOCATING THE 6                   
C6C-----SURROUNDING CELLS                                                       
      NRN=N+IDNCOL                                                              
      NRL=N-IDNCOL                                                              
      NCN=N+1                                                                   
      NCL=N-1                                                                   
      NLN=N+IDNRC                                                               
      NLL=N-IDNRC                                                               
C                                                                               
C6D-----CALCULATE 1 DIMENSIONAL SUBSCRIPTS FOR CONDUCTANCE TO THE 6             
C6D-----SURROUNDING CELLS.  THESE DEPEND ON ORDERING OF EQUATIONS.              
      IF(IDIR.LE.0)GO TO 124                                                    
      NCF=N                                                                     
      NCD=NCL                                                                   
      NRB=NRL                                                                   
      NRH=N                                                                     
      NLS=N                                                                     
      NLZ=NLL                                                                   
      GO TO 126                                                                 
  124 NCF=N                                                                     
      NCD=NCL                                                                   
      NRB=N                                                                     
      NRH=NRN                                                                   
      NLS=NLN                                                                   
      NLZ=N                                                                     
C                                                                               
C6E-----ASSIGN VARIABLES IN MATRICES A & U INVOLVING ADJACENT CELLS             
C6E1----NEIGHBOR IS 1 ROW BACK                                                  
  126 B=DZERO                                                                   
      ELNRL=DZERO                                                               
      FLNRL=DZERO                                                               
      GLNRL=DZERO                                                               
      BHNEW=DZERO                                                               
      VNRL=DZERO                                                                
      IF(I.EQ.1) GO TO 128                                                      
      B=CC(NRB)                                                                 
      ELNRL=EL(NRL)                                                             
      FLNRL=FL(NRL)                                                             
      GLNRL=GL(NRL)                                                             
      BHNEW=B*HNEW(NRL)                                                         
      VNRL=V(NRL)                                                               
C                                                                               
C6E2----NEIGHBOR IS 1 ROW AHEAD                                                 
  128 H=DZERO                                                                   
      HHNEW=DZERO                                                               
      IF(I.EQ.NROW) GO TO 130                                                   
      H=CC(NRH)                                                                 
      HHNEW=H*HNEW(NRN)                                                         
C                                                                               
C6E3----NEIGHBOR IS 1 COLUMN BACK                                               
  130 D=DZERO                                                                   
      ELNCL=DZERO                                                               
      FLNCL=DZERO                                                               
      GLNCL=DZERO                                                               
      DHNEW=DZERO                                                               
      VNCL=DZERO                                                                
      IF(J.EQ.1) GO TO 132                                                      
      D=CR(NCD)                                                                 
      ELNCL=EL(NCL)                                                             
      FLNCL=FL(NCL)                                                             
      GLNCL=GL(NCL)                                                             
      DHNEW=D*HNEW(NCL)                                                         
      VNCL=V(NCL)                                                               
C                                                                               
C6E4----NEIGHBOR IS 1 COLUMN AHEAD                                              
  132 F=DZERO                                                                   
      FHNEW=DZERO                                                               
      IF(J.EQ.NCOL) GO TO 134                                                   
      F=CR(NCF)                                                                 
      FHNEW=F*HNEW(NCN)                                                         
C                                                                               
C6E5----NEIGHBOR IS 1 LAYER BEHIND                                              
  134 Z=DZERO                                                                   
      ELNLL=DZERO                                                               
      FLNLL=DZERO                                                               
      GLNLL=DZERO                                                               
      ZHNEW=DZERO                                                               
      VNLL=DZERO                                                                
      IF(K.EQ.1) GO TO 136                                                      
      Z=CV(NLZ)                                                                 
      ELNLL=EL(NLL)                                                             
      FLNLL=FL(NLL)                                                             
      GLNLL=GL(NLL)                                                             
      ZHNEW=Z*HNEW(NLL)                                                         
      VNLL=V(NLL)                                                               
C                                                                               
C6E6----NEIGHBOR IS 1 LAYER AHEAD                                               
  136 S=DZERO                                                                   
      SHNEW=DZERO                                                               
      IF(K.EQ.NLAY) GO TO 138                                                   
      S=CV(NLS)                                                                 
      SHNEW=S*HNEW(NLN)                                                         
C                                                                               
C6E7----CALCULATE THE NEGATIVE SUM OF ALL CONDUCTANCES TO NEIGHBORING           
C6E7----CELLS                                                                   
  138 E=-Z-B-D-F-H-S                                                            
C                                                                               
C6F-----CALCULATE COMPONENTS OF THE UPPER AND LOWER MATRICES, WHICH             
C6F-----ARE THE FACTORS OF MATRIX (A+B)                                         
      AL=Z/(DONE+DITPAR*(ELNLL+FLNLL))                                          
      BL=B/(DONE+DITPAR*(ELNRL+GLNRL))                                          
      CL=D/(DONE+DITPAR*(FLNCL+GLNCL))                                          
      AP=AL*ELNLL                                                               
      CP=BL*ELNRL                                                               
      GP=CL*FLNCL                                                               
      RP=CL*GLNCL                                                               
      TP=AL*FLNLL                                                               
      UP=BL*GLNRL                                                               
      HHCOF=HCOF(N)                                                             
      DL=E+HHCOF+DITPAR*(AP+TP+CP+GP+UP+RP)-AL*GLNLL-BL*FLNRL-CL*ELNCL          
      EL(N)=(F-DITPAR*(AP+CP))/DL                                               
      FL(N)=(H-DITPAR*(TP+GP))/DL                                               
      GL(N)=(S-DITPAR*(RP+UP))/DL                                               
C                                                                               
C6G-----CALCULATE THE RESIDUAL                                                  
      RRHS=RHS(N)                                                               
      HNW=HNEW(N)                                                               
      HCFHNW=HNW*HCOF(N)                                                        
      RES=RRHS-ZHNEW-BHNEW-DHNEW-E*HNEW(N)-HCFHNW-FHNEW-HHNEW-SHNEW             
C                                                                               
C6H-----CALCULATE THE INTERMEDIATE VECTOR V                                     
      V(N)=(AC*RES-AL*VNLL-BL*VNRL-CL*VNCL)/DL                                  
C                                                                               
  150 CONTINUE                                                                  
C                                                                               
C7------STEP THROUGH EACH CELL AND SOLVE FOR HEAD CHANGE BY BACK                
C7------SUBSTITUTION                                                            
      DO 160 K=1,NLAY                                                           
      DO 160 I=1,NROW                                                           
      DO 160 J=1,NCOL                                                           
C                                                                               
C7A-----SET UP CURRENT CELL LOCATION INDEXES.  THESE ARE DEPENDENT              
C7A-----ON THE DIRECTION OF EQUATION ORDERING.                                  
      IF(IDIR.LT.0) GO TO 152                                                   
      KK=NLAY-K+1                                                               
      II=NROW-I+1                                                               
      JJ=NCOL-J+1                                                               
      GO TO 154                                                                 
  152 KK=K                                                                      
      II=I                                                                      
      JJ=NCOL-J+1                                                               
C                                                                               
C7B-----CALCULATE 1 DIMENSIONAL SUBSCRIPT OF CURRENT CELL AND                   
C7B-----SKIP CALCULATIONS IF CELL IS NOFLOW OR CONSTANT HEAD                    
  154 N=JJ+(II-1)*NCOL+(KK-1)*NRC                                               
      IF(IBOUND(N).LE.0)GO TO 160                                               
C                                                                               
C7C-----CALCULATE 1 DIMENSIONAL SUBSCRIPTS FOR THE 3 NEIGHBORING CELLS          
C7C-----BEHIND (RELATIVE TO THE DIRECTION OF THE BACK SUBSTITUTION              
C7C-----ORDERING) THE CURRRENT CELL.                                            
      NC=N+1                                                                    
      NR=N+IDNCOL                                                               
      NL=N+IDNRC                                                                
C                                                                               
C7D-----BACK SUBSTITUTE, STORING HEAD CHANGE IN ARRAY V IN PLACE OF             
C7D-----INTERMEDIATE FORWARD SUBSTITUTION VALUES.                               
      ELXI=DZERO                                                                
      FLXI=DZERO                                                                
      GLXI=DZERO                                                                
      IF(JJ.NE.NCOL) ELXI=EL(N)*V(NC)                                           
      IF(I.NE.1) FLXI=FL(N)*V(NR)                                               
      IF(K.NE.1) GLXI=GL(N)*V(NL)                                               
      VN=V(N)                                                                   
      V(N)=VN-ELXI-FLXI-GLXI                                                    
C                                                                               
C7E-----GET THE ABSOLUTE HEAD CHANGE. IF IT IS MAX OVER GRID SO FAR.            
C7E-----THEN SAVE IT ALONG WITH CELL INDICES AND HEAD CHANGE.                   
      TCHK=ABS(V(N))                                                            
      IF (TCHK.LE.BIGG) GO TO 155                                               
      BIGG=TCHK                                                                 
      BIG=V(N)                                                                  
      IB=II                                                                     
      JB=JJ                                                                     
      KB=KK                                                                     
C                                                                               
C7F-----ADD HEAD CHANGE THIS ITERATION TO HEAD FROM THE PREVIOUS                
C7F-----ITERATION TO GET A NEW ESTIMATE OF HEAD.                                
  155 XI=V(N)                                                                   
      HNEW(N)=HNEW(N)+XI                                                        
C                                                                               
  160 CONTINUE                                                                  
C                                                                               
C8------STORE THE LARGEST ABSOLUTE HEAD CHANGE (THIS ITERATION) AND             
C8------AND ITS LOCATION.                                                       
      HDCG(KITER)=BIG                                                           
      LRCH(1,KITER)=KB                                                          
      LRCH(2,KITER)=IB                                                          
      LRCH(3,KITER)=JB                                                          
      ICNVG=0                                                                   
      IF(BIGG.LE.HCLOSE) ICNVG=1                                                
C                                                                               
C9------IF END OF TIME STEP, PRINT # OF ITERATIONS THIS STEP                    
      IF(ICNVG.EQ.0 .AND. KITER.NE.MXITER) GO TO 600                            
      IF(KSTP.EQ.1) WRITE(IOUT,500)                                             
  500 FORMAT(1H0)                                                               
      WRITE(IOUT,501) KITER,KSTP,KPER                                           
  501 FORMAT(1X,I5,' ITERATIONS FOR TIME STEP',I4,' IN STRESS PERIOD',          
     1        I3)                                                               
C                                                                               
C10-----PRINT HEAD CHANGE EACH ITERATION IF PRINTOUT INTERVAL IS REACHED        
      IF(ICNVG.EQ.0 .OR. KSTP.EQ.NSTP .OR. MOD(KSTP,IPRSIP).EQ.0)               
     1      CALL SSIP1P(HDCG,LRCH,KITER,KSTP,KPER,MXITER,IOUT)                  
C                                                                               
C11-----RETURN                                                                  
600   RETURN                                                                    
C                                                                               
      END                                                                       
      SUBROUTINE SSIP1I(CR,CC,CV,IBOUND,NPARM,W,NCOL,NROW,NLAY,                 
     1          IOUT)                                                           
C                                                                               
C-----VERSION 1743 26APR1983 SSIP1I                                             
C     ******************************************************************        
C     CALCULATE AN ITERATION PARAMETER SEED AND USE IT TO CALCULATE SIP         
C     ITERATION PARAMETERS                                                      
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
      DIMENSION CR(NCOL,NROW,NLAY),CC(NCOL,NROW,NLAY)                           
     1       ,CV(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY),W(NPARM)                
C                                                                               
      DOUBLE PRECISION DWMIN,AVGSUM                                             
C     ------------------------------------------------------------------        
C                                                                               
C1------CALCULATE CONSTANTS AND INITIALIZE VARIABLES                            
      PIEPIE=9.869604                                                           
      R=NROW                                                                    
      C=NCOL                                                                    
      L=NLAY                                                                    
      CCOL=PIEPIE/(2.*C*C)                                                      
      CROW=PIEPIE/(2.*R*R)                                                      
      CLAY=PIEPIE/(2.*L*L)                                                      
      WMINMN=1.                                                                 
      AVGSUM=0.                                                                 
      NODES=0                                                                   
C                                                                               
C2------LOOP THROUGH ALL CELLS, CALCULATING A SEED FOR EACH CELL                
C2------THAT IS ACTIVE                                                          
      DO 100 K=1,NLAY                                                           
      DO 100 I=1,NROW                                                           
      DO 100 J=1,NCOL                                                           
      IF(IBOUND(J,I,K).LE.0) GO TO 100                                          
C                                                                               
C2A-----CONDUCTANCE FROM THIS CELL                                              
C2A-----TO EACH OF THE 6 ADJACENT CELLS                                         
      D=0.                                                                      
      IF(J.NE.1) D=CR(J-1,I,K)                                                  
      F=0.                                                                      
      IF(J.NE.NCOL) F=CR(J,I,K)                                                 
      B=0.                                                                      
      IF(I.NE.1) B=CC(J,I-1,K)                                                  
      H=0.                                                                      
      IF(I.NE.NROW) H=CC(J,I,K)                                                 
      Z=0.                                                                      
      IF(K.NE.1) Z=CV(J,I,K-1)                                                  
      S=0.                                                                      
      IF(K.NE.NLAY) S=CV(J,I,K)                                                 
C                                                                               
C2B-----FIND THE MAXIMUM AND MINIMUM OF THE 2 CONDUCTANCE COEFFICIENTS          
C2B-----IN EACH PRINCIPAL COORDINATE DIRECTION                                  
      DFMX=AMAX1(D,F)                                                           
      BHMX=AMAX1(B,H)                                                           
      ZSMX=AMAX1(Z,S)                                                           
      DFMN=AMIN1(D,F)                                                           
      BHMN=AMIN1(B,H)                                                           
      ZSMN=AMIN1(Z,S)                                                           
      IF(DFMN.EQ.0.) DFMN=DFMX                                                  
      IF(BHMN.EQ.0.) BHMN=BHMX                                                  
      IF(ZSMN.EQ.0.) ZSMN=ZSMX                                                  
C                                                                               
C2C-----CALCULATE A SEED IN EACH PRINCIPAL COORDINATE DIRECTION                 
      WCOL=1.                                                                   
      IF(DFMN.NE.0.) WCOL=CCOL/(1.+(BHMX+ZSMX)/DFMN)                            
      WROW=1.                                                                   
      IF(BHMN.NE.0.) WROW=CROW/(1.+(DFMX+ZSMX)/BHMN)                            
      WLAY=1.                                                                   
      IF(ZSMN.NE.0.) WLAY=CLAY/(1.+(DFMX+BHMX)/ZSMN)                            
C                                                                               
C2D-----SELECT THE CELL SEED, WHICH IS THE MINIMUM SEED OF THE 3.               
C2D-----SELECT THE MINIMUM SEED OVER THE WHOLE GRID.                            
      WMIN=AMIN1(WCOL,WROW,WLAY)                                                
      WMINMN=AMIN1(WMINMN,WMIN)                                                 
C                                                                               
C2E-----ADD THE CELL SEED TO THE ACCUMULATOR AVGSUM FOR USE                     
C2E-----IN GETTING THE AVERAGE SEED.                                            
      DWMIN=WMIN                                                                
      AVGSUM=AVGSUM+DWMIN                                                       
      NODES=NODES+1                                                             
C                                                                               
  100 CONTINUE                                                                  
C                                                                               
C3------CALCULATE THE AVERAGE SEED OF THE CELL SEEDS, AND PRINT                 
C3------THE AVERAGE AND MINIMUM SEEDS.                                          
      TMP=NODES                                                                 
      AVGMIN=AVGSUM                                                             
      AVGMIN=AVGMIN/TMP                                                         
      WRITE(IOUT,101) AVGMIN,WMINMN                                             
  101 FORMAT(1H0,'AVERAGE SEED =',F11.8/1X,'MINIMUM SEED =',F11.8)              
C                                                                               
C4------CALCULATE AND PRINT ITERATION PARAMETERS FROM THE AVERAGE SEED          
      P1=-1.                                                                    
      P2=NPARM-1                                                                
      DO 50 I=1,NPARM                                                           
      P1=P1+1.                                                                  
   50 W(I)=1.-AVGMIN**(P1/P2)                                                   
      WRITE(IOUT,150) NPARM,(W(J),J=1,NPARM)                                    
  150 FORMAT(1H0,/,I5,' ITERATION PARAMETERS CALCULATED FROM',                  
     1      ' AVERAGE SEED:'//(10X,6E15.7))                                     
C                                                                               
C5------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE SSIP1P(HDCG,LRCH,KITER,KSTP,KPER,MXITER,IOUT)                  
C                                                                               
C                                                                               
C-----VERSION 1504 08DEC1982 SSIP1P                                             
C     ******************************************************************        
C     PRINT MAXIMUM HEAD CHANGE FOR EACH ITERATION DURING A TIME STEP           
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
C                                                                               
      DIMENSION HDCG(MXITER), LRCH(3,MXITER)                                    
C     ------------------------------------------------------------------        
C                                                                               
      WRITE(IOUT,5)                                                             
    5 FORMAT(1H0,'MAXIMUM HEAD CHANGE FOR EACH ITERATION:'/                     
     1    1H0,5(' HEAD CHANGE LAYER,ROW,COL')/1X,132('-'))                      
      WRITE (IOUT,10) (HDCG(J),(LRCH(I,J),I=1,3),J=1,KITER)                     
   10 FORMAT((1X,5(G12.4,' (',I3,',',I3,',',I3,')')))                           
      WRITE(IOUT,11)                                                            
   11 FORMAT(1H0)                                                               
C                                                                               
      RETURN                                                                    
C                                                                               
      END                                                                       
C FILE SOR1.F77                                                                 
      SUBROUTINE SOR1AL(ISUM,LENX,LCA,LCRES,LCHDCG,LCLRCH,LCIEQP,               
     1       MXITER,NCOL,NROW,NLAY,NSLICE,MBW,IN,IOUT)                          
C                                                                               
C-----VERSION 1003 08DEC1983 SOR1AL                                             
C     ******************************************************************        
C     ALLOCATE STORAGE FOR SOR ARRAYS                                           
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
C     ------------------------------------------------------------------        
C                                                                               
C1------PRINT A MESSAGE IDENTIFYING SOR PACKAGE                                 
      WRITE(IOUT,1)IN                                                           
    1 FORMAT(1H0,'SOR1 -- SLICE-SUCCESSIVE OVERRELAXATION PACKAGE'              
     1,', VERSION 1, 12/08/83',' INPUT READ FROM UNIT',I3)                      
C                                                                               
C2------READ AND PRINT MXITER (MAXIMUM # OF ITERATIONS)                         
      READ(IN,2) MXITER                                                         
    2 FORMAT(I10)                                                               
      WRITE(IOUT,3) MXITER                                                      
    3 FORMAT(1X,I5,' ITERATIONS ALLOWED FOR SOR CLOSURE')                       
C                                                                               
C3------ALLOCATE SPACE FOR THE SOR ARRAYS                                       
      ISOLD=ISUM                                                                
      NSLICE=NCOL*NLAY                                                          
      MBW=NLAY+1                                                                
      LCA=ISUM                                                                  
      ISUM=ISUM+NSLICE*MBW                                                      
      LCRES=ISUM                                                                
      ISUM=ISUM+NSLICE                                                          
      LCIEQP=ISUM                                                               
      ISUM=ISUM+NSLICE                                                          
      LCHDCG=ISUM                                                               
      ISUM=ISUM+MXITER                                                          
      LCLRCH=ISUM                                                               
      ISUM=ISUM+3*MXITER                                                        
      ISP=ISUM-ISOLD                                                            
C                                                                               
C4------CALCULATE AND PRINT THE SPACE USED IN THE X ARRAY                       
      WRITE(IOUT,4) ISP                                                         
    4 FORMAT(1X,I6,' ELEMENTS IN X ARRAY ARE USED BY SOR')                      
      ISUM1=ISUM-1                                                              
      WRITE(IOUT,5) ISUM1,LENX                                                  
    5 FORMAT(1X,I6,' 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                                                                               
C5------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE SOR1RP(MXITER,ACCL,HCLOSE,IN,IPRSOR,IOUT)                      
C                                                                               
C                                                                               
C-----VERSION 1005 16MAR1983 SOR1RP                                             
C     ******************************************************************        
C     READ PARAMETERS FOR SOR                                                   
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
C     ------------------------------------------------------------------        
C                                                                               
C1------READ THE ACCELERATION PARAMETER/RELAXATION FACTOR (ACCL) THE            
C1------CLOSURE CRITERION (HCLOSE) AND THE NUMBER OF TIME STEPS                 
C1------BETWEEN PRINTOUTS OF MAXIMUM HEAD CHANGES (IPRSOR).                     
      READ(IN,1) ACCL,HCLOSE,IPRSOR                                             
    1 FORMAT(2F10.0,I10)                                                        
      IF(ACCL.EQ.0.) ACCL=1.                                                    
      IF(IPRSOR.LT.1) IPRSOR=999                                                
C                                                                               
C2------PRINT ACCL, HCLOSE, IPRSOR                                              
      WRITE(IOUT,100)                                                           
  100 FORMAT(1H0,///57X,'SOLUTION BY SLICE-SUCCESSIVE OVERRELAXATION'           
     1/57X,43('-'))                                                             
      WRITE(IOUT,115) MXITER                                                    
  115 FORMAT(1H0,47X,'MAXIMUM ITERATIONS ALLOWED FOR CLOSURE =',I9)             
      WRITE(IOUT,120) ACCL                                                      
  120 FORMAT(1H ,63X,'ACCELERATION PARAMETER =',G15.5)                          
      WRITE(IOUT,125) HCLOSE                                                    
  125 FORMAT(1H ,52X,'HEAD CHANGE CRITERION FOR CLOSURE =',E15.5)               
      WRITE(IOUT,130) IPRSOR                                                    
  130 FORMAT(1H ,52X,'SOR HEAD CHANGE PRINTOUT INTERVAL =',I9)                  
C                                                                               
C3------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE SOR1AP(HNEW,IBOUND,CR,CC,CV,HCOF,RHS,A,RES,IEQPNT,             
     1      HDCG,LRCH,KITER,HCLOSE,ACCL,ICNVG,KSTP,KPER,                        
     2      IPRSOR,MXITER,NSTP,NCOL,NROW,NLAY,NSLICE,MBW,IOUT)                  
C-----VERSION 0936 09MAY1983 SOR1AP                                             
C     ******************************************************************        
C     SOLUTION BY SLICE-SUCCESSIVE OVERRELAXATION -- 1 ITERATION                
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
      DOUBLE PRECISION HNEW,DIFF,DP,EE,R,HCFHNW,HHCOF                           
C                                                                               
      DIMENSION HNEW(NCOL,NROW,NLAY), IBOUND(NCOL,NROW,NLAY),                   
     1   CR(NCOL,NROW,NLAY), CC(NCOL,NROW,NLAY),                                
     1   CV(NCOL,NROW,NLAY), HCOF(NCOL,NROW,NLAY), RHS(NCOL,NROW,NLAY),         
     2   HDCG(MXITER), LRCH(3,MXITER),A(MBW,NSLICE),RES(NSLICE),                
     3   IEQPNT(NLAY,NCOL)                                                      
C     ------------------------------------------------------------------        
C                                                                               
C1------CALCULATE # OF ELEMENTS IN COMPRESSED MATRIX A AND                      
C1------INITIALIZE FIELDS TO SAVE LARGEST HEAD CHANGE.                          
      NA=MBW*NSLICE                                                             
      BIG=0.                                                                    
      ABSBIG=0.                                                                 
      IB=0                                                                      
      JB=0                                                                      
      KB=0                                                                      
C                                                                               
C2------PROCESS EACH SLICE                                                      
      DO 500 I=1,NROW                                                           
C                                                                               
C3------CLEAR A                                                                 
      DO 110 J=1,NSLICE                                                         
      DO 110 K=1,MBW                                                            
  110 A(K,J)=0.                                                                 
C                                                                               
C4------ASSIGN A SEQUENCE # TO EACH VARIABLE HEAD CELL.                         
      NEQT=0                                                                    
      DO 200 J=1,NCOL                                                           
      DO 200 K=1,NLAY                                                           
      IEQPNT(K,J)=0                                                             
      IF(IBOUND(J,I,K).LE.0) GO TO 200                                          
      NEQT=NEQT+1                                                               
      IEQPNT(K,J)=NEQT                                                          
  200 CONTINUE                                                                  
C                                                                               
C5------FOR EACH CELL LOAD MATRIX A AND VECTOR RES                              
      DO 300 J=1,NCOL                                                           
      DO 300 K=1,NLAY                                                           
C                                                                               
C5A-----IF SEQUENCE # IS 0 (CELL IS EXTERNAL) GO ON TO NEXT CELL                
      NEQ=IEQPNT(K,J)                                                           
      IF(NEQ.EQ.0) GO TO 300                                                    
C                                                                               
C5B-----INITIALIZE ACCUMULATORS EE AND R                                        
      EE=0.                                                                     
      R=RHS(J,I,K)                                                              
C                                                                               
C5C-----IF NODE TO LEFT SUBTRACT TERMS FROM EE AND R                            
      IF(J.EQ.1) GO TO 120                                                      
      DP=CR(J-1,I,K)                                                            
      R=R-DP*HNEW(J-1,I,K)                                                      
      EE=EE-DP                                                                  
C                                                                               
C5D-----IF NODE TO RIGHT SUBTRACT TERMS FROM EE & R, MOVE COND TO A             
  120 IF(J.EQ.NCOL) GO TO 125                                                   
      SP=CR(J,I,K)                                                              
      DP=SP                                                                     
      R=R-DP*HNEW(J+1,I,K)                                                      
      EE=EE-DP                                                                  
      NXT=IEQPNT(K,J+1)                                                         
      IF(NXT.GT.0) A(1+NXT-NEQ,NEQ)=SP                                          
C                                                                               
C5E-----IF NODE TO REAR SUBTRACT TERMS FROM EE AND R                            
  125 IF(I.EQ.1) GO TO 130                                                      
      DP=CC(J,I-1,K)                                                            
      R=R-DP*HNEW(J,I-1,K)                                                      
      EE=EE-DP                                                                  
C                                                                               
C5F-----IF NODE TO FRONT SUBTRACT TERMS FROM EE AND R                           
  130 IF(I.EQ.NROW) GO TO 132                                                   
      DP=CC(J,I,K)                                                              
      R=R-DP*HNEW(J,I+1,K)                                                      
      EE=EE-DP                                                                  
C                                                                               
C5G-----IF NODE ABOVE SUBTRACT TERMS FROM EE AND R                              
  132 IF(K.EQ.1) GO TO 134                                                      
      DP=CV(J,I,K-1)                                                            
      R=R-DP*HNEW(J,I,K-1)                                                      
      EE=EE-DP                                                                  
C                                                                               
C5H-----IF NODE BELOW SUBTRACT TERMS FROM EE & R AND MOVE COND TO A             
  134 IF(K.EQ.NLAY) GO TO 136                                                   
      SP=CV(J,I,K)                                                              
      DP=SP                                                                     
      R=R-DP*HNEW(J,I,K+1)                                                      
      EE=EE-DP                                                                  
      IF(IEQPNT(K+1,J).GT.0) A(2,NEQ)=SP                                        
C                                                                               
C5I-----MOVE EE INTO A, SUBTRACT EE TIMES LAST HEAD FROM R TO GET RES           
  136 HHCOF=HCOF(J,I,K)                                                         
      A(1,NEQ)=EE+HHCOF                                                         
      HNW=HNEW(J,I,K)                                                           
      HCFHNW=HNW*HCOF(J,I,K)                                                    
      RES(NEQ)=R-EE*HNEW(J,I,K)-HCFHNW                                          
  300 CONTINUE                                                                  
C                                                                               
C6------IF NO EQUATIONS GO TO NEXT SLICE, IF ONE EQUATION SOLVE                 
C6------DIRECTLY, IF 2 EQUATIONS CALL SSOR1B TO SOLVE FOR FIRST                 
C6------ESTIMATE OF HEAD CHANGE FOR THIS ITERATION.                             
      IF(NEQT.LT.1) GO TO 500                                                   
      IF(NEQT.EQ.1) RES(1)=RES(1)/A(1,1)                                        
      IF(NEQT.GE.2) CALL SSOR1B(A,RES,NEQT,NA,MBW)                              
C                                                                               
C7------FOR EACH CELL IN SLICE CALCULATE FINAL HEAD CHANGE THEN HEAD.           
      DO 400 J=1,NCOL                                                           
      DO 400 K=1,NLAY                                                           
      NEQ=IEQPNT(K,J)                                                           
      IF(NEQ.EQ.0) GO TO 400                                                    
C                                                                               
C7A-----MULTIPLY FIRST ESTIMATE OF HEAD CHANGE BY RELAX FACTOR TO               
C7A-----GET FINAL ESTIMATE OF HEAD CHANGE FOR THIS ITERATION.                   
      DH=RES(NEQ)*ACCL                                                          
      DIFF=DH                                                                   
C                                                                               
C7B-----ADD FINAL ESTIMATE TO HEAD FROM LAST ITERATION TO GET HEAD              
C7B-----FOR THIS ITERATION.                                                     
      HNEW(J,I,K)=HNEW(J,I,K)+DIFF                                              
C                                                                               
C7C-----SAVE FINAL HEAD CHANGE IF IT IS THE LARGEST                             
      ABSDH=ABS(DH)                                                             
      IF(ABSDH.LE.ABSBIG) GO TO 400                                             
      ABSBIG=ABSDH                                                              
      BIG=DH                                                                    
      IB=I                                                                      
      JB=J                                                                      
      KB=K                                                                      
  400 CONTINUE                                                                  
C                                                                               
C                                                                               
  500 CONTINUE                                                                  
C                                                                               
C8------SAVE LARGEST HEAD CHANGE FOR THIS ITERATION                             
      HDCG(KITER)=BIG                                                           
      LRCH(1,KITER)=KB                                                          
      LRCH(2,KITER)=IB                                                          
      LRCH(3,KITER)=JB                                                          
C                                                                               
C9------IF LARGEST HEAD CHANGE IS SMALLER THAN CLOSURE THEN SET                 
C9------CONVERGE FLAG (ICNVG) EQUAL TO 1.                                       
      ICNVG=0                                                                   
      IF(ABSBIG.LE.HCLOSE) ICNVG=1                                              
C                                                                               
C10-----IF NOT CONVERGED AND NOT EXCEDED ITERATIONS THEN RETURN                 
      IF(ICNVG.EQ.0 .AND. KITER.NE.MXITER) RETURN                               
      IF(KSTP.EQ.1) WRITE(IOUT,600)                                             
  600 FORMAT(1H0)                                                               
C                                                                               
C11-----PRINT NUMBER OF ITERATIONS                                              
      WRITE(IOUT,601) KITER,KSTP,KPER                                           
  601 FORMAT(1X,I5,' ITERATIONS FOR TIME STEP',I4,' IN STRESS PERIOD',          
     1        I3)                                                               
C                                                                               
C12-----IF FAILED TO CONVERGE OR LAST TIME STEP OR PRINTOUT                     
C12-----INTERVAL SPECIFIED BY USER IS HERE THEN PRINT MAXIMUM                   
C12-----HEAD CHANGES FOR EACH ITERATION.                                        
      IF(ICNVG.NE.0 .AND. KSTP.NE.NSTP .AND. MOD(KSTP,IPRSOR).NE.0)             
     1      GO TO 700                                                           
      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) (HDCG(J),(LRCH(I,J),I=1,3),J=1,KITER)                     
   10 FORMAT((1X,4(4X,G12.4,' (',I3,',',I3,',',I3,')')))                        
      WRITE(IOUT,11)                                                            
   11 FORMAT(1H0)                                                               
C                                                                               
C13-----RETURN                                                                  
  700 RETURN                                                                    
C                                                                               
      END                                                                       
      SUBROUTINE SSOR1B(A,B,N,NA,MBW)                                           
C                                                                               
C                                                                               
C-----VERSION 1359 31MAR1983 SSOR1B                                             
C     ******************************************************************        
C     SOLVE A SYMMETRIC SET OF EQUATIONS                                        
C        A IS COEFFICIENT MATRIX IN COMPRESSED FORM                             
C        B IS RIGHT HAND SIDE AND IS REPLACED BY SOLUTION                       
C        N IS NUMBER OF EQUATIONS TO BE SOLVED                                  
C        MBW IS BANDWIDTH OF A                                                  
C        NA IS ONE-DIMENSION SIZE OF A                                          
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
      DIMENSION A(NA),B(N)                                                      
C     ------------------------------------------------------------------        
C                                                                               
      NM1=N-1                                                                   
      MBW1=MBW-1                                                                
      ID=1-MBW                                                                  
C                                                                               
C1------SEQUENTIALLY USE EACH OF THE FIRST N-1 ROWS AS                          
C1------THE PIVOT ROW.                                                          
      DO 20 I=1,NM1                                                             
C                                                                               
C2------CALCULATE THE INVERSE OF THE PIVOT.                                     
      ID=ID+MBW                                                                 
      C1=1./A(ID)                                                               
      LD=ID                                                                     
      L=I                                                                       
C                                                                               
C3------FOR EACH ROW AFTER THE PIVOT ROW (THE TARGET ROW)                       
C3------ELIMINATE THE COLUMN CORRESPONDING TO THE PIVOT.                        
      DO 15 J=1,MBW1                                                            
      L=L+1                                                                     
      IF(L.GT.N) GO TO 20                                                       
      IB=ID+J                                                                   
C                                                                               
C4------CALCULATE THE FACTOR NEEDED TO ELIMINATE A TERM IN THE                  
C4------TARGET ROW.                                                             
      C=A(IB)*C1                                                                
      LD=LD+MBW                                                                 
      LB=LD-1                                                                   
C                                                                               
C5------MODIFY THE REST OF THE TERMS IN THE TARGET ROW.                         
      DO 10 K=J,MBW1                                                            
C                                                                               
C6------SUBTRACT THE FACTOR TIMES A TERM IN THE PIVOT ROW                       
C6------FROM THE CORRESPONDING COLUMN IN THE TARGET ROW.                        
      LB=LB+1                                                                   
      A(LB)=A(LB)-C*A(ID+K)                                                     
   10 CONTINUE                                                                  
C                                                                               
C7------MODIFY THE RIGHT SIDE OF THE EQUATION CORRESPONDING                     
C7------TO THE TARGET ROW.                                                      
      B(I+J)=B(I+J)-C*B(I)                                                      
   15 CONTINUE                                                                  
   20 CONTINUE                                                                  
      ID=ID+MBW                                                                 
C                                                                               
C8------SOLVE THE LAST EQUATION.                                                
      B(N)=B(N)/A(ID)                                                           
C                                                                               
C9------WORKING BACKWARDS SOLVE THE REST OF THE EQUATIONS.                      
      DO 70 I=1,NM1                                                             
      ID=ID-MBW                                                                 
C                                                                               
C10-----CLEAR THE ACCUMULATOR SUM.                                              
      SUM=0.0                                                                   
      L=N-I                                                                     
      MBW1M=MIN0(MBW1,I)                                                        
C                                                                               
C11-----ADD THE KNOWN TERMS IN EQUATION L TO SUM.                               
      DO 60 J=1,MBW1M                                                           
      SUM=SUM+A(ID+J)*B(L+J)                                                    
   60 CONTINUE                                                                  
C                                                                               
C12-----SOLVE FOR THE ONE UNKNOWN IN EQUATION L.                                
      B(L)=(B(L)-SUM)/A(ID)                                                     
   70 CONTINUE                                                                  
C                                                                               
C13-----RETURN                                                                  
      RETURN                                                                    
      END                                                                       
