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 HFD, 20230410: changed order of COMMON/LANDCHANGE/ variables (integer var to the end) C This allows compilation with -fdefault-real-8 flag C Did the same for the COMMON/CELLINFO/ block C HFD, 20230714: IMPORTANT CORE CHANGE: PAET calculation for NE hucs was changed C**********************************************************************C C C C *** SUBROUTINE WARMPET *** C C INPUT MONTHLY PRECIPITATION AND TEMPERATURE, AND CALCULATE C C MONTHLY PET AND POTENTIAL AET C C C C**********************************************************************C SUBROUTINE WARMPET(I, J, M, MNDAY) IMPLICIT NONE COMMON/BASIC/NGRID,NYEAR,NLC,BYEAR,IYSTART,IYEND,POP_FLAG, & REGION,CLIMATE_FLAG COMMON/CELLINFO/LADUSE(4000,200,10),LATUDE(4000),LONGI(4000), & TAREA(4000),HUCNO(4000) COMMON/LANDCHANGE/FPERD, LAI_CHG, FOR_CONV_TO 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/CLIMATE/ RAIN(4000,200,12), TEMP(4000,200, 12) COMMON/LAI/LAI(4000,200,12,10) INTEGER NLC, HUCNO, M, BYEAR, WRR INTEGER MONTHD(12),MONTHL(12), MJD(12), MMD REAL DEGLAT, LATUDE, PET,TAREA REAL APET, PE, LADUSE, TEMP REAL TPET, RAIN, FOREST REAL FPERD REAL LAI REAL PAET, TPAET, APAET CHARACTER*5 CLIMATE_FLAG CHARACTER*2 REGION !Previously undeclared INTEGER I, J, MNDAY, FOR_CONV_TO, IYEND, IYSTART, NGRID, NYEAR, & POP_FLAG, K REAL AET, AVLZFPC, AVLZFSC, AVLZTWC, AVUZFWC, AVUZTWC, DTEMP, & LAI_CHG, LONGI, RO 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 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-----JULIAN DAY FOR MID-DAY OF EACH MONTH DATA MJD/15,46,76,107,137,168,198,229,259,290,321,351/ C --- Calculate Monthly potential evapotranspiration DTEMP=TEMP(I,J,M) TPET=0. TPAET=0. DO 40 K=1, NLC C -- MMD = JULIAN DATE FOR MONTH M MMD=MJD(M) DEGLAT = LATUDE(I) C -K= LANDUSE TYPE, DTEMP = AIR TEMP, HPEC = CORRECTION PARAMETER, PE =CALCUALTED PET (MM) CALL HAMON(K,DTEMP,M,MMD,DEGLAT,PE) C ----PET (YEAR, MONTH, LANDUSE) FOR THE CURRENT CELL) PET(J,M,K) = PE PET(J,M,K) = PET(J,M,K) * MNDAY C ----CALCULATE PAET C ----PAET (YEAR, MONTH, LANDUSE) FOR THE CURRENT CELL) C ----PAET IS THE POTENTIAL AET, ASSUMING SOIL MOISTURE NOT LIMITING IF (REGION.EQ.'US') THEN WRR=INT(HUCNO(I)/1000000) C --- COMPUTE THE TOTAL FOREST LAND COVER PERCENTAGE (PRIOR TO CHANGE BY USER) C --- THIS IS NECESSARY TO PREVENT USING BOTH EQUATIONS FOR THE SAME HUC C --- FOR ORIGINAL FOREST COVER AND REDUCED FOREST COVER IF REDUCING THE FOREST C --- COVER CAUSES FOREST TO GO FROM >0.2 TO <0.2 C -- SET PAET FOR OPEN WATER LAND COVER CLASS TO PET IF(K.EQ.8) THEN PAET(J,M,K)=PET(J,M,K) ELSE C --- COMPUTE THE TOTAL FOREST LAND COVER PERCENTAGE (PRIOR TO CHANGE BY USER) C --- THIS IS NECESSARY TO PREVENT USING BOTH EQUATIONS FOR THE SAME HUC C --- FOR ORIGINAL FOREST COVER AND REDUCED FOREST COVER IF REDUCING THE FOREST C --- COVER CAUSES FOREST TO GO FROM >0.2 TO <0.2 FOREST=(LADUSE(I,J,2)+LADUSE(I,J,3) & +LADUSE(I,J,4))/(1+FPERD) IF (WRR.EQ.1.OR.(WRR.EQ.2.AND.LATUDE(I).GT.40).OR. & WRR.EQ.4.OR.WRR.EQ.5) THEN C HFD: Removing forest conditional to avoid inconsistencies over time, when C forest (now dynamic) oscillates above and below 0.2. C We are still using the special eq for this NE region, which is dominated C by forest anyways. The difference now is that the especial eq (lower PAET) C will now be used for HUCs dominated by non-forest cover (>0.8) c HFD IF (FOREST.GT.0.2) THEN PAET(J,M,K) = 0.00169*PET(J,M,K)*RAIN(I,J,M) + & 0.4*PET(J,M,K)+7.83*LAI(I,J,M,K) c HFD ELSE c HFD c HFD PAET(J,M,K) = (PET(J,M,K))* c HFD & LAI(I,J,M,K)*0.0222+0.174* c HFD & RAIN(I,J,M)+0.502*(PET(J,M,K))+5.31*LAI(I,J,M,K) c HFD c HFD ENDIF ELSE PAET(J,M,K) = (PET(J,M,K))* & LAI(I,J,M,K)*0.0222+0.174* & RAIN(I,J,M)+0.502*(PET(J,M,K))+5.31*LAI(I,J,M,K) ENDIF ENDIF ELSEIF (REGION.EQ.'MX'.OR.REGION.EQ.'RW') THEN PAET(J,M,K) = (PET(J,M,K))* & LAI(I,J,M,K)*0.0222+0.174* & RAIN(I,J,M)+0.502*(PET(J,M,K))+5.31*LAI(I,J,M,K) ENDIF C ----CALCULATE TOTAL PET AND PAET FOR THE HUC FOR A GIVEN YEAR AND MONTH TPET = TPET + PET(J,M,K) * LADUSE(I,J,K)/TAREA(I) TPAET = TPAET + PAET(J,M,K) * LADUSE(I,J,K)/TAREA(I) 40 CONTINUE C ------APET =AVERAGE PET FOR CURRENT CELL, FOR ALL YEAR, MONTH APET(M) = TPET C ------APAET =AVERAGE PAET FOR CURRENT CELL, FOR ALL YEAR, MONTH APAET(M) = TPAET C ------TEST OUTPUT C WRITE(77, 5040) HUCNO(I), J, M, RAIN(I,J,M), APET(M), APAET(M) 5040 FORMAT (3I10, 3F20.1) C --- Return RETURN END C**********************************************************************C C C C *** SUBROUTINE HAMON *** C C Calculate potential evapotranspiration by using Hamon's C C equation C C C C**********************************************************************C SUBROUTINE HAMON(K, TEMP,MONTH,J,DEGLAT,PE) IMPLICIT NONE REAL PE, TEMP, DEGLAT REAL SOLDEC, SSANG, DAY REAL ESAT, RHOSAT, PI INTEGER J !Previously undeclared INTEGER K, MONTH PI = 3.14159265 C --- ESAT = saturated vapor pressure at TEMP C --- RHOSAT = saturated vapor density at TEMP ESAT = 6.108*EXP(17.2693882*TEMP/(TEMP+237.3)) RHOSAT = 216.7*ESAT/(TEMP+273.3) C --- CALCULATE MEAN MONTHLY DAY LENGTH (DAY) BASED ON LATITUDE AND MID-MONTH JULIAN DATE C --- SHUTTLEWORTH, W.J. 1993. Evaporation, in Handbook of Hydrology, C --- MAIDMENT, D.R. ED., MCGRAW HILL, NY, PP. 4, 1-53 C --- CALCULATE THE ANGLE OF SOLAR DECLINATION IN RADIANS SOLDEC = 0.4093*SIN(2*PI*J/365.-1.405) C --- CALCULATE THE SUNSET HOUR ANGLE IN RADIANS SSANG = ACOS(-1*TAN(DEGLAT*.0174529)*TAN(SOLDEC)) C --- CALCULATE THE ADJUSTMENT IN LENGTH OF DAY FOR 12 HR PERIOD DAY = 2*SSANG/PI C --- Calculate Daily PE PE = 0.1651*DAY*RHOSAT RETURN END