Skip to content

Commit

Permalink
Channge dry/cloudy for unsat/sat. Change th/thv to \theta, \thetav.
Browse files Browse the repository at this point in the history
  • Loading branch information
@ilopezgp
ilopezgp committed Oct 21, 2021
1 parent 863d9e3 commit cc6b31c
Show file tree
Hide file tree
Showing 4 changed files with 127 additions and 127 deletions.
54 changes: 27 additions & 27 deletions src/EDMF_Environment.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,21 +47,21 @@ function update_env_precip_tendencies(en_thermo::EnvironmentThermodynamics, stat
return
end

function update_cloud_dry(en_thermo::EnvironmentThermodynamics, state, k, ts)
function update_sat_unsat(en_thermo::EnvironmentThermodynamics, state, k, ts)
q_liq = TD.liquid_specific_humidity(ts)
aux_en = center_aux_environment(state)
if q_liq > 0.0
aux_en.cloud_fraction[k] = 1.0
en_thermo.th_cloudy[k] = TD.dry_pottemp(ts)
en_thermo.thl_cloudy[k] = TD.liquid_ice_pottemp(ts)
en_thermo.t_cloudy[k] = TD.air_temperature(ts)
en_thermo.qt_cloudy[k] = TD.total_specific_humidity(ts)
en_thermo.qv_cloudy[k] = TD.vapor_specific_humidity(ts)
en_thermo.th_sat[k] = TD.dry_pottemp(ts)
en_thermo.θ_liq_ice_sat[k] = TD.liquid_ice_pottemp(ts)
en_thermo.t_sat[k] = TD.air_temperature(ts)
en_thermo.qt_sat[k] = TD.total_specific_humidity(ts)
en_thermo.qv_sat[k] = TD.vapor_specific_humidity(ts)
else
aux_en.cloud_fraction[k] = 0.0
en_thermo.th_dry[k] = TD.dry_pottemp(ts)
en_thermo.thv_dry[k] = TD.virtual_pottemp(ts)
en_thermo.qt_dry[k] = TD.total_specific_humidity(ts)
en_thermo.θ_unsat[k] = TD.dry_pottemp(ts)
en_thermo.θv_unsat[k] = TD.virtual_pottemp(ts)
en_thermo.qt_unsat[k] = TD.total_specific_humidity(ts)
end
return
end
Expand All @@ -79,7 +79,7 @@ function sgs_mean(en_thermo::EnvironmentThermodynamics, grid, state, en, rain, d
ts = thermo_state_pθq(param_set, p0_c[k], aux_en.θ_liq_ice[k], q_tot_en)
# autoconversion and accretion
mph = precipitation_formation(param_set, rain.rain_model, prog_ra.qr[k], aux_en.area[k], ρ0_c[k], dt, ts)
update_cloud_dry(en_thermo, state, k, ts)
update_sat_unsat(en_thermo, state, k, ts)
update_env_precip_tendencies(en_thermo, state, k, mph.qt_tendency, mph.θ_liq_ice_tendency)
end
return
Expand Down Expand Up @@ -117,7 +117,7 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, r
outer_env = zeros(env_len)
inner_src = zeros(src_len)
outer_src = zeros(src_len)
i_ql, i_T, i_cf, i_qt_cld, i_qt_dry, i_T_cld, i_T_dry, i_rf = 1:env_len
i_ql, i_T, i_cf, i_qt_sat, i_qt_unsat, i_T_sat, i_T_unsat, i_rf = 1:env_len
i_SH_qt, i_Sqt_H, i_SH_H, i_Sqt_qt, i_Sqt, i_SH = 1:src_len

@inbounds for k in real_center_indices(grid)
Expand Down Expand Up @@ -216,11 +216,11 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, r
# cloudy/dry categories for buoyancy in TKE
if q_liq_en > 0.0
inner_env[i_cf] += weights[m_h] * sqpi_inv
inner_env[i_qt_cld] += qt_hat * weights[m_h] * sqpi_inv
inner_env[i_T_cld] += T * weights[m_h] * sqpi_inv
inner_env[i_qt_sat] += qt_hat * weights[m_h] * sqpi_inv
inner_env[i_T_sat] += T * weights[m_h] * sqpi_inv
else
inner_env[i_qt_dry] += qt_hat * weights[m_h] * sqpi_inv
inner_env[i_T_dry] += T * weights[m_h] * sqpi_inv
inner_env[i_qt_unsat] += qt_hat * weights[m_h] * sqpi_inv
inner_env[i_T_unsat] += T * weights[m_h] * sqpi_inv
end
# products for variance and covariance source terms
inner_src[i_Sqt] += mph.qt_tendency * weights[m_h] * sqpi_inv
Expand All @@ -244,21 +244,21 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, r

# update cloudy/dry variables for buoyancy in TKE
aux_en.cloud_fraction[k] = outer_env[i_cf]
en_thermo.qt_dry[k] = outer_env[i_qt_dry]
# TD.jl cannot compute θ_dry when T=0
en_thermo.th_dry[k] = if outer_env[i_T_dry] > 0
ts_dry = TD.PhaseEquil_pTq(param_set, p0_c[k], outer_env[i_T_dry], en_thermo.qt_dry[k])
TD.dry_pottemp(ts_dry)
en_thermo.qt_unsat[k] = outer_env[i_qt_unsat]
# TD.jl cannot compute θ_unsat when T=0
en_thermo.θ_unsat[k] = if outer_env[i_T_unsat] > 0
ts_unsat = TD.PhaseEquil_pTq(param_set, p0_c[k], outer_env[i_T_unsat], en_thermo.qt_unsat[k])
TD.dry_pottemp(ts_unsat)
else
0
end

en_thermo.t_cloudy[k] = outer_env[i_T_cld]
en_thermo.qv_cloudy[k] = outer_env[i_qt_cld] - outer_env[i_ql]
en_thermo.qt_cloudy[k] = outer_env[i_qt_cld]
ts_cld = TD.PhaseEquil_pTq(param_set, p0_c[k], en_thermo.t_cloudy[k], en_thermo.qt_cloudy[k])
en_thermo.th_cloudy[k] = TD.dry_pottemp(ts_cld)
en_thermo.thl_cloudy[k] = TD.liquid_ice_pottemp(ts_cld)
en_thermo.t_sat[k] = outer_env[i_T_sat]
en_thermo.qv_sat[k] = outer_env[i_qt_sat] - outer_env[i_ql]
en_thermo.qt_sat[k] = outer_env[i_qt_sat]
ts_sat = TD.PhaseEquil_pTq(param_set, p0_c[k], en_thermo.t_sat[k], en_thermo.qt_sat[k])
en_thermo.th_sat[k] = TD.dry_pottemp(ts_sat)
en_thermo.θ_liq_ice_sat[k] = TD.liquid_ice_pottemp(ts_sat)

# update var/covar rain sources
en_thermo.Hvar_rain_dt[k] = outer_src[i_SH_H] - outer_src[i_SH] * aux_en.θ_liq_ice[k]
Expand All @@ -272,7 +272,7 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, r
ts = thermo_state_pθq(param_set, p0_c[k], aux_en.θ_liq_ice[k], aux_en.q_tot[k])
mph = precipitation_formation(param_set, rain.rain_model, prog_ra.qr[k], aux_en.area[k], ρ0_c[k], dt, ts)
update_env_precip_tendencies(en_thermo, state, k, mph.qt_tendency, mph.θ_liq_ice_tendency)
update_cloud_dry(en_thermo, state, k, ts)
update_sat_unsat(en_thermo, state, k, ts)

en_thermo.Hvar_rain_dt[k] = 0.0
en_thermo.QTvar_rain_dt[k] = 0.0
Expand Down
46 changes: 23 additions & 23 deletions src/closures/buoyancy_gradients.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,30 +19,30 @@ function buoyancy_gradients(param_set, bg_model::EnvBuoyGrad{FT, EBG}) where {FT
∂b∂θv = g * (R_d / bg_model.alpha0 / bg_model.p0) * Π

if bg_model.en_cld_frac > 0.0
ts_cloudy = thermo_state_pθq(param_set, bg_model.p0, bg_model.θ_liq_ice_cloudy, bg_model.qt_cloudy)
phase_part = TD.PhasePartition(ts_cloudy)
ts_sat = thermo_state_pθq(param_set, bg_model.p0, bg_model.θ_liq_ice_sat, bg_model.qt_sat)
phase_part = TD.PhasePartition(ts_sat)
lh = TD.latent_heat_liq_ice(param_set, phase_part)
cp_m = TD.cp_m(ts_cloudy)
∂b∂θl_cld = (
∂b∂θv * (1 + molmass_ratio * (1 + lh / R_v / bg_model.t_cloudy) * bg_model.qv_cloudy - bg_model.qt_cloudy) /
(1 + lh * lh / cp_m / R_v / bg_model.t_cloudy / bg_model.t_cloudy * bg_model.qv_cloudy)
cp_m = TD.cp_m(ts_sat)
∂b∂θl_sat = (
∂b∂θv * (1 + molmass_ratio * (1 + lh / R_v / bg_model.t_sat) * bg_model.qv_sat - bg_model.qt_sat) /
(1 + lh * lh / cp_m / R_v / bg_model.t_sat / bg_model.t_sat * bg_model.qv_sat)
)
∂b∂qt_cld = (lh / cp_m / bg_model.t_cloudy * ∂b∂θl_cld - ∂b∂θv) * bg_model.θ_cloudy
∂b∂qt_sat = (lh / cp_m / bg_model.t_sat * ∂b∂θl_sat - ∂b∂θv) * bg_model.θ_sat
else
∂b∂θl_cld = FT(0)
∂b∂qt_cld = FT(0)
∂b∂θl_sat = FT(0)
∂b∂qt_sat = FT(0)
end

∂b∂z, ∂b∂z_dry, ∂b∂z_cloudy = buoyancy_gradient_chain_rule(bg_model, ∂b∂θv, ∂b∂θl_cld, ∂b∂qt_cld)
return GradBuoy(∂b∂z, ∂b∂z_dry, ∂b∂z_cloudy)
∂b∂z, ∂b∂z_unsat, ∂b∂z_sat = buoyancy_gradient_chain_rule(bg_model, ∂b∂θv, ∂b∂θl_sat, ∂b∂qt_sat)
return GradBuoy(∂b∂z, ∂b∂z_unsat, ∂b∂z_sat)
end

"""
buoyancy_gradient_chain_rule(
bg_model::EnvBuoyGrad{FT, EBG},
∂b∂θv::FT,
∂b∂θl_cld::FT,
∂b∂qt_cld::FT,
∂b∂θl_sat::FT,
∂b∂qt_sat::FT,
) where {FT <: Real, EBG <: EnvBuoyGradClosure}
Returns the vertical buoyancy gradients in the environment, as well as in its dry and cloudy volume fractions,
Expand All @@ -51,21 +51,21 @@ from the partial derivatives with respect to thermodynamic variables in dry and
function buoyancy_gradient_chain_rule(
bg_model::EnvBuoyGrad{FT, EBG},
∂b∂θv::FT,
∂b∂θl_cld::FT,
∂b∂qt_cld::FT,
∂b∂θl_sat::FT,
∂b∂qt_sat::FT,
) where {FT <: Real, EBG <: EnvBuoyGradClosure}
if bg_model.en_cld_frac > FT(0)
∂b∂z_θl_cld = ∂b∂θl_cld * bg_model.∂θl∂z_cloudy
∂b∂z_qt_cld = ∂b∂qt_cld * bg_model.∂qt∂z_cloudy
∂b∂z_θl_sat = ∂b∂θl_sat * bg_model.∂θl∂z_sat
∂b∂z_qt_sat = ∂b∂qt_sat * bg_model.∂qt∂z_sat
else
∂b∂z_θl_cld = FT(0)
∂b∂z_qt_cld = FT(0)
∂b∂z_θl_sat = FT(0)
∂b∂z_qt_sat = FT(0)
end

∂b∂z_dry = bg_model.en_cld_frac < FT(1) ? ∂b∂θv * bg_model.∂θv∂z_dry : FT(0)
∂b∂z_unsat = bg_model.en_cld_frac < FT(1) ? ∂b∂θv * bg_model.∂θv∂z_unsat : FT(0)

∂b∂z_cloudy = ∂b∂z_θl_cld + ∂b∂z_qt_cld
∂b∂z = (1 - bg_model.en_cld_frac) * ∂b∂z_dry + bg_model.en_cld_frac * ∂b∂z_cloudy
∂b∂z_sat = ∂b∂z_θl_sat + ∂b∂z_qt_sat
∂b∂z = (1 - bg_model.en_cld_frac) * ∂b∂z_unsat + bg_model.en_cld_frac * ∂b∂z_sat

return ∂b∂z, ∂b∂z_dry, ∂b∂z_cloudy
return ∂b∂z, ∂b∂z_unsat, ∂b∂z_sat
end
94 changes: 47 additions & 47 deletions src/types.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,10 +150,10 @@ $(DocStringExtensions.FIELDS)
Base.@kwdef struct GradBuoy{FT}
"environmental vertical buoyancy gradient"
∂b∂z::FT
"vertical buoyancy gradient in the dry part of the environment"
∂b∂z_dry::FT
"vertical buoyancy gradient in the cloudy part of the environment"
∂b∂z_cloudy::FT
"vertical buoyancy gradient in the unsaturated part of the environment"
∂b∂z_unsat::FT
"vertical buoyancy gradient in the saturated part of the environment"
∂b∂z_sat::FT
end

abstract type EnvBuoyGradClosure end
Expand All @@ -168,31 +168,31 @@ Variables used in the environmental buoyancy gradient computation.
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct EnvBuoyGrad{FT, EBC <: EnvBuoyGradClosure}
"temperature in the cloudy part"
t_cloudy::FT
"vapor specific humidity in the cloudy part"
qv_cloudy::FT
"total specific humidity in the cloudy part"
qt_cloudy::FT
"potential temperature in the cloudy part"
θ_cloudy::FT
"liquid ice potential temperature in the cloudy part"
θ_liq_ice_cloudy::FT
"virtual potential temperature gradient in the non cloudy part"
∂θv∂z_dry::FT
"total specific humidity gradient in the cloudy part"
∂qt∂z_cloudy::FT
"liquid ice potential temperature gradient in the cloudy part"
∂θl∂z_cloudy::FT
"temperature in the saturated part"
t_sat::FT
"vapor specific humidity in the saturated part"
qv_sat::FT
"total specific humidity in the saturated part"
qt_sat::FT
"potential temperature in the saturated part"
θ_sat::FT
"liquid ice potential temperature in the saturated part"
θ_liq_ice_sat::FT
"virtual potential temperature gradient in the non saturated part"
∂θv∂z_unsat::FT
"total specific humidity gradient in the saturated part"
∂qt∂z_sat::FT
"liquid ice potential temperature gradient in the saturated part"
∂θl∂z_sat::FT
"reference pressure"
p0::FT
"cloud fraction"
en_cld_frac::FT
"specific volume"
alpha0::FT
end
function EnvBuoyGrad(::EBG; t_cloudy::FT, bg_kwargs...) where {FT <: Real, EBG <: EnvBuoyGradClosure}
return EnvBuoyGrad{FT, EBG}(; t_cloudy, bg_kwargs...)
function EnvBuoyGrad(::EBG; t_sat::FT, bg_kwargs...) where {FT <: Real, EBG <: EnvBuoyGradClosure}
return EnvBuoyGrad{FT, EBG}(; t_sat, bg_kwargs...)
end

Base.@kwdef mutable struct RainVariables
Expand Down Expand Up @@ -338,14 +338,14 @@ end
struct EnvironmentThermodynamics{A1}
quadrature_order::Int
quadrature_type::String
qt_dry::A1
th_dry::A1
thv_dry::A1
t_cloudy::A1
qv_cloudy::A1
qt_cloudy::A1
th_cloudy::A1
thl_cloudy::A1
qt_unsat::A1
θ_unsat::A1
θv_unsat::A1
t_sat::A1
qv_sat::A1
qt_sat::A1
th_sat::A1
θ_liq_ice_sat::A1
Hvar_rain_dt::A1
QTvar_rain_dt::A1
HQTcov_rain_dt::A1
Expand All @@ -355,34 +355,34 @@ struct EnvironmentThermodynamics{A1}
quadrature_order = parse_namelist(namelist, "thermodynamics", "quadrature_order"; default = 3)
quadrature_type = parse_namelist(namelist, "thermodynamics", "quadrature_type"; default = "gaussian")

qt_dry = center_field(grid)
th_dry = center_field(grid)
thv_dry = center_field(grid)
qt_unsat = center_field(grid)
θ_unsat = center_field(grid)
θv_unsat = center_field(grid)

t_cloudy = center_field(grid)
qv_cloudy = center_field(grid)
qt_cloudy = center_field(grid)
th_cloudy = center_field(grid)
thl_cloudy = center_field(grid)
t_sat = center_field(grid)
qv_sat = center_field(grid)
qt_sat = center_field(grid)
th_sat = center_field(grid)
θ_liq_ice_sat = center_field(grid)

Hvar_rain_dt = center_field(grid)
QTvar_rain_dt = center_field(grid)
HQTcov_rain_dt = center_field(grid)

qt_tendency_rain_formation = center_field(grid)
θ_liq_ice_tendency_rain_formation = center_field(grid)
A1 = typeof(qt_dry)
A1 = typeof(qt_unsat)
return new{A1}(
quadrature_order,
quadrature_type,
qt_dry,
th_dry,
thv_dry,
t_cloudy,
qv_cloudy,
qt_cloudy,
th_cloudy,
thl_cloudy,
qt_unsat,
θ_unsat,
θv_unsat,
t_sat,
qv_sat,
qt_sat,
th_sat,
θ_liq_ice_sat,
Hvar_rain_dt,
QTvar_rain_dt,
HQTcov_rain_dt,
Expand Down
Loading

0 comments on commit cc6b31c

Please sign in to comment.