#if defined O_mtlm_nsoil 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) # if defined O_mtlm !----------------------------------------------------------------------- ! 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 ) ! O_mtlm_kovdiff 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 # endif return end #else subroutine PENMON (POINTS, LAND_PTS, LAND_INDEX, DZ_SOIL &, HCON_SOIL, RS, Z0, LW, SWN, PSTAR, Q1, T1 &, TS1, WIND, Z1, LC, LF, EPSILON, SIGMA, TM &, ZERODEGC, E, LE, SH, G, RADNET, TSTAR, MNEG &, LW_OUT, LYING_SNOW) # if defined O_mtlm !----------------------------------------------------------------------- ! 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 ! POINTS = IN Total number of land points. ! LAND_PTS = IN Number of points on which TRIFFID may operate. ! LAND_INDEX = IN Indices of land points on which TRIFFID may operate. integer POINTS, LAND_PTS, LAND_INDEX(POINTS), I, L ! Surface variables ! DZ_SOIL = IN Soil layer thickness (m). ! HCON_SOIL = IN Soil heat capacity (W/m/K). ! RS = IN Surface resistance (s/m). ! Z0 = IN Roughness length (m). ! MNEG = IN Negative soil moisture (kg/m2). ! LYING_SNOW = IN Snow (kg/m2). ! Driving variables ! LW = IN Downward LW (W/m2). ! SWN = IN Net downward SW (W/m2). ! PSTAR = IN Surface pressure (Pa). ! Q1 = IN Specific humidity (kg/kg). ! T1 = IN Atmospheric temperature (K). ! TS1 = IN Sub-surface temperature (K). ! WIND = IN Windspeed (m/s). ! Parameters ! Z1 = IN Reference height (m). ! LC = IN Latent heat of condensation (J/kg). ! LF = IN Latent heat of fusion (J/kg). ! EPSILON = IN Ratio of molecular weights of water and dry air. ! SIGMA = IN Stefan-Boltzman constant (W/m2/K4). ! TM = IN Melting point of fresh water (K). ! ZERODEGC = IN Zero Celsius (K). ! Outputs ! E = OUT Evapotranspiration (kg/m2/s). ! LE = OUT Latent heat flux (W/m2). ! SH = OUT Sensible heat flux (W/m2). ! G = OUT Ground heat flux (W/M2). ! RADNET = OUT Surface net radiation (W/m2). ! TSTAR = OUT Surface temperature (K). ! LW_OUT = OUT net longwave radiation (W/m2). ! Work Variables ! AHAT = WORK "Available energy" (W/m2). ! AS1 = WORK 2*HCON_SOIL/DZ_SOIL (W/m2/K). ! CHN = WORK Neutral transfer coefficients. ! DENOM = WORK Denominator of PM equation. ! DQ1 = WORK Humidity deficit (kg/kg). ! DQS_DT = WORK Rate of change of saturated specific humidity with ! temperature (kg/kg/K). ! DUM = WORK Workspace variable. ! LAT = WORK Latent heat constant (J/kg). ! NUMER = WORK Numerator of PM equation. ! RESF = WORK 1/(1+RS/RA). ! QS1 = WORK Saturated specific humidity at (T1,PSTAR) (kg/kg). ! RA = WORK Aerodynamic resistance (s/m). ! RHOSTAR = WORK Surface air density (kg/m3). ! ZETAH,ZETAM = WORK Temporaries in calculation of CHN. real DZ_SOIL ,HCON_SOIL, RS(POINTS), Z0(POINTS) real MNEG(POINTS), LYING_SNOW(POINTS), LW(POINTS) real SWN(POINTS), PSTAR(POINTS), Q1(POINTS), T1(POINTS) real TS1(POINTS), WIND(POINTS), Z1, LC, LF, EPSILON, SIGMA real TM, ZERODEGC, E(POINTS), LE(POINTS), SH(POINTS) real G(POINTS), RADNET(POINTS), TSTAR(POINTS) real LW_OUT(POINTS), AHAT(POINTS), AS1, CHN, DENOM real DQ1(POINTS), DQS_DT(POINTS), DUM, LAT(POINTS), NUMER real RESF, QS1(POINTS), RA(POINTS), RHOSTAR(POINTS), ZETAH real ZETAM ! Local parameters ! R = Gas constant (J/kg/K). ! CP = Specific heat of dry air at constant pressure (J/kg/K). ! KARMAN = Von Karman's constant. real R, CP, KARMAN,KARMAN_SQ parameter (R=287.05, CP=1005.0, KARMAN=0.4, KARMAN_SQ=0.16 ) AS1 = 2*HCON_SOIL/DZ_SOIL !---------------------------------------------------------------------- ! Calculate the surface air density (use atmospheric rather than ! surface temperature). !---------------------------------------------------------------------- do I=1,LAND_PTS L = LAND_INDEX(I) RHOSTAR(L) = PSTAR(L)/(R*T1(L)) enddo !---------------------------------------------------------------------- ! Calculate the saturated specific humidity, its gradient w.r.t. ! temperature, and the humidity deficit. !---------------------------------------------------------------------- call QSAT (POINTS, LAND_PTS, LAND_INDEX, EPSILON, ZERODEGC &, QS1, T1, PSTAR) !CDIR NODEP do I=1,LAND_PTS L = LAND_INDEX(I) if (LYING_SNOW(L) .le. 50.) then LAT(L) = LC else LAT(L) = LC + LF endif DQS_DT(L) = (EPSILON*LAT(L)*QS1(L))/(R*T1(L)*T1(L)) DQ1(L) = QS1(L) - Q1(L) !----------------------------------------------------------------------- ! Calculate available energy when the surface temperature is equal to ! the atmospheric temperature (AHAT). !----------------------------------------------------------------------- AHAT(L) = SWN(L) + LW(L) - SIGMA*T1(L)*T1(L)*T1(L)*T1(L) & - AS1*(T1(L)-TS1(L)) !----------------------------------------------------------------------- ! 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)) !---------------------------------------------------------------------- ! Calculate the evaporation rate and diagnose the surface temperature. !---------------------------------------------------------------------- RESF = 1.0/(1.0+RS(L)/RA(L)) DUM = RHOSTAR(L)*CP/RA(L) + 4*SIGMA*T1(L)**3 + AS1 NUMER = (DQS_DT(L)*AHAT(L) + DUM*DQ1(L)) * RESF DENOM = RESF*LAT(L)*DQS_DT(L) + RA(L)*DUM/RHOSTAR(L) E(L) = NUMER / DENOM if (MNEG(L).lt.0.) E(L) = amin1(E(L),0.) LE(L) = LAT(L) * E(L) TSTAR(L) = T1(L) + & (AHAT(L) - LAT(L)*RHOSTAR(L)*DQ1(L)*RESF/RA(L)) & /(DUM + DQS_DT(L)*LAT(L)*RHOSTAR(L)*RESF/RA(L)) 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) G(L) = RADNET(L) - LE(L) - SH(L) enddo # endif return end #endif