

C**********************************************************************C
C                                                                      C
C     *** SUBROUTINE WATERBAL ***                                      C
C     SIMULATES MONTHLY WATER BALANCE USING 2 LAYER SOIL MOISTURE      C
C     ALGORITHM FROM NOAA NATIONAL WEATHER SERVICE SACRAMENTO SOIL     C
C     MOISTURE ACCOUNTING MODEL (SAC-SMA)                              C
C                                                                      C
C**********************************************************************C
      SUBROUTINE WATERBAL(I,J,M,MNDAY)
      
      COMMON/BASIC/NGRID,NYEAR,NLC,BYEAR,IYSTART,IYEND,POP_FLAG,REGION,
     &CLIMATE_FLAG
       
      COMMON/OUTPUT/ PET(200,12,10),APET(12),PAET(200,12,10),APAET(12),
     &AET(12), AVUZTWC(12), AVUZFWC(12), AVLZTWC(12), AVLZFPC(12),
     &AVLZFSC(12),RO(2200,200,12)
          
      COMMON/SNOWPACK/SP(12), SNOWPACK, ANUMAXSWE(200)
       
      COMMON/SUMMARY/ANURAIN(200),ANURUN(200),ANUPET(200),ANUAET(200),
     &ANUSOIL(200),ANUTEMP(200)
       
      COMMON/SOIL/LZTWM(2200), LZFPM(2200), LZFSM(2200), LZSK(2200),
     >  LZPK(2200), UZTWM(2200), UZFWM(2200), UZK(2200), ZPERC(2200),
     >  REXP(2200), PFREE(2200), SMC(12)

      COMMON/CLIMATE/ RAIN(4000,200,12), TEMP(4000,200, 12)
                                       
            
      COMMON/CELLINFO/LADUSE(4000,10),HUCNO(4000),
     >                LATUDE(4000),LONGI(4000),TAREA(4000)
     
      COMMON/IMPERVIOUS/IMPERV(4000,10)

      COMMON/ELEVATION/HUCELE(2200)
      COMMON/HUC/ HUCAREA(2200)
       
      COMMON/FLOW/ STRFLOW(2200, 200, 12)
       
      COMMON/CARBON/ GEPM(2200, 200, 12),RECOM(2200, 200, 12), 
     >               NEEM(2200,200,12)
      
      COMMON/MOMEAN/MOAVGRO(2200,12),MOAVGPET(2200,12),
     &MOAVGAET(2200,12),MOAVGPPT(2200,12),MOAVGSNW(2200,12),
     &MOAVGSOL(2200,12),MOAVGTEM(2200,12)           
C ----------------------------------------------------------------------------     
     
     
      INTEGER HUCNO(4000), MNDAY
      
      INTEGER J, IYSTART, BYEAR
      
      REAL HUCELE(2200), HUCAREA(2200), LATUDE(4000)
      
      REAL RAIN(4000,200,12), TEMP(4000,200, 12)
      
      REAL AETTEMP, RUNOFFTEMP, TBFTEMP, 
     >  IFTEMP, GEPTEMP, RECOTEMP, NEETEMP
      
      REAL LZTWM(2200), LZFPM(2200), LZFSM(2200), LZSK(2200),
     >  LZPK(2200), UZTWM(2200), UZFWM(2200), UZK(2200), ZPERC(2200),
     >  REXP(2200), PFREE(2200)

      REAL UZTWC(10), UZFWC(10), LZTWC(10), LZFSC(10), LZFPC(10)
      
      REAL ETUZTW(200,12,10), RESIDET(200,12,10), ETUZFW(200,12,10)
      
      REAL ETLZTW(200, 12, 10), RATLZT(10), RATLZ(10)
      
      REAL SNOWPACK, SNOWPPT, RAINPPT, SNOWTEMP, RAINTEMP
      
      REAL MELTMAX, SMF, SNOWW, SP(12)
      
      REAL LADUSE(4000,10), TAREA(4000)
      
      REAL ET(200,12,10), SSUR(10), GEP(200,12,10), INFIL(10),
     >     RECO(200, 12, 10), NEE(200, 12,10)      
      
      REAL PET(200,12,10)
      
      REAL UZRAT(10), TWX(10), PERCM(10), PERC(10), DEFR(10), LZDEF(10) 
      
      REAL PERCT(10), PERCF(10)
      
      REAL HPL(10), RATLP(10), RATLS(10), FRACP(10), PERCP(10),PERCS(10)
      
      REAL BF(10), SBF(10), SPBF(10), SIF(10)
      
      REAL TAUZTWC, TAUZFWC, TALZTWC, TALZFPC, TALZFSC
      
      REAL ASM, AET(12), RUNOFF(12), TOTBF(12), INTF(12), 
     &SMC(12)
      
      REAL AVUZTWC(12), AVUZFWC(12), AVLZTWC(12), AVLZFPC(12), 
     >     AVLZFSC(12)
      
       REAL STRFLOW(2200,200,12),GEPM(2200,200,12),RECOM(2200,200,12),
     > NEEM(2200,200,12)
     
      REAL SIDE, CHBF(12)
      
      INTEGER NINC, WRR
      
      REAL DINC,PINC,DUZ,DLZP,DLZS,SPERC(10),DEL(10)
      
      CHARACTER*5 CLIMATE_FLAG
      
      CHARACTER*2 REGION      
 
      REAL MOAVGRO(2200,12),MOAVGPET(2200,12),MOAVGAET(2200,12),
     &MOAVGPPT(2200,12),MOAVGSNW(2200,12),MOAVGSOL(2200,12),
     &MOAVGTEM(2200,12)  
     
      REAL IMPERV(4000,10)

C *****************************************************************************************************
C *****************************************************************************************************
C *****************************************************************************************************
C --- INITIALIZE VARIABLES FOR START OF SIMULATION

      
      AETTEMP =0.
      RUNOFFTEMP = 0.
      TBFTEMP = 0.0
      IFTEMP = 0.0
                          
      GEPTEMP = 0.
      RECOTEMP = 0.
      NEETEMP =0.
             
      IF (J .EQ. ((IYSTART-1)-BYEAR+1) .AND. M .EQ. 1) THEN
                  
C---ESTIMATE INITIAL SNOWPACK (DECEMBER OF  YEAR BEFORE FIRST YEAR OF SIMULATION) 
C---AS A FUNCTION OF LATITUDE (DECIMAL DEGREES)
C---INITIAL SNOWPACK = 0.0002*EXP(0.2644*LATUDE(I))
C---R2=0.42, N=1559
C---BASED ON MEASURED NSIDC 12/1999 SNOWPACK

        SNOWPACK = 0.0002*EXP(0.2644*LATUDE(I))

C---SET INITIAL SOIL MOISTURE

        DO 50 K=1, NLC

           UZTWC(K) = UZTWM(I)
           UZFWC(K) = 1.0
           LZTWC(K) = LZTWM(I)
           LZFSC(K) = 0.75*LZFSM(I)
           LZFPC(K) = 0.75*LZFPM(I)
           SSUR(K) = 0.
               
50      CONTINUE

      ENDIF
        
         
C *****************************************************************************************************
C *****************************************************************************************************
C----- SIMULATE SNOWPACK (SEE McCABE AND WOLOCK 2007)
C ---- (A MONTHLY WATER-BALANCE MODEL DRIVEN BY A GRAPHICAL USER INTERFACE)


C----- SET INITIAL SNOW AND RAIN PPT, SNOWW TO ZERO

        SNOWPPT = 0.
        RAINPPT = 0.
        SNOWW = 0.
      
      IF(REGION.EQ.'US') THEN

        WRR=INT(HUCNO(I)/1000000)

C----- ESTABLISH TEMPERATURE ABOVE WHICH ALL PPT FALLS AS RAIN

        IF (WRR.EQ.1.OR.WRR.EQ.2.OR.WRR.EQ.3.OR.WRR.EQ.5.OR.WRR.EQ.6.OR.
     &WRR.EQ.8) THEN
        
           RAINTEMP = 2.4
        
        ELSEIF(WRR.EQ.4.OR.WRR.EQ.7.OR.WRR.EQ.9) THEN
        
           RAINTEMP = 4.1
           
        ELSEIF(WRR.EQ.10) THEN
        
           RAINTEMP = 9.5
           
        ELSEIF(WRR.EQ.11) THEN
        
           RAINTEMP = 8.7
           
        ELSEIF(WRR.EQ.12.OR.WRR.EQ.13.OR.WRR.EQ.14.OR.WRR.EQ.15.OR.WRR.
     &EQ.16.OR.WRR.EQ.17.OR.WRR.EQ.18) THEN
        
           RAINTEMP = 8.1
        
        ENDIF
 
C----- ESTABLISH TEMPERATURE BELOW WHICH ALL PPT FALLS AS SNOW
        
        IF (WRR.EQ.1.OR.WRR.EQ.2.OR.WRR.EQ.3.OR.WRR.EQ.5.OR.WRR.EQ.6.OR.
     &WRR.EQ.8) THEN
        
           SNOWTEMP = -12.6
        
        ELSEIF(WRR.EQ.4.OR.WRR.EQ.7.OR.WRR.EQ.9) THEN
        
           SNOWTEMP = -9.0
           
        ELSEIF(WRR.EQ.10) THEN
        
           IF(HUCELE(I).GT.1000) THEN
           
              SNOWTEMP = -6.8
           
           ELSE
           
              SNOWTEMP = -8.6
              
           ENDIF
           
        ELSEIF(WRR.EQ.11) THEN
        
           SNOWTEMP = -8.2
           
        ELSEIF(WRR.EQ.12.OR.WRR.EQ.13.OR.WRR.EQ.14.OR.WRR.EQ.15.OR.WRR.
     &EQ.16.OR.WRR.EQ.17) THEN
        
           IF(HUCELE(I).GT.1000) THEN
           
              SNOWTEMP = -5.6
           
           ELSE
           
              SNOWTEMP = -8.6
           
           ENDIF
       
           
        ELSEIF(WRR.EQ.18) THEN
        
           IF(HUCELE(I).GT.1000) THEN
           
              SNOWTEMP = -2.5
              
           ELSE
           
              SNOWTEMP = -8.6
              
           ENDIF
        
        ENDIF
        
C----- ESTABLISH MAXIMUM SNOW MELT RATE

        IF (WRR.EQ.1.OR.WRR.EQ.2.OR.WRR.EQ.3.OR.WRR.EQ.5.OR.WRR.EQ.6.OR.
     &WRR.EQ.8) THEN
        
           MELTMAX = 0.85
        
        ELSEIF(WRR.EQ.4.OR.WRR.EQ.7.OR.WRR.EQ.9) THEN
        
           MELTMAX = 0.90
           
        ELSEIF(WRR.EQ.11) THEN
        
           MELTMAX = 0.75
           
        ELSEIF(WRR.EQ.10.OR.WRR.EQ.12.OR.WRR.EQ.13.OR.WRR.EQ.14.OR.WRR.
     &EQ.16.OR.WRR.EQ.17.OR.WRR.EQ.18) THEN
        
           MELTMAX = 0.60
           
        ELSEIF(WRR.EQ.15) THEN
        
           MELTMAX = 0.50
        
        ENDIF

        ELSEIF(REGION.EQ.'MX') THEN
        
           RAINTEMP = 3.3
           
           IF (HUCELE(I) .LT. 1000.0) THEN
        
              SNOWTEMP = -10.0
        
           ELSE
        
              SNOWTEMP = -1.0
        
           ENDIF
           
           MELTMAX = 0.5
           
        ENDIF


C----- FOR TEMPERATURE LESS THAN SNOW TEMP, CALCULATE SNOWPPT

        IF (TEMP(I, J, M) .LE. SNOWTEMP) THEN
        
           SNOWPPT = RAIN(I, J, M)
             
C----- FOR TEMPERATURE GREATER THAN RAIN TEMP, CALCULATE RAINPPT

        ELSEIF (TEMP(I,J, M) .GE. RAINTEMP) THEN
        
           RAINPPT = RAIN(I, J, M)
        
C----- FOR TEMPERATURE BETWEEN RAINTEMP AND SNOWTEMP, CALCULATE SNOWPPT AND RAINPPT

        ELSE
        
           SNOWPPT = RAIN(I,J, M)*((RAINTEMP - TEMP(I,J, M))/
     &     (RAINTEMP - SNOWTEMP))
           
           RAINPPT = RAIN(I,J,M) - SNOWPPT
           
        ENDIF
        
C----- ACCUMULATE SNOWPPT AS SNOWPACK
           
        SNOWPACK = SNOWPACK + SNOWPPT
        
C----- CALCULATE SNOW MELT FRACTION BASED ON MAXIMUM MELT RATE (MELTMAX) AND MONTHLY TEMPERATURE

        SMF = ((TEMP(I,J, M) - SNOWTEMP) / (RAINTEMP - SNOWTEMP)) 
     &  * MELTMAX
        
        IF (SMF .GT. MELTMAX) THEN
        
           SMF = MELTMAX
           
        ENDIF
        
        IF (SMF .LT. 0.) THEN
        
           SMF = 0.
           
        ENDIF
        
C----- CALCULATE AMOUNT OF SNOW MELTED (MM) TO CONTRIBUTE TO INFILTRATION & RUNOFF

C----- IF SNOWPACK IS LESS THAN 10.0 MM, ASSUME IT WILL ALL MELT (MCCABE AND WOLOCK, 1999)
C----- (GENERAL-CIRCULATION-MODEL SIMULATIONS OF FUTURE SNOWPACK IN THE WESTERN UNITED STATES)

        SNOWW = SNOWPACK * SMF
        
        SNOWPACK = SNOWPACK - SNOWW        

        IF (SNOWPACK .LT. 10.0) THEN
        
           SNOWW = SNOWW+SNOWPACK
        
           SNOWPACK = 0.0
        
          
        ENDIF
        

C *****************************************************************************************************
C *****************************************************************************************************
C -- LOOP THROUGH LAND COVERS IN THE HUC AND PERFORM WATER BALANCE COMPUTATIONS
         
             ASM = 0. 
             TAUZTWC=0.
             TAUZFWC=0.
             TALZTWC=0.
             TALZFPC=0.
             TALZFSC=0.           
             K=0
          
             DO 40 K=1, NLC
  
C *****************************************************************************************************
C *****************************************************************************************************
C -- COMPUTE THE INFILTRATION FOR A GIVEN MONTH FOR EACH LAND USE

                 INFIL(K) = RAINPPT + SNOWW
                 
C --- INITIALIZE FLOWS
 
                  SBF(K)=0.0
                  SSUR(K)=0.0
                  SIF(K)=0.0
                  SPERC(K)=0.0
                  SPBF(K)=0.0
               
                 
C -- FOR OPEN WATER, ET=PAET=PET, INFILTRATION IN EXCESS OF ET GOES TO SURFACE RUNOFF, GEP=RECO=NEE=0.


                   IF(REGION.EQ.'US'.AND.K.EQ.8) THEN
                   
                      ET(J, M, K) = PAET(J,M,K)
                      
                      SSUR(K)=INFIL(K)-ET(J,M,K)
                      
c                      IF (SSUR(K).LT.0) SSUR(K)=0
                      
                      GEP(J,M,K)=0.
                      
                      RECO(J,M,K)=0.
                      
                      NEE(J,M,K)=0.
                      
    
                   ELSE      
            
C *****************************************************************************************************
C --- COMPUTE AET GIVEN TOTAL WATER STORED IN UPPER SOIL LAYER STORAGES AND PAET CALCULATED IN PET.FOR

                   ET(J, M, K) = PAET(J, M, K)
                
C --- COMPUTE ET FROM UZ TENSION WATER STORAGE, RECALCULATE UZTWC, CALCULATE RESIDUAL ET DEMAND

                   ETUZTW(J, M, K) = ET(J, M, K) * (UZTWC(K)/UZTWM(I))
                   
                   RESIDET(J, M, K) = ET(J, M, K) - ETUZTW(J, M, K)
                   
                   UZTWC(K) = UZTWC(K) - ETUZTW(J, M, K)
                   
                   ETUZFW(J,M,K)=0.0
                   
                   IF (UZTWC(K).GE.0.0) GOTO 220
                   
                   ETUZTW(J, M, K) = ETUZTW(J, M, K) + UZTWC(K)
                   
                   UZTWC(K) = 0.0
                   
                   RESIDET(J, M, K) = ET(J, M, K) - ETUZTW(J, M, K)
                   
C --- COMPUTE ET FROM UZ FREE WATER STORAGE, RECALCULATE UZFWC, CALCULATE RESIDUAL ET DEMAND                   
                   
                   IF (UZFWC(K) .GE. RESIDET(J, M, K)) GO TO 221
                   
                   ETUZFW(J, M, K) = UZFWC(K)
                   
                   UZFWC(K) = 0.0
                   
                   RESIDET(J, M, K) = RESIDET(J, M, K) - ETUZFW(J, M, K)
                   
                   GO TO 225
                   
221                ETUZFW(J, M, K) = RESIDET(J, M, K)

                   UZFWC(K) = UZFWC(K) - ETUZFW(J, M, K)
                   
                   RESIDET(J, M, K) = 0.0
                   
C --- REDISTRIBUTE WATER BETWEEN UZ TENSION WATER AND FREE WATER STORAGES

220                IF((UZTWC(K)/UZTWM(I)).GE.(UZFWC(K)/UZFWM(I))) 
     &             GO TO 225

                   UZRAT(K)=(UZTWC(K)+UZFWC(K))/(UZTWM(I)+UZFWM(I))
                      
                   UZTWC(K) = UZTWM(I) * UZRAT(K)
                      
                   UZFWC(K) = UZFWM(I) * UZRAT(K)
                        
225                IF (UZTWC(K) .LT. 0.00001) UZTWC(K) = 0.0

                   IF (UZFWC(K) .LT. 0.00001) UZFWC(K) = 0.0
                   
                   
C --- COMPUTE ET FROM LZ TENSION WATER STORAGE, RECALCULATE LZTWC, CALCULATE RESIDUAL ET DEMAND

                   ETLZTW(J, M, K) = RESIDET(J, M, K) * (LZTWC(K) / 
     &             (UZTWM(I) + LZTWM(I)))
                   
                   LZTWC(K) = LZTWC(K) - ETLZTW(J, M, K)
                   
                   IF(LZTWC(K) .GE. 0.0) GO TO 226
                   
                   ETLZTW(J, M, K) = ETLZTW(J, M, K) + LZTWC(K)
                   
                   LZTWC(K) = 0.0
                   
226                RATLZT(K) = LZTWC(K) / LZTWM(I)

                   RATLZ(K) = (LZTWC(K) + LZFPC(K) + LZFSC(K)) / 
     &             (LZTWM(I) + LZFPM(I) + LZFSM(I))
     
                   IF (RATLZT(K) .GE. RATLZ(K)) GO TO 230
                  
                   LZTWC(K) = LZTWC(K) + (RATLZ(K) - RATLZT(K)) * 
     &             LZTWM(I)

                   
                   LZFSC(K) = LZFSC(K) - (RATLZ(K) - RATLZT(K)) * 
     &             LZTWM(I)
                   
                   IF(LZFSC(K) .GE. 0.0) GO TO 230
                   
                   LZFPC(K) = LZFPC(K) + LZFSC(K)
                   
                   LZFSC(K) = 0.0
                   
230                IF (LZTWC(K) .LT. 0.00001) LZTWC(K) = 0.0

C --- CALCULATE TOTAL ET SUPPLIED BY UPPER AND LOWER LAYERS

                   ET(J, M, K) = ETUZTW(J, M, K) + ETUZFW(J, M, K) + 
     &             ETLZTW(J, M, K)
                   

C *****************************************************************************************************
C *****************************************************************************************************
C --- COMPUTE PERCOLATION INTO SOIL WATER STORAGES

C     --- COMPUTE WATER IN EXCESS OF UZ TENSION WATER CAPACITY (TWX)

                  TWX(K) = INFIL(K) + UZTWC(K) - UZTWM(I)

                  IF (TWX(K).GE.0.0) THEN
             	
C     --- IF INFIL EXCEEDS UZ TENSION WATER CAPACITY, SET UZ TENSION WATER STORAGE TO CAPACITY

                     UZTWC(K) = UZTWM(I)     

                  ELSE        	

C     --- IF INFIL DOES NOT EXCEED UZ TENSION WATER CAPACITY, ALL INFIL GOES TO UZ TENSION WATER STORAGE

                     UZTWC(K) = UZTWC(K) + INFIL(K)
                     
                     TWX(K)=0.0
             	
                  ENDIF
           
C --- DETERMINE COMPUTATIONAL TIME INCREMENTS FOR BASIC TIME INTERVAL

                  NINC=INT(1.0+0.2*(UZFWC(K)+TWX(K)))

      
C     NINC=NUMBER OF TIME INCREMENTS THAT THE TIME INTERVAL
C     IS DIVIDED INTO FOR FURTHER SOIL-MOISTURE ACCOUNTING.  
C     NO ONE INCREMENT WILL EXCEED 5.0 MILLIMETERS OF UZFWC+PAV

                  DINC=(1.0/NINC)
C                 DINC IS THE LENGTH OF EACH INCREMENT IN MONTHS

                  PINC=TWX(K)/NINC
C                 PINC IS THE AMOUNT OF AVAILABLE MOISTURE FOR EACH INCREMENT


C --- COMPUTE FREE WATER DEPLETION COEFICIENTS FOR THE TIME INCREMENT BEING USED

                  DUZ=UZK(I)*DINC
      
                  DLZP=LZPK(I)*DINC
      
                  DLZS=LZSK(I)*DINC
                  

C --- SET UP INCREMENTAL DO LOOP FOR THE TIME INTERVAL

                  DO 240 F=1,NINC
      
C --  COMPUTE BASEFLOW AND KEEP TRACK OF TIME INTERVAL SUM.

                     BF(K) = LZFPC(K) * DLZP
                
                     LZFPC(K) = LZFPC(K) - BF(K)
                
                     IF (LZFPC(K) .LE. 0.0001) THEN 
                
                        BF(K) = BF(K) + LZFPC(K)
                
                        LZFPC(K) = 0.0
                
                     ENDIF

                     SBF(K) = SBF(K) + BF(K)
                 
                     SPBF(K) = SPBF(K) + BF(K)
                 
                     BF(K) = LZFSC(K) * DLZS

                     LZFSC(K) = LZFSC(K) - BF(K)
                
                     IF (LZFSC(K) .LE. 0.0001) THEN
                
                        BF(K) = BF(K) + LZFSC(K)
                
                        LZFSC(K) = 0.0
                
                     ENDIF         
                 
                     SBF(K) = SBF(K) + BF(K)       
                 

C --- COMPUTE PERCOLATION TO LZ IF FREE WATER IS AVAILABLE IN UZ

                  IF (PINC + UZFWC(K).GT.0.01) GOTO 251
		  
                  UZFWC(K) = UZFWC(K) + PINC
		  
                  GOTO 240

C     --- COMPUTE PERCOLATION DEMAND FROM LZ

251               PERCM(K) = LZFPM(I) * DLZP + LZFSM(I) * DLZS
                
                  PERC(K) = PERCM(K) * (UZFWC(K)/UZFWM(I))
                
                  DEFR(K)=1.0-((LZTWC(K)+LZFPC(K)+LZFSC(K))/
     &            (LZTWM(I)+LZFPM(I)+LZFSM(I)))
                
                  PERC(K) = PERC(K) * (1.0 + ZPERC(I) * (DEFR(K)
     &            **REXP(I)))

C     --- COMPARE LZ PERCOLATION DEMAND TO UZ FREE WATER AVAILABLE AND COMPUTE ACTUAL PERCOLATION

                  IF (PERC(K) .GE. UZFWC(K)) THEN
                
                     PERC(K) = UZFWC(K)
                       
                  ENDIF
                       
                  UZFWC(K) = UZFWC(K) - PERC(K)
                
            
C      --- CHECK TO SEE IF PERC EXCEEDS LZ TENSION AND FREE WATER DEFICIENCY, IF SO SET PERC TO LZ DEFICIENCY


                  LZDEF(K) = (LZTWC(K) + LZFPC(K) + LZFSC(K)) - 
     &            (LZTWM(I) + LZFPM(I) + LZFSM(I)) + PERC(K)
                    
                  IF (LZDEF(K) .GT. 0.0) THEN
                    
                     PERC(K) = PERC(K) - LZDEF(K)
          
                     UZFWC(K) = UZFWC(K) + LZDEF(K)
                       
                  ENDIF
                    
                  SPERC(K) = SPERC(K) + PERC(K)
C                 SPERC IS THE TIME INTERVAL SUM OF PERC

C --- COMPUTE INTERFLOW FROM UZ

                  DEL(K) = UZFWC(K) * DUZ
      
                  IF (DEL(K).GT.UZFWC(K))THEN
      
                  DEL(K)=UZFWC(K)
      
                  UZFWC(K)=0.
      
                  ELSE
      
                  UZFWC(K) = UZFWC(K) - DEL(K)
      
                  ENDIF
                 
                  SIF(K) = SIF(K) + DEL(K) 
                                 
C --- DISRIBUTE PERCOLATED WATER INTO THE LZ STORAGES AND COMPUTE THE REMAINDER IN UZ FREE WATER STORAGE AND RESIDUAL AVAIL FOR RUNOFF
   

C     --- COMPUTE PERC WATER GOING INTO LZ TENSION WATER STORAGE AND COMPARE TO AVAILABLE STORAGE

                  PERCT(K) = PERC(K) * (1.0 - PFREE(I))
                
                  IF ((PERCT(K) + LZTWC(K)) .GT. LZTWM(I)) THEN
                
C     --- WHEN PERC IS GREATER THAN AVAILABLE TENSION WATER STORAGE, SET TENSION WATER STORAGE TO MAX, REMAINDER OF PERC GETS EVALUATED AGAINST FREE WATER STORAGE

                     PERCF(K) = PERCT(K) + LZTWC(K) - LZTWM(I)
                
                     LZTWC(K) = LZTWM(I)
                
                  ELSE
                
C     --- WHEN PERC IS LESS THAN AVAILABLE TENSION WATER STORAGE, UPDATE TENSION WATER STORAGE

                     LZTWC(K) = LZTWC(K) + PERCT(K)
                
                     PERCF(K) = 0.0
                
                  ENDIF
                
C     --- COMPUTE TOTAL PERC WATER GOING INTO LZ FREE WATER STORAGE

                  PERCF(K) = PERCF(K) + PERC(K) * PFREE(I)                

                  IF(PERCF(K) .EQ. 0.0) GOTO 245
                
C     --- COMPUTE RELATIVE SIZE OF LZ PRIMARY FREE WATER STORAGE COMPARED TO LZ TOTAL FREE WATER STORAGE

                  HPL(K) = LZFPM(I) / (LZFPM(I) + LZFSM(I))
                
C     --- COMPUTE LZ PRIMARY AND SECONDARY FREE WATER CONTENT TO CAPACITY RATIOS

                  RATLP(K) = LZFPC(K) / LZFPM(I)
                
                  RATLS(K) = LZFSC(K) / LZFSM(I)
                
C     --- COMPUTE FRACTIONS AND PERCENTAGES OF FREE WATER PERC TO GO TO LZ PRIMARY STORAGE

                  FRACP(K) = (HPL(K) * 2.0 * (1.0 - RATLP(K))) 
     >            / ((1.0 - RATLP(K)) + (1.0 - RATLS(K)))
                
                  IF (FRACP(K) .GT. 1.0) FRACP(K) = 1.0

                  PERCP(K) = PERCF(K) * FRACP(K)
                
                  PERCS(K) = PERCF(K) - PERCP(K)
                
C     --- COMPUTE NEW PRIMARY AND SUPPLEMENTAL STORAGE

C         --- COMPUTE NEW SUPPLEMENTAL FREE WATER STORAGE

                  LZFSC(K) = LZFSC(K) + PERCS(K)

                  IF(LZFSC(K) .GT. LZFSM(I)) THEN
                
C         --- IF NEW SUPPLEMENTAL FREE WATER STORAGE EXCEEDS CAPACITY SET SUPPLEMENTAL STORAGE TO CAPACITY AND EXCESS GOES TO PRIMARY FREE WATER STORAGE

                     PERCS(K) = PERCS(K) - LZFSC(K) + LZFSM(I)
                
                     LZFSC(K) = LZFSM(I)
                          
                  ENDIF
                
            
C         --- IF NEW LZ SUPPLEMENTAL FREE WATER STORAGE IS LESS THAN CAPACITY MOVE ON TO COMPUTE NEW PRIMARY FREE WATER STORAGE


                  LZFPC(K) = LZFPC(K) + (PERCF(K) - PERCS(K))

                
                  IF (LZFPC(K) .LE. LZFPM(I)) GOTO 245

C             --- IF LZ FREE PRIMARY WATER STORAGE EXCEEDS CAPACITY SET PRIMARY STORAGE TO CAPACITY AND EVALUATE EXCESS AGAINST LZ TENSION WATER STORAGE

                  LZTWC(K) = LZTWC(K) + LZFPC(K) - LZFPM(I)
                
                  LZFPC(K) = LZFPM(I)
                
C            ---DISTRIBUTE PINC BETWEEN UZFWC AND SURFACE RUNOFF

                       
C                 IF UZFWC CAN ABSORB EXCESS, NO SURFACE RUNOFF

245               IF(PINC.EQ.0.0) GOTO 240

                  IF (PINC+UZFWC(K).GT.UZFWM(I)) GOTO 248
                       
                  UZFWC(K)=UZFWC(K)+PINC
                          
                  GOTO 240
                          
C                 IF UZFWC CAN NOT ABSORB EXCESS, THE COMPUTE SURFACE RUNOFF
                       
248               SSUR(K)=SSUR(K)+PINC+UZFWC(K)-UZFWM(I)

                         
                  UZFWC(K)=UZFWM(I)
                  

240               CONTINUE


C *****************************************************************************************************
C *****************************************************************************************************

C CALCULATING MONTHLY GEP G C/M2/MO
C NOTE  the following is based on New Analysis by Eddy Flux Asko Noormets(Aug 24, 2010)

C 1 = Crop; 2= Deciduous, 3=Evergreen, 4=Mixed  forest
C 5 = Grassland, 6 = SHRUBLANDS ; 7 = WETLAND ; 8= OPEN WATER; 9=URBAN; 10=BARREN

               IF (K.EQ.1) THEN 
              
                  GEP(J, M, K) = 3.13 * ET(J,M,K)
              
                  IF (TEMP(I,J,M) .LE. -1.0) THEN 
                
                     GEP(J, M, K) = 0.0
                
                  ENDIF 
              
C Reco - temp exp function                                                      
                          
                  RECO(J,M,K)= 40.59 + 0.43 * GEP(J,M,K)
                 
C -- DECIDUOUS FOREST
 
               ELSEIF (K .EQ. 2) THEN 

                  GEP(J, M, K) = 3.2 * ET(J,M,K)
             
                  IF (TEMP(I,J,M) .LE. -1.0) THEN 
                
                     GEP(J, M, K) = 0.0
                
                  ENDIF 
             
                  RECO(J,M,K)= 30.81 + 0.45* GEP(J,M,K)
          
C -- EVERGREEN FOREST

               ELSEIF (K .EQ. 3) THEN 
  
                  GEP(J, M, K) = 2.46 * ET(J,M,K)
 
                  IF (TEMP(I,J,M).LE.-1.0.AND.(M.GE.11.OR.M.LE.4)) THEN 
                
                     GEP(J, M, K) = 0.0
                
                  ENDIF 
               
                  RECO(J, M, K) = 9.9 + 0.685 * GEP(J,M,K)
               
C -- MIXED FOREST (AS ANERAGE OF DECIDUOUS AND EVERGREEN)
               
               ELSEIF (K .EQ. 4) THEN 
               
                  GEP(J, M, K) = 2.74 * ET(J,M,K)
             
                  IF (TEMP(I,J,M) .LE. -1.0) THEN 
                
                     GEP(J, M, K) = 0.0
                
                  ENDIF 

                  RECO(J,M,K)= 24.44 + 0.617 * GEP(J,M,K)

C -- GRASSLANDS               
               
               ELSEIF (K .EQ. 5) THEN                
               
                  GEP(J, M, K) = 2.12 * ET(J,M,K)
             
                  IF (TEMP(I,J,M) .LE. -1.0) THEN 
                
                     GEP(J, M, K) = 0.0
                
                  ENDIF 

                  RECO(J,M,K)= 18.9 + 0.64 * GEP(J,M,K)
               
C -- SHRUBLANDS                   
               
               ELSEIF (K .EQ. 6) THEN                

                  GEP(J, M, K) = 1.35 * ET(J,M,K)  
             
                  IF (TEMP(I,J,M) .LE. -1.0) THEN 
                
                     GEP(J, M, K) = 0.0
                
                  ENDIF 

                  RECO(J,M,K)= 9.7 + 0.56 * GEP(J,M,K)
                                           
C -- WETLAND FOR US(SAME AS MIXED????????)                   
               
               ELSEIF (REGION.EQ.'US'.AND.K .EQ. 7) THEN                

                  GEP(J, M, K) = 2.74 * ET(J,M,K)
             
                  IF (TEMP(I,J,M) .LE. -1.0) THEN 
                
                     GEP(J, M, K) = 0.0
                
                  ENDIF 

                  RECO(J,M,K)= 24.44 + 0.617 * GEP(J,M,K)
                               
C -- SAVANNA FOR MX AND RW             
               
               ELSEIF(REGION.EQ.'MX'.OR.REGION.EQ.'RW'.AND.K .EQ. 7)THEN
               
               GEP(J, M, K) = 1.26 * ET(J,M,K)
             
               IF (TEMP(I,J,M) .LE. -1.0) THEN 
                
                GEP(J, M, K) = 0.0
                
               ENDIF 
                              
               RECO(J,M,K)= 25.24 + 0.53 * GEP(J,M,K)


C -- W/U/B FOR MX AND RW 

               ELSEIF(REGION.EQ.'MX'.OR.REGION.EQ.'RW'.AND.K .EQ. 8)THEN

               GEP(J, M, K) = (1.531 * ET(J,M,K) - 11.78)
               
               IF (TEMP(I,J,M) .LE. -1.0) THEN 
                
                   GEP(J, M, K) = 0.0
                
               ENDIF   
                
                RECO(J,M,K)= 9.7 + 0.56 * GEP(J,M,K)
                        
C -- URBAN(SAME AS OPEN SHRUB???)
               
               ELSEIF (REGION.EQ.'US'.AND.K .EQ. 9) THEN                

                  GEP(J, M, K) = 1.35 * ET(J,M,K)  
             
                  IF (TEMP(I,J,M) .LE. -1.0) THEN 
                
                     GEP(J, M, K) = 0.0
                
                  ENDIF 

                  RECO(J,M,K)= 9.7 + 0.56 * GEP(J,M,K)
                                           
C -- BARREN(SAME AS OPEN SHRUB???)
               
               ELSEIF (REGION.EQ.'US'.AND.K .EQ. 10) THEN                

                  GEP(J, M, K) = 1.35 * ET(J,M,K)  
             
                  IF (TEMP(I,J,M) .LE. -1.0) THEN 
                
                     GEP(J, M, K) = 0.0
                
                  ENDIF 

                  RECO(J,M,K)= 9.7 + 0.56 * GEP(J,M,K)
                             
               ENDIF

               IF (GEP(J,M,K) .LE. 0.) THEN
               
                  GEP(J,M,K) = 0.                  
               
               ENDIF
               
               
C Calcualte ecossytem respiration and NEE per MONTH                                 
 
          
               NEE(J, M, K) = -GEP(J,M,K) + RECO(J,M,K)    


C *****************************************************************************************************
C *****************************************************************************************************
C --- COMPUTE FRACTION OF EACH WATER BALANCE COMPONENT AND GEP FOR EACH LAND COVER

C *****************************************************************************************************
      IF(REGION.EQ.'US'.AND.LADUSE(I,8).NE.TAREA(I)) THEN

C *****************************************************************************************************
C -- ACCUMULATE THE LAND COVER AREA WEIGHTED AVERAGE GEP FOR THE HUC

               GEPTEMP = GEPTEMP + GEP(J,M,K) * (LADUSE(I,K)*
     &(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,8))
     
               RECOTEMP = RECOTEMP + RECO(J,M,K)*(LADUSE(I,K)*
     &(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,8))
     
               NEETEMP = NEETEMP + NEE(J,M,K)*(LADUSE(I,K)*
     &(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,8))

C --- ACCUMULATE THE LAND COVER AREA WEIGHTED AVERAGE AET FOR THE HUC (WEIGHTED TO LAND AREA ONLY, OPEN WATER NOT INCLUDED)
               
             AETTEMP = AETTEMP+ET(J,M,K)*(LADUSE(I,K)*(1-IMPERV(I,K)))/
     &(TAREA(I)-LADUSE(I,8))
C *****************************************************************************************************
C--- ACCUMULATE THE LAND COVER AREA WEIGHTED AVERAGE FLOWS FOR THE HUC (WEIGHTED TO LAND AREA ONLY, OPEN WATER NOT INCLUDED)
               
C--- SURFACE RUNOFF FROM NON-IMPERVIOUS AREAS

            RUNOFFTEMP=RUNOFFTEMP+SSUR(K)*(LADUSE(I,K)*(1-IMPERV(I,K)))/
     &(TAREA(I)-LADUSE(I,8))
                
C--- SURFACE RUNOFF FROM IMPERVIOUS AREAS

             RUNOFFTEMP=RUNOFFTEMP+INFIL(K)*(LADUSE(I,K)*(IMPERV(I,K)))/
     &(TAREA(I)-LADUSE(I,8))

C--- INTERFLOW FROM NON-IMPERVIOUS AREAS

               IFTEMP = IFTEMP + SIF(K) *(LADUSE(I,K)*(1-IMPERV(I,K)))/
     &(TAREA(I)-LADUSE(I,8))

C--- TOTAL BASEFLOW FROM NON-IMPERVIOUS AREAS
              
               TBFTEMP = TBFTEMP + SBF(K)*(LADUSE(I,K)*(1-IMPERV(I,K)))/
     &(TAREA(I)-LADUSE(I,8))               
               
              
C *****************************************************************************************************
C--- ACCUMULATE THE LAND COVER AREA WEIGHTED AVERAGE SOIL MOISTURE FOR THE HUC
C--- NON-IMPERVIOUS AREAS ONLY (WEIGHTED TO LAND AREA ONLY, OPEN WATER NOT INCLUDED)

      TAUZTWC = TAUZTWC + UZTWC(K) * (LADUSE(I,K)*
     &(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,8))
              
      TAUZFWC = TAUZFWC + UZFWC(K) * (LADUSE(I,K)*
     &(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,8))
              
      TALZTWC = TALZTWC + LZTWC(K) * (LADUSE(I,K)*
     &(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,8))
              
      TALZFPC = TALZFPC + LZFPC(K) * (LADUSE(I,K)*
     &(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,8))
              
      TALZFSC = TALZFSC + LZFSC(K) * (LADUSE(I,K)*
     &(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,8))
              
      ASM = ASM + (UZTWC(K)+UZFWC(K)+LZTWC(K)+LZFPC(K)
     &+LZFSC(K)) * (LADUSE(I,K)*(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,8))
     
      ELSEIF(REGION.EQ.'US'.AND.LADUSE(I,8).EQ.TAREA(I)) THEN
      
      RUNOFFTEMP=SSUR(8)
      IFTEMP=0.
      TBFTEMP=0.
      TAUZTWC=0.
      TAUZFWC=0.
      TALZTWC=0.
      TALZFPC=0.
      TALZFSC=0.
      ASM=0.
      GEPTEMP=0.
      RECOTEMP=0.
      NEETEMP=0.
      
      ELSEIF(REGION.EQ.'MX'.OR.REGION.EQ.'RW') THEN
      
      GEPTEMP = GEPTEMP + GEP(J,M,K) * LADUSE(I,K)/TAREA(I)
      RECOTEMP = RECOTEMP + RECO(J,M,K)*LADUSE(I,K)/TAREA(I)               
      NEETEMP = NEETEMP + NEE(J,M,K)*LADUSE(I,K)/TAREA(I)
      AETTEMP = AETTEMP + ET(J,M,K) * LADUSE(I,K)/TAREA(I)
      RUNOFFTEMP = RUNOFFTEMP + SSUR(K) * LADUSE(I,K)/TAREA(I)
      TBFTEMP = TBFTEMP + SBF(K) * LADUSE(I,K)/TAREA(I)
      IFTEMP = IFTEMP + SIF(K) * LADUSE(I,K)/TAREA(I)
      TAUZTWC = TAUZTWC + UZTWC(K) * LADUSE(I,K)/TAREA(I)
      TAUZFWC = TAUZFWC + UZFWC(K) * LADUSE(I,K)/TAREA(I)
      TALZTWC = TALZTWC + LZTWC(K) * LADUSE(I,K)/TAREA(I)
      TALZFPC = TALZFPC + LZFPC(K) * LADUSE(I,K)/TAREA(I)
      TALZFSC = TALZFSC + LZFSC(K) * LADUSE(I,K)/TAREA(I)
      ASM = ASM + (UZTWC(K)+UZFWC(K)+LZTWC(K)+LZFPC(K)
     &+LZFSC(K)) * LADUSE(I,K)/TAREA(I)
      
      ENDIF
             
      ENDIF
      
40            CONTINUE

C -- CALCULATE MONTHLY AVERAGE WATER BALANCE COMPONENTS (LAND AREA ONLY, OPEN WATER NOT INCLUDED)

              AET(M) = AETTEMP        
              RUNOFF(M) = RUNOFFTEMP
              TOTBF(M) = TBFTEMP
              INTF(M) = IFTEMP
              SMC (M) = ASM          
              SP(M) = SNOWPACK
              AVUZTWC(M) = TAUZTWC
              AVUZFWC(M) = TAUZFWC
              AVLZTWC(M) = TALZTWC
              AVLZFPC(M) = TALZFPC
              AVLZFSC(M) = TALZFSC
                   
              GEPM(I,J, M) = GEPTEMP
              RECOM(I,J,M)  = RECOTEMP
              NEEM(I,J,M) = NEETEMP


      IF(REGION.EQ.'US') THEN

C -- CALCULATE THE TOTAL HUC AET, RUNOFF, SOIL STORAGE FOR BOTH LAND AND WATER, WEIGHTED BY LAND AND WATER AREA

      AET(M) = AET(M)*(TAREA(I)-LADUSE(I,8))/TAREA(I) + 
     &ET(J,M,8)*LADUSE(I,8)/TAREA(I)
      
      RUNOFF(M)=RUNOFF(M)*(TAREA(I)-LADUSE(I,8))/TAREA(I)+
     &SSUR(8)*LADUSE(I,8)/TAREA(I)
     
      TOTBF(M)=TOTBF(M)*(TAREA(I)-LADUSE(I,8))/TAREA(I)
      
      INTF(M)=INTF(M)*(TAREA(I)-LADUSE(I,8))/TAREA(I)
     
      AVUZTWC(M)=AVUZTWC(M)*(TAREA(I)-LADUSE(I,8))/TAREA(I)
      AVUZFWC(M)=AVUZFWC(M)*(TAREA(I)-LADUSE(I,8))/TAREA(I)
      AVLZTWC(M)=AVLZTWC(M)*(TAREA(I)-LADUSE(I,8))/TAREA(I)
      AVLZFPC(M)=AVLZFPC(M)*(TAREA(I)-LADUSE(I,8))/TAREA(I)
      AVLZFSC(M)=AVLZFSC(M)*(TAREA(I)-LADUSE(I,8))/TAREA(I)
      SMC (M)=SMC (M)*(TAREA(I)-LADUSE(I,8))/TAREA(I)

      GEPM(I,J,M) = GEPM(I,J,M)*(TAREA(I)-LADUSE(I,8))/TAREA(I) + 
     &GEP(J,M,8)*LADUSE(I,8)/TAREA(I)
      RECOM(I,J,M) = RECOM(I,J,M)*(TAREA(I)-LADUSE(I,8))/TAREA(I) + 
     &RECO(J,M,8)*LADUSE(I,8)/TAREA(I)
      NEEM(I,J,M) = NEEM(I,J,M)*(TAREA(I)-LADUSE(I,8))/TAREA(I) + 
     &NEE(J,M,8)*LADUSE(I,8)/TAREA(I)
     
      ENDIF


C -- CALCULATE AMOUNT OF BF THAT GOES TO STREAM CHANNEL
C -- REMAINDER IS LOST TO DEEP GROUNDWATER

C -- SIDE = RATIO OF DEEP RECHARGE TO CHANNEL BASEFLOW (DIMENSIONLESS)
c           SIDE = 0.85

              SIDE = 0.00
           
              CHBF(M) = (TOTBF(M))*(1-SIDE)
 
C -- CALCULATE TOTAL RUNOFF (RO) IN MM FOR EACH HUC FOR MONTH M

              RO(I,J,M) = RUNOFF(M) + CHBF(M) + INTF(M)
              
c              IF(HUCNO(I).EQ.3130001.AND.J.NE.21) THEN
c              WRITE(77,*) HUCNO(I),J,M,RUNOFF(M),CHBF(M),INTF(M)
c              ENDIF

              

C -- SET TOTAL WATER YIELD TO ZERO IF NEGATIVE.  THIS COULD HAPPEN IF PET>PPT FOR OPEN WATER DEPENDING ON THE % OF HUC IN OPEN WATER
C -- MAY WANT TO REVISIT THIS FOR ROUTING, COULD ACCOUNT FOR INSTREAM FLOW LOSSES DUE TO ET IN ARID AREAS.

              IF(RO(I,J,M).LT.0) THEN
              
              RO(I,J,M)=0.
              
              ENDIF
          
C -- CALCULATE STREAMFLOW IN MILLION M3 FOR EACH HUC FOR MONTH M

              STRFLOW(I, J, M) = RO(I,J,M) * HUCAREA(I)/1000./1000000.
      

C --  ACCUMULATE MONTHLY MEAN WATER BALANCE OUTPUT
      
         IF(J.GT.(IYSTART-BYEAR)) THEN
         
         MOAVGRO(I,M)=MOAVGRO(I,M)+RO(I,J,M)
         MOAVGPET(I,M)=MOAVGPET(I,M)+APET(M)
         MOAVGAET(I,M)=MOAVGAET(I,M)+AET(M)
         MOAVGPPT(I,M)=MOAVGPPT(I,M)+RAIN(I,J,M)
         MOAVGTEM(I,M)=MOAVGTEM(I,M)+TEMP(I,J,M)
         MOAVGSNW(I,M)=MOAVGSNW(I,M)+SP(M)
         MOAVGSOL(I,M)=MOAVGSOL(I,M)+AVUZTWC(M)+AVUZFWC(M)+
     &AVLZTWC(M)+AVLZFPC(M)+AVLZFSC(M)
     
         ENDIF


      RETURN
      END
