Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Make several variables ClimaCore fields #418

Merged
merged 2 commits into from
Oct 17, 2021
Merged

Make several variables ClimaCore fields #418

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ce8b177
Make some edmf variables fields
charleskawczynski Oct 17, 2021
82f577b
Make 2m env variables fields
charleskawczynski Oct 17, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 21 additions & 2 deletions integration_tests/utils/main.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,18 +50,37 @@ cent_aux_vars_gm(FT) = (;
W_third_m = FT(0),
QT_third_m = FT(0),
)
cent_aux_vars_en_2m(FT) = (;
dissipation = FT(0),
shear = FT(0),
entr_gain = FT(0),
detr_loss = FT(0),
press = FT(0),
buoy = FT(0),
interdomain = FT(0),
rain_src = FT(0),
)
cent_aux_vars_up(FT) = (; q_liq = FT(0), T = FT(0), RH = FT(0), buoy = FT(0))
cent_aux_vars_edmf(FT, n_up) = (;
turbconv = (;
bulk = (; area = FT(0), θ_liq_ice = FT(0), RH = FT(0), buoy = FT(0), q_tot = FT(0), q_liq = FT(0), T = FT(0)),
up = ntuple(i -> cent_aux_vars_up(FT), n_up),
en = (; area = 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),
),
KM = FT(0),
KH = FT(0),
),
)
cent_aux_vars(FT, n_up) = (; aux_vars_ref_state(FT)..., cent_aux_vars_gm(FT)..., cent_aux_vars_edmf(FT, n_up)...)

# Face only
face_aux_vars_gm(FT) = ()
face_aux_vars_edmf(FT, n_up) = (; turbconv = (; bulk = (; w = FT(0))))
face_aux_vars_edmf(FT, n_up) = (; turbconv = (; bulk = (; w = FT(0)), ρ_ae_KM = FT(0), ρ_ae_KH = FT(0)))
face_aux_vars(FT, n_up) = (; aux_vars_ref_state(FT)..., face_aux_vars_gm(FT)..., face_aux_vars_edmf(FT, n_up)...)

#####
Expand All @@ -86,7 +105,7 @@ face_diagnostic_vars(FT, n_up) = (; face_diagnostic_vars_gm(FT)..., face_diagnos
cent_prognostic_vars(FT, n_up) = (; cent_prognostic_vars_gm(FT)..., cent_prognostic_vars_edmf(FT, n_up)...)
cent_prognostic_vars_gm(FT) = (; u = FT(0), v = FT(0), θ_liq_ice = FT(0), q_tot = FT(0))
cent_prognostic_vars_up(FT) = (; area = FT(0), θ_liq_ice = FT(0), q_tot = FT(0))
cent_prognostic_vars_en(FT) = ()
cent_prognostic_vars_en(FT) = (; tke = FT(0), Hvar = FT(0), QTvar = FT(0), HQTcov = 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)))
# cent_prognostic_vars_edmf(FT, n_up) = (;) # could also use this for empty model
Expand Down
39 changes: 15 additions & 24 deletions src/EDMF_Environment.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,10 +6,6 @@ function initialize_io(en::EnvironmentVariables, Stats::NetCDFIO_Stats)
add_profile(Stats, "env_temperature")
add_profile(Stats, "env_RH")
add_profile(Stats, "env_thetal")
add_profile(Stats, "env_tke")
add_profile(Stats, "env_Hvar")
add_profile(Stats, "env_QTvar")
add_profile(Stats, "env_HQTcov")

add_profile(Stats, "env_cloud_fraction")

Expand All @@ -29,10 +25,6 @@ function io(en::EnvironmentVariables, grid, state, Stats::NetCDFIO_Stats)
write_profile(Stats, "env_temperature", en.T.values)
write_profile(Stats, "env_RH", en.RH.values)
write_profile(Stats, "env_thetal", en.H.values)
write_profile(Stats, "env_tke", en.TKE.values)
write_profile(Stats, "env_Hvar", en.Hvar.values)
write_profile(Stats, "env_QTvar", en.QTvar.values)
write_profile(Stats, "env_HQTcov", en.HQTcov.values)

write_profile(Stats, "env_cloud_fraction", en.cloud_fraction.values)

Expand Down Expand Up @@ -111,6 +103,7 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, r
a, w = FastGaussQuadrature.gausshermite(en_thermo.quadrature_order)
p0_c = center_ref_state(state).p0
ρ0_c = center_ref_state(state).ρ0
prog_en = center_prog_environment(state)

#TODO - remember you output source terms multipierd by dt (bec. of instanteneous autoconcv)
#TODO - add tendencies for gm H, QT and QR due to rain
Expand Down Expand Up @@ -139,30 +132,28 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, r

@inbounds for k in real_center_indices(grid)
if (
en.QTvar.values[k] > epsilon &&
en.Hvar.values[k] > epsilon &&
abs(en.HQTcov.values[k]) > epsilon &&
prog_en.QTvar[k] > epsilon &&
prog_en.Hvar[k] > epsilon &&
abs(prog_en.HQTcov[k]) > epsilon &&
en.QT.values[k] > epsilon &&
sqrt(en.QTvar.values[k]) < en.QT.values[k]
sqrt(prog_en.QTvar[k]) < en.QT.values[k]
)

if en_thermo.quadrature_type == "log-normal"
# Lognormal parameters (mu, sd) from mean and variance
sd_q = sqrt(log(en.QTvar.values[k] / en.QT.values[k] / en.QT.values[k] + 1.0))
sd_h = sqrt(log(en.Hvar.values[k] / en.H.values[k] / en.H.values[k] + 1.0))
sd_q = sqrt(log(prog_en.QTvar[k] / en.QT.values[k] / en.QT.values[k] + 1.0))
sd_h = sqrt(log(prog_en.Hvar[k] / en.H.values[k] / en.H.values[k] + 1.0))
# Enforce Schwarz"s inequality
corr = max(min(en.HQTcov.values[k] / sqrt(en.Hvar.values[k] * en.QTvar.values[k]), 1.0), -1.0)
sd2_hq =
log(corr * sqrt(en.Hvar.values[k] * en.QTvar.values[k]) / en.H.values[k] / en.QT.values[k] + 1.0)
corr = max(min(prog_en.HQTcov[k] / sqrt(prog_en.Hvar[k] * prog_en.QTvar[k]), 1.0), -1.0)
sd2_hq = log(corr * sqrt(prog_en.Hvar[k] * prog_en.QTvar[k]) / en.H.values[k] / en.QT.values[k] + 1.0)
sd_cond_h_q = sqrt(max(sd_h * sd_h - sd2_hq * sd2_hq / sd_q / sd_q, 0.0))
mu_q = log(
en.QT.values[k] * en.QT.values[k] / sqrt(en.QT.values[k] * en.QT.values[k] + en.QTvar.values[k]),
)
mu_h = log(en.H.values[k] * en.H.values[k] / sqrt(en.H.values[k] * en.H.values[k] + en.Hvar.values[k]))
mu_q =
log(en.QT.values[k] * en.QT.values[k] / sqrt(en.QT.values[k] * en.QT.values[k] + prog_en.QTvar[k]))
mu_h = log(en.H.values[k] * en.H.values[k] / sqrt(en.H.values[k] * en.H.values[k] + prog_en.Hvar[k]))
else
sd_q = sqrt(en.QTvar.values[k])
sd_h = sqrt(en.Hvar.values[k])
corr = max(min(en.HQTcov.values[k] / max(sd_h * sd_q, 1e-13), 1.0), -1.0)
sd_q = sqrt(prog_en.QTvar[k])
sd_h = sqrt(prog_en.Hvar[k])
corr = max(min(prog_en.HQTcov[k] / max(sd_h * sd_q, 1e-13), 1.0), -1.0)

# limit sd_q to prevent negative qt_hat
sd_q_lim = (1e-10 - en.QT.values[k]) / (sqrt2 * abscissas[1])
Expand Down
10 changes: 5 additions & 5 deletions src/Operators.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -619,10 +619,7 @@ function construct_tridiag_diffusion_en(
param_set::APS,
state,
TS,
KM,
KH,
w_en,
tke_en,
n_updrafts::Int,
minimum_area::Float64,
pressure_plume_spacing::Vector,
Expand All @@ -645,11 +642,14 @@ function construct_tridiag_diffusion_en(
aux_tc = center_aux_tc(state)
prog_up = center_prog_updrafts(state)
prog_up_f = face_prog_updrafts(state)
prog_en = center_prog_environment(state)

ae = vec(1 .- aux_tc.bulk.area)
ae = 1 .- aux_tc.bulk.area
rho_ae_K_m = face_field(grid)
w_en_c = center_field(grid)
D_env = 0.0
KM = center_aux_tc(state).KM
KH = center_aux_tc(state).KH

aeK = is_tke ? ae .* KM : ae .* KH
aeK_bcs = (; bottom = SetValue(aeK[kc_surf]), top = SetValue(aeK[kc_toa]))
Expand Down Expand Up @@ -690,7 +690,7 @@ function construct_tridiag_diffusion_en(
rho_ae_K_m[k + 1] * Δzi * Δzi +
rho_ae_K_m[k] * Δzi * Δzi +
D_env +
ρ0_c[k] * ae[k] * c_d * sqrt(max(tke_en[k], 0)) / max(mixing_length[k], 1)
ρ0_c[k] * ae[k] * c_d * sqrt(max(prog_en.tke[k], 0)) / max(mixing_length[k], 1)
)
if is_toa_center(grid, k)
c[k] = 0.0
Expand Down
Loading