From 016d165959e3171a5f3de3ce42a3cbfd5d8c8793 Mon Sep 17 00:00:00 2001 From: Zhaoyi Shen <11598433+szy21@users.noreply.github.com> Date: Thu, 19 Oct 2023 13:47:52 -0700 Subject: [PATCH] add prognostic tke equation for prognostic edmf --- .../prognostic_edmfx_bomex_box.yml | 76 +++++++++---------- .../prognostic_edmfx_dycoms_rf01_box.yml | 6 +- .../prognostic_edmfx_gabls_box.yml | 2 +- .../prognostic_edmfx_rico_box.yml | 2 +- .../prognostic_edmfx_trmm_box.yml | 2 +- src/cache/precomputed_quantities.jl | 1 + .../prognostic_edmf_precomputed_quantities.jl | 31 ++++++-- src/prognostic_equations/edmfx_closures.jl | 20 +++-- src/prognostic_equations/edmfx_entr_detr.jl | 9 ++- src/prognostic_equations/edmfx_sgs_flux.jl | 4 +- src/prognostic_equations/edmfx_tke.jl | 13 ++-- 11 files changed, 96 insertions(+), 70 deletions(-) diff --git a/config/model_configs/prognostic_edmfx_bomex_box.yml b/config/model_configs/prognostic_edmfx_bomex_box.yml index bd8a46b0967..067cf3d7ae1 100644 --- a/config/model_configs/prognostic_edmfx_bomex_box.yml +++ b/config/model_configs/prognostic_edmfx_bomex_box.yml @@ -11,7 +11,7 @@ edmfx_detr_model: "BOverW" edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true edmfx_nh_pressure: true -prognostic_tke: false +prognostic_tke: true moist: "equil" config: "box" hyperdiff: "true" @@ -24,78 +24,78 @@ y_elem: 2 z_elem: 60 z_stretch: false perturb_initstate: false -dt: "5secs" +dt: "2secs" t_end: "6hours" -dt_save_to_disk: "10mins" +dt_save_to_disk: "2mins" toml: [toml/prognostic_edmfx_box.toml] diagnostics: - short_name: ts - period: 10mins + period: 1hours - short_name: ta - period: 10mins + period: 1hours - short_name: thetaa - period: 10mins + period: 1hours - short_name: ha - period: 10mins + period: 1hours - short_name: pfull - period: 10mins + period: 1hours - short_name: rhoa - period: 10mins + period: 1hours - short_name: ua - period: 10mins + period: 1hours - short_name: va - period: 10mins + period: 1hours - short_name: wa - period: 10mins + period: 1hours - short_name: hur - period: 10mins + period: 1hours - short_name: hus - period: 10mins + period: 1hours - short_name: clw - period: 10mins + period: 1hours - short_name: cli - period: 10mins + period: 1hours - short_name: hussfc - period: 10mins + period: 1hours - short_name: evspsbl - period: 10mins + period: 1hours - short_name: arup - period: 10mins + period: 1hours - short_name: waup - period: 10mins + period: 1hours - short_name: taup - period: 10mins + period: 1hours - short_name: thetaaup - period: 10mins + period: 1hours - short_name: haup - period: 10mins + period: 1hours - short_name: husup - period: 10mins + period: 1hours - short_name: hurup - period: 10mins + period: 1hours - short_name: clwup - period: 10mins + period: 1hours - short_name: cliup - period: 10mins + period: 1hours - short_name: rhoaen - period: 10mins + period: 1hours - short_name: waen - period: 10mins + period: 1hours - short_name: taen - period: 10mins + period: 1hours - short_name: thetaaen - period: 10mins + period: 1hours - short_name: haen - period: 10mins + period: 1hours - short_name: husen - period: 10mins + period: 1hours - short_name: huren - period: 10mins + period: 1hours - short_name: clwen - period: 10mins + period: 1hours - short_name: clien - period: 10mins + period: 1hours - short_name: tke - period: 10mins + period: 1hours - short_name: lmix - period: 10mins + period: 1hours diff --git a/config/model_configs/prognostic_edmfx_dycoms_rf01_box.yml b/config/model_configs/prognostic_edmfx_dycoms_rf01_box.yml index f308d10d55c..1e293b90019 100644 --- a/config/model_configs/prognostic_edmfx_dycoms_rf01_box.yml +++ b/config/model_configs/prognostic_edmfx_dycoms_rf01_box.yml @@ -11,7 +11,7 @@ edmfx_detr_model: "BOverW" edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true edmfx_nh_pressure: true -prognostic_tke: false +prognostic_tke: true moist: equil config: box hyperdiff: "true" @@ -23,7 +23,7 @@ x_elem: 2 y_elem: 2 z_elem: 30 z_stretch: false -dt: 5secs +dt: 2secs t_end: 4hours -dt_save_to_disk: 10mins +dt_save_to_disk: 2mins toml: [toml/prognostic_edmfx_box.toml] diff --git a/config/model_configs/prognostic_edmfx_gabls_box.yml b/config/model_configs/prognostic_edmfx_gabls_box.yml index ef407a6c332..407e8f0474f 100644 --- a/config/model_configs/prognostic_edmfx_gabls_box.yml +++ b/config/model_configs/prognostic_edmfx_gabls_box.yml @@ -9,7 +9,7 @@ edmfx_detr_model: "BOverW" edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true edmfx_nh_pressure: true -prognostic_tke: false +prognostic_tke: true moist: "equil" config: "box" hyperdiff: "true" diff --git a/config/model_configs/prognostic_edmfx_rico_box.yml b/config/model_configs/prognostic_edmfx_rico_box.yml index e437c91236f..1dd10451f5c 100644 --- a/config/model_configs/prognostic_edmfx_rico_box.yml +++ b/config/model_configs/prognostic_edmfx_rico_box.yml @@ -10,7 +10,7 @@ edmfx_detr_model: "BOverW" edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true edmfx_nh_pressure: true -prognostic_tke: false +prognostic_tke: true moist: "equil" precip_model: "0M" config: "box" diff --git a/config/model_configs/prognostic_edmfx_trmm_box.yml b/config/model_configs/prognostic_edmfx_trmm_box.yml index 1cbbb0ba9dd..f6eac8fbd7f 100644 --- a/config/model_configs/prognostic_edmfx_trmm_box.yml +++ b/config/model_configs/prognostic_edmfx_trmm_box.yml @@ -9,7 +9,7 @@ edmfx_detr_model: "BOverW" edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true edmfx_nh_pressure: true -prognostic_tke: false +prognostic_tke: true moist: equil apply_limiter: false precip_model: "0M" diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index bdaf9f47c89..2adbe35a445 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -78,6 +78,7 @@ function precomputed_quantities(Y, atmos) ᶜmixing_length = similar(Y.c, FT), ᶜK_u = similar(Y.c, FT), ᶜK_h = similar(Y.c, FT), + ρatke_flux = similar(Fields.level(Y.f, half), C3{FT}), ᶜuʲs = similar(Y.c, NTuple{n, C123{FT}}), ᶠu³ʲs = similar(Y.f, NTuple{n, CT3{FT}}), ᶜKʲs = similar(Y.c, NTuple{n, FT}), diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl index a0a62358d75..15f4aa12c63 100644 --- a/src/cache/prognostic_edmf_precomputed_quantities.jl +++ b/src/cache/prognostic_edmf_precomputed_quantities.jl @@ -22,9 +22,9 @@ function set_prognostic_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t) n = n_mass_flux_subdomains(turbconv_model) thermo_args = (thermo_params, energy_form, moisture_model) - (; ᶜspecific, ᶜp, ᶜΦ, ᶜh_tot, ᶜρ_ref) = p + (; ᶜspecific, ᶜp, ᶜΦ, ᶜh_tot, ᶜρ_ref, ᶜu) = p (; ᶜtke⁰, ᶜρa⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰, ᶜh_tot⁰, ᶜq_tot⁰) = p - (; ᶜmixing_length, ᶜlinear_buoygrad, ᶜstrain_rate_norm, ᶜK_u, ᶜK_h) = p + (; ᶜmixing_length, ᶜlinear_buoygrad, ᶜstrain_rate_norm, ᶜK_u, ᶜK_h, ρatke_flux) = p (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜtsʲs, ᶜρʲs, ᶜentrʲs, ᶜdetrʲs) = p (; ustar, obukhov_length, buoyancy_flux) = p.sfc_conditions @@ -226,9 +226,7 @@ function set_prognostic_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t) @. ᶜtke_exch += Y.c.sgsʲs.:($$j).ρa * ᶜdetrʲs.:($$j) / ᶜρa⁰ * ( 1 / 2 * - ( - get_physical_w(ᶜuʲs.:($$j), ᶜlg) - get_physical_w(ᶜu⁰, ᶜlg) - )^2 - ᶜtke⁰ + norm_sqr(ᶜinterp(ᶠu³⁰) - ᶜinterp(ᶠu³ʲs.:($$j))) - ᶜtke⁰ ) end @@ -239,9 +237,9 @@ function set_prognostic_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t) ᶜz, z_sfc, ᶜdz, - sfc_tke, + max(sfc_tke, 0), ᶜlinear_buoygrad, - ᶜtke⁰, + max(ᶜtke⁰, 0), obukhov_length, ᶜstrain_rate_norm, ᶜprandtl_nvec, @@ -251,8 +249,25 @@ function set_prognostic_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t) turbconv_params = CAP.turbconv_params(params) c_m = TCP.tke_ed_coeff(turbconv_params) @. ᶜK_u = c_m * ᶜmixing_length * sqrt(max(ᶜtke⁰, 0)) - # TODO: add Prantdl number @. ᶜK_h = ᶜK_u / ᶜprandtl_nvec + ρatke_flux_values = Fields.field_values(ρatke_flux) + ρ_int_values = Fields.field_values(Fields.level(ᶜρa⁰, 1)) + u_int_values = Fields.field_values(Fields.level(ᶜu, 1)) + ustar_values = Fields.field_values(ustar) + int_local_geometry_values = + Fields.field_values(Fields.level(Fields.local_geometry_field(Y.c), 1)) + sfc_local_geometry_values = Fields.field_values( + Fields.level(Fields.local_geometry_field(Y.f), half), + ) + @. ρatke_flux_values = surface_flux_tke( + turbconv_params, + ρ_int_values, + u_int_values, + ustar_values, + int_local_geometry_values, + sfc_local_geometry_values, + ) + return nothing end diff --git a/src/prognostic_equations/edmfx_closures.jl b/src/prognostic_equations/edmfx_closures.jl index 2bd3e7e9ae2..596ff60a639 100644 --- a/src/prognostic_equations/edmfx_closures.jl +++ b/src/prognostic_equations/edmfx_closures.jl @@ -194,10 +194,10 @@ function mixing_length( # kz scale (surface layer) if obukhov_length < 0.0 #unstable l_W = - vkc * (ᶜz - z_sfc) / (sqrt(sfc_tke / ustar / ustar) * c_m) * + vkc * (ᶜz - z_sfc) / (max(sqrt(sfc_tke / ustar / ustar), eps(FT)) * c_m) * min((1 - 100 * (ᶜz - z_sfc) / obukhov_length)^FT(0.2), 1 / vkc) else # neutral or stable - l_W = vkc * (ᶜz - z_sfc) / (sqrt(sfc_tke / ustar / ustar) * c_m) + l_W = vkc * (ᶜz - z_sfc) / (max(sqrt(sfc_tke / ustar / ustar), eps(FT)) * c_m) end # compute l_TKE - the production-dissipation balanced length scale @@ -235,15 +235,19 @@ function mixing_length( l_smag = c_smag * ᶜdz end + l_smag = FT(10) # add limiters - l = SA.SVector( - (l_N < eps(FT) || l_N > l_z) ? l_z : l_N, - (l_TKE < eps(FT) || l_TKE > l_z) ? l_z : l_TKE, - (l_W < eps(FT) || l_W > l_z) ? l_z : l_W, - ) + # l = SA.SVector( + # (l_N < eps(FT) || l_N > l_z) ? l_z : l_N, + # (l_TKE < eps(FT) || l_TKE > l_z) ? l_z : l_TKE, + # (l_W < eps(FT) || l_W > l_z) ? l_z : l_W, + # ) + # l_limited = lamb_smooth_minimum(l, smin_ub, smin_rm) + l = SA.SVector(l_N, l_TKE, l_W) + l_limited = max(l_smag, min(lamb_smooth_minimum(l, smin_ub, smin_rm), l_z)) # get soft minimum # TODO: limit it with l_smag - return lamb_smooth_minimum(l, smin_ub, smin_rm) + return l_limited end """ diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl index de6620abf64..41cecbdfbab 100644 --- a/src/prognostic_equations/edmfx_entr_detr.jl +++ b/src/prognostic_equations/edmfx_entr_detr.jl @@ -106,7 +106,10 @@ function entrainment( ::ConstantCoefficientEntrainment, ) where {FT} entr_coeff = CAP.entr_coeff(params) - entr = min(entr_coeff * abs(ᶜwʲ) / (ᶜz - z_sfc), 1 / dt) + min_area_limiter = FT(2e-3) * exp(-12 * ᶜaʲ) + #entr = min(entr_coeff * abs(ᶜwʲ) / (ᶜz - z_sfc), 1 / dt) + entr = + min(abs(ᶜwʲ) * (entr_coeff / (ᶜz - z_sfc) + min_area_limiter), 1 / dt) return entr end @@ -258,11 +261,13 @@ function detrainment( ) where {FT} detr_coeff = CAP.detr_coeff(params) detr_buoy_coeff = CAP.detr_buoy_coeff(params) + max_area_limiter = FT(0.1) * exp(-10 * (0.6 - ᶜaʲ)) detr = min( abs(ᶜwʲ) * ( detr_coeff + detr_buoy_coeff * abs(min(ᶜbuoyʲ - ᶜbuoy⁰, 0)) / - max(eps(FT), (ᶜwʲ - ᶜw⁰) * (ᶜwʲ - ᶜw⁰)) + max(eps(FT), (ᶜwʲ - ᶜw⁰) * (ᶜwʲ - ᶜw⁰)) + + max_area_limiter ), 1 / dt, ) diff --git a/src/prognostic_equations/edmfx_sgs_flux.jl b/src/prognostic_equations/edmfx_sgs_flux.jl index d8c574acd3d..ed0c554100b 100644 --- a/src/prognostic_equations/edmfx_sgs_flux.jl +++ b/src/prognostic_equations/edmfx_sgs_flux.jl @@ -169,7 +169,7 @@ function edmfx_sgs_diffusive_flux_tendency!( if p.atmos.edmfx_sgs_diffusive_flux # energy ᶠρaK_h = p.ᶠtemp_scalar - @. ᶠρaK_h[colidx] = ᶠinterp(ᶜρa⁰[colidx]) * ᶠinterp(ᶜK_h[colidx]) + @. ᶠρaK_h[colidx] = ᶠinterp(ᶜρa⁰[colidx]) * max(ᶠinterp(ᶜK_h[colidx]), FT(0)) ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C( top = Operators.SetValue(C3(FT(0))), @@ -192,7 +192,7 @@ function edmfx_sgs_diffusive_flux_tendency!( # momentum ᶠρaK_u = p.ᶠtemp_scalar - @. ᶠρaK_u[colidx] = ᶠinterp(ᶜρa⁰[colidx]) * ᶠinterp(ᶜK_u[colidx]) + @. ᶠρaK_u[colidx] = ᶠinterp(ᶜρa⁰[colidx]) * max(ᶠinterp(ᶜK_u[colidx]), FT(0)) ᶠstrain_rate = p.ᶠtemp_UVWxUVW compute_strain_rate_face!(ᶠstrain_rate[colidx], ᶜu⁰[colidx]) @. Yₜ.c.uₕ[colidx] -= C12( diff --git a/src/prognostic_equations/edmfx_tke.jl b/src/prognostic_equations/edmfx_tke.jl index 07841b25db1..5aa7647213e 100644 --- a/src/prognostic_equations/edmfx_tke.jl +++ b/src/prognostic_equations/edmfx_tke.jl @@ -10,7 +10,7 @@ function edmfx_tke_tendency!( p, t, colidx, - turbconv_model::DiagnosticEDMFX, + turbconv_model::Union{PrognosticEDMFX, DiagnosticEDMFX}, ) FT = Spaces.undertype(axes(Y.c)) @@ -21,12 +21,13 @@ function edmfx_tke_tendency!( (; ᶜentrʲs, ᶜdetrʲs, ᶠu³ʲs) = p (; ᶠu³⁰, ᶜstrain_rate_norm, ᶜlinear_buoygrad, ᶜtke⁰, ᶜmixing_length) = p (; ᶜK_u, ᶜK_h, ρatke_flux) = p + ᶜρa⁰ = turbconv_model isa PrognosticEDMFX ? p.ᶜρa⁰ : Y.c.ρ ᶠgradᵥ = Operators.GradientC2F() ᶠρaK_u = p.ᶠtemp_scalar if use_prognostic_tke(turbconv_model) # turbulent transport (diffusive flux) - @. ᶠρaK_u[colidx] = ᶠinterp(Y.c.ρ[colidx]) * ᶠinterp(ᶜK_u[colidx]) + @. ᶠρaK_u[colidx] = ᶠinterp(ᶜρa⁰[colidx]) * ᶠinterp(ᶜK_u[colidx]) # boundary condition for the diffusive flux ᶜdivᵥ_ρatke = Operators.DivergenceF2C( top = Operators.SetValue(C3(FT(0))), @@ -36,15 +37,15 @@ function edmfx_tke_tendency!( ᶜdivᵥ_ρatke(-(ᶠρaK_u[colidx] * ᶠgradᵥ(ᶜtke⁰[colidx]))) # shear production @. Yₜ.c.sgs⁰.ρatke[colidx] += - 2 * Y.c.ρ[colidx] * ᶜK_u[colidx] * ᶜstrain_rate_norm[colidx] + 2 * ᶜρa⁰[colidx] * ᶜK_u[colidx] * ᶜstrain_rate_norm[colidx] # buoyancy production @. Yₜ.c.sgs⁰.ρatke[colidx] -= - Y.c.ρ[colidx] * ᶜK_h[colidx] * ᶜlinear_buoygrad[colidx] + ᶜρa⁰[colidx] * ᶜK_h[colidx] * ᶜlinear_buoygrad[colidx] # entrainment and detraiment # using ᶜu⁰ and local geometry results in allocation for j in 1:n @. Yₜ.c.sgs⁰.ρatke[colidx] += - Y.c.ρ[colidx] * ( + ᶜρa⁰[colidx] * ( ᶜentrʲs.:($$j)[colidx] * 1 / 2 * norm_sqr( ᶜinterp(ᶠu³⁰[colidx]) - ᶜinterp(ᶠu³ʲs.:($$j)[colidx]), ) - ᶜdetrʲs.:($$j)[colidx] * ᶜtke⁰[colidx] @@ -53,7 +54,7 @@ function edmfx_tke_tendency!( # pressure work # dissipation @. Yₜ.c.sgs⁰.ρatke[colidx] -= - Y.c.ρ[colidx] * c_d * max(ᶜtke⁰[colidx], 0)^(FT(3) / 2) / + ᶜρa⁰[colidx] * c_d * max(ᶜtke⁰[colidx], 0)^(FT(3) / 2) / ᶜmixing_length[colidx] end