C**********************************************************************C
C                                                                      C
C     *** SUBROUTINE FLOWROUTING ***                                   C
C     SIMULATE TOTAL FLOW FOR EACH MONTH AND EACH HUC                  C
C     WRITE ACCUMULATED FLOW TO MONTHACCFLOW.TXT                       C
C     PERFORM WATER SUPPLY/DEMAND AND WASSI CALCULATIONS               C
C     WRITE WASSI OUTPUT TO ANNUALWaSSI.TXT                            C
C     WRITE WASSI OUTPUT TO HUCWaSSI.TXT                               C
C**********************************************************************C
C        I=HUC; J= YEAR, M =MONTH, MNDAY= NO OF DAYS IN A MONTH

      SUBROUTINE FLOWROUTING_US
      
      COMMON/BASIC/NGRID,NYEAR,NLC,BYEAR,IYSTART,IYEND,POP_FLAG,REGION,
     &CLIMATE_FLAG

      COMMON/FLOW/ STRFLOW(2200, 200, 12)        

      COMMON/HUC/ HUCAREA(2200)
      
      COMMON/CELLINFO/LADUSE(4000,10),HUCNO(4000),
     >                LATUDE(4000),LONGI(4000),TAREA(4000)
         
      COMMON/WATERUSE/ WATERUSE(2200,10), HUCWUSE(2200),PERCAP(2200)
      
      COMMON/POPULATION/POPULATION(2200, 200)
      
      COMMON/GROUNDWATER/ GROUNDWATER(2200,9), GW_CHG
            
      COMMON/HUCREGION/DMC(18,4),IRRC(18,4), PPC(18,4) 

      COMMON/RETURNFLOW/ RETURNFLOW(2200,8)
      
  
C --------------------------------------------------------------
      INTEGER M
      
      REAL STRFLOW(2200, 200, 12)
      
      INTEGER X,Y,HUCID(50,2000), NFLWPATHS
  
      INTEGER HUCNO(4000), BYEAR, I, POP_FLAG
      
      INTEGER IYSTART, IYEND, J, R,NDAY,MNDAY
      
      REAL POPULATION(2200,200),WU,PERCAP(2200)
                   
      REAL HUCAREA(2200), RETURNFLOW(2200,8),WATERUSE(2200, 10)
      
      REAL GROUNDWATER (2200, 9)
      
      REAL SP1, SP2, SP3, DM

      REAL SUPPLY(2200,200,12), GROSSDEM(2200,200,12), 
     &NETDEM(2200,200,12), MWASSI(2200,200,12)
            
      REAL MONPER, DMC(18,4), IRRC(18,4), PPC(18,4)
      
      REAL ANSUPPLY(2200,200), ANGROSDEM(2200,200), ANWASSI(2200,200),
     &ANNETDEM(2200,200) 
      
      REAL AHUCSUPPLY, AHUCGROSDEM, AHUCNETDEM, AHUCFLOW, AHUCWASSI, 
     &AHUCSWS, AHUCGWS
      
      REAL MOSUPPLY, MODEMAND, CONSUSE
      
      REAL ANFLOW(2200,200)
      
      REAL GWS(2200,200,12), ANGWS(2200,200), ANSWS(2200,200)
      
      REAL DM_STORAG(2200,200,12),ANDMSTOR(2200,200),AHUCDMSTOR
      
      INTEGER YEAR1, COUNT(2200)
      
      INTEGER MONTHD(12),MONTHL(12)
      
      CHARACTER*5 CLIMATE_FLAG
      
      CHARACTER*2 REGION
      
      REAL MOAVGSWS(2200,12),MOAVGGWS(2200,12),MOAVGSUP(2200,12),
     &MOAVGGDEM(2200,12),MOAVGNDEM(2200,12),MOAVGWAS(2200,12),
     &MOAVGFLW(2200,12), SWS(2200,200,12),MOAVGDMST(2200,12)

C --- Number of days for each month during regular year
      DATA MONTHD/31,28,31,30,31,30,31,31,30,31,30,31/

C --- Number of days for each month during leap year
      DATA MONTHL/31,29,31,30,31,30,31,31,30,31,30,31/


C-----PRINT TITLE FOR MONTHLY WASSI OUTPUT MONTHWaSSI.TXT

      WRITE (100, 3000)     
3000  FORMAT   ('CELL',',','YEAR',',','MONTH',',','SWS_MGD',
     &',','GWS_MGD',',','TOTSUP_MGD',',','GROSSDEM_MGD',',',
     &'NETDEM_MGD',',','WASSI',',','CONSUSE_Mm3_MO',
     &',','FLWOUT_Mm3_MO',',','DEMSTOR_Mm3_MO')
     
C-----PRINT TITLE FOR ANNUAL WASSI OUTPUT ANNUALWaSSI.TXT

      WRITE (200, 4000)     
4000  FORMAT   ('CELL',',','YEAR',',','SWS_MGD',
     &',','GWS_MGD',',','TOTSUP_MGD',',','GROSSDEM_MGD',',',
     &'NETDEM_MGD',',','WASSI',',','CONSUSE_Mm3_YR',
     &',','FLWOUT_Mm3_YR',',','DEMSTOR_Mm3_YR') 
     
C-----PRINT TITLE FOR HUC WASSI OUTPUT AVGANNWaSSI.TXT
     
      WRITE (300, 5000)     
5000  FORMAT ('CELL',',','YRS',',','SWS_MGD',
     &',','GWS_MGD',',','TOTSUP_MGD',',','GROSSDEM_MGD',',',
     &'NETDEM_MGD',',','WASSI',',','CONSUSE_Mm3_YR',
     &',','FLWOUT_Mm3_YR',',','DEMSTOR_Mm3_YR') 
     
C-----PRINT TITLE FOR MONTHLY MEAN WASSI OUTPUT MOMEANWASSI.TXT
     
      WRITE (804, 6000)     
6000  FORMAT ('CELL',',','MO',',','SWS_MGD',
     &',','GWS_MGD',',','TOTSUP_MGD',',','GROSSDEM_MGD',',',
     &'NETDEM_MGD',',','WASSI',',','CONSUSE_Mm3_MO',
     &',','FLWOUT_Mm3_MO',',','DEMSTOR_Mm3_MO') 

C --- READ IN ROUTING NETWORK FROM ROUTEMATRIX.TXT

      DO 5062 Y=1,1000000
     
      
      READ (17,*,END=5063) (HUCID(X,Y), X=1,33)
     
C      WRITE(77,5061) HUCID(1,Y),HUCID(2,Y),HUCID(3,Y),HUCID(4,Y),
C     &HUCID(5,Y),HUCID(6,Y),HUCID(7,Y),HUCID(8,Y),HUCID(9,Y),
C     &HUCID(10,Y),HUCID(11,Y),HUCID(12,Y),HUCID(13,Y),HUCID(14,Y),
C     &HUCID(15,Y),HUCID(16,Y),HUCID(17,Y),HUCID(18,Y),HUCID(19,Y),
C     &HUCID(20,Y),HUCID(21,Y),HUCID(22,Y),HUCID(23,Y),HUCID(24,Y),
C     &HUCID(25,Y),HUCID(26,Y),HUCID(27,Y),HUCID(28,Y),HUCID(29,Y),
C     &HUCID(30,Y),HUCID(31,Y),HUCID(32,Y),HUCID(33,Y)
      
5061  FORMAT (33I10)

5062  CONTINUE

5063  NFLWPATHS=Y-1

c      DO 5064 X=1,33
c      DO 5065 Y=1,NFLWPATHS
      

c      IF(HUCID(X,Y).NE.-9999.AND.HUCID(X-1,Y).NE.-9999.AND.
c     &HUCID(X-1,Y).NE.0) THEN

c      WRITE(77,*) HUCID(X,Y),HUCID(X-1,Y)
      
c      ENDIF
      
c5065  CONTINUE
c5064  CONTINUE



C************************************************************************
C --  PERFORM SUPPLY, DEMAND, WASSI, AND FLOW ROUTING CALCULATIONS ACCOUNTING FOR CONSUMPTIVE USE 
C************************************************************************
     
      YEAR1=1960

      DO 198 I=1,NGRID
      
         DO 199 M=1,12
         
         MOAVGSWS(I,M)=0.
         MOAVGGWS(I,M)=0.
         MOAVGSUP(I,M)=0.
         MOAVGGDEM(I,M)=0.
         MOAVGNDEM(I,M)=0.
         MOAVGWAS(I,M)=0.
         MOAVGFLW(I,M)=0.
         MOAVGDMST(I,M)=0.

         
199      CONTINUE

198   CONTINUE

      DO 200 J=(IYSTART-BYEAR+1),(IYEND-BYEAR+1)
      
      YEAR = J + BYEAR - 1
      
C --- DETERMINE WHETHER YEAR IS A LEAP YEAR 
C --- http://www.timeanddate.com/date/leapyear.html

             IF (YEAR/4 .EQ. INT(YEAR/4)) THEN
             
                NDAY = 366
             
                IF (YEAR/100 .EQ. INT(YEAR/100)) THEN
                
                   NDAY = 365
                   
                   IF (YEAR/400 .EQ. INT(YEAR/400)) THEN
                   
                     NDAY = 366
                     
                   ENDIF
                   
                ENDIF
                
             ELSE
             
               NDAY = 365
               
             ENDIF

                    
         DO 100 M=1, 12
         
               IF (NDAY .EQ. 365) THEN
                 MNDAY=MONTHD(M)
               ELSE
                 MNDAY=MONTHL(M)
               ENDIF

         
            X=33
      
            DO 50 WHILE(X.GE.1)

C --- SET INITIAL COUNT FOR HUCID IN COLUMN Y TO ZERO

               DO 29 R=1,NGRID
                  
                 COUNT(R)=0.
                    
29             CONTINUE


               DO 60 Y=1,NFLWPATHS

                  IF (HUCID(X,Y).EQ.-9999) THEN
      
                     GOTO 60
      
                  ENDIF
       
                  IF(X.LT.33) THEN
      
                     IF (HUCID(X+1,Y).NE.-9999) THEN
                       
                        STRFLOW(HUCID(X,Y),J,M)=
     &                  STRFLOW(HUCID(X,Y),J,M)+
     &                  STRFLOW(HUCID(X+1,Y),J,M)   
                       
                     ELSE
      
                        STRFLOW(HUCID(X,Y),J,M)=
     &                  STRFLOW(HUCID(X,Y),J,M)
      
                     ENDIF
               
                  ELSE
            
                     STRFLOW(HUCID(X,Y),J,M)=STRFLOW(HUCID(X,Y),J,M)
               
                  ENDIF

C--CALCULATE WATER RESOURCE REGION (IR) FOR HUC 

                  IR = INT (HUCNO(HUCID(X,Y))/1000000)                             
                                          
C--SP1 --> SURFACE WATER SUPPLY IN GALLON; 1 M3 = 264.17 GALLON Streamflow is in Million m3/MONTH
                               
                  SP1=STRFLOW(HUCID(X,Y), J, M)*264.17
          
                  IF (SP1 .LE. 0.) THEN SP1=0.
          
                  SP2 = 0.
                  SP3 = 0.
                  DM = 0. 

C----CALCULATE MONTHLY PERCENT OF TOTAL USE FOR EACH SECTOR
                                 
                  DO 920 IS = 1, 8
                  
                                           
C -- FOR INDUSTRIAL(IS=2), LIVESTOCK(IS=4), MINING(IS=5), 
C    AQUACULTURE(IS=8), DISTRIBUTE ANNUAL WU EVENLY ACROSS MONTHS

                  IF (IS .EQ. 2 .OR. IS .EQ. 4 .OR. IS .EQ.
     >            5 .OR. IS .EQ. 8) THEN 
           
                     MONPER = 1.0/12.0
               
C --- FOR DOMESTIC WATER USE (IS = 1) AND PUBLIC SUPPLY(IS=7) CALCUALTE MONTHLY % USE
C----PUBLIC SUPPLY GOES TO DOMESTIC AND INDUSTRIAL USERS         
C----ASSUME ALL GOES TO DOMESTIC, AND FOLLOWS THE MONTHLY DISTRIBUTION OF DOMESTIC WATER USE
        
                  ELSEIF (IS .EQ. 1.OR. IS .EQ. 7) THEN  
               
                     MONPER=DMC(IR,1)*M**3+DMC(IR,2)*M**2+ DMC(IR,3)*M+
     >               DMC(IR,4)
                     
C -- CORRECT DOMESTIC MONPER SO THAT TOTAL ANNUAL DOMESTIC MONPER=1

                     IF(IR.EQ.1.OR.IR.EQ.2)THEN
                     MONPER=MONPER/0.95244
                     ELSEIF(IR.EQ.3)THEN
                     MONPER=MONPER/1.001
                     ELSEIF(IR.EQ.4)THEN
                     MONPER=MONPER/0.92692
                     ELSEIF(IR.EQ.5)THEN
                     MONPER=MONPER/0.98776
                     ELSEIF(IR.EQ.6)THEN
                     MONPER=MONPER/0.91312
                     ELSEIF(IR.EQ.7)THEN
                     MONPER=MONPER/1.0908
                     ELSEIF(IR.EQ.8)THEN
                     MONPER=MONPER/1.0956
                     ELSEIF(IR.EQ.9.OR.IR.EQ.10)THEN
                     MONPER=MONPER/0.951
                     ELSEIF(IR.EQ.11.OR.IR.EQ.16)THEN
                     MONPER=MONPER/0.9176
                     ELSEIF(IR.EQ.12.OR.IR.EQ.13.OR.IR.EQ.15)THEN
                     MONPER=MONPER/0.97
                     ELSEIF(IR.EQ.14)THEN
                     MONPER=MONPER/1.0636
                     ELSEIF(IR.EQ.17)THEN
                     MONPER=MONPER/0.8916
                     ELSEIF(IR.EQ.18)THEN
                     MONPER=MONPER/1.2742
                     ENDIF
                     

C----FOR IRRIGATION (IS= 3)

                     ELSEIF (IS .EQ. 3) THEN 
        
                        MONPER=IRRC(IR,1)*M**3+IRRC(IR,2)*M**2+
     &                  IRRC(IR,3)*M+IRRC(IR,4)
    
C----FOR THE NORTH STATES IRRIGATION IN WINTER MONTHS (OCT-MAY)= 0. 
      
                        IF (IR .EQ. 10 .OR. IR .EQ. 9 .OR. IR .EQ. 7 
     >                  .OR. IR .EQ. 2.OR. IR .EQ. 1 .OR. IR .EQ. 16
     >                  .OR. IR .EQ. 14) THEN
     
                           IF (M .LE. 4 .OR. M .GE. 11) THEN 
          
                              MONPER = 0.
          
                           ELSE IF (M .EQ. 7  .OR. M .EQ. 8) THEN
          
                              MONPER = 0.30
          
                           ELSEIF (M .EQ. 6 .OR. M .EQ. 9) THEN
          
                              MONPER = 0.15
                
                           ELSEIF (M .EQ. 5 .OR. M .EQ. 10) THEN      
                
                              MONPER = 0.05      
          
                           ENDIF
                              
                        ENDIF 
           
C -- CORRECT IRRIGATION MONPER SO THAT TOTAL ANNUAL IRRIGATION MONPER=1

                     IF(IR.EQ.3.OR.IR.EQ.6.OR.IR.EQ.11)THEN
                     MONPER=MONPER/0.99674
                     ELSEIF(IR.EQ.4.OR.IR.EQ.5)THEN
                     MONPER=MONPER/1.0246
                     ELSEIF(IR.EQ.8)THEN
                     MONPER=MONPER/1.2466
                     ELSEIF(IR.EQ.12.OR.IR.EQ.13.OR.IR.EQ.15)THEN
                     MONPER=MONPER/1.0646
                     ELSEIF(IR.EQ.17)THEN
                     MONPER=MONPER/0.9352
                     ELSEIF(IR.EQ.18)THEN
                     MONPER=MONPER/1.00376
                     ENDIF


C----FOR POWER PLANT (IS=6)
     
                     ELSEIF(IS .EQ. 6) THEN
                   
                        MONPER=PPC(IR,1)*M**3+PPC(IR,2)*M**2+ 
     &                  PPC(IR,3)*M + PPC(IR,4)
                    
C -- CORRECT THERMO MONPER SO THAT TOTAL ANNUAL THERMO MONPER=1

                     IF(IR.EQ.1.OR.IR.EQ.2.OR.IR.EQ.17)THEN
                     MONPER=MONPER/0.994
                     ELSEIF(IR.EQ.3)THEN
                     MONPER=MONPER/1.02402
                     ELSEIF(IR.EQ.4.OR.IR.EQ.5)THEN
                     MONPER=MONPER/1.0406
                     ELSEIF(IR.EQ.6)THEN
                     MONPER=MONPER/0.8788
                     ELSEIF(IR.EQ.7)THEN
                     MONPER=MONPER/1.0674
                     ELSEIF(IR.EQ.8)THEN
                     MONPER=MONPER/1.0006
                     ELSEIF(IR.EQ.9.OR.IR.EQ.10)THEN
                     MONPER=MONPER/1.0918
                     ELSEIF(IR.EQ.11.OR.IR.EQ.14.OR.IR.EQ.16)THEN
                     MONPER=MONPER/1.07112
                     ELSEIF(IR.EQ.12.OR.IR.EQ.13.OR.IR.EQ.15.OR.
     &               IR.EQ.18)THEN
                     MONPER=MONPER/0.913
                     ENDIF


                     ENDIF
       
                     IF (MONPER .LT. 0) THEN 
        
                        MONPER = 0.
       
                     ENDIF
        
C -- IF DOMESTIC WATER USE IS TO BE CALCULATED BASED ON POPULATION (POP_FLAG NOT EQUAL TO 9000)
C -- RATHER THAN USING THE WATER USE DATA IN WATERUSE.TXT, THEN CALCULATE DOMESTIC WATER USE

                     IF (IS.EQ.1.AND.POP_FLAG.NE.9000) THEN

C -- CALCULATE TOTAL DOMESTIC WATER USE AS A FUNCTION
C -- OF POPULATION AND ASSIGN TO THE DOMESTIC WATER USE SECTOR

                WATERUSE(HUCID(X,Y),1)= 
     &          POPULATION(HUCID(X,Y), BYEAR-YEAR1+J)*PERCAP(HUCID(X,Y))
     
                
     
                     ENDIF

C--CALCULATE MONTHLY WATER USE FOR EACH SECTOR BASED ON MONPER
         
                     WU = WATERUSE(HUCID(X,Y),IS) * NDAY * MONPER    

C--SUM WATER USE FOR ALL SECTORS TO GET TOTAL DEMAND (DM)

                     DM = DM + WU                    
                        
C--CALCULATE ADDITIONAL SUPPLY AVAILABLE FROM RETURNFLOW (FRACTION OF WU RETURNED)
                                             
                     SP2 = SP2 + RETURNFLOW(HUCID(X,Y), IS) * WU 
                            
C --CALCULATE ADDITIONAL SUPPLY AVAILABLE FROM GROUNDWATER

                     SP3 = SP3 + GROUNDWATER(HUCID(X,Y), IS)*NDAY
     &               * MONPER
                            
920               CONTINUE


C -- CALCULATE SWS, GWS, TOTAL SUPPLY, TOTAL GROSS DEMAND, AND TOTAL NET DEMAND             

                  SWS(HUCID(X,Y),J,M)= SP1
                  
                  GWS(HUCID(X,Y),J,M) = SP3
                  
                  SUPPLY(HUCID(X,Y),J,M) = SP1 + SP3
                  
                  GROSSDEM(HUCID(X,Y),J,M) = DM
                  
                  NETDEM(HUCID(X,Y),J,M) = DM - SP2

          
                  IF (COUNT(HUCID(X,Y)).LE.0.) THEN
      
C -- RECALCULATE STREAMFLOW IN Mm3/MO FOR HUC=HUCID(X,Y) ACCOUNTING FOR NET SURFACE WATER USE
C -- NET SURFACE WATER USE = TOTAL GROSS WATER USE-RETURNFLOW-GROUNDWATER WITHDRAWAL
                  
                     STRFLOW(HUCID(X,Y),J,M)=STRFLOW(HUCID(X,Y),J,M)+
     &               (SP2 + SP3 - GROSSDEM(HUCID(X,Y),J,M))/264.17
    
                     IF(STRFLOW(HUCID(X,Y),J,M).LT.0.) THEN
               
                        DM_STORAG(HUCID(X,Y),J,M)=
     &                  -STRFLOW(HUCID(X,Y),J,M)
                     
                        STRFLOW(HUCID(X,Y),J,M)=0.0
                  
                     ENDIF
     
                  ENDIF
                  
                  COUNT(HUCID(X,Y))=COUNT(HUCID(X,Y))+1.


C -- CALCULATE SUPPLY AND DEMAND IN UNITS OF MGD FOR THE MONTH

               GROSSDEM(HUCID(X,Y),J,M) = GROSSDEM(HUCID(X,Y),J,M)/MNDAY
               
               NETDEM(HUCID(X,Y),J,M) = NETDEM(HUCID(X,Y),J,M)/MNDAY
                  
               SUPPLY(HUCID(X,Y),J,M) = SUPPLY(HUCID(X,Y),J,M)/MNDAY
                  
               GWS(HUCID(X,Y),J,M) = GWS(HUCID(X,Y),J,M)/MNDAY 
               
               SWS(HUCID(X,Y),J,M) = SWS(HUCID(X,Y),J,M)/MNDAY 

C----FOR HUCS BORDERING THE GREAT LAKES, SET MONTHLY WASSI TO 0.0

         IF(HUCNO(HUCID(X,Y)).EQ.4010101.OR.HUCNO(HUCID(X,Y)).EQ.4010301
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4010302.OR.HUCNO(HUCID(X,Y)).EQ.4020101
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4020102.OR.HUCNO(HUCID(X,Y)).EQ.4020103
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4020105.OR.HUCNO(HUCID(X,Y)).EQ.4020201
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4020202.OR.HUCNO(HUCID(X,Y)).EQ.4020203
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4020300.OR.HUCNO(HUCID(X,Y)).EQ.4030101
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4030102.OR.HUCNO(HUCID(X,Y)).EQ.4030103
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4030104.OR.HUCNO(HUCID(X,Y)).EQ.4030105
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4030108.OR.HUCNO(HUCID(X,Y)).EQ.4030109
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4030110.OR.HUCNO(HUCID(X,Y)).EQ.4030111
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4030112.OR.HUCNO(HUCID(X,Y)).EQ.4030204
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4040001.OR.HUCNO(HUCID(X,Y)).EQ.4040002
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4040003.OR.HUCNO(HUCID(X,Y)).EQ.4050002
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4050006.OR.HUCNO(HUCID(X,Y)).EQ.4060101
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4060103.OR.HUCNO(HUCID(X,Y)).EQ.4060104
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4060105.OR.HUCNO(HUCID(X,Y)).EQ.4060106
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4060107.OR.HUCNO(HUCID(X,Y)).EQ.4070001
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4070002.OR.HUCNO(HUCID(X,Y)).EQ.4070003
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4070004.OR.HUCNO(HUCID(X,Y)).EQ.4070007
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4080101.OR.HUCNO(HUCID(X,Y)).EQ.4080102
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4080103.OR.HUCNO(HUCID(X,Y)).EQ.4080104
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4080206.OR.HUCNO(HUCID(X,Y)).EQ.4090001
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4090003.OR.HUCNO(HUCID(X,Y)).EQ.4120102
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4120103.OR.HUCNO(HUCID(X,Y)).EQ.4120104
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4130001.OR.HUCNO(HUCID(X,Y)).EQ.4130003
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4140101.OR.HUCNO(HUCID(X,Y)).EQ.4140102
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4140203.OR.HUCNO(HUCID(X,Y)).EQ.4150102
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4150301.OR.HUCNO(HUCID(X,Y)).EQ.4150302
     &  .OR.HUCNO(HUCID(X,Y)).EQ.7120001.OR.HUCNO(HUCID(X,Y)).EQ.7120003
     &  .OR.HUCNO(HUCID(X,Y)).EQ.7120004.OR.HUCNO(HUCID(X,Y)).EQ.7120005
     &  .OR.HUCNO(HUCID(X,Y)).EQ.7120006.OR.HUCNO(HUCID(X,Y)).EQ.4110004
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4110003.OR.HUCNO(HUCID(X,Y)).EQ.4100010
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4100009.OR.HUCNO(HUCID(X,Y)).EQ.4050001
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4110002.OR.HUCNO(HUCID(X,Y)).EQ.4120101
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4090002.OR.HUCNO(HUCID(X,Y)).EQ.4100001
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4110001.OR.HUCNO(HUCID(X,Y)).EQ.4090004
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4100011.OR.HUCNO(HUCID(X,Y)).EQ.4100012
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4090005.OR.HUCNO(HUCID(X,Y)).EQ.4100002
     &  .OR.HUCNO(HUCID(X,Y)).EQ.4050003) THEN
     
                  MWASSI(HUCID(X,Y),J,M)=0.0
                  
                  ELSE
                  
C----FOR ALL OTHER HUCS, CALCULATE MONTHLY WASSI

                     MWASSI(HUCID(X,Y),J,M) = GROSSDEM(HUCID(X,Y),J,M)/
     &               SUPPLY(HUCID(X,Y),J,M)
                  
                  ENDIF

60             CONTINUE  
      
               X=X-1
    
50          CONTINUE  

100      CONTINUE 

200   CONTINUE

C ---WRITE OUTPUT TO OUTPUT FILES

      DO 600 I=1, NGRID         
                
         DO 500 J=(IYSTART-BYEAR+1),(IYEND-BYEAR+1)
         
            IDY = J + BYEAR - 1 
            ANSUPPLY(I,J) = 0.
            ANGROSDEM(I,J) = 0.
            ANNETDEM(I,J) = 0.
            MOSUPPLY = 0.
            MODEMAND = 0.
            ANFLOW(I,J) = 0.
            ANGWS(I,J) = 0.
            ANSWS(I,J) = 0.
            ANDMSTOR(I,J) = 0.
            
C --- DETERMINE WHETHER YEAR IS A LEAP YEAR 
C --- http://www.timeanddate.com/date/leapyear.html

             IF (IDY/4 .EQ. INT(IDY/4)) THEN
             
                NDAY = 366
             
                IF (IDY/100 .EQ. INT(IDY/100)) THEN
                
                   NDAY = 365
                   
                   IF (IDY/400 .EQ. INT(IDY/400)) THEN
                   
                     NDAY = 366
                     
                   ENDIF
                   
                ENDIF
                
             ELSE
             
               NDAY = 365
               
             ENDIF


            DO 400 M=1, 12
            
               IF (NDAY .EQ. 365) THEN
                 MNDAY=MONTHD(M)
               ELSE
                 MNDAY=MONTHL(M)
               ENDIF

            
C----ONLY WRITE OUTPUT FOR YEARS AFTER THE MODEL WARMUP PERIOD

      IF(J.GT.(IYSTART-BYEAR)) THEN

C ---WRITE MONTHLY SUPPLY, DEMAND, AND WASSI TO MONTHWASSI.TXT
               
               CONSUSE = (-NETDEM(I,J,M)+GWS(I,J,M))*MNDAY/264.17
               
               WRITE (100, 2000) HUCNO(I),IDY,M,SWS(I,J,M),GWS(I,J,M),
     &         SUPPLY(I,J,M),GROSSDEM(I,J,M), NETDEM(I,J,M),
     &         MWASSI(I,J,M),CONSUSE,STRFLOW(I,J,M),DM_STORAG(I,J,M)
     
2000        FORMAT (I10,',',I10,',',I10,',',F16.2,',',F16.2,',',
     &         F16.2,',',F16.2,',',F16.2,',',F10.4,',',F16.2,
     &         ',',F16.2,',',F16.2) 
     
       ENDIF

C--- ACCUMULATE ANNUAL WATER DEMAND, SUPPLY, AND FLOW

               ANSUPPLY(I,J) = ANSUPPLY(I,J) + SUPPLY(I,J,M)
               ANGROSDEM(I,J) = ANGROSDEM(I,J) + GROSSDEM(I,J,M)  
               ANNETDEM(I,J) = ANNETDEM(I,J) + NETDEM(I,J,M) 
               ANFLOW(I,J) = ANFLOW(I,J) + STRFLOW(I,J,M)
               ANGWS(I,J) = ANGWS(I,J) + GWS(I,J,M)
               ANSWS(I,J) = ANSWS(I,J)+SWS(I,J,M)
               ANDMSTOR(I,J) = ANDMSTOR(I,J)+DM_STORAG(I,J,M)

C --  ACCUMULATE MONTHLY MEAN WASSI AND FLOWROUTING OUTPUT
      
         MOAVGSWS(I,M) = MOAVGSWS(I,M) + SWS(I,J,M)
         MOAVGGWS(I,M) = MOAVGGWS(I,M) + GWS(I,J,M)
         MOAVGSUP(I,M) = MOAVGSUP(I,M) + SUPPLY(I,J,M)
         MOAVGGDEM(I,M) = MOAVGGDEM(I,M) + GROSSDEM(I,J,M)
         MOAVGNDEM(I,M) = MOAVGNDEM(I,M) + NETDEM(I,J,M)
         MOAVGWAS(I,M) = MOAVGWAS(I,M) + MWASSI(I,J,M)
         MOAVGFLW(I,M) = MOAVGFLW(I,M) + STRFLOW(I,J,M)
         MOAVGDMST(I,M) = MOAVGDMST(I,M)+DM_STORAG(I,J,M)


400         CONTINUE

C -- CALCULATE THE ANNUAL SUPPLY AND DEMAND IN MGD FOR EACH YEAR

            ANSUPPLY(I,J) = ANSUPPLY(I,J) / 12.
            ANGROSDEM(I,J) = ANGROSDEM(I,J) / 12.
            ANNETDEM(I,J) = ANNETDEM(I,J) / 12.
            ANGWS(I,J) = ANGWS(I,J)/12.
            ANSWS(I,J) = ANSWS(I,J)/12.

C----FOR HUCS BORDERING THE GREAT LAKES, SET ANNUAL WASSI TO 0.0

                  IF(HUCNO(I).EQ.4010101.OR.HUCNO(I).EQ.4010301
     &            .OR.HUCNO(I).EQ.4010302.OR.HUCNO(I).EQ.4020101
     &            .OR.HUCNO(I).EQ.4020102.OR.HUCNO(I).EQ.4020103
     &            .OR.HUCNO(I).EQ.4020105.OR.HUCNO(I).EQ.4020201
     &            .OR.HUCNO(I).EQ.4020202.OR.HUCNO(I).EQ.4020203
     &            .OR.HUCNO(I).EQ.4020300.OR.HUCNO(I).EQ.4030101
     &            .OR.HUCNO(I).EQ.4030102.OR.HUCNO(I).EQ.4030103
     &            .OR.HUCNO(I).EQ.4030104.OR.HUCNO(I).EQ.4030105
     &            .OR.HUCNO(I).EQ.4030108.OR.HUCNO(I).EQ.4030109
     &            .OR.HUCNO(I).EQ.4030110.OR.HUCNO(I).EQ.4030111
     &            .OR.HUCNO(I).EQ.4030112.OR.HUCNO(I).EQ.4030204
     &            .OR.HUCNO(I).EQ.4040001.OR.HUCNO(I).EQ.4040002
     &            .OR.HUCNO(I).EQ.4040003.OR.HUCNO(I).EQ.4050002
     &            .OR.HUCNO(I).EQ.4050006.OR.HUCNO(I).EQ.4060101
     &            .OR.HUCNO(I).EQ.4060103.OR.HUCNO(I).EQ.4060104
     &            .OR.HUCNO(I).EQ.4060105.OR.HUCNO(I).EQ.4060106
     &            .OR.HUCNO(I).EQ.4060107.OR.HUCNO(I).EQ.4070001
     &            .OR.HUCNO(I).EQ.4070002.OR.HUCNO(I).EQ.4070003
     &            .OR.HUCNO(I).EQ.4070004.OR.HUCNO(I).EQ.4070007
     &            .OR.HUCNO(I).EQ.4080101.OR.HUCNO(I).EQ.4080102
     &            .OR.HUCNO(I).EQ.4080103.OR.HUCNO(I).EQ.4080104
     &            .OR.HUCNO(I).EQ.4080206.OR.HUCNO(I).EQ.4090001
     &            .OR.HUCNO(I).EQ.4090003.OR.HUCNO(I).EQ.4120102
     &            .OR.HUCNO(I).EQ.4120103.OR.HUCNO(I).EQ.4120104
     &            .OR.HUCNO(I).EQ.4130001.OR.HUCNO(I).EQ.4130003
     &            .OR.HUCNO(I).EQ.4140101.OR.HUCNO(I).EQ.4140102
     &            .OR.HUCNO(I).EQ.4140203.OR.HUCNO(I).EQ.4150102
     &            .OR.HUCNO(I).EQ.4150301.OR.HUCNO(I).EQ.4150302
     &            .OR.HUCNO(I).EQ.7120001.OR.HUCNO(I).EQ.7120003
     &            .OR.HUCNO(I).EQ.7120004.OR.HUCNO(I).EQ.7120005
     &            .OR.HUCNO(I).EQ.7120006.OR.HUCNO(I).EQ.4110004
     &            .OR.HUCNO(I).EQ.4110003.OR.HUCNO(I).EQ.4100010
     &            .OR.HUCNO(I).EQ.4100009.OR.HUCNO(I).EQ.4050001
     &            .OR.HUCNO(I).EQ.4110002.OR.HUCNO(I).EQ.4120101
     &            .OR.HUCNO(I).EQ.4090002.OR.HUCNO(I).EQ.4100001
     &            .OR.HUCNO(I).EQ.4110001.OR.HUCNO(I).EQ.4090004
     &            .OR.HUCNO(I).EQ.4100011.OR.HUCNO(I).EQ.4100012
     &            .OR.HUCNO(I).EQ.4090005.OR.HUCNO(I).EQ.4100002
     &            .OR.HUCNO(I).EQ.4050003) THEN
     
                 
                  ANWASSI(I,J) = 0.0
                  
                  ELSE
                  
C----FOR ALL OTHER HUCS, CALCULATE ANNUAL WASSI

                     IF (ANSUPPLY(I,J) .NE. 0.) THEN 
        
                        ANWASSI(I,J) = ANGROSDEM(I,J) / ANSUPPLY(I,J)
    
                     ELSE
        
                        ANWASSI(I,J) = 10.
                      
                     ENDIF 
                     
                  ENDIF
                  
       
C----ONLY WRITE OUTPUT FOR YEARS AFTER THE MODEL WARMUP PERIOD

      IF(J.GT.(IYSTART-BYEAR)) THEN


C -- WRITE  ANNUAL WASSI, SUPPLY, AND DEMAND 
C -- FOR EACH YEAR TO ANNUALWaSSI.TXT
       
            CONSUSE = (-ANNETDEM(I,J)+ANGWS(I,J))*NDAY/264.17
            
            WRITE (200, 7000) HUCNO(I), IDY, 
     &ANSWS(I,J), ANGWS(I,J), ANSUPPLY(I,J),ANGROSDEM(I,J), 
     &ANNETDEM(I,J),ANWASSI(I,J), CONSUSE,ANFLOW(I,J),ANDMSTOR(I,J) 
                                              
7000        FORMAT (I10,',',I10,',',F16.2,',',F16.2,',',F16.2,',',
     &      F16.2,',',F16.2,',',F10.4,',',F16.2,',',F10.1,',',F10.1)

      ENDIF

500      CONTINUE

         AHUCSUPPLY = 0.
         AHUCGROSDEM = 0.
         AHUCNETDEM = 0.
         AHUCFLOW = 0.
         AHUCWASSI = 0.
         AHUCGWS=0.
         AHUCSWS=0.
         AHUCDMSTOR=0.
         
         IY = 0
         
         DO 75 J = IYSTART,IYEND
         
            AHUCSUPPLY = AHUCSUPPLY + ANSUPPLY(I,J-(BYEAR-1))
            AHUCGROSDEM = AHUCGROSDEM + ANGROSDEM(I,J-(BYEAR-1))
            AHUCNETDEM = AHUCNETDEM + ANNETDEM(I,J-(BYEAR-1))
            AHUCFLOW = AHUCFLOW + ANFLOW(I,J-(BYEAR-1))
            AHUCWASSI = AHUCWASSI + ANWASSI(I,J-(BYEAR-1))
            AHUCGWS = AHUCGWS + ANGWS(I,J-(BYEAR-1))
            AHUCSWS = AHUCSWS + ANSWS(I,J-(BYEAR-1))
            AHUCDMSTOR=AHUCDMSTOR+ANDMSTOR(I,J-(BYEAR-1))
           
            IY = IY + 1
            
75       CONTINUE

         DO 101 R=1,12
      
            MOAVGSUP(I,R) = MOAVGSUP(I,R)/IY
            MOAVGGDEM(I,R) = MOAVGGDEM(I,R)/IY
            MOAVGNDEM(I,R) = MOAVGNDEM(I,R)/IY
            MOAVGFLW(I,R) = MOAVGFLW(I,R)/IY
            MOAVGSWS(I,R) = MOAVGSWS(I,R)/IY
            MOAVGGWS(I,R) = MOAVGGWS(I,R)/IY
            MOAVGDMST(I,R) = MOAVGDMST(I,R)/IY
            
            MOAVGWAS(I,R) = MOAVGNDEM(I,R)/MOAVGSUP(I,R)
     
101      CONTINUE 

         AHUCSUPPLY = AHUCSUPPLY / IY
         AHUCGROSDEM = AHUCGROSDEM / IY
         AHUCNETDEM = AHUCNETDEM / IY
         AHUCFLOW = AHUCFLOW / IY
         AHUCSWS = AHUCSWS/IY
         AHUCGWS = AHUCGWS/IY
         AHUCDMSTOR=AHUCDMSTOR/IY

C----FOR HUCS BORDERING THE GREAT LAKES, SET AHUCWASSI TO 0.0

                  IF(HUCNO(I).EQ.4010101.OR.HUCNO(I).EQ.4010301
     &            .OR.HUCNO(I).EQ.4010302.OR.HUCNO(I).EQ.4020101
     &            .OR.HUCNO(I).EQ.4020102.OR.HUCNO(I).EQ.4020103
     &            .OR.HUCNO(I).EQ.4020105.OR.HUCNO(I).EQ.4020201
     &            .OR.HUCNO(I).EQ.4020202.OR.HUCNO(I).EQ.4020203
     &            .OR.HUCNO(I).EQ.4020300.OR.HUCNO(I).EQ.4030101
     &            .OR.HUCNO(I).EQ.4030102.OR.HUCNO(I).EQ.4030103
     &            .OR.HUCNO(I).EQ.4030104.OR.HUCNO(I).EQ.4030105
     &            .OR.HUCNO(I).EQ.4030108.OR.HUCNO(I).EQ.4030109
     &            .OR.HUCNO(I).EQ.4030110.OR.HUCNO(I).EQ.4030111
     &            .OR.HUCNO(I).EQ.4030112.OR.HUCNO(I).EQ.4030204
     &            .OR.HUCNO(I).EQ.4040001.OR.HUCNO(I).EQ.4040002
     &            .OR.HUCNO(I).EQ.4040003.OR.HUCNO(I).EQ.4050002
     &            .OR.HUCNO(I).EQ.4050006.OR.HUCNO(I).EQ.4060101
     &            .OR.HUCNO(I).EQ.4060103.OR.HUCNO(I).EQ.4060104
     &            .OR.HUCNO(I).EQ.4060105.OR.HUCNO(I).EQ.4060106
     &            .OR.HUCNO(I).EQ.4060107.OR.HUCNO(I).EQ.4070001
     &            .OR.HUCNO(I).EQ.4070002.OR.HUCNO(I).EQ.4070003
     &            .OR.HUCNO(I).EQ.4070004.OR.HUCNO(I).EQ.4070007
     &            .OR.HUCNO(I).EQ.4080101.OR.HUCNO(I).EQ.4080102
     &            .OR.HUCNO(I).EQ.4080103.OR.HUCNO(I).EQ.4080104
     &            .OR.HUCNO(I).EQ.4080206.OR.HUCNO(I).EQ.4090001
     &            .OR.HUCNO(I).EQ.4090003.OR.HUCNO(I).EQ.4120102
     &            .OR.HUCNO(I).EQ.4120103.OR.HUCNO(I).EQ.4120104
     &            .OR.HUCNO(I).EQ.4130001.OR.HUCNO(I).EQ.4130003
     &            .OR.HUCNO(I).EQ.4140101.OR.HUCNO(I).EQ.4140102
     &            .OR.HUCNO(I).EQ.4140203.OR.HUCNO(I).EQ.4150102
     &            .OR.HUCNO(I).EQ.4150301.OR.HUCNO(I).EQ.4150302
     &            .OR.HUCNO(I).EQ.7120001.OR.HUCNO(I).EQ.7120003
     &            .OR.HUCNO(I).EQ.7120004.OR.HUCNO(I).EQ.7120005
     &            .OR.HUCNO(I).EQ.7120006.OR.HUCNO(I).EQ.4110004
     &            .OR.HUCNO(I).EQ.4110003.OR.HUCNO(I).EQ.4100010
     &            .OR.HUCNO(I).EQ.4100009.OR.HUCNO(I).EQ.4050001
     &            .OR.HUCNO(I).EQ.4110002.OR.HUCNO(I).EQ.4120101
     &            .OR.HUCNO(I).EQ.4090002.OR.HUCNO(I).EQ.4100001
     &            .OR.HUCNO(I).EQ.4110001.OR.HUCNO(I).EQ.4090004
     &            .OR.HUCNO(I).EQ.4100011.OR.HUCNO(I).EQ.4100012
     &            .OR.HUCNO(I).EQ.4090005.OR.HUCNO(I).EQ.4100002
     &            .OR.HUCNO(I).EQ.4050003) THEN
     
                  AHUCWASSI=0.0
                  
                  ELSE
                  
C----FOR ALL OTHER HUCS, CALCULATE AHUCWASSI

                  AHUCWASSI = AHUCWASSI / IY
                  
                  ENDIF

      DO 252 R=1,12
      
      MNDAY=MONTHD(R)
      
      CONSUSE = (-MOAVGNDEM(I,R)+MOAVGGWS(I,R))*MNDAY/264.17
      
      WRITE(804,253)HUCNO(I),R,MOAVGSWS(I,R),MOAVGGWS(I,R),
     &MOAVGSUP(I,R),MOAVGGDEM(I,R),MOAVGNDEM(I,R),MOAVGWAS(I,R),CONSUSE,
     &MOAVGFLW(I,R),MOAVGDMST(I,R)

253   FORMAT (I10,',',I10,',',F16.2,',',F16.2,',',F16.2,',',F16.2,',',
     &F16.2,',',F10.4,',',F16.2,',',F10.1,',',F16.2)       
     
252   CONTINUE

C -- WRITE THE AVERAGE ANNUAL WASSI, SUPPLY, AND DEMAND BETWEEN YEARS
C -- IYSTART AND IYEND FOR EACH CELL TO AVGANNWaSSI.TXT
       

C----ONLY WRITE OUTPUT FOR YEARS AFTER THE MODEL WARMUP PERIOD

      IF(J.GT.(IYSTART-BYEAR)) THEN           

         CONSUSE = (-AHUCNETDEM+AHUCGWS)*NDAY/264.17

         WRITE (300, 8000) HUCNO(I), IY, AHUCSWS,AHUCGWS,AHUCSUPPLY,
     &AHUCGROSDEM,AHUCNETDEM, AHUCWASSI, CONSUSE, AHUCFLOW,AHUCDMSTOR
                                              
8000     FORMAT (I10,',',I10,',',F16.2,',',F16.2,',',F16.2,',',F16.2,
     &',',F16.2,',',F10.4,',',F16.2,',',F10.1,',',F10.1)

      ENDIF

600   CONTINUE

      RETURN
      END
