Update turbulent Prandtl number. Limit by Prandtl and not Richardson number

Author: trontrytel

.buildkite/Manifest-v1.11.toml

@@ -399,9 +399,9 @@ version = "0.2.13"
 
 [[deps.ClimaParams]]
 deps = ["TOML"]
-git-tree-sha1 = "db2d3e069ed4c12c96a4ae071ab5fc37da1019d4"
+git-tree-sha1 = "7613d47b4f96307845cbe377dede19efa3256cfb"
 uuid = "5c42b081-d73a-476f-9059-fd94b934656c"
-version = "0.10.28"
+version = "0.10.29"
 
 [[deps.ClimaReproducibilityTests]]
 deps = ["OrderedCollections", "PrettyTables"]

Project.toml

@@ -43,7 +43,7 @@ AtmosphericProfilesLibrary = "0.1.7"
 ClimaComms = "0.6.6"
 ClimaCore = "0.14.24"
 ClimaDiagnostics = "0.2.12"
-ClimaParams = "0.10.27"
+ClimaParams = "0.10.29"
 ClimaTimeSteppers = "0.8.2"
 ClimaUtilities = "0.1.22"
 CloudMicrophysics = "0.22.9"

src/cache/diagnostic_edmf_precomputed_quantities.jl

@@ -983,12 +983,9 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_closures!
     @. ᶜstrain_rate_norm = norm_sqr(ᶜstrain_rate)
 
     ᶜprandtl_nvec = p.scratch.ᶜtemp_scalar
-    @. ᶜprandtl_nvec = turbulent_prandtl_number(
-        params,
-        obukhov_length,
-        ᶜlinear_buoygrad,
-        ᶜstrain_rate_norm,
-    )
+    @. ᶜprandtl_nvec =
+        turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm)
+
     ᶜtke_exch = p.scratch.ᶜtemp_scalar_2
     @. ᶜtke_exch = 0
     # using ᶜu⁰ would be more correct, but this is more consistent with the

src/cache/prognostic_edmf_precomputed_quantities.jl

@@ -466,7 +466,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
         )
 
         # The buoyancy term in the nonhydrostatic pressure closure is always applied
-        # for prognostic edmf. The tendency is combined with the buoyancy term in the 
+        # for prognostic edmf. The tendency is combined with the buoyancy term in the
         # updraft momentum equation in `edmfx_sgs_vertical_advection_tendency!`. This
         # term is still calculated here as it is used explicitly in the TKE equation.
         @. ᶠnh_pressure₃_buoyʲs.:($$j) = ᶠupdraft_nh_pressure_buoyancy(
@@ -501,12 +501,9 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
     @. ᶜstrain_rate_norm = norm_sqr(ᶜstrain_rate)
 
     ᶜprandtl_nvec = p.scratch.ᶜtemp_scalar
-    @. ᶜprandtl_nvec = turbulent_prandtl_number(
-        params,
-        obukhov_length,
-        ᶜlinear_buoygrad,
-        ᶜstrain_rate_norm,
-    )
+    @. ᶜprandtl_nvec =
+        turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm)
+
     ᶜtke_exch = p.scratch.ᶜtemp_scalar_2
     @. ᶜtke_exch = 0
     for j in 1:n

src/parameters/Parameters.jl

@@ -30,7 +30,7 @@ Base.@kwdef struct TurbulenceConvectionParameters{FT, VFT1, VFT2} <: ATCP
     static_stab_coeff::FT
     Prandtl_number_scale::FT
     Prandtl_number_0::FT
-    Ri_crit::FT
+    Pr_max::FT
     smin_ub::FT
     smin_rm::FT
     min_updraft_top::FT

src/parameters/create_parameters.jl

@@ -202,7 +202,6 @@ function TurbulenceConvectionParameters(
         :EDMF_max_area => :max_area,
         :mixing_length_smin_rm => :smin_rm,
         :entr_coeff => :entr_coeff,
-        :mixing_length_Ri_crit => :Ri_crit,
         :detr_coeff => :detr_coeff,
         :EDMF_surface_area => :surface_area,
         :entr_param_vec => :entr_param_vec,
@@ -219,6 +218,7 @@ function TurbulenceConvectionParameters(
         :pressure_normalmode_drag_coeff => :pressure_normalmode_drag_coeff,
         :mixing_length_Prandtl_number_scale => :Prandtl_number_scale,
         :mixing_length_Prandtl_number_0 => :Prandtl_number_0,
+        :mixing_length_Prandtl_maximum => :Pr_max,
         :mixing_length_static_stab_coeff => :static_stab_coeff,
         :pressure_normalmode_buoy_coeff1 =>
             :pressure_normalmode_buoy_coeff1,

src/prognostic_equations/edmfx_closures.jl

@@ -281,42 +281,58 @@ function mixing_length(
 end
 
 """
-    turbulent_prandtl_number(params, obukhov_length, ᶜRi_grad)
+    turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm)
 
 where:
-- `params`: set with model parameters
-- `obukhov_length`: surface Monin Obukhov length
-- `ᶜRi_grad`: gradient Richardson number
+- `params`: Parameters set
+- `ᶜlinear_buoygrad`: N^2, Brunt-Väisälä frequency squared [1/s^2].
+- `ᶜstrain_rate_norm`: Frobenius norm of strain rate tensor, |S| [1/s].
+
+Returns the turbulent Prandtl number based on the gradient Richardson number.
+
+The formula implemented is from Li et al. (JAS 2015, DOI: 10.1175/JAS-D-14-0335.1, their Eq. 39),
+with a reformulation and correction of an algebraic error in their expression:
+
+    Pr_t(Ri) = (X + sqrt(max(X^2 - 4*Pr_n*Ri, 0))) / 2
 
-Returns the turbulent Prandtl number give the obukhov length sign and
-the gradient Richardson number, which is calculated from the linearized
-buoyancy gradient and shear production.
+where X = Pr_n + ω_pr * Ri and Ri = N^2 / max(2*|S|, eps)
+using parameters Pr_n = Prandtl_number_0, ω_pr = Prandtl_number_scale.
+This formula applies in both stable (Ri > 0) and unstable (Ri < 0) conditions.
+The returned turbulent Prandtl number is limited by Pr_max parameter.
 """
-function turbulent_prandtl_number(
-    params,
-    obukhov_length,
-    ᶜlinear_buoygrad,
-    ᶜstrain_rate_norm,
-)
+function turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm)
     FT = eltype(params)
     turbconv_params = CAP.turbconv_params(params)
-    Ri_c = CAP.Ri_crit(turbconv_params)
-    ω_pr = CAP.Prandtl_number_scale(turbconv_params)
-    Pr_n = CAP.Prandtl_number_0(turbconv_params)
-    ᶜRi_grad = min(ᶜlinear_buoygrad / max(2 * ᶜstrain_rate_norm, eps(FT)), Ri_c)
-    if obukhov_length > 0 && ᶜRi_grad > 0 #stable
-        # CSB (Dan Li, 2019, eq. 75), where ω_pr = ω_1 + 1 = 53.0 / 13.0
-        prandtl_nvec =
-            Pr_n * (
-                2 * ᶜRi_grad / (
-                    1 + ω_pr * ᶜRi_grad -
-                    sqrt((1 + ω_pr * ᶜRi_grad)^2 - 4 * ᶜRi_grad)
-                )
-            )
-    else
-        prandtl_nvec = Pr_n
-    end
-    return prandtl_nvec
+    eps_FT = eps(FT)
+
+    # Parameters from CliMAParams
+    Pr_n = CAP.Prandtl_number_0(turbconv_params) # Neutral Prandtl number
+    ω_pr = CAP.Prandtl_number_scale(turbconv_params) # Prandtl number scale coefficient
+    Pr_max = CAP.Pr_max(turbconv_params) # Maximum Prandtl number limit
+
+    # Calculate the raw gradient Richardson number
+    # Using the definition Ri = N^2 / (2*|S|)
+    ᶜshear_term_safe = max(2 * ᶜstrain_rate_norm, eps_FT)
+    ᶜRi_grad = ᶜlinear_buoygrad / ᶜshear_term_safe
+
+    # --- Apply the Pr_t(Ri) formula valid for stable and unstable conditions ---
+
+    # Calculate the intermediate term X = Pr_n + ω_pr * Ri
+    X = Pr_n + ω_pr * ᶜRi_grad
+
+    # Calculate the discriminant term: (Pr_n + ω_pr*Ri)^2 - 4*Pr_n*Ri = X^2 - 4*Pr_n*Ri
+    discriminant = X * X - 4 * Pr_n * ᶜRi_grad
+    # Ensure the discriminant is non-negative before taking the square root
+    discriminant_safe = max(discriminant, FT(0))
+
+    # Calculate the Prandtl number using the positive root solution of the quadratic eq.
+    # Pr_t = ( X + sqrt(discriminant_safe) ) / 2
+    prandtl_nvec = (X + sqrt(discriminant_safe)) / 2
+
+    # Optional safety: ensure Pr_t is not excessively small or negative,
+    # though the formula should typically yield positive values if Pr_n > 0.
+    # Also ensure that it's not larger than the Pr_max parameter.
+    return min(max(prandtl_nvec, eps_FT), Pr_max)
 end
 
 edmfx_filter_tendency!(Yₜ, Y, p, t, turbconv_model) = nothing

src/prognostic_equations/gm_sgs_closures.jl

@@ -51,12 +51,8 @@ NVTX.@annotate function compute_gm_mixing_length!(ᶜmixing_length, Y, p)
     @. ᶜstrain_rate_norm = norm_sqr(ᶜstrain_rate)
 
     ᶜprandtl_nvec = p.scratch.ᶜtemp_scalar_2
-    @. ᶜprandtl_nvec = turbulent_prandtl_number(
-        params,
-        obukhov_length,
-        ᶜlinear_buoygrad,
-        ᶜstrain_rate_norm,
-    )
+    @. ᶜprandtl_nvec =
+        turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm)
 
     @. ᶜmixing_length = smagorinsky_lilly_length(
         CAP.c_smag(params),