fix: 2M tendencies, params

- rename returned NamedTuple for rain_evaporation
    tendency to (; ∂ₜρn_rai, ∂ₜq_rai)
- add 2M convenience constructors, update docs
- update tests accordingly

Author: haakon-e

docs/src/plots/RainEvapoartionSB2006.jl

@@ -7,49 +7,53 @@ import CloudMicrophysics.ThermodynamicsInterface as TDI
 import CloudMicrophysics.Parameters as CMP
 import CloudMicrophysics.Common as CO
 import CloudMicrophysics.Microphysics2M as CM2
+import CloudMicrophysics.Utilities as UT
 
 FT = Float64
 
-const tps = TDI.TD.Parameters.ThermodynamicsParameters(FT)
-const aps = CMP.AirProperties(FT)
-const SB2006 = CMP.SB2006(FT)
+tps = TDI.TD.Parameters.ThermodynamicsParameters(FT)
+aps = CMP.AirProperties(FT)
+SB2006 = CMP.SB2006(FT)
 
 function rain_evaporation_CPU(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Nᵣ, T)
+    ϵₘ, ϵₙ = UT.ϵ_numberics_2M_MN_tuple(FT)
 
-    evap_rate_0 = FT(0)
-    evap_rate_1 = FT(0)
+    ∂ₜρn_rai = zero(qᵣ)
+    ∂ₜq_rai = zero(qᵣ)
     S = TDI.supersaturation_over_liquid(tps, qₜ, qₗ + qᵣ, qᵢ + qₛ, ρ, T)
 
-    if (qᵣ > FT(0) && S < FT(0))
+    # If there are no raindrops or conditions are supersaturated, no evaporation occurs
+    (Nᵣ ≤ ϵₙ || S ≥ 0) && return (; ∂ₜρn_rai, ∂ₜq_rai)
 
-        (; ν_air, D_vapor) = aps
-        (; av, bv, α, β, ρ0) = SB2006.evap
-        ρw = SB2006.pdf_r.ρw
-        x_star = SB2006.pdf_r.xr_min
-        G = CO.G_func_liquid(aps, tps, T)
+    (; ν_air, D_vapor) = aps
+    (; av, bv, α, β, ρ0) = SB2006.evap
+    ρw = SB2006.pdf_r.ρw
+    x_star = SB2006.pdf_r.xr_min
+    G = CO.G_func_liquid(aps, tps, T)
 
-        (; xr_mean) = CM2.pdf_rain_parameters(SB2006.pdf_r, qᵣ, ρ, Nᵣ)
-        Dr = (FT(6) / FT(π) / ρw)^FT(1 / 3) * xr_mean^FT(1 / 3)
+    (; xr_mean) = CM2.pdf_rain_parameters(SB2006.pdf_r, qᵣ, ρ, Nᵣ)
+    Dr = cbrt(6 * xr_mean / (π * ρw))
 
-        t_star = (FT(6) * x_star / xr_mean)^FT(1 / 3)
-        a_vent_0 = av * FT(SF.gamma(-1, t_star)) / FT(6)^FT(-2 / 3)
-        b_vent_0 =
-            bv * FT(SF.gamma((-1 / 2) + (3 / 2) * β, t_star)) /
-            FT(6)^FT(β / 2 - 1 / 2)
+    t_star = cbrt(6 * x_star / xr_mean)
+    a_vent_0 = av * SF.gamma(-1, t_star) / FT(6)^(-2 // 3)
+    b_vent_0 = bv * SF.gamma((-1 // 2) + 3 // 2 * β, t_star) / FT(6)^(β // 2 - 1 // 2)
 
-        a_vent_1 = av * SF.gamma(FT(2)) / FT(6)^FT(1 / 3)
-        b_vent_1 =
-            bv * SF.gamma(FT(5 / 2) + FT(3 / 2) * β) / FT(6)^FT(β / 2 + 1 / 2)
+    a_vent_1 = av * SF.gamma(FT(2)) / cbrt(FT(6))
+    b_vent_1 = bv * SF.gamma(5 // 2 + 3 // 2 * β) / FT(6)^(β // 2 + 1 // 2)
 
-        N_Re = α * xr_mean^β * sqrt(ρ0 / ρ) * Dr / ν_air
-        Fv0 = a_vent_0 + b_vent_0 * (ν_air / D_vapor)^FT(1 / 3) * sqrt(N_Re)
-        Fv1 = a_vent_1 + b_vent_1 * (ν_air / D_vapor)^FT(1 / 3) * sqrt(N_Re)
+    N_Re = α * xr_mean^β * sqrt(ρ0 / ρ) * Dr / ν_air
+    Fv0 = a_vent_0 + b_vent_0 * cbrt(ν_air / D_vapor) * sqrt(N_Re)
+    Fv1 = a_vent_1 + b_vent_1 * cbrt(ν_air / D_vapor) * sqrt(N_Re)
 
-        evap_rate_0 = min(FT(0), FT(2) * FT(π) * G * S * Nᵣ * Dr * Fv0 / xr_mean)
-        evap_rate_1 = min(FT(0), FT(2) * FT(π) * G * S * Nᵣ * Dr * Fv1 / ρ)
-    end
+    ∂ₜρn_rai = min(0, 2 * FT(π) * G * S * Nᵣ * Dr * Fv0 / xr_mean)
+    ∂ₜq_rai = min(0, 2 * FT(π) * G * S * Nᵣ * Dr * Fv1 / ρ)
 
-    return (; evap_rate_0, evap_rate_1)
+    # When xr = 0 ∂ₜρn_rai becomes NaN. We replace NaN with 0 which is the limit of
+    # ∂ₜρn_rai for xr -> 0.
+    ∂ₜρn_rai = ifelse(xr_mean / x_star < eps(FT), FT(0), ∂ₜρn_rai)
+    ∂ₜq_rai = ifelse(qᵣ < ϵₘ, FT(0), ∂ₜq_rai)
+
+    return (; ∂ₜρn_rai, ∂ₜq_rai)
 end
 
 qᵥ = FT(1e-2)
@@ -67,21 +71,21 @@ Nᵣ_range = range(1e6, stop = 1e9, length = 1000)
 T_range = range(273.15, stop = 273.15 + 50, length = 1000)
 
 #! format: off
-evap_qᵣ_0 = [rain_evaporation_CPU(SB2006, aps, tps, qₜ, qₗ, qᵢ, _qᵣ, qₛ, ρ, Nᵣ, T).evap_rate_0 for _qᵣ in qᵣ_range]
-evap_Nᵣ_0 = [rain_evaporation_CPU(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, _Nᵣ, T).evap_rate_0 for _Nᵣ in Nᵣ_range]
-evap_T_0 = [rain_evaporation_CPU(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Nᵣ, _T).evap_rate_0 for _T in T_range]
+evap_qᵣ_0 = [rain_evaporation_CPU(SB2006, aps, tps, qₜ, qₗ, qᵢ, _qᵣ, qₛ, ρ, Nᵣ, T).∂ₜρn_rai for _qᵣ in qᵣ_range]
+evap_Nᵣ_0 = [rain_evaporation_CPU(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, _Nᵣ, T).∂ₜρn_rai for _Nᵣ in Nᵣ_range]
+evap_T_0 = [rain_evaporation_CPU(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Nᵣ, _T).∂ₜρn_rai for _T in T_range]
 
-evap_qᵣ_0n = [CM2.rain_evaporation(SB2006, aps, tps, qₜ, qₗ, qᵢ, _qᵣ, qₛ, ρ, Nᵣ, T).evap_rate_0 for _qᵣ in qᵣ_range]
-evap_Nᵣ_0n = [CM2.rain_evaporation(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, _Nᵣ, T).evap_rate_0 for _Nᵣ in Nᵣ_range]
-evap_T_0n = [CM2.rain_evaporation(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Nᵣ, _T).evap_rate_0 for _T in T_range]
+evap_qᵣ_0n = [CM2.rain_evaporation(SB2006, aps, tps, qₜ, qₗ, qᵢ, _qᵣ, qₛ, ρ, Nᵣ, T).∂ₜρn_rai for _qᵣ in qᵣ_range]
+evap_Nᵣ_0n = [CM2.rain_evaporation(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, _Nᵣ, T).∂ₜρn_rai for _Nᵣ in Nᵣ_range]
+evap_T_0n = [CM2.rain_evaporation(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Nᵣ, _T).∂ₜρn_rai for _T in T_range]
 
-evap_qᵣ_3 = [rain_evaporation_CPU(SB2006, aps, tps, qₜ, qₗ, qᵢ, _qᵣ, qₛ, ρ, Nᵣ, T).evap_rate_1 for _qᵣ in qᵣ_range]
-evap_Nᵣ_3 = [rain_evaporation_CPU(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, _Nᵣ, T).evap_rate_1 for _Nᵣ in Nᵣ_range]
-evap_T_3 = [rain_evaporation_CPU(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Nᵣ, _T).evap_rate_1 for _T in T_range]
+evap_qᵣ_3 = [rain_evaporation_CPU(SB2006, aps, tps, qₜ, qₗ, qᵢ, _qᵣ, qₛ, ρ, Nᵣ, T).∂ₜq_rai for _qᵣ in qᵣ_range]
+evap_Nᵣ_3 = [rain_evaporation_CPU(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, _Nᵣ, T).∂ₜq_rai for _Nᵣ in Nᵣ_range]
+evap_T_3 = [rain_evaporation_CPU(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Nᵣ, _T).∂ₜq_rai for _T in T_range]
 
-evap_qᵣ_3n = [CM2.rain_evaporation(SB2006, aps, tps, qₜ, qₗ, qᵢ, _qᵣ, qₛ, ρ, Nᵣ, T).evap_rate_1 for _qᵣ in qᵣ_range]
-evap_Nᵣ_3n = [CM2.rain_evaporation(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, _Nᵣ, T).evap_rate_1 for _Nᵣ in Nᵣ_range]
-evap_T_3n = [CM2.rain_evaporation(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Nᵣ, _T).evap_rate_1 for _T in T_range]
+evap_qᵣ_3n = [CM2.rain_evaporation(SB2006, aps, tps, qₜ, qₗ, qᵢ, _qᵣ, qₛ, ρ, Nᵣ, T).∂ₜq_rai for _qᵣ in qᵣ_range]
+evap_Nᵣ_3n = [CM2.rain_evaporation(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, _Nᵣ, T).∂ₜq_rai for _Nᵣ in Nᵣ_range]
+evap_T_3n = [CM2.rain_evaporation(SB2006, aps, tps, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Nᵣ, _T).∂ₜq_rai for _T in T_range]
 
 fig = MK.Figure(resolution = (800, 600))
 

src/BulkMicrophysicsTendencies.jl

@@ -100,16 +100,8 @@ end
 
 """
     bulk_microphysics_tendencies(
-        ::Microphysics1Moment,
-        mp,
-        tps,
-        ρ,
-        T,
-        q_tot,
-        q_lcl,
-        q_icl,
-        q_rai,
-        q_sno,
+        ::Microphysics1Moment, mp, tps,
+        ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno,
     )
 
 Compute all 1-moment microphysics tendencies in one fused call.
@@ -156,17 +148,8 @@ This is a pure function of local thermodynamic state, suitable for:
   Limiters depend on timestep `dt` and should be applied by the caller after computing tendencies.
 """
 @inline function bulk_microphysics_tendencies(
-    ::Microphysics1Moment,
-    mp::CMP.Microphysics1MParams,
-    tps,
-    ρ,
-    T,
-    q_tot,
-    q_lcl,
-    q_icl,
-    q_rai,
-    q_sno,
-    N_lcl = zero(ρ),
+    ::Microphysics1Moment, mp::CMP.Microphysics1MParams, tps,
+    ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, N_lcl = zero(ρ),
 )
     # Clamp negative inputs to zero (robustness against numerical errors)
     ρ = UT.clamp_to_nonneg(ρ)
@@ -543,12 +526,8 @@ the second form uses the supersaturation threshold `S_0 * q_vap_sat`.
 - Does NOT include geopotential (caller adds Φ for energy tendency)
 """
 @inline function bulk_microphysics_tendencies(
-    ::Microphysics0Moment,
-    mp::CMP.Microphysics0MParams,
-    tps,
-    T,
-    q_lcl,
-    q_icl,
+    ::Microphysics0Moment, mp::CMP.Microphysics0MParams, tps,
+    T, q_lcl, q_icl,
 )
     q_lcl = UT.clamp_to_nonneg(q_lcl)
     q_icl = UT.clamp_to_nonneg(q_icl)
@@ -598,7 +577,7 @@ Used by both warm-only and warm+ice dispatch methods to reduce code duplication.
 - `dn_lcl_dt`: Cloud number tendency (1/kg/s)
 - `dn_rai_dt`: Rain number tendency (1/kg/s)
 """
-@inline function warm_rain_tendencies_2m(sb, q_lcl, q_rai, ρ, n_lcl, n_rai)
+@inline function warm_rain_tendencies_2m(sb, q_tot, q_lcl, q_rai, q_ice, ρ, n_lcl, n_rai)
     # Convert to number densities for CM2 functions
     N_lcl = ρ * n_lcl
     N_rai = ρ * n_rai
@@ -610,6 +589,17 @@ Used by both warm-only and warm+ice dispatch methods to reduce code duplication.
     dn_lcl_dt = zero(FT)
     dn_rai_dt = zero(FT)
 
+    # --- Condensation of vapor / evaporation of cloud liquid water ---
+    ∂ₜq_lcl_cond = CMNonEq.conv_q_vap_to_q_lcl_icl_MM2015(lcl, tps, q_tot, q_lcl, q_ice, q_rai, zero(q_ice), ρ, T)
+    ∂ₜn_lcl_cond = zero(∂ₜq_lcl_cond)  # neglect number change from condensation/evaporation
+    dq_lcl_dt += ∂ₜq_lcl_cond
+    dn_lcl_dt += ∂ₜn_lcl_cond
+
+    # --- Evaporation of rain ---
+    evap = CM2.rain_evaporation(sb, aps, tps, q_tot, q_lcl, q_ice, q_rai, zero(q_ice), ρ, N_rai, T)
+    dq_rai_dt += evap.∂ₜq_rai
+    dn_rai_dt += evap.∂ₜρn_rai / ρ
+
     # --- Autoconversion ---
     acnv = CM2.autoconversion(sb.acnv, sb.pdf_c, q_lcl, q_rai, ρ, N_lcl)
     dq_lcl_dt += acnv.dq_lcl_dt
@@ -635,6 +625,17 @@ Used by both warm-only and warm+ice dispatch methods to reduce code duplication.
     dn_rai_br = CM2.rain_breakup(sb.pdf_r, sb.brek, q_rai, ρ, N_rai, dn_rai_sc)
     dn_rai_dt += dn_rai_br / ρ
 
+    # --- Number adjustment for mass limits ---
+    # Cloud liquid
+    dn_lcl_inc = CM2.number_increase_for_mass_limit(sb.numadj, sb.pdf_c.xc_max, q_lcl, ρ, N_lcl)
+    dn_lcl_dec = CM2.number_decrease_for_mass_limit(sb.numadj, sb.pdf_c.xc_min, q_lcl, ρ, N_lcl)
+    dn_lcl_dt += (dn_lcl_inc + dn_lcl_dec) / ρ
+
+    # Rain
+    dn_rai_inc = CM2.number_increase_for_mass_limit(sb.numadj, sb.pdf_r.xr_max, q_rai, ρ, N_rai)
+    dn_rai_dec = CM2.number_decrease_for_mass_limit(sb.numadj, sb.pdf_r.xr_min, q_rai, ρ, N_rai)
+    dn_rai_dt += (dn_rai_inc + dn_rai_dec) / ρ
+
     return (; dq_lcl_dt, dq_rai_dt, dn_lcl_dt, dn_rai_dt)
 end
 
@@ -645,7 +646,7 @@ end
     bulk_microphysics_tendencies(
         ::Microphysics2Moment,
         mp::Microphysics2MParams{FT, WR, Nothing},
-        ...
+        ρ, T, q_tot, q_lcl, n_lcl, q_rai, n_rai,
     )
 
 Compute 2-moment **warm rain only** microphysics tendencies (Seifert-Beheng 2006).
@@ -675,21 +676,9 @@ For warm rain + P3 ice, see the method that accepts `Microphysics2MParams{FT, WR
 - `db_rim_dt`: Rime volume tendency (always zero for warm-only)
 """
 @inline function bulk_microphysics_tendencies(
-    ::Microphysics2Moment,
-    mp::CMP.Microphysics2MParams{FT, WR, Nothing},
-    tps,
-    ρ,
-    T,
-    q_tot,
-    q_lcl,
-    n_lcl,
-    q_rai,
-    n_rai,
-    q_ice = zero(ρ),
-    n_ice = zero(ρ),
-    q_rim = zero(ρ),
-    b_rim = zero(ρ),
-    logλ = zero(ρ),
+    ::Microphysics2Moment, mp::CMP.Microphysics2MParams{FT, WR, Nothing}, tps,
+    ρ, T, q_tot, q_lcl, n_lcl, q_rai, n_rai,
+    q_ice = zero(ρ), n_ice = zero(ρ), q_rim = zero(ρ), b_rim = zero(ρ), logλ = zero(ρ),
 ) where {FT, WR}
     # Clamp negative inputs to zero (robustness against numerical errors)
     ρ = UT.clamp_to_nonneg(ρ)
@@ -705,48 +694,28 @@ For warm rain + P3 ice, see the method that accepts `Microphysics2MParams{FT, WR
 
     # Unpack warm rain parameters
     sb = mp.warm_rain.seifert_beheng
-    aps = mp.warm_rain.air_properties
 
     # Initialize ice-related tendencies (always zero for warm-only)
     dq_ice_dt = zero(ρ)
     dq_rim_dt = zero(ρ)
     db_rim_dt = zero(ρ)
 
     # --- Core Warm Rain Processes (shared helper) ---
-    warm = warm_rain_tendencies_2m(sb, q_lcl, q_rai, ρ, n_lcl, n_rai)
+    warm = warm_rain_tendencies_2m(sb, q_tot, q_lcl, q_rai, q_ice, ρ, n_lcl, n_rai)
     dq_lcl_dt = warm.dq_lcl_dt
     dn_lcl_dt = warm.dn_lcl_dt
     dq_rai_dt = warm.dq_rai_dt
     dn_rai_dt = warm.dn_rai_dt
 
-    # Convert to number densities for remaining functions
-    N_lcl = ρ * n_lcl
-    N_rai = ρ * n_rai
-
-    # --- Rain evaporation ---
-    evap = CM2.rain_evaporation(sb, aps, tps, zero(ρ), q_lcl, zero(ρ), q_rai, zero(ρ), ρ, N_rai, T)
-    dq_rai_dt += evap.evap_rate_1
-    dn_rai_dt += evap.evap_rate_0 / ρ
-
-    # --- Number adjustment for mass limits ---
-    # Cloud liquid
-    dn_lcl_inc = CM2.number_increase_for_mass_limit(sb.numadj, sb.pdf_c.xc_max, q_lcl, ρ, N_lcl)
-    dn_lcl_dec = CM2.number_decrease_for_mass_limit(sb.numadj, sb.pdf_c.xc_min, q_lcl, ρ, N_lcl)
-    dn_lcl_dt += (dn_lcl_inc + dn_lcl_dec) / ρ
-
-    # Rain
-    dn_rai_inc = CM2.number_increase_for_mass_limit(sb.numadj, sb.pdf_r.xr_max, q_rai, ρ, N_rai)
-    dn_rai_dec = CM2.number_decrease_for_mass_limit(sb.numadj, sb.pdf_r.xr_min, q_rai, ρ, N_rai)
-    dn_rai_dt += (dn_rai_inc + dn_rai_dec) / ρ
-
     return (; dq_lcl_dt, dn_lcl_dt, dq_rai_dt, dn_rai_dt, dq_ice_dt, dq_rim_dt, db_rim_dt)
 end
 
 """
     bulk_microphysics_tendencies(
         ::Microphysics2Moment,
         mp::Microphysics2MParams{FT, WR, <:P3IceParams},
-        ...
+        ρ, T, q_tot, q_lcl, n_lcl, q_rai, n_rai,
+        q_ice, n_ice, q_rim, b_rim, logλ,
     )
 
 Compute 2-moment **warm rain + P3 ice** microphysics tendencies.
@@ -784,21 +753,9 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch.
 - `db_rim_dt`: Rime volume tendency (m³/kg/s)
 """
 @inline function bulk_microphysics_tendencies(
-    ::Microphysics2Moment,
-    mp::CMP.Microphysics2MParams{FT, WR, ICE},
-    tps,
-    ρ,
-    T,
-    q_tot,
-    q_lcl,
-    n_lcl,
-    q_rai,
-    n_rai,
-    q_ice = zero(ρ),
-    n_ice = zero(ρ),
-    q_rim = zero(ρ),
-    b_rim = zero(ρ),
-    logλ = zero(ρ),
+    ::Microphysics2Moment, mp::CMP.Microphysics2MParams{FT, WR, ICE}, tps,
+    ρ, T, q_tot, q_lcl, n_lcl, q_rai, n_rai,
+    q_ice = zero(ρ), n_ice = zero(ρ), q_rim = zero(ρ), b_rim = zero(ρ), logλ = zero(ρ),
 ) where {FT, WR, ICE <: CMP.P3IceParams}
     # Clamp negative inputs to zero (robustness against numerical errors)
     ρ = UT.clamp_to_nonneg(ρ)
@@ -824,7 +781,7 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch.
     db_rim_dt = zero(ρ)
 
     # --- Core Warm Rain Processes (shared helper) ---
-    warm = warm_rain_tendencies_2m(sb, q_lcl, q_rai, ρ, n_lcl, n_rai)
+    warm = warm_rain_tendencies_2m(sb, q_tot, q_lcl, q_rai, q_ice, ρ, n_lcl, n_rai)
     dq_lcl_dt = warm.dq_lcl_dt
     dn_lcl_dt = warm.dn_lcl_dt
     dq_rai_dt = warm.dq_rai_dt
@@ -834,26 +791,7 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch.
     N_lcl = ρ * n_lcl
     N_rai = ρ * n_rai
 
-    # --- Rain evaporation ---
-    evap = CM2.rain_evaporation(sb, aps, tps, zero(ρ), q_lcl, zero(ρ), q_rai, zero(ρ), ρ, N_rai, T)
-    dq_rai_dt += evap.evap_rate_1
-    dn_rai_dt += evap.evap_rate_0 / ρ
-
-    # --- Number adjustment for mass limits ---
-    # Cloud liquid
-    dn_lcl_inc = CM2.number_increase_for_mass_limit(sb.numadj, sb.pdf_c.xc_max, q_lcl, ρ, N_lcl)
-    dn_lcl_dec = CM2.number_decrease_for_mass_limit(sb.numadj, sb.pdf_c.xc_min, q_lcl, ρ, N_lcl)
-    dn_lcl_dt += (dn_lcl_inc + dn_lcl_dec) / ρ
-
-    # Rain
-    dn_rai_inc = CM2.number_increase_for_mass_limit(sb.numadj, sb.pdf_r.xr_max, q_rai, ρ, N_rai)
-    dn_rai_dec = CM2.number_decrease_for_mass_limit(sb.numadj, sb.pdf_r.xr_min, q_rai, ρ, N_rai)
-    dn_rai_dt += (dn_rai_inc + dn_rai_dec) / ρ
-
     # --- P3 Ice Processes ---
-    # NOTE: P3 uses gamma_inc_inv from SpecialFunctions which is NOT GPU-compatible
-    # (it uses string formatting for errors). We must keep if-branches to skip
-    # P3 code when there is no ice, otherwise GPU compilation fails.
     p3 = mp.ice.scheme
     vel = mp.ice.terminal_velocity
     pdf_c = mp.ice.cloud_pdf

src/Common.jl

@@ -463,6 +463,22 @@ See e.g., Seifert and Beheng (2006), https://doi.org/10.1007/s00703-005-0112-4,
 
 # Returns
 - `F_v(D)`: ventilation factor function (dimensionless)
+
+# Extended help
+
+The ventilation factor is given by
+
+```math
+F_v(D) = a_v + b_v ⋅ ∛Sc ⋅ √Re(D)
+```
+
+where
+
+```math
+Sc = ν_{air} / D_{vapor}
+Re(D) = D * v_{term}(D) / ν_{air}
+```
+
 """
 @inline function ventilation_factor(vent, aps, v_term)
     (; aᵥ, bᵥ) = vent