add minimum mixing length

CliMA · Oct 23, 2023 · c0f886b · c0f886b
1 parent 2f11b1d
commit c0f886b
Show file tree

Hide file tree

Showing 2 changed files with 14 additions and 10 deletions.
diff --git a/src/prognostic_equations/edmfx_closures.jl b/src/prognostic_equations/edmfx_closures.jl
@@ -194,10 +194,10 @@ function mixing_length(
     # kz scale (surface layer)
     if obukhov_length < 0.0 #unstable
         l_W =
-            vkc * (ᶜz - z_sfc) / (sqrt(sfc_tke / ustar / ustar) * c_m) *
+            vkc * (ᶜz - z_sfc) / (max(sqrt(sfc_tke / ustar / ustar), eps(FT)) * c_m) *
             min((1 - 100 * (ᶜz - z_sfc) / obukhov_length)^FT(0.2), 1 / vkc)
     else # neutral or stable
-        l_W = vkc * (ᶜz - z_sfc) / (sqrt(sfc_tke / ustar / ustar) * c_m)
+        l_W = vkc * (ᶜz - z_sfc) / (max(sqrt(sfc_tke / ustar / ustar), eps(FT)) * c_m)
     end
 
     # compute l_TKE - the production-dissipation balanced length scale
@@ -235,15 +235,19 @@ function mixing_length(
         l_smag = c_smag * ᶜdz
     end
 
+    l_smag = FT(10)
     # add limiters
-    l = SA.SVector(
-        (l_N < eps(FT) || l_N > l_z) ? l_z : l_N,
-        (l_TKE < eps(FT) || l_TKE > l_z) ? l_z : l_TKE,
-        (l_W < eps(FT) || l_W > l_z) ? l_z : l_W,
-    )
+    # l = SA.SVector(
+    #     (l_N < eps(FT) || l_N > l_z) ? l_z : l_N,
+    #     (l_TKE < eps(FT) || l_TKE > l_z) ? l_z : l_TKE,
+    #     (l_W < eps(FT) || l_W > l_z) ? l_z : l_W,
+    # )
+    # l_limited = lamb_smooth_minimum(l, smin_ub, smin_rm)
+    l = SA.SVector(l_N, l_TKE, l_W)
+    l_limited = max(l_smag, min(lamb_smooth_minimum(l, smin_ub, smin_rm), l_z))
     # get soft minimum
     # TODO: limit it with l_smag
-    return lamb_smooth_minimum(l, smin_ub, smin_rm)
+    return l_limited
 end
 
 """

diff --git a/src/prognostic_equations/edmfx_sgs_flux.jl b/src/prognostic_equations/edmfx_sgs_flux.jl
@@ -169,7 +169,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
     if p.atmos.edmfx_sgs_diffusive_flux
         # energy
         ᶠρaK_h = p.ᶠtemp_scalar
-        @. ᶠρaK_h[colidx] = ᶠinterp(ᶜρa⁰[colidx]) * max(ᶠinterp(ᶜK_h[colidx]), FT(1))
+        @. ᶠρaK_h[colidx] = ᶠinterp(ᶜρa⁰[colidx]) * max(ᶠinterp(ᶜK_h[colidx]), FT(0))
 
         ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(
             top = Operators.SetValue(C3(FT(0))),
@@ -192,7 +192,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
 
         # momentum
         ᶠρaK_u = p.ᶠtemp_scalar
-        @. ᶠρaK_u[colidx] = ᶠinterp(ᶜρa⁰[colidx]) * max(ᶠinterp(ᶜK_u[colidx]), FT(1))
+        @. ᶠρaK_u[colidx] = ᶠinterp(ᶜρa⁰[colidx]) * max(ᶠinterp(ᶜK_u[colidx]), FT(0))
         ᶠstrain_rate = p.ᶠtemp_UVWxUVW
         compute_strain_rate_face!(ᶠstrain_rate[colidx], ᶜu⁰[colidx])
         @. Yₜ.c.uₕ[colidx] -= C12(