use cloud fraction in buoyancy gradient

Author: szy21

src/cache/cloud_fraction.jl

@@ -22,7 +22,6 @@ end
    Compute the grid scale cloud fraction based on sub-grid scale properties
 """
 NVTX.@annotate function set_cloud_fraction!(Y, p, ::DryModel, _)
-    (; turbconv_model) = p.atmos
     FT = eltype(p.params)
 
     p.precomputed.cloud_diagnostics_tuple .=
@@ -35,23 +34,9 @@ NVTX.@annotate function set_cloud_fraction!(
     ::GridScaleCloud,
 )
     (; params) = p
-    (; turbconv_model) = p.atmos
     (; ᶜts, cloud_diagnostics_tuple) = p.precomputed
     thermo_params = CAP.thermodynamics_params(params)
 
-    if isnothing(turbconv_model)
-        if p.atmos.call_cloud_diagnostics_per_stage isa
-           CallCloudDiagnosticsPerStage
-            (; ᶜgradᵥ_θ_virt, ᶜgradᵥ_q_tot, ᶜgradᵥ_θ_liq_ice) = p.precomputed
-            thermo_params = CAP.thermodynamics_params(p.params)
-            @. ᶜgradᵥ_θ_virt =
-                ᶜgradᵥ(ᶠinterp(TD.virtual_pottemp(thermo_params, ᶜts)))
-            @. ᶜgradᵥ_q_tot =
-                ᶜgradᵥ(ᶠinterp(TD.total_specific_humidity(thermo_params, ᶜts)))
-            @. ᶜgradᵥ_θ_liq_ice =
-                ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts)))
-        end
-    end
     if moist_model isa EquilMoistModel
         @. cloud_diagnostics_tuple = make_named_tuple(
             ifelse(TD.has_condensate(thermo_params, ᶜts), 1, 0),
@@ -83,10 +68,8 @@ NVTX.@annotate function set_cloud_fraction!(
     if isnothing(turbconv_model)
         if p.atmos.call_cloud_diagnostics_per_stage isa
            CallCloudDiagnosticsPerStage
-            (; ᶜgradᵥ_θ_virt, ᶜgradᵥ_q_tot, ᶜgradᵥ_θ_liq_ice) = p.precomputed
+            (; ᶜgradᵥ_q_tot, ᶜgradᵥ_θ_liq_ice) = p.precomputed
             thermo_params = CAP.thermodynamics_params(p.params)
-            @. ᶜgradᵥ_θ_virt =
-                ᶜgradᵥ(ᶠinterp(TD.virtual_pottemp(thermo_params, ᶜts)))
             @. ᶜgradᵥ_q_tot =
                 ᶜgradᵥ(ᶠinterp(TD.total_specific_humidity(thermo_params, ᶜts)))
             @. ᶜgradᵥ_θ_liq_ice =

src/cache/diagnostic_edmf_precomputed_quantities.jl

@@ -1366,10 +1366,10 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_closures!
         thermo_params,
         moisture_model,
         ᶜts,
+        p.precomputed.cloud_diagnostics_tuple.cf,
         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,8 @@ function precomputed_quantities(Y, atmos)
     )
     cloud_diagnostics_tuple =
         similar(Y.c, @NamedTuple{cf::FT, q_liq::FT, q_ice::FT})
+    # Cloud fraction is used to calculate buoyancy gradient, so we initialize it to 0 here.
+    @. 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))),
@@ -200,7 +201,6 @@ function precomputed_quantities(Y, atmos)
             ᶜentrʲs = similar(Y.c, NTuple{n, FT}),
             ᶜdetrʲs = similar(Y.c, NTuple{n, FT}),
             ᶜturb_entrʲs = similar(Y.c, NTuple{n, FT}),
-            ᶜgradᵥ_θ_virt = Fields.Field(C3{FT}, cspace),
             ᶜgradᵥ_q_tot = Fields.Field(C3{FT}, cspace),
             ᶜgradᵥ_θ_liq_ice = Fields.Field(C3{FT}, cspace),
             ᶠnh_pressure₃_buoyʲs = similar(Y.f, NTuple{n, C3{FT}}),
@@ -212,7 +212,6 @@ function precomputed_quantities(Y, atmos)
         (; ρatke_flux = similar(Fields.level(Y.f, half), C3{FT}),) : (;)
 
     sgs_quantities = (;
-        ᶜgradᵥ_θ_virt = Fields.Field(C3{FT}, cspace),
         ᶜgradᵥ_q_tot = Fields.Field(C3{FT}, cspace),
         ᶜgradᵥ_θ_liq_ice = Fields.Field(C3{FT}, cspace),
     )
@@ -548,8 +547,6 @@ NVTX.@annotate function set_explicit_precomputed_quantities_part1!(Y, p, t)
     end
 
     if turbconv_model isa AbstractEDMF
-        @. p.precomputed.ᶜgradᵥ_θ_virt =
-            ᶜgradᵥ(ᶠinterp(TD.virtual_pottemp(thermo_params, ᶜts)))
         @. p.precomputed.ᶜgradᵥ_q_tot =
             ᶜgradᵥ(ᶠinterp(TD.total_specific_humidity(thermo_params, ᶜts)))
         @. p.precomputed.ᶜgradᵥ_θ_liq_ice =
@@ -568,6 +565,15 @@ NVTX.@annotate function set_explicit_precomputed_quantities_part2!(Y, p, t)
     (; turbconv_model, moisture_model, cloud_model) = p.atmos
     (; call_cloud_diagnostics_per_stage) = p.atmos
 
+    # The buoyancy gradient depends on the cloud fraction, and the cloud fraction 
+    # depends on the mixing length, which depends on the buoyancy gradient. 
+    # We break this circular dependency by using cloud fraction from the previous time step in the 
+    # buoyancy gradient calculation. This breaks reproducible restart in general, 
+    # but we support reproducible restart with grid-scale cloud by recalculating the cloud fraction here.
+    if p.atmos.cloud_model isa GridScaleCloud
+        set_cloud_fraction!(Y, p, p.atmos.moisture_model, p.atmos.cloud_model)
+    end
+
     if turbconv_model isa PrognosticEDMFX
         set_prognostic_edmf_precomputed_quantities_explicit_closures!(Y, p, t)
         set_prognostic_edmf_precomputed_quantities_precipitation!(

src/cache/prognostic_edmf_precomputed_quantities.jl

@@ -299,14 +299,15 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
         )
     end
 
-    (; ᶜgradᵥ_θ_virt, ᶜgradᵥ_q_tot, ᶜgradᵥ_θ_liq_ice) = p.precomputed
+    (; ᶜgradᵥ_q_tot, ᶜgradᵥ_θ_liq_ice, cloud_diagnostics_tuple) =
+        p.precomputed
     @. ᶜ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/callbacks/callbacks.jl

@@ -106,16 +106,12 @@ function external_driven_single_column!(integrator)
     evaluate!(ᶜls_subsidence, wa, t)
 end
 
-
-
 NVTX.@annotate function cloud_fraction_model_callback!(integrator)
     Y = integrator.u
     p = integrator.p
-    (; ᶜts, ᶜgradᵥ_θ_virt, ᶜgradᵥ_q_tot, ᶜgradᵥ_θ_liq_ice) = p.precomputed
+    (; ᶜts, ᶜgradᵥ_q_tot, ᶜgradᵥ_θ_liq_ice) = p.precomputed
     thermo_params = CAP.thermodynamics_params(p.params)
     if isnothing(p.atmos.turbconv_model)
-        @. ᶜgradᵥ_θ_virt =
-            ᶜgradᵥ(ᶠinterp(TD.virtual_pottemp(thermo_params, ᶜts)))
         @. ᶜgradᵥ_q_tot =
             ᶜgradᵥ(ᶠinterp(TD.total_specific_humidity(thermo_params, ᶜts)))
         @. ᶜgradᵥ_θ_liq_ice =

src/prognostic_equations/eddy_diffusion_closures.jl

@@ -15,9 +15,8 @@ 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
+        ∂qt∂z::AbstractField,   # Vertical gradient of total specific humidity
+        ∂θli∂z::AbstractField,   # Vertical gradient of liquid-ice potential temperature
         local_geometry::Fields.LocalGeometry,
         # Argument for the second method (internal use with precomputed EnvBuoyGradVars):
         # bg_model::EnvBuoyGradVars
@@ -50,9 +49,8 @@ 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.
+- `∂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.
 
@@ -66,10 +64,10 @@ function buoyancy_gradients(
     thermo_params,
     moisture_model,
     ts,
+    cf,
     ::Type{C3},
-    ∂θv∂z_unsat,
-    ∂qt∂z_sat,
-    ∂θli∂z_sat,
+    ∂qt∂z,
+    ∂θli∂z,
     ᶜlg,
 ) where {C3}
     return buoyancy_gradients(
@@ -78,11 +76,11 @@ function buoyancy_gradients(
         moisture_model,
         EnvBuoyGradVars(
             ts,
+            cf,
             projected_vector_buoy_grad_vars(
                 C3,
-                ∂θv∂z_unsat,
-                ∂qt∂z_sat,
-                ∂θli∂z_sat,
+                ∂qt∂z,
+                ∂θli∂z,
                 ᶜlg,
             ),
         ),
@@ -99,36 +97,36 @@ 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) * Π
+    g_over_θv = g / TD.virtual_pottemp(thermo_params, ts)
 
-    t_sat = TD.air_temperature(thermo_params, ts)
+    T = TD.air_temperature(thermo_params, ts)
     λ = TD.liquid_fraction(thermo_params, ts)
     lh =
-        λ * TD.latent_heat_vapor(thermo_params, t_sat) +
-        (1 - λ) * TD.latent_heat_sublim(thermo_params, t_sat)
+        λ * TD.latent_heat_vapor(thermo_params, T) +
+        (1 - λ) * TD.latent_heat_sublim(thermo_params, T)
     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)
+    q_sat = TD.q_vap_saturation(thermo_params, ts)
+    q_tot = TD.total_specific_humidity(thermo_params, ts)
+    θ = TD.dry_pottemp(thermo_params, ts)
+    ∂b∂θli_unsat = g_over_θv * (1 + (Rv_over_Rd - 1) * q_tot)
+    ∂b∂qt_unsat = g_over_θ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)
+        g_over_θv *
+        (1 + Rv_over_Rd * (1 + lh / R_v / T) * q_sat - q_tot) /
+        (1 + lh^2 / cp_m / R_v / T^2 * q_sat)
     )
     ∂b∂qt_sat =
-        (lh / cp_m / t_sat * ∂b∂θli_sat - ∂b∂θv) *
-        TD.dry_pottemp(thermo_params, ts)
+        (lh / cp_m / T * ∂b∂θli_sat - g_over_θv) * θ
 
     ∂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,
     )
@@ -140,7 +138,6 @@ end
         closure::AbstractEnvBuoyGradClosure,
         bg_model::EnvBuoyGradVars,
         thermo_params,
-        ∂b∂θv::FT,
         ∂b∂θli_sat::FT,
         ∂b∂qt_sat::FT,
     ) where {FT}
@@ -165,7 +162,6 @@ Arguments:
 - `closure`: The environmental buoyancy gradient closure type.
 - `bg_model`: Precomputed environmental buoyancy gradient variables (`EnvBuoyGradVars`).
 - `thermo_params`: Thermodynamic parameters from `CLIMAParameters`.
-- `∂b∂θv`: Partial derivative of buoyancy w.r.t. virtual potential temperature (unsaturated part).
 - `∂b∂θli_sat`: Partial derivative of buoyancy w.r.t. liquid-ice potential temperature (saturated part).
 - `∂b∂qt_sat`: Partial derivative of buoyancy w.r.t. total specific humidity (saturated part).
 
@@ -176,18 +172,21 @@ function buoyancy_gradient_chain_rule(
     ::AbstractEnvBuoyGradClosure,
     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_θl_sat = ∂b∂θli_sat * bg_model.∂θli∂z_sat
-    ∂b∂z_qt_sat = ∂b∂qt_sat * bg_model.∂qt∂z_sat
+    ∂b∂z_θli_unsat = ∂b∂θli_unsat * bg_model.∂θli∂z
+    ∂b∂z_qt_unsat = ∂b∂qt_unsat * bg_model.∂qt∂z
+    ∂b∂z_unsat = ∂b∂z_θli_unsat + ∂b∂z_qt_unsat
+    ∂b∂z_θl_sat = ∂b∂θli_sat * bg_model.∂θli∂z
+    ∂b∂z_qt_sat = ∂b∂qt_sat * bg_model.∂qt∂z
     ∂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

@@ -67,17 +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) = p.precomputed
+    (; ᶜts, ᶠu³, ᶜlinear_buoygrad, ᶜstrain_rate_norm, cloud_diagnostics_tuple) =
+        p.precomputed
 
     @. ᶜlinear_buoygrad = buoyancy_gradients(
         BuoyGradMean(),
         thermo_params,
         p.atmos.moisture_model,
         ᶜts,
+        cloud_diagnostics_tuple.cf,
         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
-    ∂qt∂z_sat::FT
-    ∂θli∂z_sat::FT
+    cf::FT
+    ∂qt∂z::FT
+    ∂θli∂z::FT
 end
 
 function EnvBuoyGradVars(
     ts::TD.ThermodynamicState,
-    ∂θv∂z_unsat_∂qt∂z_sat_∂θli∂z_sat,
+    cf,
+    ∂qt∂z_∂θli∂z,
 )
-    (; ∂θ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, ∂θli∂z) = ∂qt∂z_∂θli∂z
+    return EnvBuoyGradVars(ts, cf, ∂qt∂z, ∂θli∂z)
 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 = V(v1, lg)[1] / ubvd,
+        ∂θli∂z = V(v2, lg)[1] / ubvd,
     )
 end
 

test/restart.jl

@@ -456,6 +456,7 @@ else
             "output_dir" => joinpath(amip_output_loc, amip_job_id),
             "dt_cloud_fraction" => "1secs",
             "rad" => "allskywithclear",
+            "cloud_model" => "grid_scale",
             "toml" => [
                 joinpath(@__DIR__, "../toml/longrun_aquaplanet_diagedmf.toml"),
             ],