P3 melting + F_liq

docs/src/P3Scheme.md

@@ -342,23 +342,36 @@ Melting rate is derived in the same way as in the
   [1-moment scheme](https://clima.github.io/CloudMicrophysics.jl/dev/Microphysics1M/#Snow-melt).
 We assume the same ventilation factor parameterization as in [SeifertBeheng2006](@cite),
   and use the terminal velocity parameterization from [Chen2022](@cite).
-The ``dm/dD`` derivative is computed for each P3 size regime.
-The bulk melting rate is computed by numerically integrating over the particle
-  size distribution:
+The ``dm/dD`` derivative is computed for each P3 size regime, without taking into account
+  ``L_{liq}``, since melting is dependent only on the ice core and not on the entire
+  mixed-phase particle.
+The bulk melting rates are computed by piecewise and numerically integrating over the ice core
+  size distribution using the terminal velocity of the mixed-phase particle. Particles with ``D < D_{th}``
+  are transferred directly to rain in one time step, whereas melting on particles with ``D > D_{th}``
+  is a source of ``L_{liq}``.
 ```math
 \begin{equation}
-  \left. \frac{dL}{dt} \right|_{melt} = \frac{4 \, K_{thermo}}{L_f} \left(T - T_{freeze}\right)
-    \int_{0}^{\infty} \frac{dm(D)}{dD} \frac{F_v(D) N(D)}{D}
+  \left. \frac{dL_{rai}}{dt} \right|_{melt} = \frac{4 \, K_{thermo}}{L_f} \left(T - T_{freeze}\right)
+    \int_{0}^{D_{th}} \frac{dm(D)}{dD} \frac{F_v(D) N(D)}{D}
 \end{equation}
 ```
-The melting rate for number concentration is assumed to be proportional to the
-  ice content melting rate.
 ```math
 \begin{equation}
-  \left. \frac{dN}{dt} \right|_{melt} = \frac{N}{L} \left. \frac{dL}{dt} \right|_{melt}
+  \left. \frac{dL_{liq}}{dt} \right|_{melt} = \frac{4 \, K_{thermo}}{L_f} \left(T - T_{freeze}\right)
+    \int_{D_{th}}^{\infty} \frac{dm(D)}{dD} \frac{F_v(D) N(D)}{D}
 \end{equation}
 ```
-Both rates are limited by the total available ice content and number concentration
+Sinks of ``L_{p3, tot}``, ``L_{rim}``, and ``B_{rim}`` are calculated accordingly.
+The melting rate for number concentration is assumed to be proportional to
+  the change in ``L_{ice}`` from ``\frac{dL_{rai}}{dt}`` since this represents
+  an immediate transfer to rain rather than an accumulation of ``L_{liq}``.
+```math
+\begin{equation}
+  \left. \frac{dN}{dt} \right|_{melt} = \frac{N}{L_{ice}} \left. \frac{dL_{rai}}{dt} \right|_{melt}
+\end{equation}
+```
+If ``F_{liq} > 0.99``, ``L_{p3, tot}`` is completely transferred to rain.
+All rates are limited by the total available mass contents and number concentration
   divided by model time step length.
 
 ```@example

docs/src/plots/P3Melting.jl

@@ -20,6 +20,7 @@ dt = FT(1)
 # initial ice content and number concentration
 Lᵢ = FT(1e-4)
 Nᵢ = FT(2e5)
+F_liq = FT(0)
 
 # tested temperature range
 ΔT_range = range(1e-4, stop = 0.025, length = 1000)
@@ -33,49 +34,71 @@ max_dNdt = [Nᵢ / dt for ΔT in ΔT_range]
 ρₐ1 = FT(1.2)
 Fᵣ1 = FT(0.8)
 ρᵣ1 = FT(800)
-dLdt1 = [P3.ice_melt(pps, vel.snow_ice, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ1, ρᵣ1, dt).dLdt for ΔT in ΔT_range]
-dNdt1 = [P3.ice_melt(pps, vel.snow_ice, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ1, ρᵣ1, dt).dNdt for ΔT in ΔT_range]
+dLdt1_rai = [P3.ice_melt(pps, vel, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ1, ρᵣ1, F_liq, dt).dLdt_rai for ΔT in ΔT_range]
+dLdt1_liq = [P3.ice_melt(pps, vel, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ1, ρᵣ1, F_liq, dt).dLdt_liq for ΔT in ΔT_range]
+dNdt1 = [P3.ice_melt(pps, vel, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ1, ρᵣ1, F_liq, dt).dNdt_ice for ΔT in ΔT_range]
 
 Fᵣ2 = FT(0.2)
-dLdt2 = [P3.ice_melt(pps, vel.snow_ice, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ2, ρᵣ1, dt).dLdt for ΔT in ΔT_range]
-dNdt2 = [P3.ice_melt(pps, vel.snow_ice, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ2, ρᵣ1, dt).dNdt for ΔT in ΔT_range]
+dLdt2_rai = [P3.ice_melt(pps, vel, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ2, ρᵣ1, F_liq, dt).dLdt_rai for ΔT in ΔT_range]
+dLdt2_liq = [P3.ice_melt(pps, vel, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ2, ρᵣ1, F_liq, dt).dLdt_liq for ΔT in ΔT_range]
+dNdt2 = [P3.ice_melt(pps, vel, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ2, ρᵣ1, F_liq, dt).dNdt_ice for ΔT in ΔT_range]
 
 ρᵣ2 = FT(200)
-dLdt3 = [P3.ice_melt(pps, vel.snow_ice, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ2, ρᵣ2, dt).dLdt for ΔT in ΔT_range]
-dNdt3 = [P3.ice_melt(pps, vel.snow_ice, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ2, ρᵣ2, dt).dNdt for ΔT in ΔT_range]
+dLdt3_rai = [P3.ice_melt(pps, vel, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ2, ρᵣ2, F_liq, dt).dLdt_rai for ΔT in ΔT_range]
+dLdt3_liq = [P3.ice_melt(pps, vel, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ2, ρᵣ2, F_liq, dt).dLdt_liq for ΔT in ΔT_range]
+dNdt3 = [P3.ice_melt(pps, vel, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ2, ρᵣ2, F_liq, dt).dNdt_ice for ΔT in ΔT_range]
 
 ρₐ2 = FT(0.5)
-dLdt4 = [P3.ice_melt(pps, vel.snow_ice, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ2, Fᵣ2, ρᵣ2, dt).dLdt for ΔT in ΔT_range]
-dNdt4 = [P3.ice_melt(pps, vel.snow_ice, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ2, Fᵣ2, ρᵣ2, dt).dNdt for ΔT in ΔT_range]
+dLdt4_rai = [P3.ice_melt(pps, vel, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ2, Fᵣ2, ρᵣ2, F_liq, dt).dLdt_rai for ΔT in ΔT_range]
+dLdt4_liq = [P3.ice_melt(pps, vel, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ2, Fᵣ2, ρᵣ2, F_liq, dt).dLdt_liq for ΔT in ΔT_range]
+dNdt4 = [P3.ice_melt(pps, vel, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ2, Fᵣ2, ρᵣ2, F_liq, dt).dNdt_ice for ΔT in ΔT_range]
 
 # plotting
-fig = PL.Figure(size = (1500, 500), fontsize=22, linewidth=3)
+fig = PL.Figure(size = (2000, 1500), fontsize=24, linewidth=3)
 
 ax1 = PL.Axis(fig[1, 1]; yscale = log10)
 ax2 = PL.Axis(fig[1, 2]; yscale = log10)
+ax3 = PL.Axis(fig[2, 1]; yscale = log10)
+ax4 = PL.Axis(fig[2, 2]; yscale = log10)
 
 ax1.xlabel = "T [C]"
-ax1.ylabel = "ice mass melting rate [g/m3/s]"
+ax1.ylabel = "total ice mass melting rate [g/m3/s]"
 ax2.xlabel = "T [C]"
 ax2.ylabel = "ice number melting rate [1/cm3/s]"
-
-l_max_dLdt = PL.lines!(ax1, ΔT_range,  max_dLdt * 1e3,  color = :thistle)
+ax3.xlabel = "T [C]"
+ax3.ylabel = "rain mass source [g/m3/s]"
+ax4.xlabel = "T [C]"
+ax4.ylabel = "liquid mass content source [g/m3/s]"
+
+for ax in [ax1, ax3, ax4]
+    l_max_dLdt = PL.lines!(ax, ΔT_range,  max_dLdt * 1e3,  color = :thistle)
+end
 l_max_dNdt = PL.lines!(ax2, ΔT_range,  max_dNdt * 1e-6, color = :thistle)
 
-l_dLdt1 = PL.lines!(ax1, ΔT_range,  dLdt1 * 1e3,  color = :skyblue)
+l_dLdt1 = PL.lines!(ax1, ΔT_range,  (dLdt1_rai + dLdt1_liq) * 1e3,  color = :skyblue)
 l_dNdt1 = PL.lines!(ax2, ΔT_range,  dNdt1 * 1e-6, color = :skyblue)
+l_dLdt1_rai = PL.lines!(ax3, ΔT_range,  (dLdt1_rai) * 1e3,  color = :skyblue)
+l_dLdt1_liq = PL.lines!(ax4, ΔT_range,  (dLdt1_liq) * 1e3,  color = :skyblue)
+
 
-l_dLdt2 = PL.lines!(ax1, ΔT_range,  dLdt2 * 1e3,  color = :blue3)
+l_dLdt2 = PL.lines!(ax1, ΔT_range,  (dLdt2_rai + dLdt2_liq) * 1e3,  color = :blue3)
 l_dNdt2 = PL.lines!(ax2, ΔT_range,  dNdt2 * 1e-6, color = :blue3)
+l_dLdt2_rai = PL.lines!(ax3, ΔT_range,  (dLdt2_rai) * 1e3,  color = :blue3)
+l_dLdt2_liq = PL.lines!(ax4, ΔT_range,  (dLdt2_liq) * 1e3,  color = :blue3)
 
-l_dLdt3 = PL.lines!(ax1, ΔT_range,  dLdt3 * 1e3,  color = :orchid)
+l_dLdt3 = PL.lines!(ax1, ΔT_range,  (dLdt3_rai + dLdt3_liq) * 1e3,  color = :orchid)
 l_dNdt3 = PL.lines!(ax2, ΔT_range,  dNdt3 * 1e-6, color = :orchid)
+l_dLdt3_rai = PL.lines!(ax3, ΔT_range,  (dLdt3_rai) * 1e3,  color = :orchid)
+l_dLdt3_liq = PL.lines!(ax4, ΔT_range,  (dLdt3_liq) * 1e3,  color = :orchid)
 
-l_dLdt4 = PL.lines!(ax1, ΔT_range,  dLdt4 * 1e3,  color = :purple)
+l_dLdt4 = PL.lines!(ax1, ΔT_range,  (dLdt4_rai + dLdt4_liq) * 1e3,  color = :purple)
 l_dNdt4 = PL.lines!(ax2, ΔT_range,  dNdt4 * 1e-6, color = :purple)
+l_dLdt4_rai = PL.lines!(ax3, ΔT_range,  (dLdt4_rai) * 1e3,  color = :purple)
+l_dLdt4_liq = PL.lines!(ax4, ΔT_range,  (dLdt4_liq) * 1e3,  color = :purple)
+
 
 PL.Legend(
-    fig[1, 3],
+    fig[1:2, 3],
     [l_max_dNdt, l_dNdt1, l_dNdt2, l_dNdt3, l_dNdt4],
     [
        "limit",
@@ -86,4 +109,5 @@ PL.Legend(
     ],
     framevisible = false,
 )
+PL.resize_to_layout!(fig)
 PL.save("P3_ice_melt.svg", fig)

docs/src/plots/P3SchemePlots.jl

@@ -51,11 +51,11 @@ function p3_relations_plot()
     # define plot axis
     #[row, column]
     ax1 = define_axis(fig, 1:7,  1:9,   CMK.L"m(D) regime for $ρ_r = 400 kg m^{-3}$", CMK.L"$m$ (kg)", [1e-10, 1e-8, 1e-6], 1.9)
-    ax3 = define_axis(fig, 1:7,  10:18, CMK.L"m(D) regime for $F_rim = 0.95$",          CMK.L"$m$ (kg)", [1e-10, 1e-8, 1e-6], 1.8)
+    ax3 = define_axis(fig, 1:7,  10:18, CMK.L"m(D) regime for $F_{rim} = 0.95$",          CMK.L"$m$ (kg)", [1e-10, 1e-8, 1e-6], 1.8)
     ax2 = define_axis(fig, 8:15, 1:9,   CMK.L"A(D) regime for $ρ_r = 400 kg m^{-3}$", CMK.L"$A$ ($m^2$)", [1e-8, 1e-6, 1e-4], 1.7)
-    ax4 = define_axis(fig, 8:15, 10:18, CMK.L"A(D) regime for $F_rim = 0.95$", CMK.L"$A$ ($m^2$)", [1e-8, 1e-6, 1e-4], 1.6)
+    ax4 = define_axis(fig, 8:15, 10:18, CMK.L"A(D) regime for $F_{rim} = 0.95$", CMK.L"$A$ ($m^2$)", [1e-8, 1e-6, 1e-4], 1.6)
     ax5 = define_axis(fig, 16:22, 1:9,   CMK.L"ρ(D) regime for $ρ_r = 400 kg m^{-3}$", CMK.L"$ρ$ ($kg m^{-3}$)", [100, 500, 900], 2, logscale = false)
-    ax6 = define_axis(fig, 16:22, 10:18, CMK.L"ρ(D) regime for $F_rim = 0.95$", CMK.L"$ρ$ ($kg m^{-3}$)", [100, 500, 900], 1.9, logscale = false)
+    ax6 = define_axis(fig, 16:22, 10:18, CMK.L"ρ(D) regime for $F_{rim} = 0.95$", CMK.L"$ρ$ ($kg m^{-3}$)", [100, 500, 900], 1.9, logscale = false)
 
     # Get thresholds
     sol4_0 = P3.thresholds(p3, 400.0, 0.0)

src/P3_particle_properties.jl

@@ -102,10 +102,10 @@ function thresholds(p3::PSP3{FT}, ρ_r::FT, F_rim::FT) where {FT}
     @assert F_rim >= FT(0)   # rime mass fraction must be positive ...
     @assert F_rim < FT(1)    # ... and there must always be some unrimed part
 
-    if F_rim == FT(0)
+    if F_rim == FT(0) || ρ_r == FT(0)
         return (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0))
     else
-        @assert ρ_r > FT(0)   # rime density must be positive ...
+        @assert ρ_r >= FT(0)   # rime density must be positive ...
         @assert ρ_r <= p3.ρ_l # ... and as a bulk ice density can't exceed the density of water
 
         P3_problem(ρ_d) =

src/P3_processes.jl

@@ -1,3 +1,24 @@
+"""
+A structure containing the rates of change of the specific humidities and number
+densities of liquid and rain water.
+"""
+Base.@kwdef struct P3Rates{FT}
+    "Rate of change of total ice mass content"
+    dLdt_p3_tot::FT = FT(0)
+    "Rate of change of rime mass content"
+    dLdt_rim::FT = FT(0)
+    "Rate of change of mass content of liquid on ice"
+    dLdt_liq::FT = FT(0)
+    "Rate of change of rime volume"
+    ddtB_rim::FT = FT(0)
+    "Rate of change of ice number concentration"
+    dNdt_ice::FT = FT(0)
+    "Rate of change of rain mass content"
+    dLdt_rai::FT = FT(0)
+    "Rate of change of rain number concentration"
+    dNdt_rai::FT = FT(0)
+end
+
 """
     het_ice_nucleation(pdf_c, p3, tps, q, N, T, ρₐ, p, aerosol)
 
@@ -53,64 +74,148 @@ end
  - Chen2022 - struct containing Chen 2022 velocity parameters
  - aps - air properties
  - tps - thermodynamics parameters
- - L_ice - ice content
+ - L_p3_tot- total ice content
  - N_ice - ice number concentration
  - T - temperature (K)
  - ρ_a - air density
  - F_r - rime mass fraction (q_rim/ q_i)
  - ρ_r - rime density (q_rim/B_rim)
  - dt - model time step (for limiting the tendnecy)
 
-Returns the melting rate of ice (QIMLT in Morrison and Mildbrandt (2015)).
+Returns the melting rate of ice (QIMLT in Morrison and Mildbrandt (2015))
+modified for the liquid fraction scheme described in Cholette et al (2019).
 """
 function ice_melt(
     p3::PSP3,
-    Chen2022::CMP.Chen2022VelTypeSnowIce,
+    Chen2022::CMP.Chen2022VelType,
     aps::CMP.AirProperties{FT},
     tps::TDP.ThermodynamicsParameters{FT},
-    L_ice::FT,
+    L_p3_tot::FT,
     N_ice::FT,
     Tₐ::FT,
     ρₐ::FT,
     F_rim::FT,
     ρ_rim::FT,
+    F_liq::FT,
     dt::FT,
 ) where {FT}
-    # process not dependent on F_liq
-    # (we want ice core shape params)
-    F_liq_ = FT(0)
-    # Get constants
-    (; ν_air, D_vapor, K_therm) = aps
-    L_f = TD.latent_heat_fusion(tps, Tₐ)
-    N_sc = ν_air / D_vapor
-
-    # Get the P3 diameter distribution...
-    th = thresholds(p3, ρ_rim, F_rim)
-    (λ, N_0) =
-        distribution_parameter_solver(p3, L_ice, N_ice, ρ_rim, F_rim, F_liq_)
-    N(D) = N′ice(p3, D, λ, N_0)
-    # ... and D_max for the integral
-    bound = get_ice_bound(p3, λ, FT(1e-6))
-
-    # Ice particle terminal velocity
-    v(D) = ice_particle_terminal_velocity(p3, D, Chen2022, ρₐ, F_rim, th)
-    # Reynolds number
-    N_Re(D) = D * v(D) / ν_air
-    # Ventillation factor
-    F_v(D) = p3.vent_a + p3.vent_b * N_sc^(1 / 3) * N_Re(D)^(1 / 2)
-    dmdD(D) = p3_dmdD(p3, D, F_rim, th)
-
-    f(D) = 4 * K_therm / L_f * (Tₐ - p3.T_freeze) * dmdD(D) / D * F_v(D) * N(D)
-    # Integrate
-    (dLdt, error) = QGK.quadgk(d -> f(d), 0, bound, rtol = FT(1e-6))
-
-    # only consider melting (not fusion)
-    dLdt = max(0, dLdt)
-    # compute change of N_ice proportional to change in L
-    dNdt = N_ice / L_ice * dLdt
+    # initialize rates and compute solid ice content
+    (;
+        dLdt_p3_tot,
+        dLdt_rim,
+        dLdt_liq,
+        dLdt_rai,
+        dNdt_ice,
+        dNdt_rai,
+        ddtB_rim,
+    ) = P3Rates{FT}()
+    L_ice = (1 - F_liq) * L_p3_tot
+    if F_liq > 0.99
+        # transfer total ice mass to rain
+        dLdt_p3_tot = L_p3_tot
+        dNdt_ice = N_ice
+        dLdt_rim = F_rim * L_ice
+        ddtB_rim = dLdt_rim / ρ_rim
+        dLdt_rai = L_p3_tot
+        dNdt_rai = N_ice
+    elseif L_ice > eps(FT) && N_ice > eps(FT)
+        # process not dependent on F_liq
+        # (we want ice core shape params)
+        # define temp F_liq_ = 0 for shape solver
+        F_liq_ = FT(0)
 
-    # ... and dont exceed the available number and mass of water droplets
-    dNdt = min(dNdt, N_ice / dt)
-    dLdt = min(dLdt, L_ice / dt)
-    return (; dNdt, dLdt)
+        # Get constants
+        (; ν_air, D_vapor, K_therm) = aps
+        L_f = TD.latent_heat_fusion(tps, Tₐ)
+        N_sc = ν_air / D_vapor
+
+        # Get the P3 diameter distribution...
+        D_th = D_th_helper(p3)
+        th = thresholds(p3, ρ_rim, F_rim)
+        (λ, N_0) = distribution_parameter_solver(
+            p3,
+            L_ice,
+            N_ice,
+            ρ_rim,
+            F_rim,
+            F_liq_,
+        )
+        N(D) = N′ice(p3, D, λ, N_0)
+        # ... and D_max for the integral
+        bound = get_ice_bound(p3, λ, FT(1e-6))
+
+        # check if we have a non-negligible portion of 
+        # the PSD above D_th
+        small = ifelse(bound < D_th, true, false)
+
+        # Ice particle terminal velocity
+        use_aspect_ratio = true
+        v(D) = p3_particle_terminal_velocity(
+            p3,
+            D,
+            Chen2022,
+            ρₐ,
+            F_rim,
+            F_liq,
+            th,
+            use_aspect_ratio,
+        )
+        # Reynolds number
+        N_Re(D) = D * v(D) / ν_air
+        # Ventillation factor
+        F_v(D) = p3.vent_a + p3.vent_b * N_sc^(1 / 3) * N_Re(D)^(1 / 2)
+        dmdD(D) = p3_dmdD(p3, D, F_rim, th)
+
+        f(D) =
+            4 * K_therm / L_f * (Tₐ - p3.T_freeze) * dmdD(D) / D * F_v(D) * N(D)
+        # Integrate
+        if small
+            (dLdt_rai, error) = QGK.quadgk(d -> f(d), 0, bound, rtol = FT(1e-6))
+        else
+            (dLdt_rai, error) = QGK.quadgk(d -> f(d), 0, D_th, rtol = FT(1e-6))
+            (dLdt_liq, error) =
+                QGK.quadgk(d -> f(d), D_th, bound, rtol = FT(1e-6))
+        end
+
+        # only consider melting (not fusion)
+        dLdt_rai = max(FT(0), dLdt_rai)
+        dLdt_liq = max(FT(0), dLdt_liq)
+
+        # compute sink of solid ice
+        dLdt_ice = dLdt_rai + dLdt_liq
+
+        # normalize dLdt_rai and dLdt_liq so that their sum does not
+        # exceed the available solid ice mass
+        if dLdt_ice > (L_ice / dt)
+            n = L_ice / (dLdt_ice)
+            dLdt_rai *= n
+            dLdt_liq *= n
+            # set change in solid ice mass = max
+            dLdt_ice = (L_ice / dt)
+        end
+
+        # compute change in L_rim, B_rim such that F_rim is unchanged
+        dLdt_rim = F_rim * dLdt_ice
+        ddtB_rim = dLdt_rim / ρ_rim
+
+
+        # change in total ice mass is equal to mass lost to rain
+        dLdt_p3_tot = dLdt_rai
+
+        # compute change of N_ice and N_rai proportional to the change of L_p3_tot
+        dNdt_ice = N_ice / L_ice * dLdt_p3_tot
+
+        # ... and don't exceed the available number
+        dNdt_ice = max(0, min(dNdt_ice, N_ice / dt))
+        dNdt_rai = dNdt_ice
+    end
+    return P3Rates{FT}(
+        dLdt_p3_tot = dLdt_p3_tot,
+        dLdt_rim = dLdt_rim,
+        dLdt_liq = dLdt_liq,
+        dLdt_rai = dLdt_rai,
+        dNdt_ice = dNdt_ice,
+        dNdt_rai = dNdt_rai,
+        ddtB_rim = ddtB_rim,
+    )
 end