diff --git a/CHANGELOG.md b/CHANGELOG.md index 25023cc6..8ccf7b0e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,10 +9,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed +- Corrected the units of the gravimetric soil moisture to percent instead of fractional in the FENGSHA dust scheme. + ### Added +- Additional tuning parameters for the soil moisture and drylimit calculations for application specific tuning. + ### Changed +- Correct soil moisture parameterization in FENGSHA +- Add `soil_moisture_factor` to the DU2G_instance_DU.rc (same name used in the K14 scheme) and DU2G_GridCompMod.F90 files for FENGSHA +- Add `soil_drylimit_factor` to the DU2G_instance_DU.rc and DU2G_GridCompMod.F90 files for FENGSHA + ## [v2.2.1] - 2023-05-30 ### Fixed diff --git a/ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_GridCompMod.F90 b/ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_GridCompMod.F90 index 71e78c8e..d52ccdac 100644 --- a/ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_GridCompMod.F90 +++ b/ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_GridCompMod.F90 @@ -59,6 +59,7 @@ module DU2G_GridCompMod real :: alpha ! FENGSHA scaling factor real :: gamma ! FENGSHA tuning exponent real :: kvhmax ! FENGSHA max. vertical/horizontal mass flux ratio [1] + real :: f_sdl ! FENGSHA drylimit tuning factor real :: Ch_DU_res(NHRES) ! resolutions used for Ch_DU real :: Ch_DU ! dust emission tuning coefficient [kg s2 m-5]. logical :: maringFlag=.false. ! maring settling velocity correction @@ -183,6 +184,8 @@ subroutine SetServices (GC, RC) case ('fengsha') call ESMF_ConfigGetAttribute (cfg, self%alpha, label='alpha:', __RC__) call ESMF_ConfigGetAttribute (cfg, self%gamma, label='gamma:', __RC__) + call ESMF_ConfigGetAttribute (cfg, self%f_swc, label='soil_moisture_factor:', __RC__) + call ESMF_ConfigGetAttribute (cfg, self%f_sdl, label='soil_drylimit_factor:', __RC__) call ESMF_ConfigGetAttribute (cfg, self%kvhmax, label='vertical_to_horizontal_flux_ratio_limit:', __RC__) case ('k14') call ESMF_ConfigGetAttribute (cfg, self%clayFlag, label='clayFlag:', __RC__) @@ -794,10 +797,12 @@ subroutine Run1 (GC, import, export, clock, RC) if (associated(DU_EROD)) DU_EROD = f_erod_ case ('fengsha') + call DustEmissionFENGSHA (frlake, frsnow, lwi, slc, du_clay, du_sand, du_silt, & du_ssm, du_rdrag, airdens(:,:,self%km), ustar, du_uthres, & self%alpha, self%gamma, self%kvhmax, MAPL_GRAV, & - self%rhop, self%sdist, emissions_surface, __RC__) + self%rhop, self%sdist, self%f_sdl, self%f_swc, emissions_surface, __RC__) + case ('ginoux') call DustEmissionGOCART2G(self%radius*1.e-6, frlake, wet1, lwi, u10m, v10m, & diff --git a/ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_instance_DU.rc b/ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_instance_DU.rc index 30df3064..f08b9418 100644 --- a/ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_instance_DU.rc +++ b/ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_instance_DU.rc @@ -63,3 +63,5 @@ pressure_lid_in_hPa: 0.01 alpha: 0.3 gamma: 1.3 vertical_to_horizontal_flux_ratio_limit: 2.e-04 +soil_drylimit_factor: 1.0 +# soil_moisture_factor: 1.0 diff --git a/Process_Library/GOCART2G_Process.F90 b/Process_Library/GOCART2G_Process.F90 index 6892082c..cc0b5999 100644 --- a/Process_Library/GOCART2G_Process.F90 +++ b/Process_Library/GOCART2G_Process.F90 @@ -178,7 +178,7 @@ end subroutine DustAerosolDistributionKok ! !IROUTINE: soilMoistureConvertVol2Grav - volumetric to gravimetric soil moisture ! ! !INTERFACE: - real function soilMoistureConvertVol2Grav(vsoil, sandfrac, rhop) + real function soilMoistureConvertVol2Grav(vsoil, sandfrac) ! !USES: implicit NONE @@ -186,7 +186,6 @@ real function soilMoistureConvertVol2Grav(vsoil, sandfrac, rhop) ! !INPUT PARAMETERS: real, intent(in) :: vsoil ! volumetric soil moisture fraction [1] real, intent(in) :: sandfrac ! fractional sand content [1] - real, intent(in) :: rhop ! dry dust density [kg m-3] ! !DESCRIPTION: Convert soil moisture fraction from volumetric to gravimetric. ! @@ -200,16 +199,16 @@ real function soilMoistureConvertVol2Grav(vsoil, sandfrac, rhop) ! !CONSTANTS: real, parameter :: rhow = 1000. ! density of water [kg m-3] - + real, parameter :: rhop = 1700. !EOP !------------------------------------------------------------------------- ! Begin... ! Saturated volumetric water content (sand-dependent) ! [m3 m-3] - vsat = 0.489 - 0.00126 * ( 100. * sandfrac ) + vsat = 0.489 - 0.126 * sandfrac ! Gravimetric soil content - soilMoistureConvertVol2Grav = vsoil * rhow / (rhop * (1. - vsat)) + soilMoistureConvertVol2Grav = 100. * vsoil * rhow / (rhop * (1. - vsat)) end function soilMoistureConvertVol2Grav @@ -219,7 +218,7 @@ end function soilMoistureConvertVol2Grav ! !IROUTINE: moistureCorrectionFecan - Correction factor for Fecan soil moisture ! ! !INTERFACE: - real function moistureCorrectionFecan(slc, sand, clay, rhop) + real function moistureCorrectionFecan(slc, sand, clay, b) ! !USES: implicit NONE @@ -228,7 +227,7 @@ real function moistureCorrectionFecan(slc, sand, clay, rhop) real, intent(in) :: slc ! liquid water content of top soil layer, volumetric fraction [1] real, intent(in) :: sand ! fractional sand content [1] real, intent(in) :: clay ! fractional clay content [1] - real, intent(in) :: rhop ! dry dust density [kg m-3] + real, intent(in) :: b ! drylimit factor from zender 2003 ! !DESCRIPTION: Compute correction factor to account for Fecal soil moisture ! @@ -246,10 +245,10 @@ real function moistureCorrectionFecan(slc, sand, clay, rhop) ! Begin... ! Convert soil moisture from volumetric to gravimetric - grvsoilm = soilMoistureConvertVol2Grav(slc, sand, rhop) + grvsoilm = soilMoistureConvertVol2Grav(slc, sand) ! Compute fecan dry limit - drylimit = clay * (14.0 * clay + 17.0) + drylimit = b * clay * (14.0 * clay + 17.0) ! Compute soil moisture correction moistureCorrectionFecan = sqrt(1.0 + 1.21 * max(0., grvsoilm - drylimit)**0.68) @@ -303,7 +302,7 @@ end function DustFluxV2HRatioMB95 ! !INTERFACE: subroutine DustEmissionFENGSHA(fraclake, fracsnow, oro, slc, clay, sand, silt, & ssm, rdrag, airdens, ustar, uthrs, alpha, gamma, & - kvhmax, grav, rhop, distribution, emissions, rc) + kvhmax, grav, rhop, distribution, drylimit_factor, moist_correct, emissions, rc) ! !USES: implicit NONE @@ -327,7 +326,8 @@ subroutine DustEmissionFENGSHA(fraclake, fracsnow, oro, slc, clay, sand, silt, real, intent(in) :: grav ! gravity [m/sec^2] real, dimension(:), intent(in) :: rhop ! soil class density [kg/m^3] real, dimension(:), intent(in) :: distribution ! normalized dust binned distribution [1] - + real, intent(in) :: drylimit_factor ! drylimit tuning factor from zender2003 + real, intent(in) :: moist_correct ! moisture correction factor ! !OUTPUT PARAMETERS: real, intent(out) :: emissions(:,:,:) ! binned surface emissions [kg/(m^2 sec)] integer, intent(out) :: rc ! Error return code: __SUCCESS__ or __FAIL__ @@ -352,6 +352,7 @@ subroutine DustEmissionFENGSHA(fraclake, fracsnow, oro, slc, clay, sand, silt, real :: rustar real :: total_emissions real :: u_sum, u_thresh + real :: smois ! !CONSTANTS: real, parameter :: ssm_thresh = 1.e-02 ! emit above this erodibility threshold [1] @@ -405,24 +406,25 @@ subroutine DustEmissionFENGSHA(fraclake, fracsnow, oro, slc, clay, sand, silt, ! Compute threshold wind friction velocity using drag partition ! ------------------------------------------------------------- rustar = rdrag(i,j) * ustar(i,j) - + + ! Fecan moisture correction + ! ------------------------- + smois = slc(i,j) * moist_correct + h = moistureCorrectionFecan(smois, sand(i,j), clay(i,j), drylimit_factor) + + ! Adjust threshold + ! ---------------- + u_thresh = uthrs(i,j) * h + + u_sum = rustar + u_thresh + + ! Compute Horizontal Saltation Flux according to Eq (9) in Webb et al. (2020) + ! --------------------------------------------------------------------------- + q = max(0., rustar - u_thresh) * u_sum * u_sum + ! Now compute size-dependent total emission flux ! ---------------------------------------------- do n = 1, nbins - ! Fecan moisture correction - ! ------------------------- - h = moistureCorrectionFecan(slc(i,j), sand(i,j), clay(i,j), rhop(n)) - - ! Adjust threshold - ! ---------------- - u_thresh = uthrs(i,j) * h - - u_sum = rustar + u_thresh - - ! Compute Horizontal Saltation Flux according to Eq (9) in Webb et al. (2020) - ! --------------------------------------------------------------------------- - q = max(0., rustar - u_thresh) * u_sum * u_sum - ! Distribute emissions to bins and convert to mass flux (kg s-1) ! -------------------------------------------------------------- emissions(i,j,n) = distribution(n) * total_emissions * q