C FILE WEL1.F77                                                                 
      SUBROUTINE WEL1AL(ISUM,LENX,LCWELL,MXWELL,NWELLS,IN,IOUT,                 
     1                       IWELCB)                                            
C                                                                               
C-----VERSION 1747 24APR1985                                                    
C     ******************************************************************        
C     ALLOCATE ARRAY STORAGE FOR WELL PACKAGE                                   
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
C     ------------------------------------------------------------------        
C                                                                               
C1------IDENTIFY PACKAGE AND INITIALIZE NWELLS                                  
      WRITE(IOUT,1)IN                                                           
    1 FORMAT(1H0,'WEL1 -- WELL PACKAGE, VERSION 1, 04/24/85',                   
     2' INPUT READ FROM',I3)                                                    
      NWELLS=0                                                                  
C                                                                               
C2------READ MAX NUMBER OF WELLS AND                                            
C2------UNIT OR FLAG FOR CELL-BY-CELL FLOW TERMS.                               
      READ(IN,2) MXWELL,IWELCB                                                  
    2 FORMAT(2I10)                                                              
      WRITE(IOUT,3) MXWELL                                                      
    3 FORMAT(1H ,'MAXIMUM OF',I5,' WELLS')                                      
      IF(IWELCB.GT.0) WRITE(IOUT,9) IWELCB                                      
    9 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE RECORDED ON UNIT',I3)               
      IF(IWELCB.LT.0) WRITE(IOUT,8)                                             
    8 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE PRINTED WHEN ICBCFL NOT 0')         
C                                                                               
C3------SET LCWELL EQUAL TO LOCATION OF WELL LIST IN X ARRAY.                   
      LCWELL=ISUM                                                               
C                                                                               
C4------ADD AMOUNT OF SPACE USED BY WELL LIST TO ISUM.                          
      ISP=4*MXWELL                                                              
      ISUM=ISUM+ISP                                                             
C                                                                               
C5------PRINT NUMBER OF SPACES IN X ARRAY USED BY WELL PACKAGE.                 
      WRITE(IOUT,4) ISP                                                         
    4 FORMAT(1X,I6,' ELEMENTS IN X ARRAY ARE USED FOR WELLS')                   
      ISUM1=ISUM-1                                                              
      WRITE(IOUT,5) ISUM1,LENX                                                  
    5 FORMAT(1X,I6,' ELEMENTS OF X ARRAY USED OUT OF',I7)                       
C                                                                               
C6------IF THERE ISN'T ENOUGH SPACE IN THE X ARRAY THEN PRINT                   
C6------A WARNING MESSAGE.                                                      
      IF(ISUM1.GT.LENX) WRITE(IOUT,6)                                           
    6 FORMAT(1X,'   ***X ARRAY MUST BE DIMENSIONED LARGER***')                  
C7------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE WEL1RP(WELL,NWELLS,MXWELL,IN,IOUT)                             
C                                                                               
C                                                                               
C-----VERSION 1544 22DEC1982 WEL1RP                                             
C     ******************************************************************        
C     READ NEW WELL LOCATIONS AND STRESS RATES                                  
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
      DIMENSION WELL(4,MXWELL)                                                  
C     ------------------------------------------------------------------        
C                                                                               
C1------READ ITMP(NUMBER OF WELLS OR FLAG SAYING REUSE WELL DATA)               
      READ (IN,1) ITMP                                                          
    1 FORMAT(I10)                                                               
      IF(ITMP.GE.0) GO TO 50                                                    
C                                                                               
C1A-----IF ITMP LESS THAN ZERO REUSE DATA. PRINT MESSAGE AND RETURN.            
      WRITE(IOUT,6)                                                             
    6 FORMAT(1H0,'REUSING WELLS FROM LAST STRESS PERIOD')                       
      RETURN                                                                    
C                                                                               
C1B-----ITMP=>0.  SET NWELLS EQUAL TO ITMP.                                     
   50 NWELLS=ITMP                                                               
      IF(NWELLS.LE.MXWELL) GO TO 100                                            
C                                                                               
C2------NWELLS>MXWELL.  PRINT MESSAGE. STOP.                                    
      WRITE(IOUT,99) NWELLS,MXWELL                                              
   99 FORMAT(1H0,'NWELLS(',I4,') IS GREATER THAN MXWELL(',I4,')')               
      STOP                                                                      
C                                                                               
C3------PRINT NUMBER OF WELLS IN CURRENT STRESS PERIOD.                         
  100 WRITE (IOUT,2) NWELLS                                                     
    2 FORMAT(1H0,10X,I4,' WELLS')                                               
C                                                                               
C4------IF THERE ARE NO ACTIVE WELLS IN THIS STRESS PERIOD THEN RETURN          
      IF(NWELLS.EQ.0) GO TO 260                                                 
C                                                                               
C5------READ AND PRINT LAYER,ROW,COLUMN AND RECHARGE RATE.                      
      WRITE(IOUT,3)                                                         
    3 FORMAT(1H ,47X,'LAYER    ROW    COL    STRESS RATE   WELL NO.'/           
     1,48X,45('-'))                                                             
      DO 250 II=1,NWELLS                                                        
      READ (IN,4) K,I,J,Q                                                       
    4 FORMAT(3I10,F10.0)                                                        
      WRITE (IOUT,5) K,I,J,Q,II                                             
    5 FORMAT(48X,I3,I8,I7,G16.5,I8)                                             
      WELL(1,II)=K                                                              
      WELL(2,II)=I                                                              
      WELL(3,II)=J                                                              
      WELL(4,II)=Q                                                              
  250 CONTINUE                                                                  
C                                                                               
C6------RETURN                                                                  
  260 RETURN                                                                    
      END                                                                       
      SUBROUTINE WEL1FM(NWELLS,MXWELL,RHS,WELL,IBOUND,                          
     1        NCOL,NROW,NLAY)                                                   
C                                                                               
C-----VERSION 1027 10APR1985 WEL1FM                                             
C                                                                               
C     ******************************************************************        
C     SUBSTRACT Q FROM RHS                                              
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
      DIMENSION RHS(NCOL,NROW,NLAY),WELL(4,MXWELL),                             
     1            IBOUND(NCOL,NROW,NLAY)                                        
C     ------------------------------------------------------------------        
C1------IF NUMBER OF WELLS <= 0 THEN RETURN.                                    
      IF(NWELLS.LE.0) RETURN                                                    
C                                                                               
C2------PROCESS EACH WELL IN THE WELL LIST.                                     
      DO 100 L=1,NWELLS                                                         
      IR=WELL(2,L)                                                              
      IC=WELL(3,L)                                                              
      IL=WELL(1,L)                                                              
      Q=WELL(4,L)                                                               
C                                                                               
C2A-----IF THE CELL IS INACTIVE THEN BYPASS PROCESSING.                         
      IF(IBOUND(IC,IR,IL).LE.0) GO TO 100                                       
C                                                                               
C2B-----IF THE CELL IS VARIABLE HEAD THEN ADD RECHARGE RATE                     
C       TO THE RHS ACCUMULATOR.                                                 
      RHS(IC,IR,IL)=RHS(IC,IR,IL)-Q                                             
  100 CONTINUE                                                                  
C                                                                               
C3------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE WEL1BD(NWELLS,MXWELL,VBNM,VBVL,MSUM,WELL,IBOUND,DELT,          
     1             NCOL,NROW,NLAY,KSTP,KPER,IWELCB,ICBCFL,BUFF,IOUT)            
C                                                                               
C-----VERSION 1449 20MAY1983 WEL1BD                                             
C     ******************************************************************        
C     CALCULATE VOLUMETRIC BUDGET FOR WELLS                                     
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------
      CHARACTER*4 VBNM,TEXT        
      DIMENSION VBNM(4,MSUM),VBVL(4,MSUM),WELL(4,MXWELL),                       
     1          IBOUND(NCOL,NROW,NLAY),BUFF(NCOL,NROW,NLAY)                     
      DIMENSION TEXT(4)                                                         
C                                                                               
      DATA TEXT(1),TEXT(2),TEXT(3),TEXT(4) /'    ','    ','   W','ELLS'/        
C     ------------------------------------------------------------------        
C                                                                               
C1------CLEAR RATIN AND RATOUT ACCUMULATORS.                                    
      RATIN=0.                                                                  
      RATOUT=0.                                                                 
      IBD=0                                                                     
C                                                                               
C2------IF THERE ARE NO WELLS DO NOT ACCUMULATE FLOW                            
      IF(NWELLS.EQ.0) GO TO 200                                                 
C                                                                               
C3------TEST TO SEE IF CELL-BY-CELL FLOW TERMS WILL BE RECORDED.                
      IF(ICBCFL.EQ.0 .OR. IWELCB.LE.0) GO TO 60                                 
C                                                                               
C4------IF CELL-BY-CELL FLOWS WILL BE SAVED THEN CLEAR THE BUFFER.              
      IBD=1                                                                     
      DO 50 IL=1,NLAY                                                           
      DO 50 IR=1,NROW                                                           
      DO 50 IC=1,NCOL                                                           
      BUFF(IC,IR,IL)=0.                                                         
   50 CONTINUE                                                                  
C                                                                               
C5------PROCESS WELLS ONE AT A TIME.                                            
   60 DO 100 L=1,NWELLS                                                         
      IR=WELL(2,L)                                                              
      IC=WELL(3,L)                                                              
      IL=WELL(1,L)                                                              
      Q=WELL(4,L)                                                               
C                                                                               
C5A-----IF THE CELL IS EXTERNAL IGNORE IT.                                      
      IF(IBOUND(IC,IR,IL).LE.0)GO TO 100                                        
C                                                                               
C5B-----PRINT THE INDIVIDUAL RATES IF REQUESTED(IWELCB<0).                      
      IF(IWELCB.LT.0.AND.ICBCFL.NE.0) WRITE(IOUT,900) (TEXT(N),N=1,4),          
     1    KPER,KSTP,L,IL,IR,IC,Q                                                
  900 FORMAT(1H0,4A4,'   PERIOD',I3,'   STEP',I3,'   WELL',I4,                  
     1    '   LAYER',I3,'   ROW ',I4,'   COL',I4,'   RATE',G15.7)               
C                                                                               
C5C-----IF CELL-BY-CELL FLOWS ARE TO BE SAVED THEN ADD THEM TO BUFFER.          
      IF(IBD.EQ.1) BUFF(IC,IR,IL)=BUFF(IC,IR,IL)+Q                              
      IF(Q) 90,100,80                                                           
C                                                                               
C5D-----PUMPING RATE IS POSITIVE(RECHARGE). ADD IT TO RATIN.                    
   80 RATIN=RATIN+Q                                                             
      GO TO 100                                                                 
C                                                                               
C5E-----PUMPING RATE IS NEGATIVE(DISCHARGE). ADD IT TO RATOUT.                  
   90 RATOUT=RATOUT-Q                                                           
  100 CONTINUE                                                                  
C                                                                               
C6------IF CELL-BY-CELL FLOWS WILL BE SAVED CALL UBUDSV TO RECORD THEM          
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,IWELCB,BUFF,NCOL,NROW,            
     1                          NLAY,IOUT)                                      
C                                                                               
C7------MOVE RATES INTO VBVL FOR PRINTING BY MODULE BAS1OT.                     
  200 VBVL(3,MSUM)=RATIN                                                        
      VBVL(4,MSUM)=RATOUT                                                       
C                                                                               
C8------MOVE RATES TIMES TIME STEP LENGTH INTO VBVL ACCUMULATORS.               
      VBVL(1,MSUM)=VBVL(1,MSUM)+RATIN*DELT                                      
      VBVL(2,MSUM)=VBVL(2,MSUM)+RATOUT*DELT                                     
C                                                                               
C9------MOVE BUDGET TERM LABELS INTO VBNM FOR PRINTING.                         
      VBNM(1,MSUM)=TEXT(1)                                                      
      VBNM(2,MSUM)=TEXT(2)                                                      
      VBNM(3,MSUM)=TEXT(3)                                                      
      VBNM(4,MSUM)=TEXT(4)                                                      
C                                                                               
C10-----INCREMENT BUDGET TERM COUNTER(MSUM).                                    
      MSUM=MSUM+1                                                               
C                                                                               
C11-----RETURN                                                                  
      RETURN                                                                    
      END                                                                       
C FILE RCH1.F77                                                                 
      SUBROUTINE RCH1AL(ISUM,LENX,LCIRCH,LCRECH,NRCHOP,                         
     C                  NCOL,NROW,IN,IOUT,IRCHCB)                               
C                                                                               
C-----VERSION 1744 24APR1985 RCH1AL                                             
C     ******************************************************************        
C     ALLOCATE ARRAY STORAGE FOR RECHARGE                                       
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
C     ------------------------------------------------------------------        
C                                                                               
C1------IDENTIFY PACKAGE.                                                       
      WRITE(IOUT,1)IN                                                           
    1 FORMAT(1H0,'RCH1 -- RECHARGE PACKAGE, VERSION 1, 04/24/85',               
     2' INPUT READ FROM UNIT',I3)                                               
C                                                                               
C2------READ NRCHOP AND IRCHCB.                                                 
      READ(IN,2)NRCHOP,IRCHCB                                                   
    2 FORMAT(2I10)                                                              
C                                                                               
C3------CHECK TO SEE THAT OPTION IS LEGAL.                                      
      IF(NRCHOP.GE.1.AND.NRCHOP.LE.3)GO TO 200                                  
C                                                                               
C3A-----IF ILLEGAL PRINT A MESSAGE AND ABORT SIMULATION                         
      WRITE(IOUT,8)                                                             
    8 FORMAT(1X,'ILLEGAL OPTION CODE. SIMULATION ABORTING')                     
      STOP                                                                      
C                                                                               
C4------IF OPTION IS LEGAL PRINT OPTION CODE.                                   
  200 IRK=ISUM                                                                  
      IF(NRCHOP.EQ.1) WRITE(IOUT,201)                                           
  201 FORMAT(1X,'OPTION 1 -- RECHARGE TO TOP LAYER')                            
      IF(NRCHOP.EQ.2) WRITE(IOUT,202)                                           
  202 FORMAT(1X,'OPTION 2 -- RECHARGE TO ONE SPECIFIED NODE IN EACH',           
     1        ' VERTICAL COLUMN')                                               
      IF(NRCHOP.EQ.3) WRITE(IOUT,203)                                           
  203 FORMAT(1X,'OPTION 3 -- RECHARGE TO HIGHEST ACTIVE NODE IN EACH',          
     1        ' VERTICAL COLUMN')                                               
C                                                                               
C5------IF CELL-BY-CELL FLOW  TERMS TO BE SAVED THEN PRINT UNIT #               
      IF(IRCHCB.GT.0) WRITE(IOUT,204) IRCHCB                                    
  204 FORMAT(1X,'CELL-BY-CELL FLOW TERMS WILL BE RECORDED ON UNIT',I3)          
C                                                                               
C6------ALLOCATE SPACE FOR THE RECHARGE ARRAY(RECH).                            
      LCRECH=ISUM                                                               
      ISUM=ISUM+NCOL*NROW                                                       
C                                                                               
C7------IF OPTION 2 THEN ALLOCATE SPACE FOR INDICATOR ARRAY(IRCH)               
      LCIRCH=ISUM                                                               
      IF(NRCHOP.NE.2)GO TO 300                                                  
      ISUM=ISUM+NCOL*NROW                                                       
C                                                                               
C8------CALCULATE AND PRINT AMOUNT OF SPACE USED BY RECHARGE.                   
  300 IRK=ISUM-IRK                                                              
      WRITE(IOUT,4)IRK                                                          
    4 FORMAT(1X,I6,' ELEMENTS OF X ARRAY USED FOR RECHARGE')                    
      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 MADE LARGER***')                         
C                                                                               
C9------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE RCH1RP(NRCHOP,IRCH,RECH,DELR,DELC,NROW,NCOL,                   
     C                  NLAY,IN,IOUT)                                           
C                                                                               
C-----VERSION 1513 22DEC1982 RCH1RP                                             
C     ******************************************************************        
C     READ RECHARGE RATES                                                       
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------
      CHARACTER*4 ANAME        
      DIMENSION IRCH(NCOL,NROW),RECH(NCOL,NROW),                                
     1          ANAME(6,2),DELR(NCOL),DELC(NROW)                                
C                                                                               
      DATA ANAME(1,1),ANAME(2,1),ANAME(3,1),ANAME(4,1),ANAME(5,1),              
     1  ANAME(6,1) /'    ','RECH','ARGE',' LAY','ER I','NDEX'/                  
      DATA ANAME(1,2),ANAME(2,2),ANAME(3,2),ANAME(4,2),ANAME(5,2),              
     1  ANAME(6,2) /'    ','    ','    ','    ','RECH','ARGE'/                  
C     ------------------------------------------------------------------        
C                                                                               
C1------READ FLAGS SHOWING WHETHER DATA IS TO BE REUSED.                        
      READ(IN,4)INRECH,INIRCH                                                   
    4 FORMAT(2I10)                                                              
C                                                                               
C2------TEST INRECH TO SEE WHERE RECH IS COMING FROM.                           
      IF(INRECH.GE.0)GO TO 32                                                   
C                                                                               
C2A-----IF INRECH<0 THEN REUSE RECHARGE ARRAY FROM LAST STRESS PERIOD           
      WRITE(IOUT,3)                                                             
    3 FORMAT(1H0,'REUSING RECH FROM LAST STRESS PERIOD')                        
      GO TO 55                                                                  
C                                                                               
C3------IF INRECH=>0 THEN CALL U2DREL TO READ RECHARGE RATE.                    
   32 CALL U2DREL(RECH,ANAME(1,2),NROW,NCOL,0,IN,IOUT)                          
C                                                                               
C4------MULTIPLY RECHARGE RATE BY CELL AREA TO GET VOLUMETRIC RATE.             
      DO 50 IR=1,NROW                                                           
      DO 50 IC=1,NCOL                                                           
      RECH(IC,IR)=RECH(IC,IR)*DELR(IC)*DELC(IR)                                 
   50 CONTINUE                                                                  
C                                                                               
C5------IF NRCHOP=2 THEN A LAYER INDICATOR ARRAY IS NEEDED.                     
  55  IF (NRCHOP.NE.2)GO TO 60                                                  
C                                                                               
C6------IF INIRCH<0 THEN REUSE LAYER INDICATOR ARRAY.                           
      IF(INIRCH.GE.0)GO TO 58                                                   
      WRITE(IOUT,2)                                                             
    2 FORMAT(1H0,'REUSING IRCH FROM LAST STRESS PERIOD')                        
      GO TO 60                                                                  
C                                                                               
C7------IF INIRCH=>0 CALL U2DINT TO READ LAYER IND ARRAY(IRCH)                  
   58 CALL U2DINT(IRCH,ANAME(1,1),NROW,NCOL,0,IN,IOUT)                          
C                                                                               
C8------RETURN                                                                  
   60 RETURN                                                                    
      END                                                                       
      SUBROUTINE RCH1FM(NRCHOP,IRCH,RECH,RHS,IBOUND,NCOL,                       
     1                     NROW,NLAY)                                           
C                                                                               
C-----VERSION 1025 10APR1985 RCH1FM                                             
C     ******************************************************************        
C     ADD RECHARGE TO RHS                                                       
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
      DIMENSION IRCH(NCOL,NROW),RECH(NCOL,NROW),                                
     1          RHS(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY)                      
C     ------------------------------------------------------------------        
C                                                                               
C1------IF NRCHOP IS 1 RECHARGE IS IN TOP LAYER. LAYER INDEX IS 1.              
      IF(NRCHOP.NE.1) GO TO 15                                                  
C                                                                               
      DO 10 IR=1,NROW                                                           
      DO 10 IC=1,NCOL                                                           
C                                                                               
C1A-----IF CELL IS EXTERNAL THERE IS NO RECHARGE INTO IT.                       
      IF(IBOUND(IC,IR,1).LE.0)GO TO 10                                          
C                                                                               
C1B-----SUBTRACT RECHARGE RATE FROM RIGHT-HAND-SIDE.                            
      RHS(IC,IR,1)=RHS(IC,IR,1)-RECH(IC,IR)                                     
   10 CONTINUE                                                                  
      GO TO 100                                                                 
C                                                                               
C2------IF OPTION IS 2 THEN RECHARGE IS INTO LAYER IN INDICATOR ARRAY           
   15 IF(NRCHOP.NE.2)GO TO 25                                                   
      DO 20 IR=1,NROW                                                           
      DO 20 IC=1,NCOL                                                           
C                                                                               
C2A-----LAYER INDEX IS IN INDICATOR ARRAY.                                      
      IL=IRCH(IC,IR)                                                            
C                                                                               
C2B-----IF THE CELL IS EXTERNAL THERE IS NO RECHARGE INTO IT.                   
      IF(IBOUND(IC,IR,IL).LE.0)GO TO 20                                         
C                                                                               
C2C-----SUBTRACT RECHARGE FROM RIGHT-HAND-SIDE.                                 
      RHS(IC,IR,IL)=RHS(IC,IR,IL)-RECH(IC,IR)                                   
   20 CONTINUE                                                                  
      GO TO 100                                                                 
C                                                                               
C3------IF OPTION IS 3 RECHARGE IS INTO HIGHEST INTERNAL CELL.                  
   25 IF(NRCHOP.NE.3)GO TO 100                                                  
C        CANNOT PASS THROUGH CONSTANT HEAD NODE                                 
      DO 30 IR=1,NROW                                                           
      DO 30 IC=1,NCOL                                                           
      DO 28 IL=1,NLAY                                                           
C                                                                               
C3A-----IF CELL IS CONSTANT HEAD MOVE ON TO NEXT HORIZONTAL LOCATION.           
      IF(IBOUND(IC,IR,IL).LT.0) GO TO 30                                        
C                                                                               
C3B-----IF CELL IS INACTIVE MOVE DOWN A LAYER.                                  
      IF (IBOUND(IC,IR,IL).EQ.0)GO TO 28                                        
C                                                                               
C3C-----SUBTRACT RECHARGE FROM RIGHT-HAND-SIDE.                                 
      RHS(IC,IR,IL)=RHS(IC,IR,IL)-RECH(IC,IR)                                   
      GO TO 30                                                                  
   28 CONTINUE                                                                  
   30 CONTINUE                                                                  
  100 CONTINUE                                                                  
C                                                                               
C4------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE RCH1BD(NRCHOP,IRCH,RECH,IBOUND,NROW,NCOL,NLAY,                 
     1    DELT,VBVL,VBNM,MSUM,KSTP,KPER,IRCHCB,ICBCFL,BUFF,IOUT)                
C                                                                               
C-----VERSION 1533 22DEC1982 RCH1BD                                             
C     ******************************************************************        
C     CALCULATE VOLUMETRIC BUDGET FOR RECHARGE                                  
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
      CHARACTER*4 VBNM,TEXT
      DIMENSION IRCH(NCOL,NROW),RECH(NCOL,NROW),
     1          IBOUND(NCOL,NROW,NLAY),BUFF(NCOL,NROW,NLAY),                    
     2          VBVL(4,20),VBNM(4,20)                                           
      DIMENSION TEXT(4)                                                         
      DATA TEXT(1),TEXT(2),TEXT(3),TEXT(4) /'    ','    ','RECH','ARGE'/        
C     ------------------------------------------------------------------        
C                                                                               
C1------CLEAR THE RATE ACCUMULATORS.                                            
      RATIN=0.                                                                  
      RATOUT=0.                                                                 
C                                                                               
C2------IF CELL-BY-CELL FLOW TERMS WILL BE SAVED THEN CLEAR THE BUFFER.         
      IBD=0                                                                     
      IF(ICBCFL.EQ.0 .OR. IRCHCB.LE.0) GO TO 5                                  
      IBD=1                                                                     
      DO 2 IL=1,NLAY                                                            
      DO 2 IR=1,NROW                                                            
      DO 2 IC=1,NCOL                                                            
      BUFF(IC,IR,IL)=0.                                                         
    2 CONTINUE                                                                  
C                                                                               
C3------IF NRCHOP=1 RECH GOES INTO LAYER 1. PROCESS EACH HORIZONTAL             
C3------CELL LOCATION.                                                          
    5 IF(NRCHOP.NE.1) GO TO 15                                                  
C                                                                               
C     ---RECHARGE IS APPLIED TO TOP LAYER                                       
      DO 10 IR=1,NROW                                                           
      DO 10 IC=1,NCOL                                                           
C                                                                               
C3A-----IF CELL IS EXTERNAL THEN DO NOT DO BUDGET FOR IT.                       
      IF(IBOUND(IC,IR,1).LE.0)GO TO 10                                          
      Q=RECH(IC,IR)                                                             
C                                                                               
C3B-----IF CELL-BY-CELL FLOW TERMS WILL BE SAVED THEN ADD RECH TO BUFF          
      IF(IBD.EQ.1) BUFF(IC,IR,1)=Q                                              
C                                                                               
C3C-----IF RECH POSITIVE ADD IT TO RATIN ELSE ADD IT TO RATOUT.                 
      IF(Q) 8,10,7                                                              
    7 RATIN=RATIN+Q                                                             
      GO TO 10                                                                  
    8 RATOUT=RATOUT-Q                                                           
   10 CONTINUE                                                                  
      GO TO 100                                                                 
C                                                                               
C4------IF NRCHOP=2 RECH IS IN LAYER SHOWN IN INDICATOR ARRAY(IRCH).            
C4------PROCESS HORIZONTAL CELL LOCATIONS ONE AT A TIME.                        
   15 IF(NRCHOP.NE.2)GO TO 25                                                   
      DO 20 IR=1,NROW                                                           
      DO 20 IC=1,NCOL                                                           
C                                                                               
C4A-----GET LAYER INDEX FROM INDICATOR ARRAY(IRCH).                             
      IL=IRCH(IC,IR)                                                            
C                                                                               
C4B-----IF CELL IS EXTERNAL DO NOT CALCULATE BUDGET FOR IT.                     
      IF(IBOUND(IC,IR,IL).LE.0)GO TO 20                                         
      Q=RECH(IC,IR)                                                             
C                                                                               
C4C-----IF C-B-C FLOW TERMS WILL BE SAVED THEN ADD RECHARGE TO BUFFER.          
      IF(IBD.EQ.1) BUFF(IC,IR,IL)=Q                                             
C                                                                               
C4D-----IF RECHARGE IS POSITIVE ADD TO RATIN ELSE ADD IT TO RATOUT.             
      IF(Q) 18,20,17                                                            
   17 RATIN=RATIN+Q                                                             
      GO TO 20                                                                  
   18 RATOUT=RATOUT-Q                                                           
   20 CONTINUE                                                                  
      GO TO 100                                                                 
C                                                                               
C5------IF OPTION=3 RECHARGE IS INTO HIGHEST INTERNAL CELL. IT WILL NOT         
C5------PASS THROUGH A CONSTANT HEAD CELL. PROCESS HORIZONTAL CELL              
C5------LOCATIONS ONE AT A TIME.                                                
   25 IF(NRCHOP.NE.3)GO TO 100                                                  
      DO 30 IR=1,NROW                                                           
      DO 30 IC=1,NCOL                                                           
      DO 28 IL=1,NLAY                                                           
C                                                                               
C5A-----IF CELL IS CONSTANT HEAD MOVE ON TO NEXT HORIZONTAL LOCATION.           
      IF(IBOUND(IC,IR,IL).LT.0) GO TO 30                                        
C                                                                               
C5B-----IF CELL IS INACTIVE MOVE DOWN TO NEXT CELL.                             
      IF (IBOUND(IC,IR,IL).EQ.0)GO TO 28                                        
      Q=RECH(IC,IR)                                                             
C                                                                               
C5C-----IF C-B-C FLOW TERMS TO BE SAVED THEN ADD RECHARGE TO BUFFER.            
      IF(IBD.EQ.1) BUFF(IC,IR,IL)=Q                                             
C                                                                               
C5D-----IF RECH IS POSITIVE ADD IT TO RATIN ELSE ADD IT TO RATOUT.              
      IF(Q) 27,30,26                                                            
   26 RATIN=RATIN+Q                                                             
      GO TO 30                                                                  
   27 RATOUT=RATOUT-Q                                                           
      GO TO 30                                                                  
   28 CONTINUE                                                                  
   30 CONTINUE                                                                  
C                                                                               
  100 CONTINUE                                                                  
C                                                                               
C6------IF C-B-C FLOW TERMS TO BE SAVED CALL UBUDSV TO WRITE THEM.              
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,IRCHCB,BUFF,NCOL,NROW,            
     1                          NLAY,IOUT)                                      
C                                                                               
C7------MOVE TOTAL RECHARGE RATE INTO VBVL FOR PRINTING BY BAS1OT.              
      VBVL(4,MSUM)=RATOUT                                                       
      VBVL(3,MSUM)=RATIN                                                        
C                                                                               
C8------ADD RECHARGE FOR TIME STEP TO RECHARGE ACCUMULATOR IN VBVL.             
      VBVL(2,MSUM)=VBVL(2,MSUM)+RATOUT*DELT                                     
      VBVL(1,MSUM)=VBVL(1,MSUM)+RATIN*DELT                                      
C                                                                               
C9------MOVE BUDGET TERM LABELS TO VBNM FOR PRINT BY MODULE BAS_OT.             
      VBNM(1,MSUM)=TEXT(1)                                                      
      VBNM(2,MSUM)=TEXT(2)                                                      
      VBNM(3,MSUM)=TEXT(3)                                                      
      VBNM(4,MSUM)=TEXT(4)                                                      
C                                                                               
C10-----INCREMENT BUDGET TERM COUNTER.                                          
      MSUM=MSUM+1                                                               
C                                                                               
C11-----RETURN                                                                  
      RETURN                                                                    
      END                                                                       
C FILE RIV1.F77                                                                 
      SUBROUTINE RIV1AL(ISUM,LENX,LCRIVR,MXRIVR,NRIVER,IN,IOUT,                 
     1                      IRIVCB)                                             
C                                                                               
C-----VERSION 0935 08DEC1983 RIV1AL                                             
C     ******************************************************************        
C     ALLOCATE ARRAY STORAGE FOR RIVERS                                         
C     ******************************************************************        
C                                                                               
C     SPECIFICATIONS:                                                           
C     ------------------------------------------------------------------        
C     ------------------------------------------------------------------        
C                                                                               
C1------IDENTIFY PACKAGE AND INITIALIZE NRIVER.                                 
      WRITE(IOUT,1)IN                                                           
    1 FORMAT(1H0,'RIV1 -- RIVER PACKAGE, VERSION 1,  12/08/83',                 
     2' INPUT READ FROM UNIT',I3)                                               
      NRIVER=0                                                                  
C                                                                               
C2------READ & PRINT MXRIVR & IRIVCB(UNIT OR FLAG FOR C-B-C FLOWS)              
      READ(IN,2)MXRIVR,IRIVCB                                                   
    2 FORMAT(2I10)                                                              
      WRITE(IOUT,3)MXRIVR                                                       
    3 FORMAT(1H ,'MAXIMUM OF',I5,' RIVER NODES')                                
      IF(IRIVCB.GT.0) WRITE(IOUT,9) IRIVCB                                      
    9 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE RECORDED ON UNIT',I3)               
      IF(IRIVCB.LT.0) WRITE(IOUT,8)                                             
    8 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE PRINTED')                           
C                                                                               
C3------SET LCRIVR EQUAL TO ADDRESS OF FIRST UNUSED SPACE IN X.                 
      LCRIVR=ISUM                                                               
C                                                                               
C4------CALCULATE AMOUNT OF SPACE USED BY RIVER LIST.                           
      ISP=6*MXRIVR                                                              
      ISUM=ISUM+ISP                                                             
C                                                                               
C5------PRINT AMOUNT OF SPACE USED BY RIVER PACKAGE.                            
      WRITE (IOUT,4)ISP                                                         
    4 FORMAT(1X,I6,' ELEMENTS IN X ARRAY ARE USED FOR RIVERS')                  
      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                                                                               
C7------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE RIV1RP(RIVR,NRIVER,MXRIVR,IN,IOUT)                             
C                                                                               
C                                                                               
C-----VERSION 1319 25AUG1982 RIV1RP                                             
C     ******************************************************************        
C     READ RIVER HEAD, CONDUCTANCE AND BOTTOM ELEVATION                         
C     ******************************************************************        
C                                                                               
C     SPECIFICATIONS:                                                           
C     ------------------------------------------------------------------        
      DIMENSION RIVR(6,MXRIVR)                                                  
C     ------------------------------------------------------------------        
C                                                                               
C                                                                               
C1------READ ITMP(NUMBER OF RIVER REACHES OR FLAG TO REUSE DATA)                
      READ(IN,8)ITMP                                                            
    8 FORMAT(I10)                                                               
C                                                                               
C2------TEST ITMP.                                                              
      IF(ITMP.GE.0)GO TO 50                                                     
C                                                                               
C2A-----IF ITMP <0 THEN REUSE DATA FROM LAST STRESS PERIOD.                     
      WRITE(IOUT,7)                                                             
    7 FORMAT(1H0,'REUSING RIVER REACHES FROM LAST STRESS PERIOD')               
      GO TO 260                                                                 
C                                                                               
C3------IF ITMP=> ZERO THEN IT IS THE NUMBER OF RIVER REACHES                   
   50 NRIVER=ITMP                                                               
C                                                                               
C4------IF NRIVER>MXRIVR THEN STOP.                                             
      IF(NRIVER.LE.MXRIVR)GO TO 100                                             
      WRITE(IOUT,99)NRIVER,MXRIVR                                               
   99 FORMAT(1H0,'NRIVER(',I4,') IS GREATER THAN MXRIVR(',I4,')')               
C                                                                               
C4A-----ABNORMAL STOP.                                                          
      STOP                                                                      
C                                                                               
C5------PRINT NUMBER OF RIVER REACHES IN THIS STRESS PERIOD.                    
  100 WRITE(IOUT,1)NRIVER                                                       
    1 FORMAT(1H0,//1X,I5,' RIVER REACHES')                                      
C                                                                               
C6------IF THERE ARE NO RIVER REACHES THEN RETURN.                              
      IF(NRIVER.EQ.0) GO TO 260                                                 
C                                                                               
C7------READ AND PRINT DATA FOR EACH RIVER REACH.                               
      WRITE(IOUT,3)                                                         
    3 FORMAT(1H0,15X,'LAYER',5X,'ROW',5X,'COL   '                               
     1,'  STAGE     CONDUCTANCE   BOTTOM ELEVATION   RIVER REACH'               
     2/1X,15X,80('-'))                                                          
      DO 250 II=1,NRIVER                                                        
      READ(IN,4)K,I,J,RIVR(4,II),RIVR(5,II),RIVR(6,II)                          
    4 FORMAT(3I10,3F10.0)                                                       
      WRITE(IOUT,5)K,I,J,RIVR(4,II),RIVR(5,II),RIVR(6,II),II                
    5 FORMAT(1X,15X,I4,I9,I8,G13.4,G14.4,G19.4,I10)                             
      RIVR(1,II)=K                                                              
      RIVR(2,II)=I                                                              
      RIVR(3,II)=J                                                              
  250 CONTINUE                                                                  
C                                                                               
C8------RETURN                                                                  
  260 RETURN                                                                    
      END                                                                       
      SUBROUTINE RIV1FM(NRIVER,MXRIVR,RIVR,HNEW,HCOF,RHS,IBOUND,                
     1                  NCOL,NROW,NLAY)                                         
C                                                                               
C-----VERSION 0915 27AUG1982 RIV1FM                                             
C     ******************************************************************        
C     ADD RIVER TERMS TO RHS AND HCOF                                           
C     ******************************************************************        
C                                                                               
C     SPECIFICATIONS:                                                           
C     ------------------------------------------------------------------        
C                                                                               
      DOUBLE PRECISION HNEW                                                     
      DIMENSION RIVR(6,MXRIVR),HNEW(NCOL,NROW,NLAY),                            
     1          HCOF(NCOL,NROW,NLAY),RHS(NCOL,NROW,NLAY),                       
     2          IBOUND(NCOL,NROW,NLAY)                                          
C     ------------------------------------------------------------------        
C                                                                               
C                                                                               
C1------IF NRIVER<=0 THERE ARE NO RIVERS. RETURN.                               
      IF(NRIVER.LE.0)RETURN                                                     
C                                                                               
C2------PROCESS EACH CELL IN THE RIVER LIST.                                    
      DO 100 L=1,NRIVER                                                         
C                                                                               
C3------GET COLUMN, ROW, AND LAYER OF CELL CONTAINING REACH                     
      IL=RIVR(1,L)                                                              
      IR=RIVR(2,L)                                                              
      IC=RIVR(3,L)                                                              
C                                                                               
C4------IF THE CELL IS EXTERNAL SKIP IT.                                        
      IF(IBOUND(IC,IR,IL).LE.0)GO TO 100                                        
C                                                                               
C5------SINCE THE CELL IS INTERNAL GET THE RIVER DATA.                          
      HRIV=RIVR(4,L)                                                            
      CRIV=RIVR(5,L)                                                            
      RBOT=RIVR(6,L)                                                            
      HHNEW=HNEW(IC,IR,IL)                                                      
C                                                                               
C6------COMPARE AQUIFER HEAD TO BOTTOM OF STREAM BED.                           
      IF(HHNEW.LE.RBOT)GO TO 96                                                 
C                                                                               
C7------SINCE HEAD>BOTTOM ADD TERMS TO RHS AND HCOF.                            
      RHS(IC,IR,IL)=RHS(IC,IR,IL)-CRIV*HRIV                                     
      HCOF(IC,IR,IL)=HCOF(IC,IR,IL)-CRIV                                        
      GO TO 100                                                                 
C                                                                               
C8------SINCE HEAD<BOTTOM ADD TERM ONLY TO RHS.                                 
   96 RHS(IC,IR,IL)=RHS(IC,IR,IL)-CRIV*(HRIV-RBOT)                              
  100 CONTINUE                                                                  
C                                                                               
C9------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE RIV1BD(NRIVER,MXRIVR,RIVR,IBOUND,HNEW,                         
     1      NCOL,NROW,NLAY,DELT,VBVL,VBNM,MSUM,KSTP,KPER,IRIVCB,                
     2      ICBCFL,BUFF,IOUT)                                                   
C-----VERSION 1256 28DEC1983 RIV1BD                                             
C     ******************************************************************        
C     CALCULATE VOLUMETRIC BUDGET FOR RIVERS                                    
C     ******************************************************************        
C                                                                               
C     SPECIFICATIONS:                                                           
C     ------------------------------------------------------------------        
      DOUBLE PRECISION HNEW
      CHARACTER*4 VBNM,TEXT
      DIMENSION RIVR(6,MXRIVR),IBOUND(NCOL,NROW,NLAY),                          
     1          HNEW(NCOL,NROW,NLAY),VBVL(4,20),VBNM(4,20),                     
     2          BUFF(NCOL,NROW,NLAY)                                            
      DIMENSION TEXT(4)                                                         
      DATA TEXT(1),TEXT(2),TEXT(3),TEXT(4) /'   R','IVER',' LEA','KAGE'/        
C     ------------------------------------------------------------------        
C                                                                               
C1------INITIALIZE CELL-BY-CELL FLOW TERM FLAG (IBD) AND                        
C1------ACCUMULATORS (RATIN AND RATOUT).                                        
      IBD=0                                                                     
      RATIN=0.                                                                  
      RATOUT=0.                                                                 
C                                                                               
C2------IF NO REACHES KEEP ZEROES IN ACCUMULATORS.                              
      IF(NRIVER.EQ.0)GO TO 200                                                  
C                                                                               
C3------TEST TO SEE IF CELL-BY-CELL FLOW TERMS ARE NEEDED.                      
      IF(ICBCFL.EQ.0  .OR. IRIVCB.LE.0 ) GO TO 10                               
C                                                                               
C3A-----CELL-BY-CELL FLOW TERMS ARE NEEDED SET IBD AND CLEAR BUFFER.            
      IBD=1                                                                     
      DO 5 IL=1,NLAY                                                            
      DO 5 IR=1,NROW                                                            
      DO 5 IC=1,NCOL                                                            
      BUFF(IC,IR,IL)=0.                                                         
    5 CONTINUE                                                                  
C                                                                               
C4------FOR EACH RIVER REACH ACCUMULATE RIVER FLOW (STEPS 5-15)                 
   10 DO 100 L=1,NRIVER                                                         
C                                                                               
C5------GET LAYER, ROW & COLUMN OF CELL CONTAINING REACH.                       
      IL=RIVR(1,L)                                                              
      IR=RIVR(2,L)                                                              
      IC=RIVR(3,L)                                                              
C                                                                               
C6------IF CELL IS EXTERNAL MOVE ON TO NEXT REACH.                              
      IF(IBOUND(IC,IR,IL).LE.0)GO TO 100                                        
C                                                                               
C7------GET RIVER PARAMETERS FROM RIVER LIST.                                   
      HRIV=RIVR(4,L)                                                            
      CRIV=RIVR(5,L)                                                            
      RBOT=RIVR(6,L)                                                            
      HHNEW=HNEW(IC,IR,IL)                                                      
C                                                                               
C8------COMPARE HEAD IN AQUIFER TO BOTTOM OF RIVERBED.                          
C                                                                               
C9------AQUIFER HEAD > BOTTOM THEN RATE=CRIV*(HRIV-HNEW).                       
      IF(HHNEW.GT.RBOT)RATE=CRIV*(HRIV-HHNEW)                                   
C                                                                               
C10-----AQUIFER HEAD < BOTTOM THEN RATE=CRIV*(HRIV-RBOT)                        
      IF(HHNEW.LE.RBOT)RATE=CRIV*(HRIV-RBOT)                                    
C                                                                               
C11-----PRINT THE INDIVIDUAL RATES IF REQUESTED(IRIVCB<0).                      
      IF(IRIVCB.LT.0.AND.ICBCFL.NE.0) WRITE(IOUT,900) (TEXT(N),N=1,4),          
     1    KPER,KSTP,L,IL,IR,IC,RATE                                             
  900 FORMAT(1H0,4A4,'   PERIOD',I3,'   STEP',I3,'   REACH',I4,                 
     1    '   LAYER',I3,'   ROW',I4,'   COL',I4,'   RATE',G15.7)                
C                                                                               
C12------IF C-B-C FLOW TERMS ARE TO BE SAVED THEN ADD RATE TO BUFFER.           
      IF(IBD.EQ.1) BUFF(IC,IR,IL)=BUFF(IC,IR,IL)+RATE                           
C                                                                               
C13-----SEE IF FLOW IS INTO AQUIFER OR INTO RIVER.                              
      IF(RATE)94,100,96                                                         
C                                                                               
C14-----AQUIFER IS DISCHARGING TO RIVER SUBTRACT RATE FROM RATOUT.              
   94 RATOUT=RATOUT-RATE                                                        
      GO TO 100                                                                 
C                                                                               
C15-----AQUIFER IS RECHARGED FROM RIVER ADD RATE TO RATIN.                      
   96 RATIN=RATIN+RATE                                                          
  100 CONTINUE                                                                  
C                                                                               
C16-----IF C-B-C FLOW TERMS WILL BE SAVED CALL UBUDSV TO RECORD THEM.           
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,IRIVCB,BUFF,NCOL,NROW,            
     1                          NLAY,IOUT)                                      
C                                                                               
C17-----MOVE RATES,VOLUMES & LABELS INTO ARRAYS FOR PRINTING.                   
  200 VBVL(3,MSUM)=RATIN                                                        
      VBVL(4,MSUM)=RATOUT                                                       
      VBVL(1,MSUM)=VBVL(1,MSUM)+RATIN*DELT                                      
      VBVL(2,MSUM)=VBVL(2,MSUM)+RATOUT*DELT                                     
      VBNM(1,MSUM)=TEXT(1)                                                      
      VBNM(2,MSUM)=TEXT(2)                                                      
      VBNM(3,MSUM)=TEXT(3)                                                      
      VBNM(4,MSUM)=TEXT(4)                                                      
C                                                                               
C18-----INCREMENT BUDGET TERM COUNTER                                           
      MSUM=MSUM+1                                                               
C                                                                               
C19-----RETURN                                                                  
      RETURN                                                                    
      END                                               
