From 3e87cde9ca5dc88ab351704d0cc66913eefcccb0 Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Sun, 15 Aug 2021 15:10:22 -0700 Subject: [PATCH] Use dispatch for entr-detr, + DocStringExtensions --- Project.toml | 1 + docs/Project.toml | 1 + src/TurbulenceConvection.jl | 1 + src/Turbulence_PrognosticTKE.jl | 56 ++++++------- src/closures/entr_detr.jl | 30 +++---- .../nondimensional_exchange_functions.jl | 20 ++--- src/types.jl | 84 +++++++++++-------- src/utility_functions.jl | 1 - 8 files changed, 105 insertions(+), 89 deletions(-) diff --git a/Project.toml b/Project.toml index 7c0923f4e4..ca6516dfd0 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" diff --git a/docs/Project.toml b/docs/Project.toml index fb4cfb469f..6e53463a69 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,7 @@ [deps] ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53" +DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" diff --git a/src/TurbulenceConvection.jl b/src/TurbulenceConvection.jl index 7d498a6095..6d61f3371b 100644 --- a/src/TurbulenceConvection.jl +++ b/src/TurbulenceConvection.jl @@ -3,6 +3,7 @@ module TurbulenceConvection using StaticArrays using StatsBase using LinearAlgebra +import DocStringExtensions import Dierckx import Statistics diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl index 459fb434e1..4774a75788 100755 --- a/src/Turbulence_PrognosticTKE.jl +++ b/src/Turbulence_PrognosticTKE.jl @@ -900,27 +900,17 @@ function get_env_covar_from_GMV( end function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::CasesBase) - εδ_model = entr_detr_model() - sa = eos_struct() - quadrature_order = 3 grid = get_grid(self) param_set = parameter_set(self) - - upd_cloud_diagnostics(self.UpdVar, reference_state(self)) - - # TODO: merge this with grid.cinterior - kc_surf = kc_surface(grid) - kc_toa = kc_top_of_atmos(grid) - cinterior = kc_surf:kc_toa - ref_state = reference_state(self) - l = zeros(2) + + upd_cloud_diagnostics(self.UpdVar, ref_state) @inbounds for k in real_center_indicies(grid) @inbounds for i in xrange(self.n_updrafts) # entrainment - w_upd = interpf2c(self.UpdVar.W.values, grid, k, i) if self.UpdVar.Area.values[i, k] > 0.0 + w_upd = interpf2c(self.UpdVar.W.values, grid, k, i) # compute dMdz at half levels gmv_w_k = interpf2c(GMV.W.values, grid, k) gmv_w_km = interpf2c(GMV.W.values, grid, k - 1) @@ -930,24 +920,24 @@ function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariabl dMdz = (M - Mm) * grid.dzi w_min = 0.001 - εδ_model.b_upd = self.UpdVar.B.values[i, k] - εδ_model.b_env = self.EnvVar.B.values[k] - εδ_model.w_upd = interpf2c(self.UpdVar.W.values, grid, k, i) - εδ_model.w_env = interpf2c(self.EnvVar.W.values, grid, k) - εδ_model.a_upd = self.UpdVar.Area.values[i, k] - εδ_model.a_env = self.EnvVar.Area.values[k] - εδ_model.ql_upd = self.UpdVar.QL.values[i, k] - εδ_model.ql_env = self.EnvVar.QL.values[k] - εδ_model.RH_upd = self.UpdVar.RH.values[i, k] - εδ_model.RH_env = self.EnvVar.RH.values[k] - εδ_model.M = M - εδ_model.dMdz = dMdz - εδ_model.tke = self.EnvVar.TKE.values[k] - εδ_model.n_up = self.n_updrafts - εδ_model.ρ = ref_state.rho0[k] - εδ_model.R_up = self.pressure_plume_spacing[i] - - self.entr_sc[i, k], self.detr_sc[i, k], self.frac_turb_entr[i, k], self.horiz_K_eddy[i, k] = entr_detr( + εδ_model = MoistureDeficitEntr(; + q_liq_up = self.UpdVar.QL.values[i, k], + q_liq_en = self.EnvVar.QL.values[k], + w_up = interpf2c(self.UpdVar.W.values, grid, k, i), + w_en = interpf2c(self.EnvVar.W.values, grid, k), + b_up = self.UpdVar.B.values[i, k], + b_en = self.EnvVar.B.values[k], + tke = self.EnvVar.TKE.values[k], + dMdz = dMdz, + M = M, + a_up = self.UpdVar.Area.values[i, k], + a_en = self.EnvVar.Area.values[k], + R_up = self.pressure_plume_spacing[i], + RH_up = self.UpdVar.RH.values[i, k], + RH_en = self.EnvVar.RH.values[k], + ) + + er = entr_detr( param_set, w_min, self.sorting_power, @@ -960,6 +950,10 @@ function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariabl self.turbulent_entrainment_factor, εδ_model, ) + self.entr_sc[i, k] = er.ε_dyn + self.detr_sc[i, k] = er.δ_dyn + self.frac_turb_entr[i, k] = er.ε_turb + self.horiz_K_eddy[i, k] = er.K_ε else self.entr_sc[i, k] = 0.0 self.detr_sc[i, k] = 0.0 diff --git a/src/closures/entr_detr.jl b/src/closures/entr_detr.jl index 0557c6909f..2dd360d3de 100644 --- a/src/closures/entr_detr.jl +++ b/src/closures/entr_detr.jl @@ -1,4 +1,5 @@ #### Entrainment-Detrainment kernels + """ entr_detr( param_set, @@ -21,35 +22,35 @@ Cohen et al. (JAMES, 2020), given: - `β`: sorting power for moist mixing - `c_δ`: detrainment factor - `c_div`: divergence factor for bubble case (zero otherwise) - - `c_μ`: logisitc function scale + - `c_μ`: logisitc function scale - `c_μ0`: logisitc function timescale - `χ_upd`: updraft mixing fraction - `c_λ`: tke scale factor - `c_εt`: turbulent entrainment factor - - `εδ_model`: entrainment detrainment model type + - `εδ_model`: a [`MoistureDeficitEntr`](@ref) """ -function entr_detr(param_set, w_min, β, c_δ, c_div, c_μ, c_μ0, χ_upd, c_λ, c_εt, εδ_model) +function entr_detr(param_set, w_min, β, c_δ, c_div, c_μ, c_μ0, χ_upd, c_λ, c_εt, εδ_model::MoistureDeficitEntr) - l = zeros(2) c_ε = CPEDMF.c_ε(param_set) - if (εδ_model.ql_upd + εδ_model.ql_env) == 0.0 + if (εδ_model.q_liq_up + εδ_model.q_liq_en) == 0.0 c_δ = 0.0 end - Δw = εδ_model.w_upd - εδ_model.w_env + Δw = εδ_model.w_up - εδ_model.w_en if Δw < 0.0 Δw -= w_min else Δw += w_min end - Δb = (εδ_model.b_upd - εδ_model.b_env) + Δb = (εδ_model.b_up - εδ_model.b_en) - D_ε, D_δ, M_δ, M_ε = nondimensional_exchange_functions(c_ε, c_δ, c_μ, c_μ0, β, χ_upd, Δw, Δb, εδ_model) + D_ε, D_δ, M_δ, M_ε = nondimensional_exchange_functions(param_set, c_δ, c_μ, c_μ0, β, χ_upd, Δw, Δb, εδ_model) - l[1] = c_λ * abs(Δb / sqrt(εδ_model.tke + 1e-8)) - l[2] = abs(Δb / Δw) + l_1 = c_λ * abs(Δb / sqrt(εδ_model.tke + 1e-8)) + l_2 = abs(Δb / Δw) + l = (l_1, l_2) λ = lamb_smooth_minimum(l, 0.1, 0.0005) MdMdz = max(εδ_model.dMdz / max(εδ_model.M, 1e-12), 0.0) @@ -58,9 +59,9 @@ function entr_detr(param_set, w_min, β, c_δ, c_div, c_μ, c_μ0, χ_upd, c_λ, # turbulent entrainment - k_ε = εδ_model.a_upd * c_εt * sqrt(max(εδ_model.tke, 0.0)) * εδ_model.R_up - if εδ_model.w_upd * εδ_model.a_upd > 0.0 - ε_turb = (2.0 / εδ_model.R_up^2.0) * k_ε / (εδ_model.w_upd * εδ_model.a_upd) + K_ε = εδ_model.a_up * c_εt * sqrt(max(εδ_model.tke, 0.0)) * εδ_model.R_up + if εδ_model.w_up * εδ_model.a_up > 0.0 + ε_turb = (2.0 / εδ_model.R_up^2.0) * K_ε / (εδ_model.w_up * εδ_model.a_up) else ε_turb = 0.0 end @@ -68,5 +69,6 @@ function entr_detr(param_set, w_min, β, c_δ, c_div, c_μ, c_μ0, χ_upd, c_λ, ε_dyn = λ / Δw * (c_ε * D_ε + c_δ * M_ε) + MdMdz * c_div δ_dyn = λ / Δw * (c_ε * D_δ + c_δ * M_δ) + MdMdz * c_div - return ε_dyn, δ_dyn, ε_turb, k_ε + + return EntrDetr(ε_dyn, δ_dyn, ε_turb, K_ε) end diff --git a/src/closures/nondimensional_exchange_functions.jl b/src/closures/nondimensional_exchange_functions.jl index b4f84f0edb..8da1d3c221 100644 --- a/src/closures/nondimensional_exchange_functions.jl +++ b/src/closures/nondimensional_exchange_functions.jl @@ -1,6 +1,6 @@ """ nondimensional_exchange_functions( - c_ε, + param_set, c_δ, c_μ, c_μ0, @@ -9,29 +9,29 @@ Δw, Δb, εδ_model, - ) + ) Returns the nondimensional entrainment and detrainment functions following Cohen et al. (JAMES, 2020), given: - - `c_ε`, entrainment factor - - `c_δ`, detrainment factor - - `c_μ`, logistic function factor - - `c_μ0`, logistic function timescale (sec) + - `param_set` a parameter set + - `c_δ`, detrainment factor + - `c_μ`, logistic function factor + - `c_μ0`, logistic function timescale (sec) - `β`, sorting power - `χ_up`, updraft mixing fraction - `Δw`, updraft - environment vertical velocity differnce - `Δb`, updraft - environment buoynacy differnce - `εδ_model`, entrainment detrainment model type """ -function nondimensional_exchange_functions(c_ε, c_δ, c_μ, c_μ0, β, χ_upd, Δw, Δb, εδ_model) +function nondimensional_exchange_functions(param_set, c_δ, c_μ, c_μ0, β, χ_upd, Δw, Δb, εδ_model::MoistureDeficitEntr) # should be: c_δ = sign(condensate(ts_en) + condensate(ts_up[i])) * entr.c_δ - μ_ij = (χ_upd - εδ_model.a_upd / (εδ_model.a_upd + εδ_model.a_env)) * Δb / Δw + μ_ij = (χ_upd - εδ_model.a_up / (εδ_model.a_up + εδ_model.a_en)) * Δb / Δw D_ε = 1.0 / (1.0 + exp(-c_μ / c_μ0 * μ_ij)) D_δ = 1.0 / (1.0 + exp(c_μ / c_μ0 * μ_ij)) - M_δ = (max((εδ_model.RH_upd / 100.0)^β - (εδ_model.RH_env / 100.0)^β, 0.0))^(1.0 / β) - M_ε = (max((εδ_model.RH_env / 100.0)^β - (εδ_model.RH_upd / 100.0)^β, 0.0))^(1.0 / β) + M_δ = (max((εδ_model.RH_up / 100.0)^β - (εδ_model.RH_en / 100.0)^β, 0.0))^(1.0 / β) + M_ε = (max((εδ_model.RH_en / 100.0)^β - (εδ_model.RH_up / 100.0)^β, 0.0))^(1.0 / β) return D_ε, D_δ, M_δ, M_ε end; diff --git a/src/types.jl b/src/types.jl index 560ce04119..b27ef4f212 100644 --- a/src/types.jl +++ b/src/types.jl @@ -19,41 +19,59 @@ Base.@kwdef mutable struct mph_struct thl_rain_src::Float64 = 0 qr_src::Float64 = 0 end + +""" + EntrDetr + +$(DocStringExtensions.FIELDS) +""" +Base.@kwdef struct EntrDetr{FT} + "Dynamical entrainment" + ε_dyn::FT + "Dynamical detrainment" + δ_dyn::FT + "Turbulent entrainment" + ε_turb::FT + "Horizontal eddy-diffusivity" + K_ε::FT +end + """ - Entrainment detrainment model type - - `εδ_model.b_upd`: updraft buoyancy - - `εδ_model.b_env`: environment vertical velocity - - `εδ_model.w_upd`: updraft vertical velocity - - `εδ_model.w_env`: environment area fraction - - `εδ_model.a_upd`: updraft area fraction - - `εδ_model.a_env`: environment buoyancy - - `εδ_model.ql_up`: updraft liquid water - - `εδ_model.ql_env`: environment liquid water - - `εδ_model.RH_upd`: updraft relative humidity - - `εδ_model.RH_env`: environment relative humidity - - `εδ_model.M`: updraft momentum - - `εδ_model.dMdz`: updraft momentum divergence - - `εδ_model.tke`: env TKE - - `εδ_model.N_up`: total number of updrafts - - `εδ_model.ρ`: referance density + MoistureDeficitEntr + +My entrainment detrainment model + +$(DocStringExtensions.FIELDS) """ -Base.@kwdef mutable struct entr_detr_model - b_upd::Float64 = 0 - b_env::Float64 = 0 - w_upd::Float64 = 0 - w_env::Float64 = 0 - a_upd::Float64 = 0 - a_env::Float64 = 0 - ql_upd::Float64 = 0 - ql_env::Float64 = 0 - RH_upd::Float64 = 0 - RH_env::Float64 = 0 - M::Float64 = 0 - dMdz::Float64 = 0 - tke::Float64 = 0 - n_up::Float64 = 0 - ρ::Float64 = 0 - R_up::Float64 = 0 +Base.@kwdef struct MoistureDeficitEntr{FT} + "updraft liquid water" + q_liq_up::FT + "environment liquid water" + q_liq_en::FT + "updraft vertical velocity" + w_up::FT + "environment vertical velocity" + w_en::FT + "updraft buoyancy" + b_up::FT + "environment buoyancy" + b_en::FT + "environment tke" + tke::FT + "updraft momentum divergence" + dMdz::FT + "updraft momentum" + M::FT + "updraft area fraction" + a_up::FT + "environment area fraction" + a_en::FT + "pressure plume spacing" + R_up::FT + "updraft relative humidity" + RH_up::FT + "environment relative humidity" + RH_en::FT end struct RainVariable{T} diff --git a/src/utility_functions.jl b/src/utility_functions.jl index 7011998f45..9d1fa2f59a 100644 --- a/src/utility_functions.jl +++ b/src/utility_functions.jl @@ -22,7 +22,6 @@ function logistic(x, slope, mid) end function lamb_smooth_minimum(l, lower_bound, upper_bound) - leng = size(l) x_min = minimum(l) λ_0 = max(x_min * lower_bound / real(LambertW.lambertw(2.0 / MathConstants.e)), upper_bound)