Skip to content

Commit

Permalink
Merge #514
Browse files Browse the repository at this point in the history
514: add snow advection r=charleskawczynski a=trontrytel

Another peel off from #506 . This one adds snow advection. 

Does not change RiCO.

(looks like the change is somewhere in updraft tendencies then)

Co-authored-by: Anna Jaruga <[email protected]>
  • Loading branch information
bors[bot] and trontrytel authored Nov 4, 2021
2 parents 4c2a068 + 43d22f0 commit d52b1e6
Show file tree
Hide file tree
Showing 6 changed files with 39 additions and 13 deletions.
8 changes: 6 additions & 2 deletions integration_tests/utils/main.jl
Original file line number Diff line number Diff line change
Expand Up @@ -117,13 +117,15 @@ cent_aux_vars_edmf(FT, n_up) = (;
qt_tendency_rain_evap = FT(0),
qr_tendency_rain_evap = FT(0),
qr_tendency_advection = FT(0),
qs_tendency_advection = FT(0),
en_2m = (;
tke = cent_aux_vars_en_2m(FT),
Hvar = cent_aux_vars_en_2m(FT),
QTvar = cent_aux_vars_en_2m(FT),
HQTcov = cent_aux_vars_en_2m(FT),
),
term_vel = FT(0),
term_vel_rain = FT(0),
term_vel_snow = FT(0),
KM = FT(0),
KH = FT(0),
mixing_length = FT(0),
Expand Down Expand Up @@ -184,7 +186,9 @@ cent_prognostic_vars_up(FT) = (; ρarea = FT(0), ρaθ_liq_ice = FT(0), ρaq_tot
cent_prognostic_vars_en(FT) = (; ρatke = FT(0), ρaHvar = FT(0), ρaQTvar = FT(0), ρaHQTcov = FT(0))
cent_prognostic_vars_edmf(FT, n_up) = (;
turbconv = (;
en = cent_prognostic_vars_en(FT), up = ntuple(i -> cent_prognostic_vars_up(FT), n_up), pr = (; qr = FT(0)),
en = cent_prognostic_vars_en(FT),
up = ntuple(i -> cent_prognostic_vars_up(FT), n_up),
pr = (; qr = FT(0), qs = FT(0)),
),
)
# cent_prognostic_vars_edmf(FT, n_up) = (;) # could also use this for empty model
Expand Down
3 changes: 3 additions & 0 deletions src/EDMF_Environment.jl
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@ function sgs_mean(en_thermo::EnvironmentThermodynamics, grid, state, en, precip,
param_set,
precip.precipitation_model,
prog_pr.qr[k],
prog_pr.qs[k],
aux_en.area[k],
ρ0_c[k],
dt,
Expand Down Expand Up @@ -171,6 +172,7 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, p
param_set,
precip.precipitation_model,
prog_pr.qr[k],
prog_pr.qs[k],
aux_en.area[k],
ρ0_c[k],
dt,
Expand Down Expand Up @@ -242,6 +244,7 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, p
param_set,
precip.precipitation_model,
prog_pr.qr[k],
prog_pr.qs[k],
aux_en.area[k],
ρ0_c[k],
dt,
Expand Down
33 changes: 25 additions & 8 deletions src/EDMF_Precipitation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,34 +15,51 @@ function compute_rain_advection_tendencies(::PrecipPhysics, grid, state, gm, TS:
# helper to calculate the rain velocity
# TODO: assuming gm.W = 0
# TODO: verify translation
term_vel = aux_tc.term_vel
term_vel_rain = aux_tc.term_vel_rain
term_vel_snow = aux_tc.term_vel_snow

@inbounds for k in real_center_indices(grid)
term_vel[k] = CM1.terminal_velocity(param_set, rain_type, ρ_0_c[k], prog_pr.qr[k])
term_vel_rain[k] = CM1.terminal_velocity(param_set, rain_type, ρ_0_c[k], prog_pr.qr[k])
term_vel_snow[k] = CM1.terminal_velocity(param_set, snow_type, ρ_0_c[k], prog_pr.qs[k])
end

@inbounds for k in reverse(real_center_indices(grid))
# check stability criterion
CFL_out = Δt / Δz * term_vel[k]
CFL_out_rain = Δt / Δz * term_vel_rain[k]
CFL_out_snow = Δt / Δz * term_vel_snow[k]
if is_toa_center(grid, k)
CFL_in = 0.0
CFL_in_rain = 0.0
CFL_in_snow = 0.0
else
CFL_in = Δt / Δz * term_vel[k + 1]
CFL_in_rain = Δt / Δz * term_vel_rain[k + 1]
CFL_in_snow = Δt / Δz * term_vel_snow[k + 1]
end
if max(CFL_in, CFL_out) > CFL_limit
if max(CFL_in_rain, CFL_in_snow, CFL_out_rain, CFL_out_snow) > CFL_limit
error("Time step is too large for rain fall velocity!")
end

ρ_0_cut = ccut_downwind(ρ_0_c, grid, k)

QR_cut = ccut_downwind(prog_pr.qr, grid, k)
w_cut = ccut_downwind(term_vel, grid, k)
ρQRw_cut = ρ_0_cut .* QR_cut .* w_cut
QS_cut = ccut_downwind(prog_pr.qs, grid, k)

w_cut_rain = ccut_downwind(term_vel_rain, grid, k)
w_cut_snow = ccut_downwind(term_vel_snow, grid, k)

ρQRw_cut = ρ_0_cut .* QR_cut .* w_cut_rain
ρQSw_cut = ρ_0_cut .* QS_cut .* w_cut_snow

∇ρQRw = c∇_downwind(ρQRw_cut, grid, k; bottom = FreeBoundary(), top = SetValue(0))
∇ρQSw = c∇_downwind(ρQSw_cut, grid, k; bottom = FreeBoundary(), top = SetValue(0))

ρ_0_c_k = ρ_0_c[k]

# TODO - some positivity limiters are needed
aux_tc.qr_tendency_advection[k] = ∇ρQRw / ρ_0_c_k
aux_tc.qs_tendency_advection[k] = ∇ρQSw / ρ_0_c_k

tendencies_pr.qr[k] += aux_tc.qr_tendency_advection[k]
tendencies_pr.qs[k] += aux_tc.qs_tendency_advection[k]
end
return
end
Expand Down
1 change: 1 addition & 0 deletions src/EDMF_Updrafts.jl
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@ function compute_precipitation_formation_tendencies(
param_set,
precip.precipitation_model,
prog_pr.qr[k],
prog_pr.qs[k],
aux_up[i].area[k],
ρ0_c[k],
dt,
Expand Down
6 changes: 3 additions & 3 deletions src/microphysics_coupling.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,7 @@ end
"""
Computes the tendencies to qt and θ_liq_ice due to precipitation formation
"""
function precipitation_formation(param_set::APS, precipitation_model, qr, area, ρ0, dt, ts)

function precipitation_formation(param_set::APS, precipitation_model, qr, qs, area, ρ0, dt, ts)

qr_tendency = 0.0

Expand All @@ -36,7 +35,8 @@ function precipitation_formation(param_set::APS, precipitation_model, qr, area,
end

qt_tendency = -qr_tendency
qs_tendency = 0.0
θ_liq_ice_tendency = θ_liq_ice_helper(ts, qt_tendency)

return PrecipFormation(θ_liq_ice_tendency, qt_tendency, qr_tendency)
return PrecipFormation(θ_liq_ice_tendency, qt_tendency, qr_tendency, qs_tendency)
end
1 change: 1 addition & 0 deletions src/types.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ Base.@kwdef struct PrecipFormation{FT}
θ_liq_ice_tendency::FT
qt_tendency::FT
qr_tendency::FT
qs_tendency::FT
end

"""
Expand Down

0 comments on commit d52b1e6

Please sign in to comment.