Skip to content

Commit f496d5c

Browse files
sajjadazimisajjadazimi
authored
Merge pull request #3947 from CliMA/aj/edmf_sgs_flux_1M
Redefine sedimentation velocity in Prognostic EDMF (1M/2M)
2 parents 20ef522 + bfdebde commit f496d5c

File tree

12 files changed

+500
-132
lines changed

12 files changed

+500
-132
lines changed

‎reproducibility_tests/ref_counter.jl‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
256
1+
257
22

33
# **README**
44
#
@@ -20,6 +20,10 @@
2020

2121

2222
#=
23+
257
24+
- Redefine sedimentation velocity for Prognostic EDMF with 1-moment or 2-moment microphysics
25+
on the grid scale; fix a bug in EDMFx SGS mass flux.
26+
2327
256
2428
- Fix a bug in EDMF diffusive flux
2529

‎src/cache/precipitation_precomputed_quantities.jl‎

Lines changed: 310 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -31,12 +31,12 @@ function Kin(ᶜw_precip, ᶜu_air)
3131
end
3232

3333
"""
34-
set_precipitation_velocities!(Y, p, moisture_model, microphysics_model)
34+
set_precipitation_velocities!(Y, p, moisture_model, microphysics_model, turbconv_model)
3535
3636
Updates the precipitation terminal velocity, cloud sedimentation velocity,
3737
and their contribution to total water and energy advection.
3838
"""
39-
function set_precipitation_velocities!(Y, p, _, _)
39+
function set_precipitation_velocities!(Y, p, _, _, _)
4040
(; ᶜwₜqₜ, ᶜwₕhₜ) = p.precomputed
4141
@. ᶜwₜqₜ = Geometry.WVector(0)
4242
@. ᶜwₕhₜ = Geometry.WVector(0)
@@ -47,6 +47,7 @@ function set_precipitation_velocities!(
4747
p,
4848
moisture_model::NonEquilMoistModel,
4949
microphysics_model::Microphysics1Moment,
50+
_,
5051
)
5152
(; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜts, ᶜu) = p.precomputed
5253
(; ᶜΦ) = p.core
@@ -98,13 +99,140 @@ function set_precipitation_velocities!(
9899
) / Y.c.ρ
99100
return nothing
100101
end
102+
function set_precipitation_velocities!(
103+
Y,
104+
p,
105+
moisture_model::NonEquilMoistModel,
106+
microphysics_model::Microphysics1Moment,
107+
turbconv_model::PrognosticEDMFX,
108+
)
109+
(; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwₜqₜ, ᶜwₕhₜ) = p.precomputed
110+
(; ᶜΦ) = p.core
111+
(; ᶜwₗʲs, ᶜwᵢʲs, ᶜwᵣʲs, ᶜwₛʲs, ᶜwₕʲs) = p.precomputed
112+
(; ᶜts⁰) = p.precomputed
113+
cmc = CAP.microphysics_cloud_params(p.params)
114+
cmp = CAP.microphysics_1m_params(p.params)
115+
thp = CAP.thermodynamics_params(p.params)
116+
FT = eltype(p.params)
117+
118+
ᶜρ⁰ = @. lazy(TD.air_density(thp, ᶜts⁰))
119+
ᶜρa⁰ = @. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, turbconv_model))
120+
n = n_mass_flux_subdomains(turbconv_model)
121+
122+
# scratch to compute env velocities
123+
ᶜw⁰ = p.scratch.ᶜtemp_scalar
124+
# scratch to add positive masses of subdomains
125+
# TODO use Y.c.ρq instead of ᶜρχ for computing gs velocities by averaging velocities
126+
# over subdomains once negative subdomain mass issues are resolved
127+
ᶜρχ = p.scratch.ᶜtemp_scalar_2
128+
# scratch for adding energy fluxes over subdomains
129+
ᶜρwₕhₜ = p.scratch.ᶜtemp_scalar_3
130+
131+
# Compute gs sedimentation velocity based on subdomain velocities (assuming gs flux
132+
# equals sum of sgs fluxes)
133+
134+
ᶜq_liq⁰ = ᶜspecific_env_value(@name(q_liq), Y, p)
135+
ᶜq_ice⁰ = ᶜspecific_env_value(@name(q_ice), Y, p)
136+
ᶜq_rai⁰ = ᶜspecific_env_value(@name(q_rai), Y, p)
137+
ᶜq_sno⁰ = ᶜspecific_env_value(@name(q_sno), Y, p)
138+
139+
# Cloud liquid
140+
@. ᶜw⁰ = CMNe.terminal_velocity(
141+
cmc.liquid,
142+
cmc.Ch2022.rain,
143+
ᶜρ⁰,
144+
max(zero(Y.c.ρ), ᶜq_liq⁰),
145+
)
146+
@. ᶜwₗ = ᶜρa⁰ * ᶜq_liq⁰ * ᶜw⁰
147+
@. ᶜρχ = max(zero(Y.c.ρ), ᶜρa⁰ * ᶜq_liq⁰)
148+
for j in 1:n
149+
@. ᶜwₗ += Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_liq * ᶜwₗʲs.:($$j)
150+
@. ᶜρχ += max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_liq)
151+
end
152+
@. ᶜwₗ = ifelse(ᶜρχ > FT(0), ᶜwₗ / ᶜρχ, FT(0))
153+
154+
# contribution of env cloud liquid advection to htot advection
155+
@. ᶜρwₕhₜ = ᶜρa⁰ * ᶜq_liq⁰ * ᶜw⁰ * (Iₗ(thp, ᶜts⁰) + ᶜΦ)
156+
157+
# Cloud ice
158+
@. ᶜw⁰ = CMNe.terminal_velocity(
159+
cmc.ice,
160+
cmc.Ch2022.small_ice,
161+
ᶜρ⁰,
162+
max(zero(Y.c.ρ), ᶜq_ice⁰),
163+
)
164+
@. ᶜwᵢ = ᶜρa⁰ * ᶜq_ice⁰ * ᶜw⁰
165+
@. ᶜρχ = max(zero(Y.c.ρ), ᶜρa⁰ * ᶜq_ice⁰)
166+
for j in 1:n
167+
@. ᶜwᵢ += Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_ice * ᶜwᵢʲs.:($$j)
168+
@. ᶜρχ += max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_ice)
169+
end
170+
@. ᶜwᵢ = ifelse(ᶜρχ > FT(0), ᶜwᵢ / ᶜρχ, FT(0))
171+
172+
# contribution of env cloud ice advection to htot advection
173+
@. ᶜρwₕhₜ += ᶜρa⁰ * ᶜq_ice⁰ * ᶜw⁰ * (Iᵢ(thp, ᶜts⁰) + ᶜΦ)
174+
175+
# Rain
176+
@. ᶜw⁰ = CM1.terminal_velocity(
177+
cmp.pr,
178+
cmp.tv.rain,
179+
ᶜρ⁰,
180+
max(zero(Y.c.ρ), ᶜq_rai⁰),
181+
)
182+
@. ᶜwᵣ = ᶜρa⁰ * ᶜq_rai⁰ * ᶜw⁰
183+
@. ᶜρχ = max(zero(Y.c.ρ), ᶜρa⁰ * ᶜq_rai⁰)
184+
for j in 1:n
185+
@. ᶜwᵣ += Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_rai * ᶜwᵣʲs.:($$j)
186+
@. ᶜρχ += max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_rai)
187+
end
188+
@. ᶜwᵣ = ifelse(ᶜρχ > FT(0), ᶜwᵣ / ᶜρχ, FT(0))
189+
190+
# contribution of env rain advection to htot advection
191+
@. ᶜρwₕhₜ += ᶜρa⁰ * ᶜq_rai⁰ * ᶜw⁰ * (Iₗ(thp, ᶜts⁰) + ᶜΦ)
192+
193+
# Snow
194+
@. ᶜw⁰ = CM1.terminal_velocity(
195+
cmp.ps,
196+
cmp.tv.snow,
197+
ᶜρ⁰,
198+
max(zero(Y.c.ρ), ᶜq_sno⁰),
199+
)
200+
@. ᶜwₛ = ᶜρa⁰ * ᶜq_sno⁰ * ᶜw⁰
201+
@. ᶜρχ = max(zero(Y.c.ρ), ᶜρa⁰ * ᶜq_sno⁰)
202+
for j in 1:n
203+
@. ᶜwₛ += Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_sno * ᶜwₛʲs.:($$j)
204+
@. ᶜρχ += max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_sno)
205+
end
206+
@. ᶜwₛ = ifelse(ᶜρχ > FT(0), ᶜwₛ / ᶜρχ, FT(0))
207+
208+
# contribution of env snow advection to htot advection
209+
@. ᶜρwₕhₜ += ᶜρa⁰ * ᶜq_sno⁰ * ᶜw⁰ * (Iᵢ(thp, ᶜts⁰) + ᶜΦ)
210+
211+
# Add contributions to energy and total water advection
212+
# TODO do we need to add kinetic energy of subdomains?
213+
for j in 1:n
214+
@. ᶜρwₕhₜ += Y.c.sgsʲs.:($$j).ρa * ᶜwₕʲs.:($$j) * Y.c.sgsʲs.:($$j).mse
215+
end
216+
@. ᶜwₕhₜ = Geometry.WVector(ᶜρwₕhₜ) / Y.c.ρ
217+
218+
@. ᶜwₜqₜ =
219+
Geometry.WVector(
220+
ᶜwₗ * Y.c.ρq_liq +
221+
ᶜwᵢ * Y.c.ρq_ice +
222+
ᶜwᵣ * Y.c.ρq_rai +
223+
ᶜwₛ * Y.c.ρq_sno,
224+
) / Y.c.ρ
225+
226+
return nothing
227+
end
101228
function set_precipitation_velocities!(
102229
Y,
103230
p,
104231
moisture_model::NonEquilMoistModel,
105232
microphysics_model::Microphysics2Moment,
233+
_,
106234
)
107-
(; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwnₗ, ᶜwnᵣ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜts, ᶜu) = p.precomputed
235+
(; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwₙₗ, ᶜwₙᵣ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜts, ᶜu) = p.precomputed
108236
(; ᶜΦ) = p.core
109237

110238
cmc = CAP.microphysics_cloud_params(p.params)
@@ -114,7 +242,7 @@ function set_precipitation_velocities!(
114242

115243
# compute the precipitation terminal velocity [m/s]
116244
# TODO sedimentation of snow is based on the 1M scheme
117-
@. ᶜwnᵣ = getindex(
245+
@. ᶜwₙᵣ = getindex(
118246
CM2.rain_terminal_velocity(
119247
cm2p.sb,
120248
cm2p.rtv,
@@ -142,7 +270,7 @@ function set_precipitation_velocities!(
142270
)
143271
# compute sedimentation velocity for cloud condensate [m/s]
144272
# TODO sedimentation of ice is based on the 1M scheme
145-
@. ᶜwnₗ = getindex(
273+
@. ᶜwₙₗ = getindex(
146274
CM2.cloud_terminal_velocity(
147275
cm2p.sb.pdf_c,
148276
cm2p.ctv,
@@ -186,6 +314,183 @@ function set_precipitation_velocities!(
186314
) / Y.c.ρ
187315
return nothing
188316
end
317+
function set_precipitation_velocities!(
318+
Y,
319+
p,
320+
moisture_model::NonEquilMoistModel,
321+
microphysics_model::Microphysics2Moment,
322+
turbconv_model::PrognosticEDMFX,
323+
)
324+
(; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwₙₗ, ᶜwₙᵣ, ᶜwₜqₜ, ᶜwₕhₜ) = p.precomputed
325+
(; ᶜΦ) = p.core
326+
(; ᶜwₗʲs, ᶜwᵢʲs, ᶜwᵣʲs, ᶜwₛʲs, ᶜwₙₗʲs, ᶜwₙᵣʲs, ᶜwₕʲs) = p.precomputed
327+
(; ᶜts⁰) = p.precomputed
328+
cmc = CAP.microphysics_cloud_params(p.params)
329+
cm1p = CAP.microphysics_1m_params(p.params)
330+
cm2p = CAP.microphysics_2m_params(p.params)
331+
thp = CAP.thermodynamics_params(p.params)
332+
FT = eltype(p.params)
333+
334+
ᶜρ⁰ = @. lazy(TD.air_density(thp, ᶜts⁰))
335+
ᶜρa⁰ = @. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, turbconv_model))
336+
n = n_mass_flux_subdomains(turbconv_model)
337+
338+
# scratch to compute env velocities
339+
ᶜw⁰ = p.scratch.ᶜtemp_scalar
340+
# scratch to add positive masses of subdomains
341+
# TODO use Y.c.ρq instead of ᶜρχ for computing gs velocities by averaging velocities
342+
# over subdomains once negative subdomain mass issues are resolved
343+
ᶜρχ = p.scratch.ᶜtemp_scalar_2
344+
# scratch for adding energy fluxes over subdomains
345+
ᶜρwₕhₜ = p.scratch.ᶜtemp_scalar_3
346+
347+
# Compute gs sedimentation velocity based on subdomain velocities (assuming gs flux
348+
# equals sum of sgs fluxes)
349+
350+
ᶜq_liq⁰ = ᶜspecific_env_value(@name(q_liq), Y, p)
351+
ᶜq_ice⁰ = ᶜspecific_env_value(@name(q_ice), Y, p)
352+
ᶜq_rai⁰ = ᶜspecific_env_value(@name(q_rai), Y, p)
353+
ᶜq_sno⁰ = ᶜspecific_env_value(@name(q_sno), Y, p)
354+
ᶜn_liq⁰ = ᶜspecific_env_value(@name(n_liq), Y, p)
355+
ᶜn_rai⁰ = ᶜspecific_env_value(@name(n_rai), Y, p)
356+
357+
# Cloud liquid (number)
358+
@. ᶜw⁰ = getindex(
359+
CM2.rain_terminal_velocity(
360+
cm2p.sb,
361+
cm2p.rtv,
362+
max(zero(Y.c.ρ), ᶜq_liq⁰),
363+
ᶜρ⁰,
364+
max(zero(Y.c.ρ), ᶜn_liq⁰),
365+
),
366+
1,
367+
)
368+
@. ᶜwₙₗ = ᶜρa⁰ * ᶜn_liq⁰ * ᶜw⁰
369+
@. ᶜρχ = max(zero(Y.c.ρ), ᶜρa⁰ * ᶜn_liq⁰)
370+
for j in 1:n
371+
@. ᶜwₙₗ += Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).n_liq * ᶜwₙₗʲs.:($$j)
372+
@. ᶜρχ += max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).n_liq)
373+
end
374+
@. ᶜwₙₗ = ifelse(ᶜρχ > FT(0), ᶜwₙₗ / ᶜρχ, FT(0))
375+
376+
# Cloud liquid (mass)
377+
@. ᶜw⁰ = getindex(
378+
CM2.rain_terminal_velocity(
379+
cm2p.sb,
380+
cm2p.rtv,
381+
max(zero(Y.c.ρ), ᶜq_liq⁰),
382+
ᶜρ⁰,
383+
max(zero(Y.c.ρ), ᶜn_liq⁰),
384+
),
385+
2,
386+
)
387+
@. ᶜwₗ = ᶜρa⁰ * ᶜq_liq⁰ * ᶜw⁰
388+
@. ᶜρχ = max(zero(Y.c.ρ), ᶜρa⁰ * ᶜq_liq⁰)
389+
for j in 1:n
390+
@. ᶜwₗ += Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_liq * ᶜwₗʲs.:($$j)
391+
@. ᶜρχ += max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_liq)
392+
end
393+
@. ᶜwₗ = ifelse(ᶜρχ > FT(0), ᶜwₗ / ᶜρχ, FT(0))
394+
395+
# contribution of env cloud liquid advection to htot advection
396+
@. ᶜρwₕhₜ = ᶜρa⁰ * ᶜq_liq⁰ * ᶜw⁰ * (Iₗ(thp, ᶜts⁰) + ᶜΦ)
397+
398+
# Cloud ice
399+
# TODO sedimentation of ice is based on the 1M scheme
400+
@. ᶜw⁰ = CMNe.terminal_velocity(
401+
cmc.ice,
402+
cmc.Ch2022.small_ice,
403+
ᶜρ⁰,
404+
max(zero(Y.c.ρ), ᶜq_ice⁰),
405+
)
406+
@. ᶜwᵢ = ᶜρa⁰ * ᶜq_ice⁰ * ᶜw⁰
407+
@. ᶜρχ = max(zero(Y.c.ρ), ᶜρa⁰ * ᶜq_ice⁰)
408+
for j in 1:n
409+
@. ᶜwᵢ += Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_ice * ᶜwᵢʲs.:($$j)
410+
@. ᶜρχ += max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_ice)
411+
end
412+
@. ᶜwᵢ = ifelse(ᶜρχ > FT(0), ᶜwᵢ / ᶜρχ, FT(0))
413+
414+
# contribution of env cloud ice advection to htot advection
415+
@. ᶜρwₕhₜ += ᶜρa⁰ * ᶜq_ice⁰ * ᶜw⁰ * (Iᵢ(thp, ᶜts⁰) + ᶜΦ)
416+
417+
# Rain (number)
418+
@. ᶜw⁰ = getindex(
419+
CM2.rain_terminal_velocity(
420+
cm2p.sb,
421+
cm2p.rtv,
422+
max(zero(Y.c.ρ), ᶜq_rai⁰),
423+
ᶜρ⁰,
424+
max(zero(Y.c.ρ), ᶜn_rai⁰),
425+
),
426+
1,
427+
)
428+
@. ᶜwₙᵣ = ᶜρa⁰ * ᶜn_rai⁰ * ᶜw⁰
429+
@. ᶜρχ = max(zero(Y.c.ρ), ᶜρa⁰ * ᶜn_rai⁰)
430+
for j in 1:n
431+
@. ᶜwₙᵣ += Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).n_rai * ᶜwₙᵣʲs.:($$j)
432+
@. ᶜρχ += max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).n_rai)
433+
end
434+
@. ᶜwₙᵣ = ifelse(ᶜρχ > FT(0), ᶜwₙᵣ / ᶜρχ, FT(0))
435+
436+
# Rain (mass)
437+
@. ᶜw⁰ = getindex(
438+
CM2.rain_terminal_velocity(
439+
cm2p.sb,
440+
cm2p.rtv,
441+
max(zero(Y.c.ρ), ᶜq_rai⁰),
442+
ᶜρ⁰,
443+
max(zero(Y.c.ρ), ᶜn_rai⁰),
444+
),
445+
2,
446+
)
447+
@. ᶜwᵣ = ᶜρa⁰ * ᶜq_rai⁰ * ᶜw⁰
448+
@. ᶜρχ = max(zero(Y.c.ρ), ᶜρa⁰ * ᶜq_rai⁰)
449+
for j in 1:n
450+
@. ᶜwᵣ += Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_rai * ᶜwᵣʲs.:($$j)
451+
@. ᶜρχ += max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_rai)
452+
end
453+
@. ᶜwᵣ = ifelse(ᶜρχ > FT(0), ᶜwᵣ / ᶜρχ, FT(0))
454+
455+
# contribution of env rain advection to qtot advection
456+
@. ᶜρwₕhₜ += ᶜρa⁰ * ᶜq_rai⁰ * ᶜw⁰ * (Iₗ(thp, ᶜts⁰) + ᶜΦ)
457+
458+
# Snow
459+
# TODO sedimentation of snow is based on the 1M scheme
460+
@. ᶜw⁰ = CM1.terminal_velocity(
461+
cm1p.ps,
462+
cm1p.tv.snow,
463+
ᶜρ⁰,
464+
max(zero(Y.c.ρ), ᶜq_sno⁰),
465+
)
466+
@. ᶜwₛ = ᶜρa⁰ * ᶜq_sno⁰ * ᶜw⁰
467+
@. ᶜρχ = max(zero(Y.c.ρ), ᶜρa⁰ * ᶜq_sno⁰)
468+
for j in 1:n
469+
@. ᶜwₛ += Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_sno * ᶜwₛʲs.:($$j)
470+
@. ᶜρχ += max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_sno)
471+
end
472+
@. ᶜwₛ = ifelse(ᶜρχ > FT(0), ᶜwₛ / ᶜρχ, FT(0))
473+
474+
# contribution of env snow advection to htot advection
475+
@. ᶜρwₕhₜ += ᶜρa⁰ * ᶜq_sno⁰ * ᶜw⁰ * (Iᵢ(thp, ᶜts⁰) + ᶜΦ)
476+
477+
# Add contributions to energy and total water advection
478+
# TODO do we need to add kinetic energy of subdomains?
479+
for j in 1:n
480+
@. ᶜρwₕhₜ += Y.c.sgsʲs.:($$j).ρa * ᶜwₕʲs.:($$j) * Y.c.sgsʲs.:($$j).mse
481+
end
482+
@. ᶜwₕhₜ = Geometry.WVector(ᶜρwₕhₜ) / Y.c.ρ
483+
484+
@. ᶜwₜqₜ =
485+
Geometry.WVector(
486+
ᶜwₗ * Y.c.ρq_liq +
487+
ᶜwᵢ * Y.c.ρq_ice +
488+
ᶜwᵣ * Y.c.ρq_rai +
489+
ᶜwₛ * Y.c.ρq_sno,
490+
) / Y.c.ρ
491+
492+
return nothing
493+
end
189494

190495
"""
191496
set_precipitation_cache!(Y, p, microphysics_model, turbconv_model)

0 commit comments

Comments
 (0)