add entrainment and detrainment limiter to diagnostic edmf

config/model_configs/diagnostic_edmfx_aquaplanet.yml

@@ -12,7 +12,7 @@ edmfx_sgs_diffusive_flux: true
 moist: equil 
 precip_model: 0M
 override_τ_precip: false
-dt: 10secs 
+dt: 50secs 
 t_end: 1hours 
 dt_save_to_disk: 600secs
 toml: [toml/diagnostic_edmfx_trmm_box.toml]

config/model_configs/diagnostic_edmfx_bomex_box.yml

@@ -25,7 +25,7 @@ y_elem: 2
 z_elem: 60 
 z_max: 3e3 
 z_stretch: false
-dt: "200secs" 
+dt: "100secs" 
 t_end: "6hours" 
 dt_save_to_disk: "10mins"
 toml: [toml/diagnostic_edmfx_box.toml]

config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml

@@ -25,7 +25,7 @@ y_elem: 2
 z_elem: 30 
 z_max: 1500 
 z_stretch: false
-dt: 200secs
+dt: 100secs
 t_end: 4hours 
 dt_save_to_disk: 10mins
 toml: [toml/diagnostic_edmfx_box.toml]

config/model_configs/diagnostic_edmfx_rico_box.yml

@@ -25,7 +25,7 @@ y_elem: 2
 z_elem: 100
 z_max: 4e3
 z_stretch: false
-dt: 200secs
+dt: 100secs
 t_end: 8hours
 dt_save_to_disk: 10mins
 toml: [toml/diagnostic_edmfx_trmm_box.toml]

config/model_configs/diagnostic_edmfx_trmm_box.yml

@@ -26,7 +26,7 @@ y_elem: 2
 z_elem: 82
 z_max: 16400
 z_stretch: false
-dt: 300secs
+dt: 200secs
 t_end: 6hours
 dt_save_to_disk: 10mins
 FLOAT_TYPE: "Float64"

src/cache/diagnostic_edmf_precomputed_quantities.jl

@@ -229,6 +229,7 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
     FT = eltype(Y)
     n = n_mass_flux_subdomains(turbconv_model)
     ᶜz = Fields.coordinate_field(Y.c).z
+    ᶜdz = Fields.Δz_field(axes(Y.c))
     (; params) = p
     (; dt) = p
     (; ᶜΦ, ᶜρ_ref) = p.core
@@ -295,6 +296,7 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
         ts_prev_level = Fields.field_values(Fields.level(ᶜts, i - 1))
         p_prev_level = Fields.field_values(Fields.level(ᶜp, i - 1))
         z_prev_level = Fields.field_values(Fields.level(ᶜz, i - 1))
+        dz_prev_level = Fields.field_values(Fields.level(ᶜdz, i - 1))
 
         local_geometry_prev_level = Fields.field_values(
             Fields.level(Fields.local_geometry_field(Y.c), i - 1),
@@ -364,10 +366,29 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
                 ),
                 TD.relative_humidity(thermo_params, ts_prev_level),
                 FT(0),
-                dt,
                 p.atmos.edmfx_entr_model,
             )
 
+            # We don't have an upper limit to entrainment for the first level 
+            # (calculated at i=2), as the vertical at the first level is zero
+            if i == 2
+                @. entrʲ_prev_level = limit_entrainment(
+                    entrʲ_prev_level,
+                    draft_area(ρaʲ_prev_level, ρʲ_prev_level),
+                    dt,
+                )
+            else
+                @. entrʲ_prev_level = limit_entrainment(
+                    entrʲ_prev_level,
+                    draft_area(ρaʲ_prev_level, ρʲ_prev_level),
+                    get_physical_w(
+                        u³ʲ_prev_halflevel,
+                        local_geometry_prev_halflevel,
+                    ),
+                    dz_prev_level,
+                )
+            end
+
             # TODO: use updraft top instead of scale height
             @. nh_pressure³ʲ_prev_halflevel = ᶠupdraft_nh_pressure(
                 params,
@@ -493,10 +514,15 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
                 FT(0),
                 entrʲ_prev_level,
                 vert_div_level,
-                dt,
                 p.atmos.edmfx_detr_model,
             )
 
+            @. detrʲ_prev_level = limit_detrainment(
+                detrʲ_prev_level,
+                draft_area(ρaʲ_prev_level, ρʲ_prev_level),
+                dt,
+            )
+
             ρaʲu³ʲ_data = p.scratch.temp_data_level_2
             ρaʲu³ʲ_datah_tot = ρaʲu³ʲ_dataq_tot = p.scratch.temp_data_level_3
 

src/cache/prognostic_edmf_precomputed_quantities.jl

@@ -206,9 +206,13 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             get_physical_w(ᶜu, ᶜlg),
             TD.relative_humidity(thermo_params, ᶜts⁰),
             FT(0),
-            dt,
             p.atmos.edmfx_entr_model,
         )
+        @. ᶜentrʲs.:($$j) = limit_entrainment(
+            ᶜentrʲs.:($$j),
+            draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)),
+            dt,
+        )
         @. ᶜvert_div = ᶜdivᵥ(ᶠinterp(ᶜρʲs.:($$j)) * ᶠu³ʲs.:($$j)) / ᶜρʲs.:($$j)
         @. ᶜdetrʲs.:($$j) = detrainment(
             params,
@@ -226,9 +230,13 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             FT(0),
             ᶜentrʲs.:($$j),
             ᶜvert_div,
-            dt,
             p.atmos.edmfx_detr_model,
         )
+        @. ᶜdetrʲs.:($$j) = limit_detrainment(
+            ᶜdetrʲs.:($$j),
+            draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)),
+            dt,
+        )
     end
 
     # First order approximation: Use environmental mean fields.

src/prognostic_equations/edmfx_entr_detr.jl

@@ -33,7 +33,6 @@ function entrainment(
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
-    dt::FT,
     ::NoEntrainment,
 ) where {FT}
     return FT(0)
@@ -53,7 +52,6 @@ function entrainment(
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
-    dt::FT,
     ::PiGroupsEntrainment,
 ) where {FT}
     if ᶜaʲ <= FT(0)
@@ -74,15 +72,11 @@ function entrainment(
         Π₃ = sqrt(ᶜaʲ)
         Π₄ = ᶜRHʲ - ᶜRH⁰
         Π₆ = (ᶜz - z_sfc) / ref_H
-        entr = max(
-            0,
-            min(
-                abs(ᶜwʲ) / (ᶜz - z_sfc) * (
-                    -4.013288 - 0.000968 * Π₁ + 0.356974 * Π₃ - 0.403124 * Π₄ + 1.503261 * Π₆
-                ),
-                (1 - ᶜaʲ) / max(ᶜaʲ, eps(FT)) / dt,
-            ),
-        )
+        entr =
+            abs(ᶜwʲ) / (ᶜz - z_sfc) * (
+                -4.013288 - 0.000968 * Π₁ + 0.356974 * Π₃ - 0.403124 * Π₄ +
+                1.503261 * Π₆
+            )
 
         return entr
     end
@@ -102,7 +96,6 @@ function entrainment(
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
-    dt::FT,
     ::GeneralizedEntrainment,
 ) where {FT}
     turbconv_params = CAP.turbconv_params(params)
@@ -115,15 +108,9 @@ function entrainment(
     min_area_limiter =
         min_area_limiter_scale *
         exp(-min_area_limiter_power * (max(ᶜaʲ, 0) - a_min))
-    entr = max(
-        0,
-        min(
-            entr_inv_tau +
-            entr_coeff * abs(ᶜwʲ) / (ᶜz - z_sfc) +
-            min_area_limiter,
-            (1 - ᶜaʲ) / max(ᶜaʲ, eps(FT)) / dt,
-        ),
-    )
+    entr =
+        entr_inv_tau + entr_coeff * abs(ᶜwʲ) / (ᶜz - z_sfc) + min_area_limiter
+
     return entr
 end
 
@@ -141,7 +128,6 @@ function entrainment(
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
-    dt::FT,
     ::GeneralizedHarmonicsEntrainment,
 ) where {FT}
     turbconv_params = CAP.turbconv_params(params)
@@ -154,15 +140,9 @@ function entrainment(
     min_area_limiter =
         min_area_limiter_scale *
         exp(-min_area_limiter_power * (max(ᶜaʲ, 0) - a_min))
-    entr = max(
-        0,
-        min(
-            entr_inv_tau +
-            entr_coeff * abs(ᶜwʲ) / (ᶜz - z_sfc) +
-            min_area_limiter,
-            (1 - ᶜaʲ) / max(ᶜaʲ, eps(FT)) / dt,
-        ),
-    )
+    entr =
+        entr_inv_tau + entr_coeff * abs(ᶜwʲ) / (ᶜz - z_sfc) + min_area_limiter
+
     return entr * FT(2) * hm_limiter(ᶜaʲ)
 end
 
@@ -193,7 +173,6 @@ function detrainment(
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
-    dt::FT,
     ᶜentr::FT,
     ᶜvert_div::FT,
     ::NoDetrainment,
@@ -217,7 +196,6 @@ function detrainment(
     ᶜbuoy⁰::FT,
     ᶜentr::FT,
     ᶜvert_div::FT,
-    dt::FT,
     ::PiGroupsDetrainment,
 ) where {FT}
 
@@ -242,15 +220,11 @@ function detrainment(
         Π₃ = sqrt(ᶜaʲ)
         Π₄ = ᶜRHʲ - ᶜRH⁰
         Π₆ = (ᶜz - z_sfc) / ref_H
-        detr = max(
-            0,
-            min(
-                abs(ᶜwʲ) * (
-                    3.535208 + 0.598496 * Π₁ + 1.583348 * Π₃ + 0.046275 * Π₄ - 0.344836 * Π₆ + max_area_limiter
-                ),
-                1 / dt,
-            ),
-        )
+        detr =
+            abs(ᶜwʲ) * (
+                3.535208 + 0.598496 * Π₁ + 1.583348 * Π₃ + 0.046275 * Π₄ -
+                0.344836 * Π₆ + max_area_limiter
+            )
         return detr
     end
 end
@@ -271,7 +245,6 @@ function detrainment(
     ᶜbuoy⁰::FT,
     ᶜentr::FT,
     ᶜvert_div::FT,
-    dt::FT,
     ::GeneralizedDetrainment,
 ) where {FT}
     turbconv_params = CAP.turbconv_params(params)
@@ -286,17 +259,13 @@ function detrainment(
     max_area_limiter =
         max_area_limiter_scale *
         exp(-max_area_limiter_power * (a_max - min(ᶜaʲ, 1)))
-    detr = max(
-        min(
-            detr_inv_tau +
-            detr_coeff * abs(ᶜwʲ) +
-            detr_buoy_coeff * abs(min(ᶜbuoyʲ - ᶜbuoy⁰, 0)) /
-            max(eps(FT), abs(ᶜwʲ - ᶜw⁰)) - detr_vertdiv_coeff * ᶜvert_div +
-            max_area_limiter,
-            1 / dt,
-        ),
-        0,
-    )
+    detr =
+        detr_inv_tau +
+        detr_coeff * abs(ᶜwʲ) +
+        detr_buoy_coeff * abs(min(ᶜbuoyʲ - ᶜbuoy⁰, 0)) /
+        max(eps(FT), abs(ᶜwʲ - ᶜw⁰)) - detr_vertdiv_coeff * ᶜvert_div +
+        max_area_limiter
+
     return detr
 end
 
@@ -316,7 +285,6 @@ function detrainment(
     ᶜbuoy⁰::FT,
     ᶜentr::FT,
     ᶜvert_div::FT,
-    dt::FT,
     ::GeneralizedHarmonicsDetrainment,
 ) where {FT}
     turbconv_params = CAP.turbconv_params(params)
@@ -331,17 +299,13 @@ function detrainment(
     max_area_limiter =
         max_area_limiter_scale *
         exp(-max_area_limiter_power * (a_max - min(ᶜaʲ, 1)))
-    detr = max(
-        min(
-            detr_inv_tau +
-            detr_coeff * abs(ᶜwʲ) +
-            detr_buoy_coeff * abs(min(ᶜbuoyʲ - ᶜbuoy⁰, 0)) /
-            max(eps(FT), abs(ᶜwʲ - ᶜw⁰)) - detr_vertdiv_coeff * ᶜvert_div +
-            max_area_limiter,
-            1 / dt,
-        ),
-        0,
-    )
+    detr =
+        detr_inv_tau +
+        detr_coeff * abs(ᶜwʲ) +
+        detr_buoy_coeff * abs(min(ᶜbuoyʲ - ᶜbuoy⁰, 0)) /
+        max(eps(FT), abs(ᶜwʲ - ᶜw⁰)) - detr_vertdiv_coeff * ᶜvert_div +
+        max_area_limiter
+
     return detr * FT(2) * hm_limiter(ᶜaʲ)
 end
 
@@ -361,10 +325,9 @@ function detrainment(
     ᶜbuoy⁰::FT,
     ᶜentr::FT,
     ᶜvert_div::FT,
-    dt::FT,
     ::ConstantAreaDetrainment,
 ) where {FT}
-    detr = max(min(ᶜentr - ᶜvert_div, 1 / dt), 0)
+    detr = ᶜentr - ᶜvert_div
     return detr
 end
 
@@ -403,3 +366,9 @@ function edmfx_entr_detr_tendency!(
     end
     return nothing
 end
+
+limit_entrainment(entr::FT, a, dt) where {FT} =
+    max(min(entr, (1 - a) / max(a, eps(FT)) / dt), 0)
+limit_entrainment(entr, a, w, dz) = max(min(entr, w / dz), 0)
+limit_detrainment(detr, a, dt) = max(min(detr, 1 / dt), 0)
+limit_detrainment(detr, a, w, dz) = max(min(detr, w / dz), 0)