fix: limiting behavior of 2M and P3 PSDs

Author: haakon-e

src/Microphysics2M.jl

@@ -64,6 +64,7 @@ Return the parameters of the rain drop diameter distribution
  - A `NamedTuple` with the fields `(; N₀r, Dr_mean, xr_mean)`
 """
 function pdf_rain_parameters(pdf_r::CMP.RainParticlePDF_SB2006_notlimited, qᵣ, ρₐ, Nᵣ)
+    (iszero(Nᵣ) && iszero(qᵣ)) && return (; N₀r = zero(Nᵣ), Dr_mean = zero(qᵣ), xr_mean = zero(qᵣ))
     (; ρw) = pdf_r
     Lᵣ = ρₐ * qᵣ
 
@@ -75,11 +76,12 @@ function pdf_rain_parameters(pdf_r::CMP.RainParticlePDF_SB2006_notlimited, qᵣ,
     return (; N₀r, Dr_mean, xr_mean)
 end
 function pdf_rain_parameters(pdf_r::CMP.RainParticlePDF_SB2006_limited, qᵣ, ρₐ, Nᵣ)
+    (iszero(Nᵣ) && iszero(qᵣ)) && return (; N₀r = zero(Nᵣ), Dr_mean = zero(qᵣ), xr_mean = zero(qᵣ))
     (; xr_min, xr_max, N0_min, N0_max, λ_min, λ_max, ρw) = pdf_r
     Lᵣ = ρₐ * max(0, qᵣ)
 
     # Sequence of limiting steps in Seifert and Beheng 2006:
-    x̃r = clamp(Lᵣ / Nᵣ, xr_min, xr_max)  # Eq. (94)  # TODO: Ill-defined for 0 / 0
+    x̃r = clamp(Lᵣ / Nᵣ, xr_min, xr_max)  # Eq. (94)
     N₀r = clamp(Nᵣ * ∛(π * ρw / x̃r), N0_min, N0_max)  # Eq. (95)
     λr = clamp(∜(π * ρw * N₀r / Lᵣ), λ_min, λ_max)  # Eq. (96)
     xr_mean = clamp(Lᵣ * λr / N₀r, xr_min, xr_max)  # Eq. (97)
@@ -153,6 +155,7 @@ That is,
  - `(logA, logB)`: Log of the parameters of the generalized gamma distribution
 """
 function log_pdf_cloud_parameters_mass(pdf_c, q, ρₐ, N)
+    (N < eps(one(N)) && q < eps(one(q))) && return (log(zero(N)), log(1 / zero(q)))
     (; νc, μc) = pdf_c
     L = ρₐ * q
     logx̄ = log(L / N)
@@ -242,7 +245,7 @@ Return `n(D)`, a function that computes the size distribution for rain particles
 """
 function DT.size_distribution(pdf::CMP.RainParticlePDF_SB2006, q, ρₐ, N)
     (; N₀r, Dr_mean) = pdf_rain_parameters(pdf, q, ρₐ, N)
-    return n(D) = N₀r * exp(-D / Dr_mean)
+    return n(D) = iszero(N₀r) ? zero(D) : N₀r * exp(-D / Dr_mean)
 end
 
 """
@@ -263,7 +266,7 @@ The size distribution is given by:
 """
 function DT.size_distribution(pdf::CMP.CloudParticlePDF_SB2006, q, ρₐ, N)
     (; logN₀c, λc, νcD, μcD) = pdf_cloud_parameters(pdf, q, ρₐ, N)
-    return n(D) = exp(logN₀c + νcD * log(D) - λc * D^μcD)
+    return n(D) = logN₀c == -Inf ? zero(D) : exp(logN₀c + νcD * log(D) - λc * D^μcD)
 end
 
 """
@@ -303,6 +306,7 @@ function get_size_distribution_bounds(
     q, ρₐ, N, p = eps(FT),
 ) where {FT}
     (; Dr_mean) = pdf_rain_parameters(pdf, q, ρₐ, N)
+    iszero(Dr_mean) && return (FT(0), FT(0))
     D_min = DT.exponential_quantile(Dr_mean, p)
     D_max = DT.exponential_quantile(Dr_mean, 1 - p)
     return D_min, D_max

src/P3_size_distribution.jl

@@ -225,6 +225,7 @@ where `m(D)` is the mass of a particle at diameter `D` (see [`ice_mass`](@ref)).
 function get_distribution_logλ(
     state::P3State{FT}; logλ_min = log(1e1), logλ_max = log(1e7),
 ) where {FT}
+    (iszero(state.N_ice) || iszero(state.L_ice)) && return log(zero(FT))
     target_log_LdN = log(state.L_ice) - log(state.N_ice)
 
     shape_problem(logλ) = logLdivN(state, logλ) - target_log_LdN

test/microphysics2M_tests.jl

@@ -143,6 +143,34 @@ function test_microphysics2M(FT)
     end
 
     # 2M_microphysics - Seifert and Beheng 2006 double moment scheme tests
+    TT.@testset "Seifert and Beheng 2006 - PDF parameters limiting behavior" begin
+        N = 0.0
+        q = 0.0
+        ρₐ = 1.2
+        # limited rain drop size distribution
+        params = CM2.pdf_rain_parameters(SB2006.pdf_r, q, ρₐ, N)
+        TT.@test all(iszero, params)
+        n = CM2.size_distribution(SB2006.pdf_r, q, ρₐ, N)
+        TT.@test all(iszero, (n(0), n(0.1), n(Inf)))
+        bnds = CM2.get_size_distribution_bounds(SB2006.pdf_r, q, ρₐ, N)
+        TT.@test all(iszero, bnds)
+        # not limited rain drop size distribution
+        params = CM2.pdf_rain_parameters(SB2006_no_limiters.pdf_r, q, ρₐ, N)
+        TT.@test all(iszero, params)
+        n = CM2.size_distribution(SB2006_no_limiters.pdf_r, q, ρₐ, N)
+        TT.@test all(iszero, (n(0), n(0.1), n(Inf)))
+        bnds = CM2.get_size_distribution_bounds(SB2006_no_limiters.pdf_r, q, ρₐ, N)
+        TT.@test all(iszero, bnds)
+        # cloud drop size distribution
+        logA, logB = CM2.log_pdf_cloud_parameters_mass(SB2006.pdf_c, q, ρₐ, N)
+        TT.@test logA == -Inf
+        TT.@test logB == Inf
+        A, B = CM2.pdf_cloud_parameters_mass(SB2006.pdf_c, q, ρₐ, N)
+        TT.@test A == 0
+        TT.@test B == Inf
+        n = CM2.size_distribution(SB2006.pdf_c, q, ρₐ, N)
+        TT.@test all(iszero, (n(0), n(0.1), n(Inf)))
+    end
     TT.@testset "limiting lambda_r and x_r - Seifert and Beheng 2006" begin
         #setup
         q_rai = [FT(0), FT(1e-3), FT(1e-4), FT(1e-2)]

test/p3_tests.jl

@@ -178,6 +178,11 @@ function test_shape_solver(FT)
         params = CMP.ParametersP3(FT; slope_law)
 
         @testset "Shape parameters - nonlinear solver" begin
+            # -- First, test limiting behavior: `N_ice = L_ice = 0` --
+            state = P3.get_state(params; F_rim = FT(0.5), ρ_rim = FT(500), L_ice = FT(0), N_ice = FT(0))
+            logλ = P3.get_distribution_logλ(state)
+            @test logλ == -Inf
+            # --
 
             # initialize test values:
             ep = 1 #1e4 * eps(FT)