refactor: liquid (cloud/rain) terminal velocity

Author: haakon-e

docs/src/plots/TerminalVelocity.jl

@@ -47,30 +47,6 @@ function aspect_ratio_snow_1M_prolate(snow::CMP.Snow, D::FT) where {FT <: Real}
     return 16 * ρᵢ^2 * aᵢ^3 / (9 * FT(π) * mᵢ^2)
 end
 
-"""
-    rain_terminal_velocity_individual_Chen(velo_scheme ρₐ, ρᵢ, D_r)
-
- - `velo_scheme` - structs with free parameters
- - `ρₐ` - air density
- - `ρᵢ` - apparent density of ice particles
- - `D` - particle diameter
-
-Returns the fall velocity of a raindrops or ice crystals from Chen et al 2022
-"""
-function rain_terminal_velocity_individual_Chen(
-    velo_scheme::CMP.Chen2022VelTypeRain,
-    ρₐ::FT,
-    D::FT, #in m
-) where {FT <: Real}
-    ai, bi, ci = CMO.Chen2022_vel_coeffs_B1(velo_scheme, ρₐ)
-
-    v = 0
-    for i in 1:3
-        v += ai[i] * D^bi[i] * exp(-D * ci[i])
-    end
-    return v
-end
-
 function ice_terminal_velocity_individual_Chen(
     ice::CMP.CloudIce,
     velo_scheme::CMP.Chen2022VelTypeSmallIce,
@@ -181,22 +157,23 @@ q_range = range(0, stop = 5 * 1e-3, length = 100)
 
 #! format: off
 # velocity values for cloud particle sizes
+v_term_rain = CMO.liquid_particle_terminal_velocity(Chen2022.rain, ρ_air)
 SB_rain_small = [rain_terminal_velocity_individual_SB(SB2006Vel, ρ_air, D_r)                 for D_r in D_r_range_small]
 M1_rain_small = [terminal_velocity_individual_1M(Blk1MVel.rain, ρ_air, D_r)                  for D_r in D_r_range_small]
 M1_snow_small = [terminal_velocity_individual_1M(Blk1MVel.snow, ρ_air, D_r)                  for D_r in D_r_range_small]
-Ch_liq_small =  [rain_terminal_velocity_individual_Chen(Chen2022.rain, ρ_air, D_r)           for D_r in D_r_range_small]
+Ch_liq_small = v_term_rain.(D_r_range_small)
 Ch_ice_small =  [ice_terminal_velocity_individual_Chen(ice, Chen2022.small_ice, ρ_air, D_r)        for D_r in D_r_range_small]
-Ch_rain_small = [rain_terminal_velocity_individual_Chen(Chen2022.rain, ρ_air, D_r)           for D_r in D_r_range_small]
+Ch_rain_small = v_term_rain.(D_r_range_small)
 Ch_snow_small = [snow_terminal_velocity_individual_Chen(snow, Chen2022.large_ice, ρ_air, D_r) for D_r in D_r_range_small]
 Ch_snow_small_oblate = [snow_terminal_velocity_individual_Chen_oblate(snow, Chen2022.large_ice, ρ_air, D_r) for D_r in D_r_range_small]
 Ch_snow_small_prolate = [snow_terminal_velocity_individual_Chen_prolate(snow, Chen2022.large_ice, ρ_air, D_r) for D_r in D_r_range_small]
 # velocity values for precip particle sizes
 SB_rain = [rain_terminal_velocity_individual_SB(SB2006Vel, ρ_air, D_r)                 for D_r in D_r_range]
 M1_rain = [terminal_velocity_individual_1M(Blk1MVel.rain, ρ_air, D_r)                  for D_r in D_r_range]
 M1_snow = [terminal_velocity_individual_1M(Blk1MVel.snow, ρ_air, D_r)                  for D_r in D_r_range]
-Ch_liq =  [rain_terminal_velocity_individual_Chen(Chen2022.rain, ρ_air, D_r)           for D_r in D_r_range]
+Ch_liq = v_term_rain.(D_r_range)
 Ch_ice =  [ice_terminal_velocity_individual_Chen(ice, Chen2022.small_ice, ρ_air, D_r)        for D_r in D_r_range]
-Ch_rain = [rain_terminal_velocity_individual_Chen(Chen2022.rain, ρ_air, D_r)           for D_r in D_r_range]
+Ch_rain = v_term_rain.(D_r_range)
 Ch_snow = [snow_terminal_velocity_individual_Chen(snow, Chen2022.large_ice, ρ_air, D_r) for D_r in D_r_range]
 Ch_snow_oblate  = [snow_terminal_velocity_individual_Chen_oblate(snow, Chen2022.large_ice, ρ_air, D_r) for D_r in D_r_range]
 Ch_snow_prolate = [snow_terminal_velocity_individual_Chen_prolate(snow, Chen2022.large_ice, ρ_air, D_r) for D_r in D_r_range]

src/Common.jl

@@ -336,6 +336,34 @@ function Chen2022_exponential_pdf(a::FT, b::FT, c::FT, λ_inv::FT, k::Int) where
     return a * exp(δ * log(1 / λ_inv) - (b + δ) * log(1 / λ_inv + c)) * SF.gamma(b + δ) / SF.gamma(δ)
 end
 
+"""
+    liquid_particle_terminal_velocity(velocity_params, ρₐ)
+
+Compute the terminal velocity of a liquid particle as a function of its size 
+    (maximum dimension, `D`) using the Chen 2022 parametrization.
+
+# Arguments
+- `velocity_params`: a struct with terminal velocity parameters from Chen 2022
+- `ρₐ`: air density [kg/m³]
+
+# Returns
+- The terminal velocity of a liquid particle as a function of its size (maximum dimension)
+    following Chen 2022 velocity parametrization.
+    
+Needed for numerical integrals in the P3 scheme.
+
+!!! note
+    We use the same terminal velocity parametrization for cloud and rain water.
+"""
+function liquid_particle_terminal_velocity(velocity_params::CMP.Chen2022VelTypeRain, ρₐ)
+    (ai, bi, ci) = Chen2022_vel_coeffs_B1(velocity_params, ρₐ)
+    v_term(D) = sum(@. sum(ai * D^bi * exp(-ci * D)))
+    return v_term
+end
+liquid_particle_terminal_velocity(velocity_params::CMP.Chen2022VelType, ρₐ) =
+    liquid_particle_terminal_velocity(velocity_params.rain, ρₐ)
+
+
 """
     volume_sphere_D(D)
 

src/Microphysics2M.jl

@@ -590,29 +590,6 @@ function rain_self_collection_and_breakup(
     return (; sc, br)
 end
 
-"""
-    rain_particle_terminal_velocity(Chen2022, ρₐ)
-
-Compute 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.
-
-# Arguments
- - `Chen2022`: terminal velocity parameters from Chen 2022, [`CMP.Chen2022VelTypeRain`](@ref)
- - `ρₐ`: 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(Chen2022::CMP.Chen2022VelTypeRain, ρₐ)
-    (ai, bi, ci) = CO.Chen2022_vel_coeffs_B1(Chen2022, ρₐ)
-
-    v(D) = sum(@. sum(ai * D^bi * exp(-ci * D)))
-    return v
-end
-
 """
     rain_terminal_velocity(SB2006, vel, q_rai, ρ, N_rai)
 

test/p3_tests.jl

@@ -2,6 +2,7 @@ using Test: @testset, @test, @test_throws, @test_broken
 import CloudMicrophysics.P3Scheme as P3
 import CloudMicrophysics.Parameters as CMP
 import CloudMicrophysics.Microphysics2M as CM2
+import CloudMicrophysics.Common as CO
 import Thermodynamics as TD
 import ClimaParams as CP
 
@@ -236,9 +237,9 @@ function test_particle_terminal_velocities(FT)
     @testset "Smoke tests for cloud/rain particle terminal vel from Chen 2022" begin
         Ds = range(FT(1e-6), stop = FT(1e-5), length = 5)  # TODO: Add tests for larger sizes
         expected = [0.002508, 0.009156, 0.01632, 0.02377, 0.03144]
-        vel_func = CM2.rain_particle_terminal_velocity(Chen2022.rain, ρ_a)
+        v_term = CO.liquid_particle_terminal_velocity(Chen2022.rain, ρ_a)
         for i in axes(Ds, 1)
-            vel = vel_func(Ds[i])
+            vel = v_term(Ds[i])
             @test vel >= 0
             @test vel ≈ expected[i] rtol = 1e-3
         end