Skip to content

Commit 683f35b

Browse files
szy21szy21
authored
Merge pull request #3741 from CliMA/zs/seperate_nh_pressure
Separate buoyancy and drag terms in edmf pressure closure
2 parents cab8399 + ff7438a commit 683f35b

File tree

10 files changed

+156
-77
lines changed

10 files changed

+156
-77
lines changed

‎config/default_configs/default_config.yml

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -319,7 +319,8 @@ implicit_sgs_entr_detr:
319319
help: "Whether to treat the subgrid-scale entrainment and detrainment tendency implicitly [`false` (default), `true`]. Setting it to true only works if implicit_sgs_advection is set to true."
320320
value: false
321321
implicit_sgs_nh_pressure:
322-
help: "Whether to treat the subgrid-scale nonhydrostatic pressure drag tendency implicitly [`false` (default), `true`]. Setting it to true only works if implicit_sgs_advection is set to true."
322+
help: "Whether to treat the subgrid-scale nonhydrostatic pressure closure implicitly [`false` (default), `true`]. Setting it to true only works if implicit_sgs_advection is set to true.
323+
This flag only controls whether the drag term in the pressure closure is treated implicitly. The buoyancy term is always treated explicitly."
323324
value: false
324325
implicit_sgs_mass_flux:
325326
help: "Whether to treat the subgrid-scale mass flux tendency implicitly or explicitly in grid-mean equations. Currently updraft only with Jacobian terms 0. [`false` (default), `true`]. Setting it to true only works if both implicit_sgs_advection and implicit_diffusion are set to true."

‎reproducibility_tests/ref_counter.jl

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
227
1+
228
22

33
# **README**
44
#
@@ -21,6 +21,9 @@
2121

2222
#=
2323
24+
228
25+
- Only treat the drag term in edmf pressure closure implicitly
26+
2427
227
2528
- Move nonhydrostatic pressure drag calculation to implicit precomputed quantities
2629

‎src/cache/diagnostic_edmf_precomputed_quantities.jl

Lines changed: 58 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -223,7 +223,8 @@ function compute_u³ʲ_u³ʲ(
223223
entrʲ_prev_level,
224224
turb_entrʲ_prev_level,
225225
u³⁰_data_prev_halflevel,
226-
nh_pressure³ʲ_data_prev_halflevel,
226+
nh_pressure³_buoyʲ_data_prev_halflevel,
227+
nh_pressure³_dragʲ_data_prev_halflevel,
227228
)
228229
u³ʲ_u³ʲ =
229230
(1 / (J_halflevel^2)) *
@@ -251,8 +252,14 @@ function compute_u³ʲ_u³ʲ(
251252
)
252253

253254
u³ʲ_u³ʲ -=
254-
(1 / (J_halflevel^2)) *
255-
(J_prev_level^2 * 2 * nh_pressure³ʲ_data_prev_halflevel)
255+
(1 / (J_halflevel^2)) * (
256+
J_prev_level^2 *
257+
2 *
258+
(
259+
nh_pressure³_buoyʲ_data_prev_halflevel +
260+
nh_pressure³_dragʲ_data_prev_halflevel
261+
)
262+
)
256263
return u³ʲ_u³ʲ
257264
end
258265

@@ -311,7 +318,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
311318
ᶜentrʲs,
312319
ᶜdetrʲs,
313320
ᶜturb_entrʲs,
314-
ᶠnh_pressure³ʲs,
321+
ᶠnh_pressure³_buoyʲs,
322+
ᶠnh_pressure³_dragʲs,
315323
) = p.precomputed
316324
(; ᶠu³⁰, ᶜK⁰, ᶜtke⁰) = p.precomputed
317325

@@ -396,7 +404,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
396404
ᶜentrʲ = ᶜentrʲs.:($j)
397405
ᶜdetrʲ = ᶜdetrʲs.:($j)
398406
ᶜturb_entrʲ = ᶜturb_entrʲs.:($j)
399-
ᶠnh_pressure³ʲ = ᶠnh_pressure³ʲs.:($j)
407+
ᶠnh_pressure³_buoyʲ = ᶠnh_pressure³_buoyʲs.:($j)
408+
ᶠnh_pressure³_dragʲ = ᶠnh_pressure³_dragʲs.:($j)
400409

401410
if precip_model isa Union{Microphysics0Moment, Microphysics1Moment}
402411
ᶜS_q_totʲ = p.precomputed.ᶜSqₜᵖʲs.:($j)
@@ -430,8 +439,12 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
430439
detrʲ_prev_level = Fields.field_values(Fields.level(ᶜdetrʲ, i - 1))
431440
turb_entrʲ_prev_level =
432441
Fields.field_values(Fields.level(ᶜturb_entrʲ, i - 1))
433-
nh_pressure³ʲ_prev_halflevel =
434-
Fields.field_values(Fields.level(ᶠnh_pressure³ʲ, i - 1 - half))
442+
nh_pressure³_buoyʲ_prev_halflevel = Fields.field_values(
443+
Fields.level(ᶠnh_pressure³_buoyʲ, i - 1 - half),
444+
)
445+
nh_pressure³_dragʲ_prev_halflevel = Fields.field_values(
446+
Fields.level(ᶠnh_pressure³_dragʲ, i - 1 - half),
447+
)
435448
scale_height =
436449
CAP.R_d(params) * CAP.T_surf_ref(params) / CAP.grav(params)
437450

@@ -516,21 +529,29 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
516529

517530
# TODO: use updraft top instead of scale height
518531
if p.atmos.edmfx_model.nh_pressure isa Val{true}
519-
@. nh_pressure³ʲ_prev_halflevel = ᶠupdraft_nh_pressure(
520-
params,
521-
local_geometry_prev_halflevel,
522-
-∇Φ³_prev_level * (ρʲ_prev_level - ρ_prev_level) /
523-
ρʲ_prev_level,
524-
u³ʲ_prev_halflevel,
525-
u³⁰_prev_halflevel,
526-
scale_height,
527-
)
532+
@. nh_pressure³_buoyʲ_prev_halflevel =
533+
ᶠupdraft_nh_pressure_buoyancy(
534+
params,
535+
-∇Φ³_prev_level * (ρʲ_prev_level - ρ_prev_level) /
536+
ρʲ_prev_level,
537+
)
538+
@. nh_pressure³_dragʲ_prev_halflevel =
539+
ᶠupdraft_nh_pressure_drag(
540+
params,
541+
local_geometry_prev_halflevel,
542+
u³ʲ_prev_halflevel,
543+
u³⁰_prev_halflevel,
544+
scale_height,
545+
)
528546
else
529-
@. nh_pressure³ʲ_prev_halflevel = CT3(0)
547+
@. nh_pressure³_buoyʲ_prev_halflevel = CT3(0)
548+
@. nh_pressure³_dragʲ_prev_halflevel = CT3(0)
530549
end
531550

532-
nh_pressure³ʲ_data_prev_halflevel =
533-
nh_pressure³ʲ_prev_halflevel.components.data.:1
551+
nh_pressure³_buoyʲ_data_prev_halflevel =
552+
nh_pressure³_buoyʲ_prev_halflevel.components.data.:1
553+
nh_pressure³_dragʲ_data_prev_halflevel =
554+
nh_pressure³_dragʲ_prev_halflevel.components.data.:1
534555

535556
# Updraft q_tot sources from precipitation formation
536557
# To be applied in updraft continuity, moisture and energy
@@ -574,7 +595,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
574595
entrʲ_prev_level,
575596
turb_entrʲ_prev_level,
576597
u³⁰_data_prev_halflevel,
577-
nh_pressure³ʲ_data_prev_halflevel,
598+
nh_pressure³_buoyʲ_data_prev_halflevel,
599+
nh_pressure³_dragʲ_data_prev_halflevel,
578600
)
579601

580602
# get u³ʲ to calculate divergence term for detrainment,
@@ -842,7 +864,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_top_bc!(
842864
)
843865
n = n_mass_flux_subdomains(p.atmos.turbconv_model)
844866
(; ᶜentrʲs, ᶜdetrʲs, ᶜturb_entrʲs) = p.precomputed
845-
(; ᶠu³⁰, ᶠu³ʲs, ᶜuʲs, ᶠnh_pressure³ʲs) = p.precomputed
867+
(; ᶠu³⁰, ᶠu³ʲs, ᶜuʲs, ᶠnh_pressure³_buoyʲs, ᶠnh_pressure³_dragʲs) =
868+
p.precomputed
846869
(; precip_model) = p.atmos
847870

848871
# set values for the top level
@@ -853,19 +876,26 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_top_bc!(
853876
for j in 1:n
854877
ᶜuʲ = ᶜuʲs.:($j)
855878
ᶠu³ʲ = ᶠu³ʲs.:($j)
856-
ᶠnh_pressure³ʲ = ᶠnh_pressure³ʲs.:($j)
879+
ᶠnh_pressure³_buoyʲ = ᶠnh_pressure³_buoyʲs.:($j)
880+
ᶠnh_pressure³_dragʲ = ᶠnh_pressure³_dragʲs.:($j)
857881
ᶜentrʲ = ᶜentrʲs.:($j)
858882
ᶜdetrʲ = ᶜdetrʲs.:($j)
859883
ᶜturb_entrʲ = ᶜturb_entrʲs.:($j)
860884

861885
u³ʲ_halflevel = Fields.field_values(Fields.level(ᶠu³ʲ, i_top + half))
862886
@. u³ʲ_halflevel = CT3(0)
863-
nh_pressure³ʲ_halflevel =
864-
Fields.field_values(Fields.level(ᶠnh_pressure³ʲ, i_top - half))
865-
@. nh_pressure³ʲ_halflevel = CT3(0)
866-
nh_pressure³ʲ_halflevel =
867-
Fields.field_values(Fields.level(ᶠnh_pressure³ʲ, i_top + half))
868-
@. nh_pressure³ʲ_halflevel = CT3(0)
887+
nh_pressure³_buoyʲ_halflevel =
888+
Fields.field_values(Fields.level(ᶠnh_pressure³_buoyʲ, i_top - half))
889+
@. nh_pressure³_buoyʲ_halflevel = CT3(0)
890+
nh_pressure³_buoyʲ_halflevel =
891+
Fields.field_values(Fields.level(ᶠnh_pressure³_buoyʲ, i_top + half))
892+
@. nh_pressure³_buoyʲ_halflevel = CT3(0)
893+
nh_pressure³_dragʲ_halflevel =
894+
Fields.field_values(Fields.level(ᶠnh_pressure³_dragʲ, i_top - half))
895+
@. nh_pressure³_dragʲ_halflevel = CT3(0)
896+
nh_pressure³_dragʲ_halflevel =
897+
Fields.field_values(Fields.level(ᶠnh_pressure³_dragʲ, i_top + half))
898+
@. nh_pressure³_dragʲ_halflevel = CT3(0)
869899

870900
entrʲ_level = Fields.field_values(Fields.level(ᶜentrʲ, i_top))
871901
detrʲ_level = Fields.field_values(Fields.level(ᶜdetrʲ, i_top))

‎src/cache/precomputed_quantities.jl

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -83,7 +83,7 @@ function implicit_precomputed_quantities(Y, atmos)
8383
ᶠKᵥʲs = similar(Y.f, NTuple{n, FT}),
8484
ᶜtsʲs = similar(Y.c, NTuple{n, TST}),
8585
ᶜρʲs = similar(Y.c, NTuple{n, FT}),
86-
ᶠnh_pressure₃ʲs = similar(Y.f, NTuple{n, C3{FT}}),
86+
ᶠnh_pressure₃_dragʲs = similar(Y.f, NTuple{n, C3{FT}}),
8787
moisture_sgs_quantities...,
8888
) : (;)
8989
return (; gs_quantities..., sgs_quantities..., prognostic_sgs_quantities...)
@@ -172,6 +172,7 @@ function precomputed_quantities(Y, atmos)
172172
ᶜgradᵥ_θ_virt⁰ = Fields.Field(C3{FT}, cspace),
173173
ᶜgradᵥ_q_tot⁰ = Fields.Field(C3{FT}, cspace),
174174
ᶜgradᵥ_θ_liq_ice⁰ = Fields.Field(C3{FT}, cspace),
175+
ᶠnh_pressure₃_buoyʲs = similar(Y.f, NTuple{n, C3{FT}}),
175176
precipitation_sgs_quantities...,
176177
) : (;)
177178

@@ -205,7 +206,8 @@ function precomputed_quantities(Y, atmos)
205206
ᶜentrʲs = similar(Y.c, NTuple{n, FT}),
206207
ᶜdetrʲs = similar(Y.c, NTuple{n, FT}),
207208
ᶜturb_entrʲs = similar(Y.c, NTuple{n, FT}),
208-
ᶠnh_pressure³ʲs = similar(Y.f, NTuple{n, CT3{FT}}),
209+
ᶠnh_pressure³_buoyʲs = similar(Y.f, NTuple{n, CT3{FT}}),
210+
ᶠnh_pressure³_dragʲs = similar(Y.f, NTuple{n, CT3{FT}}),
209211
ᶠu³⁰ = similar(Y.f, CT3{FT}),
210212
ᶜu⁰ = similar(Y.c, C123{FT}),
211213
ᶜK⁰ = similar(Y.c, FT),

‎src/cache/prognostic_edmf_precomputed_quantities.jl

Lines changed: 33 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -325,26 +325,22 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_implicit_clos
325325
(; params) = p
326326
n = n_mass_flux_subdomains(turbconv_model)
327327

328-
(; ᶠgradᵥ_ᶜΦ) = p.core
329-
(; ᶠu₃⁰, ᶜρʲs, ᶠnh_pressure₃ʲs) = p.precomputed
328+
(; ᶠu₃⁰, ᶠnh_pressure₃_dragʲs) = p.precomputed
330329
ᶠlg = Fields.local_geometry_field(Y.f)
331330

332331
scale_height = CAP.R_d(params) * CAP.T_surf_ref(params) / CAP.grav(params)
332+
# nonhydrostatic pressure closure drag term
333333
for j in 1:n
334-
# nonhydrostatic pressure drag
335-
for j in 1:n
336-
if p.atmos.edmfx_model.nh_pressure isa Val{true}
337-
@. ᶠnh_pressure₃ʲs.:($$j) = ᶠupdraft_nh_pressure(
338-
params,
339-
ᶠlg,
340-
ᶠbuoyancy(ᶠinterp(Y.c.ρ), ᶠinterp(ᶜρʲs.:($$j)), ᶠgradᵥ_ᶜΦ),
341-
Y.f.sgsʲs.:($$j).u₃,
342-
ᶠu₃⁰,
343-
scale_height,
344-
)
345-
else
346-
@. ᶠnh_pressure₃ʲs.:($$j) = C3(0)
347-
end
334+
if p.atmos.edmfx_model.nh_pressure isa Val{true}
335+
@. ᶠnh_pressure₃_dragʲs.:($$j) = ᶠupdraft_nh_pressure_drag(
336+
params,
337+
ᶠlg,
338+
Y.f.sgsʲs.:($$j).u₃,
339+
ᶠu₃⁰,
340+
scale_height,
341+
)
342+
else
343+
@. ᶠnh_pressure₃_dragʲs.:($$j) = C3(0)
348344
end
349345
end
350346

@@ -367,6 +363,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
367363

368364
(; params) = p
369365
(; dt) = p
366+
(; ᶠgradᵥ_ᶜΦ) = p.core
370367
thermo_params = CAP.thermodynamics_params(params)
371368
turbconv_params = CAP.turbconv_params(params)
372369

@@ -383,7 +380,16 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
383380
ᶜK_h,
384381
ρatke_flux,
385382
) = p.precomputed
386-
(; ᶜuʲs, ᶜtsʲs, ᶠu³ʲs, ᶜρʲs, ᶜentrʲs, ᶜdetrʲs, ᶜturb_entrʲs) = p.precomputed
383+
(;
384+
ᶜuʲs,
385+
ᶜtsʲs,
386+
ᶠu³ʲs,
387+
ᶜρʲs,
388+
ᶜentrʲs,
389+
ᶜdetrʲs,
390+
ᶜturb_entrʲs,
391+
ᶠnh_pressure₃_buoyʲs,
392+
) = p.precomputed
387393
(; ustar, obukhov_length) = p.precomputed.sfc_conditions
388394

389395
ᶜz = Fields.coordinate_field(Y.c).z
@@ -461,6 +467,16 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
461467
draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)),
462468
dt,
463469
)
470+
471+
# nonhydrostatic pressure closure buoyancy term
472+
if p.atmos.edmfx_model.nh_pressure isa Val{true}
473+
@. ᶠnh_pressure₃_buoyʲs.:($$j) = ᶠupdraft_nh_pressure_buoyancy(
474+
params,
475+
ᶠbuoyancy(ᶠinterp(Y.c.ρ), ᶠinterp(ᶜρʲs.:($$j)), ᶠgradᵥ_ᶜΦ),
476+
)
477+
else
478+
@. ᶠnh_pressure₃_buoyʲs.:($$j) = C3(0)
479+
end
464480
end
465481

466482
(; ᶜgradᵥ_θ_virt⁰, ᶜgradᵥ_q_tot⁰, ᶜgradᵥ_θ_liq_ice⁰) = p.precomputed

‎src/prognostic_equations/edmfx_closures.jl

Lines changed: 39 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -52,46 +52,69 @@ function surface_flux_tke(
5252
end
5353

5454
"""
55-
Return the nonhydrostatic pressure drag for updrafts [m/s2 * m]
55+
Return the virtual mass term of the pressure closure for updrafts [m/s2 * m]
5656
5757
Inputs (everything defined on cell faces):
5858
- params - set with model parameters
59-
- ᶠlg - local geometry (needed to compute the norm inside a local function)
6059
- ᶠbuoyʲ - covariant3 or contravariant3 updraft buoyancy
60+
"""
61+
function ᶠupdraft_nh_pressure_buoyancy(params, ᶠbuoyʲ)
62+
turbconv_params = CAP.turbconv_params(params)
63+
# factor multiplier for pressure buoyancy terms (effective buoyancy is (1-α_b))
64+
α_b = CAP.pressure_normalmode_buoy_coeff1(turbconv_params)
65+
return α_b * ᶠbuoyʲ
66+
end
67+
68+
"""
69+
Return the drag term of the pressure closure for updrafts [m/s2 * m]
70+
71+
Inputs (everything defined on cell faces):
72+
- params - set with model parameters
73+
- ᶠlg - local geometry (needed to compute the norm inside a local function)
6174
- ᶠu3ʲ, ᶠu3⁰ - covariant3 or contravariant3 velocity for updraft and environment.
6275
covariant3 velocity is used in prognostic edmf, and contravariant3
6376
velocity is used in diagnostic edmf.
64-
- updraft top height
77+
- scale height - an approximation for updraft top height
6578
"""
66-
function ᶠupdraft_nh_pressure(params, ᶠlg, ᶠbuoyʲ, ᶠu3ʲ, ᶠu3⁰, plume_height)
79+
function ᶠupdraft_nh_pressure_drag(params, ᶠlg, ᶠu3ʲ, ᶠu3⁰, scale_height)
6780
turbconv_params = CAP.turbconv_params(params)
68-
# factor multiplier for pressure buoyancy terms (effective buoyancy is (1-α_b))
69-
α_b = CAP.pressure_normalmode_buoy_coeff1(turbconv_params)
7081
# factor multiplier for pressure drag
7182
α_d = CAP.pressure_normalmode_drag_coeff(turbconv_params)
72-
73-
# Independence of aspect ratio hardcoded: α₂_asp_ratio² = FT(0)
74-
7583
H_up_min = CAP.min_updraft_top(turbconv_params)
7684

85+
# Independence of aspect ratio hardcoded: α₂_asp_ratio² = FT(0)
7786
# We also used to have advection term here: α_a * w_up * div_w_up
78-
return α_b * ᶠbuoyʲ +
79-
α_d * (ᶠu3ʲ - ᶠu3⁰) * CC.Geometry._norm(ᶠu3ʲ - ᶠu3⁰, ᶠlg) /
80-
max(plume_height, H_up_min)
87+
return α_d * (ᶠu3ʲ - ᶠu3⁰) * CC.Geometry._norm(ᶠu3ʲ - ᶠu3⁰, ᶠlg) /
88+
max(scale_height, H_up_min)
89+
end
90+
91+
edmfx_nh_pressure_buoyancy_tendency!(Yₜ, Y, p, t, turbconv_model) = nothing
92+
function edmfx_nh_pressure_buoyancy_tendency!(
93+
Yₜ,
94+
Y,
95+
p,
96+
t,
97+
turbconv_model::PrognosticEDMFX,
98+
)
99+
(; ᶠnh_pressure₃_buoyʲs) = p.precomputed
100+
n = n_mass_flux_subdomains(turbconv_model)
101+
for j in 1:n
102+
@. Yₜ.f.sgsʲs.:($$j).u₃ -= ᶠnh_pressure₃_buoyʲs.:($$j)
103+
end
81104
end
82105

83-
edmfx_nh_pressure_tendency!(Yₜ, Y, p, t, turbconv_model) = nothing
84-
function edmfx_nh_pressure_tendency!(
106+
edmfx_nh_pressure_drag_tendency!(Yₜ, Y, p, t, turbconv_model) = nothing
107+
function edmfx_nh_pressure_drag_tendency!(
85108
Yₜ,
86109
Y,
87110
p,
88111
t,
89112
turbconv_model::PrognosticEDMFX,
90113
)
91-
(; ᶠnh_pressure₃ʲs) = p.precomputed
114+
(; ᶠnh_pressure₃_dragʲs) = p.precomputed
92115
n = n_mass_flux_subdomains(turbconv_model)
93116
for j in 1:n
94-
@. Yₜ.f.sgsʲs.:($$j).u₃ -= ᶠnh_pressure₃ʲs.:($$j)
117+
@. Yₜ.f.sgsʲs.:($$j).u₃ -= ᶠnh_pressure₃_dragʲs.:($$j)
95118
end
96119
end
97120

0 commit comments

Comments
 (0)