read entrainment parameters from param_set

integration_tests/utils/generate_namelist.jl

@@ -80,6 +80,7 @@ function default_namelist(case_name::String)
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_sigma"] = 10.0
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_scale"] = 0.004
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["sorting_power"] = 2.0
+    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["min_upd_velocity"] = 0.001
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["aspect_ratio"] = 0.2
     # pressure 
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["min_updraft_top"] = 500.0

integration_tests/utils/parameter_set.jl

@@ -9,12 +9,19 @@ end
 
 CLIMAParameters.Planet.MSLP(ps::EarthParameterSet) = ps.nt.MSLP
 CLIMAParameters.Atmos.Microphysics.τ_cond_evap(ps::EarthParameterSet) = ps.nt.τ_cond_evap
-CLIMAParameters.Atmos.EDMF.c_ε(ps::EarthParameterSet) = ps.nt.c_ε
-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
+CLIMAParameters.Atmos.EDMF.c_ε(ps::EarthParameterSet) = ps.nt.c_ε # factor multiplyer for dry term in entrainment/detrainment
+CLIMAParameters.Atmos.EDMF.α_b(ps::EarthParameterSet) = ps.nt.α_b # factor multiplyer for pressure buoyancy terms (effective buoyancy is (1-α_b))
+CLIMAParameters.Atmos.EDMF.α_a(ps::EarthParameterSet) = ps.nt.α_a # factor multiplyer for pressure advection
+CLIMAParameters.Atmos.EDMF.α_d(ps::EarthParameterSet) = ps.nt.α_d # factor multiplyer for pressure drag
+CLIMAParameters.Atmos.EDMF.H_up_min(ps::EarthParameterSet) = ps.nt.H_up_min # minimum updarft top to avoid zero devision in pressure drag and turb-entr
+CLIMAParameters.Atmos.EDMF.ω_pr(ps::EarthParameterSet) = ps.nt.ω_pr # experiment scale factor for turbulent Prandtl number
+CLIMAParameters.Atmos.EDMF.c_δ(ps::EarthParameterSet) = ps.nt.c_δ # factor multiplyer for moist term in entrainment/detrainment
+CLIMAParameters.Atmos.EDMF.β(ps::EarthParameterSet) = ps.nt.β # sorting power for ad-hoc moisture detrainment function
+CLIMAParameters.Atmos.EDMF.χ(ps::EarthParameterSet) = ps.nt.χ # fraction of updraft air for buoyancy mixing in entrainment/detrainment (0≤χ≤1)
+CLIMAParameters.Atmos.EDMF.c_t(ps::EarthParameterSet) = ps.nt.c_t # factor multiplyer for turbulent term in entrainment/detrainment
+CLIMAParameters.Atmos.EDMF.c_λ(ps::EarthParameterSet) = ps.nt.c_λ # scaling factor for TKE in entrainment scale calculations
+CLIMAParameters.Atmos.EDMF.w_min(ps::EarthParameterSet) = ps.nt.w_min # minimum updraft velocity to aviod zero division in b/w²
+CLIMAParameters.Atmos.EDMF.μ_0(ps::EarthParameterSet) = ps.nt.μ_0 # dimentional scale logistic function in the dry term in entrainment/detrainment
 
 #! format: off
 function create_parameter_set(namelist)
@@ -29,6 +36,14 @@ function create_parameter_set(namelist)
         α_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"),
+        c_δ = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "detrainment_factor"),
+        β = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "sorting_power"),
+        χ = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "updraft_mixing_frac"),
+        c_t = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "turbulent_entrainment_factor"),
+        c_λ = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "entrainment_smin_tke_coeff"),
+        w_min = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "min_upd_velocity"),
+        μ_0 = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "entrainment_scale"),
+        μ = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "entrainment_sigma"),
     )
     return EarthParameterSet(nt)
 end

src/ClimaParams.jl

@@ -17,5 +17,6 @@ const APS = CP.AbstractEarthParameterSet
 
 """ divergence factor for bubble case (zero otherwise) """
 entrainment_massflux_div_factor(ps::APS) = ps.nt.c_div
+entrainment_sigma(ps::APS) = ps.nt.μ
 
 end

src/Turbulence_PrognosticTKE.jl

@@ -926,18 +926,7 @@ function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariabl
                     RH_en = self.EnvVar.RH.values[k],
                 )
 
-                er = entr_detr(
-                    param_set,
-                    w_min,
-                    self.sorting_power,
-                    self.detrainment_factor,
-                    self.entrainment_sigma,
-                    self.entrainment_scale,
-                    self.updraft_mixing_frac,
-                    self.entrainment_smin_tke_coeff,
-                    self.turbulent_entrainment_factor,
-                    εδ_model,
-                )
+                er = entr_detr(param_set, εδ_model)
                 self.entr_sc[i, k] = er.ε_dyn
                 self.detr_sc[i, k] = er.δ_dyn
                 self.frac_turb_entr[i, k] = er.ε_turb

src/closures/entr_detr.jl

@@ -3,35 +3,25 @@
 """
     entr_detr(
         param_set,
-        w_min,
-        β,
-        c_δ,
-        c_μ,
-        c_μ0,
-        χ_upd,
-        c_λ,
         εδ_model,
     )
-
 Returns the dynamic entrainment and detrainment rates,
 as well as the turbulent entrainment rate, following
 Cohen et al. (JAMES, 2020), given:
  - `param_set`: parameter set
- - `w_min`: minimum veritcal velocity
- - `β`: sorting power for moist mixing
- - `c_δ`: detrainment factor
- - `c_μ`: logisitc function scale
- - `c_μ0`: logisitc function timescale
- - `χ_upd`: updraft mixing fraction
- - `c_λ`: tke scale factor
- - `c_εt`: turbulent entrainment factor
  - `εδ_model`: a [`MoistureDeficitEntr`](@ref)
 """
-function entr_detr(param_set, w_min, β, c_δ, c_μ, c_μ0, χ_upd, c_λ, c_εt, εδ_model::MoistureDeficitEntr)
+function entr_detr(param_set, εδ_model::MoistureDeficitEntr)
 
+    l = zeros(2)
     c_ε = CPEDMF.c_ε(param_set)
+    c_δ = CPEDMF.c_δ(param_set)
+    c_t = CPEDMF.c_t(param_set)
+    c_λ = CPEDMF.c_λ(param_set)
+    w_min = CPEDMF.w_min(param_set)
     c_div = ICP.entrainment_massflux_div_factor(param_set)
 
+    # should be: c_δ = sign(condensate(ts_en) + condensate(ts_up[i])) * entr.c_δ
     if (εδ_model.q_liq_up + εδ_model.q_liq_en) == 0.0
         c_δ = 0.0
     end
@@ -45,7 +35,7 @@ function entr_detr(param_set, w_min, β, c_δ, c_μ, c_μ0, χ_upd, c_λ, c_εt,
 
     Δb = (εδ_model.b_up - εδ_model.b_en)
 
-    D_ε, D_δ, M_δ, M_ε = nondimensional_exchange_functions(param_set, c_δ, c_μ, c_μ0, β, χ_upd, Δw, Δb, εδ_model)
+    D_ε, D_δ, M_δ, M_ε = nondimensional_exchange_functions(param_set, Δw, Δb, εδ_model)
 
     l_1 = c_λ * abs(Δb / sqrt(εδ_model.tke + 1e-8))
     l_2 = abs(Δb / Δw)
@@ -55,10 +45,8 @@ function entr_detr(param_set, w_min, β, c_δ, c_μ, c_μ0, χ_upd, c_λ, c_εt,
     MdMdz = max(εδ_model.dMdz / max(εδ_model.M, 1e-12), 0.0)
     MdMdz = max(-εδ_model.dMdz / max(εδ_model.M, 1e-12), 0.0)
 
-
-
     # turbulent entrainment
-    K_ε = εδ_model.a_up * c_εt * sqrt(max(εδ_model.tke, 0.0)) * εδ_model.R_up
+    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
@@ -68,6 +56,5 @@ function entr_detr(param_set, w_min, β, c_δ, c_μ, c_μ0, χ_upd, c_λ, c_εt,
     ε_dyn = λ / Δw * (c_ε * D_ε + c_δ * M_ε) + MdMdz * c_div
     δ_dyn = λ / Δw * (c_ε * D_δ + c_δ * M_δ) + MdMdz * c_div
 
-
     return EntrDetr(ε_dyn, δ_dyn, ε_turb, K_ε)
 end

src/closures/nondimensional_exchange_functions.jl

@@ -1,34 +1,27 @@
 """
     nondimensional_exchange_functions(
-        param_set,
-        c_δ,
-        c_μ,
-        c_μ0,
-        β,
-        χ_up,
         Δw,
         Δb,
         εδ_model,
     )
 
 Returns the nondimensional entrainment and detrainment
 functions following Cohen et al. (JAMES, 2020), given:
- - `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(param_set, c_δ, c_μ, c_μ0, β, χ_upd, Δw, Δb, εδ_model::MoistureDeficitEntr)
+function nondimensional_exchange_functions(param_set, Δw, Δb, εδ_model)
 
-    # should be: c_δ = sign(condensate(ts_en) + condensate(ts_up[i])) * entr.c_δ
-    μ_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))
+    μ_0 = CPEDMF.μ_0(param_set)
+    μ = ICP.entrainment_sigma(param_set)
+    β = CPEDMF.β(param_set)
+    χ = CPEDMF.χ(param_set)
+
+
+    μ_ij = (χ - εδ_model.a_up / (εδ_model.a_up + εδ_model.a_en)) * Δb / Δw
+    D_ε = 1.0 / (1.0 + exp(-μ / μ_0 * μ_ij))
+    D_δ = 1.0 / (1.0 + exp(μ / μ_0 * μ_ij))
 
     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 / β)

src/types.jl

@@ -739,13 +739,6 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
     extrapolate_buoyancy::Bool
     surface_area::Float64
     max_area::Float64
-    updraft_mixing_frac::Float64
-    entrainment_sigma::Float64
-    entrainment_smin_tke_coeff::Float64
-    entrainment_scale::Float64
-    detrainment_factor::Float64
-    sorting_power::Float64
-    turbulent_entrainment_factor::Float64
     aspect_ratio::Float64
     tke_ed_coeff::Float64
     tke_diss_coeff::Float64
@@ -843,14 +836,6 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
         surface_area = namelist["turbulence"]["EDMF_PrognosticTKE"]["surface_area"]
         max_area = namelist["turbulence"]["EDMF_PrognosticTKE"]["max_area"]
         # entrainment parameters
-        entrainment_factor = namelist["turbulence"]["EDMF_PrognosticTKE"]["entrainment_factor"]
-        updraft_mixing_frac = namelist["turbulence"]["EDMF_PrognosticTKE"]["updraft_mixing_frac"]
-        entrainment_sigma = namelist["turbulence"]["EDMF_PrognosticTKE"]["entrainment_sigma"]
-        entrainment_smin_tke_coeff = namelist["turbulence"]["EDMF_PrognosticTKE"]["entrainment_smin_tke_coeff"]
-        entrainment_scale = namelist["turbulence"]["EDMF_PrognosticTKE"]["entrainment_scale"]
-        detrainment_factor = namelist["turbulence"]["EDMF_PrognosticTKE"]["detrainment_factor"]
-        sorting_power = namelist["turbulence"]["EDMF_PrognosticTKE"]["sorting_power"]
-        turbulent_entrainment_factor = namelist["turbulence"]["EDMF_PrognosticTKE"]["turbulent_entrainment_factor"]
         # pressure parameters
         aspect_ratio = namelist["turbulence"]["EDMF_PrognosticTKE"]["aspect_ratio"]
 
@@ -960,13 +945,6 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
             extrapolate_buoyancy,
             surface_area,
             max_area,
-            updraft_mixing_frac,
-            entrainment_sigma,
-            entrainment_smin_tke_coeff,
-            entrainment_scale,
-            detrainment_factor,
-            sorting_power,
-            turbulent_entrainment_factor,
             aspect_ratio,
             tke_ed_coeff,
             tke_diss_coeff,