Skip to content

Commit

Permalink
Merge #418
Browse files Browse the repository at this point in the history 
418: Make several variables ClimaCore fields r=charleskawczynski a=charleskawczynski

This PR:
 - Makes several EDMF and 2nd moment environment variables fields
 - Index bugfix into `rho_ae_KH`/`rho_ae_KM`, which are cell face fields. Since ClimaCore fields explicitly require the correct index type, I suspect that we may find more of these as we convert more fields. It's still non-behavior changing since `kc_surface(grid) == kf_surface(grid)`.
 - Completely removes the `EnvironmentVariable_2m` struct

Co-authored-by: Charles Kawczynski <kawczynski.charles@gmail.com>
  • Loading branch information
@bors @charleskawczynski
bors[bot] and charleskawczynski committed Oct 17, 2021
2 parents 646c4ca + 82f577b commit 0e339c6
Show file tree
Hide file tree
Showing 9 changed files with 242 additions and 302 deletions.
23 changes: 21 additions & 2 deletions integration_tests/utils/main.jl
View file
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
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
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

0 comments on commit 0e339c6

Please sign in to comment.