diff --git a/src/cache/precipitation_precomputed_quantities.jl b/src/cache/precipitation_precomputed_quantities.jl
index ca172846d6..3c67f5f246 100644
--- a/src/cache/precipitation_precomputed_quantities.jl
+++ b/src/cache/precipitation_precomputed_quantities.jl
@@ -4,6 +4,7 @@
 
 import CloudMicrophysics.MicrophysicsNonEq as CMNe
 import CloudMicrophysics.Microphysics1M as CM1
+import CloudMicrophysics.Microphysics2M as CM2
 
 import Thermodynamics as TD
 import ClimaCore.Operators as Operators
@@ -97,6 +98,94 @@ function set_precipitation_velocities!(
         ) / Y.c.ρ
     return nothing
 end
+function set_precipitation_velocities!(
+    Y,
+    p,
+    moisture_model::NonEquilMoistModel,
+    precip_model::Microphysics2Moment,
+)
+    (; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwNₗ, ᶜwNᵢ, ᶜwNᵣ, ᶜwNₛ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜts, ᶜu) = p.precomputed
+    (; q_liq, q_ice, q_rai, q_sno) = p.precomputed.ᶜspecific
+    (; ᶜΦ) = p.core
+
+    cm1c = CAP.microphysics_cloud_params(p.params)
+    cm1p = CAP.microphysics_1m_params(p.params)
+    cm2p = CAP.microphysics_2m_params(p.params)
+    thp = CAP.thermodynamics_params(p.params)
+
+    # compute the precipitation terminal velocity [m/s]
+    # TODO sedimentation of snow is based on the 1M scheme
+    @. ᶜwNᵣ = getindex(CM2.rain_terminal_velocity(
+        cm2p.sb,
+        cm2p.tv,
+        q_rai,
+        Y.c.ρ,
+        Y.c.N_rai,
+    ), 1)
+    @. ᶜwᵣ = getindex(CM2.rain_terminal_velocity(
+        cm2p.sb,
+        cm2p.tv,
+        q_rai,
+        Y.c.ρ,
+        Y.c.N_rai,
+    ), 2)
+    @. ᶜwNₛ = CM1.terminal_velocity(
+        cm1p.ps,
+        cm1p.tv.snow,
+        Y.c.ρ,
+        q_sno,
+    )
+    @. ᶜwₛ = CM1.terminal_velocity(
+        cm1p.ps,
+        cm1p.tv.snow,
+        Y.c.ρ,
+        q_sno,
+    )
+    # compute sedimentation velocity for cloud condensate [m/s]
+    # TODO sedimentation velocities of cloud condensates are based 
+    # on the 1M scheme.
+    @. ᶜwNₗ = CMNe.terminal_velocity(
+        cm1c.liquid,
+        cm1c.Ch2022.rain,
+        Y.c.ρ,
+        q_liq,
+    )
+    @. ᶜwₗ = CMNe.terminal_velocity(
+        cm1c.liquid,
+        cm1c.Ch2022.rain,
+        Y.c.ρ,
+        q_liq,
+    )
+    @. ᶜwNᵢ = CMNe.terminal_velocity(
+        cm1c.ice,
+        cm1c.Ch2022.small_ice,
+        Y.c.ρ,
+        q_ice,
+    )
+    @. ᶜwᵢ = CMNe.terminal_velocity(
+        cm1c.ice,
+        cm1c.Ch2022.small_ice,
+        Y.c.ρ,
+        q_ice,
+    )
+
+    # compute their contributions to energy and total water advection
+    @. ᶜwₜqₜ =
+        Geometry.WVector(
+            ᶜwₗ * Y.c.ρq_liq +
+            ᶜwᵢ * Y.c.ρq_ice +
+            ᶜwᵣ * Y.c.ρq_rai +
+            ᶜwₛ * Y.c.ρq_sno,
+        ) / Y.c.ρ
+    @. ᶜwₕhₜ =
+        Geometry.WVector(
+            ᶜwₗ * Y.c.ρq_liq * (Iₗ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwₗ, ᶜu))) +
+            ᶜwᵢ * Y.c.ρq_ice * (Iᵢ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwᵢ, ᶜu))) +
+            ᶜwᵣ * Y.c.ρq_rai * (Iₗ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwᵣ, ᶜu))) +
+            ᶜwₛ * Y.c.ρq_sno * (Iᵢ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwₛ, ᶜu))),
+        ) / Y.c.ρ
+    return nothing
+end
 
 """
     set_precipitation_cache!(Y, p, precip_model, turbconv_model)
@@ -262,6 +351,64 @@ function set_precipitation_cache!(
     # in edmf sub-domains.
     return nothing
 end
+function set_precipitation_cache!(Y, p, ::Microphysics2Moment, _)
+    (; dt) = p
+    (; ᶜts) = p.precomputed
+    (; ᶜSqₗᵖ, ᶜSqᵢᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ) = p.precomputed
+    (; ᶜSNₗᵖ, ᶜSNᵢᵖ, ᶜSNᵣᵖ, ᶜSNₛᵖ) = p.precomputed
+
+    (; q_liq, q_rai, q_sno) = p.precomputed.ᶜspecific
+
+    ᶜSᵖ = p.scratch.ᶜtemp_scalar
+    ᶜS₂ᵖ = p.scratch.ᶜtemp_scalar_2
+
+    # get thermodynamics and microphysics params
+    (; params) = p
+    cm1p = CAP.microphysics_1m_params(params)
+    cm2p = CAP.microphysics_2m_params(params)
+    thp = CAP.thermodynamics_params(params)
+
+    # compute warm precipitation sources on the grid mean (based on SB2006 2M scheme)
+    compute_warm_precipitation_sources_2M!(
+        ᶜSᵖ,
+        ᶜS₂ᵖ,
+        ᶜSNₗᵖ,
+        ᶜSNᵣᵖ,
+        ᶜSqₗᵖ,
+        ᶜSqᵣᵖ,
+        Y.c.ρ,
+        Y.c.N_liq,
+        Y.c.N_rai,
+        q_liq,
+        q_rai,
+        ᶜts,
+        dt,
+        cm2p,
+        thp,
+    )
+
+    #TODO - implement 2M cold processes!
+    
+    return nothing
+end
+function set_precipitation_cache!(
+    Y,
+    p,
+    ::Microphysics2Moment,
+    ::DiagnosticEDMFX,
+)
+    error("Not implemented yet")
+    return nothing
+end
+function set_precipitation_cache!(
+    Y,
+    p,
+    ::Microphysics2Moment,
+    ::PrognosticEDMFX,
+)
+    error("Not implemented yet")
+    return nothing
+end
 
 """
     set_precipitation_surface_fluxes!(Y, p, precipitation model)
@@ -301,7 +448,7 @@ end
 function set_precipitation_surface_fluxes!(
     Y,
     p,
-    precip_model::Microphysics1Moment,
+    precip_model::Union{Microphysics1Moment, Microphysics2Moment},
 )
     (; surface_rain_flux, surface_snow_flux) = p.precomputed
     (; col_integrated_precip_energy_tendency,) = p.conservation_check
diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
index f713f3fa34..9d8bef8b8c 100644
--- a/src/cache/precomputed_quantities.jl
+++ b/src/cache/precomputed_quantities.jl
@@ -126,20 +126,35 @@ function precomputed_quantities(Y, atmos)
     )
     sedimentation_quantities =
         atmos.moisture_model isa NonEquilMoistModel ?
-        (; ᶜwₗ = similar(Y.c, FT), ᶜwᵢ = similar(Y.c, FT)) : (;)
+        (; ᶜwₗ = similar(Y.c, FT), ᶜwᵢ = similar(Y.c, FT), ᶜwNₗ = similar(Y.c, FT), ᶜwNᵢ = similar(Y.c, FT)) : (;)
     if atmos.precip_model isa Microphysics0Moment
         precipitation_quantities =
             (; ᶜS_ρq_tot = similar(Y.c, FT), ᶜS_ρe_tot = similar(Y.c, FT))
-    elseif atmos.precip_model isa Microphysics1Moment
-        precipitation_quantities = (;
-            ᶜwᵣ = similar(Y.c, FT),
-            ᶜwₛ = similar(Y.c, FT),
-            ᶜSqₗᵖ = similar(Y.c, FT),
-            ᶜSqᵢᵖ = similar(Y.c, FT),
-            ᶜSqᵣᵖ = similar(Y.c, FT),
-            ᶜSqₛᵖ = similar(Y.c, FT),
-        )
-    else
+        elseif atmos.precip_model isa Microphysics1Moment
+            precipitation_quantities = (;
+                ᶜwᵣ = similar(Y.c, FT),
+                ᶜwₛ = similar(Y.c, FT),
+                ᶜSqₗᵖ = similar(Y.c, FT),
+                ᶜSqᵢᵖ = similar(Y.c, FT),
+                ᶜSqᵣᵖ = similar(Y.c, FT),
+                ᶜSqₛᵖ = similar(Y.c, FT),
+            )
+        elseif atmos.precip_model isa Microphysics2Moment
+            precipitation_quantities = (;
+                ᶜwᵣ = similar(Y.c, FT),
+                ᶜwₛ = similar(Y.c, FT),
+                ᶜSqₗᵖ = similar(Y.c, FT),
+                ᶜSqᵢᵖ = similar(Y.c, FT),
+                ᶜSqᵣᵖ = similar(Y.c, FT),
+                ᶜSqₛᵖ = similar(Y.c, FT),
+                ᶜwNᵣ = similar(Y.c, FT),
+                ᶜwNₛ = similar(Y.c, FT),
+                ᶜSNₗᵖ = similar(Y.c, FT),
+                ᶜSNᵢᵖ = similar(Y.c, FT),
+                ᶜSNᵣᵖ = similar(Y.c, FT),
+                ᶜSNₛᵖ = similar(Y.c, FT),
+            )
+        else
         precipitation_quantities = (;)
     end
     precipitation_sgs_quantities =
diff --git a/src/diagnostics/Diagnostics.jl b/src/diagnostics/Diagnostics.jl
index 8191e0fad4..66fb063405 100644
--- a/src/diagnostics/Diagnostics.jl
+++ b/src/diagnostics/Diagnostics.jl
@@ -28,6 +28,7 @@ import ..NonEquilMoistModel
 import ..NoPrecipitation
 import ..Microphysics0Moment
 import ..Microphysics1Moment
+import ..Microphysics2Moment
 
 # radiation
 import ClimaAtmos.RRTMGPInterface as RRTMGPI
diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl
index ed70a1e7a6..846ac20e40 100644
--- a/src/diagnostics/core_diagnostics.jl
+++ b/src/diagnostics/core_diagnostics.jl
@@ -707,6 +707,7 @@ function compute_pr!(
         NoPrecipitation,
         Microphysics0Moment,
         Microphysics1Moment,
+        Microphysics2Moment,
     },
 )
     if isnothing(out)
@@ -742,6 +743,7 @@ function compute_prra!(
         NoPrecipitation,
         Microphysics0Moment,
         Microphysics1Moment,
+        Microphysics2Moment,
     },
 )
     if isnothing(out)
@@ -774,6 +776,7 @@ function compute_prsn!(
         NoPrecipitation,
         Microphysics0Moment,
         Microphysics1Moment,
+        Microphysics2Moment,
     },
 )
     if isnothing(out)
@@ -805,7 +808,10 @@ function compute_husra!(
     state,
     cache,
     time,
-    precip_model::Microphysics1Moment,
+    precip_model::Union{
+        Microphysics1Moment,
+        Microphysics2Moment,
+    },
 )
     if isnothing(out)
         return state.c.ρq_rai ./ state.c.ρ
@@ -836,7 +842,10 @@ function compute_hussn!(
     state,
     cache,
     time,
-    precip_model::Microphysics1Moment,
+    precip_model::Union{
+        Microphysics1Moment,
+        Microphysics2Moment,
+    },
 )
     if isnothing(out)
         return state.c.ρq_sno ./ state.c.ρ
diff --git a/src/initial_conditions/InitialConditions.jl b/src/initial_conditions/InitialConditions.jl
index eb172507f3..6086d3b4a3 100644
--- a/src/initial_conditions/InitialConditions.jl
+++ b/src/initial_conditions/InitialConditions.jl
@@ -7,6 +7,7 @@ import ..NonEquilMoistModel
 import ..NoPrecipitation
 import ..Microphysics0Moment
 import ..Microphysics1Moment
+import ..Microphysics2Moment
 import ..PrescribedSurfaceTemperature
 import ..PrognosticSurfaceTemperature
 import ..ᶜinterp
diff --git a/src/initial_conditions/atmos_state.jl b/src/initial_conditions/atmos_state.jl
index b24ac4ebc2..ac3d99c282 100644
--- a/src/initial_conditions/atmos_state.jl
+++ b/src/initial_conditions/atmos_state.jl
@@ -127,6 +127,32 @@ precip_variables(ls, ::Microphysics1Moment) = (;
     ρq_rai = ls.ρ * ls.precip_state.q_rai,
     ρq_sno = ls.ρ * ls.precip_state.q_sno,
 )
+precip_variables(ls, ::Microphysics2Moment) = (;
+    ρq_rai = ls.ρ * ls.precip_state.q_rai,
+    ρq_sno = ls.ρ * ls.precip_state.q_sno,
+    N_rai = ls.precip_state.N_rai,
+    N_sno = ls.precip_state.N_sno,
+    N_liq = ls.precip_state.N_liq,
+    N_ice = ls.precip_state.N_ice,
+)
+
+"""
+    precip_variables(ls, ::MicrophysicsP3)
+
+Define prognostic variables for P3 microphysics scheme.
+"""
+precip_variables(ls, ::MicrophysicsP3) = (;
+    # Warm rain variables (from 2M)
+    ρq_rai = ls.ρ * ls.precip_state.q_rai,
+    N_rai = ls.precip_state.N_rai,
+    N_liq = ls.precip_state.N_liq,
+    # P3 ice variables
+    ρq_ice = ls.ρ * ls.precip_state.q_ice,
+    N_ice = ls.precip_state.N_ice,  # Rimed mass mixing ratio
+    L_ice = ls.precip_state.L_ice,  # Ice mass mixing ratio
+    B_rim = ls.precip_state.B_rim,  # Rimed volume mixing ratio
+    L_rim = ls.precip_state.L_rim,  # Rimed mass mixing ratio
+)
 
 # We can use paper-based cases for LES type configurations (no TKE)
 # or SGS type configurations (initial TKE needed), so we do not need to assert
diff --git a/src/parameterized_tendencies/microphysics/cloud_condensate.jl b/src/parameterized_tendencies/microphysics/cloud_condensate.jl
index 7275c10b2b..55c9d9736a 100644
--- a/src/parameterized_tendencies/microphysics/cloud_condensate.jl
+++ b/src/parameterized_tendencies/microphysics/cloud_condensate.jl
@@ -16,7 +16,7 @@ function cloud_condensate_tendency!(
     ::Union{NoPrecipitation, Microphysics0Moment},
 )
     error(
-        "NonEquilMoistModel can only be run with Microphysics1Moment precipitation",
+        "NonEquilMoistModel can only be run with Microphysics1Moment or Microphysics2Moment precipitation",
     )
 end
 
@@ -37,3 +37,24 @@ function cloud_condensate_tendency!(
     @. Yₜ.c.ρq_liq += Y.c.ρ * cloud_sources(cmc.liquid, thp, ᶜts, q_rai, dt)
     @. Yₜ.c.ρq_ice += Y.c.ρ * cloud_sources(cmc.ice, thp, ᶜts, q_sno, dt)
 end
+
+function cloud_condensate_tendency!(
+    Yₜ,
+    Y,
+    p,
+    ::NonEquilMoistModel,
+    ::Microphysics2Moment,
+)
+    (; ᶜts) = p.precomputed
+    (; params, dt) = p
+    (; q_rai, q_sno) = p.precomputed.ᶜspecific
+    FT = eltype(params)
+    thp = CAP.thermodynamics_params(params)
+    cmc = CAP.microphysics_cloud_params(params)
+
+    @. Yₜ.c.ρq_liq += Y.c.ρ * cloud_sources(cmc.liquid, thp, ᶜts, q_rai, dt)
+    @. Yₜ.c.ρq_ice += Y.c.ρ * cloud_sources(cmc.ice, thp, ᶜts, q_sno, dt)
+
+    @. Yₜ.c.N_liq += aerosol_activation_sources(cmc.liquid, thp, ᶜts, Y.c.ρ, q_rai, Y.c.N_ice, cmc.N_cloud_liquid_droplets, dt)
+    @. Yₜ.c.N_ice += aerosol_activation_sources(cmc.ice, thp, ᶜts, Y.c.ρ, q_sno, Y.c.N_ice, cmc.N_cloud_liquid_droplets, dt)
+end
diff --git a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl
index 598b71cbda..5b6fff8fc6 100644
--- a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl
+++ b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl
@@ -4,6 +4,7 @@ import Thermodynamics as TD
 import CloudMicrophysics.Microphysics0M as CM0
 import CloudMicrophysics.Microphysics1M as CM1
 import CloudMicrophysics.Microphysics2M as CM2
+import CloudMicrophysics.MicrophysicsP3 as P3
 import CloudMicrophysics.MicrophysicsNonEq as CMNe
 import CloudMicrophysics.Parameters as CMP
 
@@ -62,11 +63,12 @@ function ml_N_cloud_liquid_droplets(cmc, c_dust, c_seasalt, c_SO4, q_liq)
 end
 
 """
-    cloud_sources(cm_params, thp, ts, dt)
+    cloud_mass_sources(cm_params, thp, ts, dt)
 
  - cm_params - CloudMicrophysics parameters struct for cloud water or ice condensate
  - thp - Thermodynamics parameters struct
  - ts - thermodynamics state
+ - qᵣ or qₛ - rain or snow specific humidity 
  - dt - model time step
 
 Returns the condensation/evaporation or deposition/sublimation rate for
@@ -154,7 +156,7 @@ end
  - thp, cmp - structs with thermodynamic and microphysics parameters
 
 Returns the q source terms due to precipitation formation from the 1-moment scheme.
-The specific humidity source terms are defined as defined as Δmᵢ / (m_dry + m_tot)
+The specific humidity source terms are defined as Δmᵢ / (m_dry + m_tot)
 where i stands for total, rain or snow.
 Also returns the total energy source term due to the microphysics processes.
 """
@@ -278,7 +280,7 @@ end
  - thp, cmp - structs with thermodynamic and microphysics parameters
 
 Returns the q source terms due to precipitation sinks from the 1-moment scheme.
-The specific humidity source terms are defined as defined as Δmᵢ / (m_dry + m_tot)
+The specific humidity source terms are defined as Δmᵢ / (m_dry + m_tot)
 where i stands for total, rain or snow.
 Also returns the total energy source term due to the microphysics processes.
 """
@@ -324,3 +326,201 @@ function compute_precipitation_sinks!(
     @. Sqₛᵖ += Sᵖ
     #! format: on
 end
+
+#####
+##### 2M microphysics
+#####
+
+"""
+    aerosol_activation_sources(cm_params, thp, ts, dt)
+
+ - cm_params - CloudMicrophysics parameters struct for cloud water or ice condensate
+ - thp - Thermodynamics parameters struct
+ - ts - thermodynamics state
+ - ρ - air density
+ - qₚ - precipitation (rain or snow) specific humidity
+ - N_dp - droplets (liquid or ice) number concentration
+ _ N_dp_prescribed - droplets (liquid or ice) prescribed number concentration
+ - dt - model time step
+
+Returns the activation rate. #TODO This function temporarily computes activation rate 
+based on mass rates and a prescribed droplet mass (no activation parameterization yet).
+"""
+function aerosol_activation_sources(
+    cm_params::Union{CMP.CloudLiquid{FT}, CMP.CloudIce{FT}},
+    thp,
+    ts,
+    ρ,
+    qₚ,
+    N_dp,
+    N_dp_prescribed,
+    dt,
+) where {FT}
+    r_dp = FT(2e-6) # 2 μm
+    m_dp = 4 / 3 * π * r_dp^3
+    S = ρ * cloud_sources(cm_params, thp, ts, qₚ, dt) / m_dp
+
+    return ifelse(
+        S > FT(0),
+        triangle_inequality_limiter(S, limit((N_dp_prescribed - N_dp), dt, 2)),
+        -triangle_inequality_limiter(abs(S), limit(N_dp, dt, 2)),
+        )
+end
+
+"""
+    compute_warm_precipitation_sources_2M!(Sᵖ, S₂ᵖ, SNₗᵖ, SNᵣᵖ, Sqₗᵖ, Sqᵣᵖ, ρ, Nₗ, Nᵣ, qₗ, qᵣ, dt, sb, thp)
+
+ - Sᵖ, S₂ᵖ - temporary containters to help compute precipitation source terms
+ - SNₗᵖ, SNᵣᵖ, Sqₗᵖ, Sqᵣᵖ - cached storage for precipitation source terms
+ - ρ - air density
+ - Nₗ, Nᵣ - cloud liquid and rain number concentration
+ - qₗ, qᵣ - cloud liquid and rain specific humidity
+ - ts - thermodynamic state (see td package for details)
+ - dt - model time step
+ - thp, mp - structs with thermodynamic and microphysics parameters
+
+Computes precipitation number and mass sources due to warm precipitation processes based on the 2-moment 
+Seifert and Beheng (2006) scheme.
+"""
+function compute_warm_precipitation_sources_2M!(
+    Sᵖ,
+    S₂ᵖ,
+    SNₗᵖ,
+    SNᵣᵖ,
+    Sqₗᵖ,
+    Sqᵣᵖ,
+    ρ,
+    Nₗ,
+    Nᵣ,
+    qₗ,
+    qᵣ,
+    ts,
+    dt,
+    mp,
+    thp,
+)
+
+    FT = eltype(thp)
+    @. SNₗᵖ = ρ * FT(0)
+    @. SNᵣᵖ = ρ * FT(0)
+    @. Sqₗᵖ = ρ * FT(0)
+    @. Sqᵣᵖ = ρ * FT(0)
+
+    # auto-conversion (mass)
+    @. Sᵖ = triangle_inequality_limiter(
+        CM2.autoconversion(mp.sb.acnv, mp.sb.pdf_c, qₗ, qᵣ, ρ, Nₗ).dq_rai_dt,
+        limit(qₗ, dt, 5),
+    )
+    @. Sqₗᵖ -= Sᵖ
+    @. Sqᵣᵖ += Sᵖ
+
+    # auto-conversion (number) and liquid self-collection
+    @. Sᵖ = triangle_inequality_limiter(
+        CM2.autoconversion(mp.sb.acnv, mp.sb.pdf_c, qₗ, qᵣ, ρ, Nₗ).dN_liq_dt,
+        limit(Nₗ, dt, 10),
+    )
+    @. S₂ᵖ = -triangle_inequality_limiter(
+        -CM2.liquid_self_collection(mp.sb.acnv, mp.sb.pdf_c, qₗ, ρ, Sᵖ),
+        limit(Nₗ, dt, 5)
+    )
+    @. SNₗᵖ += Sᵖ
+    @. SNₗᵖ += S₂ᵖ
+    @. SNᵣᵖ -= 0.5*Sᵖ
+
+    # rain self-collection and breakup
+    @. Sᵖ = -triangle_inequality_limiter(
+        -CM2.rain_self_collection(mp.sb.pdf_r, mp.sb.self, qᵣ, ρ, Nᵣ),
+        limit(Nᵣ, dt, 5),
+    )
+    @. S₂ᵖ = triangle_inequality_limiter(
+        CM2.rain_breakup(mp.sb.pdf_r, mp.sb.brek, qᵣ, ρ, Nᵣ, Sᵖ),
+        limit(Nᵣ, dt, 5),
+    )
+    @. SNᵣᵖ += Sᵖ
+    @. SNᵣᵖ += S₂ᵖ
+
+    # accretion (mass)
+    @. Sᵖ = triangle_inequality_limiter(
+        CM2.accretion(mp.sb, qₗ, qᵣ, ρ, Nₗ).dq_rai_dt,
+        limit(qₗ, dt, 5),
+    )
+    @. Sqₗᵖ -= Sᵖ
+    @. Sqᵣᵖ += Sᵖ
+
+    # accretion (number)
+    @. Sᵖ = -triangle_inequality_limiter(
+        -CM2.accretion(mp.sb, qₗ, qᵣ, ρ, Nₗ).dN_liq_dt,
+        limit(Nₗ, dt, 5),
+    )
+    @. SNₗᵖ += Sᵖ
+
+    # evaporation (mass)
+    @. Sᵖ = -triangle_inequality_limiter(
+        -CM2.rain_evaporation(mp.sb, mp.aps, thp, PP(thp, ts), qᵣ, ρ, Nᵣ, Tₐ(thp, ts)).evap_rate_1,
+        limit(qᵣ, dt, 5),
+    )
+    @. Sqᵣᵖ += Sᵖ
+
+    # evaporation (number)
+    @. Sᵖ = -triangle_inequality_limiter(
+        -CM2.rain_evaporation(mp.sb, mp.aps, thp, PP(thp, ts), qᵣ, ρ, Nᵣ, Tₐ(thp, ts)).evap_rate_0,
+        limit(Nᵣ, dt, 5),
+    )
+    @. SNᵣᵖ += Sᵖ
+
+end
+
+"""
+    compute_p3_microphysics_sources!(Sᵖ,S₂ᵖ,SNₗᵖ,SNᵣᵖ,Sqₗᵖ,Sqᵣᵖ,Sqᵢᵖ,SL_rimᵖ,SB_rimᵖ,SL_iceᵖ,ρ,Nₗ,Nᵣ,qₗ,qᵣ,qᵢ,L_rim,B_rim,L_ice,ts,dt,mp,thp,)
+
+Compute source terms for P3 microphysics scheme, combining 2M warm rain with P3 ice microphysics.
+
+# Arguments
+- `Sᵖ, S₂ᵖ`: Temporary containers to help compute precipitation source terms
+- `SNᵢᵖ, SNᵣᵖ`: Cached storage for precipitation source terms
+- `ρ`: Air density
+- `Nₗ, Nᵣ`: Cloud liquid and rain number concentration
+- `qₗ, qᵣ`: Cloud liquid and rain specific humidity
+- `L_ice`: Ice mass
+- `N_ice`: Ice number concentration
+- `L_rim`: Rimed mass
+- `B_rim`: Rimed volume
+- `ts`: Thermodynamic state
+- `dt`: Model time step
+- `mp`: Microphysics parameters
+- `thp`: Thermodynamic parameters
+
+Returns the source terms for the P3 microphysics scheme.
+"""
+function compute_p3_microphysics_sources!(
+    Sᵖ, S₂ᵖ,
+    SNᵢᵖ, SNᵣᵖ, SL_rimᵖ, SB_rimᵖ, SL_iceᵖ,
+    ρ, Nₗ, Nᵣ, qₗ, qᵣ,
+    L_ice, N_ice,
+    L_rim, B_rim,
+    ts, dt, mp, thp,
+)
+    FT = eltype(thp)
+
+    # Initialize all tendencies to zero
+    @. SNᵢᵖ = ρ * FT(0)
+    @. SNᵣᵖ = ρ * FT(0)
+    @. SL_rimᵖ = ρ * FT(0)
+    @. SB_rimᵖ = ρ * FT(0)
+    @. SL_iceᵖ = ρ * FT(0)
+
+    # TODO: Add P3 ice microphysics processes:
+    # - Ice nucleation
+    # @. Sᵖ = P3.het_ice_nucleation(...)
+    # - Ice growth/sublimation
+    # - Riming
+    # - Melting
+    # @. Sᵖ = P3.ice_melt(...)
+    # - Ice-ice collection
+    # - Ice-rain collection
+    # @. Sᵖ = P3.bulk_collision_rate_with_liquid(...)
+    # - Ice sedimentation
+    # - Ice breakup
+
+    return nothing
+end
\ No newline at end of file
diff --git a/src/parameterized_tendencies/microphysics/precipitation.jl b/src/parameterized_tendencies/microphysics/precipitation.jl
index 84079fe244..dd547e9ff3 100644
--- a/src/parameterized_tendencies/microphysics/precipitation.jl
+++ b/src/parameterized_tendencies/microphysics/precipitation.jl
@@ -165,3 +165,115 @@ function precipitation_tendency!(
         @. Yₜ.c.ρq_sno += Y.c.sgsʲs.:($$j).ρa * ᶜSqₛᵖʲs.:($$j)
     end
 end
+
+#####
+##### 2-moment microphysics without sgs scheme
+#####
+
+function precipitation_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    ::DryModel,
+    precip_model::Microphysics2Moment,
+    _,
+)
+    error("Microphysics2Moment precipitation should not be used with DryModel")
+end
+function precipitation_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    ::EquilMoistModel,
+    precip_model::Microphysics2Moment,
+    _,
+)
+    error(
+        "Microphysics2Moment precipitation and EquilMoistModel precipitation_tendency is not implemented",
+    )
+end
+function precipitation_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    ::NonEquilMoistModel,
+    precip_model::Microphysics2Moment,
+    _,
+)
+    (; ᶜSqₗᵖ, ᶜSqᵢᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ) = p.precomputed
+    (; ᶜSNₗᵖ, ᶜSNᵢᵖ, ᶜSNᵣᵖ, ᶜSNₛᵖ) = p.precomputed
+
+    # Update grid mean tendencies
+    @. Yₜ.c.ρq_liq += Y.c.ρ * ᶜSqₗᵖ
+    @. Yₜ.c.ρq_ice += Y.c.ρ * ᶜSqᵢᵖ
+    @. Yₜ.c.ρq_rai += Y.c.ρ * ᶜSqᵣᵖ
+    @. Yₜ.c.ρq_sno += Y.c.ρ * ᶜSqₛᵖ
+
+    @. Yₜ.c.N_liq += ᶜSNₗᵖ
+    @. Yₜ.c.N_ice += ᶜSNᵢᵖ
+    @. Yₜ.c.N_rai += ᶜSNᵣᵖ
+    @. Yₜ.c.N_sno += ᶜSNₛᵖ
+
+    return nothing
+end
+function precipitation_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    ::NonEquilMoistModel,
+    precip_model::Microphysics2Moment,
+    turbconv_model::DiagnosticEDMFX,
+)
+    error("Not implemented yet")
+end
+function precipitation_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    ::NonEquilMoistModel,
+    precip_model::Microphysics1Moment,
+    turbconv_model::PrognosticEDMFX,
+)
+    error("Not implemented yet")
+end
+
+"""
+    precipitation_tendency!(Yₜ, Y, p, t, ::NonEquilMoistModel, ::MicrophysicsP3, _)
+
+Compute tendencies for P3 microphysics scheme, combining 2M warm rain with P3 ice microphysics.
+"""
+function precipitation_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    ::NonEquilMoistModel,
+    precip_model::MicrophysicsP3,
+    _,
+)
+    # Warm rain tendencies (from 2M)
+    (; ᶜSqₗᵖ, ᶜSqᵣᵖ) = p.precomputed
+    (; ᶜSNₗᵖ, ᶜSNᵣᵖ) = p.precomputed
+    
+    # P3 ice tendencies
+    (; ᶜSqᵢᵖ, ᶜSL_rimᵖ, ᶜSB_rimᵖ, ᶜSL_iceᵖ) = p.precomputed
+
+    # Update warm rain tendencies
+    @. Yₜ.c.ρq_liq += Y.c.ρ * ᶜSqₗᵖ
+    @. Yₜ.c.ρq_rai += Y.c.ρ * ᶜSqᵣᵖ
+    @. Yₜ.c.N_liq += ᶜSNₗᵖ
+    @. Yₜ.c.N_rai += ᶜSNᵣᵖ
+
+    # Update P3 ice tendencies
+    @. Yₜ.c.ρq_ice += Y.c.ρ * ᶜSqᵢᵖ
+    @. Yₜ.c.L_rim += ᶜSL_rimᵖ
+    @. Yₜ.c.B_rim += ᶜSB_rimᵖ
+    @. Yₜ.c.L_ice += ᶜSL_iceᵖ
+
+    return nothing
+end
diff --git a/src/parameters/Parameters.jl b/src/parameters/Parameters.jl
index 16931aa087..3423955abd 100644
--- a/src/parameters/Parameters.jl
+++ b/src/parameters/Parameters.jl
@@ -68,6 +68,8 @@ Base.@kwdef struct ClimaAtmosParameters{
     MPC,
     MP0M,
     MP1M,
+    MP2M,
+    MPP3,
     SFP,
     TCP,
     STP,
@@ -80,6 +82,8 @@ Base.@kwdef struct ClimaAtmosParameters{
     microphysics_cloud_params::MPC
     microphysics_0m_params::MP0M
     microphysics_1m_params::MP1M
+    microphysics_2m_params::MP2M
+    microphysics_p3_params::MPP3
     surface_fluxes_params::SFP
     turbconv_params::TCP
     surface_temp_params::STP
diff --git a/src/parameters/create_parameters.jl b/src/parameters/create_parameters.jl
index 22e816cb1f..deaffdb6cd 100644
--- a/src/parameters/create_parameters.jl
+++ b/src/parameters/create_parameters.jl
@@ -50,8 +50,11 @@ function ClimaAtmosParameters(toml_dict::TD) where {TD <: CP.AbstractTOMLDict}
 
     microphysics_0m_params = CM.Parameters.Parameters0M(toml_dict)
     microphysics_1m_params = microphys_1m_parameters(toml_dict)
+    microphysics_2m_params = microphys_2m_parameters(toml_dict)
+    microphysics_p3_params = microphys_p3_parameters(toml_dict)
     MP0M = typeof(microphysics_0m_params)
     MP1M = typeof(microphysics_1m_params)
+    MP2M = typeof(microphysics_2m_params)
 
     vert_diff_params = vert_diff_parameters(toml_dict)
     VDP = typeof(vert_diff_params)
@@ -69,6 +72,8 @@ function ClimaAtmosParameters(toml_dict::TD) where {TD <: CP.AbstractTOMLDict}
         MPC,
         MP0M,
         MP1M,
+        MP2M,
+        MPP3,
         SFP,
         TCP,
         STP,
@@ -82,6 +87,8 @@ function ClimaAtmosParameters(toml_dict::TD) where {TD <: CP.AbstractTOMLDict}
         microphysics_cloud_params,
         microphysics_0m_params,
         microphysics_1m_params,
+        microphysics_2m_params,
+        microphysics_p3_params,
         surface_fluxes_params,
         turbconv_params,
         surface_temp_params,
@@ -152,6 +159,24 @@ microphys_1m_parameters(toml_dict::CP.AbstractTOMLDict) = (;
     ).prescribed_cloud_droplet_number_concentration,
 )
 
+microphys_2m_parameters(::Type{FT}) where {FT <: AbstractFloat} =
+    microphys_2m_parameters(CP.create_toml_dict(FT))
+
+microphys_2m_parameters(toml_dict::CP.AbstractTOMLDict) = (;
+    sb = CM.Parameters.SB2006(toml_dict),
+    aps = CM.Parameters.AirProperties(toml_dict),
+    tv = CM.Parameters.Chen2022VelType(toml_dict),
+)
+
+"""
+    microphys_p3_parameters(toml_dict)
+
+Create parameters for P3 microphysics scheme. This combines 2M warm rain with P3 ice microphysics.
+"""
+microphys_p3_parameters(toml_dict::CP.AbstractTOMLDict) = (;
+    p3 = CM.Parameters.ParametersP3(toml_dict),
+)
+
 function vert_diff_parameters(toml_dict)
     name_map = (; :C_E => :C_E, :H_diffusion => :H, :D_0_diffusion => :D₀)
     return CP.get_parameter_values(toml_dict, name_map, "ClimaAtmos")
diff --git a/src/prognostic_equations/implicit/manual_sparse_jacobian.jl b/src/prognostic_equations/implicit/manual_sparse_jacobian.jl
index fb6db02ae7..a338713d4a 100644
--- a/src/prognostic_equations/implicit/manual_sparse_jacobian.jl
+++ b/src/prognostic_equations/implicit/manual_sparse_jacobian.jl
@@ -151,7 +151,7 @@ function jacobian_cache(alg::ManualSparseJacobian, Y, atmos)
     sfc_if_available = is_in_Y(@name(sfc)) ? (@name(sfc),) : ()
 
     condensate_names =
-        (@name(c.ρq_liq), @name(c.ρq_ice), @name(c.ρq_rai), @name(c.ρq_sno))
+        (@name(c.ρq_liq), @name(c.ρq_ice), @name(c.ρq_rai), @name(c.ρq_sno), @name(c.N_liq), @name(c.N_ice), @name(c.N_rai), @name(c.N_sno))
     available_condensate_names =
         MatrixFields.unrolled_filter(is_in_Y, condensate_names)
     available_tracer_names =
@@ -592,6 +592,18 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
                 DiagonalMatrixRow(-Geometry.WVector(ᶜwₚ) / ᶜρ) - (I,)
         end
 
+        # ∂ᶜN_liq_err_∂ᶜN_liq = matrix[@name(c.N_liq), @name(c.N_liq)]
+        # @. ∂ᶜN_liq_err_∂ᶜN_liq = zero(typeof(∂ᶜN_liq_err_∂ᶜN_liq)) - (I,)
+
+        # ∂ᶜN_ice_err_∂ᶜN_ice = matrix[@name(c.N_ice), @name(c.N_ice)]
+        # @. ∂ᶜN_ice_err_∂ᶜN_ice = zero(typeof(∂ᶜN_ice_err_∂ᶜN_ice)) - (I,)
+
+        # ∂ᶜN_rai_err_∂ᶜN_rai = matrix[@name(c.N_rai), @name(c.N_rai)]
+        # @. ∂ᶜN_rai_err_∂ᶜN_rai = zero(typeof(∂ᶜN_rai_err_∂ᶜN_rai)) - (I,)
+
+        # ∂ᶜN_sno_err_∂ᶜN_sno = matrix[@name(c.N_sno), @name(c.N_sno)]
+        # @. ∂ᶜN_sno_err_∂ᶜN_sno = zero(typeof(∂ᶜN_sno_err_∂ᶜN_sno)) - (I,)
+
     end
 
     if use_derivative(diffusion_flag)
diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl
index fe069a3376..ff14a96340 100644
--- a/src/solver/model_getters.jl
+++ b/src/solver/model_getters.jl
@@ -305,6 +305,8 @@ function get_precipitation_model(parsed_args)
         Microphysics0Moment()
     elseif precip_model == "1M"
         Microphysics1Moment()
+    elseif precip_model == "2M"
+        Microphysics2Moment()
     else
         error("Invalid precip_model $(precip_model)")
     end
diff --git a/src/solver/types.jl b/src/solver/types.jl
index 804fd8a3c9..bb7451e064 100644
--- a/src/solver/types.jl
+++ b/src/solver/types.jl
@@ -16,6 +16,8 @@ abstract type AbstractPrecipitationModel end
 struct NoPrecipitation <: AbstractPrecipitationModel end
 struct Microphysics0Moment <: AbstractPrecipitationModel end
 struct Microphysics1Moment <: AbstractPrecipitationModel end
+struct Microphysics2Moment <: AbstractPrecipitationModel end
+struct MicrophysicsP3 <: AbstractPrecipitationModel end
 
 """