Merge pull request ESCOMP#2803 from dmleung/dustemis_tune4cesm3

Tuning Leung_2023 dust emissions for clm6_0_cam7.0

bld/unit_testers/build-namelist_test.pl

@@ -1359,8 +1359,8 @@ sub cat_and_create_namelistinfile {
                                      namelst=>"fsurdat='build-namelist_test.pl'",
                                      phys=>"clm6_0",
                                    },
-     "hillslope with init_interp"=>{ options=>"-bgc bgc -envxml_dir .",
-                                     namelst=>"use_init_interp=.true.,use_hillslope=.true.",
+     "hillslope with init_interp"=>{ options=>"--res 10x15 --bgc bgc --envxml_dir .",
+                                     namelst=>"use_init_interp=.true.,use_hillslope=.true.,hillslope_file='/dev/null'",
                                      phys=>"clm6_0",
                                    },
                );

cime_config/testdefs/ExpectedTestFails.xml

@@ -135,6 +135,17 @@
     </phase>
   </test>
 
+  <test name="SMS_Ld10_D_Mmpi-serial.CLM_USRDAT.I1PtClm60Bgc.izumi_nag.clm-NEON-MOAB--clm-PRISM">
+    <phase name="SHAREDLIB_BUILD">
+      <status>FAIL</status>
+      <issue>#2310</issue>
+    </phase>
+    <phase name="RUN">
+      <status>FAIL</status>
+      <issue>#2310</issue>
+    </phase>
+  </test>
+
   <test name="SMS_Ld10_D_Mmpi-serial.CLM_USRDAT.I1PtClm60Fates.derecho_gnu.clm-FatesPRISM--clm-NEON-FATES-YELL">
     <phase name="SHAREDLIB_BUILD">
       <status>FAIL</status>
@@ -171,6 +182,13 @@
     </phase>
   </test>
 
+  <test name="SMS_Lm3_D_Mmpi-serial.1x1_brazil.I2000Clm50FatesCruRsGs.izumi_intel.clm-FatesColdHydro">
+    <phase name="CREATE_NEWCASE">
+      <status>FAIL</status>
+      <issue>#2454</issue>
+    </phase>
+  </test>
+
   <test name="SMS_Ld10_D_Mmpi-serial.CLM_USRDAT.I1PtClm60Fates.izumi_nag.clm-FatesFireLightningPopDens--clm-NEON-FATES-NIWO">
     <phase name="SHAREDLIB_BUILD">
       <status>FAIL</status>

doc/ChangeLog

@@ -1,4 +1,95 @@
 ===============================================================
+Tag name: ctsm5.3.005
+Originator(s): dmleung (Danny Leung)
+Date: Thu 10 Oct 2024 03:15:52 AM MDT
+One-line Summary: Hardcoded tuning adjustments for Leung_2024 dust emissions
+
+Purpose and description of changes
+----------------------------------
+
+Changes to dust emissions when Leung_2024 method is being used.
+
+Tuning was needed since we saw some high biases in dust over semiarid regions given all the updates in CTSM and CAM. The biggest
+changes in CTSM that affects dust is an increase in friction velocity (ustar; fv in CTSM) over vegetated, semiarid regions, mainly
+due to a switch in the roughness scheme from 'ZengWang2007' to 'Meier2022'. Since dust emission is very sensitive to ustar, the dust
+emission scheme magnifies this increase and caused strong high biases over semiarid regions, including Australia and Patagonia (see
+plots in issue #2732). To enhance the robustness of Leung_2023 and reduce the likelihood to see huge changes in answers in the
+future, we tried to limit the sensitivity of dust emissions to ustar. We also tried different methods (see specific notes below) to
+reduce dust emissions from semiarid regions given the high biases over there.
+
+
+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.]
+
+[x] clm6_0
+
+[ ] clm5_1
+
+[ ] clm5_0
+
+[ ] ctsm5_0-nwp
+
+[ ] clm4_5
+
+
+Bugs fixed
+----------
+List of CTSM issues fixed (include CTSM Issue # and description) [one per line]:
+
+   Start of work on tuning for #2732
+   Fixes #2826 -- Fix failing warning test
+
+Notes of particular relevance for developers:
+---------------------------------------------
+
+Caveats for developers (e.g., code that is duplicated that requires double maintenance):
+   Dust tuning parameters are hardcoded in FORTRAN and need to be moved to namelists
+
+Testing summary: Regular
+----------------
+ [PASS means all tests PASS; OK means tests PASS other than expected fails.]
+
+  build-namelist tests (if CLMBuildNamelist.pm has changed):
+
+    derecho - PASS
+
+  python testing (if python code has changed; see instructions in python/README.md; document testing done):
+
+    derecho - PASS
+
+  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 for clm6_0
+
+  Summarize any changes to answers, i.e.,
+    - what code configurations: clm6_0
+    - what platforms/compilers: all
+    - nature of change (roundoff; larger than roundoff/same climate; new climate):
+        Just updates dust emissions when Leung_2024 method is used
+
+    Only the dust emission fields are changed all other fields are untouched
+
+Other details
+-------------
+
+Pull Requests that document the changes (include PR ids):
+(https://github.com/ESCOMP/ctsm/pull)
+
+   https://github.com/ESCOMP/CTSM/pull/2803 -- Tuning Leung_2023 dust emissions for clm6_0_cam7.0
+
+===============================================================
+===============================================================
 Tag name: ctsm5.3.004
 Originator(s): samrabin (Sam Rabin)
 Date: Mon Oct  7 21:25:20 MDT 2024

doc/ChangeSum

@@ -1,5 +1,6 @@
 Tag                   Who      Date  Summary
 ============================================================================================================================
+       ctsm5.3.005     erik 10/10/2024 Hardcoded tuning adjustments for Leung_2024 dust emissions
        ctsm5.3.004 samrabin 10/07/2024 Move hillslope data off surface datasets
        ctsm5.3.003 multiple 10/07/2024 FATES default parameter file update
        ctsm5.3.002 multiple 09/26/2024 Duplicate tag (same as 5.3.001)

src/biogeochem/DustEmisLeung2023.F90

@@ -358,10 +358,6 @@ subroutine DustEmission (this, bounds, &
     ! !LOCAL VARIABLES
     integer  :: fp,p,c,l,g,m,n      ! indices
     real(r8) :: liqfrac             ! fraction of total water that is liquid
-    real(r8) :: wnd_frc_rat         ! [frc] Wind friction threshold over wind friction
-    real(r8) :: wnd_frc_slt_dlt     ! [m s-1] Friction velocity increase from saltatn
-    real(r8) :: wnd_rfr_dlt         ! [m s-1] Reference windspeed excess over threshld
-    real(r8) :: dst_slt_flx_rat_ttl
     real(r8) :: flx_mss_hrz_slt_ttl
     real(r8) :: flx_mss_vrt_dst_ttl(bounds%begp:bounds%endp)
     real(r8) :: frc_thr_wet_fct
@@ -384,14 +380,14 @@ subroutine DustEmission (this, bounds, &
     real(r8) :: wnd_frc_thr_slt_it  ! [m/s] created for impact threshold friction velocity
     real(r8) :: wnd_frc_thr_slt     ! [m/s] used for wet fluid threshold friction velocity
     real(r8) :: K_length            ! [dimless] normalized mean interobstacle distance, or called gap length (Okin, 2008)
+    ! dmleung has these variables and will change them into pointers and prepare for their history outputs. 30 Sep 2024
     real(r8) :: bare_frc            ! LUH2 bare soil land cover fraction
-    real(r8) :: veg_frc             ! LUH2 natural vegetation + crop land cover fraction
+    real(r8) :: natveg_frc          ! LUH2 natural vegetation cover fraction
+    real(r8) :: crop_frc            ! LUH2 crop cover fraction.
     !
     ! constants
     !
-    real(r8), parameter :: cst_slt = 2.61_r8           ! [frc] Saltation constant
-    real(r8), parameter :: flx_mss_fdg_fct = 5.0e-4_r8 ! [frc] Empir. mass flx tuning eflx_lh_vegt
-    real(r8), parameter :: vai_mbl_thr = 1.0_r8        ! [m2 m-2] new VAI threshold; Danny M. Leung suggests 1, and Zender's scheme uses 0.3
+    real(r8), parameter :: vai_mbl_thr = 0.6_r8        ! [m2 m-2] new VAI threshold; Danny M. Leung suggests something between 0.6 and 1 for tuning. Zender's scheme uses 0.3. Simone Tilmes might want this as a namelist variable for easier CESM tuning. dmleung 30 Sep 2024.
 
     real(r8), parameter :: Cd0 = 4.4e-5_r8             ! [dimless] proportionality constant in calculation of dust emission coefficient
     real(r8), parameter :: Ca = 2.7_r8                 ! [dimless] proportionality constant in scaling of dust emission exponent
@@ -407,11 +403,12 @@ subroutine DustEmission (this, bounds, &
     real(r8), parameter :: D_p = 130e-6_r8             ! [m] Medium soil particle diameter, assuming a global constant of ~130 um following Leung et al. (2023). dmleung 16 Feb 2024
     real(r8), parameter :: gamma_Shao = 1.65e-4_r8     ! [kg s-2] interparticle cohesion: fitting parameter in Shao and Lu (2000) (S&L00). dmleung 16 Feb 2024
     real(r8), parameter :: A_Shao = 0.0123_r8          ! [dimless] coefficient for aerodynamic force: fitting parameter in Shao and Lu (2000). dmleung 16 Feb 2024
-    real(r8), parameter :: frag_expt_thr = 5.0_r8      ! [dimless] threshold for fragmentation exponent defined in Leung et al. (2023), somewhere within 3 to 5. It is used to prevent a local AOD blowup (over Patagonia, Argentina), but one can test larger values and relax the threshold if wanted. dmleung 16 Feb 2024
+    real(r8), parameter :: frag_expt_thr = 2.5_r8      ! [dimless] Maximum value or threshold for fragmentation exponent defined in Leung et al. (2023). Danny M. Leung suggested it to be somewhere between 3 and 5 for tuning. It is used to prevent a local AOD blowup (over Patagonia, Argentina), but one can test larger values and relax the threshold if wanted. dmleung 16 Feb 2024. Update: Simone Tilmes might want this as a namelist variable for easier CESM tuning. 30 Sep 2024.
     real(r8), parameter :: z0a_glob = 1e-4_r8          ! [m] assumed globally constant aeolian roughness length value in Leung et al. (2023), for the log law of the wall for Comola et al. (2019) intermittency scheme. dmleung 20 Feb 2024
     real(r8), parameter :: hgt_sal = 0.1_r8            ! [m] saltation height used by Comola et al. (2019) intermittency scheme for the log law of the wall. dmleung 20 Feb 2024
     real(r8), parameter :: vai0_Okin = 0.1_r8          ! [m2/m2] minimum VAI needed for Okin-Pierre's vegetation drag partition equation. lai=0 in the equation will lead to infinity, so a small value is added into this lai dmleung defined.
     real(r8), parameter :: zii = 1000.0_r8             ! [m] convective boundary layer height added by dmleung 20 Feb 2024, following other CTSM modules (e.g., CanopyFluxesMod). Should we transfer PBL height (PBLH) from CAM?
+    real(r8), parameter :: dust_veg_drag_fact = 0.7_r8 ! [dimless] dmleung added a tuning factor for Greg Okin's vegetation drag partition effect. dmleung suggested a smaller vegetation drag partition effect given an increase in vegetation roughness after CTSM switched from using ZengWang2007 to Meier2022. This is simply because the drag partition effect should decrease with increasing roughness, but Okin's scheme is only a function of LAI. One might want to change this factor to 1_r8 when using ZengWang2007. dmleung 30 Sep 2024
     real(r8) :: numer                                  ! Numerator term for threshold crossing rate
     real(r8) :: denom                                  ! Denominator term for threshold crossing rate
     !------------------------------------------------------------------------
@@ -650,8 +647,8 @@ subroutine DustEmission (this, bounds, &
 
 
             ! calculate Okin's shear stress ratio (SSR, which is vegetation drag partition factor) using Pierre's equation
-            K_length = 2.0_r8 * (1.0_r8/vai_Okin(p) - 1.0_r8)   ! Here LAI has to be non-zero to avoid blowup, and < 1 to avoid -ve K_length. See this equation in Leung et al. (2023). This line is Okin's formulation
-            ssr(p) = (K_length+f_0*c_e)/(K_length+c_e) ! see this equation in Caroline Pierre et al. (2014) or Leung et al. (2023). This line is Pierre's formulation.
+            K_length = 2.0_r8 * (1.0_r8/vai_Okin(p) - 1.0_r8)   ! Here VAI has to be non-zero to avoid blowup, and < 1 to avoid -ve K_length. See this equation in Leung et al. (2023). This line is Okin's formulation
+            ssr(p) = dust_veg_drag_fact * (K_length+f_0*c_e)/(K_length+c_e) ! see this equation in Caroline Pierre et al. (2014) or Leung et al. (2023). This line is Pierre's formulation. dmleung added a tuning factor for Okin's vegetation drag partition effect (SSR) on 30 Sep 2024.
 
             ! calculation of the hybrid/total drag partition effect considering both rock and vegetation drag partitioning using LUH2 bare and veg fractions within a grid
             if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then
@@ -886,4 +883,4 @@ end subroutine SetDragPartition
 
  !==============================================================================
 
-end module DustEmisLeung2023
+end module DustEmisLeung2023

src/biogeophys/SoilStateInitTimeConstMod.F90

@@ -220,6 +220,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
     real(r8)           :: perturbed_sand                ! temporary for paramfile implementation of +/- sand percentage
     real(r8)           :: residual_clay_frac            ! temporary for paramfile implementation of +/- residual clay percentage
     real(r8)           :: perturbed_residual_clay_frac  ! temporary for paramfile implementation of +/- residual clay percentage
+    real(r8)           :: dust_moist_fact               ! tuning factor for soil moisture effect on limiting dust emissions, used by Charlie Zender. Simone Tilmes suggested to change this parameter into a namelist variable for easier CESM tuning. dmleung added 30 Sep 2024
     integer            :: dimid                         ! dimension id
     logical            :: readvar 
     type(file_desc_t)  :: ncid                          ! netcdf id
@@ -708,17 +709,24 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
     ! Initialize threshold soil moisture, and mass fraction of clay as
     ! scaling coefficient of dust emission flux (kg/m2/s) in each DustEmisType
     ! module. See the comments in each function.
+    ! Zender suggested that threshold soil moisture is tunable (see comment
+    ! inside ThresholdSoilMoistZender2003). dmleung further add dust_moist_fact
+    ! ofr modelers to tune the threshold soil moisture. The resulting tuning 
+    ! factor is thus a = dust_moist_fact / (clay3d). dmleung 30 Sep 2024
     ! --------------------------------------------------------------------
 
     do c = begc,endc
        g = col%gridcell(c)
 
-       soilstate_inst%gwc_thr_col(c) = ThresholdSoilMoistZender2003( clay3d(g,1) )
+       !soilstate_inst%gwc_thr_col(c) = ThresholdSoilMoistZender2003( clay3d(g,1) )
        if ( is_dust_emis_leung() )then
           soilstate_inst%mss_frc_cly_vld_col(c) = MassFracClayLeung2023( clay3d(g,1) )
+          dust_moist_fact = 1.0_r8   ! change this into a namelist variable later., currrently not used but could be in the future
        else
           soilstate_inst%mss_frc_cly_vld_col(c) = MassFracClay( clay3d(g,1) )
+          dust_moist_fact = 1.0_r8
        end if
+       soilstate_inst%gwc_thr_col(c) = dust_moist_fact * ThresholdSoilMoistZender2003( clay3d(g,1) )
 
     end do
 
@@ -739,9 +747,11 @@ real(r8) function ThresholdSoilMoistZender2003( clay )
   ! Calculate the threshold gravimetric water content needed for dust emission, based on clay content
   ! This was the original equation with a = 1 / (%clay) being the tuning factor for soil
   ! moisture effect in Zender's 2003 dust emission scheme (only for top layer).
+  ! dmleung further added dust_moist_fact for more flexibility in tuning, so the tuning factor here
+  ! is a = dust_moist_fact / (%clay). dmleung added dust_moist_fact on 30 Sep 2024.
   !
   ! 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.
+  ! convert surface clay from percentage to fraction.
   ! The equation comes from Eq. 14 of Fecan et al. (1999; https://doi.org/10.1007/s00585-999-0149-7).
   !
   ! NOTE: dmleung 19 Feb 2024.