Merge pull request #2588 from HuiWangWanderInGitHub/CLM5_MEGAN_ISOP

New updates on MEGANv2.1

doc/ChangeLog

@@ -1,4 +1,68 @@
 ===============================================================
+Tag name: ctsm5.2.020
+Originator(s): HuiWangWanderInGitHub (Hui Wang,huiw16@uci.edu)
+Date: Mon 12 Aug 2024 10:51:04 AM MDT
+One-line Summary: MEGAN updates (MEGAN-CLM6)
+
+Purpose and description of changes
+----------------------------------
+ We add new features to MEGANv2.1 for simulating isoprene emissions based on three recent studies conducted at the BAI lab in the University of California, Irvine. The first one is about the drought impact on isoprene (Wang et al., 2022). The second study investigates the effect of temperature history on the emission factors of boreal broadleaf deciduous shrubs (Wang et al., 2024a). The third study explores a different temperature response curve for C3 Arctic grass (Wang et al., 2024b, under press). These changes improved the model's representation of isoprene emissions during drought and in high-latitude ecosystems.
+
+ This work relates to issue #1323 and is based on the old MEGANv2.1 framework, but incorporates new scientific insights.
+
+Significant changes to scientifically-supported configurations
+--------------------------------------------------------------
+
+Does this tag change answers significantly for any of the following physics configurations?
+(Details of any changes will be given in the "Answer changes" section below.)
+
+    [Put an [X] in the box for any configuration with significant answer changes.]
+
+[ ] clm6_0
+
+[ ] clm5_1
+
+[ ] clm5_0
+
+[ ] ctsm5_0-nwp
+
+[ ] clm4_5
+
+
+Notes of particular relevance for users
+---------------------------------------
+Changes to documentation:
+ Hui Wang has offered to provide the modifications to the documentation that we need.
+ The PR https://github.com/ESCOMP/CTSM/pull/2588 includes a list of references.
+
+
+Testing summary:
+----------------
+
+ [PASS means all tests PASS; OK means tests PASS other than expected fails.]
+
+  regular tests (aux_clm: https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#pre-merge-system-testing):
+
+    derecho ----- OK
+    izumi ------- OK
+
+Answer changes
+--------------
+
+Changes answers relative to baseline: Yes
+
+  Summarize any changes to answers, i.e.,
+    - what code configurations: All that calculate BVOCs
+    - what platforms/compilers: All
+    - nature of change: Larger than roundoff, only in BVOC output and specifically limited to MEG_ and _voc variables.
+
+Other details
+-------------
+Pull Requests that document the changes (include PR ids):
+ https://github.com/ESCOMP/CTSM/pull/2588 New updates on MEGANv2.1
+
+===============================================================
+===============================================================
 Tag name: ctsm5.2.019
 Originator(s): erik (Erik Kluzek,UCAR/TSS,303-497-1326)
 Date: Sun 11 Aug 2024 09:00:43 AM MDT

doc/ChangeSum

@@ -1,5 +1,6 @@
 Tag                   Who      Date  Summary
 ============================================================================================================================
+       ctsm5.2.020   slevis 08/12/2024 MEGAN updates (MEGAN-CLM6)
        ctsm5.2.019     erik 08/11/2024 Add in an additional dust emission method Leung_2023, by default off
        ctsm5.2.018 mvdebols 08/02/2024 Fix/excess ice cold start
        ctsm5.2.017     erik 07/30/2024 Dust emissions control moved to cmeps

src/biogeochem/VOCEmissionMod.F90

@@ -30,6 +30,7 @@ module VOCEmissionMod
   use SolarAbsorbedType  , only : solarabs_type
   use TemperatureType    , only : temperature_type
   use PatchType          , only : patch                
+  use EnergyFluxType     , only : energyflux_type
   !
   implicit none
   private 
@@ -378,7 +379,7 @@ end subroutine InitCold
   !-----------------------------------------------------------------------
   subroutine VOCEmission (bounds, num_soilp, filter_soilp, &
        atm2lnd_inst, canopystate_inst, photosyns_inst, temperature_inst, &
-       vocemis_inst)
+       vocemis_inst, energyflux_inst)
     !
     ! ! NEW DESCRIPTION
     ! Volatile organic compound emission
@@ -422,6 +423,7 @@ subroutine VOCEmission (bounds, num_soilp, filter_soilp, &
     type(photosyns_type)   , intent(in)    :: photosyns_inst
     type(temperature_type) , intent(in)    :: temperature_inst
     type(vocemis_type)     , intent(inout) :: vocemis_inst
+    type(energyflux_type)  , intent(in)    :: energyflux_inst
     !
     ! !REVISION HISTORY:
     ! 4/29/11: Colette L. Heald: expand MEGAN to 20 compound classes
@@ -476,14 +478,7 @@ subroutine VOCEmission (bounds, num_soilp, filter_soilp, &
     end if
 
     associate(                                                    & 
-         !dz           => col%dz                                , & ! Input:  [real(r8) (:,:) ]  depth of layer (m)                              
-         !bsw          => soilstate_inst%bsw_col                , & ! Input:  [real(r8) (:,:) ]  Clapp and Hornberger "b" (nlevgrnd)             
-         !clayfrac     => soilstate_inst%clayfrac_col           , & ! Input:  [real(r8) (:)   ]  fraction of soil that is clay                     
-         !sandfrac     => soilstate_inst%sandfrac_col           , & ! Input:  [real(r8) (:)   ]  fraction of soil that is sand                     
-         !watsat       => soilstate_inst%watsat_col             , & ! Input:  [real(r8) (:,:) ]  volumetric soil water at saturation (porosity) (nlevgrnd)
-         !sucsat       => soilstate_inst%sucsat_col             , & ! Input:  [real(r8) (:,:) ]  minimum soil suction (mm) (nlevgrnd)            
-         !h2osoi_vol   => waterstate_inst%h2osoi_vol_col        , & ! Input:  [real(r8) (:,:) ]  volumetric soil water (m3/m3)                   
-         !h2osoi_ice   => waterstate_inst%h2osoi_ice_col        , & ! Input:  [real(r8) (:,:) ]  ice soil content (kg/m3)                        
+         btran         => energyflux_inst%btran_patch           , & ! Input:  [real(r8) (:)   ]  transpiration wetness factor (0 to 1)
 
          forc_solad    => atm2lnd_inst%forc_solad_downscaled_col, & ! Input:  [real(r8) (:,:) ]  direct beam radiation (visible only)            
          forc_solai    => atm2lnd_inst%forc_solai_grc           , & ! Input:  [real(r8) (:,:) ]  diffuse radiation     (visible only)            
@@ -569,10 +564,8 @@ subroutine VOCEmission (bounds, num_soilp, filter_soilp, &
           ! Activity factor for LAI (Guenther et al., 2006): all species
           gamma_l = get_gamma_L(fsun240(p), elai(p))
 
-          ! Activity factor for soil moisture: all species (commented out for now)
-          !          gamma_sm = get_gamma_SM(clayfrac(p), sandfrac(p), h2osoi_vol(c,:), h2osoi_ice(c,:), &
-          !               col%dz(c,:), soilstate_inst%bsw_col(c,:), watsat(c,:), sucsat(c,:), root_depth(patch%itype(p)))
-          gamma_sm = 1.0_r8
+          ! Impact of soil moisture on isoprene emission
+          gamma_sm = get_gamma_SM(btran(p))
 
           ! Loop through VOCs for light, temperature and leaf age activity factor & apply
           ! all final activity factors to baseline emission factors
@@ -599,7 +592,7 @@ subroutine VOCEmission (bounds, num_soilp, filter_soilp, &
 
              ! Activity factor for T
              gamma_t = get_gamma_T(t_veg240(p), t_veg24(p),t_veg(p), ct1(class_num), ct2(class_num),&
-                                   betaT(class_num),LDF(class_num), Ceo(class_num), Eopt, topt)
+                                   betaT(class_num),LDF(class_num), Ceo(class_num), Eopt, topt, patch%itype(p))
 
              ! Activity factor for Leaf Age
              gamma_a = get_gamma_A(patch%itype(p), elai240(p),elai(p),class_num)
@@ -818,90 +811,53 @@ function get_gamma_L(fsun240_in,elai_in)
   end function get_gamma_L
 
   !-----------------------------------------------------------------------
-  function get_gamma_SM(clayfrac_in, sandfrac_in, h2osoi_vol_in, h2osoi_ice_in, dz_in, &
-       bsw_in, watsat_in, sucsat_in, root_depth_in)
-    !
-    ! Activity factor for soil moisture (Guenther et al., 2006): all species
-    !----------------------------------
-    ! 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
-    !
-    ! !USES:
-    use clm_varcon , only : denice
-    use clm_varpar , only : nlevsoi
-    !
-    ! !ARGUMENTS:
-    implicit none
-    real(r8),intent(in) :: clayfrac_in
-    real(r8),intent(in) :: sandfrac_in
-    real(r8),intent(in) :: h2osoi_vol_in(nlevsoi)
-    real(r8),intent(in) :: h2osoi_ice_in(nlevsoi)
-    real(r8),intent(in) :: dz_in(nlevsoi)
-    real(r8),intent(in) :: bsw_in(nlevsoi)
-    real(r8),intent(in) :: watsat_in(nlevsoi)
-    real(r8),intent(in) :: sucsat_in(nlevsoi)
-    real(r8),intent(in) :: root_depth_in
-    !
-    ! !LOCAL VARIABLES:
-    real(r8)            :: get_gamma_SM
-    integer  :: j
-    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
-    real(r8), parameter :: deltheta1=0.06_r8               ! empirical coefficient
-    real(r8), parameter :: smpmax = 2.57e5_r8              ! maximum soil matrix potential
-    !-----------------------------------------------------------------------
+  function get_gamma_SM(btran_in)
 
-    if ((clayfrac_in > 0) .and. (sandfrac_in > 0)) then 
-       get_gamma_SM = 0._r8
-       nl=0._r8
-
-       do j = 1,nlevsoi
-          if  (sum(dz_in(1:j)) < root_depth_in)  then
-             theta_ice = h2osoi_ice_in(j)/(dz_in(j)*denice)
-             wilt = ((smpmax/sucsat_in(j))**(-1._r8/bsw_in(j))) * (watsat_in(j) - theta_ice)
-             theta1 = wilt + deltheta1
-             if (h2osoi_vol_in(j) >= theta1) then 
-                get_gamma_SM = get_gamma_SM + 1._r8
-             else if ( (h2osoi_vol_in(j) > wilt) .and. (h2osoi_vol_in(j) < theta1) ) then
-                get_gamma_SM = get_gamma_SM + ( h2osoi_vol_in(j) - wilt ) / deltheta1
-             else
-                get_gamma_SM = get_gamma_SM + 0._r8
-             end if
-             nl=nl+1._r8
-          end if
-       end do
-
-       if (nl > 0._r8) then
-          get_gamma_SM = get_gamma_SM/nl
-       endif
+    !---------------------------------------
+    ! May 22, 2024
+    ! Activity factor for soil moisture of Isoprene (Wang et al., 2022, JAMES)
+    ! It is based on eq. (11) in the paper. Because the temperature response
+    ! of isoprene has been explicitly included in CLM;
 
-       if (get_gamma_SM > 1.0_r8) then 
-          write(iulog,*) 'healdSM > 1: gamma_SM, nl', get_gamma_SM, nl
-          get_gamma_SM=1.0_r8
-       endif
+    !ARGUMENTS:
+    implicit none
+    real(r8),intent(in) :: btran_in
+
+    !!!------- the drought algorithm--------
+    real(r8), parameter :: a1 = -7.4463_r8
+    real(r8), parameter :: b1 = 3.2552_r8
+    real(r8), parameter :: btran_threshold = 0.2_r8
+    real(r8)            :: get_gamma_SM
+    !---------------------------------------
 
+    if (btran_in >= 1._r8) then
+       get_gamma_SM = 1._r8
     else
-       get_gamma_SM = 1.0_r8
-    end if
+       get_gamma_SM = 1._r8 / (1._r8 + b1 * exp(a1 * (btran_in - btran_threshold)))
+    endif
 
   end function get_gamma_SM
-  
+
   !-----------------------------------------------------------------------
-  function get_gamma_T(t_veg240_in, t_veg24_in,t_veg_in, ct1_in, ct2_in, betaT_in, LDF_in, Ceo_in, Eopt, topt)
+  function get_gamma_T(t_veg240_in, t_veg24_in,t_veg_in, ct1_in, ct2_in, betaT_in, LDF_in, Ceo_in, Eopt, topt, ivt_in)
 
-    ! Activity factor for temperature 
+    ! Activity factor for temperature
+    !--------------------------------
+    ! May 24, 2024 Hui updated the temperature response curves of isoprene for
+    ! Boreal Broadleaf Deciduous Shrub and Arctic C3 grass based on
+    ! Wang et al., 2024 (GRL) and Wang et al., 2024 (Nature Communications)
     !--------------------------------
     ! Calculate both a light-dependent fraction as in Guenther et al., 2006 for isoprene
     ! of a max saturation type form. Also caculate a light-independent fraction of the
     ! form of an exponential. Final activity factor depends on light dependent fraction
     ! of compound type.
     !
+    ! !USES:
+    use clm_varcon, only: tfrz
+
     ! !ARGUMENTS:
     implicit none
+    integer,intent(in)  :: ivt_in
     real(r8),intent(in) :: t_veg240_in
     real(r8),intent(in) :: t_veg24_in
     real(r8),intent(in) :: t_veg_in
@@ -917,7 +873,9 @@ function get_gamma_T(t_veg240_in, t_veg24_in,t_veg_in, ct1_in, ct2_in, betaT_in,
     real(r8) :: get_gamma_T
     real(r8) :: gamma_t_LDF             ! activity factor for temperature
     real(r8) :: gamma_t_LIF             ! activity factor for temperature
-    real(r8) :: x                       ! temporary 
+    real(r8) :: x                       ! temporary i
+    real(r8) :: bet_arc_c3  ! activity factor for temperature for arctic C3 grass
+    real(r8), parameter :: bet_arc_c3_max = 300._r8  ! max value, activity factor for temperature for arctic C3 graass
     real(r8), parameter :: co1 = 313._r8                   ! empirical coefficient
     real(r8), parameter :: co2 = 0.6_r8                    ! empirical coefficient
     real(r8), parameter :: co4 = 0.05_r8                   ! empirical coefficient
@@ -927,19 +885,40 @@ function get_gamma_T(t_veg240_in, t_veg24_in,t_veg_in, ct1_in, ct2_in, betaT_in,
     real(r8), parameter :: ct3 = 0.00831_r8                ! empirical coefficient (0.0083 in User's Guide)
     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 :: std_act_energy_isopr = 95._r8  ! standard activation energy for isoprene
+    real(r8), parameter :: empirical_param_1 = 9.49_r8  ! empirical param for the activation energy in response to 10-day temperature change
+    real(r8), parameter :: empirical_param_2 = 0.53_r8  ! empirical param for the activation energy in response to 10-day temperature change
+    real(r8), parameter :: empirical_param_3 = 0.12_r8  ! empirical param for the emission factors of arctic C3 grass in response to 10-day temperature change
+    real(r8), parameter :: empirical_param_4 = 7.9_r8  ! empirical param for the emission factors of broadleaf deciduous boreal shrubs in response to 10-day temperature change
+    real(r8), parameter :: empirical_param_5 = 0.217_r8  ! empirical param for the emission factors of broadleaf deciduous boreal shrubs in response to 10-day temperature change
     !-----------------------------------------------------------------------
 
     ! Light dependent fraction (Guenther et al., 2006)
     if ( (t_veg240_in > 0.0_r8) .and. (t_veg240_in < 1.e30_r8) ) then 
        ! topt and Eopt from eq 8 and 9:
        topt = co1 + (co2 * (t_veg240_in-tstd0))
-       Eopt = Ceo_in * exp (co4 * (t_veg24_in-tstd0)) * exp(co4 * (t_veg240_in -tstd0))
+       if ( (ivt_in == nbrdlf_dcd_brl_shrub) ) then  ! boreal-deciduous-shrub
+       ! coming from BEAR-oNS campaign willows results
+          Eopt = empirical_param_4 * exp (empirical_param_5 * (t_veg24_in - tfrz - 24.0_r8))
+       else if ( (ivt_in == nc3_arctic_grass ) ) then  ! boreal-grass
+          Eopt = exp(empirical_param_3 * (t_veg240_in - tfrz - 15._r8))
+       else
+          Eopt = Ceo_in * exp (co4 * (t_veg24_in - tstd0)) * exp(co4 * (t_veg240_in - tstd0))
+       endif
+
     else
        topt = topt_fix
        Eopt = Eopt_fix
     endif
     x = ( (1._r8/topt) - (1._r8/(t_veg_in)) ) / ct3
-    gamma_t_LDF = Eopt * ( ct2_in * exp(ct1_in * x)/(ct2_in - ct1_in * (1._r8 - exp(ct2_in * x))) )
+    ! for the boreal grass from BEAR-oNS campaign
+    if ( (ivt_in == nc3_arctic_grass ) ) then  ! boreal-grass
+        bet_arc_c3 = std_act_energy_isopr + empirical_param_1 * exp(empirical_param_2 * (tfrz + 15._r8 - t_veg240_in))
+        bet_arc_c3 = min(bet_arc_c3, bet_arc_c3_max)
+        gamma_t_LDF = Eopt * exp(bet_arc_c3 * ((1._r8 / (tfrz + 30._r8) - 1._r8 / t_veg_in) / ct3))
+    else
+        gamma_t_LDF = Eopt * ( ct2_in * exp(ct1_in * x) / (ct2_in - ct1_in * (1._r8 - exp(ct2_in * x))) )
+    endif
 
 
     ! Light independent fraction (of exp(beta T) form)

src/main/clm_driver.F90

@@ -815,11 +815,11 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro
        call dust_emis_inst%DustDryDep(bounds_clump, &
             atm2lnd_inst, frictionvel_inst)
 
-       ! VOC emission (A. Guenther's MEGAN (2006) model)
+       ! VOC emission (A. Guenther's MEGAN (2006) model; Wang et al., 2022, 2024a, 2024b)
        call VOCEmission(bounds_clump,                                         &
                filter(nc)%num_soilp, filter(nc)%soilp,                           &
                atm2lnd_inst, canopystate_inst, photosyns_inst, temperature_inst, &
-               vocemis_inst)
+               vocemis_inst, energyflux_inst)
 
        call t_stopf('bgc')