Add P3 terminal velocity functions

docs/src/API.md

@@ -66,12 +66,7 @@ Microphysics2M.conv_q_liq_to_q_rai
 P3Scheme
 P3Scheme.thresholds
 P3Scheme.distribution_parameter_solver
-```
-
-# Terminal Velocity
-```@docs
-TerminalVelocity
-TerminalVelocity.velocity_chen
+P3Scheme.ice_terminal_velocity
 ```
 
 # Aerosol model

docs/src/plots/P3TerminalVelocityPlots.jl

@@ -31,7 +31,7 @@ function get_values(
             ρ_r = ρ_rs[j]
 
             V_m[i, j] =
-                P3.terminal_velocity(p3, Chen2022, q, N, ρ_r, F_r, ρ_a)[2]
+                P3.ice_terminal_velocity(p3, Chen2022, q, N, ρ_r, F_r, ρ_a)[2]
             # get D_m in mm for plots
             D_m[i, j] = 1e3 * P3.D_m(p3, q, N, ρ_r, F_r)
         end

src/CloudMicrophysics.jl

@@ -14,13 +14,12 @@ include("Common.jl")
 include("Microphysics0M.jl")
 include("Microphysics1M.jl")
 include("Microphysics2M.jl")
-include("TerminalVelocity.jl")
+include("IceNucleation.jl")
 include("P3.jl")
 include("MicrophysicsNonEq.jl")
 include("AerosolModel.jl")
 include("AerosolActivation.jl")
 include("Nucleation.jl")
-include("IceNucleation.jl")
 include("PrecipitationSusceptibility.jl")
 include("ArtifactCalling.jl")
 

src/Microphysics2M.jl

@@ -504,6 +504,28 @@ function rain_self_collection_and_breakup(
     return (; sc, br)
 end
 
+"""
+    rain_particle_terminal_velocity(D, Chen2022, ρₐ)
+
+ - D - maximum particle dimension
+ - Chen2022 - a struct with terminal velocity parameters from Chen 2022
+ - ρₐ - air density
+
+Returns the terminal velocity of an individual rain drop as a function
+of its size (maximum dimension) following Chen 2022 velocity parametrization.
+Needed for numerical integrals in the P3 scheme.
+"""
+function rain_particle_terminal_velocity(
+    D::FT,
+    Chen2022::CMP.Chen2022VelTypeRain,
+    ρₐ::FT,
+) where {FT}
+    (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρₐ)
+
+    v = sum(@. sum(ai * D^bi * exp(-ci * D)))
+    return v
+end
+
 """
     rain_terminal_velocity(SB2006, vel, q_rai, ρ, N_rai)
 

src/P3.jl

@@ -14,7 +14,11 @@ import HCubature as HC
 import ClimaParams as CP
 import CloudMicrophysics.Parameters as CMP
 import CloudMicrophysics.Common as CO
-import CloudMicrophysics.TerminalVelocity as TV
+import CloudMicrophysics.HetIceNucleation as CM_HetIce
+import CloudMicrophysics.Microphysics2M as CM2
+
+import Thermodynamics as TD
+import Thermodynamics.Parameters as TDP
 
 const PSP3 = CMP.ParametersP3
 

src/P3_terminal_velocity.jl

@@ -1,5 +1,53 @@
 """
-    terminal_velocity(p3, Chen2022, q, N, ρ_r, F_r, ρ_a)
+    ϕᵢ(mᵢ, aᵢ, ρᵢ)
+
+ - mᵢ - particle mass
+ - aᵢ - particle area
+ - ρᵢ - ice density
+
+Returns the aspect ratio (ϕ) for an ice particle with given mass, area, and ice density.
+"""
+function ϕᵢ(mᵢ::FT, aᵢ::FT, ρᵢ::FT) where {FT}
+    # TODO - add some notes on how we derived it
+    # TODO - should we make use of other P3 properties like rimed fraction and volume?
+    return 16 * ρᵢ^2 * aᵢ^3 / (9 * π * mᵢ^2)
+end
+
+"""
+    ice_particle_terminal_velocity(D, Chen2022, ρₐ, mᵢ, aᵢ, ρᵢ)
+
+ - D - maximum particle dimension
+ - Chen2022 - a struct with terminal velocity parameters from Chen 2022
+ - ρₐ - air density
+ - mᵢ - particle mass
+ - aᵢ - particle area
+ - ρᵢ - ice density
+
+Returns the terminal velocity of a single ice particle as a function
+of its size (maximum dimension) using Chen 2022 parametrization.
+Needed for numerical integrals in the P3 scheme.
+"""
+function ice_particle_terminal_velocity(
+    D::FT,
+    Chen2022::CMP.Chen2022VelTypeSnowIce,
+    ρₐ::FT,
+    mᵢ::FT,
+    aᵢ::FT,
+    ρᵢ::FT,
+) where {FT}
+    if D <= Chen2022.cutoff
+        (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρₐ)
+    else
+        (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρₐ)
+    end
+    v = sum(@. sum(ai * D^bi * exp(-ci * D)))
+
+    κ = FT(-1 / 6)
+    return ϕᵢ(mᵢ, aᵢ, ρᵢ)^κ * v
+end
+
+"""
+    ice_terminal_velocity(p3, Chen2022, q, N, ρ_r, F_r, ρ_a)
 
  - p3 - a struct with P3 scheme parameters
  - Chen2022 - a struch with terminal velocity parameters as in Chen(2022)
@@ -9,10 +57,10 @@
  - F_r - rime mass fraction (q_rim/q_i)
  - ρ_a - density of air
 
- Returns the mass and number weighted fall speeds
- Eq C10 of Morrison and Milbrandt (2015)
+Returns the mass and number weighted fall speeds for ice following
+eq C10 of Morrison and Milbrandt (2015).
 """
-function terminal_velocity(
+function ice_terminal_velocity(
     p3::PSP3,
     Chen2022::CMP.Chen2022VelTypeSnowIce,
     q::FT,
@@ -52,6 +100,7 @@ function terminal_velocity(
     spheres_n(a, b, c) = (aₛ(a), bₛ(b), cₛ(c))
     spheres_m(a, b, c) = (aₛ_m(a), bₛ_m(b), cₛ(c))
 
+    # TODO - combine different aspect ratio computations
     aₙₛ(a) = aₛ(a) * (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ
     bₙₛ(b) = bₛ(b) + κ * (3 * p3.σ - 2 * p3.β_va)
 
@@ -183,3 +232,67 @@ function terminal_velocity(
     end
     return (v_n / N, v_m / q)
 end
+
+"""
+    velocity_difference(type, Dₗ, Dᵢ, p3, Chen2022, ρ_a, F_r, th)
+
+ - type - a struct containing the size distribution parameters of the particle colliding with ice
+ - Dₗ - maximum dimension of the particle colliding with ice
+ - Dᵢ - maximum dimension of ice particle
+ - p3 - a struct containing P3 parameters
+ - Chen2022 - a struct containing Chen 2022 velocity parameters
+ - ρ_a - density of air
+ - F_r - rime mass fraction (q_rim/ q_i)
+ - th - P3 particle properties thresholds
+
+Returns the absolute value of the velocity difference between an ice particle and
+cloud or rain drop as a function of their sizes. It uses Chen 2022 velocity
+parameterization for ice and rain and assumes no sedimentation of cloud droplets.
+Needed for P3 numerical integrals
+"""
+function velocity_difference(
+    pdf_r::Union{
+        CMP.RainParticlePDF_SB2006{FT},
+        CMP.RainParticlePDF_SB2006_limited{FT},
+    },
+    Dₗ::FT,
+    Dᵢ::FT,
+    p3::PSP3,
+    Chen2022::CMP.Chen2022VelType,
+    ρ_a::FT,
+    F_r::FT,
+    th,
+) where {FT}
+    return abs(
+        ice_particle_terminal_velocity(
+            Dᵢ,
+            Chen2022.snow_ice,
+            ρ_a,
+            p3_mass(p3, Dᵢ, F_r, th),
+            p3_area(p3, Dᵢ, F_r, th),
+            p3.ρ_i,
+        ) - CM2.rain_particle_terminal_velocity(Dₗ, Chen2022.rain, ρ_a),
+    )
+end
+function velocity_difference(
+    pdf_c::CMP.CloudParticlePDF_SB2006{FT},
+    Dₗ::FT,
+    Dᵢ::FT,
+    p3::PSP3,
+    Chen2022::CMP.Chen2022VelType,
+    ρ_a::FT,
+    F_r::FT,
+    th,
+) where {FT}
+    # velocity difference for cloud collisions
+    return abs(
+        ice_particle_terminal_velocity(
+            Dᵢ,
+            Chen2022.snow_ice,
+            ρ_a,
+            p3_mass(p3, Dᵢ, F_r, th),
+            p3_area(p3, Dᵢ, F_r, th),
+            p3.ρ_i,
+        ),
+    )
+end