Re-add diagnostic-Nc activation for the 2M scheme

docs/src/API.md

@@ -67,7 +67,7 @@ Microphysics2M
 
 ```@autodocs
 Modules = [Parameters]
-Pages   = ["Microphysics2M.jl", "Microphysics2MParams.jl"]
+Pages   = ["Microphysics2M.jl", "Microphysics2MParams.jl", "ActivationSchemes.jl"]
 ```
 
 ## Size distributions
@@ -136,6 +136,7 @@ BulkMicrophysicsTendencies.Instantaneous
 BulkMicrophysicsTendencies.InstantaneousVerbose
 BulkMicrophysicsTendencies.LinearizedAverage
 BulkMicrophysicsTendencies.bulk_microphysics_tendencies
+BulkMicrophysicsTendencies.activation_source
 ```
 
 # P3 scheme

src/BulkMicrophysicsTendencies.jl

@@ -663,6 +663,38 @@ end
 
 # --- 2-Moment Microphysics Helper Functions ---
 
+"""
+    activation_source(scheme, tps, ρ, T, q_tot, q_lcl, q_ice, n_lcl, w, p)
+
+Return the cloud-droplet number source `dn_lcl_activation_dt` [kg⁻¹ air / s]
+implied by the activation `scheme` at the given local state.
+
+All arguments are passed by value so the function is pure and suitable for
+both broadcast and GPU kernels; inputs not needed by a particular scheme
+(e.g. `w`, `p` for [`CMP.DiagnosticNc`](@ref)) are simply ignored.
+
+# Arguments
+- `scheme`: activation scheme; dispatch target.
+- `tps`: thermodynamics parameters.
+- `ρ`, `T`, `q_tot`, `q_lcl`, `q_ice`: local thermodynamic state.
+- `n_lcl`: cloud-droplet number per mass of air [kg⁻¹].
+- `w`: vertical velocity [m/s] (for supersaturation-driven schemes).
+- `p`: pressure [Pa] (for supersaturation-driven schemes).
+"""
+@inline function activation_source(
+    ::CMP.NoActivation, tps, ρ, T, q_tot, q_lcl, q_ice, n_lcl, w, p,
+)
+    return zero(n_lcl)
+end
+
+@inline function activation_source(
+    scheme::CMP.DiagnosticNc, tps, ρ, T, q_tot, q_lcl, q_ice, n_lcl, w, p,
+)
+    FT = UT.promote_typeof(n_lcl, q_lcl)
+    target = ifelse(q_lcl > FT(scheme.q_thresh), FT(scheme.N_c), zero(FT))
+    return (target - n_lcl) / FT(scheme.τ_relax)
+end
+
 """
     warm_rain_tendencies_2m(sb, q_lcl, q_rai, ρ, n_lcl, n_rai)
 
@@ -708,7 +740,13 @@ Used by both warm-only and warm+ice dispatch methods to reduce code duplication.
     dn_rai_dt = zero(FT)
 
     # --- Aerosol activation ---
-    dn_lcl_activation_dt = zero(FT)
+    # `w`, `p` are the per-cell ambient inputs that supersaturation-driven
+    # schemes will need; positional so the call site can `@.`-broadcast.
+    # `NoActivation`/`DiagnosticNc` ignore them.
+    dn_lcl_activation_dt = activation_source(
+        warm_rain.activation_scheme, tps, ρ, T, q_tot, q_lcl, q_ice, n_lcl, w, p,
+    )
+    dn_lcl_dt += dn_lcl_activation_dt
 
     # --- Condensation of vapor / evaporation of cloud liquid water ---
     micro_mock = (; q_tot, q_lcl, q_icl = q_ice, q_rai, q_sno = zero(q_ice))

src/parameters/ActivationSchemes.jl

@@ -0,0 +1,63 @@
+export AbstractActivationScheme, NoActivation, DiagnosticNc
+
+"""
+    AbstractActivationScheme
+
+Abstract super-type for aerosol-activation closures used inside
+`BulkMicrophysicsTendencies`.
+
+Every concrete subtype is a pure, immutable parameter struct that fully
+specifies one activation closure. BMT dispatches on the scheme type to
+compute the cloud-droplet number source `dn_lcl_activation_dt` [kg⁻¹ s⁻¹]
+given the local thermodynamic state — the same options-pattern the 1-moment
+scheme uses for its prescribed-droplet-number autoconversion
+([`PrescribedNd`](@ref)).
+
+Currently provided: [`NoActivation`](@ref) (the null source, default) and
+[`DiagnosticNc`](@ref) (relaxation toward a prescribed droplet number).
+Supersaturation-driven closures (Twomey, Abdul-Razzak–Ghan) are deferred to
+a broader activation API.
+"""
+abstract type AbstractActivationScheme <: ParametersType end
+
+"""
+    NoActivation()
+
+Singleton scheme that suppresses aerosol activation: the activation source
+is zero regardless of state. The default — used when the host model supplies
+its own activation, and for tests that isolate other processes.
+"""
+struct NoActivation <: AbstractActivationScheme end
+
+"""
+    DiagnosticNc(; N_c, q_thresh = 1e-7, τ_relax = 60)
+
+Diagnostic-`N_c` activation: relax the droplet number toward a prescribed
+target. The RCEMIP-I default and the pre-ARG GCM standard.
+
+- While `q_lcl > q_thresh`, the droplet target is `N_c` [kg⁻¹ air] and the
+  tendency is `(N_c − n_lcl) / τ_relax`.
+- Otherwise the target is zero (cloud has evaporated) and the tendency is
+  `-n_lcl / τ_relax`.
+
+# Design notes
+- `τ_relax` is a FIXED physical timescale, NOT the integrator timestep:
+  relaxation by `1/dt` produces an infinitely-stiff source as adaptive
+  solvers refine the step. The 60 s default keeps `n_lcl` tracking cloud
+  mass within one macro-step while remaining integrable.
+- `q_thresh` avoids activating into essentially vapor-only cells, where
+  numerical noise in `q_lcl` would otherwise create phantom droplets.
+
+# Fields
+$(DocStringExtensions.FIELDS)
+"""
+@kwdef struct DiagnosticNc{FT} <: AbstractActivationScheme
+    "Droplet number target [kg⁻¹ air] while cloud mass is present"
+    N_c::FT
+    "Cloud-existence threshold on `q_lcl` [kg/kg]. Default: 1e-7."
+    q_thresh::FT = oftype(N_c, 1e-7)
+    "Fixed relaxation timescale [s]. Default: 60."
+    τ_relax::FT = oftype(N_c, 60)
+end
+
+Base.broadcastable(x::AbstractActivationScheme) = tuple(x)

src/parameters/Microphysics2MParams.jl

@@ -10,20 +10,27 @@ Parameters for 2-moment warm rain processes (Seifert-Beheng 2006).
 - `air_properties::AP`: AirProperties — air properties for evaporation
 - `condevap::CE`: MM2015 cond-evap relaxation timescale
 - `subdep::SD`: MM2015 sub-dep relaxation timescale
+- `activation_scheme::AS`: aerosol → cloud-droplet activation closure
+  (defaults to [`NoActivation`](@ref); set to e.g. [`DiagnosticNc`](@ref)
+  to enable activation in the fused 2M tendency).
 """
-@kwdef struct WarmRainParams2M{SB, AP, CE, SD} <: ParametersType
+@kwdef struct WarmRainParams2M{SB, AP, CE, SD, AS <: AbstractActivationScheme} <: ParametersType
     seifert_beheng::SB
     air_properties::AP
     condevap::CE
     subdep::SD
+    activation_scheme::AS = NoActivation()
 end
 # Construct WarmRainParams2M from a ClimaParams TOML dictionary
-WarmRainParams2M(toml_dict::CP.ParamDict; is_limited = true) =
+WarmRainParams2M(toml_dict::CP.ParamDict; is_limited = true,
+    activation_scheme::AbstractActivationScheme = NoActivation(),
+) =
     WarmRainParams2M(;
         seifert_beheng = SB2006(toml_dict; is_limited),
         air_properties = AirProperties(toml_dict),
         condevap = CondEvap2M(toml_dict),
         subdep = SubDep2M(toml_dict),
+        activation_scheme,
     )
 
 Base.show(io::IO, mime::MIME"text/plain", x::WarmRainParams2M) =
@@ -146,10 +153,11 @@ Create a `Microphysics2MParams` object from a ClimaParams TOML dictionary.
 Microphysics2MParams(toml_dict::CP.ParamDict;
     with_ice = false, is_limited = true,
     quadrature = QUAD.ChebyshevGauss(100),
+    activation_scheme = NoActivation(),
     inp_depletion_model = NIceProxyDepletion(τ_act = 300),
 ) = Microphysics2MParams(;
     # Warm rain parameters (always present)
-    warm_rain = WarmRainParams2M(toml_dict; is_limited),
+    warm_rain = WarmRainParams2M(toml_dict; is_limited, activation_scheme),
     # Optional ice phase parameters
     ice = with_ice ?
           P3IceParams(toml_dict; is_limited, quadrature, inp_depletion_model) :

src/parameters/Parameters.jl

@@ -49,6 +49,7 @@ include("MicrophysicsP3.jl")
 include("TerminalVelocity.jl")
 
 # Unified parameter containers
+include("ActivationSchemes.jl")
 include("Microphysics0MParams.jl")
 include("Microphysics1MOptions.jl")
 include("Microphysics1MParams.jl")

test/activation_schemes_tests.jl

@@ -0,0 +1,74 @@
+using Test
+
+import ClimaParams as CP
+import CloudMicrophysics as CM
+import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.BulkMicrophysicsTendencies as BMT
+import CloudMicrophysics.ThermodynamicsInterface as TDI
+import ForwardDiff as FD
+
+function test_activation_schemes(FT)
+    tps = TDI.TD.Parameters.ThermodynamicsParameters(FT)
+    args = (tps, FT(1.05), FT(288), FT(0.012), FT(2e-3), FT(0), FT(5e7), FT(0), FT(0))
+    splice(a; q_lcl = a[5], n_lcl = a[7]) = (a[1:4]..., q_lcl, a[6], n_lcl, a[8:9]...)
+
+    @testset "NoActivation is the null source ($FT)" begin
+        @test BMT.activation_source(CMP.NoActivation(), args...) === zero(FT)
+    end
+
+    @testset "DiagnosticNc relaxation ($FT)" begin
+        scheme = CMP.DiagnosticNc(; N_c = FT(1e8))
+        # cloudy: relax up toward the target
+        s = BMT.activation_source(scheme, args...)
+        @test s ≈ (FT(1e8) - FT(5e7)) / FT(60)
+        # above target: relax down
+        s_dn = BMT.activation_source(scheme, splice(args; n_lcl = FT(2e8))...)
+        @test s_dn ≈ (FT(1e8) - FT(2e8)) / FT(60)
+        # cloud-free: target zero (decay of leftover droplets)
+        s0 = BMT.activation_source(scheme, splice(args; q_lcl = FT(0))...)
+        @test s0 ≈ -FT(5e7) / FT(60)
+        # custom timescale and threshold
+        s2 = BMT.activation_source(
+            CMP.DiagnosticNc(; N_c = FT(1e8), q_thresh = FT(1e-2), τ_relax = FT(30)), args...,
+        )
+        @test s2 ≈ -FT(5e7) / FT(30)  # q_lcl below the raised threshold
+    end
+
+    @testset "activation in the fused 2M tendency ($FT)" begin
+        mp_on = CMP.Microphysics2MParams(
+            FT; activation_scheme = CMP.DiagnosticNc(; N_c = FT(1e8)),
+        )
+        mp_off = CMP.Microphysics2MParams(FT)
+        @test mp_off.warm_rain.activation_scheme isa CMP.NoActivation
+        t_on = BMT.bulk_microphysics_tendencies(
+            BMT.Microphysics2Moment(), mp_on, tps,
+            FT(1.05), FT(288), FT(0.012), FT(2e-3), FT(5e7), FT(1e-4), FT(4e4),
+        )
+        t_off = BMT.bulk_microphysics_tendencies(
+            BMT.Microphysics2Moment(), mp_off, tps,
+            FT(1.05), FT(288), FT(0.012), FT(2e-3), FT(5e7), FT(1e-4), FT(4e4),
+        )
+        @test t_off.dn_lcl_activation_dt === zero(FT)
+        @test t_on.dn_lcl_activation_dt ≈ (FT(1e8) - FT(5e7)) / FT(60)
+        # the activation source is folded into the total droplet tendency
+        @test t_on.dn_lcl_dt - t_off.dn_lcl_dt ≈ t_on.dn_lcl_activation_dt
+    end
+
+    @testset "activation source is concretely typed under mixed args ($FT)" begin
+        D = FD.Dual{Nothing, FT, 8}
+        scheme = CMP.DiagnosticNc(; N_c = FT(1e8))
+        for sig in (
+            (FT, FT, FT, D, FT, D, FT, FT),  # Dual q_lcl, n_lcl
+            (FT, FT, FT, FT, FT, D, FT, FT), # Dual n_lcl only
+            (FT, FT, FT, D, FT, FT, FT, FT), # Dual q_lcl only
+        )
+            rts = Base.return_types(
+                BMT.activation_source, Tuple{typeof(scheme), typeof(tps), FT, sig...},
+            )
+            @test all(isconcretetype, rts)
+        end
+    end
+end
+
+test_activation_schemes(Float64)
+test_activation_schemes(Float32)

test/runtests.jl

@@ -27,6 +27,7 @@ TT.@testset "All tests" begin
     include("ventilation_tests.jl")
     include("ad_compat_tests.jl")
     include("return_type_tests.jl")
+    include("activation_schemes_tests.jl")
     include("DistributionTools_tests.jl")
     include("unrolled_logsumexp.jl")
 end