From a832aeaed9633f25d2f4c83409eef772ab08e50b Mon Sep 17 00:00:00 2001 From: Zhaoyi Shen <11598433+szy21@users.noreply.github.com> Date: Tue, 25 Nov 2025 23:21:58 -0800 Subject: [PATCH] use cloud fraction in buoyancy gradient --- .../diagnostic_edmf_precomputed_quantities.jl | 5 ++-- src/cache/precomputed_quantities.jl | 2 +- .../prognostic_edmf_precomputed_quantities.jl | 10 +++---- .../eddy_diffusion_closures.jl | 30 +++++++++++-------- src/prognostic_equations/gm_sgs_closures.jl | 5 ++-- src/solver/types.jl | 9 +++--- src/utils/utilities.jl | 7 ++--- 7 files changed, 35 insertions(+), 33 deletions(-) diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl index 8aac50a7ab..54c0b1c76b 100644 --- a/src/cache/diagnostic_edmf_precomputed_quantities.jl +++ b/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, ) diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index 0f9d00fb85..0f03cd9a84 100644 --- a/src/cache/precomputed_quantities.jl +++ b/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))), diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl index 9322b93a86..c2c2dec8bb 100644 --- a/src/cache/prognostic_edmf_precomputed_quantities.jl +++ b/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, diff --git a/src/prognostic_equations/eddy_diffusion_closures.jl b/src/prognostic_equations/eddy_diffusion_closures.jl index af65232463..d1b8bcddd7 100644 --- a/src/prognostic_equations/eddy_diffusion_closures.jl +++ b/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 diff --git a/src/prognostic_equations/gm_sgs_closures.jl b/src/prognostic_equations/gm_sgs_closures.jl index 274a27169b..be65438f84 100644 --- a/src/prognostic_equations/gm_sgs_closures.jl +++ b/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, ) diff --git a/src/solver/types.jl b/src/solver/types.jl index decb1efa7c..368538f382 100644 --- a/src/solver/types.jl +++ b/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 diff --git a/src/utils/utilities.jl b/src/utils/utilities.jl index fe0b74a10a..be8b42e106 100644 --- a/src/utils/utilities.jl +++ b/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