diff --git a/integration_tests/utils/generate_namelist.jl b/integration_tests/utils/generate_namelist.jl index 4acd536168..3e7ef2fa87 100644 --- a/integration_tests/utils/generate_namelist.jl +++ b/integration_tests/utils/generate_namelist.jl @@ -57,7 +57,7 @@ function default_namelist(case_name::String) namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["max_area"] = 0.9 namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_factor"] = 0.13 namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["detrainment_factor"] = 0.51 - namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_massflux_div_factor"] = 0.4 + namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_massflux_div_factor"] = 0.0 namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["turbulent_entrainment_factor"] = 0.015 namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_ed_mf_sigma"] = 50.0 namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_smin_tke_coeff"] = 0.3 @@ -345,7 +345,7 @@ function DryBubble(namelist_defaults) namelist["meta"]["simname"] = "DryBubble" namelist["meta"]["casename"] = "DryBubble" - namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment"] = "moisture_deficit_div" + namelist["turbulence"]["EDMF_PrognosticTKE"]["entrainment_massflux_div_factor"] = 0.4 return namelist end diff --git a/src/TurbulenceConvection.jl b/src/TurbulenceConvection.jl index d3bb3ba8dd..7d498a6095 100644 --- a/src/TurbulenceConvection.jl +++ b/src/TurbulenceConvection.jl @@ -73,5 +73,7 @@ include("forcing_functions.jl") include("Surface.jl") include("surface_functions.jl") include("closures/perturbation_pressure.jl") +include("closures/entr_detr.jl") +include("closures/nondimensional_exchange_functions.jl") end diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl index 58bd0159af..7a5db9a252 100755 --- a/src/Turbulence_PrognosticTKE.jl +++ b/src/Turbulence_PrognosticTKE.jl @@ -146,8 +146,6 @@ function io(self::EDMF_PrognosticTKE, Stats::NetCDFIO_Stats, TS::TimeStepping) write_profile(Stats, "horiz_K_eddy", mean_horiz_K_eddy[cinterior]) write_profile(Stats, "entrainment_sc", mean_entr_sc[cinterior]) write_profile(Stats, "detrainment_sc", mean_detr_sc[cinterior]) - write_profile(Stats, "sorting_function", mean_sorting_function[cinterior]) - write_profile(Stats, "b_mix", mean_b_mix[cinterior]) write_profile(Stats, "nh_pressure", mean_nh_pressure[cinterior]) write_profile(Stats, "nh_pressure_adv", mean_nh_pressure_adv[cinterior]) write_profile(Stats, "nh_pressure_drag", mean_nh_pressure_drag[cinterior]) @@ -902,8 +900,7 @@ function get_env_covar_from_GMV( end function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::CasesBase) - ret = entr_struct() - input = entr_in_struct() + εδ_model = entr_detr_model() sa = eos_struct() quadrature_order = 3 grid = get_grid(self) @@ -911,94 +908,87 @@ function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariabl upd_cloud_diagnostics(self.UpdVar, reference_state(self)) - input.wstar = self.wstar - - input.dz = get_grid(self).dz - input.zbl = compute_zbl_qt_grad(GMV) - input.sort_pow = self.sorting_power - input.c_det = self.detrainment_factor - input.c_div = self.entrainment_Mdiv_factor - input.c_mu = self.entrainment_sigma - input.c_mu0 = self.entrainment_scale - input.c_ed_mf = self.entrainment_ed_mf_sigma - input.chi_upd = self.updraft_mixing_frac - input.tke_coef = self.entrainment_smin_tke_coeff - # 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) - tau = get_mixing_tau(self.zi, self.wstar) - - ae = 1 .- self.UpdVar.Area.bulkvalues l = zeros(2) @inbounds for k in real_center_indicies(grid) @inbounds for i in xrange(self.n_updrafts) - alen = max(length(argwhere(self.UpdVar.Area.values[i, cinterior])), 1) - avals = self.UpdVar.Area.values[i, cinterior] - input.zi = self.UpdVar.cloud_base[i] # entrainment - input.buoy_ed_flux = self.EnvVar.TKE.buoy[k] - if self.UpdVar.Area.values[i, k] > 0.0 - input.b_upd = self.UpdVar.B.values[i, k] - input.w_upd = interpf2c(self.UpdVar.W.values, grid, k, i) - input.z = grid.z_half[k] - input.a_upd = self.UpdVar.Area.values[i, k] - input.a_env = 1.0 - self.UpdVar.Area.bulkvalues[k] - input.tke = self.EnvVar.TKE.values[k] - input.ql_env = self.EnvVar.QL.values[k] - input.b_env = self.EnvVar.B.values[k] - input.w_env = interpf2c(self.EnvVar.W.values, grid, k) - input.ql_up = self.UpdVar.QL.values[i, k] - input.ql_env = self.EnvVar.QL.values[k] - input.RH_upd = self.UpdVar.RH.values[i, k] - input.RH_env = self.EnvVar.RH.values[k] - + w_upd = interpf2c(self.UpdVar.W.values, grid, k, i) + a_upd = interpf2c(self.UpdVar.Area.values, grid, k, i) + if w_upd * a_upd > 0.0 # compute dMdz at half levels gmv_w_k = interpf2c(GMV.W.values, grid, k) gmv_w_km = interpf2c(GMV.W.values, grid, k - 1) upd_w_km = interpf2c(self.UpdVar.W.values, grid, k - 1, i) - input.M = input.a_upd * (input.w_upd - gmv_w_k) + M = a_upd * (w_upd - gmv_w_k) Mm = self.UpdVar.Area.values[i, k - 1] * (upd_w_km - gmv_w_km) - input.dMdz = (input.M - Mm) * grid.dzi - input.tke = self.EnvVar.TKE.values[k] - - ret = self.entr_detr_fp(param_set, input) - self.entr_sc[i, k] = ret.entr_sc - self.detr_sc[i, k] = ret.detr_sc - self.sorting_function[i, k] = ret.sorting_function - self.b_mix[i, k] = ret.b_mix - + 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( + param_set, + w_min, + self.sorting_power, + self.detrainment_factor, + self.entrainment_Mdiv_factor, + self.entrainment_sigma, + self.entrainment_scale, + self.entrainment_ed_mf_sigma, + self.updraft_mixing_frac, + self.entrainment_smin_tke_coeff, + self.turbulent_entrainment_factor, + εδ_model, + ) else self.entr_sc[i, k] = 0.0 self.detr_sc[i, k] = 0.0 - self.sorting_function[i, k] = 0.0 - self.b_mix[i, k] = self.EnvVar.B.values[k] - end - - # horizontal eddy diffusivities - if self.UpdVar.Area.values[i, k] > 0.0 - self.horiz_K_eddy[i, k] = - self.UpdVar.Area.values[i, k] * - self.turbulent_entrainment_factor * - sqrt(fmax(self.EnvVar.TKE.values[k], 0.0)) * - self.pressure_plume_spacing[i] - else + self.frac_turb_entr[i, k] = 0.0 self.horiz_K_eddy[i, k] = 0.0 end - # turbulent entrainment - wu_half = interp2pt(self.UpdVar.W.values[i, k], self.UpdVar.W.values[i, k - 1]) - if self.UpdVar.Area.values[i, k] * wu_half > 0.0 - self.frac_turb_entr[i, k] = - (2.0 / self.pressure_plume_spacing[i]^2.0) * self.horiz_K_eddy[i, k] / wu_half / - self.UpdVar.Area.values[i, k] - else - self.frac_turb_entr[i, k] = 0.0 - end + # # horizontal eddy diffusivities + # if self.UpdVar.Area.values[i, k] > 0.0 + # self.horiz_K_eddy[i, k] = + # self.UpdVar.Area.values[i, k] * + # self.turbulent_entrainment_factor * + # sqrt(fmax(self.EnvVar.TKE.values[k], 0.0)) * + # self.pressure_plume_spacing[i] + # else + # self.horiz_K_eddy[i, k] = 0.0 + # end + + # # turbulent entrainment + # wu_half = interp2pt(self.UpdVar.W.values[i, k], self.UpdVar.W.values[i, k - 1]) + # if self.UpdVar.Area.values[i, k] * wu_half > 0.0 + # self.frac_turb_entr[i, k] = + # (2.0 / self.pressure_plume_spacing[i]^2.0) * self.horiz_K_eddy[i, k] / wu_half / + # self.UpdVar.Area.values[i, k] + # else + # self.frac_turb_entr[i, k] = 0.0 + # end # pressure a_kfull = interpc2f(self.UpdVar.Area.values, grid, k, i) @@ -1033,16 +1023,6 @@ function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariabl return end -function compute_zbl_qt_grad(GMV::GridMeanVariables) - grid = get_grid(GMV) - # computes inversion height as z with max gradient of qt - z∇q_tot = map(real_center_indicies(grid)) do k - (grid.z_half[k], abs(∇c2f(GMV.QT.values, grid, k))) - end - k_star = argmax(last.(z∇q_tot)) - return first.(z∇q_tot)[k_star] -end - function compute_pressure_plume_spacing(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::CasesBase) @inbounds for i in xrange(self.n_updrafts) @@ -1785,7 +1765,6 @@ function compute_covariance(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, Ca self.UpdVar.H, self.UpdVar.QT, ) - cleanup_covariance(self, GMV) return end @@ -1872,49 +1851,6 @@ function initialize_covariance(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, return end -function cleanup_covariance(self::EDMF_PrognosticTKE, GMV::GridMeanVariables) - - tmp_eps = 1e-18 - grid = get_grid(self) - - @inbounds for k in real_center_indicies(grid) - if GMV.TKE.values[k] < tmp_eps - GMV.TKE.values[k] = 0.0 - end - if GMV.Hvar.values[k] < tmp_eps - GMV.Hvar.values[k] = 0.0 - end - if GMV.QTvar.values[k] < tmp_eps - GMV.QTvar.values[k] = 0.0 - end - if fabs(GMV.HQTcov.values[k]) < tmp_eps - GMV.HQTcov.values[k] = 0.0 - end - if fabs(GMV.W_third_m.values[k]) < tmp_eps - GMV.W_third_m.values[k] = 0.0 - end - if fabs(GMV.H_third_m.values[k]) < tmp_eps - GMV.H_third_m.values[k] = 0.0 - end - if fabs(GMV.QT_third_m.values[k]) < tmp_eps - GMV.QT_third_m.values[k] = 0.0 - end - if self.EnvVar.Hvar.values[k] < tmp_eps - self.EnvVar.Hvar.values[k] = 0.0 - end - if self.EnvVar.TKE.values[k] < tmp_eps - self.EnvVar.TKE.values[k] = 0.0 - end - if self.EnvVar.QTvar.values[k] < tmp_eps - self.EnvVar.QTvar.values[k] = 0.0 - end - if fabs(self.EnvVar.HQTcov.values[k]) < tmp_eps - self.EnvVar.HQTcov.values[k] = 0.0 - end - end -end - - function compute_covariance_shear( self::EDMF_PrognosticTKE, GMV::GridMeanVariables, diff --git a/src/closures/entr_detr.jl b/src/closures/entr_detr.jl new file mode 100644 index 0000000000..7b83648e78 --- /dev/null +++ b/src/closures/entr_detr.jl @@ -0,0 +1,85 @@ +#### Entrainment-Detrainment kernels +""" + entr_detr( + param_set, + w_min, + β, + c_δ, + c_div, + c_μ, + c_μ0, + c_ed_mf, + χ_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_div`: divergence factor for bubble case (zero otherwise) + - `c_μ`: logisitc function scale + - `c_μ0`: logisitc function timescale + - `c_ed_mf`: ED-MF coeeficeicent for relating BL tke and velocity scales + - `χ_upd`: updraft mixing fraction + - `c_λ`: + - `εδ_model`: updraft buoyancy + -- Where: + - `εδ_model.b_upd`: updraft buoyancy + - `εδ_model.w_upd`: updraft vertical velocity + - `εδ_model.a_upd`: updraft area fraction + - `εδ_model.a_env`: environment buoyancy + - `εδ_model.b_env`: environment vertical velocity + - `εδ_model.w_env`: environment area fraction + - `εδ_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 +""" +function entr_detr(param_set, w_min, β, c_δ, c_div, c_μ, c_μ0, c_ed_mf, χ_upd, c_λ, ε_t, εδ_model) + + l = zeros(2) + c_ε = CPEDMF.c_ε(param_set) + + if (εδ_model.ql_upd + εδ_model.ql_env) == 0.0 + c_δ = 0.0 + end + + Δw = εδ_model.w_upd - εδ_model.w_env + if Δw < 0.0 + Δw -= w_min + else + Δw += w_min + end + + Δb = (εδ_model.b_upd - εδ_model.b_env) + + D_ε, D_δ, M_δ, M_ε = nondimensional_exchange_functions(c_ε, c_δ, c_μ, c_μ0, β, χ_upd, Δw, Δb, εδ_model) + + l[1] = c_λ * fabs(Δb / sqrt(εδ_model.tke + 1e-8)) + l[2] = fabs(Δb / Δw) + λ = lamb_smooth_minimum(l, 0.1, 0.0005) + + MdMdz = fmax(εδ_model.dMdz / fmax(εδ_model.M, 1e-12), 0.0) + MdMdz = fmax(-εδ_model.dMdz / fmax(εδ_model.M, 1e-12), 0.0) + + + + # turbulent entrainment + k_ε = εδ_model.a_upd * ε_t * sqrt(fmax(εδ_model.tke, 0.0)) * εδ_model.R_up + ε_turb = (2.0 / εδ_model.R_up^2.0) * k_ε / (εδ_model.w_upd * εδ_model.a_upd) + + ε_dyn = λ / Δw * (c_ε * D_ε + c_δ * M_ε) + MdMdz * c_div + δ_dyn = λ / Δw * (c_ε * D_δ + c_δ * M_δ) + MdMdz * c_div + + return ε_dyn, δ_dyn, ε_turb, k_ε +end diff --git a/src/closures/nondimensional_exchange_functions.jl b/src/closures/nondimensional_exchange_functions.jl new file mode 100644 index 0000000000..76210ca1a0 --- /dev/null +++ b/src/closures/nondimensional_exchange_functions.jl @@ -0,0 +1,42 @@ +""" + nondimensional_exchange_functions( + c_ε, + c_δ, + c_μ, + c_μ0, + β, + χ_up, + Δ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) + - `β`, sorting power + - `χ_up`, updraft mixing fraction + - `Δw`, updraft - environment vertical velocity differnce + - `Δb`, updraft - environment buoynacy differnce + - `εδ_model`, updraft mixing fraction + -- Where + - `εδ_model.RH_upd`, updraft relative humidity + - `εδ_model.RH_env`, environment relative humidity + - `εδ_model.a_upd`, updraft area fraction + - `εδ_model.a_env`, environment area fraction +""" +function nondimensional_exchange_functions(c_ε, c_δ, c_μ, c_μ0, β, χ_upd, Δw, Δb, εδ_model) + + # 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 + D_ε = 1.0 / (1.0 + exp(-c_μ / c_μ0 * μ_ij)) + D_δ = 1.0 / (1.0 + exp(c_μ / c_μ0 * μ_ij)) + + M_δ = (fmax((εδ_model.RH_upd / 100.0)^β - (εδ_model.RH_env / 100.0)^β, 0.0))^(1.0 / β) + M_ε = (fmax((εδ_model.RH_env / 100.0)^β - (εδ_model.RH_upd / 100.0)^β, 0.0))^(1.0 / β) + + return D_ε, D_δ, M_δ, M_ε +end; diff --git a/src/turbulence_functions.jl b/src/turbulence_functions.jl index 27754264c4..845eeef276 100644 --- a/src/turbulence_functions.jl +++ b/src/turbulence_functions.jl @@ -1,147 +1,3 @@ -# Entrainment Rates - -function entr_detr_env_moisture_deficit_b_ED_MF(param_set::APS, entr_in::entr_struct) - _ret = entr_struct() - - moisture_deficit_d = - (fmax( - (entr_in.RH_upd / 100.0)^entr_in.sort_pow - (entr_in.RH_env / 100.0)^entr_in.sort_pow, - 0.0, - ))^(1.0 / entr_in.sort_pow) - moisture_deficit_e = - (fmax( - (entr_in.RH_env / 100.0)^entr_in.sort_pow - (entr_in.RH_upd / 100.0)^entr_in.sort_pow, - 0.0, - ))^(1.0 / entr_in.sort_pow) - _ret.sorting_function = moisture_deficit_e - c_det = entr_in.c_det - if (entr_in.ql_up + entr_in.ql_env) ≈ 0.0 - c_det = 0.0 - end - - dw = entr_in.w_upd - entr_in.w_env - if dw < 0.0 - dw -= 0.001 - else - dw += 0.001 - end - - db = (entr_in.b_upd - entr_in.b_env) - mu = entr_in.c_mu / entr_in.c_mu0 - - inv_timescale = fabs(db / dw) - logistic_e = 1.0 / (1.0 + exp(-mu * db / dw * (entr_in.chi_upd - entr_in.a_upd / (entr_in.a_upd + entr_in.a_env)))) - logistic_d = 1.0 / (1.0 + exp(mu * db / dw * (entr_in.chi_upd - entr_in.a_upd / (entr_in.a_upd + entr_in.a_env)))) - - #Logistic of buoyancy fluxes - inv_timescale = fabs(db / dw) - ed_mf_ratio = - fabs(entr_in.buoy_ed_flux) / - (fabs(entr_in.a_upd * entr_in.a_env * (entr_in.w_upd - entr_in.w_env) * (entr_in.b_upd - entr_in.b_env)) + 1e-8) - logistic_e *= (1.0 / (1.0 + exp(entr_in.c_ed_mf * (ed_mf_ratio - 1.0)))) - c_ε = CPEDMF.c_ε(param_set) - _ret.entr_sc = inv_timescale / dw * (c_ε * logistic_e + c_det * moisture_deficit_e) - _ret.detr_sc = inv_timescale / dw * (c_ε * logistic_d + c_det * moisture_deficit_d) - - return _ret -end - -function entr_detr_env_moisture_deficit(param_set::APS, entr_in::entr_in_struct) - _ret = entr_struct() - l = zeros(2) - - moisture_deficit_d = - (fmax( - (entr_in.RH_upd / 100.0)^entr_in.sort_pow - (entr_in.RH_env / 100.0)^entr_in.sort_pow, - 0.0, - ))^(1.0 / entr_in.sort_pow) - moisture_deficit_e = - (fmax( - (entr_in.RH_env / 100.0)^entr_in.sort_pow - (entr_in.RH_upd / 100.0)^entr_in.sort_pow, - 0.0, - ))^(1.0 / entr_in.sort_pow) - _ret.sorting_function = moisture_deficit_e - c_det = entr_in.c_det - if (entr_in.ql_up + entr_in.ql_env) == 0.0 - c_det = 0.0 - end - - dw = entr_in.w_upd - entr_in.w_env - if dw < 0.0 - dw -= 0.001 - else - dw += 0.001 - end - - db = (entr_in.b_upd - entr_in.b_env) - mu = entr_in.c_mu / entr_in.c_mu0 - - inv_timescale = fabs(db / dw) - logistic_e = 1.0 / (1.0 + exp(-mu * db / dw * (entr_in.chi_upd - entr_in.a_upd / (entr_in.a_upd + entr_in.a_env)))) - logistic_d = 1.0 / (1.0 + exp(mu * db / dw * (entr_in.chi_upd - entr_in.a_upd / (entr_in.a_upd + entr_in.a_env)))) - - #smooth min - l[1] = entr_in.tke_coef * fabs(db / sqrt(entr_in.tke + 1e-8)) - l[2] = fabs(db / dw) - inv_timescale = lamb_smooth_minimum(l, 0.1, 0.0005) - c_ε = CPEDMF.c_ε(param_set) - _ret.entr_sc = inv_timescale / dw * (c_ε * logistic_e + c_det * moisture_deficit_e) - _ret.detr_sc = inv_timescale / dw * (c_ε * logistic_d + c_det * moisture_deficit_d) - - return _ret -end - -function entr_detr_env_moisture_deficit_div(param_set::APS, entr_in::entr_struct) - _ret = entr_struct() - l = zeros(2) - - moisture_deficit_d = - (fmax( - (entr_in.RH_upd / 100.0)^entr_in.sort_pow - (entr_in.RH_env / 100.0)^entr_in.sort_pow, - 0.0, - ))^(1.0 / entr_in.sort_pow) - moisture_deficit_e = - (fmax( - (entr_in.RH_env / 100.0)^entr_in.sort_pow - (entr_in.RH_upd / 100.0)^entr_in.sort_pow, - 0.0, - ))^(1.0 / entr_in.sort_pow) - _ret.sorting_function = moisture_deficit_e - c_det = entr_in.c_det - if (entr_in.ql_up + entr_in.ql_env) == 0.0 - c_det = 0.0 - end - - dw = entr_in.w_upd - entr_in.w_env - if dw < 0.0 - dw -= 0.001 - else - dw += 0.001 - end - - db = (entr_in.b_upd - entr_in.b_env) - mu = entr_in.c_mu / entr_in.c_mu0 - - inv_timescale = fabs(db / dw) - logistic_e = 1.0 / (1.0 + exp(-mu * db / dw * (entr_in.chi_upd - entr_in.a_upd / (entr_in.a_upd + entr_in.a_env)))) - logistic_d = 1.0 / (1.0 + exp(mu * db / dw * (entr_in.chi_upd - entr_in.a_upd / (entr_in.a_upd + entr_in.a_env)))) - - entr_MdMdz = fmax(entr_in.dMdz / fmax(entr_in.M, 1e-12), 0.0) - detr_MdMdz = fmax(-entr_in.dMdz / fmax(entr_in.M, 1e-12), 0.0) - - - #smooth min - l[1] = entr_in.tke_coef * fabs(db / sqrt(entr_in.tke + 1e-8)) - l[2] = fabs(db / dw) - inv_timescale = lamb_smooth_minimum(l, 0.1, 0.0005) - - c_ε = CPEDMF.c_ε(param_set) - _ret.entr_sc = inv_timescale / dw * (c_ε * logistic_e + c_det * moisture_deficit_e) + entr_MdMdz * entr_in.c_div - _ret.detr_sc = inv_timescale / dw * (c_ε * logistic_d + c_det * moisture_deficit_d) + detr_MdMdz * entr_in.c_div - - - return _ret -end - # convective velocity scale function get_wstar(bflux, zi) return cbrt(fmax(bflux * zi, 0.0)) diff --git a/src/types.jl b/src/types.jl index 29ee0adf17..3cabf0db1d 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,10 +1,3 @@ -Base.@kwdef mutable struct entr_struct - entr_sc::Float64 = 0 - detr_sc::Float64 = 0 - sorting_function::Float64 = 0 - b_mix::Float64 = 0 -end - Base.@kwdef mutable struct eos_struct T::Float64 = 0 ql::Float64 = 0 @@ -27,41 +20,23 @@ Base.@kwdef mutable struct mph_struct qr_src::Float64 = 0 end -Base.@kwdef mutable struct entr_in_struct - c_div::Float64 = 0 - M::Float64 = 0 - dMdz::Float64 = 0 - zi::Float64 = 0 - wstar::Float64 = 0 - z::Float64 = 0 - sort_pow::Float64 = 0 - c_det::Float64 = 0 - chi_upd::Float64 = 0 - tke_coef::Float64 = 0 - c_mu::Float64 = 0 - c_ed_mf::Float64 = 0 - c_mu0::Float64 = 0 - dz::Float64 = 0 - w_upd::Float64 = 0 +Base.@kwdef mutable struct entr_detr_model b_upd::Float64 = 0 - dt::Float64 = 0 b_env::Float64 = 0 + w_upd::Float64 = 0 + w_env::Float64 = 0 a_upd::Float64 = 0 a_env::Float64 = 0 - tke::Float64 = 0 + ql_upd::Float64 = 0 + ql_env::Float64 = 0 RH_upd::Float64 = 0 RH_env::Float64 = 0 - qt_up::Float64 = 0 - ql_up::Float64 = 0 - T_env::Float64 = 0 - qt_env::Float64 = 0 - ql_env::Float64 = 0 - w_env::Float64 = 0 - nh_pressure::Float64 = 0 - L::Float64 = 0 - zbl::Float64 = 0 - poisson::Float64 = 0 - buoy_ed_flux::Float64 = 0 + M::Float64 = 0 + dMdz::Float64 = 0 + tke::Float64 = 0 + n_up::Float64 = 0 + ρ::Float64 = 0 + R_up::Float64 = 0 end struct RainVariable{T} @@ -723,7 +698,6 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2} zi::Float64 n_updrafts::Int use_const_plume_spacing::Bool - entr_detr_fp::Function drag_sign::Int asp_label extrapolate_buoyancy::Bool @@ -813,24 +787,24 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2} use_const_plume_spacing = parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "use_constant_plume_spacing"; default = false) - entr_detr_fp_str = parse_namelist( - namelist, - "turbulence", - "EDMF_PrognosticTKE", - "entrainment"; - default = "entr_detr_env_moisture_deficit", - ) - entr_detr_fp = if entr_detr_fp_str == "moisture_deficit" - entr_detr_env_moisture_deficit - elseif entr_detr_fp_str == "moisture_deficit_b_ED_MF" - entr_detr_env_moisture_deficit_b_ED_MF - elseif entr_detr_fp_str == "moisture_deficit_div" - entr_detr_env_moisture_deficit_div - elseif entr_detr_fp_str == "none" - entr_detr_none - else - error("Turbulence--EDMF_PrognosticTKE: Entrainment rate namelist option is not recognized") - end + # entr_detr_fp_str = parse_namelist( + # namelist, + # "turbulence", + # "EDMF_PrognosticTKE", + # "entrainment"; + # default = "entr_detr_env_moisture_deficit", + # ) + # entr_detr_fp = if entr_detr_fp_str == "moisture_deficit" + # entr_detr_env_moisture_deficit + # elseif entr_detr_fp_str == "moisture_deficit_b_ED_MF" + # entr_detr_env_moisture_deficit_b_ED_MF + # elseif entr_detr_fp_str == "moisture_deficit_div" + # entr_detr_env_moisture_deficit_div + # elseif entr_detr_fp_str == "none" + # entr_detr_none + # else + # error("Turbulence--EDMF_PrognosticTKE: Entrainment rate namelist option is not recognized") + # end pressure_func_drag_str = parse_namelist( namelist, @@ -863,13 +837,7 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2} surface_area = namelist["turbulence"]["EDMF_PrognosticTKE"]["surface_area"] max_area = namelist["turbulence"]["EDMF_PrognosticTKE"]["max_area"] entrainment_factor = namelist["turbulence"]["EDMF_PrognosticTKE"]["entrainment_factor"] - entrainment_Mdiv_factor = parse_namelist( - namelist, - "turbulence", - "EDMF_PrognosticTKE", - "entrainment_massflux_div_factor"; - default = 0.0, - ) + entrainment_Mdiv_factor = namelist["turbulence"]["EDMF_PrognosticTKE"]["entrainment_massflux_div_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"] @@ -1020,7 +988,6 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2} zi, n_updrafts, use_const_plume_spacing, - entr_detr_fp, drag_sign, asp_label, extrapolate_buoyancy,