Skip to content

Commit 7d34bfb

Browse files
trontryteltrontrytel
committed
TEsting simplified hyperdiffusion
1 parent d248230 commit 7d34bfb

File tree

2 files changed

+85
-206
lines changed

2 files changed

+85
-206
lines changed

‎src/prognostic_equations/hyperdiffusion.jl‎

Lines changed: 82 additions & 202 deletions
Original file line numberDiff line numberDiff line change
@@ -50,40 +50,20 @@ function hyperdiffusion_cache(
5050
gs_quantities = (;
5151
ᶜ∇²u = similar(Y.c, C123{FT}),
5252
ᶜ∇²specific_energy = similar(Y.c, FT),
53-
ᶜ∇²specific_tracers = Base.materialize(ᶜspecific_gs_tracers(Y)),
53+
ᶜ∇²q_tot = similar(Y.c, FT),
5454
)
5555

5656
# Sub-grid scale quantities
5757
ᶜ∇²uʲs =
5858
turbconv_model isa PrognosticEDMFX ? similar(Y.c, NTuple{n, C123{FT}}) :
5959
(;)
60-
moisture_sgs_quantities =
61-
moisture_model isa NonEquilMoistModel &&
62-
microphysics_model isa Microphysics1Moment ?
63-
(;
64-
ᶜ∇²q_liqʲs = similar(Y.c, NTuple{n, FT}),
65-
ᶜ∇²q_iceʲs = similar(Y.c, NTuple{n, FT}),
66-
ᶜ∇²q_raiʲs = similar(Y.c, NTuple{n, FT}),
67-
ᶜ∇²q_snoʲs = similar(Y.c, NTuple{n, FT}),
68-
) :
69-
moisture_model isa NonEquilMoistModel &&
70-
microphysics_model isa Microphysics2Moment ?
71-
(;
72-
ᶜ∇²q_liqʲs = similar(Y.c, NTuple{n, FT}),
73-
ᶜ∇²q_iceʲs = similar(Y.c, NTuple{n, FT}),
74-
ᶜ∇²q_raiʲs = similar(Y.c, NTuple{n, FT}),
75-
ᶜ∇²q_snoʲs = similar(Y.c, NTuple{n, FT}),
76-
ᶜ∇²n_liqʲs = similar(Y.c, NTuple{n, FT}),
77-
ᶜ∇²n_raiʲs = similar(Y.c, NTuple{n, FT}),
78-
) : (;)
7960
sgs_quantities =
8061
turbconv_model isa PrognosticEDMFX ?
8162
(;
8263
ᶜ∇²uₕʲs = similar(Y.c, NTuple{n, C12{FT}}),
8364
ᶜ∇²uᵥʲs = similar(Y.c, NTuple{n, C3{FT}}),
8465
ᶜ∇²mseʲs = similar(Y.c, NTuple{n, FT}),
8566
ᶜ∇²q_totʲs = similar(Y.c, NTuple{n, FT}),
86-
moisture_sgs_quantities...,
8767
) : (;)
8868
maybe_ᶜ∇²tke⁰ =
8969
use_prognostic_tke(turbconv_model) ? (; ᶜ∇²tke⁰ = similar(Y.c, FT)) :
@@ -112,9 +92,9 @@ NVTX.@annotate function prep_hyperdiffusion_tendency!(Yₜ, Y, p, t)
11292
diffuse_tke = use_prognostic_tke(turbconv_model)
11393
(; ᶜp) = p.precomputed
11494
(; ᶜh_ref) = p.core
115-
(; ᶜ∇²u, ᶜ∇²specific_energy) = p.hyperdiff
95+
(; ᶜ∇²u, ᶜ∇²specific_energy, ᶜ∇²q_tot) = p.hyperdiff
11696
if turbconv_model isa PrognosticEDMFX
117-
(; ᶜ∇²uₕʲs, ᶜ∇²uᵥʲs, ᶜ∇²uʲs, ᶜ∇²mseʲs) = p.hyperdiff
97+
(; ᶜ∇²uₕʲs, ᶜ∇²uᵥʲs, ᶜ∇²uʲs, ᶜ∇²mseʲs, ᶜ∇²q_totʲs) = p.hyperdiff
11898
end
11999

120100
# Grid scale hyperdiffusion
@@ -125,6 +105,8 @@ NVTX.@annotate function prep_hyperdiffusion_tendency!(Yₜ, Y, p, t)
125105
@. ᶜ∇²specific_energy =
126106
wdivₕ(gradₕ(specific(Y.c.ρe_tot, Y.c.ρ) + ᶜp / Y.c.ρ - ᶜh_ref))
127107

108+
@. ᶜ∇²q_tot = wdivₕ(gradₕ(specific(Y.c.ρq_tot, Y.c.ρ)))
109+
128110
if diffuse_tke
129111
ᶜρa⁰ =
130112
turbconv_model isa PrognosticEDMFX ?
@@ -142,6 +124,7 @@ NVTX.@annotate function prep_hyperdiffusion_tendency!(Yₜ, Y, p, t)
142124
C123(wgradₕ(divₕ(p.precomputed.ᶜuʲs.:($$j)))) -
143125
C123(wcurlₕ(C123(curlₕ(p.precomputed.ᶜuʲs.:($$j)))))
144126
@. ᶜ∇²mseʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).mse))
127+
@. ᶜ∇²q_totʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).q_tot))
145128
@. ᶜ∇²uₕʲs.:($$j) = C12(ᶜ∇²uʲs.:($$j))
146129
@. ᶜ∇²uᵥʲs.:($$j) = C3(ᶜ∇²uʲs.:($$j))
147130
end
@@ -156,41 +139,71 @@ NVTX.@annotate function apply_hyperdiffusion_tendency!(Yₜ, Y, p, t)
156139

157140
(; divergence_damping_factor) = hyperdiff
158141
(; ν₄_scalar, ν₄_vorticity) = ν₄(hyperdiff, Y)
142+
# Rescale the hyperdiffusivity for precipitating species.
143+
ν₄_scalar_for_precip = CAP.α_hyperdiff_tracer(p.params) * ν₄_scalar
159144

160145
n = n_mass_flux_subdomains(turbconv_model)
161146
diffuse_tke = use_prognostic_tke(turbconv_model)
162147
ᶜJ = Fields.local_geometry_field(Y.c).J
163148
point_type = eltype(Fields.coordinate_field(Y.c))
164-
(; ᶜ∇²u, ᶜ∇²specific_energy) = p.hyperdiff
149+
150+
(; ᶜ∇²u, ᶜ∇²specific_energy, ᶜ∇²q_tot) = p.hyperdiff
151+
165152
if turbconv_model isa PrognosticEDMFX
166153
ᶜρa⁰ = @. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, turbconv_model))
167-
(; ᶜ∇²uₕʲs, ᶜ∇²uᵥʲs, ᶜ∇²uʲs, ᶜ∇²mseʲs) = p.hyperdiff
154+
(; ᶜ∇²uₕʲs, ᶜ∇²uᵥʲs, ᶜ∇²uʲs, ᶜ∇²mseʲs, ᶜ∇²q_totʲs) = p.hyperdiff
168155
end
169156
if use_prognostic_tke(turbconv_model)
170157
(; ᶜ∇²tke⁰) = p.hyperdiff
171158
end
172159

173-
if turbconv_model isa PrognosticEDMFX
174-
for j in 1:n
175-
@. ᶜ∇²uʲs.:($$j) = C123(ᶜ∇²uₕʲs.:($$j)) + C123(ᶜ∇²uᵥʲs.:($$j))
176-
end
177-
end
178-
160+
###
161+
### GS velocities
162+
###
179163
# re-use to store the curl-curl part
180164
@. ᶜ∇²u =
181165
divergence_damping_factor * C123(wgradₕ(divₕ(ᶜ∇²u))) -
182166
C123(wcurlₕ(C123(curlₕ(ᶜ∇²u))))
183167
@. Yₜ.c.uₕ -= ν₄_vorticity * C12(ᶜ∇²u)
184168
@. Yₜ.f.u₃ -= ν₄_vorticity * ᶠwinterp(ᶜJ * Y.c.ρ, C3(ᶜ∇²u))
185169

170+
###
171+
### GS energy, density and total moisture
172+
###
186173
@. Yₜ.c.ρe_tot -= ν₄_scalar * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²specific_energy))
174+
@. Yₜ.c.ρq_tot -= ν₄_scalar * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²q_tot))
175+
@. Yₜ.c.ρ -= ν₄_scalar * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²q_tot))
187176

188-
# Sub-grid scale hyperdiffusion continued
189-
if (turbconv_model isa PrognosticEDMFX) && diffuse_tke
190-
@. Yₜ.c.sgs⁰.ρatke -= ν₄_vorticity * wdivₕ(ᶜρa⁰ * gradₕ(ᶜ∇²tke⁰))
177+
# TODO: Since we are not applying the limiter to density (or area-weighted
178+
# density), the mass redistributed by hyperdiffusion will not be conserved
179+
# by the limiter. Is this a significant problem?
180+
181+
###
182+
### GS moisture tracers
183+
###
184+
FT = eltype(Y.c.ρ)
185+
if p.atmos.moisture_model isa NonEquilMoistModel && (p.atmos.microphysics_model isa Microphysics1Moment || p.atmos.microphysics_model isa Microphysics2Moment)
186+
# cloud condensate mass
187+
@. Yₜ.c.ρq_liq -= ν₄_scalar * Y.c.ρq_liq / max(eps(FT), Y.c.ρq_tot) * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²q_tot))
188+
@. Yₜ.c.ρq_ice -= ν₄_scalar * Y.c.ρq_ice / max(eps(FT), Y.c.ρq_tot) * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²q_tot))
189+
# precipitation mass
190+
@. Yₜ.c.ρq_rai -= ν₄_scalar_for_precip * Y.c.ρq_rai / max(eps(FT), Y.c.ρq_tot) * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²q_tot))
191+
@. Yₜ.c.ρq_sno -= ν₄_scalar_for_precip * Y.c.ρq_sno / max(eps(FT), Y.c.ρq_tot) * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²q_tot))
191192
end
193+
if p.atmos.moisture_model isa NonEquilMoistModel && p.atmos.microphysics_model isa Microphysics2Moment
194+
# number concnetrations
195+
# TODO - should I multiply by som reference number concentration?
196+
@. Yₜ.c.ρn_liq -= ν₄_scalar * Y.c.ρq_liq / max(eps(FT), Y.c.ρq_tot) * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²q_tot))
197+
@. Yₜ.c.ρn_rai -= ν₄_scalar_for_precip * Y.c.ρq_rai / max(eps(FT), Y.c.ρq_tot) * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²q_tot))
198+
end
199+
200+
###
201+
### SGS
202+
###
192203
if turbconv_model isa PrognosticEDMFX
193204
for j in 1:n
205+
@. ᶜ∇²uʲs.:($$j) = C123(ᶜ∇²uₕʲs.:($$j)) + C123(ᶜ∇²uᵥʲs.:($$j))
206+
194207
if point_type <: Geometry.Abstract3DPoint
195208
# only need curl-curl part
196209
@. ᶜ∇²uᵥʲs.:($$j) = C3(wcurlₕ(C123(curlₕ(ᶜ∇²uʲs.:($$j)))))
@@ -200,9 +213,39 @@ NVTX.@annotate function apply_hyperdiffusion_tendency!(Yₜ, Y, p, t)
200213
# Note: It is more correct to have ρa inside and outside the divergence
201214
@. Yₜ.c.sgsʲs.:($$j).mse -=
202215
ν₄_scalar * wdivₕ(gradₕ(ᶜ∇²mseʲs.:($$j)))
216+
217+
@. Yₜ.c.sgsʲs.:($$j).ρa -=
218+
ν₄_scalar *
219+
wdivₕ(Y.c.sgsʲs.:($$j).ρa * gradₕ(ᶜ∇²q_totʲs.:($$j)))
220+
@. Yₜ.c.sgsʲs.:($$j).q_tot -=
221+
ν₄_scalar * wdivₕ(gradₕ(ᶜ∇²q_totʲs.:($$j)))
222+
223+
if p.atmos.moisture_model isa NonEquilMoistModel &&
224+
(p.atmos.microphysics_model isa Microphysics1Moment || p.atmos.microphysics_model isa Microphysics2Moment)
225+
@. Yₜ.c.sgsʲs.:($$j).q_liq -=
226+
ν₄_scalar * Y.c.sgsʲs.q_liq / max(FT(0), Y.c.sgsʲs.q_tot) * wdivₕ(gradₕ(ᶜ∇²q_totʲs.:($$j)))
227+
@. Yₜ.c.sgsʲs.:($$j).q_ice -=
228+
ν₄_scalar * Y.c.sgsʲs.q_ice / max(FT(0), Y.c.sgsʲs.q_tot) * wdivₕ(gradₕ(ᶜ∇²q_totʲs.:($$j)))
229+
@. Yₜ.c.sgsʲs.:($$j).q_rai -=
230+
ν₄_scalar_for_precip * Y.c.sgsʲs.q_rai / max(FT(0), Y.c.sgsʲs.q_tot) * wdivₕ(gradₕ(ᶜ∇²q_totʲs.:($$j)))
231+
@. Yₜ.c.sgsʲs.:($$j).q_sno -=
232+
ν₄_scalar_for_precip * Y.c.sgsʲs.q_sno / max(FT(0), Y.c.sgsʲs.q_tot) * wdivₕ(gradₕ(ᶜ∇²q_totʲs.:($$j)))
233+
end
234+
if p.atmos.moisture_model isa NonEquilMoistModel && p.atmos.microphysics_model isa Microphysics2Moment
235+
@. Yₜ.c.sgsʲs.:($$j).n_liq -=
236+
ν₄_scalar * Y.c.sgsʲs.q_liq / max(FT(0), Y.c.sgsʲs.q_tot) * wdivₕ(gradₕ(ᶜ∇²q_totʲs.:($$j)))
237+
@. Yₜ.c.sgsʲs.:($$j).n_rai -=
238+
ν₄_scalar_for_precip * Y.c.sgsʲs.q_rai / max(FT(0), Y.c.sgsʲs.q_tot) * wdivₕ(gradₕ(ᶜ∇²q_totʲs.:($$j)))
239+
end
203240
end
204241
end
205242

243+
###
244+
### SGS TKE
245+
###
246+
if (turbconv_model isa PrognosticEDMFX) && diffuse_tke
247+
@. Yₜ.c.sgs⁰.ρatke -= ν₄_vorticity * wdivₕ(ᶜρa⁰ * gradₕ(ᶜ∇²tke⁰))
248+
end
206249
if turbconv_model isa DiagnosticEDMFX && diffuse_tke
207250
@. Yₜ.c.sgs⁰.ρatke -= ν₄_vorticity * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²tke⁰))
208251
end
@@ -211,10 +254,10 @@ end
211254
function dss_hyperdiffusion_tendency_pairs(p)
212255
(; hyperdiff, turbconv_model) = p.atmos
213256
buffer = p.hyperdiff.hyperdiffusion_ghost_buffer
214-
(; ᶜ∇²u, ᶜ∇²specific_energy) = p.hyperdiff
257+
(; ᶜ∇²u, ᶜ∇²specific_energy, ᶜ∇²q_tot) = p.hyperdiff
215258
diffuse_tke = use_prognostic_tke(turbconv_model)
216259
if turbconv_model isa PrognosticEDMFX
217-
(; ᶜ∇²uₕʲs, ᶜ∇²uᵥʲs, ᶜ∇²mseʲs) = p.hyperdiff
260+
(; ᶜ∇²uₕʲs, ᶜ∇²uᵥʲs, ᶜ∇²mseʲs, ᶜ∇²q_totʲs) = p.hyperdiff
218261
end
219262
if use_prognostic_tke(turbconv_model)
220263
(; ᶜ∇²tke⁰) = p.hyperdiff
@@ -223,6 +266,7 @@ function dss_hyperdiffusion_tendency_pairs(p)
223266
core_dynamics_pairs = (
224267
ᶜ∇²u => buffer.ᶜ∇²u,
225268
ᶜ∇²specific_energy => buffer.ᶜ∇²specific_energy,
269+
ᶜ∇²q_tot => buffer.ᶜ∇²q_tot,
226270
(diffuse_tke ? (ᶜ∇²tke⁰ => buffer.ᶜ∇²tke⁰,) : ())...,
227271
)
228272
tc_dynamics_pairs =
@@ -231,171 +275,7 @@ function dss_hyperdiffusion_tendency_pairs(p)
231275
ᶜ∇²uₕʲs => buffer.ᶜ∇²uₕʲs,
232276
ᶜ∇²uᵥʲs => buffer.ᶜ∇²uᵥʲs,
233277
ᶜ∇²mseʲs => buffer.ᶜ∇²mseʲs,
278+
ᶜ∇²q_totʲs => buffer.ᶜ∇²q_totʲs,
234279
) : ()
235-
dynamics_pairs = (core_dynamics_pairs..., tc_dynamics_pairs...)
236-
237-
(; ᶜ∇²specific_tracers) = p.hyperdiff
238-
core_tracer_pairs =
239-
!isempty(propertynames(ᶜ∇²specific_tracers)) ?
240-
(ᶜ∇²specific_tracers => buffer.ᶜ∇²specific_tracers,) : ()
241-
tc_tracer_pairs =
242-
turbconv_model isa PrognosticEDMFX ?
243-
(p.hyperdiff.ᶜ∇²q_totʲs => buffer.ᶜ∇²q_totʲs,) : ()
244-
tc_moisture_pairs =
245-
turbconv_model isa PrognosticEDMFX &&
246-
p.atmos.moisture_model isa NonEquilMoistModel &&
247-
p.atmos.microphysics_model isa Microphysics1Moment ?
248-
(
249-
p.hyperdiff.ᶜ∇²q_liqʲs => buffer.ᶜ∇²q_liqʲs,
250-
p.hyperdiff.ᶜ∇²q_iceʲs => buffer.ᶜ∇²q_iceʲs,
251-
p.hyperdiff.ᶜ∇²q_raiʲs => buffer.ᶜ∇²q_raiʲs,
252-
p.hyperdiff.ᶜ∇²q_snoʲs => buffer.ᶜ∇²q_snoʲs,
253-
) :
254-
turbconv_model isa PrognosticEDMFX &&
255-
p.atmos.moisture_model isa NonEquilMoistModel &&
256-
p.atmos.microphysics_model isa Microphysics2Moment ?
257-
(
258-
p.hyperdiff.ᶜ∇²q_liqʲs => buffer.ᶜ∇²q_liqʲs,
259-
p.hyperdiff.ᶜ∇²q_iceʲs => buffer.ᶜ∇²q_iceʲs,
260-
p.hyperdiff.ᶜ∇²q_raiʲs => buffer.ᶜ∇²q_raiʲs,
261-
p.hyperdiff.ᶜ∇²q_snoʲs => buffer.ᶜ∇²q_snoʲs,
262-
p.hyperdiff.ᶜ∇²q_liqʲs => buffer.ᶜ∇²q_liqʲs,
263-
p.hyperdiff.ᶜ∇²q_raiʲs => buffer.ᶜ∇²q_raiʲs,
264-
) : ()
265-
tracer_pairs =
266-
(core_tracer_pairs..., tc_tracer_pairs..., tc_moisture_pairs...)
267-
return (dynamics_pairs..., tracer_pairs...)
268-
end
269-
270-
# This should prep variables that we will dss in
271-
# dss_hyperdiffusion_tendency_pairs
272-
NVTX.@annotate function prep_tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t)
273-
(; hyperdiff, turbconv_model) = p.atmos
274-
isnothing(hyperdiff) && return nothing
275-
276-
(; ᶜ∇²specific_tracers) = p.hyperdiff
277-
278-
# TODO: Fix RecursiveApply bug in gradₕ to fuse this operation.
279-
# ᶜ∇²specific_tracers .= wdivₕ.(gradₕ.(ᶜspecific_gs_tracers(Y)))
280-
foreach_gs_tracer(Y, ᶜ∇²specific_tracers) do ᶜρχ, ᶜ∇²χ, _
281-
@. ᶜ∇²χ = wdivₕ(gradₕ(specific(ᶜρχ, Y.c.ρ)))
282-
end
283-
284-
if turbconv_model isa PrognosticEDMFX
285-
n = n_mass_flux_subdomains(turbconv_model)
286-
(; ᶜ∇²q_totʲs) = p.hyperdiff
287-
for j in 1:n
288-
# Note: It is more correct to have ρa inside and outside the divergence
289-
@. ᶜ∇²q_totʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).q_tot))
290-
end
291-
if p.atmos.moisture_model isa NonEquilMoistModel &&
292-
p.atmos.microphysics_model isa Microphysics1Moment
293-
(; ᶜ∇²q_liqʲs, ᶜ∇²q_iceʲs, ᶜ∇²q_raiʲs, ᶜ∇²q_snoʲs) = p.hyperdiff
294-
for j in 1:n
295-
# Note: It is more correct to have ρa inside and outside the divergence
296-
@. ᶜ∇²q_liqʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).q_liq))
297-
@. ᶜ∇²q_iceʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).q_ice))
298-
@. ᶜ∇²q_raiʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).q_rai))
299-
@. ᶜ∇²q_snoʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).q_sno))
300-
end
301-
elseif p.atmos.moisture_model isa NonEquilMoistModel &&
302-
p.atmos.microphysics_model isa Microphysics2Moment
303-
(;
304-
ᶜ∇²q_liqʲs,
305-
ᶜ∇²q_iceʲs,
306-
ᶜ∇²q_raiʲs,
307-
ᶜ∇²q_snoʲs,
308-
ᶜ∇²n_liqʲs,
309-
ᶜ∇²n_raiʲs,
310-
) = p.hyperdiff
311-
for j in 1:n
312-
# Note: It is more correct to have ρa inside and outside the divergence
313-
@. ᶜ∇²q_liqʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).q_liq))
314-
@. ᶜ∇²q_iceʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).q_ice))
315-
@. ᶜ∇²q_raiʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).q_rai))
316-
@. ᶜ∇²q_snoʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).q_sno))
317-
@. ᶜ∇²n_liqʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).n_liq))
318-
@. ᶜ∇²n_raiʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).n_rai))
319-
end
320-
end
321-
end
322-
return nothing
323-
end
324-
325-
# This requires dss to have been called on
326-
# variables in dss_hyperdiffusion_tendency_pairs
327-
NVTX.@annotate function apply_tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t)
328-
(; hyperdiff, turbconv_model) = p.atmos
329-
isnothing(hyperdiff) && return nothing
330-
331-
# Rescale the hyperdiffusivity for precipitating species.
332-
(; ν₄_scalar) = ν₄(hyperdiff, Y)
333-
ν₄_scalar_for_precip = CAP.α_hyperdiff_tracer(p.params) * ν₄_scalar
334-
335-
n = n_mass_flux_subdomains(turbconv_model)
336-
(; ᶜ∇²specific_tracers) = p.hyperdiff
337-
338-
# TODO: Since we are not applying the limiter to density (or area-weighted
339-
# density), the mass redistributed by hyperdiffusion will not be conserved
340-
# by the limiter. Is this a significant problem?
341-
foreach_gs_tracer(Yₜ, ᶜ∇²specific_tracers) do ᶜρχₜ, ᶜ∇²χ, ρχ_name
342-
ν₄_scalar_for_χ =
343-
ρχ_name in (@name(ρq_rai), @name(ρq_sno), @name(ρn_rai)) ?
344-
ν₄_scalar_for_precip : ν₄_scalar
345-
@. ᶜρχₜ -= ν₄_scalar_for_χ * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²χ))
346-
347-
# Take into account the effect of total water diffusion on density.
348-
if ρχ_name == @name(ρq_tot)
349-
@. Yₜ.c.ρ -= ν₄_scalar * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²χ))
350-
end
351-
end
352-
if turbconv_model isa PrognosticEDMFX
353-
(; ᶜ∇²q_totʲs) = p.hyperdiff
354-
for j in 1:n
355-
@. Yₜ.c.sgsʲs.:($$j).ρa -=
356-
ν₄_scalar *
357-
wdivₕ(Y.c.sgsʲs.:($$j).ρa * gradₕ(ᶜ∇²q_totʲs.:($$j)))
358-
@. Yₜ.c.sgsʲs.:($$j).q_tot -=
359-
ν₄_scalar * wdivₕ(gradₕ(ᶜ∇²q_totʲs.:($$j)))
360-
end
361-
if p.atmos.moisture_model isa NonEquilMoistModel &&
362-
p.atmos.microphysics_model isa Microphysics1Moment
363-
(; ᶜ∇²q_liqʲs, ᶜ∇²q_iceʲs, ᶜ∇²q_raiʲs, ᶜ∇²q_snoʲs) = p.hyperdiff
364-
for j in 1:n
365-
@. Yₜ.c.sgsʲs.:($$j).q_liq -=
366-
ν₄_scalar * wdivₕ(gradₕ(ᶜ∇²q_liqʲs.:($$j)))
367-
@. Yₜ.c.sgsʲs.:($$j).q_ice -=
368-
ν₄_scalar * wdivₕ(gradₕ(ᶜ∇²q_iceʲs.:($$j)))
369-
@. Yₜ.c.sgsʲs.:($$j).q_rai -=
370-
ν₄_scalar_for_precip * wdivₕ(gradₕ(ᶜ∇²q_raiʲs.:($$j)))
371-
@. Yₜ.c.sgsʲs.:($$j).q_sno -=
372-
ν₄_scalar_for_precip * wdivₕ(gradₕ(ᶜ∇²q_snoʲs.:($$j)))
373-
end
374-
elseif p.atmos.moisture_model isa NonEquilMoistModel &&
375-
p.atmos.microphysics_model isa Microphysics2Moment
376-
(;
377-
ᶜ∇²q_liqʲs,
378-
ᶜ∇²q_iceʲs,
379-
ᶜ∇²q_raiʲs,
380-
ᶜ∇²q_snoʲs,
381-
ᶜ∇²n_liqʲs,
382-
ᶜ∇²n_raiʲs,
383-
) = p.hyperdiff
384-
for j in 1:n
385-
@. Yₜ.c.sgsʲs.:($$j).q_liq -=
386-
ν₄_scalar * wdivₕ(gradₕ(ᶜ∇²q_liqʲs.:($$j)))
387-
@. Yₜ.c.sgsʲs.:($$j).q_ice -=
388-
ν₄_scalar * wdivₕ(gradₕ(ᶜ∇²q_iceʲs.:($$j)))
389-
@. Yₜ.c.sgsʲs.:($$j).n_liq -=
390-
ν₄_scalar * wdivₕ(gradₕ(ᶜ∇²n_liqʲs.:($$j)))
391-
@. Yₜ.c.sgsʲs.:($$j).q_rai -=
392-
ν₄_scalar_for_precip * wdivₕ(gradₕ(ᶜ∇²q_raiʲs.:($$j)))
393-
@. Yₜ.c.sgsʲs.:($$j).q_sno -=
394-
ν₄_scalar_for_precip * wdivₕ(gradₕ(ᶜ∇²q_snoʲs.:($$j)))
395-
@. Yₜ.c.sgsʲs.:($$j).n_rai -=
396-
ν₄_scalar_for_precip * wdivₕ(gradₕ(ᶜ∇²n_raiʲs.:($$j)))
397-
end
398-
end
399-
end
400-
return nothing
280+
return (core_dynamics_pairs..., tc_dynamics_pairs...)
401281
end

0 commit comments

Comments
 (0)