From 250845c062e6175caf95bd00a24e2c108b4bf6d3 Mon Sep 17 00:00:00 2001
From: Zhaoyi Shen <11598433+szy21@users.noreply.github.com>
Date: Wed, 25 Oct 2023 11:45:10 -0700
Subject: [PATCH] limit area fraction in entrainment and detrainment limiter

---
 .../diagnostic_edmfx_dycoms_rf01_box.yml         |  2 +-
 src/prognostic_equations/edmfx_entr_detr.jl      | 16 ++++++++++------
 2 files changed, 11 insertions(+), 7 deletions(-)

diff --git a/config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml b/config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml
index 5d216e618f..bcca068b72 100644
--- a/config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml
+++ b/config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml
@@ -23,7 +23,7 @@ y_elem: 2
 z_elem: 30 
 z_max: 1500 
 z_stretch: false
-dt: 6secs 
+dt: 10secs 
 t_end: 4hours 
 dt_save_to_disk: 2mins
 toml: [toml/diagnostic_edmfx_box.toml]
diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl
index 2a5c7f12c5..e968841bb6 100644
--- a/src/prognostic_equations/edmfx_entr_detr.jl
+++ b/src/prognostic_equations/edmfx_entr_detr.jl
@@ -112,7 +112,8 @@ function entrainment(
     turbconv_params = CAP.turbconv_params(params)
     a_min = TCP.min_area(turbconv_params)
     min_area_limiter =
-        min_area_limiter_scale * exp(-min_area_limiter_power * (ᶜaʲ - a_min))
+        min_area_limiter_scale *
+        exp(-min_area_limiter_power * (max(ᶜaʲ, 0) - a_min))
     entr = min(
         entr_inv_tau + entr_coeff * abs(ᶜwʲ) / (ᶜz - z_sfc) + min_area_limiter,
         1 / dt,
@@ -144,7 +145,8 @@ function entrainment(
     turbconv_params = CAP.turbconv_params(params)
     a_min = TCP.min_area(turbconv_params)
     min_area_limiter =
-        min_area_limiter_scale * exp(-min_area_limiter_power * (ᶜaʲ - a_min))
+        min_area_limiter_scale *
+        exp(-min_area_limiter_power * (max(ᶜaʲ, 0) - a_min))
     entr = min(
         entr_inv_tau + entr_coeff * abs(ᶜwʲ) / (ᶜz - z_sfc) + min_area_limiter,
         1 / dt,
@@ -262,12 +264,13 @@ function detrainment(
     turbconv_params = CAP.turbconv_params(params)
     a_max = TCP.max_area(turbconv_params)
     max_area_limiter =
-        max_area_limiter_scale * exp(-max_area_limiter_power * (a_max - ᶜaʲ))
+        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(eps(FT), abs(ᶜwʲ - ᶜw⁰)) +
         max_area_limiter,
         1 / dt,
     )
@@ -299,12 +302,13 @@ function detrainment(
     turbconv_params = CAP.turbconv_params(params)
     a_max = TCP.max_area(turbconv_params)
     max_area_limiter =
-        max_area_limiter_scale * exp(-max_area_limiter_power * (a_max - ᶜaʲ))
+        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(eps(FT), abs(ᶜwʲ - ᶜw⁰)) +
         max_area_limiter,
         1 / dt,
     )