C FILE BCF1.F77                                                                 
      SUBROUTINE BCF1AL(ISUM,LENX,LCSC1,LCHY,LCBOT,                             
     1     LCTOP,LCSC2,LCTRPY,IN,ISS,NCOL,NROW,NLAY,IOUT,IBCFCB)                
C                                                                               
C-----VERSION 1738 24APR1985 BCF1AL                                             
C                                                                               
C     ******************************************************************        
C     ALLOCATE ARRAY STORAGE FOR BLOCK-CENTERED FLOW PACKAGE                    
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
      COMMON /FLWCOM/LAYCON(80)                                                 
C     ------------------------------------------------------------------        
C                                                                               
C1------IDENTIFY PACKAGE                                                        
      WRITE(IOUT,1)IN                                                           
    1 FORMAT(1H0,'BCF1 -- BLOCK-CENTERED FLOW PACKAGE, VERSION 1',              
     1', 04/24/85',' INPUT READ FROM UNIT',I3)                                  
C                                                                               
C2------READ AND PRINT ISS (STEADY-STATE FLAG) AND IBCFCB (FLAG FOR             
C2------PRINTING OR UNIT# FOR RECORDING CELL-BY-CELL FLOW TERMS)                
      READ(IN,2) ISS,IBCFCB                                                     
    2 FORMAT(2I10)                                                              
      IF(ISS.EQ.0) WRITE(IOUT,3)                                                
    3 FORMAT(1X,'TRANSIENT SIMULATION')                                         
      IF(ISS.NE.0) WRITE(IOUT,4)                                                
    4 FORMAT(1X,'STEADY-STATE SIMULATION')                                      
      IF(IBCFCB.GT.0) WRITE(IOUT,9) IBCFCB                                      
    9 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE RECORDED ON UNIT',I3)               
      IF(IBCFCB.LT.0) WRITE(IOUT,88)                                            
   88 FORMAT(1X,'CONSTANT HEAD CELL-BY-CELL FLOWS WILL BE PRINTED')             
C                                                                               
C3------READ TYPE CODE FOR EACH LAYER AND COUNT TOPS AND BOTTOMS                
      IF(NLAY.LE.80) GO TO 50                                                   
      WRITE(IOUT,11)                                                            
   11 FORMAT(1H0,'YOU HAVE SPECIFIED MORE THAN 80 MODEL LAYERS'/1X,             
     1  'SPACE IS RESERVED FOR A MAXIMUM OF 80 LAYERS IN ARRAY LAYCON')         
      STOP                                                                      
C                                                                               
C3A-----READ LAYER TYPE CODES.                                                  
   50 READ(IN,51) (LAYCON(I),I=1,NLAY)                                          
   51 FORMAT(40I2)                                                              
C        BOTTOM IS READ FOR TYPES 1,3    TOP IS READ FOR TYPES 2,3              
      WRITE(IOUT,52)                                                            
   52 FORMAT(1X,5X,'LAYER  AQUIFER TYPE',/1X,5X,19('-'))                        
C                                                                               
C3B-----INITIALIZE TOP AND BOTTOM COUNTERS.                                     
      NBOT=0                                                                    
      NTOP=0                                                                    
C                                                                               
C3C------PRINT LAYER TYPE AND COUNT TOPS AND BOTTOMS NEEDED.                    
      DO 100 I=1,NLAY                                                           
C                                                                               
C3C1----PRINT LAYER NUMBER AND LAYER TYPE CODE.                                 
      L=LAYCON(I)                                                               
      WRITE(IOUT,7) I,L                                                         
    7 FORMAT(1X,I9,I10)                                                         
C                                                                               
C3C2----ONLY THE TOP LAYER CAN BE UNCONFINED(LAYCON=1).                         
      IF(L.NE.1 .OR. I.EQ.1) GO TO 70                                           
      WRITE(IOUT,8)                                                             
    8 FORMAT(1H0,'AQUIFER TYPE 1 IS ONLY ALLOWED IN TOP LAYER')                 
      STOP                                                                      
C                                                                               
C3C3----LAYER TYPES 1 AND 3 NEED A BOTTOM. ADD 1 TO KB.                         
   70 IF(L.EQ.1 .OR. L.EQ.3) NBOT=NBOT+1                                        
C                                                                               
C3C4----LAYER TYPES 2 AND 3 NEED A TOP. ADD 1 TO KT.                            
      IF(L.EQ.2 .OR. L.EQ.3) NTOP=NTOP+1                                        
  100 CONTINUE                                                                  
C                                                                               
C                                                                               
C                                                                               
C4------COMPUTE DIMENSIONS FOR ARRAYS.                                          
      NRC=NROW*NCOL                                                             
      ISIZ=NRC*NLAY                                                             
C                                                                               
C5------ALLOCATE SPACE FOR ARRAYS. IF RUN IS TRANSIENT(ISS=0)                   
C5------THEN SPACE MUST BE ALLOCATED FOR STORAGE.                               
      ISOLD=ISUM                                                                
      LCSC1=ISUM                                                                
      IF(ISS.EQ.0) ISUM=ISUM+ISIZ                                               
      LCSC2=ISUM                                                                
      IF(ISS.EQ.0) ISUM=ISUM+NRC*NTOP                                           
      LCTRPY=ISUM                                                               
      ISUM=ISUM+NLAY                                                            
      LCBOT=ISUM                                                                
      ISUM=ISUM+NRC*NBOT                                                        
      LCHY=ISUM                                                                 
      ISUM=ISUM+NRC*NBOT                                                        
      LCTOP=ISUM                                                                
      ISUM=ISUM+NRC*NTOP                                                        
C                                                                               
C6------PRINT THE AMOUNT OF SPACE USED BY THE BCF PACKAGE.                      
      ISP=ISUM-ISOLD                                                            
      WRITE(IOUT,101) ISP                                                       
  101 FORMAT(1X,I6,' ELEMENTS IN X ARRAY ARE USED BY BCF')                      
      ISUM1=ISUM-1                                                              
      WRITE(IOUT,102) ISUM1,LENX                                                
  102 FORMAT(1X,I6,' ELEMENTS OF X ARRAY USED OUT OF',I7)                       
      IF(ISUM1.GT.LENX) WRITE(IOUT,103)                                         
  103 FORMAT(1X,'   ***X ARRAY MUST BE DIMENSIONED LARGER***')                  
C                                                                               
C7------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE BCF1RP(IBOUND,HNEW,SC1,HY,CR,CC,CV,DELR,DELC,                  
     1        BOT,TOP,SC2,TRPY,IN,ISS,NCOL,NROW,NLAY,NODES,IOUT)                
C                                                                               
C-----VERSION 1003 03MAY1983 BCF1RP                                             
C                                                                               
C     ******************************************************************        
C     READ AND INITIALIZE DATA FOR BLOCK-CENTERED FLOW PACKAGE                  
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------
      CHARACTER*4 ANAME        
      DOUBLE PRECISION HNEW                                                     
C                                                                               
      DIMENSION HNEW(NODES),SC1(NODES),HY(NODES),CR(NODES),CC(NODES),           
     1          CV(NODES),ANAME(6,10),DELR(NCOL),DELC(NROW),BOT(NODES),         
     1          TOP(NODES),SC2(NODES),TRPY(NLAY),IBOUND(NODES)                  
C                                                                               
      COMMON /FLWCOM/LAYCON(80)                                                 
C                                                                               
      DATA ANAME(1,1),ANAME(2,1),ANAME(3,1),ANAME(4,1),ANAME(5,1),              
     1  ANAME(6,1) /'    ','PRIM','ARY ','STOR','AGE ','COEF'/                 
      DATA ANAME(1,2),ANAME(2,2),ANAME(3,2),ANAME(4,2),ANAME(5,2),              
     1  ANAME(6,2) /'    ','TRAN','SMIS','. AL','ONG ','ROWS'/                 
      DATA ANAME(1,3),ANAME(2,3),ANAME(3,3),ANAME(4,3),ANAME(5,3),              
     1  ANAME(6,3) /'   H','YD. ','COND','. AL','ONG ','ROWS'/                 
      DATA ANAME(1,4),ANAME(2,4),ANAME(3,4),ANAME(4,4),ANAME(5,4),              
     1  ANAME(6,4) /'VERT',' HYD',' CON','D /T','HICK','NESS'/                 
      DATA ANAME(1,5),ANAME(2,5),ANAME(3,5),ANAME(4,5),ANAME(5,5),              
     1  ANAME(6,5) /'    ','    ','    ','    ','  BO','TTOM'/                  
      DATA ANAME(1,6),ANAME(2,6),ANAME(3,6),ANAME(4,6),ANAME(5,6),              
     1  ANAME(6,6) /'    ','    ','    ','    ','    ',' TOP'/                  
      DATA ANAME(1,7),ANAME(2,7),ANAME(3,7),ANAME(4,7),ANAME(5,7),              
     1  ANAME(6,7) /'  SE','COND','ARY ','STOR','AGE ','COEF'/                  
      DATA ANAME(1,8),ANAME(2,8),ANAME(3,8),ANAME(4,8),ANAME(5,8),              
     1  ANAME(6,8) /'COLU','MN T','O RO','W AN','ISOT','ROPY'/                  
      DATA ANAME(1,9),ANAME(2,9),ANAME(3,9),ANAME(4,9),ANAME(5,9),              
     1  ANAME(6,9) /'    ','    ','    ','    ','    ','DELR'/                  
      DATA ANAME(1,10),ANAME(2,10),ANAME(3,10),ANAME(4,10),ANAME(5,10),         
     1  ANAME(6,10) /'    ','    ','    ','    ','    ','DELC'/                 
C     ------------------------------------------------------------------        
C                                                                               
C1------CALCULATE NUMBER OF NODES IN A LAYER AND READ TRPY,DELR,DELC            
      NIJ=NCOL*NROW                                                             
C                                                                               
      CALL U1DREL(TRPY,ANAME(1,8),NLAY,IN,IOUT)                                 
      CALL U1DREL(DELR,ANAME(1,9),NCOL,IN,IOUT)                                 
      CALL U1DREL(DELC,ANAME(1,10),NROW,IN,IOUT)                                
C                                                                               
C2------READ ALL PARAMETERS FOR EACH LAYER                                      
      KT=0                                                                      
      KB=0                                                                      
      DO 200 K=1,NLAY                                                           
C                                                                               
C2A-----FIND ADDRESS OF EACH LAYER IN THREE DIMENSION ARRAYS.                   
      IF(LAYCON(K).EQ.1 .OR. LAYCON(K).EQ.3) KB=KB+1                            
      IF(LAYCON(K).EQ.2 .OR. LAYCON(K).EQ.3) KT=KT+1                            
      LOC=1+(K-1)*NIJ                                                           
      LOCB=1+(KB-1)*NIJ                                                         
      LOCT=1+(KT-1)*NIJ                                                         
C                                                                               
C2B-----READ PRIMARY STORAGE COEFFICIENT INTO ARRAY SC1 IF TRANSIENT            
      IF(ISS.EQ.0)CALL U2DREL(SC1(LOC),ANAME(1,1),NROW,NCOL,KK,IN,IOUT)   
C                                                                               
C2C-----READ TRANSMISSIVITY INTO ARRAY CC IF LAYER TYPE IS 0 OR 2               
      IF(LAYCON(K).EQ.3 .OR. LAYCON(K).EQ.1) GO TO 100                          
      CALL U2DREL(CC(LOC),ANAME(1,2),NROW,NCOL,KK,IN,IOUT)               
      GO TO 110                                                                 
C                                                                               
C2D-----READ HYDRAULIC CONDUCTIVITY(HY) AND BOTTOM ELEVATION(BOT)               
C2D-----IF LAYER TYPE IS 1 OR 3                                                 
  100 CALL U2DREL(HY(LOCB),ANAME(1,3),NROW,NCOL,KK,IN,IOUT)               
      CALL U2DREL(BOT(LOCB),ANAME(1,5),NROW,NCOL,KK,IN,IOUT)              
C                                                                               
C2E-----READ VERTICAL HYCOND/THICK INTO ARRAY CV IF NOT BOTTOM LAYER            
C2E----- READ AS HYCOND/THICKNESS -- CONVERTED TO CONDUCTANCE LATER             
  110 IF(K.EQ.NLAY) GO TO 120                                                   
      CALL U2DREL(CV(LOC),ANAME(1,4),NROW,NCOL,KK,IN,IOUT)               
C                                                                               
C2F-----READ SECONDARY STORAGE COEFFICIENT INTO ARRAY SC2 IF TRANSIENT          
C2F-----AND LAYER TYPE IS 2 OR 3                                                
  120 IF(LAYCON(K).NE.3 .AND. LAYCON(K).NE.2) GO TO 200                         
      IF(ISS.EQ.0)CALL U2DREL(SC2(LOCT),ANAME(1,7),NROW,NCOL,KK,IN,IOUT)       
C                                                                               
C2G-----READ TOP ELEVATION(TOP) IF LAYER TYPE IS 2 OR 3                         
      CALL U2DREL(TOP(LOCT),ANAME(1,6),NROW,NCOL,KK,IN,IOUT)              
  200 CONTINUE                                                                  
C                                                                               
C3------PREPARE AND CHECK BCF DATA                                              
      CALL SBCF1N(HNEW,IBOUND,SC1,SC2,CR,CC,CV,HY,TRPY,DELR,DELC,ISS,           
     1         NCOL,NROW,NLAY,IOUT)                                             
C                                                                               
C4------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE BCF1FM(HCOF,RHS,HOLD,SC1,HNEW,IBOUND,CR,CC,CV,HY,TRPY,         
     1                BOT,TOP,SC2,DELR,DELC,DELT,ISS,KITER,KSTP,KPER,           
     2                NCOL,NROW,NLAY,IOUT)                                      
C-----VERSION 1003 10APR1985 BCF1FM                                             
C                                                                               
C     ******************************************************************        
C     ADD LEAKAGE CORRECTION AND STORAGE TO HCOF AND RHS, AND CALCULATE         
C     CONDUCTANCE AS REQUIRED                                                   
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
      DOUBLE PRECISION HNEW                                                     
C                                                                               
      DIMENSION HCOF(NCOL,NROW,NLAY),RHS(NCOL,NROW,NLAY),                       
     1    HOLD(NCOL,NROW,NLAY),SC1(NCOL,NROW,NLAY),HNEW(NCOL,NROW,NLAY),        
     2    IBOUND(NCOL,NROW,NLAY),CR(NCOL,NROW,NLAY),                            
     3    CC(NCOL,NROW,NLAY),CV(NCOL,NROW,NLAY),HY(NCOL,NROW,NLAY),             
     4    TRPY(NLAY),BOT(NCOL,NROW,NLAY),TOP(NCOL,NROW,NLAY),DELR(NCOL),        
     5    DELC(NROW),SC2(NCOL,NROW,NLAY)                                        
C                                                                               
      COMMON /FLWCOM/LAYCON(80)                                                 
C     ------------------------------------------------------------------        
      KB=0                                                                      
      KT=0                                                                      
C                                                                               
C1------FOR EACH LAYER: IF T VARIES CALCULATE HORIZONTAL CONDUCTANCES           
      DO 100 K=1,NLAY                                                           
      IF(LAYCON(K).EQ.3 .OR. LAYCON(K).EQ.2) KT=KT+1                            
C                                                                               
C1A-----IF LAYER TYPE IS NOT 1 OR 3 THEN SKIP THIS LAYER.                       
      IF(LAYCON(K).NE.3 .AND. LAYCON(K).NE.1) GO TO 100                         
      KB=KB+1                                                                   
C                                                                               
C1B-----FOR LAYER TYPES 1 & 3 CALL SBCFH1 TO CALCULATE                          
C1B-----HORIZONTAL CONDUCTANCES.                                                
      CALL SBCF1H(HNEW,IBOUND,CR,CC,CV,HY,TRPY,DELR,DELC,BOT,TOP,               
     1         K,KB,KT,KITER,KSTP,KPER,NCOL,NROW,NLAY,IOUT)                     
  100 CONTINUE                                                                  
C                                                                               
C2------IF THE SIMULATION IS TRANSIENT ADD STORAGE TO HCOF AND RHS              
      IF(ISS.NE.0) GO TO 201                                                    
      TLED=1./DELT                                                              
      KT=0                                                                      
      DO 200 K=1,NLAY                                                           
C                                                                               
C3------SEE IF THIS LAYER IS CONVERTIBLE OR NON-CONVERTIBLE.                    
      IF(LAYCON(K).EQ.3 .OR. LAYCON(K).EQ.2) GO TO 150                          
C4------NON-CONVERTIBLE LAYER, SO USE PRIMARY STORAGE                           
      DO 140 I=1,NROW                                                           
      DO 140 J=1,NCOL                                                           
      RHO=SC1(J,I,K)*TLED                                                       
      HCOF(J,I,K)=HCOF(J,I,K)-RHO                                               
      RHS(J,I,K)=RHS(J,I,K)-RHO*HOLD(J,I,K)                                     
  140 CONTINUE                                                                  
      GO TO 200                                                                 
C                                                                               
C5------A CONVERTIBLE LAYER, SO CHECK OLD AND NEW HEADS TO DETERMINE            
C5------WHEN TO USE PRIMARY AND SECONDARY STORAGE                               
  150 KT=KT+1                                                                   
      DO 180 I=1,NROW                                                           
      DO 180 J=1,NCOL                                                           
C                                                                               
C5A-----IF THE CELL IS EXTERNAL THEN SKIP IT.                                   
      IF(IBOUND(J,I,K).LE.0) GO TO 180                                          
      TP=TOP(J,I,KT)                                                            
      RHO2=SC2(J,I,KT)*TLED                                                     
      RHO1=SC1(J,I,K)*TLED                                                      
C                                                                               
C5B-----FIND STORAGE FACTOR AT START OF TIME STEP.                              
      SOLD=RHO2                                                                 
      IF(HOLD(J,I,K).GT.TP) SOLD=RHO1                                           
C                                                                               
C5C-----FIND STORAGE FACTOR AT END OF TIME STEP.                                
      HTMP=HNEW(J,I,K)                                                          
      SNEW=RHO2                                                                 
      IF(HTMP.GT.TP) SNEW=RHO1                                                  
C                                                                               
C5D-----ADD STORAGE TERMS TO RHS AND HCOF.                                      
      HCOF(J,I,K)=HCOF(J,I,K)-SNEW                                              
      RHS(J,I,K)=RHS(J,I,K) - SOLD*(HOLD(J,I,K)-TP) - SNEW*TP                   
C                                                                               
  180 CONTINUE                                                                  
C                                                                               
  200 CONTINUE                                                                  
C                                                                               
C6------FOR EACH LAYER DETERMINE IF CORRECTION TERMS ARE NEEDED FOR             
C6------FLOW DOWN INTO PARTIALLY SATURATED LAYERS.                              
  201 KT=0                                                                      
      DO 300 K=1,NLAY                                                           
C                                                                               
C7------SEE IF CORRECTION IS NEEDED FOR LEAKAGE FROM ABOVE.                     
      IF(LAYCON(K).NE.3 .AND. LAYCON(K).NE.2) GO TO 250                         
      KT=KT+1                                                                   
      IF(K.EQ.1) GO TO 250                                                      
C                                                                               
C7A-----FOR EACH CELL MAKE THE CORRECTION IF NEEDED.                            
      DO 220 I=1,NROW                                                           
      DO 220 J=1,NCOL                                                           
C                                                                               
C7B-----IF THE CELL IS EXTERNAL(IBOUND<=0) THEN SKIP IT.                        
      IF(IBOUND(J,I,K).LE.0) GO TO 220                                          
      HTMP=HNEW(J,I,K)                                                          
C                                                                               
C7C-----IF HEAD IS ABOVE TOP THEN CORRECTION NOT NEEDED                         
      IF(HTMP.GE.TOP(J,I,KT)) GO TO 220                                         
C                                                                               
C7D-----WITH HEAD BELOW TOP ADD CORRECTION TERMS TO RHS AND HCOF.               
      RHS(J,I,K)=RHS(J,I,K) + CV(J,I,K-1)*TOP(J,I,KT)                           
      HCOF(J,I,K)=HCOF(J,I,K) + CV(J,I,K-1)                                     
  220 CONTINUE                                                                  
C                                                                               
C8------SEE IF THIS LAYER MAY NEED CORRECTION FOR LEAKAGE TO BELOW.             
  250 IF(K.EQ.NLAY) GO TO 300                                                   
      IF(LAYCON(K+1).NE.3 .AND. LAYCON(K+1).NE.2) GO TO 300                     
      KTT=KT+1                                                                  
C                                                                               
C8A-----FOR EACH CELL MAKE THE CORRECTION IF NEEDED.                            
      DO 280 I=1,NROW                                                           
      DO 280 J=1,NCOL                                                           
C                                                                               
C8B-----IF CELL IS EXTERNAL (IBOUND<=0) THEN SKIP IT.                           
      IF(IBOUND(J,I,K).LE.0) GO TO 280                                          
C                                                                               
C8C-----IF HEAD IN THE LOWER CELL IS LESS THAN TOP ADD CORRECTION               
C8C-----TERM TO RHS.                                                            
      HTMP=HNEW(J,I,K+1)                                                        
      IF(HTMP.LT.TOP(J,I,KTT)) RHS(J,I,K)=RHS(J,I,K)                            
     1                        - CV(J,I,K)*(TOP(J,I,KTT)-HTMP)                   
  280 CONTINUE                                                                  
  300 CONTINUE                                                                  
C                                                                               
C9------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE BCF1BD(VBNM,VBVL,MSUM,HNEW,IBOUND,HOLD,SC1,CR,CC,CV,           
     1   TOP,SC2,DELT,ISS,NCOL,NROW,NLAY,KSTP,KPER,IBCFCB,                      
     2   ICBCFL,BUFF,IOUT)                                                      
C-----VERSION 1250 28DEC1983 BCF1BD                                             
C                                                                               
C     ******************************************************************        
C     COMPUTE BUDGET FLOW TERMS FOR BCF -- STORAGE, CONSTANT HEAD, AND          
C     FLOW ACROSS CELL WALLS                                                    
C     ******************************************************************        
C                                                                               
C     SPECIFICATIONS:                                                           
C     ------------------------------------------------------------------
      CHARACTER*4 VBNM,TEXT        
      DOUBLE PRECISION HNEW                                                     
C                                                                               
      DIMENSION HNEW(NCOL,NROW,NLAY), IBOUND(NCOL,NROW,NLAY),                   
     1   HOLD(NCOL,NROW,NLAY), SC1(NCOL,NROW,NLAY),                             
     2   CR(NCOL,NROW,NLAY), CC(NCOL,NROW,NLAY),                                
     3   CV(NCOL,NROW,NLAY), VBNM(4,20), VBVL(4,20),                            
     4   SC2(NCOL,NROW,NLAY),                                                   
     5   TOP(NCOL,NROW,NLAY),BUFF(NCOL,NROW,NLAY)                               
C                                                                               
      COMMON /FLWCOM/LAYCON(80)                                                 
C                                                                               
      DIMENSION TEXT(4)                                                         
C                                                                               
      DATA TEXT(1),TEXT(2),TEXT(3),TEXT(4) /'    ','    ',' STO','RAGE'/        
C     ------------------------------------------------------------------        
C                                                                               
C1------INITIALIZE BUDGET ACCUMULATORS                                          
      STOIN=0.                                                                  
      STOUT=0.                                                                  
C                                                                               
C2------IF CELL-BY-CELL FLOWS ARE NEEDED THEN SET FLAG IBD.                     
      IBD=0                                                                     
      IF(ICBCFL.NE.0 .AND. IBCFCB.GT.0) IBD=1                                   
C                                                                               
C3------IF STEADY STATE THEN SKIP ALL STORAGE CALCULATIONS                      
      IF(ISS.NE.0) GO TO 305                                                    
C                                                                               
C4------IF CELL-BY-CELL FLOWS ARE NEEDED (IBD IS SET) CLEAR BUFFER              
      IF(IBD.EQ.0) GO TO 220                                                    
      DO 210 K=1,NLAY                                                           
      DO 210 I=1,NROW                                                           
      DO 210 J=1,NCOL                                                           
      BUFF(J,I,K)=0.                                                            
  210 CONTINUE                                                                  
C                                                                               
C5------RUN THROUGH EVERY CELL IN THE GRID                                      
  220 KT=0                                                                      
      DO 300 K=1,NLAY                                                           
      LC=LAYCON(K)                                                              
      IF(LC.EQ.3 .OR. LC.EQ.2) KT=KT+1                                          
      DO 300 I=1,NROW                                                           
      DO 300 J=1,NCOL                                                           
C                                                                               
C6------CALCULATE FLOW FROM STORAGE (VARIABLE HEAD CELLS ONLY)                  
      IF(IBOUND(J,I,K).LE.0) GO TO 300                                          
      HSING=HNEW(J,I,K)                                                         
C                                                                               
C6A----CHECK LAYER TYPE TO SEE IF ONE STORAGE FACTOR OR TWO                     
      IF(LC.NE.3 .AND. LC.NE.2) GO TO 285                                       
C                                                                               
C6B----TWO STORAGE FACTORS                                                      
      TP=TOP(J,I,KT)                                                            
      SYA=SC2(J,I,KT)                                                           
      SCFA=SC1(J,I,K)                                                           
      SOLD=SYA                                                                  
      IF(HOLD(J,I,K).GT.TP) SOLD=SCFA                                           
      SNEW=SYA                                                                  
      IF(HSING.GT.TP) SNEW=SCFA                                                 
      STRG=SOLD*(HOLD(J,I,K)-TP) + SNEW*TP - SNEW*HSING                         
      GO TO 288                                                                 
C                                                                               
C6C----ONE STORAGE FACTOR                                                       
  285 SC=SC1(J,I,K)                                                             
      STRG=SC*HOLD(J,I,K) - SC*HSING                                            
C                                                                               
C7-----STORE CELL-BY-CELL FLOW IN BUFFER AND ADD TO ACCUMULATORS                
  288 IF(IBD.EQ.1) BUFF(J,I,K)=STRG/DELT                                        
      IF(STRG) 292,300,294                                                      
  292 STOUT=STOUT-STRG                                                          
      GO TO 300                                                                 
  294 STOIN=STOIN+STRG                                                          
C                                                                               
  300 CONTINUE                                                                  
C                                                                               
C8-----IF IBD FLAG IS SET RECORD THE CONTENTS OF THE BUFFER                     
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,                                  
     1                       IBCFCB,BUFF,NCOL,NROW,NLAY,IOUT)                   
C                                                                               
C9------ADD TOTAL RATES AND VOLUMES TO VBVL & PUT TITLES IN VBNM                
  305 VBVL(1,MSUM)=VBVL(1,MSUM)+STOIN                                           
      VBVL(2,MSUM)=VBVL(2,MSUM)+STOUT                                           
      VBVL(3,MSUM)=STOIN/DELT                                                   
      VBVL(4,MSUM)=STOUT/DELT                                                   
      VBNM(1,MSUM)=TEXT(1)                                                      
      VBNM(2,MSUM)=TEXT(2)                                                      
      VBNM(3,MSUM)=TEXT(3)                                                      
      VBNM(4,MSUM)=TEXT(4)                                                      
      MSUM=MSUM+1                                                               
C                                                                               
C10-----CALCULATE FLOW FROM CONSTANT HEAD NODES                                 
      CALL SBCF1F(VBNM,VBVL,MSUM,HNEW,IBOUND,CR,CC,CV,TOP,DELT,                 
     1        NCOL,NROW,NLAY,KSTP,KPER,IBD,IBCFCB,ICBCFL,BUFF,IOUT)             
C                                                                               
C11-----CALCULATE AND SAVE FLOW ACROSS CELL BOUNDARIES IF C-B-C                 
C11-----FLOW TERMS ARE REQUESTED.                                               
      IF(IBD.NE.0) CALL SBCF1B(HNEW,IBOUND,CR,CC,CV,TOP,NCOL,NROW,NLAY,         
     1         KSTP,KPER,IBCFCB,BUFF,IOUT)                                      
C                                                                               
C12----RETURN                                                                   
      RETURN                                                                    
      END                                                                       
      SUBROUTINE SBCF1C(CR,CC,TRPY,DELR,DELC,K,NCOL,NROW,NLAY)                  
C                                                                               
C                                                                               
C-----VERSION 1010 16NOV1982 SBCF1C                                             
C     ******************************************************************        
C     COMPUTE BRANCH CONDUCTANCE USING HARMONIC MEAN OF BLOCK                   
C     CONDUCTANCES -- BLOCK TRANSMISSIVITY IS IN CC UPON ENTRY                  
C     ******************************************************************        
C                                                                               
C      SPECIFICATIONS:                                                          
C     ------------------------------------------------------------------        
C                                                                               
      DIMENSION CR(NCOL,NROW,NLAY), CC(NCOL,NROW,NLAY)                          
     2   , TRPY(NLAY), DELR(NCOL), DELC(NROW)                                   
C                                                                               
C     ------------------------------------------------------------------        
      YX=TRPY(K)*2.                                                             
C                                                                               
C1------FOR EACH CELL CALCULATE BRANCH CONDUCTANCES FROM THAT CELL              
C1------TO THE ONE ON THE RIGHT AND THE IN FRONT.                               
      DO 40 I=1,NROW                                                            
      DO 40 J=1,NCOL                                                            
      T1=CC(J,I,K)                                                              
C                                                                               
C2------IF T=0 THEN SET CONDUCTANCE EQUAL TO 0. GO ON TO NEXT CELL.             
      IF(T1.NE.0.) GO TO 10                                                     
      CR(J,I,K)=0.                                                              
      GO TO 40                                                                  
C                                                                               
C3------IF THIS IS NOT THE LAST COLUMN(RIGHTMOST) THEN CALCULATE                
C3------BRANCH CONDUCTANCE IN THE ROW DIRECTION (CR) TO THE RIGHT.              
   10 IF(J.EQ.NCOL) GO TO 30                                                    
      T2=CC(J+1,I,K)                                                            
      CR(J,I,K)=2.*T2*T1*DELC(I)/(T1*DELR(J+1)+T2*DELR(J))                      
C                                                                               
C4------IF THIS IS NOT THE LAST ROW(FRONTMOST) THEN CALCULATE                   
C4------BRANCH CONDUCTANCE IN THE COLUMN DIRECTION (CC) TO THE FRONT.           
   30 IF(I.EQ.NROW) GO TO 40                                                    
      T2=CC(J,I+1,K)                                                            
      CC(J,I,K)=YX*T2*T1*DELR(J)/(T1*DELC(I+1)+T2*DELC(I))                      
   40 CONTINUE                                                                  
C                                                                               
C5------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE SBCF1B(HNEW,IBOUND,CR,CC,CV,TOP,NCOL,NROW,NLAY,                
     1      KSTP,KPER,IBCFCB,BUFF,IOUT)                                         
C                                                                               
C-----VERSION 1004 03MAY1983 SBCF1B                                             
C                                                                               
C     ******************************************************************        
C     COMPUTE FLOW ACROSS EACH CELL WALL                                        
C     ******************************************************************        
C                                                                               
C     SPECIFICATIONS:                                                           
C     ------------------------------------------------------------------
      CHARACTER*4 TEXT        
      DOUBLE PRECISION HNEW,HD                                                  
C                                                                               
      DIMENSION HNEW(NCOL,NROW,NLAY), IBOUND(NCOL,NROW,NLAY),                   
     1     CR(NCOL,NROW,NLAY), CC(NCOL,NROW,NLAY),                              
     2     CV(NCOL,NROW,NLAY), TOP(NCOL,NROW,NLAY),                             
     3     BUFF(NCOL,NROW,NLAY)                                                 
C                                                                               
      COMMON /FLWCOM/LAYCON(80)                                                 
C                                                                               
      DIMENSION TEXT(12)                                                        
C                                                                               
      DATA TEXT(1),TEXT(2),TEXT(3),TEXT(4),TEXT(5),TEXT(6),TEXT(7),             
     1   TEXT(8),TEXT(9),TEXT(10),TEXT(11),TEXT(12)                             
     2   /'FLOW',' RIG','HT F','ACE ',                                          
     2    'FLOW',' FRO','NT F','ACE ','FLOW',' LOW','ER F','ACE '/              
C     ------------------------------------------------------------------        
C                                                                               
      NCM1=NCOL-1                                                               
      IF(NCM1.LT.1) GO TO 405                                                   
C                                                                               
C1-----CLEAR THE BUFFER                                                         
      DO 310 K=1,NLAY                                                           
      DO 310 I=1,NROW                                                           
      DO 310 J=1,NCOL                                                           
      BUFF(J,I,K)=0.                                                            
  310 CONTINUE                                                                  
C                                                                               
C2-----FOR EACH CELL CALCULATE FLOW THRU RIGHT FACE & STORE IN BUFFER           
      DO 400 K=1,NLAY                                                           
      DO 400 I=1,NROW                                                           
      DO 400 J=1,NCM1                                                           
      IF((IBOUND(J,I,K).LE.0) .AND. (IBOUND(J+1,I,K).LE.0)) GO TO 400           
      HDIFF=HNEW(J,I,K)-HNEW(J+1,I,K)                                           
      BUFF(J,I,K)=HDIFF*CR(J,I,K)                                               
  400 CONTINUE                                                                  
C                                                                               
C3-----RECORD CONTENTS OF BUFFER                                                
      CALL UBUDSV(KSTP,KPER,TEXT(1),IBCFCB,BUFF,NCOL,NROW,NLAY,IOUT)            
C                                                                               
C4-----CLEAR THE BUFFER                                                         
  405 NRM1=NROW-1                                                               
      IF(NRM1.LT.1) GO TO 505                                                   
      DO 410 K=1,NLAY                                                           
      DO 410 I=1,NROW                                                           
      DO 410 J=1,NCOL                                                           
      BUFF(J,I,K)=0.                                                            
  410 CONTINUE                                                                  
C                                                                               
C5-----FOR EACH CELL CALCULATE FLOW THRU FRONT FACE & STORE IN BUFFER           
      DO 500 K=1,NLAY                                                           
      DO 500 I=1,NRM1                                                           
      DO 500 J=1,NCOL                                                           
      IF((IBOUND(J,I,K).LE.0) .AND. (IBOUND(J,I+1,K).LE.0)) GO TO 500           
      HDIFF=HNEW(J,I,K)-HNEW(J,I+1,K)                                           
      BUFF(J,I,K)=HDIFF*CC(J,I,K)                                               
  500 CONTINUE                                                                  
C                                                                               
C6-----RECORD CONTENTS OF BUFFER.                                               
      CALL UBUDSV(KSTP,KPER,TEXT(5),IBCFCB,BUFF,NCOL,NROW,NLAY,IOUT)            
  505 NLM1=NLAY-1                                                               
      IF(NLM1.LT.1) GO TO 1000                                                  
C                                                                               
C7-----CLEAR THE BUFFER                                                         
      DO 510 K=1,NLAY                                                           
      DO 510 I=1,NROW                                                           
      DO 510 J=1,NCOL                                                           
      BUFF(J,I,K)=0.                                                            
  510 CONTINUE                                                                  
C                                                                               
C8-----FOR EACH CELL CALCULATE FLOW THRU LOWER FACE & STORE IN BUFFER           
      KT=0                                                                      
      DO 600 K=1,NLM1                                                           
      IF(LAYCON(K).EQ.3 .OR. LAYCON(K).EQ.2) KT=KT+1                            
      DO 600 I=1,NROW                                                           
      DO 600 J=1,NCOL                                                           
      IF((IBOUND(J,I,K).LE.0) .AND. (IBOUND(J,I,K+1).LE.0)) GO TO 600           
      HD=HNEW(J,I,K+1)                                                          
      IF(LAYCON(K+1).NE.3 .AND. LAYCON(K+1).NE.2) GO TO 580                     
      TMP=HD                                                                    
      IF(TMP.LT.TOP(J,I,KT+1)) HD=TOP(J,I,KT+1)                                 
  580 HDIFF=HNEW(J,I,K)-HD                                                      
      BUFF(J,I,K)=HDIFF*CV(J,I,K)                                               
  600 CONTINUE                                                                  
C                                                                               
C9-----RECORD CONTENTS OF BUFFER.                                               
      CALL UBUDSV(KSTP,KPER,TEXT(9),IBCFCB,BUFF,NCOL,NROW,NLAY,IOUT)            
C                                                                               
C10----RETURN                                                                   
 1000 RETURN                                                                    
      END                                                                       
      SUBROUTINE SBCF1F(VBNM,VBVL,MSUM,HNEW,IBOUND,CR,CC,CV,                    
     1   TOP,DELT,NCOL,NROW,NLAY,KSTP,KPER,IBD,IBCFCB,ICBCFL,                   
     2   BUFF,IOUT)                                                             
C-----VERSION 1123 29MAY1983 SBCF1F                                             
C                                                                               
C     ******************************************************************        
C     COMPUTE FLOW FROM CONSTANT HEAD NODES                                     
C     ******************************************************************        
C                                                                               
C     SPECIFICATIONS:                                                           
C     ------------------------------------------------------------------
      CHARACTER*4 VBNM,TEXT        
      DOUBLE PRECISION HNEW,HD                                                  
C                                                                               
      DIMENSION HNEW(NCOL,NROW,NLAY), IBOUND(NCOL,NROW,NLAY),                   
     1     CR(NCOL,NROW,NLAY), CC(NCOL,NROW,NLAY),                              
     2     CV(NCOL,NROW,NLAY), VBNM(4,20), VBVL(4,20),                          
     3     TOP(NCOL,NROW,NLAY),BUFF(NCOL,NROW,NLAY)                             
C                                                                               
      COMMON /FLWCOM/LAYCON(80)                                                 
C                                                                               
      DIMENSION TEXT(4)                                                         
C                                                                               
      DATA TEXT(1),TEXT(2),TEXT(3),TEXT(4) /'   C','ONST','ANT ','HEAD'/        
C     ------------------------------------------------------------------        
C                                                                               
C1------CLEAR BUDGET ACCUMULATORS                                               
      CHIN=0.                                                                   
      CHOUT=0.                                                                  
C                                                                               
C2------CLEAR BUFFER IF CELL-BY-CELL FLOW TERM FLAG(IBD) IS SET                 
      IF(IBD.EQ.0) GO TO 8                                                      
      DO 5 K=1,NLAY                                                             
      DO 5 I=1,NROW                                                             
      DO 5 J=1,NCOL                                                             
      BUFF(J,I,K)=0.                                                            
    5 CONTINUE                                                                  
C                                                                               
C3------FOR EACH CELL IF IT IS CONSTANT HEAD COMPUTE FLOW ACROSS 6              
C3-----FACES.                                                                   
    8 KT=0                                                                      
      DO 200 K=1,NLAY                                                           
      LC=LAYCON(K)                                                              
      IF(LC.EQ.3 .OR. LC.EQ.2) KT=KT+1                                          
      DO 200 I=1,NROW                                                           
      DO 200 J=1,NCOL                                                           
C                                                                               
C4-----IF CELL IS NOT CONSTANT HEAD SKIP IT & GO ON TO NEXT CELL.               
      IF (IBOUND(J,I,K).GE.0)GO TO 200                                          
C                                                                               
C5-----CLEAR FIELDS FOR SIX FLOW RATES.                                         
      X1=0.                                                                     
      X2=0.                                                                     
      X3=0.                                                                     
      X4=0.                                                                     
      X5=0.                                                                     
      X6=0.                                                                     
C6-----FOR EACH FACE OF THE CELL CALCULATE FLOW THROUGH THAT FACE               
C6-----OUT OF THE CONSTANT HEAD CELL AND INTO THE FLOW DOMAIN.                  
C6-----COMMENTS 7-11 APPEAR ONLY IN THE SECTION HEADED BY COMMENT 6A            
C6-----BUT THEY APPLY IN A SIMILAR MANNER TO THE SECTIONS HEADED                
C6-----BY COMMENTS 6B-6F.                                                       
C                                                                               
C6A----CALCULATE FLOW THROUGH THE LEFT FACE                                     
C                                                                               
C7-----IF THERE IS NOT A VARIABLE HEAD CELL ON THE OTHER SIDE OF THIS           
C7-----FACE THEN GO ON TO THE NEXT FACE.                                        
      IF(J.EQ.1) GO TO 30                                                       
      IF(IBOUND(J-1,I,K).LE.0)GO TO 30                                          
      HDIFF=HNEW(J,I,K)-HNEW(J-1,I,K)                                           
C                                                                               
C8-----CALCULATE FLOW THROUGH THIS FACE INTO THE ADJACENT CELL.                 
      X1=HDIFF*CR(J-1,I,K)                                                      
C                                                                               
C9-----TEST TO SEE IF FLOW IS POSITIVE OR NEGATIVE                              
      IF (X1) 10,30,20                                                          
C                                                                               
C10----IF NEGATIVE ADD TO CHOUT(FLOW OUT OF DOMAIN TO CONSTANT HEAD).           
   10 CHOUT=CHOUT-X1                                                            
      GO TO 30                                                                  
C                                                                               
C11----IF POSITIVE ADD TO CHIN(FLOW INTO DOMAIN FROM CONSTANT HEAD).            
   20 CHIN=CHIN+X1                                                              
C                                                                               
C6B----CALCULATE FLOW THROUGH THE RIGHT FACE                                    
   30 IF(J.EQ.NCOL) GO TO 60                                                    
      IF(IBOUND(J+1,I,K).LE.0) GO TO 60                                         
      HDIFF=HNEW(J,I,K)-HNEW(J+1,I,K)                                           
      X2=HDIFF*CR(J,I,K)                                                        
      IF(X2)40,60,50                                                            
   40 CHOUT=CHOUT-X2                                                            
      GO TO 60                                                                  
   50 CHIN=CHIN+X2                                                              
C                                                                               
C6C----CALCULATE FLOW THROUGH THE BACK FACE.                                    
   60 IF(I.EQ.1) GO TO 90                                                       
      IF (IBOUND(J,I-1,K).LE.0) GO TO 90                                        
      HDIFF=HNEW(J,I,K)-HNEW(J,I-1,K)                                           
      X3=HDIFF*CC(J,I-1,K)                                                      
      IF(X3) 70,90,80                                                           
   70 CHOUT=CHOUT-X3                                                            
      GO TO 90                                                                  
   80 CHIN=CHIN+X3                                                              
C                                                                               
C6D----CALCULATE FLOW THROUGH THE FRONT FACE.                                   
   90 IF(I.EQ.NROW) GO TO 120                                                   
      IF(IBOUND(J,I+1,K).LE.0) GO TO 120                                        
      HDIFF=HNEW(J,I,K)-HNEW(J,I+1,K)                                           
      X4=HDIFF*CC(J,I,K)                                                        
      IF (X4) 100,120,110                                                       
  100 CHOUT=CHOUT-X4                                                            
      GO TO 120                                                                 
  110 CHIN=CHIN+X4                                                              
C                                                                               
C6E----CALCULATE FLOW THROUGH THE UPPER FACE                                    
  120 IF(K.EQ.1) GO TO 150                                                      
      IF (IBOUND(J,I,K-1).LE.0) GO TO 150                                       
      HD=HNEW(J,I,K)                                                            
      IF(LC.NE.3 .AND. LC.NE.2) GO TO 122                                       
      TMP=HD                                                                    
      IF(TMP.LT.TOP(J,I,KT)) HD=TOP(J,I,KT)                                     
  122 HDIFF=HD-HNEW(J,I,K-1)                                                    
      X5=HDIFF*CV(J,I,K-1)                                                      
      IF(X5) 130,150,140                                                        
  130 CHOUT=CHOUT-X5                                                            
      GO TO 150                                                                 
  140 CHIN=CHIN+X5                                                              
C                                                                               
C6F----CALCULATE FLOW THROUGH THE LOWER FACE.                                   
  150 IF(K.EQ.NLAY) GO TO 180                                                   
      IF(IBOUND(J,I,K+1).LE.0) GO TO 180                                        
      HD=HNEW(J,I,K+1)                                                          
      IF(LAYCON(K+1).NE.3 .AND. LAYCON(K+1).NE.2) GO TO 152                     
      TMP=HD                                                                    
      IF(TMP.LT.TOP(J,I,KT+1)) HD=TOP(J,I,KT+1)                                 
  152 HDIFF=HNEW(J,I,K)-HD                                                      
      X6=HDIFF*CV(J,I,K)                                                        
      IF(X6) 160,180,170                                                        
  160 CHOUT=CHOUT-X6                                                            
      GO TO 180                                                                 
  170 CHIN=CHIN+X6                                                              
C                                                                               
C12-----SUM UP FLOWS THROUGH SIX SIDES OF CONSTANT HEAD CELL.                   
 180  RATE=X1+X2+X3+X4+X5+X6                                                    
C                                                                               
C13-----PRINT THE INDIVIDUAL RATES IF REQUESTED(IBCFCB<0).                      
      IF(IBCFCB.LT.0.AND.ICBCFL.NE.0) WRITE(IOUT,900) (TEXT(N),N=1,4),          
     1    KPER,KSTP,K,I,J,RATE                                                  
  900 FORMAT(1H0,4A4,'   PERIOD',I3,'   STEP',I3,'   LAYER',I3,                 
     1    '   ROW',I4,'   COL',I4,'   RATE ',G15.7)                             
C                                                                               
C14----IF CELL-BY-CELL FLAG SET STORE SUM OF FLOWS FOR CELL IN BUFFER           
      IF(IBD.EQ.1) BUFF(J,I,K)=RATE                                             
C                                                                               
  200 CONTINUE                                                                  
C                                                                               
C15----IF CELL-BY-CELL FLAG SET THEN RECORD CONTENTS OF BUFFER                  
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT(1),                               
     1                   IBCFCB,BUFF,NCOL,NROW,NLAY,IOUT)                       
C                                                                               
C                                                                               
C16----SAVE TOTAL CONSTANT HEAD FLOWS AND VOLUMES IN VBVL TABLE                 
C16----FOR INCLUSION IN BUDGET. PUT LABELS IN VBNM TABLE.                       
      VBVL(1,MSUM)=VBVL(1,MSUM)+CHIN*DELT                                       
      VBVL(2,MSUM)=VBVL(2,MSUM)+CHOUT*DELT                                      
      VBVL(3,MSUM)=CHIN                                                         
      VBVL(4,MSUM)=CHOUT                                                        
C                                                                               
C     ---SETUP VOLUMETRIC BUDGET NAMES                                          
      VBNM(1,MSUM)=TEXT(1)                                                      
      VBNM(2,MSUM)=TEXT(2)                                                      
      VBNM(3,MSUM)=TEXT(3)                                                      
      VBNM(4,MSUM)=TEXT(4)                                                      
C                                                                               
      MSUM=MSUM+1                                                               
C                                                                               
C                                                                               
C17----RETURN                                                                   
      RETURN                                                                    
      END                                                                       
      SUBROUTINE SBCF1H(HNEW,IBOUND,CR,CC,CV,HY,TRPY,DELR,DELC                  
     1,BOT,TOP,K,KB,KT,KITER,KSTP,KPER,NCOL,NROW,NLAY,IOUT)                     
C                                                                               
C-----VERSION 1006 03MAY1983 SBCF1H                                             
C                                                                               
C     ******************************************************************        
C     COMPUTE CONDUCTANCE FROM SATURATED THICKNESS AND HYDRAULIC                
C     CONDUCTIVITY                                                              
C     ******************************************************************        
C                                                                               
C      SPECIFICATIONS:                                                          
C     ------------------------------------------------------------------        
      DOUBLE PRECISION HNEW                                                     
C                                                                               
      DIMENSION HNEW(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY)                     
     1, CR(NCOL,NROW,NLAY), CC(NCOL,NROW,NLAY), CV(NCOL,NROW,NLAY)              
     2,HY(NCOL,NROW,NLAY), TRPY(NLAY), DELR(NCOL), DELC(NROW)                   
     3, BOT(NCOL,NROW,NLAY),TOP(NCOL,NROW,NLAY)                                 
C                                                                               
      COMMON /FLWCOM/LAYCON(80)                                                 
C     ------------------------------------------------------------------        
C                                                                               
C1------CALCULATE TRANSMISSIVITY AT EACH ACTIVE CELL. TRANSMISSIVITY            
C1------WILL BE STORED TEMPORARILY IN THE CC ARRAY.                             
      DO 200 I=1,NROW                                                           
      DO 200 J=1,NCOL                                                           
C                                                                               
C2------IF CELL IS INACTIVE THEN SET T=0 & MOVE ON TO NEXT CELL.                
      IF(IBOUND(J,I,K).NE.0) GO TO 10                                           
      CC(J,I,K)=0.                                                              
      GO TO 200                                                                 
C                                                                               
C3------CALCULATE SATURATED THICKNESS.                                          
   10 HD=HNEW(J,I,K)                                                            
      IF(LAYCON(K).EQ.1) GO TO 50                                               
      IF(HD.GT.TOP(J,I,KT)) HD=TOP(J,I,KT)                                      
   50 THCK=HD-BOT(J,I,KB)                                                       
C                                                                               
C4------CHECK TO SEE IF SATURATED THICKNESS IS GREATER THAN ZERO.               
      IF(THCK.LE.0.) GO TO 100                                                  
C                                                                               
C5------IF SATURATED THICKNESS>0 THEN T=K*THICKNESS.                            
      CC(J,I,K)=THCK*HY(J,I,KB)                                                 
      GO TO 200                                                                 
C                                                                               
C6------WHEN SATURATED THICKNESS < 0, PRINT A MESSAGE AND SET                   
C6------TRANSMISSIVITY, IBOUND, AND VERTICAL CONDUCTANCE =0                     
  100 WRITE(IOUT,150) J,I,K,KITER,KSTP,KPER                                     
  150 FORMAT(1H0,10('*'),'NODE',3I4,' (COL,ROW,LAYER) WENT DRY'                 
     1      ,' AT ITERATION =',I3,'  TIME STEP =',I3                            
     2      ,'  STRESS PERIOD =',I3)                                            
      HNEW(J,I,K)=1.E30                                                         
      CC(J,I,K)=0.                                                              
      IBOUND(J,I,K)=0                                                           
      IF(K.LT.NLAY) CV(J,I,K)=0.                                                
      IF(K.GT.1) CV(J,I,K-1)=0.                                                 
      GO TO 200                                                                 
  200 CONTINUE                                                                  
C                                                                               
C7------COMPUTE HORIZONTAL BRANCH CONDUCTANCES FROM TRANSMISSIVITY              
      CALL SBCF1C(CR,CC,TRPY,DELR,DELC,K,NCOL,NROW,NLAY)                        
C                                                                               
C8------RETURN                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE SBCF1N(HNEW,IBOUND,SC1,SC2,CR,CC,CV,HY,TRPY,DELR,DELC,         
     1         ISS,NCOL,NROW,NLAY,IOUT)                                         
C                                                                               
C-----VERSION 1007 03MAY1983 SBCF1N                                             
C                                                                               
C     ******************************************************************        
C     INITIALIZE AND CHECK BCF DATA                                             
C     ******************************************************************        
C                                                                               
C        SPECIFICATIONS:                                                        
C     ------------------------------------------------------------------        
C                                                                               
      DOUBLE PRECISION HNEW,HCNV                                                
C                                                                               
      DIMENSION HNEW(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY)                     
     1    ,SC1(NCOL,NROW,NLAY),CR(NCOL,NROW,NLAY)                               
     2    ,CC(NCOL,NROW,NLAY),CV(NCOL,NROW,NLAY)                                
     3    ,HY(NCOL,NROW,NLAY),TRPY(NLAY),DELR(NCOL),DELC(NROW)                  
     4    ,SC2(NCOL,NROW,NLAY)                                                  
C                                                                               
      COMMON /FLWCOM/LAYCON(80)                                                 
C     ------------------------------------------------------------------        
C                                                                               
C1------IF IBOUND=0, SET CV=0., CC=0., AND HY=0.                                
      KB=0                                                                      
      DO 30 K=1,NLAY                                                            
      IF(LAYCON(K).EQ.3 .OR. LAYCON(K).EQ.1) KB=KB+1                            
      DO 30 I=1,NROW                                                            
      DO 30 J=1,NCOL                                                            
      IF(IBOUND(J,I,K).NE.0) GO TO 30                                           
      IF(K.NE.NLAY) CV(J,I,K)=0.                                                
      IF(K.NE.1) CV(J,I,K-1)=0.                                                 
      CC(J,I,K)=0.                                                              
      IF(LAYCON(K).EQ.3 .OR. LAYCON(K).EQ.1) HY(J,I,KB)=0.                      
   30 CONTINUE                                                                  
C                                                                               
C2------INSURE THAT EACH ACTIVE CELL HAS AT LEAST ONE NON-ZERO                  
C2------TRANSMISSIVE PARAMETER. IF NOT, CONVERT CELL TO NOFLOW.                 
      HCNV=888.88                                                               
      KB=0                                                                      
      DO 60 K=1,NLAY                                                            
      IF(LAYCON(K).EQ.1 .OR. LAYCON(K).EQ.3) GO TO 55                           
C2A-----WHEN LAYER TYPE 0 OR 2, TRANSMISSIVITY OR CV MUST BE NONZERO            
      DO 54 I=1,NROW                                                            
      DO 54 J=1,NCOL                                                            
      IF(IBOUND(J,I,K).EQ.0) GO TO 54                                           
      IF(CC(J,I,K).NE.0.) GO TO 54                                              
      IF(K.EQ.NLAY) GO TO 51                                                    
      IF(CV(J,I,K).NE.0.) GO TO 54                                              
   51 IF(K.EQ.1) GO TO 53                                                       
      IF(CV(J,I,K-1).NE.0.) GO TO 54                                            
   53 IBOUND(J,I,K)=0                                                           
      HNEW(J,I,K)=HCNV                                                          
      WRITE(IOUT,52) K,I,J                                                      
   52 FORMAT(1X,'NODE (LAYER,ROW,COL)',3I4,                                     
     1      ' ELIMINATED BECAUSE ALL CONDUCTANCES TO NODE ARE 0')               
   54 CONTINUE                                                                  
      GO TO 60                                                                  
C                                                                               
C2B-----WHEN LAYER TYPE IS 1 OR 3, HY OR CV MUST BE NONZERO                     
   55 KB=KB+1                                                                   
      DO 59 I=1,NROW                                                            
      DO 59 J=1,NCOL                                                            
      IF(IBOUND(J,I,K).EQ.0) GO TO 59                                           
      IF(HY(J,I,KB).NE.0.) GO TO 59                                             
      IF(K.EQ.NLAY) GO TO 56                                                    
      IF(CV(J,I,K).NE.0.) GO TO 59                                              
   56 IF(K.EQ.1) GO TO 57                                                       
      IF(CV(J,I,K-1).NE.0.) GO TO 59                                            
   57 IBOUND(J,I,K)=0                                                           
      HNEW(J,I,K)=HCNV                                                          
      CC(J,I,K)=0.                                                              
      WRITE(IOUT,52) K,I,J                                                      
   59 CONTINUE                                                                  
   60 CONTINUE                                                                  
C                                                                               
C3------CALCULATE HOR. CONDUCTANCE(CR AND CC) FOR CONSTANT T LAYERS             
      DO 65 K=1,NLAY                                                            
      IF(LAYCON(K).EQ.3 .OR. LAYCON(K).EQ.1) GO TO 65                           
      CALL SBCF1C(CR,CC,TRPY,DELR,DELC,K,NCOL,NROW,NLAY)                        
   65 CONTINUE                                                                  
C                                                                               
C4------MULTIPLY VERTICAL LEAKANCE BY AREA TO MAKE CONDUCTANCE                  
      IF(NLAY.EQ.1) GO TO 69                                                    
      K1=NLAY-1                                                                 
      DO 68 K=1,K1                                                              
      DO 68 I=1,NROW                                                            
      DO 68 J=1,NCOL                                                            
      CV(J,I,K)=CV(J,I,K)*DELR(J)*DELC(I)                                       
   68 CONTINUE                                                                  
C                                                                               
C5------IF TRANSIENT MULTIPLY PRIMARY STORAGE FACTOR BY DELR &                  
C5------DELC TO GET PRIMARY STORAGE CAPACITY(SC1).                              
   69 IF(ISS.NE.0) GO TO 100                                                    
      KT=0                                                                      
      DO 80 K=1,NLAY                                                            
      DO 70 I=1,NROW                                                            
      DO 70 J=1,NCOL                                                            
      SC1(J,I,K)=SC1(J,I,K)*DELR(J)*DELC(I)                                     
   70 CONTINUE                                                                  
C                                                                               
C6------IF LAYER IS CONF/UNCONF MULTIPLY SECONDARY STORAGE FACTOR               
C6------BY DELR AND DELC TO GET SECONDARY STORAGE CAPACITY(SC2).                
      IF(LAYCON(K).NE.3 .AND. LAYCON(K).NE.2) GO TO 80                          
      KT=KT+1                                                                   
      DO 75 I=1,NROW                                                            
      DO 75 J=1,NCOL                                                            
      SC2(J,I,KT)=SC2(J,I,KT)*DELR(J)*DELC(I)                                   
   75 CONTINUE                                                                  
C                                                                               
   80 CONTINUE                                                                  
C                                                                               
C7------RETURN                                                                  
  100 RETURN                                                                    
      END                                                                       
