From cc6b31cd1fe2312c29bba63edf587502dcabb93e Mon Sep 17 00:00:00 2001 From: ilopezgp Date: Thu, 21 Oct 2021 15:57:02 -0700 Subject: [PATCH] Channge dry/cloudy for unsat/sat. Change th/thv to \theta, \thetav. --- src/EDMF_Environment.jl | 54 ++++++++--------- src/closures/buoyancy_gradients.jl | 46 +++++++-------- src/types.jl | 94 +++++++++++++++--------------- src/update_aux.jl | 60 +++++++++---------- 4 files changed, 127 insertions(+), 127 deletions(-) diff --git a/src/EDMF_Environment.jl b/src/EDMF_Environment.jl index df2735b597..cb4db6b8e6 100644 --- a/src/EDMF_Environment.jl +++ b/src/EDMF_Environment.jl @@ -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 @@ -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 @@ -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) @@ -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 @@ -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] @@ -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 diff --git a/src/closures/buoyancy_gradients.jl b/src/closures/buoyancy_gradients.jl index f8d74ca68d..454c153e69 100644 --- a/src/closures/buoyancy_gradients.jl +++ b/src/closures/buoyancy_gradients.jl @@ -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, @@ -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 diff --git a/src/types.jl b/src/types.jl index 2b6e7f2410..f1153a12ee 100644 --- a/src/types.jl +++ b/src/types.jl @@ -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 @@ -168,22 +168,22 @@ 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" @@ -191,8 +191,8 @@ Base.@kwdef struct EnvBuoyGrad{FT, EBC <: EnvBuoyGradClosure} "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 @@ -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 @@ -355,15 +355,15 @@ 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) @@ -371,18 +371,18 @@ struct EnvironmentThermodynamics{A1} 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, diff --git a/src/update_aux.jl b/src/update_aux.jl index f8fd3fbdcb..9cacec4a7e 100644 --- a/src/update_aux.jl +++ b/src/update_aux.jl @@ -115,7 +115,7 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) rho = TD.air_density(ts_en) aux_en.buoy[k] = buoyancy_c(param_set, ρ0_c[k], rho) - update_cloud_dry(en_thermo, state, k, ts_en) + update_cloud_unsat(en_thermo, state, k, ts_en) aux_en.RH[k] = TD.relative_humidity(ts_en) ##### @@ -298,34 +298,34 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) QT_cut = ccut(aux_en.q_tot, grid, k) ∂qt∂z = c∇(QT_cut, grid, k; bottom = SetGradient(0), top = SetGradient(0)) - THL_cut = ccut(aux_en.θ_liq_ice, grid, k) - ∂θl∂z = c∇(THL_cut, grid, k; bottom = SetGradient(0), top = SetGradient(0)) + θ_liq_ice_cut = ccut(aux_en.θ_liq_ice, grid, k) + ∂θl∂z = c∇(θ_liq_ice_cut, grid, k; bottom = SetGradient(0), top = SetGradient(0)) p0_cut = ccut(p0_c, grid, k) T_cut = ccut(aux_en.T, grid, k) QT_cut = ccut(aux_en.q_tot, grid, k) QL_cut = ccut(aux_en.q_liq, grid, k) ts_cut = TD.PhaseEquil_pTq.(param_set, p0_cut, T_cut, QT_cut) - thv_cut = TD.virtual_pottemp.(ts_cut) + θv_cut = TD.virtual_pottemp.(ts_cut) ts = thermo_state_pθq(param_set, p0_c[k], aux_en.θ_liq_ice[k], aux_en.q_tot[k]) θ = TD.dry_pottemp(ts) θv = TD.virtual_pottemp(ts) - ∂θv∂z = c∇(thv_cut, grid, k; bottom = SetGradient(0), top = Extrapolate()) + ∂θv∂z = c∇(θv_cut, grid, k; bottom = SetGradient(0), top = Extrapolate()) # buoyancy_gradients if edmf.bg_closure == BuoyGradMean() # First order approximation: Use environmental mean fields. bg_kwargs = (; - t_cloudy = aux_en.T[k], + t_sat = aux_en.T[k], # WARNING: ICE MISSING - qv_cloudy = (aux_en.q_tot - aux_en.q_liq)[k], - qt_cloudy = aux_en.q_tot[k], - θ_cloudy = θ, - θ_liq_ice_cloudy = aux_en.θ_liq_ice[k], - ∂θv∂z_dry = ∂θv∂z, - ∂qt∂z_cloudy = ∂qt∂z, - ∂θl∂z_cloudy = ∂θl∂z, + qv_sat = (aux_en.q_tot - aux_en.q_liq)[k], + qt_sat = aux_en.q_tot[k], + θ_sat = θ, + θ_liq_ice_sat = aux_en.θ_liq_ice[k], + ∂θv∂z_unsat = ∂θv∂z, + ∂qt∂z_sat = ∂qt∂z, + ∂θl∂z_sat = ∂θl∂z, p0 = p0_c[k], en_cld_frac = aux_en.cloud_fraction[k], alpha0 = α0_c[k], @@ -335,9 +335,9 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) elseif edmf.bg_closure == BuoyGradQuadratures() # Second order approximation: Use dry and cloudy environmental fields. cf_cut = ccut(aux_en.cloud_fraction, grid, k) - QT_cloudy_cut = ccut(en_thermo.qt_cloudy, grid, k) - ∂qt∂z_cloudy = c∇_vanishing_subdomain( - QT_cloudy_cut, + QT_sat_cut = ccut(en_thermo.qt_sat, grid, k) + ∂qt∂z_sat = c∇_vanishing_subdomain( + QT_sat_cut, cf_cut, ∂qt∂z, grid, @@ -345,9 +345,9 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) bottom = SetGradient(0), top = SetGradient(0), ) - THL_cloudy_cut = ccut(en_thermo.thl_cloudy, grid, k) - ∂θl∂z_cloudy = c∇_vanishing_subdomain( - THL_cloudy_cut, + θ_liq_ice_sat_cut = ccut(en_thermo.θ_liq_ice_sat, grid, k) + ∂θl∂z_sat = c∇_vanishing_subdomain( + θ_liq_ice_sat_cut, cf_cut, ∂θl∂z, grid, @@ -355,9 +355,9 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) bottom = SetGradient(0), top = SetGradient(0), ) - THV_dry_cut = ccut(en_thermo.thv_dry, grid, k) - ∂θv∂z_dry = c∇_vanishing_subdomain( - THV_dry_cut, + θv_unsat_cut = ccut(en_thermo.θv_unsat, grid, k) + ∂θv∂z_unsat = c∇_vanishing_subdomain( + θv_unsat_cut, cf_cut, ∂θv∂z, grid, @@ -367,14 +367,14 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) ) bg_kwargs = (; - t_cloudy = en_thermo.t_cloudy[k], - qv_cloudy = en_thermo.qv_cloudy[k], - qt_cloudy = en_thermo.qt_cloudy[k], - θ_cloudy = en_thermo.th_cloudy[k], - θ_liq_ice_cloudy = en_thermo.thl_cloudy[k], - ∂θv∂z_dry = ∂θv∂z_dry, - ∂qt∂z_cloudy = ∂qt∂z_cloudy, - ∂θl∂z_cloudy = ∂θl∂z_cloudy, + t_sat = en_thermo.t_sat[k], + qv_sat = en_thermo.qv_sat[k], + qt_sat = en_thermo.qt_sat[k], + θ_sat = en_thermo.th_sat[k], + θ_liq_ice_sat = en_thermo.θ_liq_ice_sat[k], + ∂θv∂z_unsat = ∂θv∂z_unsat, + ∂qt∂z_sat = ∂qt∂z_sat, + ∂θl∂z_sat = ∂θl∂z_sat, p0 = p0_c[k], en_cld_frac = aux_en.cloud_fraction[k], alpha0 = α0_c[k],