add vertical divergence to detrainment closure

CliMA · Nov 10, 2023 · 2a1c314 · 2a1c314
1 parent e64ca3c
commit 2a1c314
Show file tree

Hide file tree

Showing 5 changed files with 91 additions and 17 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: "ConstantArea"
+edmfx_detr_model: "Generalized"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
@@ -27,4 +27,4 @@ perturb_initstate: false
 dt: "10secs"
 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/parameters/Parameters.jl b/src/parameters/Parameters.jl
@@ -36,6 +36,7 @@ Base.@kwdef struct TurbulenceConvectionParameters{FT} <: ATCP
     detr_tau::FT
     detr_coeff::FT
     detr_buoy_coeff::FT
+    detr_vertdiv_coeff::FT
     min_area_limiter_scale::FT
     min_area_limiter_power::FT
     max_area_limiter_scale::FT

diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl
@@ -142,6 +142,7 @@ function entrainment(
     turbconv_params = CAP.turbconv_params(params)
     entr_inv_tau = CAP.entr_tau(turbconv_params)
     entr_coeff = CAP.entr_coeff(turbconv_params)
+    entr_vertdiv_coeff = CAP.entr_vertdiv_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)
@@ -268,20 +269,24 @@ function detrainment(
     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 = min(
-        detr_inv_tau +
-        detr_coeff * abs(ᶜwʲ) +
-        detr_buoy_coeff * abs(min(ᶜbuoyʲ - ᶜbuoy⁰, 0)) /
-        max(eps(FT), abs(ᶜwʲ - ᶜw⁰)) +
-        max_area_limiter,
-        1 / dt,
+    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,
     )
     return detr
 end
@@ -309,20 +314,24 @@ function detrainment(
     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 = min(
-        detr_inv_tau +
-        detr_coeff * abs(ᶜwʲ) +
-        detr_buoy_coeff * abs(min(ᶜbuoyʲ - ᶜbuoy⁰, 0)) /
-        max(eps(FT), abs(ᶜwʲ - ᶜw⁰)) +
-        max_area_limiter,
-        1 / dt,
+    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,
     )
     return detr * FT(2) * hm_limiter(ᶜaʲ)
 end

diff --git a/toml/prognostic_edmfx_bomex_box.toml b/toml/prognostic_edmfx_bomex_box.toml
@@ -0,0 +1,59 @@
+[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.7
+type = "float"
+
+[entr_tau]
+alias = "entr_tau"
+value = 0
+type = "float"
+
+[entr_coeff]
+alias = "entr_coeff"
+value = 0.1
+type = "float"
+
+[min_area_limiter_scale]
+alias = "min_area_limiter_scale"
+value = 0.001
+type = "float"
+
+[min_area_limiter_power]
+alias = "min_area_limiter_power"
+value = 10
+type = "float"
+
+[detr_tau]
+alias = "detr_tau"
+value = 0
+type = "float"
+
+[detr_coeff]
+alias = "detr_coeff"
+value = 2e-3
+type = "float"
+
+[detr_buoy_coeff]
+alias = "detr_buoy_coeff"
+value = 0.3
+type = "float"
+
+[detr_vertdiv_coeff]
+alias = "detr_vertdiv_coeff"
+value = 0.6
+type = "float"
+
+[max_area_limiter_scale]
+alias = "max_area_limiter_scale"
+value = 0.01
+type = "float"
diff --git a/toml/prognostic_edmfx_dycoms_rf01_box.toml b/toml/prognostic_edmfx_dycoms_rf01_box.toml
@@ -40,14 +40,19 @@ type = "float"
 
 [detr_coeff]
 alias = "detr_coeff"
-value = 2e-3
+value = 1e-3
 type = "float"
 
 [detr_buoy_coeff]
 alias = "detr_buoy_coeff"
 value = 0.3
 type = "float"
 
+[detr_vertdiv_coeff]
+alias = "detr_vertdiv_coeff"
+value = 0.5
+type = "float"
+
 [max_area_limiter_scale]
 alias = "max_area_limiter_scale"
 value = 0.01