From b3b8e08acdcad393844207ac9e556bdee8b3ef08 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Mon, 3 Jun 2024 15:42:23 -0600 Subject: [PATCH] Add reference for ThresholdSoilMoist and add extended comments from Danny Leung from the Leung PR, actually use it in the code and make sure the unit test uses the new name --- .../DustEmis_test/test_DustEmisZender2003.pf | 4 +- src/biogeophys/SoilStateInitTimeConstMod.F90 | 49 ++++++++++++++++--- 2 files changed, 45 insertions(+), 8 deletions(-) diff --git a/src/biogeochem/test/DustEmis_test/test_DustEmisZender2003.pf b/src/biogeochem/test/DustEmis_test/test_DustEmisZender2003.pf index c99f4bd97d..a7fea904c9 100644 --- a/src/biogeochem/test/DustEmis_test/test_DustEmisZender2003.pf +++ b/src/biogeochem/test/DustEmis_test/test_DustEmisZender2003.pf @@ -19,7 +19,7 @@ module test_DustEmisZender2003 use WaterType, only : water_type use FrictionVelocityMod, only : frictionvel_type use unittestWaterTypeFactory, only : unittest_water_type_factory_type - use SoilStateInitTimeConstMod, only : ThresholdSoilMoist, ThresholdSoilMoistKok2014, MassFracClay + use SoilStateInitTimeConstMod, only : ThresholdSoilMoistZender2003, ThresholdSoilMoistKok2014, MassFracClay implicit none @@ -156,7 +156,7 @@ contains ! considered bedrock col%nbedrock(c) = nlevsoi - this%soilstate_inst%gwc_thr_col(c) = ThresholdSoilMoist( clay ) + this%soilstate_inst%gwc_thr_col(c) = ThresholdSoilMoistZender2003( clay ) this%soilstate_inst%mss_frc_cly_vld_col(c) = MassFracClay( clay ) end do diff --git a/src/biogeophys/SoilStateInitTimeConstMod.F90 b/src/biogeophys/SoilStateInitTimeConstMod.F90 index 3b71efe103..0be54571e5 100644 --- a/src/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/biogeophys/SoilStateInitTimeConstMod.F90 @@ -19,7 +19,7 @@ module SoilStateInitTimeConstMod public :: readParams ! PRIVATE FUNCTIONS MADE PUBLIC Just for unit-testing: - public :: ThresholdSoilMoist + public :: ThresholdSoilMoistZender2003 public :: ThresholdSoilMoistKok2014 public :: MassFracClay ! @@ -707,8 +707,8 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) do c = begc,endc g = col%gridcell(c) - soilstate_inst%gwc_thr_col(c) = 0.17_r8 + 0.14_r8 * clay3d(g,1) * 0.01_r8 - soilstate_inst%mss_frc_cly_vld_col(c) = min(clay3d(g,1) * 0.01_r8, 0.20_r8) + soilstate_inst%gwc_thr_col(c) = ThresholdSoilMoistZender2003( clay3d(g,1) ) + soilstate_inst%mss_frc_cly_vld_col(c) = MassFracClay( clay3d(g,1) ) end do ! -------------------------------------------------------------------- @@ -722,9 +722,31 @@ end subroutine SoilStateInitTimeConst !------------------------------------------------------------------------------ - real(r8) function ThresholdSoilMoist( clay ) + real(r8) function ThresholdSoilMoistZender2003( clay ) + !------------------------------------------------------------------------------ + ! ! Calculate the threshold soil moisture needed for dust emission, based on clay content - use abortUtils , only : endrun + ! This was the original equation with a = 1 / (%clay) being the tuning factor for soil + ! moisture effect in Zender's 2003 dust emission scheme. + ! + ! 0.17 and 0.14 are fitting coefficients in Fecan et al. (1999), and 0.01 is used to + ! convert surface clay fraction from percentage to fraction. + ! + !------------------------------------------------------------------------------ + ! For future developments Danny M. Leung decided (Dec, 2023) that the Leung et al. (2023) o + ! dust emission scheme in the CESM will use Zender's tuning as well, which overall + ! encourages more dust emissions from seriamid and more marginal dust sources. + ! Another advantage of using this tuning factor instead of a = 1 is that the dust emission + ! threshold is linearly dependent on the clay fraction instead of parabolically dependent + ! on clay fraction as in the above line. This means that dust emission becomes a little + ! less sensitive to clay content (soil texture). + ! + ! Also see the notes below for ThresholdSoilMoistKok2014. + ! + ! Notes from: dmleung 19 Feb 2024. + ! + !------------------------------------------------------------------------------ + use abortUtils , only : endrun use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) real(r8), intent(IN) :: clay ! Fraction of clay in the soil (%) @@ -734,12 +756,27 @@ real(r8) function ThresholdSoilMoist( clay ) return end if ThresholdSoilMoist = 0.17_r8 + 0.14_r8 * clay * 0.01_r8 - end function ThresholdSoilMoist + end function ThresholdSoilMoistZender2003 !------------------------------------------------------------------------------ real(r8) function ThresholdSoilMoistKok2014( clay ) + !------------------------------------------------------------------------------ ! Calculate the threshold soil moisture needed for dust emission, based on clay content + ! + ! The below calculates the threshold gravimetric water content for the dust emission + ! calculation in DustEmis. The equation comes from Eq. 14 of Fecan et al. + ! (1999; https://doi.org/10.1007/s00585-999-0149-7). + ! gwc_thr_col = 0.17*clay3d + 0.0014*(clay3d**2), and we only concern the topmost + ! soil layer. Charlie Zender later on added a tuning factor (a) such that the + ! equation becomes gwc_thr_col = a*[0.17*clay3d + 0.0014*(clay3d**2)]. + ! (Zender et al., 2003a; https://doi.org/10.1029/2002JD002775) + ! Kok et al. (2014a, b) chose to use a = 1. Resulting in this function + ! Charlie Zender (2003a) chose: a = 1/clay3d, which gives the ThresholdSoilMoistZender2003 + ! function above. + ! + ! Notes from: dmleung 24 May 2024. + !------------------------------------------------------------------------------ real(r8), intent(IN) :: clay ! Fraction of clay in the soil (%) ThresholdSoilMoistKok2014 = 0.01_r8*(0.17_r8*clay + 0.0014_r8*clay*clay)