CliMA · haakon-e · Jul 2, 2026 · Jul 2, 2026 · Jul 2, 2026 · Jul 2, 2026
diff --git a/docs/bibliography.bib b/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}
+}
diff --git a/docs/make.jl b/docs/make.jl
@@ -72,6 +72,7 @@ pages = Any[
     "Home" => "index.md",
     "Parameterizations" => Parameterizations,
     "Bulk tendencies" => "BulkTendencies.md",
+    "Rosenbrock substepping" => "RosenbrockNumerics.md",
     "Thermodynamics interface" => "Thermodynamics.md",
     "How to guides" => Guides,
     "Models" => Models,

diff --git a/docs/src/API.md b/docs/src/API.md
@@ -133,9 +133,50 @@ BulkMicrophysicsTendencies.TendencyMode
 BulkMicrophysicsTendencies.Instantaneous
 BulkMicrophysicsTendencies.InstantaneousVerbose
 BulkMicrophysicsTendencies.LinearizedAverage
+BulkMicrophysicsTendencies.RosenbrockAverage
+BulkMicrophysicsTendencies.Verbose
+BulkMicrophysicsTendencies.Jacobian
+BulkMicrophysicsTendencies.DonorJacobian
+BulkMicrophysicsTendencies.CoupledDonorJacobian
+BulkMicrophysicsTendencies.ExactJacobian
+BulkMicrophysicsTendencies.GrowthTreatment
+BulkMicrophysicsTendencies.ImplicitGrowth
+BulkMicrophysicsTendencies.ExplicitGrowthDiagonal
+BulkMicrophysicsTendencies.TendencyLimiter
+BulkMicrophysicsTendencies.NoLimiter
+BulkMicrophysicsTendencies.EndStateSaturationAdjustment
+BulkMicrophysicsTendencies.rosenbrock_donor
+BulkMicrophysicsTendencies.rosenbrock_coupled
+BulkMicrophysicsTendencies.rosenbrock_exact
 BulkMicrophysicsTendencies.bulk_microphysics_tendencies
 ```
 
+Internal Rosenbrock substep helpers:
+
+```@docs
+BulkMicrophysicsTendencies.Instantaneous2MP3Tendency
+BulkMicrophysicsTendencies.ManualJacobian
+BulkMicrophysicsTendencies.MicroState1M
+BulkMicrophysicsTendencies.MicroState2MP3
+BulkMicrophysicsTendencies.Raw1MTendency
+BulkMicrophysicsTendencies._aggregate_tendencies
+BulkMicrophysicsTendencies._euler_update
+BulkMicrophysicsTendencies._full_species_mask
+BulkMicrophysicsTendencies._instantaneous_1m_tendency
+BulkMicrophysicsTendencies._instantaneous_2mp3_tendency
+BulkMicrophysicsTendencies._jacobian_2mp3_manual
+BulkMicrophysicsTendencies._linearize
+BulkMicrophysicsTendencies._per_process_1m
+BulkMicrophysicsTendencies._per_process_2mp3
+BulkMicrophysicsTendencies._rosenbrock_solve
+BulkMicrophysicsTendencies._rosenbrock_species_mask
+BulkMicrophysicsTendencies._rosenbrock_substep_verbose
+BulkMicrophysicsTendencies._rosenbrock_system
+BulkMicrophysicsTendencies._rosenbrock_update
+BulkMicrophysicsTendencies._species_mask
+BulkMicrophysicsTendencies.rosenbrock_manual
+```
+
 # P3 scheme
 
 ```@docs
@@ -221,6 +262,7 @@ These methods integrate over the particle size distribution.
 P3Scheme.D_m
 P3Scheme.ice_particle_terminal_velocity
 P3Scheme.ice_terminal_velocity_number_weighted
+P3Scheme.ice_terminal_velocity_number_weighted_from_prognostic
 P3Scheme.ice_terminal_velocity_mass_weighted
 ```
 
@@ -247,10 +289,19 @@ P3Scheme.bulk_liquid_ice_collision_sources
 Supporting methods:
 
 ```@docs
+P3Scheme.VolumetricCollisionRate
 P3Scheme.volumetric_collision_rate_integrand
 P3Scheme.compute_max_freeze_rate
 P3Scheme.compute_local_rime_density
 P3Scheme.get_liquid_integrals
+P3Scheme.crossing_integral_bounds
+P3Scheme.crossover_diameter
+P3Scheme.closed_rain_inner_NM
+P3Scheme.collision_cross_section_ice_liquid
+P3Scheme.collision_cross_section_ice_liquid_coeffs
+P3Scheme.collision_cross_section_ice_ice
+P3Scheme.ice_self_collection
+P3Scheme.wet_growth_onset_diameter
 P3Scheme.∫liquid_ice_collisions
 ```
 
@@ -262,8 +313,9 @@ P3Scheme.integrate
 P3Scheme.subintervals
 P3Scheme.ChebyshevGauss
 Quadrature.GaussLegendre
-Quadrature.build_quadrature
 P3Scheme.integral_bounds
+P3Scheme.velocity_integral_bounds
+P3Scheme.velocity_breakpoints
 ```
 
 # Aerosol model
@@ -406,6 +458,8 @@ Utilities.fac
 
 ```@docs
 Common
+Common.Chen2022VelocityCurve
+Common.particle_terminal_velocity
 Common.G_func_liquid
 Common.G_func_ice
 Common.logistic_function

diff --git a/docs/src/BulkTendencies.md b/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.
diff --git a/docs/src/DevelopersGuide.md b/docs/src/DevelopersGuide.md
@@ -206,3 +206,25 @@ Add a comment of what you are testing and use `@test` to create your test.
 The GPU tests are ran twice: for `Float64` and `Float32`.
 Similar as with performance tests, some trial and error is needed
   to find good tolerances for both options.
+
+### P3 quadrature error study
+
+The P3 scheme evaluates its size-distribution integrals with a Gauss-Legendre
+rule whose default order is recorded in `src/parameters/Microphysics2MParams.jl`.
+The study that supports the default lives in `test/p3_quadrature_error_study.jl`:
+it measures the relative error of every quadrature-consuming P3 quantity against
+a high-order reference, across a set of physically plausible column states, and
+attributes it to the transport components (sedimentation velocities, melt), the
+components proportional to the collision efficiency, and the full bulk tendency
+vector.
+
+Run it from the repository root with
+
+```
+julia --project=test -e 'include("test/p3_quadrature_error_study.jl");
+                         run_quadrature_error_study()'
+```
+
+Rerun the study when changing the quadrature rule, the integral bounds or
+breakpoints, or the P3 process integrands, and revisit the default order and the
+tolerances in `test/bulk_tendencies_quadrature_tests.jl` accordingly.