! source file: /Users/nmengis/Documents/UVic_ESCM/2.10/updates/02_CE_permafrost_merge_CMIP6forcing/source/mtlm/penmon.F subroutine PENMON(POINTS, LAND_PTS, LAND_INDEX, NGND, DT, RHO_S &, DZ, HCON_SOIL, RSU, RSF, Z0, LW, SWN, PSTAR &, Q1, T1, TGND, WIND, LHC, LHF, EPSILON &, SIGMA, ZERODEGC, ET, SUBSNOW, SUBSOIL &, LE, SH, G, TSTAR, LW_OUT, SU, SF &, MH, SNOW_HT, HT, HCON_SNOW) !----------------------------------------------------------------------- ! Routine to calculate the evaporation using an extended version of ! the Penman-Monteith equation. !********************************************************************** ! this file is based on code that may have had the following copyright: ! (c) CROWN COPYRIGHT 1997, U.K. METEOROLOGICAL OFFICE. ! Permission has been granted by the authors to the public to copy ! and use this software without charge, provided that this Notice and ! any statement of authorship are reproduced on all copies. Neither the ! Crown nor the U.K. Meteorological Office makes any warranty, express ! or implied, or assumes any liability or responsibility for the use of ! this software. !********************************************************************** !----------------------------------------------------------------------- implicit none !~ Constants ~ ! CP = Specific heat of dry air at constant pressure (J/kg/K) ! DT = Timestep (s) ! EPSILON = Ratio of molecular weights of water and dry air. ! HCON_SNOW = Thermal conductivity of snow (W/m/K) ! KARMAN = Von Karman's constant ! KARMAN_SQ = Square of Von Karman's constant ! LHC = Latent heat of condensation (J/kg). ! LHF = Latent heat of fusion (J/kg). ! R = Gas constant, J/kg/K ! RHO_S = Density of lying snow (Kg/m3) ! SIGMA = Stefan-Boltzman constant (W/m2/K4). ! Z1 = Reference height for determening turbulent fluxes (m). ! ZERODEGC = Zero Celsius (K). real EPSILON, LHC, LHF, SIGMA, Z1, ZERODEGC real R, CP, KARMAN, KARMAN_SQ, HCON_SNOW, DT, RHO_S parameter (R=287.05, CP=1005.0, KARMAN=0.4, KARMAN_SQ=0.16 ) ! 1 change made here to make HCON_SNOW general variable parameter (Z1 = 10.) !~ Spatial indexing variables ~ ! I, L = Loop counters ! LAND_PTS = Number of points on which TRIFFID may operate. ! LAND_INDEX = Indices of land points on which TRIFFID may operate. ! POINTS = Total number of land points. ! NGND = Number of ground levels integer POINTS, LAND_PTS, LAND_INDEX(POINTS), I, L, NGND !~ Soil variables~ ! DZ = Ground layer thickness (m). ! HCON_SOIL = Soil surface thermal conductivity (W/m/K) ! SF = Soil frozen moisture content as a fraction of saturation ! SU = Soil unfrozen moisture content as a fraction of saturation ! TGND = Sub-surface temperature (K). real DZ(NGND) real SU(POINTS, NGND), SF(POINTS, NGND) real TGND(POINTS, NGND), HCON_SOIL(POINTS) !~Surface characteristics~ ! HT = Soil / vegetation height (m) ! LW = Downward LW (W/m2). ! LYING_SNOW = Snow (kg/m2). ! PSTAR = Surface pressure (Pa). ! Q1 = Specific humidity (kg/kg). ! RSU = Surface resistance - evaporation (s/m). ! RSF = Surface resistance - sublimation (s/m) ! SNOW_HT = Height of lying snow (m) ! SWN = Net downward SW (W/m2). ! T1 = Atmospheric temperature (K). ! TSTAR = Skin temperature (L) ! WIND = Windspeed (m/s). ! Z0 = Roughness length (m). real LW(POINTS), LYING_SNOW(POINTS) real PSTAR(POINTS), Q1(POINTS), RSU(POINTS), RSF(POINTS) real SWN(POINTS), SNOW_HT(POINTS), TSTAR(POINTS), HT(POINTS) real T1(POINTS), WIND(POINTS), Z0(POINTS) !~ Flux Outputs~ ! ET = Evapotranspiration (kg/m2/s). ! G = Ground heat flux (W/M2). ! LE = Latent heat flux (W/m2). ! LW_OUT = Net longwave radiation (W/m2) ! MH = Energy to melt lying snow from below (W/M2) ! RADNET = Surface net radiation (W/m2). ! SH = Sensible heat flux (W/m2). ! SUBSOIL = Sublimation from bare soil (kg/m2/s) ! SUBSNOW = Sublimation from soil (kg/m2/s) ! TSTAR = Skin temperature (K). real ET(POINTS), G(POINTS), LE(POINTS), LW_OUT(POINTS) real SUBSOIL(POINTS), SUBSNOW(POINTS) real RADNET(POINTS), SH(POINTS), MH(POINTS) !~ Other variables ~ ! AS1 = Thermal conductivity divided by DZsoil+DZsnow ! C1, C2, C3 = Switches to control sublimation / evaporation ! CHN = Neutral transfer coefficients. ! dENdT = Derivative of the energy balance eqn. w.r.t. temp (W/m2T) ! DQ1 = Humidity deficit (kg/kg). ! DQS_DT = Rate of change of saturated specific humidity with ! temperature (kg/kg/K). ! ENBAL = Surface energy budget balance (in theory should be 0) (W/m2) ! LAT = Latent heat constant (J/kg). ! MELTFLUX = Flux to melt all snow (W/m2) ! RESF = 1/(1+RS/RA). ! QS = Saturated specific humidity at (T1,PSTAR) (kg/kg). ! RA = Aerodynamic resistance (s/m). ! RHOSTAR = Surface air density (kg/m3). ! ZETAH,ZETAM = Temporaries in calculation of CHN. real AS1(POINTS), CHN real DQ1(POINTS), DQS_DT(POINTS), DUM, LAT(POINTS) real RESF, QS(POINTS), RA(POINTS), RHOSTAR(POINTS) real ZETAH, ZETAM, MELTFLUX real C1(POINTS), C2(POINTS), C3(POINTS), QS2(POINTS) real ENBAL(POINTS), dEBdT(POINTS), TSTAR0 !--------------------------------------------------------------------------------- ! Initializations - set up all the variables for solving the surface energy balance !--------------------------------------------------------------------------------- ! Calculate the saturated specific humidity call QSAT (POINTS, LAND_PTS, LAND_INDEX, EPSILON, ZERODEGC &, QS, T1, PSTAR) do I=1,LAND_PTS L = LAND_INDEX(I) !---------------------------------------------------------------------- ! First figure out what kind of calculation penmon must do: ! a) Snow covered terrain (no visible veg or soil) ! b) Vegetation w/ lying snow ! c) Vegetation (transpiration) or bare soil (evaporation) ! d) Bare soil (sublimation) !--------------------------------------------------------------------- if (SNOW_HT(L) .gt. 0.) then C2(L) = 0. C3(L) = 1. if (SNOW_HT(L) .gt. HT(L)) then ! Case a) Snow only C1(L) = 0. else ! Case b) Snow + veg C1(L) = 1. endif elseif ((HT(L).gt. 0.)) then ! Case c) Vegetation C1(L) = 1. C2(L) = 0. C3(L) = 0. else ! Case d) Bare soil C1(L) = SU(L,1) / (SU(L,1) + SF(L,1)) C2(L) = SF(L,1) / (SU(L,1) + SF(L,1)) C3(L) = 0. if (SU(L,1) .lt. 0.0001) C1(L) = 0. if (SF(L,1) .lt. 0.0001) C2(L) = 0. endif ! Calculate the surface air density and humidity deficit RHOSTAR(L) = PSTAR(L)/(R*T1(L)) DQ1(L) = QS(L) - Q1(L) ! Calculate the thermal conductivity of the soil, divided thickness of soil if (SNOW_HT(L) .gt. 0.) then if (SNOW_HT(L) .gt. (DZ(1)/2)) then AS1(L) = 2*HCON_SNOW/DZ(1) else AS1(L) = 2*HCON_SOIL(L)/(DZ(1)+2*SNOW_HT(L) & *(HCON_SOIL(L)/HCON_SNOW-1)) endif else AS1(L) = 2*HCON_SOIL(L)/DZ(1) endif ! Calculate the neutral bulk transfer coefficient and aerodynamic resistance. ZETAM = LOG((Z1 + Z0(L)) / Z0(L)) ZETAH = LOG((Z1 + Z0(L)) / (0.1*Z0(L))) CHN = KARMAN_SQ / (ZETAH * ZETAM) RA(L) = 1.0 / (CHN * WIND(L)) !---------------------------------------------------------------------- ! Now start Newton's method to solve the surface energy budget and T* !---------------------------------------------------------------------- ! Evaluate the surface energy budget at T* = T1 ENBAL(L) = SWN(L) + LW(L) - SIGMA*(T1(L)*T1(L)*T1(L)*T1(L)) & - (C1(L) * LHC /(RA(L)+RSU(L)) & + C2(L) * (LHC+LHF)/(RA(L)+RSF(L)) & + C3(L) * (LHC+LHF)/RA(L)) * RHOSTAR(L) * DQ1(L) & - AS1(L)*(T1(L)-TGND(L,1)) ! Now evaluate the derivative of the energy budget with respect to temperature (evaluated at T1) DQS_DT(L) = (EPSILON*LHC*QS(L))/(R*T1(L)*T1(L)) dEBdT(L) = -4*SIGMA*(T1(L)*T1(L)*T1(L))-RHOSTAR(L)*CP/RA(L) & - AS1(L) -(C1(L) * LHC/(RA(L)+RSU(L)) & + C2(L) *((LHC+LHF)*(LHC+LHF))/LHC/(RA(L)+RSF(L)) & + C3(L)*((LHC+LHF)*(LHC+LHF))/LHC/RA(L))*RHOSTAR(L) & * DQS_DT(L) ! T = To - f(To)/f'(To) TSTAR(L) = T1(L) - ENBAL(L)/dEBdT(L) RHOSTAR(L) = PSTAR(L)/(R*TSTAR(L)) ! Reduce TSTAR(L) to 0 oC if there's lying snow to produce snowmelt if ((SNOW_HT(L) .gt. 0) .and. (TSTAR(L) .gt. ZERODEGC)) then MELTFLUX = SNOW_HT(L) * RHO_S *LHF / DT TSTAR0 = TSTAR(L) TSTAR(L) = MAX(ZERODEGC, TSTAR(L) + MELTFLUX/dEBdT(L)) G(L) = AS1(L)*(TSTAR(L)-TGND(L,1)) SH(L) = RHOSTAR(L)*CP/RA(L)*(TSTAR(L)-T1(L)) LW_OUT(L) = LW(L) - SIGMA*TSTAR(L)*TSTAR(L)*TSTAR(L)*TSTAR(L) RADNET(L) = SWN(L) + LW_OUT(L) !Calculat QS at the new TSTAR call QSAT (1, 1, 1, EPSILON,ZERODEGC,QS(L),TSTAR(L),PSTAR(L)) ET(L) = C1(L) * RHOSTAR(L) / (RA(L)+RSU(L))*(QS(L)-Q1(L)) SUBSOIL(L) = C2(L) * RHOSTAR(L)/(RA(L)+ RSF(L))*(QS(L)-Q1(L)) SUBSNOW(L) = C3(L) * RHOSTAR(L) / RA(L) * (QS(L) - Q1(L)) LE(L) = LHC * ET(L) + (LHC+LHF)*(SUBSOIL(L) + SUBSNOW(L)) MH(L) = RADNET(L)-LE(L)-SH(L)-G(L) if (MH(L) .lt. 0.) then ! Energy balance strongly non-linear. Go back to first guess MH(L) = 0. TSTAR(L) = TSTAR0 SH(L) = RHOSTAR(L)*CP/RA(L)*(TSTAR(L)-T1(L)) LW_OUT(L) =LW(L)-SIGMA*TSTAR(L)*TSTAR(L)*TSTAR(L)*TSTAR(L) RADNET(L) = SWN(L) + LW_OUT(L) !Calculate QS at the new TSTAR call QSAT (1,1,1,EPSILON,ZERODEGC,QS(L),TSTAR(L),PSTAR(L)) ET(L) = C1(L) * RHOSTAR(L) / (RA(L)+RSU(L))*(QS(L)-Q1(L)) SUBSOIL(L)=C2(L)*RHOSTAR(L)/(RA(L)+ RSF(L))*(QS(L)-Q1(L)) SUBSNOW(L)=C3(L)*RHOSTAR(L) / RA(L) * (QS(L) - Q1(L)) LE(L) = LHC * ET(L) + (LHC+LHF)*(SUBSOIL(L) + SUBSNOW(L)) G(L) = RADNET(L) - LE(L) - SH(L) endif else SH(L) = RHOSTAR(L)*CP/RA(L)*(TSTAR(L)-T1(L)) LW_OUT(L) = LW(L) - SIGMA*TSTAR(L)*TSTAR(L)*TSTAR(L)*TSTAR(L) RADNET(L) = SWN(L) + LW_OUT(L) !Calculae QS at the new TSTAR call QSAT (1, 1, 1,EPSILON,ZERODEGC,QS(L),TSTAR(L),PSTAR(L)) ET(L) = C1(L) * RHOSTAR(L) / (RA(L)+RSU(L))*(QS(L)-Q1(L)) SUBSOIL(L) = C2(L) * RHOSTAR(L)/(RA(L)+ RSF(L))*(QS(L)-Q1(L)) SUBSNOW(L) = C3(L) * RHOSTAR(L) / RA(L) * (QS(L) - Q1(L)) LE(L) = LHC * ET(L) + (LHC+LHF)*(SUBSOIL(L) + SUBSNOW(L)) G(L) = RADNET(L) - LE(L) - SH(L) MH(L) = 0. endif enddo return end