GitHub

P3 internals

This page is a reference for the building blocks of the P3 functions on the API page, grouped by what they compute. Their names, arguments, and behavior can change without notice, so models should call the API functions instead.

Mass relation

P3Scheme.P3Particles evaluates these relations once at construction, from the rime mass fraction and the rime density. The graupel density follows from the density of the unrimed part, and the size thresholds of the mass regimes follow from the densities. The aspect ratio of a particle uses the material density of its regime.

CloudMicrophysics.P3Scheme.get_ρ_d — Function
get_ρ_d(mass, F_rim, ρ_rim)

Return the density of the unrimed portion of the particle [kg/m³].

The result is the exact solution as a function of the rime mass fraction F_rim, the mass power law parameters mass, and the rime density ρ_rim.

For the derivation of the numerically stable form used here, see the (P3 scheme documentation).

Arguments

Returns

  • ρ_d: density of the unrimed portion of the particle [kg/m³]

Examples

julia> import CloudMicrophysics.Parameters as CMP,
              ClimaParams as CP,
              CloudMicrophysics.P3Scheme as P3

julia> FT = Float64;

julia> mass = CMP.MassPowerLaw(CP.create_toml_dict(FT));

julia> F_rim, ρ_rim = FT(0.5), FT(916.7);

julia> ρ_d = P3.get_ρ_d(mass, F_rim, ρ_rim)
488.9120789986414
source
CloudMicrophysics.P3Scheme.get_ρ_g — Function
get_ρ_g(mass, F_rim, ρ_rim)

Return the density of total (deposition + rime) ice mass for graupel [kg/m³]

Arguments

  • mass: CMP.MassPowerLaw parameters
  • F_rim: rime mass fraction (L_rim / L_ice) [-]
  • ρ_rim: rime density (L_rim / B_rim) [kg/m³]

Returns

  • ρ_g: density of total (deposition + rime) ice mass for graupel [kg/m³]

Notes:

The density of the unrimed portion of the particle is get_ρ_d. See Eq. 16 in [25].

source
CloudMicrophysics.P3Scheme.ϕ_material_density — Function
ϕ_material_density(particles, D)

Return the material density of the solid of an ice particle at diameter D [kg/m³].

It is used in the aspect-ratio closure ϕᵢ. This is ρ_i in every mass regime except graupel (D_gr ≤ D < D_cr), where it is the graupel density ρ_g.

This is the density of the actual solid material, not the size-dependent effective density mᵢ / (π D³ / 6) returned by ice_density; see the P3 documentation for the distinction.

Arguments

  • particles: P3Particles
  • D: maximum particle dimension [m]
source

Slope solve

P3Scheme.P3Distribution solves for the slope with get_distribution_logλ. The solve matches log(L/N) from logLdivN, built from the mass moments of logmass_gamma_moment and the gamma moments of loggamma_moment, to the ratio of the ice mass and the bounded number. The shape parameter follows from the slope with get_μ, and the intercept from the number with get_logN₀. logN′ice evaluates the log of the number concentration at a diameter.

CloudMicrophysics.P3Scheme.get_distribution_logλ — Function
get_distribution_logλ(particles, ρq_ice, ρn_ice)

Solve for logλ, the log of the slope parameter of the P3 size distribution.

The solve uses the ice mass ρq_ice, the bounded ice number ρn_ice, and the mass relation of the P3Particles particles.

The assumed distribution is of the form

\[N′(D) = N₀ D^μ e^{-λD}\]

where N′(D) is the number concentration at diameter D and μ is the shape parameter. The shape parameter is parameterized, e.g. CMP.SlopePowerLaw or CMP.SlopeConstant.

This algorithm solves for logλ = log(λ) given ρq_ice and ρn_ice by solving the equations:

\[\begin{align*} \log(L) &= \log ∫_0^∞ m(D) N′(D)\ \mathrm{d}D, \\ \log(N) &= \log ∫_0^∞ N′(D)\ \mathrm{d}D, \\ \end{align*}\]

where m(D) is the mass of a particle at diameter D (see ice_mass). The procedure is described in detail in the P3 docs.

The solve uses the primal (as opposed to ForwardDiff.Dual) values of its arguments and returns a primal value, so logλ is constant under differentiation with respect to the state.

Arguments

  • particles: The P3Particles
  • ρq_ice: The ice mass concentration [kg/m³]
  • ρn_ice: The ice number concentration [1/m³], bounded by bounded_ice_number

The search bounds are logλ = 2 and the slope at which the mean particle mass equals the mass of a newly nucleated crystal. Without ice mass, the result is the upper bound, the smallest particles. With mass but without number, it is the lower bound, the largest particles.

source
CloudMicrophysics.P3Scheme.logmass_gamma_moment — Function
logmass_gamma_moment(particles, μ, logλ; n = 0)

Compute log(∫_0^∞ Dⁿ m(D) N′(D) dD) given the particles, μ, and logλ.

This is the log of the n-th moment of the mass-weighted PSD.

Arguments

  • particles: P3Particles
  • μ: The shape parameter [-]
  • logλ: The log of the slope parameter [log(1/m)]

Keyword arguments

  • n: The order of the moment, default is 0

Note:

  • For n = 0, this evaluates to log(L/N₀)
  • For n = 1, this evaluates to the (unnormalized) mass-weighted mean particle size, see D_m
source
CloudMicrophysics.P3Scheme.loggamma_moment — Function
loggamma_moment(μ, logλ; [k = 0], [scale = 1])

Compute log(scale ⋅ ∫_0^∞ G(D) D^k dD).

G(D) ≡ D^μ e^{-λD} is the (unnormalized) gamma kernel, k is an arbitrary exponent, and scale is a scale factor.

Arguments

  • μ: The PSD shape parameter [-]
  • logλ: The log of the PSD slope parameter [log(1/m)]

Keyword arguments

  • k: An arbitrary exponent [-], default is 0
  • scale: The scale factor [-], default is 1.

The implementation follows the same logic as loggamma_inc_moment, but with D₁ = 0 and D₂ = ∞, which implies q_D₁ = 1 and q_D₂ = 0.

source
CloudMicrophysics.P3Scheme.loggamma_inc_moment — Function
loggamma_inc_moment(D₁, D₂, μ, logλ, [k = 0], [scale = 1])

Compute log(Iᵏ) where Iᵏ is the following integral:

``I^k = ∫_{D₁}^{D₂} G(D) D^k dD``

$G(D) ≡ D^μ e^{-λD}$ is the (unnormalized) gamma kernel, and k is an arbitrary exponent.

If scale is provided, log(scale ⋅ Iᵏ) is returned.

With appropriate scaling, we can compute useful quantities like:

  • the k-th moment of the ice PSD, $M^k = N₀ I^k$
  • combined power law and moment weighted integrals, $∫_{D₁}^{D₂} (aD^b) D^n K(D) dD ≡ a I^(b + n)$

Arguments

  • D₁: The minimum diameter [m]
  • D₂: The maximum diameter [m]
  • μ: The PSD shape parameter [-]
  • logλ: The log of the PSD slope parameter [log(1/m)]
  • k: An arbitrary exponent [-], default is 0
  • scale: The scale factor [-], default is 1

Extended help

Implementation details

We can write ∫_D₁^D₂ G(D) D^k dD, where G(D) = D^μ e^{-λD} as: ∫_D₁^∞ G(D) D^k dD - ∫_D₂^∞ G(D) D^k dD with the transformation x = λD, and z = μ+k+1, each term can be written as: ∫_{Dᵢ}^∞ G(D) D^k dD = ∫_{λDᵢ}^∞ x^z e^{-x} dx / λ^z = Γ(z, λDᵢ) / λ^z where Γ(z, λDᵢ) = q ⋅ Γ(z) and q is the incomplete gamma function ratio given by (_, q) = UT.gamma_inc(z, x). This means that the integral ∫_{Dᵢ}^∞ G(D) D^k dD is computed as: Γ(z) ⋅ q / λ^z The full integral from D₁ to D₂ is then: Γ(z) ⋅ (q_D₁ - q_D₂) / λ^z In log-space, this is: - z log(λ) + logΓ(z) + log(q_D₁ - q_D₂)

See also gamma_inc_moment

source
CloudMicrophysics.P3Scheme.gamma_inc_moment — Function
gamma_inc_moment(D₁, D₂, p, α)

Compute the incomplete gamma moment ∫_{D₁}^{D₂} D^p e^{-α D} dD.

It equals α^{-(p+1)} Γ(p+1) [Q(p+1,αD₁) - Q(p+1,αD₂)], with Q the regularized upper incomplete gamma.

Returns 0 if D₂ ≤ D₁, and NaN if α ≤ 0

See also loggamma_inc_moment

source
CloudMicrophysics.P3Scheme.get_logN₀ — Function
get_logN₀(N_ice, μ, logλ)

Compute log(N₀) given N_ice, μ, and logλ,

    N  = N₀ ∫ G(D) dD
log N₀ = log N - log(∫G(D) dD)
       = log(N) - log( ∫D^μ e^{-λD} dD )
       = log(N) - log(M⁰)

Arguments

  • N_ice: The number concentration [1/m³]
  • μ: The shape parameter [-]
  • logλ: The log of the slope parameter [log(1/m)]

Without ice number, log(N₀) = -Inf, with zero derivatives.

source

Liquid-ice collisions

P3Scheme.bulk_liquid_ice_collision_sources builds the collision kernel, the local rime density, and the maximum freezing rate for the ice particles of the state. For each ice diameter, get_liquid_integrals integrates the kernel over the liquid size distribution, and ∫liquid_ice_collisions integrates the result over the ice size distribution and partitions the collected mass between freezing and shedding.

CloudMicrophysics.P3Scheme.volumetric_collision_rate_integrand — Function
volumetric_collision_rate_integrand(velocity_params, particles, ρₐ)

Return a function that computes the volumetric collision rate integrand for ice-liquid collisions [m³/s].

The returned function takes ice and liquid particle diameters as arguments.

Arguments

Returns

A function (D_ice, D_liq) -> E * K * |vᵢ - vₗ| where:

  • D_ice and D_liq are the (maximum) diameters of the ice and liquid particles
  • E is the collision efficiency
  • K is the collision cross section
  • vᵢ and vₗ are the terminal velocities of ice and liquid particles

Note that E, K, vᵢ and vₗ are all, in general, functions of D_ice and D_liq.

This function is a component of integrals like

\[∫ ∫ E * K * |vᵢ - vₗ| * N'_i * N'_l dD_i dD_l\]

source
CloudMicrophysics.P3Scheme.compute_max_freeze_rate — Function
compute_max_freeze_rate(aps, tps, velocity_params, particles, ρₐ, Tₐ)

Return a function max_freeze_rate(Dᵢ) of the maximum possible freezing rate [kg/s].

The rate is that of an ice particle of diameter Dᵢ [m]. It evaluates to 0 if T ≥ T_freeze.

Arguments

This rate represents the thermodynamic upper limit to collisional freezing, which occurs when the heat transfer from the ice particle to the environment is balanced by the latent heat of fusion.

From Eq (A7) in Musil (1970), [33].

source
CloudMicrophysics.P3Scheme.compute_local_rime_density — Function
compute_local_rime_density(velocity_params, particles, ρₐ, T)

Return a function ρ′_rim(Dᵢ, Dₗ) that computes the local rime density [kg/m³].

The arguments are the ice particle diameter Dᵢ [m] and the liquid particle diameter Dₗ [m].

Arguments

Returns

A function that computes the local rime density [kg/m³] using the equation:

\[ρ'_{rim} = a + b R_i + c R_i^2\]

where

\[R_i = -\frac{ 10^6 ⋅ D_{liq} ⋅ |v_{liq} - v_{ice}| }{ 2 T_{sfc} }\]

and $T_{sfc} < 0$ is the sub-zero surface temperature [°C], $D_{liq}$ is the liquid particle diameter [m], $v_{liq/ice}$ is the particle terminal velocity [m/s]. With the $10^6$ factor converting $D_{liq}$ from [m] to [μm], the units of $R_i$ are [μm m s⁻¹ °C⁻¹]. The units of $ρ'_{rim}$ are [kg/m³].

We assume for simplicity that $T_{sfc}$ equals $T$, the ambient air temperature. For real graupel, $T_{sfc}$ is slightly higher than $T$ due to latent heat release of freezing liquid particles onto the ice particle. Morrison & Milbrandt (2013) found little sensitivity to "realistic" increases in $T_{sfc}$.

See also LocalRimeDensity.

Extended help

Implementation follows Cober and List (1993), Eq. 16 and 17. See also the P3 fortran code, microphy_p3.f90, Line 3315-3323, which extends the range of the calculation to $R_i ≤ 12$, the upper limit of which then equals the solid bulk ice density, $ρ_ice = 916.7 kg/m^3$. At and above $T_{sfc} = 0$, the local rime density is $ρ_ice$.

Note that Morrison & Milbrandt (2015) [25] only uses this parameterization for collisions with cloud droplets. For rain drops, they use a value near the solid bulk ice density, $ρ^* = 900 kg/m^3$. We do not consider this distinction, and use this parameterization for all liquid particles.

source
CloudMicrophysics.P3Scheme.get_liquid_integrals — Function
get_liquid_integrals(n, ∂ₜV, m_liq, ρ′_rim, liq_bounds, quad)

Return a function liquid_integrals(Dᵢ) of the liquid particle integrals at ice diameter Dᵢ.

Arguments

  • n: liquid particle size distribution function n(D)
  • ∂ₜV: volumetric collision rate integrand function ∂ₜV(Dᵢ, D)
  • m_liq: liquid particle mass function m_liq(D)
  • ρ′_rim: local rime density function ρ′_rim(Dᵢ, D)
  • liq_bounds: integration bounds for liquid particles
  • quad: quadrature rule, a QuadratureRule

Notes

The function liquid_integrals(Dᵢ) returns a tuple (∂ₜN_col, ∂ₜM_col, ∂ₜB_col) of collision rates at Dᵢ, where:

  • ∂ₜN_col: number collision rate [1/s]
  • ∂ₜM_col: mass collision rate [kg/s]
  • ∂ₜB_col: rime volume collision rate [m³/s]
source
CloudMicrophysics.P3Scheme.∫liquid_ice_collisions — Function
∫liquid_ice_collisions(
    n_i, ∂ₜM_max, cloud_integrals, rain_integrals, ice_bounds, quad,
)

Computes the bulk collision rate integrands between ice and liquid particles.

Arguments

  • n_i: ice particle size distribution function n_i(D)
  • ∂ₜM_max: maximum freezing rate function ∂ₜM_max(Dᵢ)
  • cloud_integrals: inner liquid integrals for cloud particles, e.g. from get_liquid_integrals
  • rain_integrals: inner liquid integrals for rain particles, e.g. from get_liquid_integrals or get_liquid_integrals_rain_closed
  • ice_bounds: integration bounds for ice particles, from integral_bounds
  • quad: quadrature rule, a QuadratureRule

Returns

A 10-element vector of integrated rates, (QCFRZ, QCSHD, NCCOL, QRFRZ, QRSHD, NRCOL, ∫M_col, BCCOL, BRCOL, ∫𝟙_wet_M_col); see the method below for the definition of each entry.

source
∫liquid_ice_collisions(
    psd_c, psd_r, aps, tps, vel, m_liq, state,
    L_c, N_c, L_r, N_r, ρₐ, T,
)

Compute key liquid-ice collision rates and quantities. Used by bulk_liquid_ice_collision_sources.

Arguments

Returns

A tuple (QCFRZ, QCSHD, NCCOL, QRFRZ, QRSHD, NRCOL, ∫M_col, BCCOL, BRCOL, ∫𝟙_wet_M_col), where:

  1. QCFRZ - Cloud mass collision rate due to freezing [kg/s]
  2. QCSHD - Cloud mass collision rate due to shedding [kg/s]
  3. NCCOL - Cloud number collision rate [1/s]
  4. QRFRZ - Rain mass collision rate due to freezing [kg/s]
  5. QRSHD - Rain mass collision rate due to shedding [kg/s]
  6. NRCOL - Rain number collision rate [1/s]
  7. ∫M_col - Total collision rate [kg/s]
  8. BCCOL - Cloud rime volume source [m³/m³/s]
  9. BRCOL - Rain rime volume source [m³/m³/s]
  10. ∫𝟙_wet_M_col - Wet growth indicator [kg/s]
source

Quadrature

The P3 integrals split the integration range at the regime thresholds and at the kinks of the integrand, and sum a quadrature over each subinterval.

CloudMicrophysics.Quadrature.integrate — Function
integrate(f, a, b, quad)

Approximate the definite integral ∫ₐᵇ f(x) dx using the quadrature rule quad.

Mathematical Background

The integral is transformed to the standard interval [-1, 1]:

 y = (2x - (a+b)) / (b-a)    →    x = (b-a)y/2 + (a+b)/2

with Jacobian dx/dy = (b-a)/2, so that

 ∫ₐᵇ f(x) dx = (b-a)/2 ∫₋₁¹ f(x(y)) dy
             ≈ (b-a)/2 ∑ᵢ f(xᵢ) · inv_weight_fun(quad, yᵢ) · weight(quad, i)

where yᵢ = node(quad, i), xᵢ = (b-a)yᵢ/2 + (a+b)/2, and the inv_weight_fun/weight factors are rule-specific (see ChebyshevGauss and GaussLegendre).

Arguments

  • f: Function to integrate
  • a, b: Integration bounds. Note: if a ≥ b, or a or b is NaN, zero(f(a)) is returned.
  • quad: Quadrature scheme (a QuadratureRule).

Returns

Approximation to the definite integral ∫ₐᵇ f(x) dx

source
integrate(f, bnds, quad)

Integrate the function f over each subinterval of the integration bounds, bnds.

Arguments

  • f: Function to integrate
  • bnds: A tuple of bounds, (a, b, c, d, ...)
  • quad: Quadrature scheme (a QuadratureRule).

The integral is computed as the sum of the integrals over each subinterval, (a, b), (b, c), (c, d), ... — see subintervals.

source
CloudMicrophysics.Quadrature.subintervals — Function
subintervals(bnds::NTuple)

Pair adjacent elements of a flat boundary tuple into the consecutive subintervals they delimit:

subintervals((a, b, c, d)) ≡ ((a, b), (b, c), (c, d))

Base.front and Base.tail are zero-cost on fixed-length tuples, so the helper compiles to a direct tuple construction with no runtime dispatch.

source
CloudMicrophysics.P3Scheme.velocity_integral_bounds — Function
velocity_integral_bounds(particles, dist, v_term; p, moment_order = 0)

Compute the integration bounds for a velocity-weighted P3 integral.

This function returns the same integration bounds as integral_bounds, but with additional breakpoints that a terminal-velocity closure v_term may introduce, see velocity_breakpoints.

Arguments

See integral_bounds for the keyword arguments and the return value.

source