C    HFD, 20230329, fixed duplicated DIMENSION attributes
C    HFD, 20230405: code adapted to process dynamic land cover and LAI
C                   At this point, only US region, US_XX clim scenario are supported
C    HFD, 20230407: all variables are now explicitly declared

C**********************************************************************C
C                                                                      C
C     *** SUBROUTINE FLOWROUTING ***                                   C
C     SIMULATE TOTAL FLOW FOR EACH MONTH AND EACH HUC                  C
C**********************************************************************C
C        I=HUC; J= YEAR, M =MONTH, MNDAY= NO OF DAYS IN A MONTH

      SUBROUTINE FLOWROUTING_RW
      
      IMPLICIT NONE
      
      COMMON/BASIC/NGRID,NYEAR,NLC,BYEAR,IYSTART,IYEND,POP_FLAG,
     & REGION,CLIMATE_FLAG

      COMMON/FLOW/ STRFLOW(2200, 200, 12)        
         
     
  
C --------------------------------------------------------------
      INTEGER NGRID,NYEAR,NLC,BYEAR,IYSTART,IYEND,POP_FLAG
      
      REAL STRFLOW
      
      CHARACTER*4 HEADER(20)
      
      INTEGER MAX_LVL,A
      
      INTEGER F_WS(1000),T_WS(1000), T_WS_LVL(1000)
      
      INTEGER I,J,M,L,R
      
      INTEGER LEVEL
      
      INTEGER IDY,IY
      
      REAL ANFLOW(2200,200),AHUCFLOW,MOAVGFLW(2200,12)

      CHARACTER*5 CLIMATE_FLAG
      
      CHARACTER*2 REGION


C --------------------------------------------------------------

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 TABLE, DETERMINE NUMBER OF NODES AND MAX LEVEL
     
      READ(17,1) HEADER
1     FORMAT(20A4)
      
      MAX_LVL=0.
      
      DO 2 A=1,1000000
           
      READ(17,*,END=3) F_WS(A),T_WS(A), T_WS_LVL(A)
      
      IF(T_WS_LVL(A).GT.MAX_LVL) THEN
      
      MAX_LVL=T_WS_LVL(A)
      
      ENDIF
      
2     CONTINUE

3     A=A-1
     

C************************************************************************
C --  PERFORM FLOW ROUTING CALCULATIONS 
C************************************************************************
     
      DO 200 J=(IYSTART-BYEAR+1),(IYEND-BYEAR+1)
      
      DO 100 M=1, 12

C -- STARTING AT HIGHEST TO-HUC LEVEL, LOOP THROUGH LEVELS FOR FLOW ACCUMULATION

      LEVEL=MAX_LVL
      
      DO 150 L=1,MAX_LVL
      
C -- LOOK THROUGH ROUTING TABLE, IF TO-HUC AT LEVEL IS AT LEVEL
C -- L, ADD THE FROM-HUC FLOW TO FLOW AT OUTLET OF TO HUC

      DO 400 R=1,A
      
      IF (T_WS_LVL(R).EQ.LEVEL) THEN
      
      STRFLOW(T_WS(R),J,M)=STRFLOW(T_WS(R),J,M)+STRFLOW(F_WS(R),J,M)
      
      ENDIF
      
400   CONTINUE

      LEVEL=MAX_LVL-L
      
150   CONTINUE

100   CONTINUE 

200   CONTINUE


C************************************************************************
C ---WRITE OUTPUT TO OUTPUT FILES
C************************************************************************

      DO 700 I=1, NGRID  
      
      DO 199 M=1,12
         
      MOAVGFLW(I,M)=0.
         
199   CONTINUE
       
                
      DO 600 J=(IYSTART-BYEAR+1),(IYEND-BYEAR+1)
         
      IDY = J + BYEAR - 1 
      
      ANFLOW(I,J) = 0.
            
      DO 500 M=1, 12
           
C----ONLY WRITE OUTPUT AND CALCULATE MONTHLY MEANS FOR YEARS AFTER THE MODEL WARMUP PERIOD

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

C ---WRITE MONTHLY FLOW MONTHWASSI.TXT
              
      WRITE (100, 2000) I,IDY,M,STRFLOW(I,J,M)
     
2000  FORMAT (I10,',',I10,',',I10,
     &',0.0,0.0,0.0,0.0,0.0,0.0,0.0,',F16.2,',0.0') 
     
      ENDIF
     
      ANFLOW(I,J) = ANFLOW(I,J) + STRFLOW(I,J,M)

      MOAVGFLW(I,M) = MOAVGFLW(I,M) + STRFLOW(I,J,M)

500   CONTINUE

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

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

C -- WRITE  ANNUAL FLOW FOR EACH YEAR TO ANNUALFLOW.TXT
       
      WRITE (200, 7000) I, IDY,ANFLOW(I,J)
                                              
7000  FORMAT (I10,',',I10,',0.0,0.0,0.0,0.0,0.0,0.0,0.0,',F16.2,',0.0')

      ENDIF

600   CONTINUE

      AHUCFLOW = 0.
         
      IY = 0
         
      DO 75 J = IYSTART,IYEND
         
      IF(J.GT.(IYSTART-BYEAR)) THEN
      
      AHUCFLOW = AHUCFLOW + ANFLOW(I,J-(BYEAR-1))
           
      IY = IY + 1
      
      ENDIF
            
75    CONTINUE

      DO 101 R=1,12
      
      MOAVGFLW(I,R) = MOAVGFLW(I,R)/IY
      
      WRITE (804,253)I,R,MOAVGFLW(I,R)
      
253   FORMAT(I10,',',I10,',0.0,0.0,0.0,0.0,0.0,0.0,0.0,',F16.2,',0.0')
     
101   CONTINUE 

      AHUCFLOW = AHUCFLOW / IY

C -- WRITE THE AVERAGE ANNUAL FLOW BETWEEN YEARS
C -- IYSTART AND IYEND FOR EACH CELL TO AVGANNFLOW.TXT
       

      WRITE (300, 8000) I, IY, AHUCFLOW
                                              
8000  FORMAT (I10,',',I10,',0.0,0.0,0.0,0.0,0.0,0.0,0.0,',F16.2,',0.0')


700   CONTINUE

      RETURN
      END