module CanopyFluxesMod !------------------------------------------------------------------------------ !BOP ! ! !MODULE: CanopyFluxesMod ! ! !DESCRIPTION: ! Calculates the leaf temperature and the leaf fluxes, ! transpiration, photosynthesis and updates the dew ! accumulation due to evaporation. ! ! !USES: use abortutils, only: endrun use clm_varctl, only: iulog use shr_sys_mod, only: shr_sys_flush ! ! !PUBLIC TYPES: implicit none save ! ! !PUBLIC MEMBER FUNCTIONS: public :: CanopyFluxes !Calculates the leaf temperature and leaf fluxes ! ! !PRIVATE MEMBER FUNCTIONS: private :: Stomata !Leaf stomatal resistance and leaf photosynthesis ! ! !REVISION HISTORY: ! Created by Mariana Vertenstein ! 4/25/05, Peter Thornton: replaced old Stomata subroutine with what ! used to be called StomataCN, as part of migration to new sun/shade ! algorithms. ! !EOP !------------------------------------------------------------------------------ contains !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: CanopyFluxes ! ! !INTERFACE: subroutine CanopyFluxes(lbg, ubg, lbc, ubc, lbp, ubp, & num_nolakep, filter_nolakep) ! ! !DESCRIPTION: ! 1. Calculates the leaf temperature: ! 2. Calculates the leaf fluxes, transpiration, photosynthesis and ! updates the dew accumulation due to evaporation. ! ! Method: ! Use the Newton-Raphson iteration to solve for the foliage ! temperature that balances the surface energy budget: ! ! f(t_veg) = Net radiation - Sensible - Latent = 0 ! f(t_veg) + d(f)/d(t_veg) * dt_veg = 0 (*) ! ! Note: ! (1) In solving for t_veg, t_grnd is given from the previous timestep. ! (2) The partial derivatives of aerodynamical resistances, which cannot ! be determined analytically, are ignored for d(H)/dT and d(LE)/dT ! (3) The weighted stomatal resistance of sunlit and shaded foliage is used ! (4) Canopy air temperature and humidity are derived from => Hc + Hg = Ha ! => Ec + Eg = Ea ! (5) Energy loss is due to: numerical truncation of energy budget equation ! (*); and "ecidif" (see the code) which is dropped into the sensible ! heat ! (6) The convergence criteria: the difference, del = t_veg(n+1)-t_veg(n) ! and del2 = t_veg(n)-t_veg(n-1) less than 0.01 K, and the difference ! of water flux from the leaf between the iteration step (n+1) and (n) ! less than 0.1 W/m2; or the iterative steps over 40. ! ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8 use clmtype use clm_atmlnd , only : clm_a2l use clm_time_manager , only : get_step_size, get_prev_date use clm_varpar , only : nlevgrnd, nlevsno use clm_varcon , only : sb, cpair, hvap, vkc, grav, denice, & denh2o, tfrz, csoilc, tlsai_crit, alpha_aero, & isecspday, degpsec use shr_const_mod , only : SHR_CONST_TKFRZ use pftvarcon , only : nirrig use QSatMod , only : QSat use FrictionVelocityMod, only : FrictionVelocity, MoninObukIni use spmdMod , only : masterproc ! ! !ARGUMENTS: implicit none integer, intent(in) :: lbg, ubg ! gridcell bounds integer, intent(in) :: lbc, ubc ! column bounds integer, intent(in) :: lbp, ubp ! pft bounds integer, intent(in) :: num_nolakep ! number of column non-lake points in pft filter integer, intent(in) :: filter_nolakep(ubp-lbp+1) ! pft filter for non-lake points ! ! !CALLED FROM: ! subroutine Biogeophysics1 in module Biogeophysics1Mod ! ! !REVISION HISTORY: ! 15 September 1999: Yongjiu Dai; Initial code ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision ! 12/19/01, Peter Thornton ! Changed tg to t_grnd for consistency with other routines ! 1/29/02, Peter Thornton ! Migrate to new data structures, new calling protocol. For now co2 and ! o2 partial pressures are hardwired, but they should be coming in from ! forc_pco2 and forc_po2. Keeping the same hardwired values as in CLM2 to ! assure bit-for-bit results in the first comparisons. ! 27 February 2008: Keith Oleson; Sparse/dense aerodynamic parameters from ! X. Zeng ! 6 March 2009: Peter Thornton; Daylength control on Vcmax, from Bill Bauerle ! ! !LOCAL VARIABLES: ! ! local pointers to implicit in variables ! integer , pointer :: frac_veg_nosno(:) ! frac of veg not covered by snow (0 OR 1 now) [-] integer , pointer :: ivt(:) ! pft vegetation type integer , pointer :: pcolumn(:) ! pft's column index integer , pointer :: plandunit(:) ! pft's landunit index integer , pointer :: pgridcell(:) ! pft's gridcell index real(r8), pointer :: forc_th(:) ! atmospheric potential temperature (Kelvin) real(r8), pointer :: t_grnd(:) ! ground surface temperature [K] real(r8), pointer :: thm(:) ! intermediate variable (forc_t+0.0098*forc_hgt_t_pft) real(r8), pointer :: qg(:) ! specific humidity at ground surface [kg/kg] real(r8), pointer :: thv(:) ! virtual potential temperature (kelvin) real(r8), pointer :: z0mv(:) ! roughness length over vegetation, momentum [m] real(r8), pointer :: z0hv(:) ! roughness length over vegetation, sensible heat [m] real(r8), pointer :: z0qv(:) ! roughness length over vegetation, latent heat [m] real(r8), pointer :: z0mg(:) ! roughness length of ground, momentum [m] real(r8), pointer :: dqgdT(:) ! temperature derivative of "qg" real(r8), pointer :: htvp(:) ! latent heat of evaporation (/sublimation) [J/kg] real(r8), pointer :: emv(:) ! ground emissivity real(r8), pointer :: emg(:) ! vegetation emissivity real(r8), pointer :: forc_pbot(:) ! atmospheric pressure (Pa) real(r8), pointer :: forc_pco2(:) ! partial pressure co2 (Pa) #if (defined C13) ! 4/14/05: PET ! Adding isotope code real(r8), pointer :: forc_pc13o2(:) ! partial pressure c13o2 (Pa) #endif real(r8), pointer :: forc_po2(:) ! partial pressure o2 (Pa) real(r8), pointer :: forc_q(:) ! atmospheric specific humidity (kg/kg) real(r8), pointer :: forc_u(:) ! atmospheric wind speed in east direction (m/s) real(r8), pointer :: forc_v(:) ! atmospheric wind speed in north direction (m/s) real(r8), pointer :: forc_hgt_u_pft(:) !observational height of wind at pft level [m] real(r8), pointer :: forc_rho(:) ! density (kg/m**3) real(r8), pointer :: forc_lwrad(:) ! downward infrared (longwave) radiation (W/m**2) real(r8), pointer :: displa(:) ! displacement height (m) real(r8), pointer :: elai(:) ! one-sided leaf area index with burying by snow real(r8), pointer :: esai(:) ! one-sided stem area index with burying by snow real(r8), pointer :: fdry(:) ! fraction of foliage that is green and dry [-] real(r8), pointer :: fwet(:) ! fraction of canopy that is wet (0 to 1) real(r8), pointer :: laisun(:) ! sunlit leaf area real(r8), pointer :: laisha(:) ! shaded leaf area real(r8), pointer :: sabv(:) ! solar radiation absorbed by vegetation (W/m**2) real(r8), pointer :: watsat(:,:) ! volumetric soil water at saturation (porosity) real(r8), pointer :: watdry(:,:) ! btran parameter for btran=0 real(r8), pointer :: watopt(:,:) ! btran parameter for btran = 1 real(r8), pointer :: h2osoi_ice(:,:)! ice lens (kg/m2) real(r8), pointer :: h2osoi_liq(:,:)! liquid water (kg/m2) real(r8), pointer :: dz(:,:) ! layer depth (m) real(r8), pointer :: t_soisno(:,:) ! soil temperature (Kelvin) real(r8), pointer :: sucsat(:,:) ! minimum soil suction (mm) real(r8), pointer :: bsw(:,:) ! Clapp and Hornberger "b" real(r8), pointer :: rootfr(:,:) ! fraction of roots in each soil layer real(r8), pointer :: dleaf(:) ! characteristic leaf dimension (m) real(r8), pointer :: smpso(:) ! soil water potential at full stomatal opening (mm) real(r8), pointer :: smpsc(:) ! soil water potential at full stomatal closure (mm) real(r8), pointer :: frac_sno(:) ! fraction of ground covered by snow (0 to 1) real(r8), pointer :: htop(:) ! canopy top(m) real(r8), pointer :: snowdp(:) ! snow height (m) real(r8), pointer :: soilbeta(:) ! soil wetness relative to field capacity real(r8), pointer :: lat(:) ! latitude (radians) real(r8), pointer :: decl(:) ! declination angle (radians) real(r8), pointer :: max_dayl(:) !maximum daylength for this column (s) real(r8), pointer :: londeg(:) ! longitude (degrees) (for calculation of local time) ! ! local pointers to implicit inout arguments ! real(r8), pointer :: cgrnds(:) ! deriv. of soil sensible heat flux wrt soil temp [w/m2/k] real(r8), pointer :: cgrndl(:) ! deriv. of soil latent heat flux wrt soil temp [w/m**2/k] real(r8), pointer :: t_veg(:) ! vegetation temperature (Kelvin) real(r8), pointer :: t_ref2m(:) ! 2 m height surface air temperature (Kelvin) real(r8), pointer :: q_ref2m(:) ! 2 m height surface specific humidity (kg/kg) real(r8), pointer :: t_ref2m_r(:) ! Rural 2 m height surface air temperature (Kelvin) real(r8), pointer :: rh_ref2m(:) ! 2 m height surface relative humidity (%) real(r8), pointer :: rh_ref2m_r(:) ! Rural 2 m height surface relative humidity (%) real(r8), pointer :: h2ocan(:) ! canopy water (mm H2O) real(r8), pointer :: cisun(:) !sunlit intracellular CO2 (Pa) real(r8), pointer :: cisha(:) !shaded intracellular CO2 (Pa) ! ! local pointers to implicit out arguments ! real(r8), pointer :: rb1(:) ! boundary layer resistance (s/m) real(r8), pointer :: cgrnd(:) ! deriv. of soil energy flux wrt to soil temp [w/m2/k] real(r8), pointer :: dlrad(:) ! downward longwave radiation below the canopy [W/m2] real(r8), pointer :: ulrad(:) ! upward longwave radiation above the canopy [W/m2] real(r8), pointer :: ram1(:) ! aerodynamical resistance (s/m) real(r8), pointer :: btran(:) ! transpiration wetness factor (0 to 1) real(r8), pointer :: rssun(:) ! sunlit stomatal resistance (s/m) real(r8), pointer :: rssha(:) ! shaded stomatal resistance (s/m) real(r8), pointer :: psnsun(:) ! sunlit leaf photosynthesis (umol CO2 /m**2/ s) real(r8), pointer :: psnsha(:) ! shaded leaf photosynthesis (umol CO2 /m**2/ s) #if (defined C13) ! 4/14/05: PET ! Adding isotope code real(r8), pointer :: c13_psnsun(:) ! sunlit leaf photosynthesis (umol 13CO2 /m**2/ s) real(r8), pointer :: c13_psnsha(:) ! shaded leaf photosynthesis (umol 13CO2 /m**2/ s) ! 4/21/05: PET ! Adding isotope code real(r8), pointer :: rc13_canair(:) !C13O2/C12O2 in canopy air real(r8), pointer :: rc13_psnsun(:) !C13O2/C12O2 in sunlit canopy psn flux real(r8), pointer :: rc13_psnsha(:) !C13O2/C12O2 in shaded canopy psn flux real(r8), pointer :: alphapsnsun(:) !fractionation factor in sunlit canopy psn flux real(r8), pointer :: alphapsnsha(:) !fractionation factor in shaded canopy psn flux #endif real(r8), pointer :: qflx_tran_veg(:) ! vegetation transpiration (mm H2O/s) (+ = to atm) real(r8), pointer :: dt_veg(:) ! change in t_veg, last iteration (Kelvin) real(r8), pointer :: qflx_evap_veg(:) ! vegetation evaporation (mm H2O/s) (+ = to atm) real(r8), pointer :: eflx_sh_veg(:) ! sensible heat flux from leaves (W/m**2) [+ to atm] real(r8), pointer :: taux(:) ! wind (shear) stress: e-w (kg/m/s**2) real(r8), pointer :: tauy(:) ! wind (shear) stress: n-s (kg/m/s**2) real(r8), pointer :: eflx_sh_grnd(:) ! sensible heat flux from ground (W/m**2) [+ to atm] real(r8), pointer :: qflx_evap_soi(:) ! soil evaporation (mm H2O/s) (+ = to atm) real(r8), pointer :: fpsn(:) ! photosynthesis (umol CO2 /m**2 /s) real(r8), pointer :: rootr(:,:) ! effective fraction of roots in each soil layer real(r8), pointer :: rresis(:,:) ! root resistance by layer (0-1) (nlevgrnd) real(r8), pointer :: irrig_rate(:) ! current irrigation rate [mm/s] integer, pointer :: n_irrig_steps_left(:) ! # of time steps for which we still need to irrigate today ! ! ! !OTHER LOCAL VARIABLES: !EOP ! real(r8), parameter :: btran0 = 0.0_r8 ! initial value real(r8), parameter :: zii = 1000.0_r8 ! convective boundary layer height [m] real(r8), parameter :: beta = 1.0_r8 ! coefficient of conective velocity [-] real(r8), parameter :: delmax = 1.0_r8 ! maxchange in leaf temperature [K] real(r8), parameter :: dlemin = 0.1_r8 ! max limit for energy flux convergence [w/m2] real(r8), parameter :: dtmin = 0.01_r8 ! max limit for temperature convergence [K] integer , parameter :: itmax = 40 ! maximum number of iteration [-] integer , parameter :: itmin = 2 ! minimum number of iteration [-] real(r8), parameter :: irrig_min_lai = 0.0_r8 ! Minimum LAI for irrigation real(r8), parameter :: irrig_btran_thresh = 0.999999_r8 ! Irrigate when btran falls below 0.999999 rather than 1 to allow for round-off error integer , parameter :: irrig_start_time = isecspday/4 ! Time of day to check whether we need irrigation, seconds (0 = midnight). ! We start applying the irrigation in the time step FOLLOWING this time, ! since we won't begin irrigating until the next call to Hydrology1 integer , parameter :: irrig_length = isecspday/6 ! Desired amount of time to irrigate per day (sec). Actual time may ! differ if this is not a multiple of dtime. Irrigation won't work properly ! if dtime > secsperday real(r8), parameter :: irrig_factor = 0.7_r8 ! Determines target soil moisture level for irrigation. If h2osoi_liq_so ! is the soil moisture level at which stomata are fully open and ! h2osoi_liq_sat is the soil moisture level at saturation (eff_porosity), ! then the target soil moisture level is ! (h2osoi_liq_so + irrig_factor*(h2osoi_liq_sat - h2osoi_liq_so)). ! A value of 0 means that the target soil moisture level is h2osoi_liq_so. ! A value of 1 means that the target soil moisture level is h2osoi_liq_sat !added by K.Sakaguchi for litter resistance real(r8), parameter :: lai_dl = 0.5_r8 ! placeholder for (dry) plant litter area index (m2/m2) real(r8), parameter :: z_dl = 0.05_r8 ! placeholder for (dry) litter layer thickness (m) !added by K.Sakaguchi for stability formulation real(r8), parameter :: ria = 0.5_r8 ! free parameter for stable formulation (currently = 0.5, "gamma" in Sakaguchi&Zeng,2008) real(r8) :: dtime ! land model time step (sec) real(r8) :: zldis(lbp:ubp) ! reference height "minus" zero displacement height [m] real(r8) :: zeta ! dimensionless height used in Monin-Obukhov theory real(r8) :: wc ! convective velocity [m/s] real(r8) :: dth(lbp:ubp) ! diff of virtual temp. between ref. height and surface real(r8) :: dthv(lbp:ubp) ! diff of vir. poten. temp. between ref. height and surface real(r8) :: dqh(lbp:ubp) ! diff of humidity between ref. height and surface real(r8) :: obu(lbp:ubp) ! Monin-Obukhov length (m) real(r8) :: um(lbp:ubp) ! wind speed including the stablity effect [m/s] real(r8) :: ur(lbp:ubp) ! wind speed at reference height [m/s] real(r8) :: uaf(lbp:ubp) ! velocity of air within foliage [m/s] real(r8) :: temp1(lbp:ubp) ! relation for potential temperature profile real(r8) :: temp12m(lbp:ubp) ! relation for potential temperature profile applied at 2-m real(r8) :: temp2(lbp:ubp) ! relation for specific humidity profile real(r8) :: temp22m(lbp:ubp) ! relation for specific humidity profile applied at 2-m real(r8) :: ustar(lbp:ubp) ! friction velocity [m/s] real(r8) :: tstar ! temperature scaling parameter real(r8) :: qstar ! moisture scaling parameter real(r8) :: thvstar ! virtual potential temperature scaling parameter real(r8) :: taf(lbp:ubp) ! air temperature within canopy space [K] real(r8) :: qaf(lbp:ubp) ! humidity of canopy air [kg/kg] real(r8) :: rpp ! fraction of potential evaporation from leaf [-] real(r8) :: rppdry ! fraction of potential evaporation through transp [-] real(r8) :: cf ! heat transfer coefficient from leaves [-] real(r8) :: rb(lbp:ubp) ! leaf boundary layer resistance [s/m] real(r8) :: rah(lbp:ubp,2) ! thermal resistance [s/m] real(r8) :: raw(lbp:ubp,2) ! moisture resistance [s/m] real(r8) :: wta ! heat conductance for air [m/s] real(r8) :: wtg(lbp:ubp) ! heat conductance for ground [m/s] real(r8) :: wtl ! heat conductance for leaf [m/s] real(r8) :: wta0(lbp:ubp) ! normalized heat conductance for air [-] real(r8) :: wtl0(lbp:ubp) ! normalized heat conductance for leaf [-] real(r8) :: wtg0 ! normalized heat conductance for ground [-] real(r8) :: wtal(lbp:ubp) ! normalized heat conductance for air and leaf [-] real(r8) :: wtga ! normalized heat cond. for air and ground [-] real(r8) :: wtaq ! latent heat conductance for air [m/s] real(r8) :: wtlq ! latent heat conductance for leaf [m/s] real(r8) :: wtgq(lbp:ubp) ! latent heat conductance for ground [m/s] real(r8) :: wtaq0(lbp:ubp) ! normalized latent heat conductance for air [-] real(r8) :: wtlq0(lbp:ubp) ! normalized latent heat conductance for leaf [-] real(r8) :: wtgq0 ! normalized heat conductance for ground [-] real(r8) :: wtalq(lbp:ubp) ! normalized latent heat cond. for air and leaf [-] real(r8) :: wtgaq ! normalized latent heat cond. for air and ground [-] real(r8) :: el(lbp:ubp) ! vapor pressure on leaf surface [pa] real(r8) :: deldT ! derivative of "el" on "t_veg" [pa/K] real(r8) :: qsatl(lbp:ubp) ! leaf specific humidity [kg/kg] real(r8) :: qsatldT(lbp:ubp) ! derivative of "qsatl" on "t_veg" real(r8) :: e_ref2m ! 2 m height surface saturated vapor pressure [Pa] real(r8) :: de2mdT ! derivative of 2 m height surface saturated vapor pressure on t_ref2m real(r8) :: qsat_ref2m ! 2 m height surface saturated specific humidity [kg/kg] real(r8) :: dqsat2mdT ! derivative of 2 m height surface saturated specific humidity on t_ref2m real(r8) :: air(lbp:ubp),bir(lbp:ubp),cir(lbp:ubp) ! atmos. radiation temporay set real(r8) :: dc1,dc2 ! derivative of energy flux [W/m2/K] real(r8) :: delt ! temporary real(r8) :: delq(lbp:ubp) ! temporary real(r8) :: del(lbp:ubp) ! absolute change in leaf temp in current iteration [K] real(r8) :: del2(lbp:ubp) ! change in leaf temperature in previous iteration [K] real(r8) :: dele(lbp:ubp) ! change in latent heat flux from leaf [K] real(r8) :: dels ! change in leaf temperature in current iteration [K] real(r8) :: det(lbp:ubp) ! maximum leaf temp. change in two consecutive iter [K] real(r8) :: efeb(lbp:ubp) ! latent heat flux from leaf (previous iter) [mm/s] real(r8) :: efeold ! latent heat flux from leaf (previous iter) [mm/s] real(r8) :: efpot ! potential latent energy flux [kg/m2/s] real(r8) :: efe(lbp:ubp) ! water flux from leaf [mm/s] real(r8) :: efsh ! sensible heat from leaf [mm/s] real(r8) :: obuold(lbp:ubp) ! monin-obukhov length from previous iteration real(r8) :: tlbef(lbp:ubp) ! leaf temperature from previous iteration [K] real(r8) :: ecidif ! excess energies [W/m2] real(r8) :: err(lbp:ubp) ! balance error real(r8) :: erre ! balance error real(r8) :: co2(lbp:ubp) ! atmospheric co2 partial pressure (pa) #if (defined C13) ! 4/14/05: PET ! Adding isotope code real(r8) :: c13o2(lbp:ubp) ! atmospheric c13o2 partial pressure (pa) #endif real(r8) :: o2(lbp:ubp) ! atmospheric o2 partial pressure (pa) real(r8) :: svpts(lbp:ubp) ! saturation vapor pressure at t_veg (pa) real(r8) :: eah(lbp:ubp) ! canopy air vapor pressure (pa) real(r8) :: s_node ! vol_liq/eff_porosity real(r8) :: smp_node ! matrix potential real(r8) :: vol_ice ! partial volume of ice lens in layer real(r8) :: eff_porosity ! effective porosity in layer real(r8) :: vol_liq ! partial volume of liquid water in layer integer :: itlef ! counter for leaf temperature iteration [-] integer :: nmozsgn(lbp:ubp) ! number of times stability changes sign real(r8) :: w ! exp(-LSAI) real(r8) :: csoilcn ! interpolated csoilc for less than dense canopies real(r8) :: fm(lbp:ubp) ! needed for BGC only to diagnose 10m wind speed real(r8) :: wtshi ! sensible heat resistance for air, grnd and leaf [-] real(r8) :: wtsqi ! latent heat resistance for air, grnd and leaf [-] integer :: j ! soil/snow level index integer :: p ! pft index integer :: c ! column index integer :: l ! landunit index integer :: g ! gridcell index integer :: fp ! lake filter pft index integer :: fn ! number of values in pft filter integer :: fnorig ! number of values in pft filter copy integer :: fnold ! temporary copy of pft count integer :: f ! filter index integer :: filterp(ubp-lbp+1) ! temporary filter integer :: fporig(ubp-lbp+1) ! temporary filter real(r8) :: displa_loc(lbp:ubp) ! temporary copy real(r8) :: z0mv_loc(lbp:ubp) ! temporary copy real(r8) :: z0hv_loc(lbp:ubp) ! temporary copy real(r8) :: z0qv_loc(lbp:ubp) ! temporary copy logical :: found ! error flag for canopy above forcing hgt integer :: index ! pft index for error real(r8) :: egvf ! effective green vegetation fraction real(r8) :: lt ! elai+esai real(r8) :: ri ! stability parameter for under canopy air (unitless) real(r8) :: csoilb ! turbulent transfer coefficient over bare soil (unitless) real(r8) :: ricsoilc ! modified transfer coefficient under dense canopy (unitless) real(r8) :: snowdp_c ! critical snow depth to cover plant litter (m) real(r8) :: rdl ! dry litter layer resistance for water vapor (s/m) real(r8) :: elai_dl ! exposed (dry) plant litter area index real(r8) :: fsno_dl ! effective snow cover over plant litter real(r8) :: dayl ! daylength (s) real(r8) :: temp ! temporary, for daylength calculation real(r8) :: dayl_factor(lbp:ubp) ! scalar (0-1) for daylength effect on Vcmax integer :: yr ! year at start of time step integer :: mon ! month at start of time step integer :: day ! day at start of time step integer :: time ! time at start of time step (seconds after 0Z) integer :: local_time ! local time at start of time step (seconds after solar midnight) integer :: irrig_nsteps_per_day ! number of time steps per day in which we irrigate logical :: check_for_irrig(lbp:ubp) ! where do we need to check soil moisture to see if we need to irrigate? logical :: frozen_soil(lbc:ubc) ! set to true if we have encountered a frozen soil layer real(r8) :: vol_liq_so ! partial volume of liquid water in layer for which smp_node = smpso real(r8) :: h2osoi_liq_so ! liquid water corresponding to vol_liq_so for this layer [kg/m2] real(r8) :: h2osoi_liq_sat ! liquid water corresponding to eff_porosity for this layer [kg/m2] real(r8) :: deficit ! difference between desired soil moisture level for this layer and current soil moisture level [kg/m2] !------------------------------------------------------------------------------ ! Assign local pointers to derived type members (gridcell-level) forc_lwrad => clm_a2l%forc_lwrad forc_pco2 => clm_a2l%forc_pco2 #if (defined C13) forc_pc13o2 => clm_a2l%forc_pc13o2 #endif forc_po2 => clm_a2l%forc_po2 forc_q => clm_a2l%forc_q forc_pbot => clm_a2l%forc_pbot forc_u => clm_a2l%forc_u forc_v => clm_a2l%forc_v forc_th => clm_a2l%forc_th forc_rho => clm_a2l%forc_rho lat => clm3%g%lat londeg => clm3%g%londeg ! Assign local pointers to derived type members (column-level) t_soisno => clm3%g%l%c%ces%t_soisno watsat => clm3%g%l%c%cps%watsat watdry => clm3%g%l%c%cps%watdry watopt => clm3%g%l%c%cps%watopt h2osoi_ice => clm3%g%l%c%cws%h2osoi_ice dz => clm3%g%l%c%cps%dz h2osoi_liq => clm3%g%l%c%cws%h2osoi_liq sucsat => clm3%g%l%c%cps%sucsat bsw => clm3%g%l%c%cps%bsw emg => clm3%g%l%c%cps%emg t_grnd => clm3%g%l%c%ces%t_grnd qg => clm3%g%l%c%cws%qg thv => clm3%g%l%c%ces%thv dqgdT => clm3%g%l%c%cws%dqgdT htvp => clm3%g%l%c%cps%htvp z0mg => clm3%g%l%c%cps%z0mg frac_sno => clm3%g%l%c%cps%frac_sno snowdp => clm3%g%l%c%cps%snowdp soilbeta => clm3%g%l%c%cws%soilbeta decl => clm3%g%l%c%cps%decl max_dayl => clm3%g%l%c%cps%max_dayl ! Assign local pointers to derived type members (pft-level) rb1 => clm3%g%l%c%p%pps%rb1 ivt => clm3%g%l%c%p%itype pcolumn => clm3%g%l%c%p%column plandunit => clm3%g%l%c%p%landunit pgridcell => clm3%g%l%c%p%gridcell frac_veg_nosno => clm3%g%l%c%p%pps%frac_veg_nosno btran => clm3%g%l%c%p%pps%btran rootfr => clm3%g%l%c%p%pps%rootfr rootr => clm3%g%l%c%p%pps%rootr rresis => clm3%g%l%c%p%pps%rresis emv => clm3%g%l%c%p%pps%emv t_veg => clm3%g%l%c%p%pes%t_veg displa => clm3%g%l%c%p%pps%displa z0mv => clm3%g%l%c%p%pps%z0mv z0hv => clm3%g%l%c%p%pps%z0hv z0qv => clm3%g%l%c%p%pps%z0qv ram1 => clm3%g%l%c%p%pps%ram1 htop => clm3%g%l%c%p%pps%htop rssun => clm3%g%l%c%p%pps%rssun rssha => clm3%g%l%c%p%pps%rssha cisun => clm3%g%l%c%p%pps%cisun cisha => clm3%g%l%c%p%pps%cisha psnsun => clm3%g%l%c%p%pcf%psnsun psnsha => clm3%g%l%c%p%pcf%psnsha #if (defined C13) ! 4/14/05: PET ! Adding isotope code c13_psnsun => clm3%g%l%c%p%pc13f%psnsun c13_psnsha => clm3%g%l%c%p%pc13f%psnsha ! 4/21/05: PET ! Adding isotope code rc13_canair => clm3%g%l%c%p%pepv%rc13_canair rc13_psnsun => clm3%g%l%c%p%pepv%rc13_psnsun rc13_psnsha => clm3%g%l%c%p%pepv%rc13_psnsha alphapsnsun => clm3%g%l%c%p%pps%alphapsnsun alphapsnsha => clm3%g%l%c%p%pps%alphapsnsha #endif elai => clm3%g%l%c%p%pps%elai esai => clm3%g%l%c%p%pps%esai fdry => clm3%g%l%c%p%pps%fdry laisun => clm3%g%l%c%p%pps%laisun laisha => clm3%g%l%c%p%pps%laisha qflx_tran_veg => clm3%g%l%c%p%pwf%qflx_tran_veg fwet => clm3%g%l%c%p%pps%fwet h2ocan => clm3%g%l%c%p%pws%h2ocan dt_veg => clm3%g%l%c%p%pps%dt_veg sabv => clm3%g%l%c%p%pef%sabv qflx_evap_veg => clm3%g%l%c%p%pwf%qflx_evap_veg eflx_sh_veg => clm3%g%l%c%p%pef%eflx_sh_veg taux => clm3%g%l%c%p%pmf%taux tauy => clm3%g%l%c%p%pmf%tauy eflx_sh_grnd => clm3%g%l%c%p%pef%eflx_sh_grnd qflx_evap_soi => clm3%g%l%c%p%pwf%qflx_evap_soi t_ref2m => clm3%g%l%c%p%pes%t_ref2m q_ref2m => clm3%g%l%c%p%pes%q_ref2m t_ref2m_r => clm3%g%l%c%p%pes%t_ref2m_r rh_ref2m_r => clm3%g%l%c%p%pes%rh_ref2m_r rh_ref2m => clm3%g%l%c%p%pes%rh_ref2m dlrad => clm3%g%l%c%p%pef%dlrad ulrad => clm3%g%l%c%p%pef%ulrad cgrnds => clm3%g%l%c%p%pef%cgrnds cgrndl => clm3%g%l%c%p%pef%cgrndl cgrnd => clm3%g%l%c%p%pef%cgrnd fpsn => clm3%g%l%c%p%pcf%fpsn forc_hgt_u_pft => clm3%g%l%c%p%pps%forc_hgt_u_pft thm => clm3%g%l%c%p%pes%thm irrig_rate => clm3%g%l%c%cps%irrig_rate n_irrig_steps_left => clm3%g%l%c%cps%n_irrig_steps_left ! Assign local pointers to derived type members (ecophysiological) dleaf => pftcon%dleaf smpso => pftcon%smpso smpsc => pftcon%smpsc ! Determine step size dtime = get_step_size() irrig_nsteps_per_day = ((irrig_length + (dtime - 1))/dtime) ! round up ! Filter pfts where frac_veg_nosno is non-zero fn = 0 do fp = 1,num_nolakep p = filter_nolakep(fp) if (frac_veg_nosno(p) /= 0) then fn = fn + 1 filterp(fn) = p end if end do ! Initialize do f = 1, fn p = filterp(f) del(p) = 0._r8 ! change in leaf temperature from previous iteration efeb(p) = 0._r8 ! latent head flux from leaf for previous iteration wtlq0(p) = 0._r8 wtalq(p) = 0._r8 wtgq(p) = 0._r8 wtaq0(p) = 0._r8 obuold(p) = 0._r8 btran(p) = btran0 end do ! calculate daylength control for Vcmax do f = 1, fn p=filterp(f) c=pcolumn(p) g=pgridcell(p) ! calculate daylength temp = -(sin(lat(g))*sin(decl(c)))/(cos(lat(g)) * cos(decl(c))) temp = min(1._r8,max(-1._r8,temp)) dayl = 2.0_r8 * 13750.9871_r8 * acos(temp) ! calculate dayl_factor as the ratio of (current:max dayl)^2 ! set a minimum of 0.01 (1%) for the dayl_factor dayl_factor(p)=min(1._r8,max(0.01_r8,(dayl*dayl)/(max_dayl(c)*max_dayl(c)))) end do rb1(lbp:ubp) = 0._r8 ! Effective porosity of soil, partial volume of ice and liquid (needed for btran) ! and root resistance factors do j = 1,nlevgrnd do f = 1, fn p = filterp(f) c = pcolumn(p) l = plandunit(p) ! Root resistance factors vol_ice = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice)) eff_porosity = watsat(c,j)-vol_ice vol_liq = min(eff_porosity, h2osoi_liq(c,j)/(dz(c,j)*denh2o)) if (vol_liq .le. 0._r8 .or. t_soisno(c,j) .le. tfrz-2._r8) then rootr(p,j) = 0._r8 else s_node = max(vol_liq/eff_porosity,0.01_r8) smp_node = max(smpsc(ivt(p)), -sucsat(c,j)*s_node**(-bsw(c,j))) rresis(p,j) = min( (eff_porosity/watsat(c,j))* & (smp_node - smpsc(ivt(p))) / (smpso(ivt(p)) - smpsc(ivt(p))), 1._r8) rootr(p,j) = rootfr(p,j)*rresis(p,j) btran(p) = btran(p) + rootr(p,j) endif end do end do ! Normalize root resistances to get layer contribution to ET do j = 1,nlevgrnd do f = 1, fn p = filterp(f) if (btran(p) .gt. btran0) then rootr(p,j) = rootr(p,j)/btran(p) else rootr(p,j) = 0._r8 end if end do end do ! Determine if irrigation is needed (over irrigated soil columns) ! First, determine in what grid cells we need to bother 'measuring' soil water, to see if we need irrigation ! Also set n_irrig_steps_left for these grid cells ! n_irrig_steps_left(p) > 0 is ok even if irrig_rate(p) ends up = 0 ! in this case, we'll irrigate by 0 for the given number of time steps call get_prev_date(yr, mon, day, time) ! get time as of beginning of time step do f = 1, fn p = filterp(f) c = pcolumn(p) g = pgridcell(p) if (ivt(p) == nirrig .and. elai(p) > irrig_min_lai .and. btran(p) < irrig_btran_thresh) then ! see if it's the right time of day to start irrigating: local_time = modulo(time + nint(londeg(g)/degpsec), isecspday) if (modulo(local_time - irrig_start_time, isecspday) < dtime) then ! it's time to start irrigating check_for_irrig(p) = .true. n_irrig_steps_left(c) = irrig_nsteps_per_day irrig_rate(c) = 0._r8 ! reset; we'll add to this later else check_for_irrig(p) = .false. end if else ! non-irrig pft or elai<=irrig_min_lai or btran>irrig_btran_thresh check_for_irrig(p) = .false. end if end do ! Now 'measure' soil water for the grid cells identified above and see if the soil is dry enough to warrant irrigation frozen_soil(:) = .false. do j = 1,nlevgrnd do f = 1, fn p = filterp(f) c = pcolumn(p) if (check_for_irrig(p) .and. .not. frozen_soil(c)) then ! if level L was frozen, then we don't look at any levels below L if (t_soisno(c,j) <= SHR_CONST_TKFRZ) then frozen_soil(c) = .true. else if (rootfr(p,j) > 0._r8) then ! determine soil water deficit in this layer: ! Calculate vol_liq_so - i.e., vol_liq at which smp_node = smpso - by inverting the above equations ! for the root resistance factors vol_ice = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice)) ! this duplicates the above equation for vol_ice eff_porosity = watsat(c,j)-vol_ice ! this duplicates the above equation for eff_porosity vol_liq_so = eff_porosity * (-smpso(ivt(p))/sucsat(c,j))**(-1/bsw(c,j)) ! Translate vol_liq_so and eff_porosity into h2osoi_liq_so and h2osoi_liq_sat and calculate deficit h2osoi_liq_so = vol_liq_so * denh2o * dz(c,j) h2osoi_liq_sat = eff_porosity * denh2o * dz(c,j) deficit = max((h2osoi_liq_so + irrig_factor*(h2osoi_liq_sat - h2osoi_liq_so)) - h2osoi_liq(c,j), 0._r8) ! Add deficit to irrig_rate, converting units from mm to mm/sec irrig_rate(c) = irrig_rate(c) + deficit/(dtime*irrig_nsteps_per_day) end if ! else if (rootfr(p,j) .gt. 0) end if ! if (check_for_irrig(p) .and. .not. frozen_soil(c)) end do ! do f end do ! do j ! Modify aerodynamic parameters for sparse/dense canopy (X. Zeng) do f = 1, fn p = filterp(f) c = pcolumn(p) lt = min(elai(p)+esai(p), tlsai_crit) egvf =(1._r8 - alpha_aero * exp(-lt)) / (1._r8 - alpha_aero * exp(-tlsai_crit)) displa(p) = egvf * displa(p) z0mv(p) = exp(egvf * log(z0mv(p)) + (1._r8 - egvf) * log(z0mg(c))) z0hv(p) = z0mv(p) z0qv(p) = z0mv(p) end do found = .false. do f = 1, fn p = filterp(f) c = pcolumn(p) g = pgridcell(p) ! Net absorbed longwave radiation by canopy and ground ! =air+bir*t_veg**4+cir*t_grnd(c)**4 air(p) = emv(p) * (1._r8+(1._r8-emv(p))*(1._r8-emg(c))) * forc_lwrad(g) bir(p) = - (2._r8-emv(p)*(1._r8-emg(c))) * emv(p) * sb cir(p) = emv(p)*emg(c)*sb ! Saturated vapor pressure, specific humidity, and their derivatives ! at the leaf surface call QSat (t_veg(p), forc_pbot(g), el(p), deldT, qsatl(p), qsatldT(p)) ! Determine atmospheric co2 and o2 co2(p) = forc_pco2(g) o2(p) = forc_po2(g) #if (defined C13) ! 4/14/05: PET ! Adding isotope code c13o2(p) = forc_pc13o2(g) #endif ! Initialize flux profile nmozsgn(p) = 0 taf(p) = (t_grnd(c) + thm(p))/2._r8 qaf(p) = (forc_q(g)+qg(c))/2._r8 ur(p) = max(1.0_r8,sqrt(forc_u(g)*forc_u(g)+forc_v(g)*forc_v(g))) dth(p) = thm(p)-taf(p) dqh(p) = forc_q(g)-qaf(p) delq(p) = qg(c) - qaf(p) dthv(p) = dth(p)*(1._r8+0.61_r8*forc_q(g))+0.61_r8*forc_th(g)*dqh(p) zldis(p) = forc_hgt_u_pft(p) - displa(p) ! Check to see if the forcing height is below the canopy height if (zldis(p) < 0._r8) then found = .true. index = p end if end do if (found) then write(iulog,*)'Error: Forcing height is below canopy height for pft index ',index call endrun() end if do f = 1, fn p = filterp(f) c = pcolumn(p) ! Initialize Monin-Obukhov length and wind speed call MoninObukIni(ur(p), thv(c), dthv(p), zldis(p), z0mv(p), um(p), obu(p)) end do ! Set counter for leaf temperature iteration (itlef) itlef = 0 fnorig = fn fporig(1:fn) = filterp(1:fn) ! Make copies so that array sections are not passed in function calls to friction velocity do f = 1, fn p = filterp(f) displa_loc(p) = displa(p) z0mv_loc(p) = z0mv(p) z0hv_loc(p) = z0hv(p) z0qv_loc(p) = z0qv(p) end do ! Begin stability iteration ITERATION : do while (itlef <= itmax .and. fn > 0) ! Determine friction velocity, and potential temperature and humidity ! profiles of the surface boundary layer call FrictionVelocity (lbp, ubp, fn, filterp, & displa_loc, z0mv_loc, z0hv_loc, z0qv_loc, & obu, itlef+1, ur, um, ustar, & temp1, temp2, temp12m, temp22m, fm) do f = 1, fn p = filterp(f) c = pcolumn(p) g = pgridcell(p) tlbef(p) = t_veg(p) del2(p) = del(p) ! Determine aerodynamic resistances ram1(p) = 1._r8/(ustar(p)*ustar(p)/um(p)) rah(p,1) = 1._r8/(temp1(p)*ustar(p)) raw(p,1) = 1._r8/(temp2(p)*ustar(p)) ! Bulk boundary layer resistance of leaves uaf(p) = um(p)*sqrt( 1._r8/(ram1(p)*um(p)) ) cf = 0.01_r8/(sqrt(uaf(p))*sqrt(dleaf(ivt(p)))) rb(p) = 1._r8/(cf*uaf(p)) rb1(p) = rb(p) ! Parameterization for variation of csoilc with canopy density from ! X. Zeng, University of Arizona w = exp(-(elai(p)+esai(p))) ! changed by K.Sakaguchi from here ! transfer coefficient over bare soil is changed to a local variable ! just for readability of the code (from line 680) csoilb = (vkc/(0.13_r8*(z0mg(c)*uaf(p)/1.5e-5_r8)**0.45_r8)) !compute the stability parameter for ricsoilc ("S" in Sakaguchi&Zeng,2008) ri = ( grav*htop(p) * (taf(p) - t_grnd(c)) ) / (taf(p) * uaf(p) **2.00_r8) !! modify csoilc value (0.004) if the under-canopy is in stable condition if ( (taf(p) - t_grnd(c) ) > 0._r8) then ! decrease the value of csoilc by dividing it with (1+gamma*min(S, 10.0)) ! ria ("gmanna" in Sakaguchi&Zeng, 2008) is a constant (=0.5) ricsoilc = csoilc / (1.00_r8 + ria*min( ri, 10.0_r8) ) csoilcn = csoilb*w + ricsoilc*(1._r8-w) else csoilcn = csoilb*w + csoilc*(1._r8-w) end if !! Sakaguchi changes for stability formulation ends here rah(p,2) = 1._r8/(csoilcn*uaf(p)) raw(p,2) = rah(p,2) ! Stomatal resistances for sunlit and shaded fractions of canopy. ! Done each iteration to account for differences in eah, tv. svpts(p) = el(p) ! pa eah(p) = forc_pbot(g) * qaf(p) / 0.622_r8 ! pa end do ! 4/25/05, PET: Now calling the sun/shade version of Stomata by default call Stomata (fn, filterp, lbp, ubp, svpts, eah, o2, co2, rb, dayl_factor, phase='sun') call Stomata (fn, filterp, lbp, ubp, svpts, eah, o2, co2, rb, dayl_factor, phase='sha') do f = 1, fn p = filterp(f) c = pcolumn(p) g = pgridcell(p) ! Sensible heat conductance for air, leaf and ground ! Moved the original subroutine in-line... wta = 1._r8/rah(p,1) ! air wtl = (elai(p)+esai(p))/rb(p) ! leaf wtg(p) = 1._r8/rah(p,2) ! ground wtshi = 1._r8/(wta+wtl+wtg(p)) wtl0(p) = wtl*wtshi ! leaf wtg0 = wtg(p)*wtshi ! ground wta0(p) = wta*wtshi ! air wtga = wta0(p)+wtg0 ! ground + air wtal(p) = wta0(p)+wtl0(p) ! air + leaf ! Fraction of potential evaporation from leaf if (fdry(p) .gt. 0._r8) then rppdry = fdry(p)*rb(p)*(laisun(p)/(rb(p)+rssun(p)) + & laisha(p)/(rb(p)+rssha(p)))/elai(p) else rppdry = 0._r8 end if efpot = forc_rho(g)*wtl*(qsatl(p)-qaf(p)) if (efpot > 0._r8) then if (btran(p) > btran0) then qflx_tran_veg(p) = efpot*rppdry rpp = rppdry + fwet(p) else !No transpiration if btran below 1.e-10 rpp = fwet(p) qflx_tran_veg(p) = 0._r8 end if !Check total evapotranspiration from leaves rpp = min(rpp, (qflx_tran_veg(p)+h2ocan(p)/dtime)/efpot) else !No transpiration if potential evaporation less than zero rpp = 1._r8 qflx_tran_veg(p) = 0._r8 end if ! Update conductances for changes in rpp ! Latent heat conductances for ground and leaf. ! Air has same conductance for both sensible and latent heat. ! Moved the original subroutine in-line... wtaq = frac_veg_nosno(p)/raw(p,1) ! air wtlq = frac_veg_nosno(p)*(elai(p)+esai(p))/rb(p) * rpp ! leaf !Litter layer resistance. Added by K.Sakaguchi snowdp_c = z_dl ! critical depth for 100% litter burial by snow (=litter thickness) fsno_dl = snowdp(c)/snowdp_c ! effective snow cover for (dry)plant litter elai_dl = lai_dl*(1._r8 - min(fsno_dl,1._r8)) ! exposed (dry)litter area index rdl = ( 1._r8 - exp(-elai_dl) ) / ( 0.004_r8*uaf(p)) ! dry litter layer resistance ! add litter resistance and Lee and Pielke 1992 beta if (delq(p) .lt. 0._r8) then !dew. Do not apply beta for negative flux (follow old rsoil) wtgq(p) = frac_veg_nosno(p)/(raw(p,2)+rdl) else wtgq(p) = soilbeta(c)*frac_veg_nosno(p)/(raw(p,2)+rdl) end if wtsqi = 1._r8/(wtaq+wtlq+wtgq(p)) wtgq0 = wtgq(p)*wtsqi ! ground wtlq0(p) = wtlq*wtsqi ! leaf wtaq0(p) = wtaq*wtsqi ! air wtgaq = wtaq0(p)+wtgq0 ! air + ground wtalq(p) = wtaq0(p)+wtlq0(p) ! air + leaf dc1 = forc_rho(g)*cpair*wtl dc2 = hvap*forc_rho(g)*wtlq efsh = dc1*(wtga*t_veg(p)-wtg0*t_grnd(c)-wta0(p)*thm(p)) efe(p) = dc2*(wtgaq*qsatl(p)-wtgq0*qg(c)-wtaq0(p)*forc_q(g)) ! Evaporation flux from foliage erre = 0._r8 if (efe(p)*efeb(p) < 0._r8) then efeold = efe(p) efe(p) = 0.1_r8*efeold erre = efe(p) - efeold end if dt_veg(p) = (sabv(p) + air(p) + bir(p)*t_veg(p)**4 + & cir(p)*t_grnd(c)**4 - efsh - efe(p)) / & (- 4._r8*bir(p)*t_veg(p)**3 +dc1*wtga +dc2*wtgaq*qsatldT(p)) t_veg(p) = tlbef(p) + dt_veg(p) dels = dt_veg(p) del(p) = abs(dels) err(p) = 0._r8 if (del(p) > delmax) then dt_veg(p) = delmax*dels/del(p) t_veg(p) = tlbef(p) + dt_veg(p) err(p) = sabv(p) + air(p) + bir(p)*tlbef(p)**3*(tlbef(p) + & 4._r8*dt_veg(p)) + cir(p)*t_grnd(c)**4 - & (efsh + dc1*wtga*dt_veg(p)) - (efe(p) + & dc2*wtgaq*qsatldT(p)*dt_veg(p)) end if ! Fluxes from leaves to canopy space ! "efe" was limited as its sign changes frequently. This limit may ! result in an imbalance in "hvap*qflx_evap_veg" and ! "efe + dc2*wtgaq*qsatdt_veg" efpot = forc_rho(g)*wtl*(wtgaq*(qsatl(p)+qsatldT(p)*dt_veg(p)) & -wtgq0*qg(c)-wtaq0(p)*forc_q(g)) qflx_evap_veg(p) = rpp*efpot ! Calculation of evaporative potentials (efpot) and ! interception losses; flux in kg m**-2 s-1. ecidif ! holds the excess energy if all intercepted water is evaporated ! during the timestep. This energy is later added to the ! sensible heat flux. ecidif = 0._r8 if (efpot > 0._r8 .and. btran(p) > btran0) then qflx_tran_veg(p) = efpot*rppdry else qflx_tran_veg(p) = 0._r8 end if ecidif = max(0._r8, qflx_evap_veg(p)-qflx_tran_veg(p)-h2ocan(p)/dtime) qflx_evap_veg(p) = min(qflx_evap_veg(p),qflx_tran_veg(p)+h2ocan(p)/dtime) ! The energy loss due to above two limits is added to ! the sensible heat flux. eflx_sh_veg(p) = efsh + dc1*wtga*dt_veg(p) + err(p) + erre + hvap*ecidif ! Re-calculate saturated vapor pressure, specific humidity, and their ! derivatives at the leaf surface call QSat(t_veg(p), forc_pbot(g), el(p), deldT, qsatl(p), qsatldT(p)) ! Update vegetation/ground surface temperature, canopy air ! temperature, canopy vapor pressure, aerodynamic temperature, and ! Monin-Obukhov stability parameter for next iteration. taf(p) = wtg0*t_grnd(c) + wta0(p)*thm(p) + wtl0(p)*t_veg(p) qaf(p) = wtlq0(p)*qsatl(p) + wtgq0*qg(c) + forc_q(g)*wtaq0(p) ! Update Monin-Obukhov length and wind speed including the ! stability effect dth(p) = thm(p)-taf(p) dqh(p) = forc_q(g)-qaf(p) delq(p) = wtalq(p)*qg(c)-wtlq0(p)*qsatl(p)-wtaq0(p)*forc_q(g) tstar = temp1(p)*dth(p) qstar = temp2(p)*dqh(p) thvstar = tstar*(1._r8+0.61_r8*forc_q(g)) + 0.61_r8*forc_th(g)*qstar zeta = zldis(p)*vkc*grav*thvstar/(ustar(p)**2*thv(c)) if (zeta >= 0._r8) then !stable zeta = min(2._r8,max(zeta,0.01_r8)) um(p) = max(ur(p),0.1_r8) else !unstable zeta = max(-100._r8,min(zeta,-0.01_r8)) wc = beta*(-grav*ustar(p)*thvstar*zii/thv(c))**0.333_r8 um(p) = sqrt(ur(p)*ur(p)+wc*wc) end if obu(p) = zldis(p)/zeta if (obuold(p)*obu(p) < 0._r8) nmozsgn(p) = nmozsgn(p)+1 if (nmozsgn(p) >= 4) obu(p) = zldis(p)/(-0.01_r8) obuold(p) = obu(p) end do ! end of filtered pft loop ! Test for convergence itlef = itlef+1 if (itlef > itmin) then do f = 1, fn p = filterp(f) dele(p) = abs(efe(p)-efeb(p)) efeb(p) = efe(p) det(p) = max(del(p),del2(p)) end do fnold = fn fn = 0 do f = 1, fnold p = filterp(f) if (.not. (det(p) < dtmin .and. dele(p) < dlemin)) then fn = fn + 1 filterp(fn) = p end if end do end if end do ITERATION ! End stability iteration fn = fnorig filterp(1:fn) = fporig(1:fn) do f = 1, fn p = filterp(f) c = pcolumn(p) g = pgridcell(p) ! Energy balance check in canopy err(p) = sabv(p) + air(p) + bir(p)*tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p)) & + cir(p)*t_grnd(c)**4 - eflx_sh_veg(p) - hvap*qflx_evap_veg(p) ! Fluxes from ground to canopy space delt = wtal(p)*t_grnd(c)-wtl0(p)*t_veg(p)-wta0(p)*thm(p) taux(p) = -forc_rho(g)*forc_u(g)/ram1(p) tauy(p) = -forc_rho(g)*forc_v(g)/ram1(p) eflx_sh_grnd(p) = cpair*forc_rho(g)*wtg(p)*delt qflx_evap_soi(p) = forc_rho(g)*wtgq(p)*delq(p) ! 2 m height air temperature t_ref2m(p) = thm(p) + temp1(p)*dth(p)*(1._r8/temp12m(p) - 1._r8/temp1(p)) t_ref2m_r(p) = t_ref2m(p) ! 2 m height specific humidity q_ref2m(p) = forc_q(g) + temp2(p)*dqh(p)*(1._r8/temp22m(p) - 1._r8/temp2(p)) ! 2 m height relative humidity call QSat(t_ref2m(p), forc_pbot(g), e_ref2m, de2mdT, qsat_ref2m, dqsat2mdT) rh_ref2m(p) = min(100._r8, q_ref2m(p) / qsat_ref2m * 100._r8) rh_ref2m_r(p) = rh_ref2m(p) ! Downward longwave radiation below the canopy dlrad(p) = (1._r8-emv(p))*emg(c)*forc_lwrad(g) + & emv(p)*emg(c)*sb*tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p)) ! Upward longwave radiation above the canopy ulrad(p) = ((1._r8-emg(c))*(1._r8-emv(p))*(1._r8-emv(p))*forc_lwrad(g) & + emv(p)*(1._r8+(1._r8-emg(c))*(1._r8-emv(p)))*sb*tlbef(p)**3*(tlbef(p) + & 4._r8*dt_veg(p)) + emg(c)*(1._r8-emv(p))*sb*t_grnd(c)**4) ! Derivative of soil energy flux with respect to soil temperature cgrnds(p) = cgrnds(p) + cpair*forc_rho(g)*wtg(p)*wtal(p) cgrndl(p) = cgrndl(p) + forc_rho(g)*wtgq(p)*wtalq(p)*dqgdT(c) cgrnd(p) = cgrnds(p) + cgrndl(p)*htvp(c) ! Update dew accumulation (kg/m2) h2ocan(p) = max(0._r8,h2ocan(p)+(qflx_tran_veg(p)-qflx_evap_veg(p))*dtime) ! total photosynthesis fpsn(p) = psnsun(p)*laisun(p) + psnsha(p)*laisha(p) #if (defined CN) && (defined C13) ! 4/14/05: PET ! Adding isotope code rc13_canair(p) = c13o2(p)/(co2(p)-c13o2(p)) rc13_psnsun(p) = rc13_canair(p)/alphapsnsun(p) rc13_psnsha(p) = rc13_canair(p)/alphapsnsha(p) c13_psnsun(p) = psnsun(p) * (rc13_psnsun(p)/(1._r8+rc13_psnsun(p))) c13_psnsha(p) = psnsha(p) * (rc13_psnsha(p)/(1._r8+rc13_psnsha(p))) !write(iulog,*) p,ivt(p),btran(p),psnsun(p),psnsha(p),alphapsnsun(p),alphapsnsha(p) #endif end do ! Filter out pfts which have small energy balance errors; report others fnold = fn fn = 0 do f = 1, fnold p = filterp(f) if (abs(err(p)) > 0.1_r8) then fn = fn + 1 filterp(fn) = p end if end do do f = 1, fn p = filterp(f) write(iulog,*) 'energy balance in canopy ',p,', err=',err(p) end do end subroutine CanopyFluxes !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: Stomata ! ! !INTERFACE: subroutine Stomata (fn, filterp, lbp, ubp, ei, ea, o2, co2, rb, dayl_factor, phase) ! ! !DESCRIPTION: ! Leaf stomatal resistance and leaf photosynthesis. Modifications for CN code. ! !REVISION HISTORY: ! 22 January 2004: Created by Peter Thornton ! 4/14/05: Peter Thornton: Converted Ci from local variable to pps struct member ! now returns cisun or cisha per pft as implicit output argument. ! Also sets alphapsnsun and alphapsnsha. ! 4/25/05, Peter Thornton: Adopted as the default code for CLM, together with ! modifications for sun/shade canopy. Renamed from StomataCN to Stomata, ! and eliminating the older Stomata subroutine ! 3/6/09: Peter Thornton; added dayl_factor control on Vcmax, from Bill Bauerle ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8 use shr_const_mod, only : SHR_CONST_TKFRZ, SHR_CONST_RGAS use clmtype use clm_atmlnd , only : clm_a2l use spmdMod , only: masterproc use pftvarcon , only : nbrdlf_dcd_tmp_shrub use pftvarcon , only : nsoybean, npcropmin ! ! !ARGUMENTS: implicit none integer , intent(in) :: fn ! size of pft filter integer , intent(in) :: filterp(fn) ! pft filter integer , intent(in) :: lbp, ubp ! pft bounds real(r8), intent(in) :: ei(lbp:ubp) ! vapor pressure inside leaf (sat vapor press at tl) (pa) real(r8), intent(in) :: ea(lbp:ubp) ! vapor pressure of canopy air (pa) real(r8), intent(in) :: o2(lbp:ubp) ! atmospheric o2 concentration (pa) real(r8), intent(in) :: co2(lbp:ubp) ! atmospheric co2 concentration (pa) real(r8), intent(inout) :: rb(lbp:ubp) ! boundary layer resistance (s/m) real(r8), intent(in) :: dayl_factor(lbp:ubp) ! scalar (0-1) for daylength character(len=*), intent(in) :: phase ! 'sun' or 'sha' ! ! !CALLED FROM: ! subroutine CanopyFluxes in this module ! ! !LOCAL VARIABLES: ! ! local pointers to implicit in variables ! new ecophys variables (leafcn, flnr) added 1/26/04 ! integer , pointer :: pcolumn(:) ! pft's column index integer , pointer :: pgridcell(:) ! pft's gridcell index integer , pointer :: ivt(:) ! pft vegetation type real(r8), pointer :: qe25(:) ! quantum efficiency at 25C (umol CO2 / umol photon) real(r8), pointer :: c3psn(:) ! photosynthetic pathway: 0. = c4, 1. = c3 real(r8), pointer :: mp(:) ! slope of conductance-to-photosynthesis relationship real(r8), pointer :: tgcm(:) ! air temperature at agcm reference height (kelvin) real(r8), pointer :: forc_pbot(:) ! atmospheric pressure (Pa) real(r8), pointer :: tl(:) ! leaf temperature (Kelvin) real(r8), pointer :: btran(:) ! soil water transpiration factor (0 to 1) real(r8), pointer :: apar(:) ! par absorbed per unit lai (w/m**2) real(r8), pointer :: leafcn(:) ! leaf C:N (gC/gN) real(r8), pointer :: flnr(:) ! fraction of leaf N in the Rubisco enzyme (gN Rubisco / gN leaf) real(r8), pointer :: sla(:) ! specific leaf area, projected area basis (m^2/gC) real(r8), pointer :: fnitr(:) ! foliage nitrogen limitation factor (-) ! ! local pointers to implicit inout variables ! real(r8), pointer :: rs(:) ! leaf stomatal resistance (s/m) real(r8), pointer :: psn(:) ! foliage photosynthesis (umol co2 /m**2/ s) [always +] real(r8), pointer :: ci(:) ! intracellular leaf CO2 (Pa) #if (defined C13) real(r8), pointer :: alphapsn(:) ! 13C fractionation factor for PSN () #endif ! ! local pointers to implicit out variables ! real(r8), pointer :: lnc(:) ! leaf N concentration per unit projected LAI (gN leaf/m^2) real(r8), pointer :: vcmx(:) ! maximum rate of carboxylation (umol co2/m**2/s) ! ! ! !LOCAL VARIABLES: !EOP ! real(r8), parameter :: mpe = 1.e-6_r8 ! prevents overflow error if division by zero integer , parameter :: niter = 3 ! number of iterations integer :: f,p,c,g ! indices integer :: iter ! iteration index real(r8) :: ab ! used in statement functions real(r8) :: bc ! used in statement functions real(r8) :: f1 ! generic temperature response (statement function) real(r8) :: f2 ! generic temperature inhibition (statement function) real(r8) :: tc ! leaf temperature (degree celsius) real(r8) :: cs ! co2 concentration at leaf surface (pa) real(r8) :: kc ! co2 michaelis-menten constant (pa) real(r8) :: ko ! o2 michaelis-menten constant (pa) real(r8) :: atmp ! intermediate calculations for rs real(r8) :: btmp ! intermediate calculations for rs real(r8) :: ctmp ! intermediate calculations for rs real(r8) :: q ! intermediate calculations for rs real(r8) :: r1,r2 ! roots for rs real(r8) :: ppf ! absorb photosynthetic photon flux (umol photons/m**2/s) real(r8) :: wc ! rubisco limited photosynthesis (umol co2/m**2/s) real(r8) :: wj ! light limited photosynthesis (umol co2/m**2/s) real(r8) :: we ! export limited photosynthesis (umol co2/m**2/s) real(r8) :: cp ! co2 compensation point (pa) real(r8) :: awc ! intermediate calcuation for wc real(r8) :: j ! electron transport (umol co2/m**2/s) real(r8) :: cea ! constrain ea or else model blows up real(r8) :: cf ! s m**2/umol -> s/m real(r8) :: rsmax0 ! maximum stomatal resistance [s/m] real(r8) :: kc25 ! co2 michaelis-menten constant at 25c (pa) real(r8) :: akc ! q10 for kc25 real(r8) :: ko25 ! o2 michaelis-menten constant at 25c (pa) real(r8) :: ako ! q10 for ko25 real(r8) :: bp ! minimum leaf conductance (umol/m**2/s) ! additional variables for new treatment of Vcmax, Peter Thornton, 1/26/04 real(r8) :: act25 ! (umol/mgRubisco/min) Rubisco activity at 25 C real(r8) :: act ! (umol/mgRubisco/min) Rubisco activity real(r8) :: q10act ! (DIM) Q_10 for Rubisco activity real(r8) :: fnr ! (gRubisco/gN in Rubisco) !------------------------------------------------------------------------------ ! Set statement functions f1(ab,bc) = ab**((bc-25._r8)/10._r8) f2(ab) = 1._r8 + exp((-2.2e05_r8+710._r8*(ab+SHR_CONST_TKFRZ))/(SHR_CONST_RGAS*0.001_r8*(ab+SHR_CONST_TKFRZ))) ! Assign local pointers to derived type members (pft-level) pcolumn => clm3%g%l%c%p%column pgridcell => clm3%g%l%c%p%gridcell ivt => clm3%g%l%c%p%itype tl => clm3%g%l%c%p%pes%t_veg btran => clm3%g%l%c%p%pps%btran if (phase == 'sun') then apar => clm3%g%l%c%p%pef%parsun rs => clm3%g%l%c%p%pps%rssun psn => clm3%g%l%c%p%pcf%psnsun ci => clm3%g%l%c%p%pps%cisun #if (defined C13) alphapsn => clm3%g%l%c%p%pps%alphapsnsun #endif sla => clm3%g%l%c%p%pps%slasun lnc => clm3%g%l%c%p%pps%lncsun vcmx => clm3%g%l%c%p%pps%vcmxsun else if (phase == 'sha') then apar => clm3%g%l%c%p%pef%parsha rs => clm3%g%l%c%p%pps%rssha psn => clm3%g%l%c%p%pcf%psnsha ci => clm3%g%l%c%p%pps%cisha sla => clm3%g%l%c%p%pps%slasha #if (defined C13) alphapsn => clm3%g%l%c%p%pps%alphapsnsha #endif lnc => clm3%g%l%c%p%pps%lncsha vcmx => clm3%g%l%c%p%pps%vcmxsha end if ! Assign local pointers to derived type members (gridcell-level) forc_pbot => clm_a2l%forc_pbot ! Assign local pointers to derived type members (column-level) tgcm => clm3%g%l%c%p%pes%thm ! Assign local pointers to pft constants ! new ecophys constants added 1/26/04 qe25 => pftcon%qe25 c3psn => pftcon%c3psn mp => pftcon%mp leafcn => pftcon%leafcn flnr => pftcon%flnr fnitr => pftcon%fnitr ! Set constant values kc25 = 30._r8 akc = 2.1_r8 ko25 = 30000._r8 ako = 1.2_r8 bp = 2000._r8 ! New constants for CN code, added 1/26/04 act25 = 3.6_r8 q10act = 2.4_r8 fnr = 7.16_r8 ! Convert rubisco activity units from umol/mgRubisco/min -> umol/gRubisco/s act25 = act25 * 1000.0_r8 / 60.0_r8 do f = 1, fn p = filterp(f) c = pcolumn(p) g = pgridcell(p) ! Initialize rs=rsmax and psn=0 because calculations are performed only ! when apar > 0, in which case rs <= rsmax and psn >= 0 ! Set constants rsmax0 = 2.e4_r8 cf = forc_pbot(g)/(SHR_CONST_RGAS*0.001_r8*tgcm(p))*1.e06_r8 if (apar(p) <= 0._r8) then ! night time rs(p) = min(rsmax0, 1._r8/bp * cf) psn(p) = 0._r8 lnc(p) = 0._r8 vcmx(p) = 0._r8 #if (defined C13) alphapsn(p) = 1._r8 #endif else ! day time tc = tl(p) - SHR_CONST_TKFRZ ppf = 4.6_r8 * apar(p) j = ppf * qe25(ivt(p)) kc = kc25 * f1(akc,tc) ko = ko25 * f1(ako,tc) awc = kc * (1._r8+o2(p)/ko) cp = 0.5_r8*kc/ko*o2(p)*0.21_r8 ! Modification for shrubs proposed by X.D.Z ! Why does he prefer this line here instead of in subr. ! CanopyFluxes? (slevis) ! Equivalent modification for soy following AgroIBIS #if (defined CNDV) if (ivt(p) == nbrdlf_dcd_tmp_shrub) btran(p) = min(1._r8, btran(p) * 3.33_r8) #endif if (ivt(p) == nsoybean) btran(p) = min(1._r8, btran(p) * 1.25_r8) ! new calculations for vcmax, 1/26/04 lnc(p) = 1._r8 / (sla(p) * leafcn(ivt(p))) act = act25 * f1(q10act,tc) #if (defined CN) if ( ivt(p) < npcropmin )then vcmx(p) = lnc(p) * flnr(ivt(p)) * fnr * act / f2(tc) * btran(p) * & dayl_factor(p) else vcmx(p) = 101._r8 * f1(q10act,tc) / f2(tc) * btran(p) * dayl_factor(p) end if #else vcmx(p) = lnc(p) * flnr(ivt(p)) * fnr * act / f2(tc) * btran(p) * & dayl_factor(p) * fnitr(ivt(p)) #endif ! First guess ci ci(p) = 0.7_r8*co2(p)*c3psn(ivt(p)) + 0.4_r8*co2(p)*(1._r8-c3psn(ivt(p))) ! rb: s/m -> s m**2 / umol rb(p) = rb(p)/cf ! Constrain ea cea = max(0.25_r8*ei(p)*c3psn(ivt(p))+0.40_r8*ei(p)*(1._r8-c3psn(ivt(p))), min(ea(p),ei(p)) ) ! ci iteration for 'actual' photosynthesis do iter = 1, niter wj = max(ci(p)-cp,0._r8)*j/(ci(p)+2._r8*cp)*c3psn(ivt(p)) + j*(1._r8-c3psn(ivt(p))) wc = max(ci(p)-cp,0._r8)*vcmx(p)/(ci(p)+awc)*c3psn(ivt(p)) + vcmx(p)*(1._r8-c3psn(ivt(p))) we = 0.5_r8*vcmx(p)*c3psn(ivt(p)) + 4000._r8*vcmx(p)*ci(p)/forc_pbot(g)*(1._r8-c3psn(ivt(p))) psn(p) = min(wj,wc,we) cs = max( co2(p)-1.37_r8*rb(p)*forc_pbot(g)*psn(p), mpe ) atmp = mp(ivt(p))*psn(p)*forc_pbot(g)*cea / (cs*ei(p)) + bp btmp = ( mp(ivt(p))*psn(p)*forc_pbot(g)/cs + bp ) * rb(p) - 1._r8 ctmp = -rb(p) if (btmp >= 0._r8) then q = -0.5_r8*( btmp + sqrt(btmp*btmp-4._r8*atmp*ctmp) ) else q = -0.5_r8*( btmp - sqrt(btmp*btmp-4._r8*atmp*ctmp) ) end if r1 = q/atmp r2 = ctmp/q rs(p) = max(r1,r2) ci(p) = max( cs-psn(p)*forc_pbot(g)*1.65_r8*rs(p), 0._r8 ) end do ! rs, rb: s m**2 / umol -> s/m rs(p) = min(rsmax0, rs(p)*cf) rb(p) = rb(p) * cf #if (defined C13) ! 4/14/05: PET ! Adding isotope code alphapsn(p) = 1._r8 + (((c3psn(ivt(p)) * (4.4_r8 + (22.6_r8*(ci(p)/co2(p))))) + & ((1._r8 - c3psn(ivt(p))) * 4.4_r8))/1000._r8) !alphapsn(p) = 1._r8 !write(iulog,*) 'in StomataCN ',p,ivt(p),c3psn(ivt(p)),ci(p),co2(p),alphapsn(p) #endif end if end do end subroutine Stomata end module CanopyFluxesMod