feat(Rosenbrock): add the 1M donor/coupled/exact modes on the framework

Extend the unified RosenbrockAverage framework to the one-moment model:
add the DonorCellJacobian and CoupledDonorJacobian providers and channel
masks, the rosenbrock_donor() and rosenbrock_coupled() presets, and the
Microphysics1Moment substep driver. LinearizedAverage now forwards to
rosenbrock_donor(), the donor-cell configuration of the framework. The
one-moment kernels are made ForwardDiff-able so rosenbrock_exact() runs
on the 1M model. Documentation and Float64/Float32 tests cover the
donor/coupled/exact presets, the keyword constructor, donor-equals-
LinearizedAverage, and zero-allocation hot calls.

Author: haakon-e

docs/bibliography.bib

@@ -819,4 +819,15 @@ @article{Horn2012
   pages={2635--2660},
   year={2011},
   publisher={G{\"o}ttingen, Germany}
-}
+}
+
+@article{Wan2020,
+  author = {Wan, Hui and Zhang, Shixuan and Rasch, Philip J. and Larson, Vincent E. and Zeng, Xubin and Yan, Huiping},
+  title = {Quantifying and attributing time step sensitivities in present-day climate simulations conducted with EAMv1},
+  journal = {Journal of Advances in Modeling Earth Systems},
+  year = {2020},
+  volume = {12},
+  number = {12},
+  doi = {10.1029/2019MS001982},
+  pages = {e2019MS001982}
+}

docs/src/API.md

@@ -146,6 +146,8 @@ BulkMicrophysicsTendencies.ExplicitGrowthDiagonal
 BulkMicrophysicsTendencies.TendencyLimiter
 BulkMicrophysicsTendencies.NoLimiter
 BulkMicrophysicsTendencies.EndStateSaturationAdjustment
+BulkMicrophysicsTendencies.rosenbrock_donor
+BulkMicrophysicsTendencies.rosenbrock_coupled
 BulkMicrophysicsTendencies.rosenbrock_exact
 BulkMicrophysicsTendencies.bulk_microphysics_tendencies
 ```

docs/src/BulkTendencies.md

@@ -153,3 +153,23 @@ This demonstrates that the linearized implicit substepping method provides a con
 
 - Average (implicit) bulk tendencies are currently implemented **only for the one-moment microphysics scheme**.
 - For other microphysics schemes, only **instantaneous bulk tendencies** are available at present.
+
+## Rosenbrock-averaged tendencies (2M+P3)
+
+For the 2-moment + P3 configuration, `RosenbrockAverage` replaces the hand-built linearization above with the exact Jacobian of the fused tendency, obtained by forward-mode automatic differentiation. The interval ``\Delta t`` is divided into `nsub` substeps of length ``h``, and each substep performs one linearized-implicit (Rosenbrock–Euler) update of the eight prognostic species ``x = (q_{\mathrm{lcl}}, n_{\mathrm{lcl}}, q_{\mathrm{rai}}, n_{\mathrm{rai}}, q_{\mathrm{ice}}, n_{\mathrm{ice}}, q_{\mathrm{rim}}, b_{\mathrm{rim}})``:
+
+```math
+\left(\frac{I}{h} - P J P\right) \Delta x = f(x), \qquad x \leftarrow \max(x + \Delta x,\, 0)
+```
+
+where ``f`` is the **raw instantaneous tendency** — the unmodified `Microphysics2Moment` process rates, with no timestep-dependent clipping — and ``J = \partial f / \partial x`` is its exact 8×8 `ForwardDiff` Jacobian.
+
+The differentiated tendency is the model physics, not a stabilized variant of it. The P3 condensation/deposition scheme is an analytic time-averaged relaxation [MorrisonMilbrandt2015](@cite) with no tendency clip; a supersaturation cap and ``1/h`` sink limits are explicit-Euler stabilization devices, not physical terms. They are unnecessary here because the one-stage Rosenbrock update is L-stable: it damps the stiff vapor-exchange subsystem monotonically, so the saturation overshoot and oscillation those limiters suppress cannot occur. They also degrade the solution: ``1/h`` tendency clips inject ``h``-independent error and break convergence under refinement [Wan2020](@cite), and a saturation clip structurally forbids the mixed-phase quasi-steady vapor pressure, which lies between liquid and ice saturation [KorolevMazin2003](@cite).
+
+The discrete stabilization steps are therefore conditioning and projection devices, all ``h``-free and applied to the linear solve rather than to the physics:
+
+- **Species projection** ``P = \mathrm{Diag}(z)``: species whose condensed mass is below ``10^{-10}`` are projected out of the Jacobian. Their rows of ``I/h - PJP`` reduce to the identity, so the solve returns exactly a forward-Euler update for those species while active species stay implicit — an IMEX-style splitting at species granularity. Near-empty species otherwise produce finite but very large Jacobian entries whose linearized steady state produces spurious number concentrations that substep refinement cannot remove.
+- **Equilibration** ``S = \mathrm{Diag}(|x| + h|f| + \epsilon)``: the linear system is solved as ``S^{-1} A S`` so the rows, which span roughly nine orders of magnitude across number and mass species, become O(1)-conditioned. This keeps single-precision roundoff relative to each species' own scale; an unscaled Float32 factorization deposits roundoff from the large rows into empty species as spurious mass.
+- **Positivity clamp** ``x \leftarrow \max(x + \Delta x, 0)``: a projection onto the physical nonnegative orthant after each substep.
+
+The local temperature is advanced each substep from the latent heating of the realized increments. `logλ` and `q_tot` are held fixed across the interval, matching the explicit-substepping semantics; non-finite states or Jacobians fall back to forward-Euler substeps of the raw tendency. The implicit update makes the stiff ice-process path insensitive to the substep length; at very large substeps (``h \gtrsim 100`` s) the single linearization carries the usual first-order error of a one-stage method, so increase `nsub` to refine.

docs/src/RosenbrockNumerics.md

@@ -17,6 +17,14 @@ substeps from the latent heat of the realized increment.
 `RosenbrockAverage` is parameterized by three independent option families:
 
 - **`Jacobian`** — the matrix `J` used in the substep solve.
+  - `DonorJacobian` — the donor-based linearization `M`: each transfer is linearized in its donor species,
+    vapor sources enter as a constant, and rates are floored by `max(q_min, q_donor)`. This is the matrix the
+    operational `LinearizedAverage` mode uses.
+  - `CoupledDonorJacobian` — the donor-based matrix with the vapor-competition (Wegener–Bergeron–Findeisen)
+    coupling added. The donor-based linearization keeps only donor-species slopes; restoring the dependence of
+    each rate on the shared vapor specific content recovers the cross-species coupling and corrects the sign of
+    the snow-from-cloud-liquid entry. The direct condensate dependence of the rates (rain ventilation, the
+    availability terms) is not recovered; use `ExactJacobian` for the full derivative.
   - `ExactJacobian` — the exact tendency derivative, formed with `ForwardDiff`.
 
 - **`GrowthTreatment`** — how the positive (growth) diagonal of `J` enters the implicit operator.
@@ -29,21 +37,25 @@ substeps from the latent heat of the realized increment.
   - `EndStateSaturationAdjustment` — scale the increment so the latent-heated end state does not cross ice
     saturation (see below).
 
-The supported preset configuration is:
+Three preset configurations are supported:
 
 | preset | Jacobian | growth | limiter |
 |---|---|---|---|
+| `rosenbrock_donor()` | `DonorJacobian` | `ImplicitGrowth` | `NoLimiter` |
+| `rosenbrock_coupled()` | `CoupledDonorJacobian` | `ImplicitGrowth` | `NoLimiter` |
 | `rosenbrock_exact()` | `ExactJacobian` | `ExplicitGrowthDiagonal` | `EndStateSaturationAdjustment` |
 
-On the two-moment + P3 model only `ExactJacobian` is available; use `rosenbrock_exact()`.
+`rosenbrock_donor()` reproduces `LinearizedAverage` (the operational donor-based scheme), now expressed within
+the unified framework; in `Float64` the two agree to round-off. On the two-moment + P3 model only
+`ExactJacobian` is available (there is no donor-based matrix there); use `rosenbrock_exact()`.
 
 ### Extending the framework
 
-To add a new Jacobian, define `struct MyJacobian <: Jacobian end` and a `_species_mask(::MyJacobian, ::GrowthTreatment)`
-method (returning a `x -> z` species projection). A new growth treatment is a `GrowthTreatment` subtype plus an
-`_apply_growth(::MyGrowth, J)` method; a new limiter is a `TendencyLimiter` subtype plus an
-`_apply_limiter(::MyLimiter, x, Δ, ...)` method. The substep driver dispatches on the option types at compile time,
-so a configured mode resolves with no run-time branch.
+To add a new Jacobian, define `struct MyJacobian <: Jacobian end` and the methods `_jacobian_provider(::MyJacobian)`
+(returning a `(g, x) -> J` provider) and `_species_mask(::MyJacobian, ::GrowthTreatment)`. A new growth treatment
+is a `GrowthTreatment` subtype plus an `_apply_growth(::MyGrowth, J)` method; a new limiter is a `TendencyLimiter`
+subtype plus an `_apply_limiter(::MyLimiter, x, Δ, ...)` method. The substep driver dispatches on the option types
+at compile time, so a configured mode resolves with no run-time branch.
 
 ## The coarse-step deposition instability
 
@@ -56,14 +68,14 @@ linear operator does not see): snow is over-deposited far past the available vap
 spurious temperature excursion, and the state goes non-physical. With the exact Jacobian this crash occurs at
 every coarse time step in the single-column convective test.
 
-### What cures it
+### What resolves it
 
 The exact preset removes the growth mode from the implicit operator and bounds the now-explicit growth by the
 physical saturation limit:
 
 - **`ExplicitGrowthDiagonal`** zeros the positive diagonal, so the implicit operator carries only non-positive
-  modes and is well-conditioned at any substep size. The exact off-diagonal couplings are retained, so accuracy
-  at cold, supersaturated cells is good.
+  modes and is well-conditioned at any substep size. The exact off-diagonal couplings (which the donor-based
+  matrix drops) are retained, so accuracy at cold, supersaturated cells is better than the donor scheme.
 - **`EndStateSaturationAdjustment`** scales the substep increment by the largest `s ∈ [0, 1]` for which the
   latent-heated end state stays at or above ice saturation. It acts only on cells that begin at or above
   saturation (a subsaturated, evaporating or sublimating cell cannot over-deposit, so its increment is returned
@@ -79,7 +91,7 @@ exact scheme robust across the resolved time-step envelope.
     precipitation at coarse time steps. Two or more substeps recover accurate precipitation; the saturation
     adjustment is then rarely active.
 
-## Approaches that did not cure the instability
+## Approaches that did not resolve the instability
 
 These were tried and are not part of the supported framework; they are recorded here because the failure modes
 are instructive.

docs/src/plots/BulkTendencies_plots.jl

@@ -115,7 +115,8 @@ function integrate_bulk_microphysics_linearized_one_step(
     Ls_over_cp = TDI.TD.Parameters.LH_s0(tps) / TDI.TD.Parameters.cp_d(tps)
 
     for i in 1:nsub
-        rates = BMT._linearized_implicit_step(
+        rates = BMT.bulk_microphysics_tendencies(
+            BMT.LinearizedAverage(),
             cm,
             mp,
             tps,