From aa82966bca932c61ac025fc13738c5d9c2c5fa67 Mon Sep 17 00:00:00 2001 From: yairchn Date: Mon, 16 Aug 2021 11:02:47 -0500 Subject: [PATCH] add RI and Pr functions WIP move pr scale to paramset WIP WIP debugging fix bug --- integration_tests/utils/generate_namelist.jl | 1 + integration_tests/utils/parameter_set.jl | 2 ++ src/Turbulence_PrognosticTKE.jl | 34 ++++++++------------ src/turbulence_functions.jl | 14 ++++++++ 4 files changed, 30 insertions(+), 21 deletions(-) diff --git a/integration_tests/utils/generate_namelist.jl b/integration_tests/utils/generate_namelist.jl index 74d842f47..b762bb304 100644 --- a/integration_tests/utils/generate_namelist.jl +++ b/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 diff --git a/integration_tests/utils/parameter_set.jl b/integration_tests/utils/parameter_set.jl index 5fb942524..f6eebef7b 100644 --- a/integration_tests/utils/parameter_set.jl +++ b/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 diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl index 66f9f6bd1..5883409a5 100755 --- a/src/Turbulence_PrognosticTKE.jl +++ b/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 diff --git a/src/turbulence_functions.jl b/src/turbulence_functions.jl index 5fc552a52..dc9a32d13 100644 --- a/src/turbulence_functions.jl +++ b/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