add precipitation fraction

src/cache/microphysics_cache.jl

@@ -936,15 +936,17 @@ function set_microphysics_tendency_cache!(
             ᶜT, ᶜw_air, cmp, thp, dt, nsubs,
         )
     else
-        (; ᶜT′T′, ᶜq′q′, ᶜcorr_Tq, ᶜsgs_moments) = p.precomputed
-        α = sgs_variance_fidelity(CAP.cloud_fraction_steepness_scale(p.params))
-        ξ_liq = CAP.sgs_liquid_uniform_fraction(p.params)
-        ξ_ice = CAP.sgs_ice_uniform_fraction(p.params)
+        (; ᶜT′T′, ᶜq′q′, ᶜcorr_Tq, ᶜsgs_moments, ᶜprecip_frac) = p.precomputed
+        # One packed scalar-options argument keeps the broadcast short enough
+        # for the GPU kernel (see `sgs_microphysics_options`).
+        opts = sgs_microphysics_options(p.params)
         @. ᶜmp_tendency = microphysics_tendencies_1m(
             BMT.Microphysics1Moment(), sgs_quad, cmp, thp, Y.c.ρ, ᶜT, ᶜw_air,
             ᶜq_tot_nonneg, ᶜq_lcl, ᶜq_icl, ᶜq_rai, ᶜq_sno,
-            ᶜT′T′, ᶜq′q′, ᶜcorr_Tq, ᶜsgs_moments.λ_lagrange, α, ξ_liq, ξ_ice,
-            dt, nsubs_quad,
+            ᶜT′T′, ᶜq′q′, ᶜcorr_Tq, ᶜsgs_moments.λ_lagrange, dt, nsubs_quad,
+            TD.liquid_fraction(thp, ᶜT, max(0, ᶜq_lcl), max(0, ᶜq_icl)),
+            ᶜq_tot_nonneg - TD.q_vap_saturation(thp, ᶜT, Y.c.ρ),
+            ᶜprecip_frac, ᶜsgs_moments.sigma_S, ᶜsgs_moments.CF_d, opts,
         )
     end
 
@@ -1002,10 +1004,10 @@ function set_microphysics_tendency_cache!(
             ᶜT⁰, ᶜw⁰_air, cmp, thp, dt, nsubs,
         )
     else
-        (; ᶜT′T′, ᶜq′q′, ᶜcorr_Tq, ᶜsgs_moments) = p.precomputed
-        α = sgs_variance_fidelity(CAP.cloud_fraction_steepness_scale(p.params))
-        ξ_liq = CAP.sgs_liquid_uniform_fraction(p.params)
-        ξ_ice = CAP.sgs_ice_uniform_fraction(p.params)
+        (; ᶜT′T′, ᶜq′q′, ᶜcorr_Tq, ᶜsgs_moments, ᶜprecip_frac) = p.precomputed
+        # One packed scalar-options argument keeps the broadcast short enough
+        # for the GPU kernel (see `sgs_microphysics_options`).
+        opts = sgs_microphysics_options(p.params)
         # The liquid fraction `λ` and the linearized SGS saturation-excess mean
         # `mu_S` are held fixed across the quadrature (they depend only on the mean
         # state), so compute them once here and pass them in, instead of recomputing
@@ -1017,8 +1019,9 @@ function set_microphysics_tendency_cache!(
         @. ᶜmp_tendency⁰ = microphysics_tendencies_1m(
             BMT.Microphysics1Moment(), sgs_quad, cmp, thp, ᶜρ⁰, ᶜT⁰, ᶜw⁰_air,
             ᶜq_tot_nonneg⁰, ᶜq_lcl⁰, ᶜq_icl⁰, ᶜq_rai⁰, ᶜq_sno⁰,
-            ᶜT′T′, ᶜq′q′, ᶜcorr_Tq, ᶜsgs_moments.λ_lagrange, α, ξ_liq, ξ_ice,
-            dt, nsubs_quad, ᶜλ⁰, ᶜmu_S⁰,
+            ᶜT′T′, ᶜq′q′, ᶜcorr_Tq, ᶜsgs_moments.λ_lagrange, dt, nsubs_quad,
+            ᶜλ⁰, ᶜmu_S⁰, ᶜprecip_frac, ᶜsgs_moments.sigma_S, ᶜsgs_moments.CF_d,
+            opts,
         )
     end
 

src/cache/precomputed_quantities.jl

@@ -280,12 +280,16 @@ function precomputed_quantities(Y, atmos)
     uses_microphysics_quadrature_moments =
         atmos.microphysics_model isa
         Union{NonEquilibriumMicrophysics1M, NonEquilibriumMicrophysics2M}
-    # `ᶜsgs_moments` caches `(sigma_S, λ_lagrange)` — the SGS standard
-    # deviation and the Lagrange multiplier used by `Microphysics1MEvaluator`.
-    #  Allocated only for 1M/2M schemes.
+    # `ᶜsgs_moments` caches `(sigma_S, λ_lagrange, CF_d)` — the SGS standard
+    # deviation, the Lagrange multiplier used by `Microphysics1MEvaluator`,
+    # and the discrete cloudy mass of the quadrature measure
+    # (`_discrete_cloud_fraction`). Allocated only for 1M/2M schemes, together
+    # with `ᶜprecip_frac`, the overlap precipitation fraction
+    # (`set_precip_fraction!`).
     SGSMomentsNT = @NamedTuple{
         sigma_S::FT,
         λ_lagrange::FT,
+        CF_d::FT,
     }
     covariance_quantities = if uses_sgs_quadrature
         base = (;
@@ -312,8 +316,11 @@ function precomputed_quantities(Y, atmos)
             )...,
         )
         uses_microphysics_quadrature_moments ?
-        (; base..., ᶜsgs_moments = similar(Y.c, SGSMomentsNT)) :
-        base
+        (;
+            base...,
+            ᶜsgs_moments = similar(Y.c, SGSMomentsNT),
+            ᶜprecip_frac = zeros(axes(Y.c)),
+        ) : base
     else
         (;)
     end

src/parameterized_tendencies/microphysics/cloud_fraction.jl

@@ -1062,10 +1062,73 @@ uses the same quadrature points.
     return λ
 end
 
+"""
+    discrete_cloudy_weight_width_coeff(FT)
+
+Smoothing width of the discrete cloudy weight `sᵢ`, in units of the
+equilibrium PDF width `α·σ_S` (see `_discrete_cloud_fraction`). A hard
+indicator `1[shifted_excess > 0]` would make each point's cloudy/clear
+assignment jump as it crosses the cloud threshold between steps; `0.25` keeps
+the transition narrow compared with the Gauss–Hermite point spacing
+(≈ 1.7 σ_S at order 3), so `CF_d` stays close to the hard count.
+"""
+@inline discrete_cloudy_weight_width_coeff(::Type{FT}) where {FT} = FT(0.25)
+
+"""
+    discrete_cloudy_weight_width(α, sigma_S)
+
+Smoothing width `ε_w = c_w·α·σ_S` of the discrete cloudy weight, floored at
+`ϵ_numerics` so the zero-variance limit stays finite.
+"""
+@inline discrete_cloudy_weight_width(α, sigma_S) = max(
+    discrete_cloudy_weight_width_coeff(typeof(sigma_S)) * α * sigma_S,
+    ϵ_numerics(typeof(sigma_S)),
+)
+
+"""
+    discrete_cloudy_weight(shifted_excess, ε_w)
+
+Smooth cloudy weight `s = sigmoid(shifted_excess / ε_w) ∈ [0, 1]` of one
+quadrature point, where `shifted_excess = λ_lagrange + α·S′` is the signed
+quantity whose positive part is the local condensate. Evaluated in `tanh`
+form, which cannot overflow.
+"""
+@inline discrete_cloudy_weight(shifted_excess, ε_w) =
+    (1 + tanh(shifted_excess / (2 * ε_w))) / 2
+
+"""
+    _discrete_cloud_fraction(q_c, λ_lagrange, α, sigma_S, S′s, ws)
+
+Discrete cloudy mass of the quadrature measure,
+
+    shifted_excessᵢ = λ_lagrange + α·S′ᵢ
+    sᵢ              = sigmoid(shifted_excessᵢ / ε_w),  ε_w = c_w·α·σ_S
+    CF_d            = Σᵢ wᵢ·sᵢ
+
+a smoothed count of the quadrature points that carry condensate under exactly
+the measure `Microphysics1MEvaluator` integrates over. It is the cloud cover of
+the sampled PDF itself, without the augmented-σ floor of `ᶜcloud_fraction`, and
+is used as a diagnostic of that PDF. Condensate-free cells (`q_c ≤ 0`) return
+`0`: there the fit parks `λ_lagrange` exactly on the largest kink, so the
+moistest point sits on the cloud threshold and the smooth weight would
+otherwise report half its quadrature weight as cloudy in every clear cell.
+"""
+@inline function _discrete_cloud_fraction(q_c, λ_lagrange, α, sigma_S, S′s, ws)
+    FT = typeof(λ_lagrange)
+    q_c <= zero(FT) && return zero(FT)
+    ε_w = discrete_cloudy_weight_width(α, sigma_S)
+    CF_d = zero(FT)
+    @inbounds for i in eachindex(S′s)
+        shifted_excess = λ_lagrange + α * S′s[i]
+        CF_d += ws[i] * discrete_cloudy_weight(shifted_excess, ε_w)
+    end
+    return CF_d
+end
+
 """
     _compute_sgs_moments(thp, ρ, T, q_tot, q_c, sgs_quad, T′T′, q′q′, corr_Tq, α)
 
-Single quadrature pass returning `(sigma_S, λ_lagrange)`:
+Single quadrature pass returning `(sigma_S, λ_lagrange, CF_d)`:
 
   - `sigma_S = sqrt(Σᵢ wᵢ·S′ᵢ²)`: SGS standard deviation of the sampled
     centred excess, clipped at `ϵ_numerics(FT)`.
@@ -1074,22 +1137,29 @@ Single quadrature pass returning `(sigma_S, λ_lagrange)`:
     quadrature measure — the same points and weights the microphysics
     evaluator integrates over (see `_fit_discrete_lagrange`; the analytic
     truncated-Gaussian inverse `_compute_z` provides the seed).
+  - `CF_d = Σᵢ wᵢ·sᵢ`: the discrete cloudy mass of the same measure, with `sᵢ`
+    a smoothed cloudy indicator (see `_discrete_cloud_fraction`), accumulated
+    over points the pass already visits.
 
 The SGS mean `μ_S = q_tot − q_sat(T, ρ)` is analytic under the closure's
 linearization (see `_sgs_saturation_moments`) and is recomputed on demand
 wherever it is needed downstream.
 
 Without SGS sampling (`nothing`, `GridMeanSGS`), all mass sits at the mean:
 `sigma_S = ϵ_numerics(FT)` and the constraint gives `λ_lagrange = q_c`
-directly, matching the σ_S → 0 limit of the sampled branch.
+directly, matching the σ_S → 0 limit of the sampled branch; `CF_d` is then
+the all-mass-at-the-mean limit `q_c > 0 ? 1 : 0`.
 """
 @inline function _compute_sgs_moments(
     thp, ρ, T, q_tot, q_c,
     sgs_quad, T′T′, q′q′, corr_Tq, α,
 )
     FT = typeof(ρ)
-    not_quadrature(sgs_quad) &&
-        return (; sigma_S = ϵ_numerics(FT), λ_lagrange = q_c)
+    not_quadrature(sgs_quad) && return (;
+        sigma_S = ϵ_numerics(FT),
+        λ_lagrange = q_c,
+        CF_d = ifelse(q_c > zero(FT), one(FT), zero(FT)),
+    )
 
     mu_S = q_tot - TD.q_vap_saturation(thp, T, ρ)
     transform =
@@ -1104,7 +1174,8 @@ directly, matching the σ_S → 0 limit of the sampled branch.
     σ_S_eff = α * sigma_S
     λ0 = _compute_z(q_c / σ_S_eff) * σ_S_eff
     λ_lagrange = _fit_discrete_lagrange(λ0, q_c, α, S′s, ws)
-    return (; sigma_S, λ_lagrange)
+    CF_d = _discrete_cloud_fraction(q_c, λ_lagrange, α, sigma_S, S′s, ws)
+    return (; sigma_S, λ_lagrange, CF_d)
 end
 
 """
@@ -1113,10 +1184,12 @@ end
 Final post-Aitken update. No-op when `ᶜsgs_moments` is not allocated (dry / 0M).
 
 Uses ONE quadrature pass via `_compute_sgs_moments` to fill
-`ᶜsgs_moments = (sigma_S, λ_lagrange)`, then computes
+`ᶜsgs_moments = (sigma_S, λ_lagrange, CF_d)`, then computes
 `ᶜcloud_fraction` consistently with the augmented `σ_aug` closure (see
 `_compute_cloud_fraction`) from the grid-mean cloud condensate
-(`_grid_mean_cloud_condensate`).
+(`_grid_mean_cloud_condensate`). Finally runs `set_precip_fraction!`, the
+column overlap sweep that turns the fresh cover into the precipitation
+fraction `ᶜprecip_frac`.
 
 Overwrites `p.scratch.ᶜtemp_scalar`, `ᶜtemp_scalar_2`, `ᶜtemp_scalar_3`,
 `ᶜtemp_scalar_5`, and `ᶜtemp_scalar_6`. It must not touch `ᶜtemp_scalar_4` or
@@ -1173,7 +1246,7 @@ NVTX.@annotate function set_sgs_moments_and_cloud_fraction!(Y, p)
     ᶜq′q′,
     ᶜcorr_Tq
 
-        # ONE quadrature pass → (sigma_S, λ_lagrange).
+        # ONE quadrature pass → (sigma_S, λ_lagrange, CF_d).
         @. ᶜsgs_moments = _compute_sgs_moments(
             thermo_params, ᶜρ_env, ᶜT_mean, ᶜq_mean, ᶜq_lcl + ᶜq_icl,
             $(sgs_quad), ᶜT′T′, ᶜq′q′, ᶜcorr_Tq, α_ft,
@@ -1194,6 +1267,122 @@ NVTX.@annotate function set_sgs_moments_and_cloud_fraction!(Y, p)
             $(floor),
         )
     end
+    set_precip_fraction!(Y, p)
+    return nothing
+end
+
+"""
+    set_precip_fraction!(Y, p)
+
+Fill `p.precomputed.ᶜprecip_frac` with the precipitation fraction `a_p`, the
+area fraction of the cell that precipitation falls through, by maximum-random
+overlap of the cloud cover swept from the model top down:
+
+    a_p(k) = max(CF(k), f_decay · a_p(k+1))   where q_rai + q_sno > q_min
+    a_p(k) = 0                                 otherwise
+
+with `f_decay = sgs_precip_overlap_decay` (`1` is pure maximum overlap; the
+shrink is applied per level). A precipitation-free level resets the recursion,
+so a shaft that evaporates completely does not seed the layers below it. The
+presence threshold is Thermodynamics' `q_min`, the one `TD.has_condensate`
+uses, so the mask closes where real precipitation ends rather than on
+numerical residue.
+
+The recursion is seeded with `ᶜcloud_fraction`, the single-domain cover of the
+cell (environment PDF plus updraft condensate, with the cover floor), because
+`a_p` describes the geometry of the cloud the precipitation came from, which
+includes the updraft; the precipitation mask uses the precipitation of the
+domain the quadrature integrates (the environment under `PrognosticEDMFX`).
+`a_p` enters the quadrature only as the fraction of the PDF that carries the
+shaft (`sgs_precip_shaft_threshold`); the per-point assignment is normalized
+by the discrete weight of the selected points, so no measure consistency
+between the cover and the quadrature is required for conservation.
+
+A negative `sgs_precip_overlap_decay` disables the overlap fraction:
+`ᶜprecip_frac` is set to `1` and the node placement runs in its moist-half
+mode (`Microphysics1MEvaluator`).
+
+Mutates `p.precomputed.ᶜprecip_frac`; the return value is unused.
+"""
+NVTX.@annotate function set_precip_fraction!(Y, p)
+    hasproperty(p.precomputed, :ᶜprecip_frac) || return nothing
+    FT = eltype(p.params)
+    f_decay = FT(CAP.sgs_precip_overlap_decay(p.params))
+    if f_decay < zero(FT)
+        @. p.precomputed.ᶜprecip_frac = one(FT)
+        return nothing
+    end
+    thermo_params = CAP.thermodynamics_params(p.params)
+    _precip_fraction_sweep!(
+        p.precomputed.ᶜprecip_frac,
+        p.precomputed.ᶜcloud_fraction,
+        _get_precip_mean(Y, p, p.atmos.turbconv_model),
+        f_decay,
+        FT(TD.Parameters.q_min(thermo_params)),
+    )
+    return nothing
+end
+
+"""
+    _precip_fraction_sweep!(ᶜprecip_frac, ᶜcf, ᶜq_precip, f_decay, q_precip_min)
+
+Run the maximum-random overlap recursion of `set_precip_fraction!` on plain
+fields, writing `a_p` into `ᶜprecip_frac`. Levels with
+`ᶜq_precip ≤ q_precip_min` carry no shaft.
+
+Implemented as a loop of level broadcasts rather than with
+`Operators.column_accumulate!`, whose CPU path materializes a level `Field`
+per level per column; the loop allocates per level independent of the number
+of columns and gives the same result. Its cost on GPU is one kernel launch per
+level.
+
+The precipitation mask is folded into `ᶜprecip_frac` first, as a negative
+sentinel, so the recursion reads and writes a single field: a level holding no
+precipitation must both report `a_p = 0` and stop the shaft above it from
+being inherited further down, so "no precipitation" stays distinguishable from
+"precipitating, but with no cloud of its own" (`cf = 0`, which does inherit).
+"""
+function _precip_fraction_sweep!(ᶜprecip_frac, ᶜcf, ᶜq_precip, f_decay, q_precip_min)
+    FT = eltype(ᶜprecip_frac)
+    @. ᶜprecip_frac = ifelse(ᶜq_precip > q_precip_min, ᶜcf, -one(FT))
+    # Level 1 is the bottom model level, so the recursion runs from the last
+    # level down. The top level has no shaft above it, so its sentinel just
+    # resolves to zero.
+    nz = Spaces.nlevels(axes(ᶜprecip_frac))
+    ᶜa_p_above = Fields.level(ᶜprecip_frac, nz)
+    @. ᶜa_p_above = max(ᶜa_p_above, zero(FT))
+    for level in (nz - 1):-1:1
+        ᶜa_p = Fields.level(ᶜprecip_frac, level)
+        @. ᶜa_p = ifelse(
+            ᶜa_p < zero(FT),
+            zero(FT),
+            min(one(FT), max(ᶜa_p, f_decay * ᶜa_p_above)),
+        )
+        ᶜa_p_above = ᶜa_p
+    end
+    return nothing
+end
+
+"""
+    _get_precip_mean(Y, p, turbconv_model)
+
+Mean precipitation specific humidity `q_rai + q_sno` of the domain that carries
+the SGS closure: the environment for `PrognosticEDMFX`, the grid mean
+otherwise. Matches the domain of `_get_condensate_means`, so the precipitation
+mask of `set_precip_fraction!` describes the air the quadrature integrates.
+"""
+function _get_precip_mean(Y, p, turbconv_model)
+    if turbconv_model isa PrognosticEDMFX
+        return @. lazy(
+            max(0, $(ᶜspecific_env_value(@name(q_rai), Y, p))) +
+            max(0, $(ᶜspecific_env_value(@name(q_sno), Y, p))),
+        )
+    else
+        return @. lazy(
+            max(0, specific(Y.c.ρq_rai, Y.c.ρ)) +
+            max(0, specific(Y.c.ρq_sno, Y.c.ρ)),
+        )
+    end
 end