#include <misc.h>
#include <preproc.h>
module CanopyFluxesMod 1,3
!------------------------------------------------------------------------------
!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, & 1,18
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
use clm_varpar , only : nlevgrnd, nlevsno
use clm_varcon , only : sb, cpair, hvap, vkc, grav, denice, &
denh2o, tfrz, csoilc, tlsai_crit, alpha_aero
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)
!
! 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)
!
!
! !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 [-]
!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
!------------------------------------------------------------------------------
! 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
! 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
! Assign local pointers to derived type members (ecophysiological)
dleaf => pftcon%dleaf
smpso => pftcon%smpso
smpsc => pftcon%smpsc
! Determine step size
dtime = get_step_size()
! 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))))
#if (defined NO_DAYLEN_VCMAX)
dayl_factor(p) = 1.0_r8
#endif
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
!dir$ concurrent
!cdir nodep
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
!dir$ concurrent
!cdir nodep
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
! 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.
!dir$ concurrent
!cdir nodep
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
!dir$ concurrent
!cdir nodep
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
!dir$ concurrent
!cdir nodep
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)
!dir$ concurrent
!cdir nodep
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')
!dir$ concurrent
!cdir nodep
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
!dir$ concurrent
!cdir nodep
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)
!dir$ concurrent
!cdir nodep
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) 2,6
!
! !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
!
! !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 :: vcmx25(:) ! max rate of carboxylation at 25C (umol CO2/m**2/s)
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) :: avcmx ! q10 for vcmx25
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
vcmx25 => pftcon%vcmx25
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
avcmx = 2.4_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
!dir$ concurrent
!cdir nodep
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)
#if (defined CNDV)
if (ivt(p) == nbrdlf_dcd_tmp_shrub) btran(p) = min(1._r8, btran(p) * 3.33_r8)
#endif
! new calculations for vcmax, 1/26/04
lnc(p) = 1._r8 / (sla(p) * leafcn(ivt(p)))
act = act25 * f1(q10act,tc)
#if (defined CN)
vcmx(p) = lnc(p) * flnr(ivt(p)) * fnr * act / f2(tc) * btran(p) * dayl_factor(p)
#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
#if (defined NEC_SX)
!NEC NEC NEC
!
! ITER = 1
!
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 )
!
! ITER = 2
!
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 )
!
! ITER = 3
!
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 )
!NEC NEC NEC
#else
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
#endif
! 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