Merge #2331

2331: add a detrainment closure that produces constant area fraction r=szy21 a=szy21 Co-authored-by: Zhaoyi Shen <11598433+szy21@users.noreply.github.com>
CliMA · Nov 3, 2023 · 66f8b05 · 66f8b05
2 parents 2d6eb60 + 65f1fca
commit 66f8b05
Show file tree

Hide file tree

Showing 11 changed files with 99 additions and 8 deletions.
diff --git a/config/model_configs/prognostic_edmfx_bomex_box.yml b/config/model_configs/prognostic_edmfx_bomex_box.yml
@@ -7,7 +7,7 @@ surface_setup: "Bomex"
 turbconv: "prognostic_edmfx"
 edmfx_upwinding: first_order
 edmfx_entr_model: "Generalized"
-edmfx_detr_model: "Generalized"
+edmfx_detr_model: "ConstantArea"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
@@ -27,7 +27,7 @@ perturb_initstate: false
 dt: "5secs"
 t_end: "6hours"
 dt_save_to_disk: "10mins"
-toml: [toml/prognostic_edmfx_box.toml]
+toml: [toml/prognostic_edmfx_bomex_box.toml]
 diagnostics:
   - short_name: ts
     period: 10mins

diff --git a/config/model_configs/prognostic_edmfx_bomex_column.yml b/config/model_configs/prognostic_edmfx_bomex_column.yml
@@ -7,7 +7,7 @@ surface_setup: "Bomex"
 turbconv: "prognostic_edmfx"
 edmfx_upwinding: first_order
 edmfx_entr_model: "Generalized"
-edmfx_detr_model: "Generalized"
+edmfx_detr_model: "ConstantArea"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
@@ -22,4 +22,4 @@ perturb_initstate: false
 dt: "5secs"
 t_end: "6hours"
 dt_save_to_disk: "10mins"
-toml: [toml/prognostic_edmfx_box.toml]
+toml: [toml/prognostic_edmfx_bomex_box.toml]
diff --git a/config/model_configs/prognostic_edmfx_bomex_tke_box.yml b/config/model_configs/prognostic_edmfx_bomex_tke_box.yml
@@ -7,7 +7,7 @@ surface_setup: "Bomex"
 turbconv: "prognostic_edmfx"
 edmfx_upwinding: first_order
 edmfx_entr_model: "Generalized"
-edmfx_detr_model: "Generalized"
+edmfx_detr_model: "ConstantArea"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
@@ -27,4 +27,4 @@ perturb_initstate: false
 dt: "5secs"
 t_end: "6hours"
 dt_save_to_disk: "10mins"
-toml: [toml/prognostic_edmfx_box.toml]
+toml: [toml/prognostic_edmfx_bomex_box.toml]
diff --git a/config/model_configs/prognostic_edmfx_bomex_tke_column.yml b/config/model_configs/prognostic_edmfx_bomex_tke_column.yml
@@ -7,7 +7,7 @@ surface_setup: "Bomex"
 turbconv: "prognostic_edmfx"
 edmfx_upwinding: first_order
 edmfx_entr_model: "Generalized"
-edmfx_detr_model: "Generalized"
+edmfx_detr_model: "ConstantArea"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
@@ -22,4 +22,4 @@ perturb_initstate: false
 dt: "5secs"
 t_end: "6hours"
 dt_save_to_disk: "10mins"
-toml: [toml/prognostic_edmfx_box.toml]
+toml: [toml/prognostic_edmfx_bomex_box.toml]
diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl
@@ -392,6 +392,8 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
                 ),
                 TD.relative_humidity(thermo_params, ts_prev_level),
                 FT(0),
+                FT(0), # ᶜentr, not implemented
+                FT(0), # ᶜvert_div, not implemented
                 dt,
                 p.atmos.edmfx_detr_model,
             )

diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
@@ -44,6 +44,8 @@ function precomputed_quantities(Y, atmos)
         !(atmos.moisture_model isa DryModel) &&
         atmos.energy_form isa TotalEnergy
     ) || !(atmos.turbconv_model isa PrognosticEDMFX)
+    @assert !(atmos.edmfx_detr_model isa ConstantAreaDetrainment) ||
+            !(atmos.turbconv_model isa DiagnosticEDMFX)
     TST = thermo_state_type(atmos.moisture_model, FT)
     SCT = SurfaceConditions.surface_conditions_type(atmos, FT)
     n = n_mass_flux_subdomains(atmos.turbconv_model)

diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl
@@ -192,6 +192,7 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
     ᶜdz = Fields.Δz_field(axes(Y.c))
     ᶜlg = Fields.local_geometry_field(Y.c)
 
+    ᶜvert_div = p.scratch.ᶜtemp_scalar
     for j in 1:n
         @. ᶜentrʲs.:($$j) = entrainment(
             params,
@@ -210,6 +211,7 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             dt,
             p.atmos.edmfx_entr_model,
         )
+        @. ᶜvert_div = ᶜdivᵥ(ᶠinterp(ᶜρʲs.:($$j)) * ᶠu³ʲs.:($$j)) / ᶜρʲs.:($$j)
         @. ᶜdetrʲs.:($$j) = detrainment(
             params,
             ᶜz,
@@ -224,6 +226,8 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             get_physical_w(ᶜu, ᶜlg),
             TD.relative_humidity(thermo_params, ᶜts⁰),
             FT(0),
+            ᶜentrʲs.:($$j),
+            ᶜvert_div,
             dt,
             p.atmos.edmfx_detr_model,
         )

diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl
@@ -184,6 +184,8 @@ function detrainment(
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
     dt::FT,
+    ᶜentr::FT,
+    ᶜvert_div::FT,
     ::NoDetrainment,
 ) where {FT}
     return FT(0)
@@ -203,6 +205,8 @@ function detrainment(
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
+    ᶜentr::FT,
+    ᶜvert_div::FT,
     dt::FT,
     ::PiGroupsDetrainment,
 ) where {FT}
@@ -255,6 +259,8 @@ function detrainment(
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
+    ᶜentr::FT,
+    ᶜvert_div::FT,
     dt::FT,
     ::GeneralizedDetrainment,
 ) where {FT}
@@ -294,6 +300,8 @@ function detrainment(
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
+    ᶜentr::FT,
+    ᶜvert_div::FT,
     dt::FT,
     ::GeneralizedHarmonicsDetrainment,
 ) where {FT}
@@ -319,6 +327,29 @@ function detrainment(
     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,
+    dt::FT,
+    ::ConstantAreaDetrainment,
+) where {FT}
+    detr = max(min(ᶜentr - ᶜvert_div, 1 / dt), 0)
+    return detr
+end
+
 edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model) = nothing
 
 function edmfx_entr_detr_tendency!(

diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl
@@ -376,6 +376,8 @@ function get_detrainment_model(parsed_args)
         GeneralizedDetrainment()
     elseif detr_model == "GeneralizedHarmonics"
         GeneralizedHarmonicsDetrainment()
+    elseif detr_model == "ConstantArea"
+        ConstantAreaDetrainment()
     else
         error("Invalid detr_model $(detr_model)")
     end

diff --git a/src/solver/types.jl b/src/solver/types.jl
@@ -188,6 +188,7 @@ struct NoDetrainment <: AbstractDetrainmentModel end
 struct PiGroupsDetrainment <: AbstractDetrainmentModel end
 struct GeneralizedDetrainment <: AbstractDetrainmentModel end
 struct GeneralizedHarmonicsDetrainment <: AbstractDetrainmentModel end
+struct ConstantAreaDetrainment <: AbstractDetrainmentModel end
 
 abstract type AbstractQuadratureType end
 struct LogNormalQuad <: AbstractQuadratureType end

diff --git a/toml/prognostic_edmfx_bomex_box.toml b/toml/prognostic_edmfx_bomex_box.toml
@@ -0,0 +1,49 @@
+[EDMF_surface_area]
+alias = "surface_area"
+value = 0.1
+type = "float"
+
+[EDMF_min_area]
+alias = "min_area"
+value = 1.0e-5
+type = "float"
+
+[EDMF_max_area]
+alias = "max_area"
+value = 0.6
+type = "float"
+
+[entr_tau]
+alias = "entr_tau"
+value = 0.002
+type = "float"
+
+[entr_coeff]
+alias = "entr_coeff"
+value = 0
+type = "float"
+
+[min_area_limiter_scale]
+alias = "min_area_limiter_scale"
+value = 0
+type = "float"
+
+[detr_tau]
+alias = "detr_tau"
+value = 0
+type = "float"
+
+[detr_coeff]
+alias = "detr_coeff"
+value = 0
+type = "float"
+
+[detr_buoy_coeff]
+alias = "detr_buoy_coeff"
+value = 0
+type = "float"
+
+[max_area_limiter_scale]
+alias = "max_area_limiter_scale"
+value = 0
+type = "float"