From 5703a85408b9665d53c9afe662d12e7855a41d96 Mon Sep 17 00:00:00 2001 From: costachris Date: Fri, 13 Sep 2024 15:15:41 -0700 Subject: [PATCH] Add turb entr to TKE equation, simplify entr/detr limiters --- .../diagnostic_edmf_precomputed_quantities.jl | 38 ++-------- .../prognostic_edmf_precomputed_quantities.jl | 30 ++------ src/prognostic_equations/edmfx_entr_detr.jl | 75 ++++--------------- src/prognostic_equations/edmfx_tke.jl | 10 ++- 4 files changed, 37 insertions(+), 116 deletions(-) diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl index ebc08c6c170..ddb1f98af2a 100644 --- a/src/cache/diagnostic_edmf_precomputed_quantities.jl +++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl @@ -411,36 +411,24 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( draft_area(ρaʲ_prev_level, ρʲ_prev_level) ) - @. turb_entrʲ_prev_level = limit_turb_entrainment( - entrʲ_prev_level, - turb_entrʲ_prev_level, - draft_area(ρaʲ_prev_level, ρʲ_prev_level), - get_physical_w( - u³ʲ_prev_halflevel, - local_geometry_prev_halflevel, - )/dz_prev_level, - dt, - ) - # We don't have an upper limit to entrainment for the first level # (calculated at i=2), as the vertical velocity at the first level is zero if i > 2 @. entrʲ_prev_level = limit_entrainment( entrʲ_prev_level, - turb_entrʲ_prev_level, draft_area(ρaʲ_prev_level, ρʲ_prev_level), get_physical_w( u³ʲ_prev_halflevel, local_geometry_prev_halflevel, - )/dz_prev_level, - dt, + ), + dz_prev_level, ) end - # @. entrʲ_prev_level = limit_entrainment( - # entrʲ_prev_level, - # draft_area(ρaʲ_prev_level, ρʲ_prev_level), - # dt, - # ) + @. entrʲ_prev_level = limit_entrainment( + entrʲ_prev_level, + draft_area(ρaʲ_prev_level, ρʲ_prev_level), + dt, + ) # TODO: use updraft top instead of scale height @. nh_pressure³ʲ_prev_halflevel = ᶠupdraft_nh_pressure( @@ -586,19 +574,9 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( @. detrʲ_prev_level = limit_detrainment( detrʲ_prev_level, - turb_entrʲ_prev_level, draft_area(ρaʲ_prev_level, ρʲ_prev_level), - get_physical_w( - u³ʲ_prev_halflevel, - local_geometry_prev_halflevel, - )/dz_prev_level, - dt, + dt, ) - # @. detrʲ_prev_level = limit_detrainment( - # detrʲ_prev_level, - # draft_area(ρaʲ_prev_level, ρʲ_prev_level), - # dt, - # ) ρaʲu³ʲ_data = p.scratch.temp_data_level_2 ρaʲu³ʲ_datamse = ρaʲu³ʲ_dataq_tot = p.scratch.temp_data_level_3 diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl index f1ce248239d..449b0cbf53e 100644 --- a/src/cache/prognostic_edmf_precomputed_quantities.jl +++ b/src/cache/prognostic_edmf_precomputed_quantities.jl @@ -233,33 +233,23 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_closures!( p.atmos.edmfx_entr_model, ) - @. ᶜturb_entrʲs.:($$j) = turbulent_entrainment( - params, - draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)) - ) - - @. ᶜturb_entrʲs.:($$j) = limit_turb_entrainment( + @. ᶜentrʲs.:($$j) = limit_entrainment( ᶜentrʲs.:($$j), - ᶜturb_entrʲs.:($$j), draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)), - FT(1.0), dt, ) + @. ᶜturb_entrʲs.:($$j) = turbulent_entrainment( + params, + draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)) + ) - @. ᶜentrʲs.:($$j) = limit_entrainment( + @. ᶜturb_entrʲs.:($$j) = limit_turb_entrainment( ᶜentrʲs.:($$j), ᶜturb_entrʲs.:($$j), - draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)), - FT(1.0), dt, ) - # @. ᶜentrʲs.:($$j) = limit_entrainment( - # ᶜentrʲs.:($$j), - # draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)), - # dt, - # ) @. ᶜvert_div = ᶜdivᵥ(ᶠinterp(ᶜρʲs.:($$j)) * ᶠu³ʲs.:($$j)) / ᶜρʲs.:($$j) @. ᶜmassflux_vert_div = ᶜdivᵥ(ᶠinterp(Y.c.sgsʲs.:($$j).ρa) * ᶠu³ʲs.:($$j)) @@ -286,17 +276,9 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_closures!( @. ᶜdetrʲs.:($$j) = limit_detrainment( ᶜdetrʲs.:($$j), - ᶜturb_entrʲs.:($$j), draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)), - FT(1.0), dt, ) - - # @. ᶜdetrʲs.:($$j) = limit_detrainment( - # ᶜdetrʲs.:($$j), - # draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)), - # dt, - # ) end (; ᶜgradᵥ_θ_virt⁰, ᶜgradᵥ_q_tot⁰, ᶜgradᵥ_θ_liq_ice⁰) = p.precomputed diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl index 9af1886928a..e2983fbe709 100644 --- a/src/prognostic_equations/edmfx_entr_detr.jl +++ b/src/prognostic_equations/edmfx_entr_detr.jl @@ -354,65 +354,20 @@ function edmfx_entr_detr_tendency!(Yₜ, Y, p, t, turbconv_model::PrognosticEDMF return nothing end -# limit_entrainment(entr::FT, a, dt) where {FT} = max( -# min(entr, FT(0.9) * (1 - a) / max(a, eps(FT)) / dt, FT(0.9) * 1 / dt), -# 0, -# ) -# limit_entrainment(entr::FT, a, w, dz) where {FT} = -# max(min(entr, FT(0.9) * w / dz), 0) -# limit_detrainment(detr::FT, a, dt) where {FT} = -# max(min(detr, FT(0.9) * 1 / dt), 0) -# limit_detrainment(detr::FT, a, w, dz) where {FT} = -# max(min(detr, FT(0.9) * w / dz), 0) - - -function limit_entrainment(dyn_entr::FT, turb_entr::FT, a, inv_timescale_lim, dt) where {FT} - - upper_limit = FT(0.9) * min(inv_timescale_lim, FT(1) / dt) - - total_entr = dyn_entr + turb_entr - - if total_entr > upper_limit - dyn_entr, turb_entr = rescale_entr_detr(dyn_entr, turb_entr, upper_limit) - return dyn_entr - else - return dyn_entr - end - -end - -function limit_turb_entrainment(dyn_entr::FT, turb_entr::FT, a, inv_timescale_lim::FT, dt) where {FT} - - upper_limit = FT(0.9) * min(inv_timescale_lim, FT(1) / dt) #w / dz - - total_entr = dyn_entr + turb_entr - - if total_entr > upper_limit - dyn_entr, turb_entr = rescale_entr_detr(dyn_entr, turb_entr, upper_limit) - return turb_entr - else - return turb_entr - end -end - -function limit_detrainment(dyn_detr::FT, turb_entr::FT, a, inv_timescale_lim, dt) where {FT} - - upper_limit = FT(0.9) * min(inv_timescale_lim, FT(1) / dt) - - total_entr = dyn_detr + turb_entr - - if total_entr > upper_limit - dyn_entr, turb_entr = rescale_entr_detr(dyn_detr, turb_entr, upper_limit) - return dyn_detr - else - return dyn_detr - end - +limit_entrainment(entr::FT, a, dt) where {FT} = max( + min(entr, FT(0.9) * (1 - a) / max(a, eps(FT)) / dt, FT(0.9) * 1 / dt), + 0, +) +limit_entrainment(entr::FT, a, w, dz) where {FT} = + max(min(entr, FT(0.9) * w / dz), 0) +limit_detrainment(detr::FT, a, dt) where {FT} = + max(min(detr, FT(0.9) * 1 / dt), 0) +limit_detrainment(detr::FT, a, w, dz) where {FT} = + max(min(detr, FT(0.9) * w / dz), 0) + + +function limit_turb_entrainment(dyn_entr::FT, turb_entr::FT, dt) where {FT} + # Ensure entrainment + turbulent entrainment < 1/dt to promote stability + return min((FT(0.9) * 1 / dt) - dyn_entr, turb_entr) end -function rescale_entr_detr(dyn_entr::FT, turb_entr::FT, upper_limit::FT) where {FT} - total_entr = dyn_entr + turb_entr - if total_entr > upper_limit - return upper_limit * dyn_entr / total_entr, upper_limit * turb_entr / total_entr - end -end diff --git a/src/prognostic_equations/edmfx_tke.jl b/src/prognostic_equations/edmfx_tke.jl index 6aa696d1222..f9d44d95a19 100644 --- a/src/prognostic_equations/edmfx_tke.jl +++ b/src/prognostic_equations/edmfx_tke.jl @@ -13,8 +13,8 @@ function edmfx_tke_tendency!( ) n = n_mass_flux_subdomains(turbconv_model) - (; ᶜentrʲs, ᶜdetrʲs, ᶠu³ʲs) = p.precomputed - (; ᶠu³⁰, ᶜstrain_rate_norm, ᶜlinear_buoygrad, ᶜtke⁰) = p.precomputed + (; ᶜturb_entrʲs, ᶜentrʲs, ᶜdetrʲs, ᶠu³ʲs) = p.precomputed + (; ᶠu³⁰, ᶠu³, ᶜstrain_rate_norm, ᶜlinear_buoygrad, ᶜtke⁰) = p.precomputed (; ᶜK_u, ᶜK_h) = p.precomputed ᶜρa⁰ = turbconv_model isa PrognosticEDMFX ? p.precomputed.ᶜρa⁰ : Y.c.ρ nh_pressure3ʲs = @@ -44,12 +44,18 @@ function edmfx_tke_tendency!( ᶜρaʲ = turbconv_model isa PrognosticEDMFX ? Y.c.sgsʲs.:($j).ρa : p.precomputed.ᶜρaʲs.:($j) + # dynamical entr/detr @. Yₜ.c.sgs⁰.ρatke += ᶜρaʲ * ( ᶜdetrʲs.:($$j) * 1 / 2 * norm_sqr(ᶜinterp(ᶠu³⁰) - ᶜinterp(ᶠu³ʲs.:($$j))) - ᶜentrʲs.:($$j) * ᶜtke⁰ ) + # turbulent entr + @. Yₜ.c.sgs⁰.ρatke += + ᶜρaʲ * ᶜturb_entrʲs.:($$j) * ( + norm(ᶜinterp(ᶠu³⁰) - ᶜinterp(ᶠu³ʲs.:($$j))) * norm(ᶜinterp(ᶠu³⁰) - ᶜinterp(ᶠu³)) - ᶜtke⁰ + ) end # pressure work @. Yₜ.c.sgs⁰.ρatke += ᶜtke_press