Skip to content

Commit 6bfcdb8

Browse files
sajjadazimisajjadazimi
committed
add edmf 2M microphysics sedimentation terms
1 parent d94fd43 commit 6bfcdb8

File tree

4 files changed

+176
-31
lines changed

4 files changed

+176
-31
lines changed

‎src/cache/precipitation_precomputed_quantities.jl

Lines changed: 22 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -107,7 +107,7 @@ function set_precipitation_velocities!(
107107
(; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwnₗ, ᶜwnᵣ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜts, ᶜu) = p.precomputed
108108
(; ᶜΦ) = p.core
109109

110-
cm1c = CAP.microphysics_cloud_params(p.params)
110+
cmc = CAP.microphysics_cloud_params(p.params)
111111
cm1p = CAP.microphysics_1m_params(p.params)
112112
cm2p = CAP.microphysics_2m_params(p.params)
113113
thp = CAP.thermodynamics_params(p.params)
@@ -141,24 +141,30 @@ function set_precipitation_velocities!(
141141
specific(Y.c.ρq_sno, Y.c.ρ),
142142
)
143143
# compute sedimentation velocity for cloud condensate [m/s]
144-
# TODO sedimentation velocities of cloud condensates are based
145-
# on the 1M scheme. Sedimentation velocity of cloud number concentration
146-
# is equal to that of the mass.
147-
@. ᶜwnₗ = CMNe.terminal_velocity(
148-
cm1c.liquid,
149-
cm1c.Ch2022.rain,
150-
Y.c.ρ,
151-
specific(Y.c.ρq_liq, Y.c.ρ),
144+
# TODO sedimentation of ice is based on the 1M scheme
145+
@. ᶜwnₗ = getindex(
146+
CM2.cloud_terminal_velocity(
147+
cm2p.sb.pdf_c,
148+
cm2p.tv,
149+
specific(Y.c.ρq_liq, Y.c.ρ),
150+
Y.c.ρ,
151+
Y.c.ρn_liq,
152+
),
153+
1,
152154
)
153-
@. ᶜwₗ = CMNe.terminal_velocity(
154-
cm1c.liquid,
155-
cm1c.Ch2022.rain,
156-
Y.c.ρ,
157-
specific(Y.c.ρq_liq, Y.c.ρ),
155+
@. ᶜwₗ = getindex(
156+
CM2.cloud_terminal_velocity(
157+
cm2p.sb.pdf_c,
158+
cm2p.tv,
159+
specific(Y.c.ρq_liq, Y.c.ρ),
160+
Y.c.ρ,
161+
Y.c.ρn_liq,
162+
),
163+
2,
158164
)
159165
@. ᶜwᵢ = CMNe.terminal_velocity(
160-
cm1c.ice,
161-
cm1c.Ch2022.small_ice,
166+
cmc.ice,
167+
cmc.Ch2022.small_ice,
162168
Y.c.ρ,
163169
specific(Y.c.ρq_ice, Y.c.ρ),
164170
)

‎src/cache/precomputed_quantities.jl

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -157,6 +157,14 @@ function precomputed_quantities(Y, atmos)
157157
ᶜSqₛᵖʲs = similar(Y.c, NTuple{n, FT}),
158158
ᶜSnₗᵖʲs = similar(Y.c, NTuple{n, FT}),
159159
ᶜSnᵣᵖʲs = similar(Y.c, NTuple{n, FT}),
160+
ᶜwₗʲs = similar(Y.c, NTuple{n, FT}),
161+
ᶜwᵢʲs = similar(Y.c, NTuple{n, FT}),
162+
ᶜwᵣʲs = similar(Y.c, NTuple{n, FT}),
163+
ᶜwₛʲs = similar(Y.c, NTuple{n, FT}),
164+
ᶜwₜʲs = similar(Y.c, NTuple{n, FT}),
165+
ᶜwₕʲs = similar(Y.c, NTuple{n, FT}),
166+
ᶜwₙₗʲs = similar(Y.c, NTuple{n, FT}),
167+
ᶜwₙᵣʲs = similar(Y.c, NTuple{n, FT}),
160168
ᶜSqₗᵖ⁰ = similar(Y.c, FT),
161169
ᶜSqᵢᵖ⁰ = similar(Y.c, FT),
162170
ᶜSqᵣᵖ⁰ = similar(Y.c, FT),

‎src/cache/prognostic_edmf_precomputed_quantities.jl

Lines changed: 96 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -770,7 +770,8 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation
770770

771771
(; params, dt) = p
772772
thp = CAP.thermodynamics_params(params)
773-
cmp = CAP.microphysics_2m_params(params)
773+
cm1p = CAP.microphysics_1m_params(p.params)
774+
cm2p = CAP.microphysics_2m_params(p.params)
774775
cmc = CAP.microphysics_cloud_params(params)
775776

776777
(;
@@ -786,6 +787,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation
786787
) = p.precomputed
787788
(; ᶜSqₗᵖ⁰, ᶜSqᵢᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰, ᶜSnₗᵖ⁰, ᶜSnᵣᵖ⁰, ᶜts⁰, ᶜu⁰) =
788789
p.precomputed
790+
(; ᶜwₗʲs, ᶜwᵢʲs, ᶜwᵣʲs, ᶜwₛʲs, ᶜwₙₗʲs, ᶜwₙᵣʲs, ᶜwₜʲs, ᶜwₕʲs) = p.precomputed
789791

790792
ᶜSᵖ = p.scratch.ᶜtemp_scalar
791793
ᶜS₂ᵖ = p.scratch.ᶜtemp_scalar_2
@@ -801,7 +803,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation
801803
seasalt_mean_radius,
802804
sulfate_num,
803805
p.tracers.prescribed_aerosols_field,
804-
cmp.aerosol,
806+
cm2p.aerosol,
805807
)
806808
else
807809
@. seasalt_num = 0
@@ -811,7 +813,95 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation
811813

812814
# Compute sources
813815
n = n_mass_flux_subdomains(p.atmos.turbconv_model)
816+
FT = eltype(params)
814817
for j in 1:n
818+
819+
# compute terminal velocity for precipitation
820+
# TODO sedimentation of snow is based on the 1M scheme
821+
@. ᶜwₙᵣʲs.:($$j) = getindex(
822+
CM2.rain_terminal_velocity(
823+
cm2p.sb,
824+
cm2p.tv,
825+
max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).q_rai),
826+
ᶜρʲs.:($$j),
827+
max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).n_rai),
828+
),
829+
1,
830+
)
831+
@. ᶜwᵣʲs.:($$j) = getindex(
832+
CM2.rain_terminal_velocity(
833+
cm2p.sb,
834+
cm2p.tv,
835+
max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).q_rai),
836+
ᶜρʲs.:($$j),
837+
max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).n_rai),
838+
),
839+
2,
840+
)
841+
@. ᶜwₛʲs.:($$j) = CM1.terminal_velocity(
842+
cm1p.ps,
843+
cm1p.tv.snow,
844+
ᶜρʲs.:($$j),
845+
max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).q_sno),
846+
)
847+
# compute sedimentation velocity for cloud condensate [m/s]
848+
# TODO sedimentation of ice is based on the 1M scheme
849+
@. ᶜwₙₗʲs.:($$j) = getindex(
850+
CM2.cloud_terminal_velocity(
851+
cm2p.sb.pdf_c,
852+
cm2p.tv,
853+
max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).q_liq),
854+
ᶜρʲs.:($$j),
855+
max(zero(Y.c.ρ), ᶜρʲs.:($$j) * Y.c.sgsʲs.:($$j).n_liq),
856+
),
857+
1,
858+
)
859+
@. ᶜwₗʲs.:($$j) = getindex(
860+
CM2.cloud_terminal_velocity(
861+
cm2p.sb.pdf_c,
862+
cm2p.tv,
863+
max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).q_liq),
864+
ᶜρʲs.:($$j),
865+
max(zero(Y.c.ρ), ᶜρʲs.:($$j) * Y.c.sgsʲs.:($$j).n_liq),
866+
),
867+
2,
868+
)
869+
@. ᶜwᵢʲs.:($$j) = CMNe.terminal_velocity(
870+
cmc.ice,
871+
cmc.Ch2022.small_ice,
872+
ᶜρʲs.:($$j),
873+
max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).q_ice),
874+
)
875+
# compute their contirbutions to energy and total water advection
876+
@. ᶜwₜʲs.:($$j) = ifelse(
877+
Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_tot > FT(0),
878+
(
879+
ᶜwₗʲs.:($$j) * Y.c.sgsʲs.:($$j).q_liq +
880+
ᶜwᵢʲs.:($$j) * Y.c.sgsʲs.:($$j).q_ice +
881+
ᶜwᵣʲs.:($$j) * Y.c.sgsʲs.:($$j).q_rai +
882+
ᶜwₛʲs.:($$j) * Y.c.sgsʲs.:($$j).q_sno
883+
) / Y.c.sgsʲs.:($$j).q_tot,
884+
FT(0),
885+
)
886+
@. ᶜwₕʲs.:($$j) = ifelse(
887+
Y.c.sgsʲs.:($$j).ρa * abs(Y.c.sgsʲs.:($$j).mse) > FT(0),
888+
(
889+
ᶜwₗʲs.:($$j) *
890+
Y.c.sgsʲs.:($$j).q_liq *
891+
(Iₗ(thp, ᶜtsʲs.:($$j)) + ᶜΦ) +
892+
ᶜwᵢʲs.:($$j) *
893+
Y.c.sgsʲs.:($$j).q_ice *
894+
(Iᵢ(thp, ᶜtsʲs.:($$j)) + ᶜΦ) +
895+
ᶜwᵣʲs.:($$j) *
896+
Y.c.sgsʲs.:($$j).q_rai *
897+
(Iₗ(thp, ᶜtsʲs.:($$j)) + ᶜΦ) +
898+
ᶜwₛʲs.:($$j) *
899+
Y.c.sgsʲs.:($$j).q_sno *
900+
(Iᵢ(thp, ᶜtsʲs.:($$j)) + ᶜΦ)
901+
) / abs(Y.c.sgsʲs.:($$j).mse),
902+
FT(0),
903+
)
904+
815905
# Precipitation sources and sinks from the updrafts
816906
compute_warm_precipitation_sources_2M!(
817907
ᶜSᵖ,
@@ -830,7 +920,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation
830920
Y.c.sgsʲs.:($j).q_sno,
831921
ᶜtsʲs.:($j),
832922
dt,
833-
cmp,
923+
cm2p,
834924
thp,
835925
)
836926
@. ᶜSqᵢᵖʲs.:($$j) = 0
@@ -870,7 +960,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation
870960
Y.c.sgsʲs.:($$j).n_liq + Y.c.sgsʲs.:($$j).n_rai,
871961
ᶜρʲs.:($$j),
872962
max(0, w_component.(Geometry.WVector.(ᶜuʲs.:($$j)))),
873-
(cmp,),
963+
(cm2p,),
874964
thp,
875965
ᶜtsʲs.:($$j),
876966
dt,
@@ -903,7 +993,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation
903993
ᶜq_sno⁰,
904994
ᶜts⁰,
905995
dt,
906-
cmp,
996+
cm2p,
907997
thp,
908998
)
909999
@. ᶜSqᵢᵖ⁰ = 0
@@ -943,7 +1033,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation
9431033
ᶜn_liq⁰ + ᶜn_rai⁰,
9441034
ᶜρ⁰,
9451035
w_component.(Geometry.WVector.(ᶜu⁰)),
946-
(cmp,),
1036+
(cm2p,),
9471037
thp,
9481038
ᶜts⁰,
9491039
dt,

‎src/prognostic_equations/advection.jl

Lines changed: 50 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -404,8 +404,10 @@ function edmfx_sgs_vertical_advection_tendency!(
404404
@. ᶜχʲₜ += va
405405
end
406406

407-
if p.atmos.moisture_model isa NonEquilMoistModel &&
408-
p.atmos.microphysics_model isa Microphysics1Moment
407+
if p.atmos.moisture_model isa NonEquilMoistModel && (
408+
p.atmos.microphysics_model isa Microphysics1Moment ||
409+
p.atmos.microphysics_model isa Microphysics2Moment
410+
)
409411
# TODO - add contibutions to sgs mass flux from tracer sedimentation
410412
# TODO - add precipitation and cloud sedimentation in implicit solver/tendency with if/else
411413

@@ -506,16 +508,55 @@ function edmfx_sgs_vertical_advection_tendency!(
506508
end
507509
if p.atmos.moisture_model isa NonEquilMoistModel &&
508510
p.atmos.microphysics_model isa Microphysics2Moment
509-
# TODO add sedimentation for number concentrations
510511

511-
sgs_microphysics_tracers =
512-
(:q_liq, :q_ice, :q_rai, :q_sno, :n_liq, :n_rai)
512+
FT = eltype(params)
513+
514+
# Sedimentation velocities for microphysics number concentrations
515+
# (or any tracers that does not directly participate in variations of q_tot and mse)
516+
sgs_microphysics_tracers = (
517+
(@name(c.sgsʲs.:(1).n_liq), @name(ᶜwₙₗʲs.:(1))),
518+
(@name(c.sgsʲs.:(1).n_rai), @name(ᶜwₙᵣʲs.:(1))),
519+
)
520+
521+
MatrixFields.unrolled_foreach(
522+
sgs_microphysics_tracers,
523+
) do (χʲ_name, χʲ_name)
524+
MatrixFields.has_field(Y, χʲ_name) || return
513525

514-
for χ_name in sgs_microphysics_tracers
515-
ᶜχʲ = getproperty(Y.c.sgsʲs.:($j), χ_name)
516-
va = vertical_advection(ᶠu³ʲs.:($j), ᶜχʲ, edmfx_upwinding)
517-
ᶜχʲₜ = getproperty(Yₜ.c.sgsʲs.:($j), χ_name)
526+
ᶜχʲ = MatrixFields.get_field(Y, χʲ_name)
527+
ᶜχʲₜ = MatrixFields.get_field(Yₜ, χʲ_name)
528+
ᶜwʲ = MatrixFields.get_field(p.precomputed, wʲ_name)
529+
ᶠw³ʲ = (@. lazy(CT3(ᶠinterp(Geometry.WVector(-1 * ᶜwʲ)))))
530+
ᶜw³ʲ = (@. lazy(CT3(Geometry.WVector(-1 * ᶜwʲ))))
531+
532+
# Advective form advection of moisture tracers with the grid mean velocity
533+
va = vertical_advection(ᶠu³ʲs.:($j), ᶜχʲ, edmf_upwnd)
534+
@. ᶜχʲₜ += va
535+
536+
# Advective form advection of moisture tracers with sedimentation velocities
537+
va = vertical_advection(ᶠw³ʲ, ᶜχʲ, edmf_upwnd)
518538
@. ᶜχʲₜ += va
539+
540+
# moisture tracers terms proportional to 1/ρ̂ ∂zρ̂
541+
ᶜinv_ρ̂ = (@. lazy(
542+
specific(
543+
FT(1),
544+
Y.c.sgsʲs.:($$j).ρa,
545+
FT(0),
546+
Y.c.ρ,
547+
turbconv_model,
548+
),
549+
))
550+
ᶜ∂ρ̂∂z = (@. lazy(
551+
upwind_biased_grad(
552+
-1 * Geometry.WVector(ᶜwʲ),
553+
Y.c.sgsʲs.:($$j).ρa,
554+
),
555+
))
556+
@. ᶜχʲₜ -= dot(ᶜinv_ρ̂ * ᶜ∂ρ̂∂z, ᶜw³ʲ) * ᶜχʲ
557+
558+
# moisture tracer term proportional to velocity gradients
559+
@. ᶜχʲₜ -= ᶜdivᵥ(ᶠw³ʲ) * ᶜχʲ
519560
end
520561

521562
end

0 commit comments

Comments
 (0)