add Ri and Pr functions

integration_tests/utils/generate_namelist.jl

@@ -69,6 +69,7 @@ function default_namelist(case_name::String)
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["tke_diss_coeff"] = 0.22
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["static_stab_coeff"] = 0.4
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["lambda_stab"] = 0.9
+    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["Prandtl_number_scale"] = 53.0 / 13.0
     # entrainment
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_factor"] = 0.13
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["detrainment_factor"] = 0.51

integration_tests/utils/parameter_set.jl

@@ -14,6 +14,7 @@ CLIMAParameters.Atmos.EDMF.α_b(ps::EarthParameterSet) = ps.nt.α_b
 CLIMAParameters.Atmos.EDMF.α_a(ps::EarthParameterSet) = ps.nt.α_a
 CLIMAParameters.Atmos.EDMF.α_d(ps::EarthParameterSet) = ps.nt.α_d
 CLIMAParameters.Atmos.EDMF.H_up_min(ps::EarthParameterSet) = ps.nt.H_up_min
+CLIMAParameters.Atmos.EDMF.ω_pr(ps::EarthParameterSet) = ps.nt.ω_pr
 
 #! format: off
 function create_parameter_set(namelist)
@@ -27,6 +28,7 @@ function create_parameter_set(namelist)
         α_a = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "pressure_normalmode_adv_coeff"),
         α_d = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "pressure_normalmode_drag_coeff"),
         H_up_min = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "min_updraft_top"),
+        ω_pr = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "Prandtl_number_scale"),
     )
     return EarthParameterSet(nt)
 end

src/Turbulence_PrognosticTKE.jl

@@ -348,7 +348,9 @@ function compute_mixing_length(self, obukhov_length, ustar, GMV::GridMeanVariabl
     kc_surf = kc_surface(grid)
     tau = get_mixing_tau(self.zi, self.wstar)
     l = zeros(3)
-    m_eps = 1.0e-9 # Epsilon to avoid zero
+    ϵ = 1.0e-9 # Epsilon to avoid zero
+    ω_pr = CPEDMF.ω_pr(param_set)
+
     @inbounds for k in real_center_indicies(grid)
         z_ = grid.z_half[k]
         # kz scale (surface layer)
@@ -367,7 +369,7 @@ function compute_mixing_length(self, obukhov_length, ustar, GMV::GridMeanVariabl
         ∇U = c∇(U_cut, grid, k; bottom = SetGradient(0), top = SetGradient(0))
         ∇V = c∇(V_cut, grid, k; bottom = SetGradient(0), top = SetGradient(0))
         ∇w = ∇f2c(w_dual, grid, k; bottom = SetGradient(0), top = SetGradient(0))
-        shear2 = ∇U^2 + ∇V^2 + ∇w^2
+        Shear² = ∇U^2 + ∇V^2 + ∇w^2
         qt_dry = self.EnvThermo.qt_dry[k]
         th_dry = self.EnvThermo.th_dry[k]
         t_cloudy = self.EnvThermo.t_cloudy[k]
@@ -407,26 +409,16 @@ function compute_mixing_length(self, obukhov_length, ustar, GMV::GridMeanVariabl
         )
 
         # Partial buoyancy gradients
-        grad_b_thl = grad_thl * d_buoy_thetal_total
-        grad_b_qt = grad_qt * d_buoy_qt_total
-        ri_grad = min(grad_b_thl / max(shear2, m_eps) + grad_b_qt / max(shear2, m_eps), 0.25)
-
-        # Turbulent Prandtl number
-        if obukhov_length > 0.0 && ri_grad > 0.0 #stable
-            # CSB (Dan Li, 2019), with Pr_neutral=0.74 and w1=40.0/13.0
-            self.prandtl_nvec[k] =
-                prandtl_number(self) * (
-                    2.0 * ri_grad /
-                    (1.0 + (53.0 / 13.0) * ri_grad - sqrt((1.0 + (53.0 / 13.0) * ri_grad)^2.0 - 4.0 * ri_grad))
-                )
-        else
-            self.prandtl_nvec[k] = prandtl_number(self)
-        end
+        ∂b∂θ_l = grad_thl * d_buoy_thetal_total
+        ∂b∂q_tot = grad_qt * d_buoy_qt_total
+
+        ∇Ri = gradient_Richardson_number(∂b∂θ_l, ∂b∂q_tot, Shear², ϵ)
+        self.prandtl_nvec[k] = turbulent_Prandtl_number(obukhov_length, ∇Ri, prandtl_number(self), ω_pr)
 
         # Production/destruction terms
         a =
             self.tke_ed_coeff *
-            (shear2 - grad_b_thl / self.prandtl_nvec[k] - grad_b_qt / self.prandtl_nvec[k]) *
+            (Shear² - ∂b∂θ_l / self.prandtl_nvec[k] - ∂b∂q_tot / self.prandtl_nvec[k]) *
             sqrt(self.EnvVar.TKE.values[k])
         # Dissipation term
         c_neg = self.tke_diss_coeff * self.EnvVar.TKE.values[k] * sqrt(self.EnvVar.TKE.values[k])
@@ -446,9 +438,9 @@ function compute_mixing_length(self, obukhov_length, ustar, GMV::GridMeanVariabl
                 wc_env / (1.0 - self.UpdVar.Area.bulkvalues[k])
         end
 
-        if abs(a) > m_eps && 4.0 * a * c_neg > -self.b[k] * self.b[k]
+        if abs(a) > ϵ && 4.0 * a * c_neg > -self.b[k] * self.b[k]
             self.l_entdet[k] = max(-self.b[k] / 2.0 / a + sqrt(self.b[k] * self.b[k] + 4.0 * a * c_neg) / 2.0 / a, 0.0)
-        elseif abs(a) < m_eps && abs(self.b[k]) > m_eps
+        elseif abs(a) < ϵ && abs(self.b[k]) > ϵ
             self.l_entdet[k] = c_neg / self.b[k]
         end
         l3 = self.l_entdet[k]
@@ -493,7 +485,7 @@ function compute_mixing_length(self, obukhov_length, ustar, GMV::GridMeanVariabl
 
         j = 1
         while (j <= length(l))
-            if l[j] < m_eps || l[j] > 1.0e6
+            if l[j] < ϵ || l[j] > 1.0e6
                 l[j] = 1.0e6
             end
             j += 1

src/turbulence_functions.jl

@@ -101,3 +101,17 @@ function set_cloudbase_flag(ql, current_flag)
     end
     return new_flag
 end
+
+function gradient_Richardson_number(∂b∂θ_l, ∂b∂q_tot, Shear², ϵ)
+    return min(∂b∂θ_l / max(Shear², ϵ) + ∂b∂q_tot / max(Shear², ϵ), 0.25)
+end
+
+function turbulent_Prandtl_number(obukhov_length, ∇Ri, Pr, ω_pr)
+    if obukhov_length > 0.0 && ∇Ri > 0.0 #stable
+        # CSB (Dan Li, 2019), with Pr_neutral=0.74 and w1=40.0/13.0
+        prandtl_nvec = Pr * (2 * ∇Ri / (1 + ω_pr * ∇Ri - sqrt((1 + ω_pr * ∇Ri)^2 - 4 * ∇Ri)))
+    else
+        prandtl_nvec = Pr
+    end
+    return prandtl_nvec
+end