add mass flux divergence to detrainment

config/model_configs/prognostic_edmfx_bomex_box.yml

@@ -23,7 +23,7 @@ y_elem: 2
 z_elem: 60
 z_stretch: false
 perturb_initstate: false
-dt: "20secs"
+dt: "10secs"
 t_end: "6hours"
 dt_save_state_to_disk: "30mins"
 toml: [toml/prognostic_edmfx_bomex_box.toml]

post_processing/ci_plots.jl

@@ -955,6 +955,7 @@ EDMFBoxPlots = Union{
     Val{:prognostic_edmfx_gabls_box},
     Val{:prognostic_edmfx_bomex_fixtke_box},
     Val{:prognostic_edmfx_bomex_box},
+    Val{:prognostic_edmfx_bomex_column},
     Val{:prognostic_edmfx_bomex_stretched_box},
     Val{:prognostic_edmfx_dycoms_rf01_box},
     Val{:prognostic_edmfx_rico_column},

src/cache/diagnostic_edmf_precomputed_quantities.jl

@@ -483,6 +483,7 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
                 p_prev_level,
                 ρ_prev_level,
                 buoyancy_flux_sfc_halflevel,
+                ρaʲ_prev_level,
                 draft_area(ρaʲ_prev_level, ρʲ_prev_level),
                 get_physical_w(
                     u³ʲ_prev_halflevel,
@@ -498,6 +499,7 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
                 FT(0),
                 entrʲ_prev_level,
                 vert_div_level,
+                FT(0), # mass flux divergence is not implemented for diagnostic edmf
                 p.atmos.edmfx_detr_model,
             )
 

src/cache/prognostic_edmf_precomputed_quantities.jl

@@ -64,8 +64,6 @@ function set_prognostic_edmf_precomputed_quantities_draft_and_bc!(Y, p, ᶠuₕ
     (; ᶜspecific, ᶜp, ᶜh_tot, ᶜK) = p.precomputed
     (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶠKᵥʲs, ᶜtsʲs, ᶜρʲs) = p.precomputed
     (; ustar, obukhov_length, buoyancy_flux) = p.precomputed.sfc_conditions
-    ᶜinterp_lb = Operators.LeftBiasedF2C()
-    ᶜinterp_rb = Operators.RightBiasedF2C()
 
     for j in 1:n
         ᶜuʲ = ᶜuʲs.:($j)
@@ -195,6 +193,7 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
     ᶜlg = Fields.local_geometry_field(Y.c)
 
     ᶜvert_div = p.scratch.ᶜtemp_scalar
+    ᶜmassflux_vert_div = p.scratch.ᶜtemp_scalar_2
     for j in 1:n
         # entrainment/detrainment
         @. ᶜentrʲs.:($$j) = entrainment(
@@ -219,13 +218,15 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             dt,
         )
         @. ᶜvert_div = ᶜdivᵥ(ᶠinterp(ᶜρʲs.:($$j)) * ᶠu³ʲs.:($$j)) / ᶜρʲs.:($$j)
+        @. ᶜmassflux_vert_div = ᶜdivᵥ(ᶠinterp(Y.c.sgsʲs.:($$j).ρa) * ᶠu³ʲs.:($$j))
         @. ᶜdetrʲs.:($$j) = detrainment(
             params,
             ᶜz,
             z_sfc,
             ᶜp,
             Y.c.ρ,
             buoyancy_flux,
+            Y.c.sgsʲs.:($$j).ρa,
             draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)),
             get_physical_w(ᶜuʲs.:($$j), ᶜlg),
             TD.relative_humidity(thermo_params, ᶜtsʲs.:($$j)),
@@ -235,6 +236,7 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             FT(0),
             ᶜentrʲs.:($$j),
             ᶜvert_div,
+            ᶜmassflux_vert_div,
             p.atmos.edmfx_detr_model,
         )
         @. ᶜdetrʲs.:($$j) = limit_detrainment(

src/prognostic_equations/edmfx_entr_detr.jl

@@ -114,50 +114,19 @@ function entrainment(
     return entr
 end
 
-function entrainment(
-    params,
-    ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    buoy_flux_surface::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    ::GeneralizedHarmonicsEntrainment,
-) where {FT}
-    turbconv_params = CAP.turbconv_params(params)
-    entr_inv_tau = CAP.entr_tau(turbconv_params)
-    entr_coeff = CAP.entr_coeff(turbconv_params)
-    min_area_limiter_scale = CAP.min_area_limiter_scale(turbconv_params)
-    min_area_limiter_power = CAP.min_area_limiter_power(turbconv_params)
-    a_min = CAP.min_area(turbconv_params)
-
-    min_area_limiter =
-        min_area_limiter_scale *
-        exp(-min_area_limiter_power * (max(ᶜaʲ, 0) - a_min))
-    entr =
-        entr_inv_tau + entr_coeff * abs(ᶜwʲ) / (ᶜz - z_sfc) + min_area_limiter
-
-    return entr * FT(2) * hm_limiter(ᶜaʲ)
-end
-
 """
    Return detrainment rate [1/s].
 
    Inputs (everything defined on cell centers):
    - params set with model parameters
    - ᶜz, z_sfc, ᶜp, ᶜρ, - grid-scale height, surface height, grid-scale pressure and density
    - buoy_flux_surface - buoyancy flux at the surface
-   - ᶜaʲ, ᶜwʲ, ᶜRHʲ, ᶜbuoyʲ - updraft area, physical vertical velocity,
+   - ᶜρaʲ, ᶜaʲ, ᶜwʲ, ᶜRHʲ, ᶜbuoyʲ - updraft effective density, updraft area, physical vertical velocity,
                                    relative humidity and buoyancy
    - ᶜw⁰, ᶜRH⁰, ᶜbuoy⁰ - environment physical vertical velocity,
                               relative humidity and buoyancy
-   - dt - timestep
+   - ᶜentr - entrainment rate
+   - ᶜvert_div,ᶜmassflux_vert_div - vertical divergence, vertical mass flux divergence
 """
 function detrainment(
     params,
@@ -166,6 +135,7 @@ function detrainment(
     ᶜp::FT,
     ᶜρ::FT,
     buoy_flux_surface::FT,
+    ᶜρaʲ::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
     ᶜRHʲ::FT,
@@ -175,6 +145,7 @@ function detrainment(
     ᶜbuoy⁰::FT,
     ᶜentr::FT,
     ᶜvert_div::FT,
+    ᶜmassflux_vert_div::FT,
     ::NoDetrainment,
 ) where {FT}
     return FT(0)
@@ -187,6 +158,7 @@ function detrainment(
     ᶜp::FT,
     ᶜρ::FT,
     buoy_flux_surface::FT,
+    ᶜρaʲ::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
     ᶜRHʲ::FT,
@@ -196,6 +168,7 @@ function detrainment(
     ᶜbuoy⁰::FT,
     ᶜentr::FT,
     ᶜvert_div::FT,
+    ᶜmassflux_vert_div::FT,
     ::PiGroupsDetrainment,
 ) where {FT}
 
@@ -236,6 +209,7 @@ function detrainment(
     ᶜp::FT,
     ᶜρ::FT,
     buoy_flux_surface::FT,
+    ᶜρaʲ::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
     ᶜRHʲ::FT,
@@ -245,6 +219,7 @@ function detrainment(
     ᶜbuoy⁰::FT,
     ᶜentr::FT,
     ᶜvert_div::FT,
+    ᶜmassflux_vert_div::FT,
     ::GeneralizedDetrainment,
 ) where {FT}
     turbconv_params = CAP.turbconv_params(params)
@@ -259,12 +234,17 @@ function detrainment(
     max_area_limiter =
         max_area_limiter_scale *
         exp(-max_area_limiter_power * (a_max - min(ᶜaʲ, 1)))
-    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
+
+    if ᶜρaʲ <= 0
+        detr = 0
+    else
+        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 * min(ᶜmassflux_vert_div, 0) / ᶜρaʲ +
+            max_area_limiter
+    end
 
     return detr
 end
@@ -276,6 +256,7 @@ function detrainment(
     ᶜp::FT,
     ᶜρ::FT,
     buoy_flux_surface::FT,
+    ᶜρaʲ::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
     ᶜRHʲ::FT,
@@ -285,46 +266,7 @@ function detrainment(
     ᶜbuoy⁰::FT,
     ᶜentr::FT,
     ᶜvert_div::FT,
-    ::GeneralizedHarmonicsDetrainment,
-) where {FT}
-    turbconv_params = CAP.turbconv_params(params)
-    detr_inv_tau = CAP.detr_tau(turbconv_params)
-    detr_coeff = CAP.detr_coeff(turbconv_params)
-    detr_buoy_coeff = CAP.detr_buoy_coeff(turbconv_params)
-    detr_vertdiv_coeff = CAP.detr_vertdiv_coeff(turbconv_params)
-    max_area_limiter_scale = CAP.max_area_limiter_scale(turbconv_params)
-    max_area_limiter_power = CAP.max_area_limiter_power(turbconv_params)
-    a_max = CAP.max_area(turbconv_params)
-
-    max_area_limiter =
-        max_area_limiter_scale *
-        exp(-max_area_limiter_power * (a_max - min(ᶜaʲ, 1)))
-    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
-
-function detrainment(
-    params,
-    ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    buoy_flux_surface::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    ᶜentr::FT,
-    ᶜvert_div::FT,
+    ᶜmassflux_vert_div::FT,
     ::ConstantAreaDetrainment,
 ) where {FT}
     detr = ᶜentr - ᶜvert_div

toml/prognostic_edmfx_bomex_box.toml

@@ -8,13 +8,13 @@ value = 1.0e-5
 value = 0.7
 
 [entr_inv_tau]
-value = 0
+value = 0.0002
 
 [entr_coeff]
-value = 0.1
+value = 0.2
 
 [min_area_limiter_scale]
-value = 0.001
+value = 0
 
 [min_area_limiter_power]
 value = 10
@@ -23,13 +23,13 @@ value = 10
 value = 0
 
 [detr_coeff]
-value = 2e-3
+value = 0
 
 [detr_buoy_coeff]
-value = 0.3
+value = 0
 
 [detr_vertdiv_coeff]
-value = 0.6
+value = 1
 
 [max_area_limiter_scale]
-value = 0.01
+value = 0

toml/prognostic_edmfx_simpleplume.toml

@@ -35,13 +35,13 @@ value = 0
 type = "float"
 
 [detr_buoy_coeff]
-value = 0.3
+value = 0
 type = "float"
 
 [detr_vertdiv_coeff]
-value = 0
+value = 1
 type = "float"
 
 [max_area_limiter_scale]
-value = 0.01
+value = 0
 type = "float"