#include <misc.h> #include <preproc.h> module VOCEmissionMod 1,2 !----------------------------------------------------------------------- !BOP ! ! !MODULE: VOCEmissionMod ! ! !DESCRIPTION: ! Volatile organic compound emission ! ! !USES: use clm_varctl, only: iulog use abortutils, only: endrun ! ! !PUBLIC TYPES: implicit none save ! ! !PUBLIC MEMBER FUNCTIONS: public :: VOCEmission ! ! !REVISION HISTORY: ! Created by Mariana Vertenstein ! !EOP !----------------------------------------------------------------------- contains !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: VOCEmission ! ! !INTERFACE: subroutine VOCEmission (lbp, ubp, num_soilp, filter_soilp ) 1,21 ! ! ! NEW DESCRIPTION ! Volatile organic compound emission ! This code simulates volatile organic compound emissions following: ! 1. Isoprene: Guenther et al., 2006 description of MEGAN emissions ! following equations 2-9, 16-17, 20 ! 2. Monoterpenes/OVOCs/ORVOCs/CO: algorithm presented in Guenther, A., ! 1999: Modeling Biogenic Volatile Organic Compound Emissions to the ! Atmosphere. In Reactive Hydrocarbons in the Atmosphere, Ch. 3 ! With updates from MEGAN online user's guide ! ( http://acd.ucar.edu/~guenther/MEGAN/MEGANusersguide.pdf) ! This model relies on the assumption that 90% of isoprene and monoterpene ! emissions originate from canopy foliage: ! E= epsilon * gamma * rho ! VOC flux (E) [ugC m-2 h-1] is calculated from baseline emission ! factors (epsilon) [ugC m-2 h-1] which are mapped for each PFT (isoprene) ! or constant for each PFT (others). Note that for constant EFs the units ! of [ugC g-1 h-1] must be multiplied by the source density factor. ! The emission activity factor (gamma) [unitless] for isoprene includes ! dependence on PPFT, temperature, LAI, leaf age and soil moisture. ! The canopy environment constant was calculated offline for CLM+CAM at ! standard conditions. ! The emission activity factor for the other emissions depends on temperature. ! We assume that the escape efficiency (rho) here is unity following ! Guenther et al., 2006. ! Subroutine written to operate at the patch level. ! IN FINAL IMPLEMENTATION, REMEMBER: ! 1. may wish to call this routine only as freq. as rad. calculations ! 2. may wish to place epsilon values directly in pft-physiology file ! Output: vocflx(nvoc) !VOC flux [ug C m-2 h-1] ! ! ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8 use clm_atmlnd , only : clm_a2l use clmtype use clm_varpar , only : nvoc, numpft use clm_varpar , only : nvoc, numpft use clm_atmlnd , only : clm_a2l use shr_const_mod, only : SHR_CONST_RGAS use clm_varcon , only : denice use clm_varpar , only : nlevsoi use pftvarcon , only : ndllf_evr_tmp_tree, ndllf_evr_brl_tree, & ndllf_dcd_brl_tree, nbrdlf_evr_trp_tree, & nbrdlf_evr_tmp_tree, nbrdlf_dcd_brl_shrub, & nbrdlf_dcd_trp_tree, nbrdlf_dcd_tmp_tree, & nbrdlf_dcd_brl_tree, nbrdlf_evr_shrub, & nc3_arctic_grass, nwheat, & ncorn, nc4_grass, & noveg ! ! !ARGUMENTS: implicit none integer, intent(in) :: lbp, ubp ! pft bounds integer, intent(in) :: num_soilp ! number of columns in soil pft filter integer, intent(in) :: filter_soilp(num_soilp) ! pft filter for soil ! ! !CALLED FROM: ! ! !REVISION HISTORY: ! Author: Sam Levis ! 2/1/02, Peter Thornton: migration to new data structure ! 4/15/06, Colette L. Heald: modify for updated MEGAN model (Guenther et al., 2006) ! ! !LOCAL VARIABLES: ! ! local pointers to implicit in arguments ! integer , pointer :: pgridcell(:) ! gridcell index of corresponding pft integer , pointer :: pcolumn(:) ! column index of corresponding pft integer , pointer :: ivt(:) ! pft vegetation type for current real(r8), pointer :: t_veg(:) ! pft vegetation temperature (Kelvin) real(r8), pointer :: fsun(:) ! sunlit fraction of canopy real(r8), pointer :: elai(:) ! one-sided leaf area index with burying by snow real(r8), pointer :: clayfrac(:) ! fraction of soil that is clay real(r8), pointer :: sandfrac(:) ! fraction of soil that is sand real(r8), pointer :: forc_solad(:,:) ! direct beam radiation (visible only) real(r8), pointer :: forc_solai(:,:) ! diffuse radiation (visible only) real(r8), pointer :: sla(:) ! specific leaf area [m2 leaf g-1 C] real(r8), pointer :: h2osoi_vol(:,:) ! volumetric soil water (m3/m3) real(r8), pointer :: h2osoi_ice(:,:) ! ice soil content (kg/m3) real(r8), pointer :: dz(:,:) ! depth of layer (m) real(r8), pointer :: coszen(:) ! cosine of solar zenith angle real(r8), pointer :: efisop(:,:) ! emission factors for isoprene for each pft [ug C m-2 h-1] real(r8), pointer :: elai_p(:) ! one-sided leaf area index from previous timestep real(r8), pointer :: t_veg24(:) ! avg pft vegetation temperature for last 24 hrs real(r8), pointer :: t_veg240(:) ! avg pft vegetation temperature for last 240 hrs real(r8), pointer :: fsun24(:) ! sunlit fraction of canopy last 24 hrs real(r8), pointer :: fsun240(:) ! sunlit fraction of canopy last 240 hrs real(r8), pointer :: forc_solad24(:) ! direct beam radiation last 24hrs (visible only) real(r8), pointer :: forc_solai24(:) ! diffuse radiation last 24hrs (visible only) real(r8), pointer :: forc_solad240(:) ! direct beam radiation last 240hrs (visible only) real(r8), pointer :: forc_solai240(:) ! diffuse radiation last 240hrs (visible only) real(r8), pointer :: bsw(:,:) ! Clapp and Hornberger "b" (nlevgrnd) real(r8), pointer :: watsat(:,:) ! volumetric soil water at saturation (porosity) (nlevgrnd) real(r8), pointer :: sucsat(:,:) ! minimum soil suction (mm) (nlevgrnd) real(r8), parameter :: smpmax = 2.57e5_r8 ! maximum soil matrix potential ! ! local pointers to original implicit out arrays ! real(r8), pointer :: vocflx(:,:) ! VOC flux [ug C m-2 h-1] real(r8), pointer :: vocflx_tot(:) ! VOC flux [ug C m-2 h-1] real(r8), pointer :: vocflx_1(:) ! VOC flux(1) [ug C m-2 h-1] real(r8), pointer :: vocflx_2(:) ! VOC flux(2) [ug C m-2 h-1] real(r8), pointer :: vocflx_3(:) ! VOC flux(3) [ug C m-2 h-1] real(r8), pointer :: vocflx_4(:) ! VOC flux(4) [ug C m-2 h-1] real(r8), pointer :: vocflx_5(:) ! VOC flux(5) [ug C m-2 h-1] real(r8), pointer :: Eopt_out(:) real(r8), pointer :: topt_out(:) real(r8), pointer :: alpha_out(:) real(r8), pointer :: cp_out(:) real(r8), pointer :: paru_out(:) real(r8), pointer :: par24u_out(:) real(r8), pointer :: par240u_out(:) real(r8), pointer :: para_out(:) real(r8), pointer :: par24a_out(:) real(r8), pointer :: par240a_out(:) real(r8), pointer :: gamma_out(:) real(r8), pointer :: gammaT_out(:) real(r8), pointer :: gammaP_out(:) real(r8), pointer :: gammaL_out(:) real(r8), pointer :: gammaA_out(:) real(r8), pointer :: gammaS_out(:) ! ! ! !OTHER LOCAL VARIABLES: !EOP ! integer :: fp,p,g,c,n,j ! indices integer :: ct_bad real(r8) :: epsilon(lbp:ubp) ! emission factor [ugC m-2 h-1] real(r8) :: par ! temporary real(r8) :: par24 ! temporary real(r8) :: par240 ! temporary real(r8) :: density ! source density factor [g dry wgt foliar mass/m2 ground] real(r8) :: gamma(lbp:ubp) ! activity factor (accounting for light, T, age, LAI conditions) real(r8) :: gamma_p ! activity factor for PPFD real(r8) :: gamma_l ! activity factor for PPFD & LAI real(r8) :: gamma_t ! activity factor for temperature real(r8) :: gamma_a ! activity factor for leaf age real(r8) :: gamma_sm ! activity factor for soil moisture real(r8) :: x ! temporary real(r8) :: Eopt ! temporary real(r8) :: topt ! temporary real(r8) :: cp ! temporary real(r8) :: alpha ! temporary real(r8) :: elai_prev ! lai for previous timestep real(r8) :: fnew, fgro, fmat, fsen ! fractions of leaves at different phenological stages real(r8) :: nl ! temporary number of soil levels real(r8) :: theta_ice ! water content in ice in m3/m3 real(r8) :: wilt ! wilting point in m3/m3 real(r8) :: theta1 ! temporary ! ! Constants ! real(r8), parameter :: R = SHR_CONST_RGAS*0.001_r8 ! univ. gas constant [J K-1 mol-1] real(r8), parameter :: scale_mw =0.882_r8 ! conversion factor for isoprene -> carbon real(r8), parameter :: alpha_fix = 0.001_r8 ! empirical coefficient real(r8), parameter :: cp_fix = 1.21_r8 ! empirical coefficient real(r8), parameter :: ct1 = 95.0_r8 ! empirical coefficient (70 in User's Guide) real(r8), parameter :: ct2 = 230.0_r8 ! empirical coefficient (200 in User's Guide) real(r8), parameter :: ct3 = 0.00831_r8 ! empirical coefficient (0.0083 in User's Guide) real(r8), parameter :: topt_fix = 317._r8 ! std temperature [K] real(r8), parameter :: Eopt_fix = 2.26_r8 ! empirical coefficient real(r8), parameter :: tstd = 303.15_r8 ! std temperature [K] real(r8), parameter :: bet = 0.09_r8 ! beta empirical coefficient [K-1] real(r8), parameter :: clai1 = 0.49_r8 ! empirical coefficient real(r8), parameter :: clai2 = 0.2_r8 ! empirical coefficient real(r8), parameter :: clai3 = 5.0_r8 ! empirical coefficient real(r8), parameter :: Anew = 0.01_r8 ! relative emission factor for new plants real(r8), parameter :: Agro = 0.5_r8 ! relative emission factor for new plants real(r8), parameter :: Amat = 1.0_r8 ! relative emission factor for new plants real(r8), parameter :: Asen = 0.33_r8 ! relative emission factor for new plants real(r8), parameter :: cce = 0.40_r8 ! factor to set emissions to unity @ std real(r8), parameter :: cce1 = 0.47_r8 ! same as Cce but for non-accumulated vars real(r8), parameter :: ca1 = 0.004_r8 ! empirical coefficent for alpha real(r8), parameter :: ca2 = 0.0005_r8 ! empirical coefficent for alpha real(r8), parameter :: ca3 = 0.0468_r8 ! empirical coefficent for cp real(r8), parameter :: par0_sun = 200._r8 ! std conditions for past 24 hrs [umol/m2/s] real(r8), parameter :: par0_shade = 50._r8 ! std conditions for past 24 hrs [umol/m2/s] real(r8), parameter :: co1 = 313._r8 ! empirical coefficient real(r8), parameter :: co2 = 0.6_r8 ! empirical coefficient real(r8), parameter :: co3 = 2.034_r8 ! empirical coefficient real(r8), parameter :: co4 = 0.05_r8 ! empirical coefficient real(r8), parameter :: tstd0 = 297_r8 ! std temperature [K] real(r8), parameter :: deltheta1=0.06_r8 ! empirical coefficient real(r8), parameter :: scaling_to_500_Tg = 5._r8/7._r8 ! J-F Larmaque's empirical scaling factor ! ! These are the values from version of genesis-ibis / 1000. ! CN calculates its own sla [m2 leaf g-1 C] ! Divide by 2 in the equation to get dry weight foliar mass from grams carbon ! real(r8) :: hardwire_sla(0:numpft) real(r8) :: slarea(lbp:ubp) ! Specific leaf areas [m2 leaf g-1 C] real(r8) :: hardwire_droot(0:numpft) ! Root depth [m] !----------------------------------------------------------------------- ! Assign local pointers to derived type members (gridcell-level) forc_solad => clm_a2l%forc_solad forc_solai => clm_a2l%forc_solai efisop => clm3%g%gve%efisop ! Assign local pointers to derived subtypes components (column-level) h2osoi_vol => clm3%g%l%c%cws%h2osoi_vol h2osoi_ice => clm3%g%l%c%cws%h2osoi_ice dz => clm3%g%l%c%cps%dz bsw => clm3%g%l%c%cps%bsw watsat => clm3%g%l%c%cps%watsat sucsat => clm3%g%l%c%cps%sucsat ! Assign local pointers to derived subtypes components (pft-level) pgridcell => clm3%g%l%c%p%gridcell pcolumn => clm3%g%l%c%p%column ivt => clm3%g%l%c%p%itype t_veg => clm3%g%l%c%p%pes%t_veg fsun => clm3%g%l%c%p%pps%fsun elai => clm3%g%l%c%p%pps%elai clayfrac => clm3%g%l%c%p%pps%clayfrac sandfrac => clm3%g%l%c%p%pps%sandfrac vocflx => clm3%g%l%c%p%pvf%vocflx vocflx_tot => clm3%g%l%c%p%pvf%vocflx_tot vocflx_1 => clm3%g%l%c%p%pvf%vocflx_1 vocflx_2 => clm3%g%l%c%p%pvf%vocflx_2 vocflx_3 => clm3%g%l%c%p%pvf%vocflx_3 vocflx_4 => clm3%g%l%c%p%pvf%vocflx_4 vocflx_5 => clm3%g%l%c%p%pvf%vocflx_5 Eopt_out => clm3%g%l%c%p%pvf%Eopt_out topt_out => clm3%g%l%c%p%pvf%topt_out alpha_out => clm3%g%l%c%p%pvf%alpha_out cp_out => clm3%g%l%c%p%pvf%cp_out paru_out => clm3%g%l%c%p%pvf%paru_out par24u_out => clm3%g%l%c%p%pvf%par24u_out par240u_out => clm3%g%l%c%p%pvf%par240u_out para_out => clm3%g%l%c%p%pvf%para_out par24a_out => clm3%g%l%c%p%pvf%par24a_out par240a_out => clm3%g%l%c%p%pvf%par240a_out gammaL_out => clm3%g%l%c%p%pvf%gammaL_out gammaT_out => clm3%g%l%c%p%pvf%gammaT_out gammaP_out => clm3%g%l%c%p%pvf%gammaP_out gammaA_out => clm3%g%l%c%p%pvf%gammaA_out gammaS_out => clm3%g%l%c%p%pvf%gammaS_out gamma_out => clm3%g%l%c%p%pvf%gamma_out sla => clm3%g%l%c%p%pps%slasha t_veg24 => clm3%g%l%c%p%pvs%t_veg24 t_veg240 => clm3%g%l%c%p%pvs%t_veg240 forc_solad24 => clm3%g%l%c%p%pvs%fsd24 forc_solad240 => clm3%g%l%c%p%pvs%fsd240 forc_solai24 => clm3%g%l%c%p%pvs%fsi24 forc_solai240 => clm3%g%l%c%p%pvs%fsi240 fsun24 => clm3%g%l%c%p%pvs%fsun24 fsun240 => clm3%g%l%c%p%pvs%fsun240 elai_p => clm3%g%l%c%p%pvs%elai_p hardwire_sla(noveg) = 0._r8 ! bare-soil hardwire_sla(ndllf_evr_tmp_tree) = 0.0125_r8 !needleleaf hardwire_sla(ndllf_evr_brl_tree) = 0.0125_r8 !Gordon Bonan suggests NET = 0.0076 hardwire_sla(ndllf_dcd_brl_tree) = 0.0125_r8 !Gordon Bonan suggests NDT = 0.0200 hardwire_sla(nbrdlf_evr_trp_tree) = 0.0250_r8 !broadleaf hardwire_sla(nbrdlf_evr_tmp_tree) = 0.0250_r8 !Gordon Bonan suggests BET = 0.0178 hardwire_sla(nbrdlf_dcd_trp_tree) = 0.0250_r8 !Gordon Bonan suggests BDT = 0.0274 hardwire_sla(nbrdlf_dcd_tmp_tree:nbrdlf_dcd_brl_shrub) = 0.0250_r8 hardwire_sla(nc3_arctic_grass:numpft) = 0.0200_r8 !grass/crop ! root depth (m) (defined based on Zeng et al., 2001, cf Guenther 2006) hardwire_droot(noveg) = 0._r8 ! bare-soil hardwire_droot(ndllf_evr_tmp_tree:ndllf_evr_brl_tree) = 1.8_r8 ! evergreen tree hardwire_droot(ndllf_dcd_brl_tree) = 2.0_r8 ! needleleaf deciduous boreal tree hardwire_droot(nbrdlf_evr_trp_tree:nbrdlf_evr_tmp_tree) = 3.0_r8 ! broadleaf evergreen tree hardwire_droot(nbrdlf_dcd_trp_tree:nbrdlf_dcd_brl_tree) = 2.0_r8 ! broadleaf deciduous tree hardwire_droot(nbrdlf_evr_shrub:nbrdlf_dcd_brl_shrub) = 2.5_r8 ! shrub hardwire_droot(nc3_arctic_grass:numpft) = 1.5_r8 ! grass/crop ! initialize variables which get passed to the atmosphere vocflx(lbp:ubp, :)=0._r8 ! Determine specific leaf array do fp = 1,num_soilp p = filter_soilp(fp) slarea(p) = hardwire_sla(ivt(p)) end do ! Begin loop through voc species !_______________________________________________________________________________ do n = 1, nvoc select case (n) case(1) do fp = 1,num_soilp p = filter_soilp(fp) g = pgridcell(p) ! epsilon: use gridded values for 6 PFTs specified by MEGAN following ! ------- Guenther et al. (2006). Map the numpft CLM PFTs to these 6. ! Units: [ug C m-2 h-1] (convert input files from units of ! [ug isop m-2 h-1]) epsilon(p) = 0._r8 ! isoprenes: if ( ivt(p) == ndllf_evr_tmp_tree & .or. ivt(p) == ndllf_evr_brl_tree) then !fineleaf evergreen epsilon(p) = efisop(2,g)*scale_mw else if (ivt(p) == ndllf_dcd_brl_tree) then !fineleaf deciduous epsilon(p) = efisop(3,g)*scale_mw else if (ivt(p) >= nbrdlf_evr_trp_tree & .and. ivt(p) <= nbrdlf_dcd_brl_tree) then !broadleaf trees epsilon(p) = efisop(1,g)*scale_mw else if (ivt(p) >= nbrdlf_evr_shrub & .and. ivt(p) <= nbrdlf_dcd_brl_shrub) then !shrubs epsilon(p) = efisop(4,g)*scale_mw else if (ivt(p) >= nc3_arctic_grass & .and. ivt(p) <= nc4_grass) then !grass epsilon(p) = efisop(5,g)*scale_mw else if (ivt(p) >= ncorn) then !crops epsilon(p) =efisop(6,g)*scale_mw end if end do case(2) do fp = 1,num_soilp p = filter_soilp(fp) g = pgridcell(p) ! epsilon: use values from table 3 in Guenther (1997) which originate in ! ------- Guenther et al. (1995). In the comments below, I mention the pft ! category as described in table 3. Some values were taken directly ! from Guenther et al. (1995). Units: [ugC g-1 h-1] ! Values were updated on 1/2002 (Guenther, personal communication) ! monoterpenes: epsilon(p) = 0._r8 ! monoterpenes: if ( ivt(p) >= ndllf_evr_tmp_tree & .and. ivt(p) <= ndllf_evr_brl_tree) then !needleleaf evergreen epsilon(p) = 2.0_r8 else if (ivt(p) == ndllf_dcd_brl_tree) then !needleleaf deciduous epsilon(p) = 1.6_r8 else if (ivt(p) >= nbrdlf_evr_trp_tree & .and. ivt(p) <= nbrdlf_dcd_brl_tree) then !broadleaf everg trop epsilon(p) = 0.4_r8 else if (ivt(p) >= nbrdlf_evr_shrub & .and. ivt(p) <= nbrdlf_dcd_brl_shrub) then !other woody veg epsilon(p) = 0.8_r8 else if (ivt(p) >= nc3_arctic_grass & .and. ivt(p) <= numpft) then !grass & crop epsilon(p) = 0.1_r8 end if end do case (3) do fp = 1,num_soilp p = filter_soilp(fp) g = pgridcell(p) ! other VOCs (OVOCs) epsilon(p) = 1.0_r8 !Guenther (personal communication) end do case (4) do fp = 1,num_soilp p = filter_soilp(fp) g = pgridcell(p) ! other reactive VOCs (ORVOCs) epsilon(p) = 1.0_r8 !Guenther (personal communication) end do case (5) do fp = 1,num_soilp p = filter_soilp(fp) g = pgridcell(p) ! CO epsilon(p) = 0.3_r8 !Guenther (personal communication) end do case default write(iulog,*)'only nvocs up to index 5 are currently supported' call endrun() end select ct_bad=0 select case (n) case (1) do fp = 1,num_soilp p = filter_soilp(fp) g = pgridcell(p) c = pcolumn(p) ! gamma: Activity factor. Units [dimensionless] ! ===== For isoprene include activity factors for LAI,PPFD, T, leaf age, and soil moisture ! Activity factor for LAI (Guenther et al., 2006) !------------------------ ! Guenther et al., 2006 eq 3 if ( (fsun240(p) > 0.0_r8) .and. (fsun240(p) < 1.e30_r8) ) then gamma_l = cce * elai(p) else gamma_l = cce1 * elai(p) end if gammaL_out(p)=gamma_l ! Activity factor for PPFD (Guenther et al., 2006) !------------------------- ! With distinction between sunlit and shaded leafs, weight scalings by ! fsun and fshade ! Scale total incident par by fraction of sunlit leaves (added on 1/2002) ! multiply w/m2 by 4.6 to get umol/m2/s for par (added 8/14/02) ! fvitt -- forc_solad240, forc_solai240 can be zero when CLM finidat is specified ! which will cause par240 to be zero and produce NaNs via log(par240) ! dml -- fsun240 can be equal to or greater than one before 10 day averages are ! set on startup or if a new pft comes online during land cover change. ! Avoid this problem by only doing calculations with fsun240 when fsun240 is ! between 0 and 1 if ( (fsun240(p) > 0._r8) .and. (fsun240(p) < 1._r8) .and. (forc_solad240(p) > 0._r8) & .and. (forc_solai240(p) > 0._r8)) then ! With alpha and cp calculated based on eq 6 and 7: ! Note indexing for accumulated variables is all at pft level ! SUN: par = (forc_solad(g,1) + fsun(p) * forc_solai(g,1)) * 4.6_r8 par24 = (forc_solad24(p) + fsun24(p) * forc_solai24(p)) * 4.6_r8 par240 = (forc_solad240(p) + fsun240(p) * forc_solai240(p)) * 4.6_r8 alpha = ca1 - ca2 * log(par240) cp = ca3 * exp(ca2 * (par24-par0_sun))*par240**(0.6_r8) gamma_p = fsun(p) * ( cp * alpha*par * (1._r8 + alpha*alpha*par*par)**(-0.5_r8) ) paru_out(p)=par par24u_out(p)=par24 par240u_out(p)=par240 ! SHADE: par = ((1._r8 - fsun(p)) * forc_solai(g,1)) * 4.6_r8 par24 = ((1._r8 - fsun24(p)) * forc_solai24(p)) * 4.6_r8 par240 = ((1._r8 - fsun240(p)) * forc_solai240(p)) * 4.6_r8 alpha = ca1 - ca2 * log(par240) cp = ca3 * exp(ca2 * (par24-par0_shade))*par240**(0.6_r8) par = ((1._r8 - fsun(p)) * forc_solai(g,1)) * 4.6_r8 gamma_p = gamma_p + (1-fsun(p)) * (cp*alpha*par*(1._r8 + alpha*alpha*par*par)**(-0.5_r8)) para_out(p)=par par24a_out(p)=par24 par240a_out(p)=par240 else ! With fixed alpha and cp (from MEGAN User's Guide): ! SUN: direct + diffuse par = (forc_solad(g,1) + fsun(p) * forc_solai(g,1)) * 4.6_r8 alpha = alpha_fix cp = cp_fix gamma_p = fsun(p) * ( cp * alpha*par * (1._r8 + alpha*alpha*par*par)**(-0.5_r8) ) paru_out(p)=par par24u_out(p)=-999 par240u_out(p)=-999 ! SHADE: diffuse par = ((1._r8 - fsun(p)) * forc_solai(g,1)) * 4.6_r8 gamma_p = gamma_p + (1-fsun(p)) * (cp*alpha*par*(1._r8 + alpha*alpha*par*par)**(-0.5_r8)) para_out(p)=par par24a_out(p)=-999 par240a_out(p)=-999 end if alpha_out(p)=alpha cp_out(p)=cp gammaP_out(p)=gamma_p ! Activity factor for temperature (Guenther et al., 2006) !-------------------------------- if ( (t_veg240(p) > 0.0_r8) .and. (t_veg240(p) < 1.e30_r8) ) then ! topt and Eopt from eq 8 and 9: topt = co1 + (co2 * (t_veg240(p)-tstd0)) Eopt = co3 * exp (co4 * (t_veg24(p)-tstd0)) * exp(co4 * (t_veg240(p) -tstd0)) else topt = topt_fix Eopt = Eopt_fix endif x = ( (1._r8/topt) - (1._r8/(t_veg(p))) ) / ct3 gamma_t = Eopt * ( ct2 * exp(ct1 * x)/(ct2 - ct1 * (1._r8 - exp(ct2 * x))) ) topt_out(p)=topt Eopt_out(p)=Eopt gammaT_out(p)=gamma_t ! Activity factor for leaf age (Guenther et al., 2006) !----------------------------- ! If not CNDV elai is constant therefore gamma_a=1.0 ! gamma_a set to unity for evergreens (PFTs 1, 2, 4, 5) ! Note that we assume here that the time step is shorter than the number of !days after budbreak required to induce isoprene emissions (ti=12 days) and ! the number of days after budbreak to reach peak emission (tm=28 days) if ( (ivt(p) == ndllf_dcd_brl_tree) .or. (ivt(p) >= nbrdlf_dcd_trp_tree) ) then ! non-evergreen if ( (elai_p(p) > 0.0_r8) .and. (elai_p(p) < 1.e30_r8) )then elai_prev = 2._r8*elai_p(p)-elai(p) ! have accumulated average lai over last timestep if (elai_prev == elai(p)) then fnew = 0.0_r8 fgro = 0.0_r8 fmat = 1.0_r8 fsen = 0.0_r8 else if (elai_prev > elai(p)) then fnew = 0.0_r8 fgro = 0.0_r8 fmat = 1.0_r8 - (elai_prev - elai(p))/elai_prev fsen = (elai_prev - elai(p))/elai_prev else if (elai_prev < elai(p)) then fnew = 1 - (elai_prev / elai(p)) fgro = 0.0_r8 fmat = (elai_prev / elai(p)) fsen = 0.0_r8 end if gamma_a = fnew * Anew + fgro * Agro + fmat * Amat + fsen * Asen else gamma_a = 1.0_r8 end if else gamma_a = 1.0_r8 end if gammaA_out(p)=gamma_a ! Activity factor for soil moisture (Guenther et al., 2006) !---------------------------------- ! Calculate the mean scaling factor throughout the root depth. ! wilting point potential is in units of matric potential (mm) ! (1 J/Kg = 0.001 MPa, approx = 0.1 m) ! convert to volumetric soil water using equation 7.118 of the CLM4 Technical Note if ((clayfrac(p) > 0) .and. (sandfrac(p) > 0)) then gamma_sm = 0._r8 nl=0._r8 do j = 1,nlevsoi if (sum(dz(c,1:j)) < hardwire_droot(ivt(p))) then theta_ice = h2osoi_ice(c,j)/(dz(c,j)*denice) wilt = ((smpmax/sucsat(c,j))**(-1._r8/bsw(c,j))) * (watsat(c,j) - theta_ice) theta1 = wilt + deltheta1 if (h2osoi_vol(c,j) >= theta1) then gamma_sm = gamma_sm + 1._r8 else if ( (h2osoi_vol(c,j) > wilt) .and. (h2osoi_vol(c,j) < theta1) ) then gamma_sm = gamma_sm + ( h2osoi_vol(c,j) - wilt ) / deltheta1 else gamma_sm = gamma_sm + 0._r8 end if nl=nl+1._r8 end if end do if (nl > 0) then gamma_sm = gamma_sm/nl endif else gamma_sm = 1.0_r8 end if gammaS_out(p)=gamma_sm ! Calculate total scaling factor !-------------------------------- gamma(p) = gamma_l * gamma_p * gamma_t * gamma_a * gamma_sm if ( (gamma(p) >=0.0_r8) .and. (gamma(p)< 100._r8) ) then gamma_out(p)=gamma(p) else gamma_out(p)=gamma(p) write(6,*) 'clh GAMMA: ',gamma(p),gamma_l,gamma_p,gamma_t,gamma_a,gamma_sm end if end do case (2,3,4,5) do fp = 1,num_soilp p = filter_soilp(fp) g = pgridcell(p) ! gamma: Activity factor. Units [dimensionless] ! ----- For monoterpenes, OVOCs, ORVOCs, CO include simple activity factors ! for LAI and T only (Guenther et al., 1995) gamma_t = exp(bet * (t_veg(p) - tstd)) gamma(p)=gamma_t end do end select do fp = 1,num_soilp p = filter_soilp(fp) g = pgridcell(p) ! density: Source density factor [g dry weight foliar mass m-2 ground] ! ------- Other than isoprene, need to convert EF units from ! [ug g-1 h-1] to [ug m-2 h-1] if (ivt(p) > noveg) then density = elai(p) / (slarea(p) * 0.5_r8) else density = 0._r8 end if ! calculate the voc flux ! ---------------------- select case (n) case(1) vocflx(p,n) = epsilon(p) * gamma(p) * scaling_to_500_Tg case(2,3,4,5) vocflx(p,n) = epsilon(p) * gamma(p) * density end select end do ! end pft loop end do ! end voc species loop !_______________________________________________________________________________ ! Calculate total voc flux and individual components for history output do fp = 1,num_soilp p = filter_soilp(fp) vocflx_tot(p) = 0._r8 end do do n = 1, nvoc do fp = 1,num_soilp p = filter_soilp(fp) vocflx_tot(p) = vocflx_tot(p) + vocflx(p,n) end do end do do fp = 1,num_soilp p = filter_soilp(fp) g = pgridcell(p) vocflx_1(p) = vocflx(p,1) vocflx_2(p) = vocflx(p,2) vocflx_3(p) = vocflx(p,3) vocflx_4(p) = vocflx(p,4) vocflx_5(p) = vocflx(p,5) end do end subroutine VOCEmission end module VOCEmissionMod