buoyancy gradient

Author: szy21

src/cache/diagnostic_edmf_precomputed_quantities.jl

@@ -1340,7 +1340,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_closures!
     t,
 )
     (; params) = p
-    (; ᶜu, ᶜts, ᶠu³) = p.precomputed
+    (; ᶜu, ᶜT, ᶜq_liq_rai, ᶜq_ice_sno, ᶠu³) = p.precomputed
     (; ustar) = p.precomputed.sfc_conditions
     (; ᶜlinear_buoygrad, ᶜstrain_rate_norm) = p.precomputed
     (; ρtke_flux) = p.precomputed
@@ -1353,10 +1353,15 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_closures!
         @. ᶜu⁰ = C123(Y.c.uₕ) + ᶜinterp(C123(ᶠu³⁰))
     end  # Set here, but used in a different function
 
+    ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
     @. ᶜlinear_buoygrad = buoyancy_gradients(
         BuoyGradMean(),
         thermo_params,
-        ᶜts,
+        ᶜT,
+        Y.c.ρ,
+        ᶜq_tot,
+        ᶜq_liq_rai,
+        ᶜq_ice_sno,
         p.precomputed.cloud_diagnostics_tuple.cf,
         C3,
         p.precomputed.ᶜgradᵥ_q_tot,

src/cache/prognostic_edmf_precomputed_quantities.jl

@@ -175,7 +175,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
     FT = eltype(params)
     n = n_mass_flux_subdomains(turbconv_model)
 
-    (; ᶜu, ᶜp, ᶠu³, ᶜts, ᶠu³⁰, ᶜts⁰) = p.precomputed
+    (; ᶜu, ᶜp, ᶠu³, ᶜT, ᶜq_liq_rai, ᶜq_ice_sno, ᶠu³⁰, ᶜts⁰) = p.precomputed
     (; ᶜlinear_buoygrad, ᶜstrain_rate_norm, ρtke_flux) = p.precomputed
     (;
         ᶜuʲs,
@@ -301,10 +301,15 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
 
     (; ᶜgradᵥ_q_tot, ᶜgradᵥ_θ_liq_ice, cloud_diagnostics_tuple) =
         p.precomputed
+    ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
     @. ᶜlinear_buoygrad = buoyancy_gradients( # TODO - do we need to modify buoyancy gradients based on NonEq + 1M tracers?
         BuoyGradMean(),
         thermo_params,
-        ᶜts,
+        ᶜT,
+        Y.c.ρ,
+        ᶜq_tot,
+        ᶜq_liq_rai,
+        ᶜq_ice_sno,
         cloud_diagnostics_tuple.cf,
         C3,
         ᶜgradᵥ_q_tot,

src/parameterized_tendencies/radiation/radiation.jl

@@ -462,21 +462,19 @@ function radiation_tendency!(Yₜ, Y, p, t, radiation_mode::RadiationDYCOMS)
     (; params) = p
     (; ᶜκρq, ∫_0_∞_κρq, ᶠ∫_0_z_κρq, isoline_z_ρ_ρq, ᶠradiation_flux) =
         p.radiation
-    (; ᶜts) = p.precomputed
+    (; ᶜts, ᶜq_liq_rai) = p.precomputed
     thermo_params = CAP.thermodynamics_params(params)
     cp_d = CAP.cp_d(params)
     FT = Spaces.undertype(axes(Y.c))
     NT = NamedTuple{(:z, :ρ, :ρq_tot), NTuple{3, FT}}
     ᶜz = Fields.coordinate_field(Y.c).z
     ᶠz = Fields.coordinate_field(Y.f).z
 
-    # TODO: According to the paper, we should replace liquid_specific_humidity
+    # TODO: According to the paper, we should replace ᶜq_liq_rai
     # with TD.mixing_ratios(thermo_params, ᶜts).liq, but this wouldn't
     # match the original code from TurbulenceConvection.
     @. ᶜκρq =
-        radiation_mode.kappa *
-        Y.c.ρ *
-        TD.liquid_specific_humidity(thermo_params, ᶜts)
+        radiation_mode.kappa * Y.c.ρ * ᶜq_liq_rai
 
     Operators.column_integral_definite!(∫_0_∞_κρq, ᶜκρq)
 
@@ -555,11 +553,11 @@ radiation_model_cache(Y, radiation_mode::RadiationISDAC; args...) = (;)  # Don't
 function radiation_tendency!(Yₜ, Y, p, t, radiation_mode::RadiationISDAC)
     (; F₀, F₁, κ) = radiation_mode
     (; params, precomputed) = p
-    (; ᶜts) = precomputed
+    (; ᶜts, ᶜq_liq_rai) = precomputed
     thermo_params = CAP.thermodynamics_params(params)
 
     ᶜρq = p.scratch.ᶜtemp_scalar
-    @. ᶜρq = Y.c.ρ * TD.liquid_specific_humidity(thermo_params, ᶜts)
+    @. ᶜρq = Y.c.ρ * ᶜq_liq_rai
 
     LWP_zₜ = p.scratch.temp_field_level  # column integral of LWP (zₜ = top-of-domain)
     Operators.column_integral_definite!(LWP_zₜ, ᶜρq)

src/prognostic_equations/eddy_diffusion_closures.jl

@@ -14,7 +14,12 @@ import SurfaceFluxes.UniversalFunctions as UF
         thermo_params,
 
         # Arguments for the first method (most commonly called):
-        ts::TD.ThermodynamicState,
+        T,      # Air temperature [K]
+        ρ,      # Air density [kg/m³]
+        q_tot,  # Total specific humidity [kg/kg]
+        q_liq,  # Liquid specific humidity [kg/kg]
+        q_ice,  # Ice specific humidity [kg/kg]
+        cf,     # Cloud fraction
         ::Type{C3}, # Covariant3 vector type, for projecting gradients
         ∂qt∂z::AbstractField,   # Vertical gradient of total specific humidity
         ∂θli∂z::AbstractField,   # Vertical gradient of liquid-ice potential temperature
@@ -36,23 +41,22 @@ fraction. The calculation involves:
    thermodynamic variables (`∂θᵥ/∂z`, `∂θₗᵢ/∂z`, `∂qₜ/∂z`), obtained from
    the input fields after projection.
 3. Blending the resulting unsaturated and saturated buoyancy gradients based on
-   the environmental cloud fraction derived from the thermodynamic state `ts`.
-
-The dispatch on `EnvBuoyGradVars` (used internally by the first method after
-constructing it from `ts` and projected gradients) occurs during its construction.
-The analytical solutions employed are consistent for both mean fields and
-conditional fields derived from assumed distributions over conserved thermodynamic
-variables.
+   the environmental cloud fraction.
 
 Arguments:
 - `closure`: The environmental buoyancy gradient closure type (e.g., `BuoyGradMean`).
 - `thermo_params`: Thermodynamic parameters from `CLIMAParameters`.
-- `ts`: Center-level thermodynamic state of the environment.
+- `T`: Air temperature [K]
+- `ρ`: Air density [kg/m³]
+- `q_tot`: Total specific humidity [kg/kg]
+- `q_liq`: Liquid specific humidity [kg/kg]
+- `q_ice`: Ice specific humidity [kg/kg]
+- `cf`: Cloud fraction
 - `C3`: The `ClimaCore.Geometry.Covariant3Vector` type, used for projecting input vertical gradients.
 - `∂qt∂z`: Field of vertical gradients of total specific humidity.
 - `∂θli∂z`: Field of vertical gradients of liquid-ice potential temperature.
 - `local_geometry`: Field of local geometry at cell centers, used for gradient projection.
-The second method takes a precomputed `EnvBuoyGradVars` object instead of `ts` and gradient fields.
+The second method takes a precomputed `EnvBuoyGradVars` object instead of T, ρ, q_tot, q_liq, q_ice and gradient fields.
 
 Returns:
 - `∂b∂z`: The mean vertical buoyancy gradient [s⁻²], as a field of scalars.
@@ -62,7 +66,11 @@ function buoyancy_gradients end
 function buoyancy_gradients(
     ebgc::AbstractEnvBuoyGradClosure,
     thermo_params,
-    ts,
+    T,
+    ρ,
+    q_tot,
+    q_liq,
+    q_ice,
     cf,
     ::Type{C3},
     ∂qt∂z,
@@ -73,7 +81,11 @@ function buoyancy_gradients(
         ebgc,
         thermo_params,
         EnvBuoyGradVars(
-            ts,
+            T,
+            ρ,
+            q_tot,
+            q_liq,
+            q_ice,
             cf,
             projected_vector_buoy_grad_vars(
                 C3,
@@ -96,18 +108,18 @@ function buoyancy_gradients(
     Rv_over_Rd = TDP.Rv_over_Rd(thermo_params)
     R_v = TDP.R_v(thermo_params)
 
-    ts = bg_model.ts
-    ∂b∂θv = g / TD.virtual_pottemp(thermo_params, ts)
+    (; T, ρ, q_tot, q_liq, q_ice) = bg_model
+    ∂b∂θv = g / TD.virtual_pottemp(thermo_params, T, ρ, q_tot, q_liq, q_ice)
 
-    T = TD.air_temperature(thermo_params, ts)
-    λ = TD.liquid_fraction(thermo_params, ts)
+    # TODO: is this correct for equilmoist?
+    λ = TD.liquid_fraction(thermo_params, T, q_liq, q_ice)
     lh =
         λ * TD.latent_heat_vapor(thermo_params, T) +
         (1 - λ) * TD.latent_heat_sublim(thermo_params, T)
-    cp_m = TD.cp_m(thermo_params, ts)
-    q_sat = TD.q_vap_saturation(thermo_params, ts)
-    q_tot = TD.total_specific_humidity(thermo_params, ts)
-    θ = TD.dry_pottemp(thermo_params, ts)
+    cp_m = TD.cp_m(thermo_params, q_tot, q_liq, q_ice)
+    # TODO: is this correct for equilmoist?
+    q_sat = TD.q_vap_saturation(thermo_params, T, ρ, q_liq, q_ice)
+    θ = TD.potential_temperature(thermo_params, T, ρ, q_tot, q_liq, q_ice)
     ∂b∂θli_unsat = ∂b∂θv * (1 + (Rv_over_Rd - 1) * q_tot)
     ∂b∂qt_unsat = ∂b∂θv * (Rv_over_Rd - 1) * θ
     ∂b∂θli_sat = (

src/prognostic_equations/gm_sgs_closures.jl

@@ -67,13 +67,18 @@ NVTX.@annotate function compute_gm_mixing_length(Y, p)
 
     ᶜdz = Fields.Δz_field(axes(Y.c))
     ᶜlg = Fields.local_geometry_field(Y.c)
-    (; ᶜts, ᶠu³, ᶜlinear_buoygrad, ᶜstrain_rate_norm, cloud_diagnostics_tuple) =
+    (; ᶜT, ᶜq_liq_rai, ᶜq_ice_sno, ᶠu³, ᶜlinear_buoygrad, ᶜstrain_rate_norm, cloud_diagnostics_tuple) =
         p.precomputed
 
+    ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
     @. ᶜlinear_buoygrad = buoyancy_gradients(
         BuoyGradMean(),
         thermo_params,
-        ᶜts,
+        ᶜT,
+        Y.c.ρ,
+        ᶜq_tot,
+        ᶜq_liq_rai,
+        ᶜq_ice_sno,
         cloud_diagnostics_tuple.cf,
         C3,
         p.precomputed.ᶜgradᵥ_q_tot,

src/solver/types.jl

@@ -318,20 +318,28 @@ Base.broadcastable(x::BuoyGradMean) = tuple(x)
 
 Variables used in the environmental buoyancy gradient computation.
 """
-Base.@kwdef struct EnvBuoyGradVars{FT, TS}
-    ts::TS
+Base.@kwdef struct EnvBuoyGradVars{FT}
+    T::FT
+    ρ::FT
+    q_tot::FT
+    q_liq::FT
+    q_ice::FT
     cf::FT
     ∂qt∂z::FT
     ∂θli∂z::FT
 end
 
 function EnvBuoyGradVars(
-    ts::TD.ThermodynamicState,
+    T,
+    ρ,
+    q_tot,
+    q_liq,
+    q_ice,
     cf,
     ∂qt∂z_∂θli∂z,
 )
     (; ∂qt∂z, ∂θli∂z) = ∂qt∂z_∂θli∂z
-    return EnvBuoyGradVars(ts, cf, ∂qt∂z, ∂θli∂z)
+    return EnvBuoyGradVars(T, ρ, q_tot, q_liq, q_ice, cf, ∂qt∂z, ∂θli∂z)
 end
 
 Base.eltype(::EnvBuoyGradVars{FT}) where {FT} = FT