C HFD, 20230329, fixed duplicated DIMENSION attributes C HFD, 20230330, 'F' declared as integer 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 HFD, 20230410: MOAVG* variables now explicitly initialized C HFD, 20230410: changed order of COMMON/CELLINFO/ vars (integers to the end) C This is to avoid potential problems when compiling with C -fdefault-real-8 flag C HFD, 20230413: creation of COMMON/MOISTURE/ (soil moisture init bug fix) 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) IMPLICIT NONE 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,200,10),LATUDE(4000),LONGI(4000), & TAREA(4000),HUCNO(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-----debug------------------------------------------------------------- COMMON/MOISTURE/UZTWC(10),UZFWC(10),LZTWC(10),LZFSC(10), & LZFPC(10),SSUR(10) c-----debug------------------------------------------------------------- C ---------------------------------------------------------------------------- INTEGER HUCNO, MNDAY INTEGER J, IYSTART, BYEAR REAL HUCELE, HUCAREA, LATUDE REAL RAIN, TEMP REAL AETTEMP, RUNOFFTEMP, TBFTEMP, > IFTEMP, GEPTEMP, RECOTEMP, NEETEMP REAL LZTWM, LZFPM, LZFSM, LZSK, > LZPK, UZTWM, UZFWM, UZK, ZPERC, > REXP, PFREE ! debug REAL UZTWC(10), UZFWC(10), LZTWC(10), LZFSC(10), LZFPC(10) REAL UZTWC, UZFWC, LZTWC, LZFSC, LZFPC !debug 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 REAL LADUSE, TAREA ! debug REAL ET(200,12,10), SSUR(10), GEP(200,12,10), INFIL(10), REAL ET(200,12,10), SSUR, GEP(200,12,10), INFIL(10), > RECO(200, 12, 10), NEE(200, 12,10) !debug REAL PET 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, RUNOFF(12), TOTBF(12), INTF(12), &SMC REAL AVUZTWC, AVUZFWC, AVLZTWC, AVLZFPC, > AVLZFSC REAL STRFLOW,GEPM,RECOM, > NEEM 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,MOAVGPET,MOAVGAET, &MOAVGPPT,MOAVGSNW,MOAVGSOL, &MOAVGTEM REAL IMPERV INTEGER F !Previously undeclared INTEGER I, M, IYEND, NGRID, NLC, NYEAR, POP_FLAG, K REAL ANUAET, ANUMAXSWE, ANUPET, ANURAIN, ANURUN, ANUSOIL, & ANUTEMP, APAET, APET, PAET, LONGI, RO 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--debug---------------------------------------------------------------- c if(I.eq.658) c & write(*,*) c & J,M,UZTWC(1),LZTWC(1),RAIN(I,J,M) c--debug---------------------------------------------------------------- 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,J,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,J,K)* &(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,J,8)) RECOTEMP = RECOTEMP + RECO(J,M,K)*(LADUSE(I,J,K)* &(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,J,8)) NEETEMP = NEETEMP + NEE(J,M,K)*(LADUSE(I,J,K)* &(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,J,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,J,K)* & (1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,J,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,J,K)* & (1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,J,8)) C--- SURFACE RUNOFF FROM IMPERVIOUS AREAS RUNOFFTEMP=RUNOFFTEMP+INFIL(K)*(LADUSE(I,J,K)* & (IMPERV(I,K)))/(TAREA(I)-LADUSE(I,J,8)) C--- INTERFLOW FROM NON-IMPERVIOUS AREAS IFTEMP = IFTEMP + SIF(K) *(LADUSE(I,J,K)* & (1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,J,8)) C--- TOTAL BASEFLOW FROM NON-IMPERVIOUS AREAS TBFTEMP = TBFTEMP + SBF(K)*(LADUSE(I,J,K)* & (1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,J,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,J,K)* &(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,J,8)) TAUZFWC = TAUZFWC + UZFWC(K) * (LADUSE(I,J,K)* &(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,J,8)) TALZTWC = TALZTWC + LZTWC(K) * (LADUSE(I,J,K)* &(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,J,8)) TALZFPC = TALZFPC + LZFPC(K) * (LADUSE(I,J,K)* &(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,J,8)) TALZFSC = TALZFSC + LZFSC(K) * (LADUSE(I,J,K)* &(1-IMPERV(I,K)))/(TAREA(I)-LADUSE(I,J,8)) ASM = ASM + (UZTWC(K)+UZFWC(K)+LZTWC(K)+LZFPC(K) &+LZFSC(K)) * (LADUSE(I,J,K)*(1-IMPERV(I,K)))/ &(TAREA(I)-LADUSE(I,J,8)) ELSEIF(REGION.EQ.'US'.AND.LADUSE(I,J,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,J,K)/TAREA(I) RECOTEMP = RECOTEMP + RECO(J,M,K)*LADUSE(I,J,K)/TAREA(I) NEETEMP = NEETEMP + NEE(J,M,K)*LADUSE(I,J,K)/TAREA(I) AETTEMP = AETTEMP + ET(J,M,K) * LADUSE(I,J,K)/TAREA(I) RUNOFFTEMP = RUNOFFTEMP + SSUR(K) * LADUSE(I,J,K)/TAREA(I) TBFTEMP = TBFTEMP + SBF(K) * LADUSE(I,J,K)/TAREA(I) IFTEMP = IFTEMP + SIF(K) * LADUSE(I,J,K)/TAREA(I) TAUZTWC = TAUZTWC + UZTWC(K) * LADUSE(I,J,K)/TAREA(I) TAUZFWC = TAUZFWC + UZFWC(K) * LADUSE(I,J,K)/TAREA(I) TALZTWC = TALZTWC + LZTWC(K) * LADUSE(I,J,K)/TAREA(I) TALZFPC = TALZFPC + LZFPC(K) * LADUSE(I,J,K)/TAREA(I) TALZFSC = TALZFSC + LZFSC(K) * LADUSE(I,J,K)/TAREA(I) ASM = ASM + (UZTWC(K)+UZFWC(K)+LZTWC(K)+LZFPC(K) &+LZFSC(K)) * LADUSE(I,J,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,J,8))/TAREA(I) + &ET(J,M,8)*LADUSE(I,J,8)/TAREA(I) RUNOFF(M)=RUNOFF(M)*(TAREA(I)-LADUSE(I,J,8))/TAREA(I)+ &SSUR(8)*LADUSE(I,J,8)/TAREA(I) TOTBF(M)=TOTBF(M)*(TAREA(I)-LADUSE(I,J,8))/TAREA(I) INTF(M)=INTF(M)*(TAREA(I)-LADUSE(I,J,8))/TAREA(I) AVUZTWC(M)=AVUZTWC(M)*(TAREA(I)-LADUSE(I,J,8))/TAREA(I) AVUZFWC(M)=AVUZFWC(M)*(TAREA(I)-LADUSE(I,J,8))/TAREA(I) AVLZTWC(M)=AVLZTWC(M)*(TAREA(I)-LADUSE(I,J,8))/TAREA(I) AVLZFPC(M)=AVLZFPC(M)*(TAREA(I)-LADUSE(I,J,8))/TAREA(I) AVLZFSC(M)=AVLZFSC(M)*(TAREA(I)-LADUSE(I,J,8))/TAREA(I) SMC (M)=SMC (M)*(TAREA(I)-LADUSE(I,J,8))/TAREA(I) GEPM(I,J,M) = GEPM(I,J,M)*(TAREA(I)-LADUSE(I,J,8))/TAREA(I) + &GEP(J,M,8)*LADUSE(I,J,8)/TAREA(I) RECOM(I,J,M) = RECOM(I,J,M)*(TAREA(I)-LADUSE(I,J,8))/TAREA(I) + &RECO(J,M,8)*LADUSE(I,J,8)/TAREA(I) NEEM(I,J,M) = NEEM(I,J,M)*(TAREA(I)-LADUSE(I,J,8))/TAREA(I) + &NEE(J,M,8)*LADUSE(I,J,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--debug---------------------------------------------------------------- c if(I.eq.658 .and. J.eq.4) c & write(*,*) M,temp(I,J,M),snowppt,rainppt,snowpack,snoww,smf c--debug---------------------------------------------------------------- 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 !Initialization, previously absent IF(J.EQ.(IYSTART-BYEAR+1)) THEN MOAVGRO(I,M) =0. MOAVGPET(I,M) =0. MOAVGAET(I,M) =0. MOAVGPPT(I,M) =0. MOAVGTEM(I,M) =0. MOAVGSNW(I,M) =0. MOAVGSOL(I,M) =0. ENDIF 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