use cloud fraction in buoyancy gradient

Author: szy21

src/cache/diagnostic_edmf_precomputed_quantities.jl

@@ -1367,9 +1367,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_closures!
         moisture_model,
         ᶜts,
         C3,
-        p.precomputed.ᶜgradᵥ_θ_virt,    # ∂θv∂z_unsat
-        p.precomputed.ᶜgradᵥ_q_tot,     # ∂qt∂z_sat
-        p.precomputed.ᶜgradᵥ_θ_liq_ice, # ∂θl∂z_sat
+        p.precomputed.ᶜgradᵥ_q_tot,
+        p.precomputed.ᶜgradᵥ_θ_liq_ice,
         ᶜlg,
     )
 

src/cache/precomputed_quantities.jl

@@ -85,7 +85,6 @@ function precomputed_quantities(Y, atmos)
     TST = thermo_state_type(atmos.moisture_model, FT)
     SCT = SurfaceConditions.surface_conditions_type(atmos, FT)
     cspace = axes(Y.c)
-    fspace = axes(Y.f)
     n = n_mass_flux_subdomains(atmos.turbconv_model)
     n_prog = n_prognostic_mass_flux_subdomains(atmos.turbconv_model)
     @assert !(atmos.turbconv_model isa PrognosticEDMFX) || n_prog == 1
@@ -98,6 +97,7 @@ function precomputed_quantities(Y, atmos)
     )
     cloud_diagnostics_tuple =
         similar(Y.c, @NamedTuple{cf::FT, q_liq::FT, q_ice::FT})
+    @. cloud_diagnostics_tuple.cf = FT(0)
     surface_precip_fluxes = (;
         surface_rain_flux = zeros(axes(Fields.level(Y.f, half))),
         surface_snow_flux = zeros(axes(Fields.level(Y.f, half))),

src/cache/prognostic_edmf_precomputed_quantities.jl

@@ -421,20 +421,20 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
         )
     end
 
-    (; ᶜgradᵥ_θ_virt, ᶜgradᵥ_q_tot, ᶜgradᵥ_θ_liq_ice) = p.precomputed
+    (; ᶜgradᵥ_θ_virt, ᶜgradᵥ_q_tot, ᶜgradᵥ_θ_liq_ice, cloud_diagnostics_tuple) = p.precomputed
     # First order approximation: Use environmental mean fields.
-    @. ᶜgradᵥ_θ_virt = ᶜgradᵥ(ᶠinterp(TD.virtual_pottemp(thermo_params, ᶜts)))       # ∂θv∂z_unsat
+    @. ᶜgradᵥ_θ_virt = ᶜgradᵥ(ᶠinterp(TD.virtual_pottemp(thermo_params, ᶜts)))
     ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
-    @. ᶜgradᵥ_q_tot = ᶜgradᵥ(ᶠinterp(ᶜq_tot))                                        # ∂qt∂z_sat
+    @. ᶜgradᵥ_q_tot = ᶜgradᵥ(ᶠinterp(ᶜq_tot))
     @. ᶜgradᵥ_θ_liq_ice =
-        ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts)))                    # ∂θl∂z_sat
+        ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts)))
     @. ᶜlinear_buoygrad = buoyancy_gradients( # TODO - do we need to modify buoyancy gradients based on NonEq + 1M tracers?
         BuoyGradMean(),
         thermo_params,
         moisture_model,
         ᶜts,
+        cloud_diagnostics_tuple.cf,
         C3,
-        ᶜgradᵥ_θ_virt,
         ᶜgradᵥ_q_tot,
         ᶜgradᵥ_θ_liq_ice,
         ᶜlg,

src/prognostic_equations/eddy_diffusion_closures.jl

@@ -15,7 +15,6 @@ import ClimaCore.Fields as Fields
         # Arguments for the first method (most commonly called):
         ts::TD.ThermodynamicState,
         ::Type{C3}, # Covariant3 vector type, for projecting gradients
-        ∂θv∂z_unsat::AbstractField, # Vertical gradient of virtual potential temperature in unsaturated part
         ∂qt∂z_sat::AbstractField,   # Vertical gradient of total specific humidity in saturated part
         ∂θli∂z_sat::AbstractField,   # Vertical gradient of liquid-ice potential temperature in saturated part
         local_geometry::Fields.LocalGeometry,
@@ -50,7 +49,6 @@ Arguments:
 - `moisture_model`: Moisture model (e.g., `EquilMoistModel`, `NonEquilMoistModel`).
 - `ts`: Center-level thermodynamic state of the environment.
 - `C3`: The `ClimaCore.Geometry.Covariant3Vector` type, used for projecting input vertical gradients.
-- `∂θv∂z_unsat`: Field of vertical gradients of virtual potential temperature in the unsaturated part.
 - `∂qt∂z_sat`: Field of vertical gradients of total specific humidity in the saturated part.
 - `∂θli∂z_sat`: Field of vertical gradients of liquid-ice potential temperature in the saturated part.
 - `local_geometry`: Field of local geometry at cell centers, used for gradient projection.
@@ -66,8 +64,8 @@ function buoyancy_gradients(
     thermo_params,
     moisture_model,
     ts,
+    cf,
     ::Type{C3},
-    ∂θv∂z_unsat,
     ∂qt∂z_sat,
     ∂θli∂z_sat,
     ᶜlg,
@@ -78,9 +76,9 @@ function buoyancy_gradients(
         moisture_model,
         EnvBuoyGradVars(
             ts,
+            cf,
             projected_vector_buoy_grad_vars(
                 C3,
-                ∂θv∂z_unsat,
                 ∂qt∂z_sat,
                 ∂θli∂z_sat,
                 ᶜlg,
@@ -99,13 +97,10 @@ function buoyancy_gradients(
 
     g = TDP.grav(thermo_params)
     Rv_over_Rd = TDP.Rv_over_Rd(thermo_params)
-    R_d = TDP.R_d(thermo_params)
     R_v = TDP.R_v(thermo_params)
 
     ts = bg_model.ts
-    p = TD.air_pressure(thermo_params, ts)
-    Π = TD.exner_given_pressure(thermo_params, p)
-    ∂b∂θv = g * (R_d * TD.air_density(thermo_params, ts) / p) * Π
+    ∂b∂θv = g / TD.virtual_pottemp(thermo_params, ts)
 
     t_sat = TD.air_temperature(thermo_params, ts)
     λ = TD.liquid_fraction(thermo_params, ts)
@@ -115,20 +110,25 @@ function buoyancy_gradients(
     cp_m = TD.cp_m(thermo_params, ts)
     qv_sat = TD.vapor_specific_humidity(thermo_params, ts)
     qt_sat = TD.total_specific_humidity(thermo_params, ts)
+    θ = TD.dry_pottemp(thermo_params, ts)
+    ∂b∂θli_unsat = ∂b∂θv * (1 + (Rv_over_Rd - 1) * qt_sat)
+    ∂b∂qt_unsat = ∂b∂θv * (Rv_over_Rd - 1) * θ
     ∂b∂θli_sat = (
         ∂b∂θv *
         (1 + Rv_over_Rd * (1 + lh / R_v / t_sat) * qv_sat - qt_sat) /
-        (1 + lh * lh / cp_m / R_v / t_sat / t_sat * qv_sat)
+        (1 + lh^2 / cp_m / R_v / t_sat^2 * qv_sat)
     )
     ∂b∂qt_sat =
         (lh / cp_m / t_sat * ∂b∂θli_sat - ∂b∂θv) *
-        TD.dry_pottemp(thermo_params, ts)
+        θ
 
     ∂b∂z = buoyancy_gradient_chain_rule(
         ebgc,
         bg_model,
         thermo_params,
         ∂b∂θv,
+        ∂b∂θli_unsat,
+        ∂b∂qt_unsat,
         ∂b∂θli_sat,
         ∂b∂qt_sat,
     )
@@ -177,17 +177,21 @@ function buoyancy_gradient_chain_rule(
     bg_model::EnvBuoyGradVars,
     thermo_params,
     ∂b∂θv,
+    ∂b∂θli_unsat,
+    ∂b∂qt_unsat,
     ∂b∂θli_sat,
     ∂b∂qt_sat,
 )
-    en_cld_frac = ifelse(TD.has_condensate(thermo_params, bg_model.ts), 1, 0)
+    #en_cld_frac = ifelse(TD.has_condensate(thermo_params, bg_model.ts), 1, 0)
 
+    ∂b∂z_θli_unsat = ∂b∂θli_unsat * bg_model.∂θli∂z_sat
+    ∂b∂z_qt_unsat = ∂b∂qt_unsat * bg_model.∂qt∂z_sat
+    ∂b∂z_unsat = ∂b∂z_θli_unsat + ∂b∂z_qt_unsat
     ∂b∂z_θl_sat = ∂b∂θli_sat * bg_model.∂θli∂z_sat
     ∂b∂z_qt_sat = ∂b∂qt_sat * bg_model.∂qt∂z_sat
     ∂b∂z_sat = ∂b∂z_θl_sat + ∂b∂z_qt_sat
-    ∂b∂z_unsat = ∂b∂θv * bg_model.∂θv∂z_unsat
 
-    ∂b∂z = (1 - en_cld_frac) * ∂b∂z_unsat + en_cld_frac * ∂b∂z_sat
+    ∂b∂z = (1 - bg_model.cf) * ∂b∂z_unsat + bg_model.cf * ∂b∂z_sat
 
     return ∂b∂z
 end

src/prognostic_equations/gm_sgs_closures.jl

@@ -75,9 +75,8 @@ NVTX.@annotate function compute_gm_mixing_length(Y, p)
         p.atmos.moisture_model,
         ᶜts,
         C3,
-        p.precomputed.ᶜgradᵥ_θ_virt,    # ∂θv∂z_unsat
-        p.precomputed.ᶜgradᵥ_q_tot,     # ∂qt∂z_sat
-        p.precomputed.ᶜgradᵥ_θ_liq_ice, # ∂θl∂z_sat
+        p.precomputed.ᶜgradᵥ_q_tot,
+        p.precomputed.ᶜgradᵥ_θ_liq_ice,
         ᶜlg,
     )
 

src/solver/types.jl

@@ -376,17 +376,18 @@ Variables used in the environmental buoyancy gradient computation.
 """
 Base.@kwdef struct EnvBuoyGradVars{FT, TS}
     ts::TS
-    ∂θv∂z_unsat::FT
+    cf::FT
     ∂qt∂z_sat::FT
     ∂θli∂z_sat::FT
 end
 
 function EnvBuoyGradVars(
     ts::TD.ThermodynamicState,
-    ∂θv∂z_unsat_∂qt∂z_sat_∂θli∂z_sat,
+    cf,
+    ∂qt∂z_sat_∂θli∂z_sat,
 )
-    (; ∂θv∂z_unsat, ∂qt∂z_sat, ∂θli∂z_sat) = ∂θv∂z_unsat_∂qt∂z_sat_∂θli∂z_sat
-    return EnvBuoyGradVars(ts, ∂θv∂z_unsat, ∂qt∂z_sat, ∂θli∂z_sat)
+    (; ∂qt∂z_sat, ∂θli∂z_sat) = ∂qt∂z_sat_∂θli∂z_sat
+    return EnvBuoyGradVars(ts, cf, ∂qt∂z_sat, ∂θli∂z_sat)
 end
 
 Base.eltype(::EnvBuoyGradVars{FT}) where {FT} = FT

src/utils/utilities.jl

@@ -279,12 +279,11 @@ The type should correspond to a vector with only one component, i.e., a basis ve
 projected_vector_data(::Type{V}, vector, local_geometry) where {V} =
     V(vector, local_geometry)[1] / unit_basis_vector_data(V, local_geometry)
 
-function projected_vector_buoy_grad_vars(::Type{V}, v1, v2, v3, lg) where {V}
+function projected_vector_buoy_grad_vars(::Type{V}, v1, v2, lg) where {V}
     ubvd = unit_basis_vector_data(V, lg)
     return (;
-        ∂θv∂z_unsat = V(v1, lg)[1] / ubvd,
-        ∂qt∂z_sat = V(v2, lg)[1] / ubvd,
-        ∂θli∂z_sat = V(v3, lg)[1] / ubvd,
+        ∂qt∂z_sat = V(v1, lg)[1] / ubvd,
+        ∂θli∂z_sat = V(v2, lg)[1] / ubvd,
     )
 end