Add (1 - a)^c limiter to PiGroups and Generalized entrainment functions

Author: costachris

src/parameters/Parameters.jl

@@ -44,6 +44,7 @@ Base.@kwdef struct TurbulenceConvectionParameters{FT, VFT1, VFT2} <: ATCP
     detr_vertdiv_coeff::FT
     entr_param_vec::VFT1
     turb_entr_param_vec::VFT2
+    entr_mult_limiter_coeff::FT
     detr_massflux_vertdiv_coeff::FT
     min_area_limiter_scale::FT
     min_area_limiter_power::FT

src/parameters/create_parameters.jl

@@ -177,6 +177,7 @@ function TurbulenceConvectionParameters(
         :EDMF_surface_area => :surface_area,
         :entr_param_vec => :entr_param_vec,
         :turb_entr_param_vec => :turb_entr_param_vec,
+        :entr_mult_limiter_coeff => :entr_mult_limiter_coeff,
         :minimum_updraft_top => :min_updraft_top,
         :mixing_length_eddy_viscosity_coefficient => :tke_ed_coeff,
         :mixing_length_smin_ub => :smin_ub,

src/prognostic_equations/edmfx_entr_detr.jl

@@ -57,6 +57,8 @@ function entrainment(
     ::PiGroupsEntrainment,
 )
     FT = eltype(thermo_params)
+    entr_mult_limiter_coeff = CAP.entr_mult_limiter_coeff(turbconv_params)
+
     if ᶜaʲ <= FT(0)
         return FT(0)
     else
@@ -68,14 +70,15 @@ function entrainment(
         # non-dimensional pi-groups
         Π₁ = (ᶜz - z_sfc) * (ᶜbuoyʲ - ᶜbuoy⁰) / ((ᶜwʲ - ᶜw⁰)^2 + eps(FT)) / 100
         Π₂ = max(ᶜtke⁰, 0) / ((ᶜwʲ - ᶜw⁰)^2 + eps(FT)) / 2
-        Π₃ = sqrt(ᶜaʲ)
+        Π₃ = sqrt(max(ᶜaʲ, 0))
         Π₄ = ᶜRHʲ - ᶜRH⁰
         Π₅ = (ᶜz - z_sfc) / ref_H
         # Π₁, Π₂ are unbounded, so clip values that blow up
         Π₁ = min(max(Π₁, -1), 1)
         Π₂ = min(max(Π₂, -1), 1)
 
         entr =
+            ((FT(1) - min(ᶜaʲ, FT(1)))^entr_mult_limiter_coeff) *
             abs(ᶜwʲ - ᶜw⁰) / (ᶜz - z_sfc) * (
                 entr_param_vec[1] * abs(Π₁) +
                 entr_param_vec[2] * abs(Π₂) +
@@ -107,6 +110,7 @@ function entrainment(
     ::GeneralizedEntrainment,
 )
     FT = eltype(thermo_params)
+    entr_mult_limiter_coeff = CAP.entr_mult_limiter_coeff(turbconv_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)
@@ -117,9 +121,11 @@ function entrainment(
         min_area_limiter_scale *
         exp(-min_area_limiter_power * (max(ᶜaʲ, 0) - a_min))
     entr =
-        entr_inv_tau +
-        entr_coeff * abs(ᶜwʲ - ᶜw⁰) / (ᶜz - z_sfc) +
-        min_area_limiter
+        ((FT(1) - min(ᶜaʲ, FT(1)))^entr_mult_limiter_coeff) * (
+            entr_inv_tau +
+            entr_coeff * abs(ᶜwʲ - ᶜw⁰) / (ᶜz - z_sfc) +
+            min_area_limiter
+        )
 
     return max(entr, 0)
 end
@@ -245,7 +251,7 @@ function detrainment(
         # non-dimensional pi-groups
         Π₁ = (ᶜz - z_sfc) * (ᶜbuoyʲ - ᶜbuoy⁰) / ((ᶜwʲ - ᶜw⁰)^2 + eps(FT)) / 100
         Π₂ = max(ᶜtke⁰, 0) / ((ᶜwʲ - ᶜw⁰)^2 + eps(FT)) / 2
-        Π₃ = sqrt(ᶜaʲ)
+        Π₃ = sqrt(max(ᶜaʲ, 0))
         Π₄ = ᶜRHʲ - ᶜRH⁰
         Π₅ = (ᶜz - z_sfc) / ref_H
         # Π₁, Π₂ are unbounded, so clip values that blow up