C FILE DRN1.F77                                                                 
      SUBROUTINE DRN1AL(ISUM,LENX,LCDRAI,NDRAIN,MXDRN,IN,IOUT,                  
     1                           IDRNCB)                                        
C                                                                               
C-----VERSION 1740 24APR1985 DRN1AL                                             
C     ******************************************************************        
C     ALLOCATE ARRAY STORAGE FOR DRAIN PACKAGE                                  
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
C     ------------------------------------------------------------------        
C                                                                               
C1------IDENTIFY PACKAGE AND INITIALIZE NDRAIN.                                 
      WRITE(IOUT,1)IN                                                           
    1 FORMAT(1H0,'DRN1 -- DRAIN PACKAGE, VERSION 1, 04/24/85',                  
     2' INPUT READ FROM UNIT',I3)                                               
      NDRAIN=0                                                                  
C                                                                               
C2------READ & PRINT MXDRN & IDRNCB(UNIT & FLAG FOR CELL-BY-CELL FLOW)          
      READ(IN,2) MXDRN,IDRNCB                                                   
    2 FORMAT(2I10)                                                              
      WRITE(IOUT,3) MXDRN                                                       
    3 FORMAT(1H ,'MAXIMUM OF',I5,' DRAINS')                                     
      IF(IDRNCB.GT.0) WRITE(IOUT,9) IDRNCB                                      
    9 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE RECORDED ON UNIT',I3)               
      IF(IDRNCB.LT.0) WRITE(IOUT,8)                                             
    8 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE PRINTED WHEN ICBCFL NOT 0')         
C                                                                               
C3------SET LCDRAI EQUAL TO ADDRESS OF FIRST UNUSED SPACE IN X.                 
      LCDRAI=ISUM                                                               
C                                                                               
C4------CALCULATE AMOUNT OF SPACE USED BY THE DRAIN PACKAGE.                    
      ISP=5*MXDRN                                                               
      ISUM=ISUM+ISP                                                             
C                                                                               
C5------PRINT AMOUNT OF SPACE USED BY DRAIN PACKAGE.                            
      WRITE(IOUT,4) ISP                                                         
    4 FORMAT(1X,I6,' ELEMENTS IN X ARRAY ARE USED FOR DRAINS')                  
      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                                                                               
C6------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE DRN1RP(DRAI,NDRAIN,MXDRN,IN,IOUT)                              
C                                                                               
C                                                                               
C-----VERSION 1603 25APR1983 DRN1RP                                             
C     ******************************************************************        
C     READ DRAIN LOCATIONS, ELEVATIONS, AND CONDUCTANCES                        
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
      DIMENSION DRAI(5,MXDRN)                                                   
C     ------------------------------------------------------------------        
C                                                                               
C1------READ ITMP(NUMBER OF DRAIN CELLS 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 DRAINS FROM LAST STRESS PERIOD')                      
      RETURN                                                                    
C                                                                               
C3------IF ITMP=>0 THEN IT IS THE NUMBER OF DRAINS.                             
   50 NDRAIN=ITMP                                                               
      IF(NDRAIN.LE.MXDRN) GO TO 100                                             
C                                                                               
C4------IF NDRAIN>MXDRN THEN STOP                                               
      WRITE(IOUT,99) NDRAIN,MXDRN                                               
   99 FORMAT(1H0,'NDRAIN(',I4,') IS GREATER THAN MXDRN(',I4,')')                
      STOP                                                                      
C                                                                               
C5------PRINT NUMBER OF DRAINS IN THIS STRESS PERIOD.                           
  100 WRITE(IOUT,1) NDRAIN                                                      
    1 FORMAT(1H0,//1X,I5,' DRAINS')                                             
C                                                                               
C6------IF THERE ARE NO DRAINS THEN RETURN.                                     
      IF(NDRAIN.EQ.0) GO TO 260                                                 
C                                                                               
C7------READ AND PRINT DATA FOR EACH DRAIN.                                     
      WRITE(IOUT,3)                                                         
    3 FORMAT(1H0,15X,'LAYER',5X,'ROW',5X                                        
     1,'COL   ELEVATION   CONDUCTANCE   DRAIN NO.'/1X,15X,60('-'))              
      DO 250 II=1,NDRAIN                                                        
      READ (IN,4) K,I,J,DRAI(4,II),DRAI(5,II)                                   
    4 FORMAT(3I10,2F10.0)                                                       
      WRITE (IOUT,5) K,I,J,DRAI(4,II),DRAI(5,II),II                         
    5 FORMAT(1X,15X,I4,I9,I8,G13.4,G14.4,I8)                                    
      DRAI(1,II)=K                                                              
      DRAI(2,II)=I                                                              
      DRAI(3,II)=J                                                              
  250 CONTINUE                                                                  
C                                                                               
C8------RETURN                                                                  
  260 RETURN                                                                    
C                                                                               
      END                                                                       
      SUBROUTINE DRN1FM(NDRAIN,MXDRN,DRAI,HNEW,HCOF,RHS,IBOUND,                 
     1              NCOL,NROW,NLAY)                                             
C                                                                               
C-----VERSION 1030 10APR1985 DRN1FM                                             
C                                                                               
C     ******************************************************************        
C     ADD DRAIN FLOW TO SOURCE TERM                                             
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
      DOUBLE PRECISION HNEW                                                     
C                                                                               
      DIMENSION DRAI(5,MXDRN),HNEW(NCOL,NROW,NLAY),                             
     1         RHS(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY),                      
     1         HCOF(NCOL,NROW,NLAY)                                             
C     ------------------------------------------------------------------        
C                                                                               
C1------IF NDRAIN<=0 THERE ARE NO DRAINS. RETURN                                
      IF(NDRAIN.LE.0) RETURN                                                    
C                                                                               
C2------PROCESS EACH CELL IN THE DRAIN LIST                                     
      DO 100 L=1,NDRAIN                                                         
C                                                                               
C3------GET COLUMN, ROW AND LAYER OF CELL CONTAINING DRAIN.                     
      IL=DRAI(1,L)                                                              
      IR=DRAI(2,L)                                                              
      IC=DRAI(3,L)                                                              
C                                                                               
C4-------IF THE CELL IS EXTERNAL SKIP IT.                                       
      IF(IBOUND(IC,IR,IL).LE.0) GO TO 100                                       
C                                                                               
C5-------IF THE CELL IS INTERNAL GET THE DRAIN DATA.                            
      EL=DRAI(4,L)                                                              
      HHNEW=HNEW(IC,IR,IL)                                                      
C                                                                               
C6------IF HEAD IS LOWER THAN DRAIN THEN SKIP THIS CELL.                        
      IF(HHNEW.LE.EL) GO TO 100                                                 
C                                                                               
C7------HEAD IS HIGHER THAN DRAIN. ADD TERMS TO RHS AND HCOF.                   
      C=DRAI(5,L)                                                               
      HCOF(IC,IR,IL)=HCOF(IC,IR,IL)-C                                           
      RHS(IC,IR,IL)=RHS(IC,IR,IL)-C*EL                                          
  100 CONTINUE                                                                  
C                                                                               
C8------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE DRN1BD(NDRAIN,MXDRN,VBNM,VBVL,MSUM,DRAI,DELT,HNEW,             
     1        NCOL,NROW,NLAY,IBOUND,KSTP,KPER,IDRNCB,ICBCFL,BUFF,IOUT)          
C                                                                               
C-----VERSION 1301 28DEC1983 DRN1BD                                             
C                                                                               
C     ******************************************************************        
C     CALCULATE VOLUMETRIC BUDGET FOR DRAINS                                    
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------
      CHARACTER*4 VBNM,TEXT        
      DOUBLE PRECISION HNEW                                                     
C                                                                               
      DIMENSION VBNM(4,MSUM),VBVL(4,MSUM),DRAI(5,MXDRN),                        
     1           HNEW(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY),                   
     2           BUFF(NCOL,NROW,NLAY)                                           
      DIMENSION TEXT(4)                                                         
C                                                                               
      DATA TEXT(1),TEXT(2),TEXT(3),TEXT(4) /'    ','    ','  DR','AINS'/        
C     ------------------------------------------------------------------        
C                                                                               
C1------INITIALIZE CELL-BY-CELL FLOW TERM FLAG (IBD) AND                        
C1------ACCUMULATORS (RATIN AND RATOUT).                                        
      RATOUT=0.                                                                 
      IBD=0                                                                     
C                                                                               
C2------IF THERE ARE NO DRAINS THEN DO NOT ACCUMULATE DRAIN FLOW                
      IF(NDRAIN.LE.0) GO TO 200                                                 
C                                                                               
C3------TEST TO SEE IF CELL-BY-CELL FLOW TERMS ARE NEEDED.                      
      IF(ICBCFL.EQ.0 .OR. IDRNCB.LE.0) GO TO 60                                 
C                                                                               
C3B-----CELL-BY-CELL FLOW TERMS ARE NEEDED SET IBD AND CLEAR 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                                                                               
C4------FOR EACH DRAIN ACCUMULATE DRAIN FLOW                                    
   60 DO 100 L=1,NDRAIN                                                         
C                                                                               
C5------GET LAYER, ROW & COLUMN OF CELL CONTAINING REACH.                       
      IL=DRAI(1,L)                                                              
      IR=DRAI(2,L)                                                              
      IC=DRAI(3,L)                                                              
C                                                                               
C6------IF CELL IS EXTERNAL IGNORE IT.                                          
      IF(IBOUND(IC,IR,IL).LE.0) GO TO 100                                       
C                                                                               
C7------GET DRAIN PARAMETERS FROM DRAIN LIST.                                   
      EL=DRAI(4,L)                                                              
      C=DRAI(5,L)                                                               
      HHNEW=HNEW(IC,IR,IL)                                                      
C                                                                               
C8------IF HEAD LOWER THAN DRAIN THEN FORGET THIS CELL.                         
      IF(HHNEW.LE.EL) GO TO 100                                                 
C                                                                               
C9------HEAD HIGHER THAN DRAIN. CALCULATE Q=C*(EL-HHNEW) ADD Q TO RATOUT        
      Q=C*(EL-HHNEW)                                                            
      RATOUT=RATOUT-Q                                                           
C                                                                               
C10-----PRINT THE INDIVIDUAL RATES IF REQUESTED(IDRNCB<0).                      
      IF(IDRNCB.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,'   DRAIN',I4,                 
     1    '   LAYER',I3,'   ROW',I4,'   COL',I4,'   RATE',G15.7)                
C                                                                               
C11-----IF C-B-C FLOW TERMS ARE TO BE SAVED THEN ADD Q TO BUFFER.               
      IF(IBD.EQ.1) BUFF(IC,IR,IL)=BUFF(IC,IR,IL)+Q                              
  100 CONTINUE                                                                  
C                                                                               
C12-----IF C-B-C FLOW TERMS WILL BE SAVED CALL UBUDSV TO RECORD THEM.           
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,IDRNCB,BUFF,NCOL,NROW,            
     1                          NLAY,IOUT)                                      
C                                                                               
C13-----MOVE RATES,VOLUMES & LABELS INTO ARRAYS FOR PRINTING.                   
  200 VBVL(3,MSUM)=0.                                                           
      VBVL(4,MSUM)=RATOUT                                                       
      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                                                                               
C14-----INCREMENT BUDGET TERM COUNTER                                           
      MSUM=MSUM+1                                                               
C                                                                               
C15-----RETURN                                                                  
      RETURN                                                                    
      END                                                                       
C FILE EVT1.F77                                                                 
      SUBROUTINE EVT1AL(ISUM,LENX,LCIEVT,LCEVTR,LCEXDP,LCSURF,                  
     1                  NCOL,NROW,NEVTOP,IN,IOUT,IEVTCB)                        
C                                                                               
C-----VERSION 1741 24APR1985 EVT1AL                                             
C     ******************************************************************        
C     ALLOCATE ARRAY STORAGE FOR EVAPOTRANSPIRATION                             
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
C     ------------------------------------------------------------------        
C                                                                               
C1------IDENTIFY PACKAGE.                                                       
      WRITE(IOUT,1)IN                                                           
    1 FORMAT(1H0,'EVT1 -- EVAPOTRANSPIRATION PACKAGE, VERSION 1,',              
     1     ' 04/24/85',' INPUT READ FROM UNIT',I3)                              
C                                                                               
C2------READ NEVTOP AND IEVTCB.                                                 
      READ(IN,3)NEVTOP,IEVTCB                                                   
    3 FORMAT(2I10)                                                              
C                                                                               
C3------CHECK TO SEE THAT ET OPTION IS LEGAL.                                   
      IF(NEVTOP.GE.1.AND.NEVTOP.LE.2)GO TO 200                                  
C                                                                               
C3A-----IF ILLEGAL PRINT A MESSAGE & ABORT SIMULATION.                          
      WRITE(IOUT,8)                                                             
    8 FORMAT(1X,'ILLEGAL ET OPTION CODE. SIMULATION ABORTING')                  
      STOP                                                                      
C                                                                               
C4------IF THE OPTION IS LEGAL THEN PRINT THE OPTION CODE.                      
  200 IF(NEVTOP.EQ.1) WRITE(IOUT,201)                                           
  201 FORMAT(1X,'OPTION 1 -- EVAPOTRANSPIRATION FROM TOP LAYER')                
      IF(NEVTOP.EQ.2) WRITE(IOUT,202)                                           
  202 FORMAT(1X,'OPTION 2 -- EVAPOTRANSPIRATION FROM ONE SPECIFIED',            
     1      ' NODE IN EACH VERTICAL COLUMN')                                    
      IRK=ISUM                                                                  
C                                                                               
C5------IF CELL-BY-CELL TERMS TO BE SAVED THEN PRINT UNIT NUMBER.               
      IF(IEVTCB.GT.0) WRITE(IOUT,203) IEVTCB                                    
  203 FORMAT(1X,'CELL-BY-CELL FLOW TERMS WILL BE SAVED ON UNIT',I3)             
C                                                                               
C6------ALLOCATE SPACE FOR THE ARRAYS EVTR, EXDP AND SURF.                      
      LCEVTR=ISUM                                                               
      ISUM=ISUM+NCOL*NROW                                                       
      LCEXDP=ISUM                                                               
      ISUM=ISUM+NCOL*NROW                                                       
      LCSURF=ISUM                                                               
      ISUM=ISUM+NCOL*NROW                                                       
C                                                                               
C7------IF OPTION 2 THEN ALLOCATE SPACE FOR THE INDICATOR ARRAY(IEVT)           
      LCIEVT=ISUM                                                               
      IF(NEVTOP.NE.2)GO TO 300                                                  
      ISUM=ISUM+NCOL*NROW                                                       
C                                                                               
C8------CALCULATE & PRINT AMOUNT OF SPACE USED BY ET PACKAGE.                   
  300 IRK=ISUM-IRK                                                              
      WRITE(IOUT,4)IRK                                                          
    4 FORMAT(1X,I6,' ELEMENTS OF X ARRAY USED FOR EVAPOTRANSPIRATION')          
      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 EVT1RP(NEVTOP,IEVT,EVTR,EXDP,SURF,DELR,DELC,                   
     1                  NCOL,NROW,NLAY,IN,IOUT)                                 
C                                                                               
C-----VERSION 1631 08FEB1983 EVT1RP                                             
C     ******************************************************************        
C     READ EVAPOTRANSPIRATION DATA                                              
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------
      CHARACTER*4 ANAME        
      DIMENSION IEVT(NCOL,NROW),EVTR(NCOL,NROW),EXDP(NCOL,NROW),                
     1          SURF(NCOL,NROW),ANAME(6,4),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) /'    ','    ','  ET',' LAY','ER I','NDEX'/                  
      DATA ANAME(1,2),ANAME(2,2),ANAME(3,2),ANAME(4,2),ANAME(5,2),              
     1  ANAME(6,2) /'    ','    ','    ','  ET',' SUR','FACE'/                  
      DATA ANAME(1,3),ANAME(2,3),ANAME(3,3),ANAME(4,3),ANAME(5,3),              
     1  ANAME(6,3) /' EVA','POTR','ANSP','IRAT','ION ','RATE'/                  
      DATA ANAME(1,4),ANAME(2,4),ANAME(3,4),ANAME(4,4),ANAME(5,4),              
     1  ANAME(6,4) /'    ','    ','EXTI','NCTI','ON D','EPTH'/                  
C     ------------------------------------------------------------------        
C                                                                               
C1------READ FLAGS SHOWING WHETHER DATA IS TO BE REUSED.                        
      READ(IN,6)INSURF,INEVTR,INEXDP,INIEVT                                     
    6 FORMAT(4I10)                                                              
C                                                                               
C2------TEST INSURF TO SEE WHERE SURFACE ELEVATION COMES FROM.                  
      IF(INSURF.GE.0)GO TO 32                                                   
C                                                                               
C2A------IF INSURF<0 THEN REUSE SURFACE ARRAY FROM LAST STRESS PERIOD           
      WRITE(IOUT,3)                                                             
    3 FORMAT(1H0,'REUSING SURF FROM LAST STRESS PERIOD')                        
      GO TO 35                                                                  
C                                                                               
C3-------IF INSURF=>0 THEN CALL MODULE U2DREL TO READ SURFACE.                  
   32 CALL U2DREL(SURF,ANAME(1,2),NROW,NCOL,0,IN,IOUT)                          
C                                                                               
C4------TEST INEVTR TO SEE WHERE MAX ET RATE COMES FROM.                        
   35 IF(INEVTR.GE.0)GO TO 37                                                   
C                                                                               
C4A-----IF INEVTR<0 THEN REUSE MAX ET RATE.                                     
      WRITE(IOUT,4)                                                             
    4 FORMAT(1H0,'REUSING EVTR FROM LAST STRESS PERIOD')                        
      GO TO 45                                                                  
C                                                                               
C5------IF INEVTR=>0 CALL MODULE U2DREL TO READ MAX ET RATE.                    
   37 CALL U2DREL(EVTR,ANAME(1,3),NROW,NCOL,0,IN,IOUT)                          
C                                                                               
C6------MULTIPLY MAX ET RATE BY CELL AREA TO GET VOLUMETRIC RATE                
      DO 40 IR=1,NROW                                                           
      DO 40 IC=1,NCOL                                                           
      EVTR(IC,IR)=EVTR(IC,IR)*DELR(IC)*DELC(IR)                                 
   40 CONTINUE                                                                  
C                                                                               
C7------TEST INEXDP TO SEE WHERE EXTINCTION DEPTH COMES FROM                    
   45 IF(INEXDP.GE.0)GO TO 47                                                   
C                                                                               
C7A------IF INEXDP<0 REUSE EXTINCTION DEPTH FROM LAST STRESS PERIOD             
      WRITE(IOUT,5)                                                             
    5 FORMAT(1H0,'REUSING EXDP FROM LAST STRESS PERIOD')                        
      GO TO 48                                                                  
C                                                                               
C8-------IF INEXDP=>0 CALL MODULE U2DREL TO READ EXTINCTION DEPTH               
   47 CALL U2DREL(EXDP,ANAME(1,4),NROW,NCOL,0,IN,IOUT)                          
C                                                                               
C9------IF OPTION(NEVTOP) IS 2 THEN WE NEED AN INDICATOR ARRAY.                 
  48  IF(NEVTOP.NE.2)GO TO 50                                                   
C                                                                               
C10------IF INIEVT<0 THEN REUSE LAYER INDICATOR ARRAY.                          
      IF(INIEVT.GE.0)GO TO 49                                                   
      WRITE(IOUT,2)                                                             
    2 FORMAT(1H0,'REUSING IEVT FROM LAST STRESS PERIOD')                        
      GO TO 50                                                                  
C                                                                               
C11------IF INIEVT=>0 THEN CALL MODULE U2DINT TO READ INDICATOR ARRAY.          
   49 CALL U2DINT(IEVT,ANAME(1,1),NROW,NCOL,0,IN,IOUT)                          
C                                                                               
C12-----RETURN                                                                  
   50 RETURN                                                                    
      END                                                                       
      SUBROUTINE EVT1FM(NEVTOP,IEVT,EVTR,EXDP,SURF,RHS,HCOF,                    
     1                  IBOUND,HNEW,NCOL,NROW,NLAY)                             
C                                                                               
C-----VERSION 1031 10APR1985 EVT1FM                                             
C     ******************************************************************        
C        ADD EVAPOTRANSPIRATION TO RHS AND HCOF                                 
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
      DOUBLE PRECISION HNEW                                                     
      DIMENSION IEVT(NCOL,NROW),EVTR(NCOL,NROW),EXDP(NCOL,NROW),                
     1          SURF(NCOL,NROW),RHS(NCOL,NROW,NLAY),                            
     2          HCOF(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY),                    
     3          HNEW(NCOL,NROW,NLAY)                                            
C     ------------------------------------------------------------------        
C                                                                               
C1------PROCESS EACH HORIZONTAL CELL LOCATION                                   
      DO 10 IR=1,NROW                                                           
      DO 10 IC=1,NCOL                                                           
C                                                                               
C2------SET THE LAYER INDEX EQUAL TO 1      .                                   
      IL=1                                                                      
C                                                                               
C3------IF OPTION 2 IS SPECIFIED THEN GET LAYER INDEX FROM IEVT ARRAY           
      IF(NEVTOP.EQ.2)IL=IEVT(IC,IR)                                             
C                                                                               
C4------IF THE CELL IS EXTERNAL IGNORE IT.                                      
      IF(IBOUND(IC,IR,IL).LE.0)GO TO 10                                         
      C=EVTR(IC,IR)                                                             
      S=SURF(IC,IR)                                                             
      H=HNEW(IC,IR,IL)                                                          
C                                                                               
C5------IF AQUIFER HEAD IS GREATER THAN OR EQUAL TO SURF, ET IS CONSTANT        
      IF(H.LT.S) GO TO 5                                                        
C                                                                               
C5A-----SUBTRACT -EVTR FROM RHS                                                 
      RHS(IC,IR,IL)=RHS(IC,IR,IL) + C                                           
      GO TO 10                                                                  
C                                                                               
C6------IF DEPTH TO WATER>=EXTINCTION DEPTH THEN ET IS 0                        
    5 D=S-H                                                                     
      X=EXDP(IC,IR)                                                             
      IF(D.GE.X)GO TO 10                                                        
C                                                                               
C7------LINEAR RANGE. ADD ET TERMS TO BOTH RHS AND HCOF.                        
      RHS(IC,IR,IL)=RHS(IC,IR,IL)+C-C*S/X                                       
      HCOF(IC,IR,IL)=HCOF(IC,IR,IL)-C/X                                         
   10 CONTINUE                                                                  
C                                                                               
C8------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE EVT1BD(NEVTOP,IEVT,EVTR,EXDP,SURF,IBOUND,HNEW,                 
     1           NCOL,NROW,NLAY,DELT,VBVL,VBNM,MSUM,KSTP,KPER,                  
     2           IEVTCB,ICBCFL,BUFF,IOUT)                                       
C-----VERSION 1405 10FEB1983 EVT1BD                                             
C     ******************************************************************        
C     CALCULATE VOLUMETRIC BUDGET FOR EVAPOTRANSPIRATION                        
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------
      CHARACTER*4 VBNM,TEXT        
      DOUBLE PRECISION HNEW                                                     
      DIMENSION IEVT(NCOL,NROW),EVTR(NCOL,NROW),EXDP(NCOL,NROW),                
     1          SURF(NCOL,NROW),IBOUND(NCOL,NROW,NLAY),                         
     2          VBVL(4,20),VBNM(4,20),HNEW(NCOL,NROW,NLAY),                     
     3          BUFF(NCOL,NROW,NLAY)                                            
      DIMENSION TEXT(4)                                                         
      DATA TEXT(1),TEXT(2),TEXT(3),TEXT(4) /'    ','    ','    ','  ET'/        
C     ------------------------------------------------------------------        
C                                                                               
C1------CLEAR THE RATE ACCUMULATOR.                                             
      RATOUT=0                                                                  
C                                                                               
C2------IF CELL-BY-CELL FLOW TERMS WILL BE SAVED THEN CLEAR THE BUFFER.         
      IBD=0                                                                     
      IF(IEVTCB.LE.0 .OR. ICBCFL.EQ.0) GO TO 5                                  
      IBD=1                                                                     
      DO 4 IL=1,NLAY                                                            
      DO 4 IR=1,NROW                                                            
      DO 4 IC=1,NCOL                                                            
      BUFF(IC,IR,IL)=0.                                                         
    4 CONTINUE                                                                  
C                                                                               
C3------PROCESS EACH HORIZONTAL CELL LOCATION                                   
    5 DO 10 IR=1,NROW                                                           
      DO 10 IC=1,NCOL                                                           
C                                                                               
C4------SET THE LAYER INDEX EQUAL TO 1      .                                   
      IL=1                                                                      
C                                                                               
C5------IF OPTION 2 IS SPECIFIED THEN GET LAYER INDEX FROM IEVT ARRAY           
      IF(NEVTOP.EQ.2)IL=IEVT(IC,IR)                                             
C                                                                               
C6------IF CELL IS EXTERNAL THEN IGNORE IT.                                     
      IF(IBOUND(IC,IR,IL).LE.0)GO TO 10                                         
      C=EVTR(IC,IR)                                                             
      S=SURF(IC,IR)                                                             
      H=HNEW(IC,IR,IL)                                                          
C                                                                               
C7------IF AQUIFER HEAD => SURF,SET Q=MAX ET RATE                               
      IF(H.LT.S) GO TO 7                                                        
      Q=-C                                                                      
      GO TO 9                                                                   
C                                                                               
C8------IF DEPTH=>EXTINCTION DEPTH, ET IS 0                                     
    7 X=EXDP(IC,IR)                                                             
      D=S-H                                                                     
      IF(D.GE.X)GO TO 10                                                        
C                                                                               
C9------LINEAR RANGE . Q=-EVTR(H-EXEL)/EXDP                                     
      Q=C*D/X-C                                                                 
C                                                                               
C10-----ACCUMULATE TOTAL FLOW RATE                                              
    9 RATOUT=RATOUT-Q                                                           
C                                                                               
C11-----IF CELL-BY-CELL FLOW TERMS TO BE SAVED THE ADD Q TO BUFFER.             
      IF(IBD.EQ.1) BUFF(IC,IR,IL)=Q                                             
   10 CONTINUE                                                                  
C                                                                               
C12-----IF C-B-C TO BE SAVED CALL MODULE UBUDSV TO RECORD THEM.                 
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,IEVTCB,BUFF,NCOL,NROW,            
     1                          NLAY,IOUT)                                      
C                                                                               
C13-----MOVE TOTAL ET RATE INTO VBVL FOR PRINTING BY BAS1OT.                    
      VBVL(3,MSUM)=0.                                                           
      VBVL(4,MSUM)=RATOUT                                                       
C                                                                               
C14-----ADD ET(ET_RATE TIMES STEP LENGTH) TO VBVL                               
      VBVL(1,MSUM)=0.                                                           
      VBVL(2,MSUM)=VBVL(2,MSUM)+RATOUT*DELT                                     
C                                                                               
C15-----MOVE BUDGET TERM LABELS TO VBNM FOR PRINT BY MODULE BAS1OT              
      VBNM(1,MSUM)=TEXT(1)                                                      
      VBNM(2,MSUM)=TEXT(2)                                                      
      VBNM(3,MSUM)=TEXT(3)                                                      
      VBNM(4,MSUM)=TEXT(4)                                                      
C                                                                               
C16-----INCREMENT BUDGET TERM COUNTER                                           
      MSUM=MSUM+1                                                               
C                                                                               
C17-----RETURN                                                                  
      RETURN                                                                    
      END                                                                       
CWES INSERT EVT2.F77 HERE                                                       
C FILE GHB1.F77                                                                 
      SUBROUTINE GHB1AL(ISUM,LENX,LCBNDS,NBOUND,MXBND,IN,IOUT,                  
     1                       IGHBCB)                                            
C                                                                               
C-----VERSION 1742 24APR1985 GHB1AL                                             
C     ******************************************************************        
C     ALLOCATE ARRAY STORAGE FOR HEAD-DEPENDENT BOUNDARIES                      
C     ******************************************************************        
C                                                                               
C     SPECIFICATIONS:                                                           
C     ------------------------------------------------------------------        
C     ------------------------------------------------------------------        
C                                                                               
C1------IDENTIFY PACKAGE AND INITIALIZE # OF GENERAL HEAD BOUNDS                
      WRITE(IOUT,1)IN                                                           
    1 FORMAT(1H0,'GHB1 -- GHB PACKAGE, VERSION 1, 04/24/85',                    
     2' INPUT READ FROM UNIT',I3)                                               
      NBOUND=0                                                                  
C                                                                               
C2------READ AND PRINT MXBND AND IGHBCB (MAX # OF BOUNDS AND UNIT               
C2------FOR CELL-BY-CELL FLOW TERMS FOR GHB)                                    
      READ(IN,2) MXBND,IGHBCB                                                   
    2 FORMAT(2I10)                                                              
      WRITE(IOUT,3) MXBND                                                       
    3 FORMAT(1H ,'MAXIMUM OF',I5,' HEAD-DEPENDENT BOUNDARY NODES')              
      IF(IGHBCB.GT.0) WRITE(IOUT,9) IGHBCB                                      
    9 FORMAT(1X,'CELL-BY-CELL FLOW WILL BE RECORDED ON UNIT',I3)                
      IF(IGHBCB.LT.0) WRITE(IOUT,8)                                             
    8 FORMAT(1X,'CELL-BY-CELL FLOW WILL BE PRINTED WHEN ICBCFL NOT 0')          
C                                                                               
C3------SET LCBNDS EQUAL TO ADDRESS OF FIRST UNUSED SPACE IN X.                 
      LCBNDS=ISUM                                                               
C                                                                               
C4------CALCULATE AMOUNT OF SPACE USED BY THE GENERAL HEAD LIST.                
      ISP=5*MXBND                                                               
      ISUM=ISUM+ISP                                                             
C                                                                               
C5------PRINT AMOUNT OF SPACE USED BY THE GHB PACKAGE                           
      WRITE(IOUT,4) ISP                                                         
    4 FORMAT(1X,I6,' ELEMENTS IN X ARRAY ARE USED FOR HEAD',                    
     1      '-DEPENDENT BOUNDARIES')                                            
      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                                                                               
C6------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE GHB1RP(BNDS,NBOUND,MXBND,IN,IOUT)                              
C                                                                               
C                                                                               
C-----VERSION 1651 02FEB1983 GHB1RP                                             
C     ******************************************************************        
C     READ DATA FOR GHB                                                         
C     ******************************************************************        
C                                                                               
C     SPECIFICATIONS:                                                           
C     ------------------------------------------------------------------        
      DIMENSION BNDS(5,MXBND)                                                   
C     ------------------------------------------------------------------        
C                                                                               
C1------READ ITMP(# OF GENERAL HEAD BOUNDS 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 HEAD-DEPENDENT BOUNDS FROM LAST STRESS',              
     1      ' PERIOD')                                                          
      GO TO 260                                                                 
C                                                                               
C3------IF ITMP=>0 THEN IT IS THE # OF GENERAL HEAD BOUNDS.                     
   50 NBOUND=ITMP                                                               
C                                                                               
C4------IF MAX NUMBER OF BOUNDS IS EXCEEDED THEN STOP                           
      IF(NBOUND.LE.MXBND) GO TO 100                                             
      WRITE(IOUT,99) NBOUND,MXBND                                               
   99 FORMAT(1H0,'NBOUND(',I4,') IS GREATER THAN MXBND(',I4,')')                
C                                                                               
C4A-----ABNORMAL STOP                                                           
      STOP                                                                      
C                                                                               
C5------PRINT # OF GENERAL HEAD BOUNDS THIS STRESS PERIOD                       
  100 WRITE(IOUT,1) NBOUND                                                      
    1 FORMAT(1H0,//1X,I5,' HEAD-DEPENDENT BOUNDARY NODES')                      
C                                                                               
C6------IF THERE ARE NO GENERAL HEAD BOUNDS THEN RETURN.                        
      IF(NBOUND.EQ.0) GO TO 260                                                 
C                                                                               
C7------READ & PRINT DATA FOR EACH GENERAL HEAD BOUNDARY.                       
      WRITE(IOUT,3)                                                         
    3 FORMAT(1H0,15X,'LAYER',5X,'ROW',5X                                        
     1,'COL   ELEVATION   CONDUCTANCE   BOUND NO.'/1X,15X,60('-'))              
      DO 250 II=1,NBOUND                                                        
      READ (IN,4) K,I,J,BNDS(4,II),BNDS(5,II)                                   
    4 FORMAT(3I10,2F10.0)                                                       
      WRITE (IOUT,5) K,I,J,BNDS(4,II),BNDS(5,II),II                         
    5 FORMAT(1X,15X,I4,I9,I8,G13.4,G14.4,I8)                                    
      BNDS(1,II)=K                                                              
      BNDS(2,II)=I                                                              
      BNDS(3,II)=J                                                              
  250 CONTINUE                                                                  
C                                                                               
C8------RETURN                                                                  
  260 RETURN                                                                    
      END                                                                       
      SUBROUTINE GHB1FM(NBOUND,MXBND,BNDS,HCOF,RHS,IBOUND,                      
     1              NCOL,NROW,NLAY)                                             
C                                                                               
C-----VERSION 1037 10APR1985 GHB1FM                                             
C     ******************************************************************        
C     ADD GHB TERMS TO RHS AND HCOF                                             
C     ******************************************************************        
C                                                                               
C     SPECIFICATIONS:                                                           
C     ------------------------------------------------------------------        
      DIMENSION BNDS(5,MXBND),HCOF(NCOL,NROW,NLAY),                             
     1         RHS(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY)                       
C     ------------------------------------------------------------------        
C                                                                               
C1------IF NBOUND<=0 THEN THERE ARE NO GENERAL HEAD BOUNDS. RETURN.             
      IF(NBOUND.LE.0) RETURN                                                    
C                                                                               
C2------PROCESS EACH ENTRY IN THE GENERAL HEAD BOUND LIST (BNDS)                
      DO 100 L=1,NBOUND                                                         
C                                                                               
C3------GET COLUMN, ROW AND LAYER OF CELL CONTAINING BOUNDARY                   
      IL=BNDS(1,L)                                                              
      IR=BNDS(2,L)                                                              
      IC=BNDS(3,L)                                                              
C                                                                               
C4------IF THE CELL IS EXTERNAL THEN SKIP IT.                                   
      IF(IBOUND(IC,IR,IL).LE.0) GO TO 100                                       
C                                                                               
C5------SINCE THE CELL IS INTERNAL GET THE BOUNDARY DATA.                       
      HB=BNDS(4,L)                                                              
      C=BNDS(5,L)                                                               
C                                                                               
C6------ADD TERMS TO RHS AND HCOF                                               
      HCOF(IC,IR,IL)=HCOF(IC,IR,IL)-C                                           
      RHS(IC,IR,IL)=RHS(IC,IR,IL)-C*HB                                          
  100 CONTINUE                                                                  
C                                                                               
C7------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE GHB1BD(NBOUND,MXBND,VBNM,VBVL,MSUM,BNDS,DELT,HNEW,             
     1   NCOL,NROW,NLAY,IBOUND,KSTP,KPER,IGHBCB,ICBCFL,BUFF,IOUT)               
C                                                                               
C-----VERSION 1152 20MAY1983 GHB1BD                                             
C     ******************************************************************        
C     CALCULATE VOLUMETRIC BUDGET FOR GHB                                       
C     ******************************************************************        
C                                                                               
C     SPECIFICATIONS:                                                           
C     ------------------------------------------------------------------
      CHARACTER*4 VBNM,TEXT        
      DOUBLE PRECISION HNEW                                                     
      DIMENSION VBNM(4,MSUM),VBVL(4,MSUM),BNDS(5,MXBND),                        
     1           HNEW(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY),                   
     2           BUFF(NCOL,NROW,NLAY)                                           
      DIMENSION TEXT(4)                                                         
      DATA TEXT(1),TEXT(2),TEXT(3),TEXT(4) /' HEA','D DE','P BO','UNDS'/        
C     ------------------------------------------------------------------        
C                                                                               
C1------INITIALIZE CELL-BY-CELL FLOW TERM FLAG (IBD) AND                        
C1------ACCUMULATORS (RATIN AND RATOUT)                                         
      IBD=0                                                                     
      RATOUT=0.                                                                 
      RATIN=0.                                                                  
C                                                                               
C2------IF NO BOUNDARIES THEN KEEP ZEROES IN ACCUMULATORS.                      
      IF(NBOUND.EQ.0) GO TO 200                                                 
C                                                                               
C3------TEST TO SEE IF CELL-BY-CELL FLOW TERMS ARE NEEDED.                      
      IF(ICBCFL.EQ.0 .OR. IGHBCB.LE.0) GO TO 10                                 
C                                                                               
C3A-----SINCE CELL-BY-CELL FLOW TERMS ARE NEEDED CLEAR BUFFER & SET             
C3A-----THE FLAG IBD.                                                           
      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 GENERAL HEAD BOUND ACCUMULATE FLOW INTO AQUIFER                
   10 DO 100 L=1,NBOUND                                                         
C                                                                               
C5------GET LAYER, ROW AND COLUMN OF EACH GENERAL HEAD BOUNDARY.                
      IL=BNDS(1,L)                                                              
      IR=BNDS(2,L)                                                              
      IC=BNDS(3,L)                                                              
C                                                                               
C6------IF CELL IS EXTERNAL THEN IGNORE IT.                                     
      IF(IBOUND(IC,IR,IL).LE.0) GO TO 100                                       
C                                                                               
C7------GET PARAMETERS FROM BOUNDARY LIST.                                      
      HHNEW=HNEW(IC,IR,IL)                                                      
      HB=BNDS(4,L)                                                              
      C=BNDS(5,L)                                                               
C                                                                               
C8------CALCULATE THE FOW RATE INTO THE CELL                                    
      RATE=C*(HB-HHNEW)                                                         
C                                                                               
C9------PRINT THE INDIVIDUAL RATES IF REQUESTED(IGHBCB<0).                      
      IF(IGHBCB.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,'   BOUNDARY',I4,              
     1    '   LAYER',I3,'   ROW',I4,'   COL',I4,'   RATE',G15.7)                
C                                                                               
C10-----IF CELL-BY-CELL TERMS ARE TO BE SAVED THEN PUT RATE IN BUFFER           
      IF(IBD.EQ.1) BUFF(IC,IR,IL)=BUFF(IC,IR,IL)+RATE                           
C                                                                               
C11-----SEE IF FLOW IS INTO AQUIFER OR OUT OF AQUIFER.                          
      IF(RATE)94,100,96                                                         
C                                                                               
C12------FLOW IS OUT OF AQUIFER SUBTRACT RATE FROM RATOUT                       
   94 RATOUT=RATOUT-RATE                                                        
      GO TO 100                                                                 
C                                                                               
C13-----FLOW IS INTO AQIFER ADD RATE TO RATIN                                   
   96 RATIN=RATIN+RATE                                                          
  100 CONTINUE                                                                  
C                                                                               
C14-----IF CELL-BY-CELL TERMS ARE TO BE SAVED THEN CALL                         
C14-----UTILITY MODULE UBUDSV                                                   
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,IGHBCB,BUFF,NCOL,NROW,            
     1                          NLAY,IOUT)                                      
C                                                                               
C15-----MOVE RATES, VOLUMES AND LABELS INTO ARRAYS FOR PRINTING                 
  200 VBVL(3,MSUM)=RATIN                                                        
      VBVL(1,MSUM)=VBVL(1,MSUM)+RATIN*DELT                                      
      VBVL(4,MSUM)=RATOUT                                                       
      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                                                                               
C16-----INCREMENT THE BUDGET TERM COUNTER                                       
      MSUM=MSUM+1                                                               
C                                                                               
C17-----RETURN                                                                  
      RETURN                                                                    
      END                                                                       
