Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

use z-z_sfc in edmf closure #2041

Merged
merged 1 commit into from
Sep 1, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions src/cache/diagnostic_edmf_precomputed_quantities.jl
Original file line number Diff line number Diff line change
Expand Up @@ -356,6 +356,7 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
params,
p.atmos.edmfx_entr_detr,
z_prev_level,
z_sfc_halflevel,
p_prev_level,
ρ_prev_level,
buoyancy_flux_sfc_halflevel,
Expand Down Expand Up @@ -674,10 +675,12 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
end

sfc_tke = Fields.level(ᶜtke⁰, 1)
z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, half)
@. ᶜmixing_length = mixing_length(
params,
ustar,
ᶜz,
z_sfc,
ᶜdz,
max(sfc_tke, 0),
ᶜlinear_buoygrad,
Expand Down
3 changes: 3 additions & 0 deletions src/cache/edmf_precomputed_quantities.jl
Original file line number Diff line number Diff line change
Expand Up @@ -145,6 +145,7 @@ function set_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
end

ᶜz = Fields.coordinate_field(Y.c).z
z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
ᶜdz = Fields.Δz_field(axes(Y.c))
ᶜlg = Fields.local_geometry_field(Y.c)

Expand All @@ -153,6 +154,7 @@ function set_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
params,
p.atmos.edmfx_entr_detr,
ᶜz,
z_sfc,
ᶜp,
Y.c.ρ,
buoyancy_flux,
Expand Down Expand Up @@ -215,6 +217,7 @@ function set_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
p.params,
ustar,
ᶜz,
z_sfc,
ᶜdz,
sfc_tke,
ᶜlinear_buoygrad,
Expand Down
20 changes: 13 additions & 7 deletions src/prognostic_equations/edmfx_closures.jl
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,7 @@ function edmfx_nh_pressure_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::EDMF
(; params, ᶜρʲs, ᶜρ_ref, ᶠgradᵥ_ᶜΦ, ᶜuʲs, ᶜu⁰, ᶠu₃⁰) = p
FT = eltype(Y)
ᶜz = Fields.coordinate_field(Y.c).z
z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
ᶠlg = Fields.local_geometry_field(Y.f)

turbconv_params = CAP.turbconv_params(params)
Expand All @@ -90,6 +91,7 @@ function edmfx_nh_pressure_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::EDMF
updraft_top = Spaces.level(ᶜz[colidx], level)[]
end
end
updraft_top = updraft_top - z_sfc[colidx][]

@. Yₜ.f.sgsʲs.:($$j).u₃[colidx] -= ᶠupdraft_nh_pressure(
params,
Expand Down Expand Up @@ -130,6 +132,7 @@ function pi_groups_entr_detr(
params,
entr_detr_flag::Bool,
ᶜz::FT,
z_sfc::FT,
ᶜp::FT,
ᶜρ::FT,
buoy_flux_surface::FT,
Expand Down Expand Up @@ -162,13 +165,15 @@ function pi_groups_entr_detr(

# non-dimensional pi-groups
# TODO - using Π₁ blows things up
Π₁ = ᶜz * (ᶜbuoyʲ - ᶜbuoy⁰) / ((ᶜwʲ - ᶜw⁰)^2 + w_star^2 + eps(FT))
Π₁ =
(ᶜz - z_sfc) * (ᶜbuoyʲ - ᶜbuoy⁰) /
((ᶜwʲ - ᶜw⁰)^2 + w_star^2 + eps(FT))
Π₃ = sqrt(ᶜaʲ)
Π₄ = ᶜRHʲ - ᶜRH⁰
Π₆ = ᶜz / ref_H
Π₆ = (ᶜz - z_sfc) / ref_H
entr = max(
0,
ᶜwʲ / ᶜz * (
ᶜwʲ / (ᶜz - z_sfc) * (
-4.013288 - 0.000968 * Π₁ + 0.356974 * Π₃ - 0.403124 * Π₄ +
1.503261 * Π₆
),
Expand All @@ -184,7 +189,7 @@ function pi_groups_entr_detr(
# TODO - Temporary: Switch to Π groups after simple tests are done
# (kinematic, bubble, Bomex)
# and/or we can calibrate things in ClimaAtmos
entr = max(0, min(entr_coeff * ᶜwʲ / ᶜz, 1 / dt))
entr = max(0, min(entr_coeff * ᶜwʲ / (ᶜz - z_sfc), 1 / dt))
detr = max(0, min(detr_coeff * ᶜwʲ, 1 / dt))

return (; entr, detr)
Expand Down Expand Up @@ -269,6 +274,7 @@ function mixing_length(
params,
ustar::FT,
ᶜz::FT,
z_sfc::FT,
ᶜdz::FT,
sfc_tke::FT,
ᶜlinear_buoygrad::FT,
Expand All @@ -293,10 +299,10 @@ function mixing_length(
# kz scale (surface layer)
if obukhov_length < 0.0 #unstable
l_W =
vkc * ᶜz / (sqrt(sfc_tke / ustar / ustar) * c_m) *
min((1 - 100 * ᶜz / obukhov_length)^FT(0.2), 1 / vkc)
vkc * (ᶜz - z_sfc) / (sqrt(sfc_tke / ustar / ustar) * c_m) *
min((1 - 100 * (ᶜz - z_sfc) / obukhov_length)^FT(0.2), 1 / vkc)
else # neutral or stable
l_W = vkc * ᶜz / (sqrt(sfc_tke / ustar / ustar) * c_m)
l_W = vkc * (ᶜz - z_sfc) / (sqrt(sfc_tke / ustar / ustar) * c_m)
end

# compute l_TKE - the production-dissipation balanced length scale
Expand Down
Loading