Skip to content

Commit

Permalink
Add turb entr to TKE equation, simplify entr/detr limiters
Browse files Browse the repository at this point in the history
  • Loading branch information
@costachris
costachris committed Sep 13, 2024
1 parent 5831acc commit 5703a85
Show file tree
Hide file tree
Showing 4 changed files with 37 additions and 116 deletions.
38 changes: 8 additions & 30 deletions src/cache/diagnostic_edmf_precomputed_quantities.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down Expand Up @@ -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
Expand Down
30 changes: 6 additions & 24 deletions src/cache/prognostic_edmf_precomputed_quantities.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand All @@ -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
Expand Down
75 changes: 15 additions & 60 deletions src/prognostic_equations/edmfx_entr_detr.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
10 changes: 8 additions & 2 deletions src/prognostic_equations/edmfx_tke.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 =
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 5703a85

Please sign in to comment.