Cleanup in limiters

src/Microphysics2M.jl

@@ -66,11 +66,11 @@ function pdf_cloud(
 ) where {FT}
 
     (; νc, μc, ρw) = pdf_c
-    ϵ_N = FT(10)     # TODO - should it be a parameter
+    #ϵ_N = FT(10)     # TODO - should it be a parameter
 
     χ = 9 # convert mean droplet mass from kg to μg
 
-    if qₗ < eps(FT) || Nₗ < ϵ_N
+    if qₗ < eps(FT) || Nₗ < eps(FT)
         return (
             xc = FT(0),
             χ,
@@ -116,9 +116,9 @@ function pdf_rain(
     Nᵣ::FT,
 ) where {FT}
     (; νr, μr, ρw) = pdf_r
-    ϵ_N = FT(10)            # TODO - should it be a parameter?
+    #ϵ_N = FT(10)            # TODO - should it be a parameter?
 
-    if qᵣ < eps(FT) || Nᵣ < ϵ_N
+    if qᵣ < eps(FT) || Nᵣ < eps(FT)
         return (λr = FT(0), xr = FT(0), Ar = FT(0), Br = FT(0))
     else
         λr = (Nᵣ * FT(π) * ρw / ρₐ / qᵣ)^FT(1 / 3)
@@ -133,24 +133,24 @@ function pdf_rain(
 end
 function pdf_rain(
     pdf_r::CMP.RainParticlePDF_SB2006_limited{FT},
-    q_rai::FT,
-    ρ::FT,
-    N_rai::FT,
+    qᵣ::FT,
+    ρₐ::FT,
+    Nᵣ::FT,
 ) where {FT}
     (; νr, μr, xr_min, xr_max, N0_min, N0_max, λ_min, λ_max, ρw) = pdf_r
 
-    L_rai = ρ * q_rai
-    xr_0 = L_rai / N_rai
+    qᵣ = qᵣ > FT(0) ? qᵣ : FT(0)
+
+    Lᵣ = ρₐ * qᵣ
+    xr_0 = Nᵣ > eps(FT) ? Lᵣ / Nᵣ : FT(0)
     xr_hat = max(xr_min, min(xr_max, xr_0))
-    N0 = max(N0_min, min(N0_max, N_rai * (FT(π) * ρw / xr_hat)^FT(1 / 3)))
-    λr = max(λ_min, min(λ_max, (FT(π) * ρw * N0 / L_rai)^FT(1 / 4)))
-    xr = max(xr_min, min(xr_max, L_rai * λr / N0))
-    αr = N_rai * λr
+    N0 = max(N0_min, min(N0_max, Nᵣ * (FT(π) * ρw / xr_hat)^FT(1 / 3)))
+    λr = max(λ_min, min(λ_max, (FT(π) * ρw * N0 / Lᵣ)^FT(1 / 4)))
+    xr = max(xr_min, min(xr_max, Lᵣ * λr / N0))
+    αr = Nᵣ * λr
 
-    Br =
-        (xr == 0) ? FT(0) :
-        (SF.gamma(FT(νr + 1) / μr) * xr / SF.gamma(FT(νr + 2) / μr))^(-μr)
-    Ar = μr * N_rai * Br^(FT(νr + 1) / μr) / SF.gamma(FT(νr + 1) / μr)
+    Br = (SF.gamma(FT(νr + 1) / μr) * xr / SF.gamma(FT(νr + 2) / μr))^(-μr)
+    Ar = μr * Nᵣ * Br^(FT(νr + 1) / μr) / SF.gamma(FT(νr + 1) / μr)
 
     return (; λr, αr, xr, Ar, Br)
 end
@@ -176,7 +176,7 @@ function autoconversion(
     N_liq,
 ) where {FT}
 
-    if q_liq < eps(FT)
+    if q_liq < eps(FT) || N_liq < eps(FT)
         return LiqRaiRates{FT}()
     end
 
@@ -221,7 +221,7 @@ collisions between raindrops and cloud droplets (accretion) for `scheme == SB200
 """
 function accretion((; accr)::CMP.SB2006{FT}, q_liq, q_rai, ρ, N_liq) where {FT}
 
-    if q_liq < eps(FT) || q_rai < eps(FT)
+    if q_liq < eps(FT) || q_rai < eps(FT) || N_liq < eps(FT)
         return LiqRaiRates{FT}()
     end
 
@@ -326,19 +326,15 @@ function rain_self_collection(
     N_rai::FT,
 ) where {FT}
 
-    if q_rai < eps(FT)
+    if q_rai < eps(FT) || N_rai < eps(FT)
         return FT(0)
     end
 
     (; krr, κrr, d) = self
     (; ρ0, ρw) = pdf
 
     L_rai = ρ * q_rai
-    #λr =
-    #    pdf_rain(pdf, q_rai, ρ, N_rai).λr *
-    #    (SF.gamma(FT(4)) / FT(π) / ρw)^FT(1 / 3)
     λr = pdf_rain(pdf, q_rai, ρ, N_rai).Br
-
     dN_rai_dt_sc = -krr * N_rai * L_rai * sqrt(ρ0 / ρ) * (1 + κrr / λr)^d
 
     return dN_rai_dt_sc
@@ -368,7 +364,7 @@ function rain_breakup(
     dN_rai_dt_sc::FT,
 ) where {FT}
 
-    if q_rai < eps(FT)
+    if q_rai < eps(FT) || N_rai < eps(FT)
         return FT(0)
     end
     (; Deq, Dr_th, kbr, κbr) = brek
@@ -431,12 +427,15 @@ function rain_terminal_velocity(
     N_rai::FT,
 ) where {FT}
     # TODO: Input argument list needs to be redesigned
-    if q_rai < eps(FT)
-        return (FT(0), FT(0))
-    end
+
     λr = pdf_rain(pdf_r, q_rai, ρ, N_rai).λr
-    vt0 = max(FT(0), sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr)))
-    vt1 = max(FT(0), sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr)^FT(4)))
+
+    vt0 =
+        N_rai < eps(FT) ? FT(0) :
+        max(FT(0), sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr)))
+    vt1 =
+        q_rai < eps(FT) ? FT(0) :
+        max(FT(0), sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr)^FT(4)))
     return (vt0, vt1)
 end
 function rain_terminal_velocity(
@@ -446,9 +445,6 @@ function rain_terminal_velocity(
     ρ::FT,
     N_rai::FT,
 ) where {FT}
-    if q_rai < eps(FT)
-        return (FT(0), FT(0))
-    end
     # coefficients from Table B1 from Chen et. al. 2022
     aiu, bi, ciu = CO.Chen2022_vel_coeffs_small(vel, ρ)
     # size distribution parameter
@@ -458,8 +454,8 @@ function rain_terminal_velocity(
     vt0 = sum(CO.Chen2022_vel_add.(aiu, bi, ciu, λ, 0))
     vt3 = sum(CO.Chen2022_vel_add.(aiu, bi, ciu, λ, 3))
 
-    vt0 = max(FT(0), vt0)
-    vt3 = max(FT(0), vt3)
+    vt0 = N_rai < eps(FT) ? FT(0) : max(FT(0), vt0)
+    vt3 = q_rai < eps(FT) ? FT(0) : max(FT(0), vt3)
     # It should be (ϕ^κ * vt0, ϕ^κ * vt3), but for rain drops ϕ = 1 and κ = 0
     return (vt0, vt3)
 end
@@ -507,7 +503,7 @@ function rain_evaporation(
     evap_rate_1 = FT(0)
     S = TD.supersaturation(tps, q, ρ, T, TD.Liquid())
 
-    if (q_rai > FT(0) && S < FT(0))
+    if ((q_rai > eps(FT) || N_rai > eps(FT)) && S < FT(0))
 
         (; ν_air, D_vapor) = aps
         (; av, bv, α, β, ρ0) = evap
@@ -534,6 +530,9 @@ function rain_evaporation(
 
         evap_rate_0 = min(FT(0), FT(2) * FT(π) * G * S * N_rai * Dr * Fv0 / xr)
         evap_rate_1 = min(FT(0), FT(2) * FT(π) * G * S * N_rai * Dr * Fv1 / ρ)
+
+        evap_rate_0 = N_rai < eps(FT) ? FT(0) : evap_rate_0
+        evap_rate_1 = q_rai < eps(FT) ? FT(0) : evap_rate_1
     end
 
     return (; evap_rate_0, evap_rate_1)
@@ -578,19 +577,19 @@ function radar_reflectivity(
     χ = pdf_cloud(pdf_c, q_liq, ρ_air, N_liq).χ
 
     Zc =
-        Bc == 0 ? FT(0) :
+        Bc < eps(FT) ? FT(0) :
         Ac *
         C^(νc + 1) *
         (Bc * C^μc)^(-(3 + νc) / μc) *
         SF.gamma((3 + νc) / μc) / μc * FT(10)^(-2 * χ)
     Zr =
-        Br == 0 ? FT(0) :
+        Br < eps(FT) ? FT(0) :
         Ar *
         C^(νr + 1) *
         (Br * C^μr)^(-(3 + νr) / μr) *
         SF.gamma((3 + νr) / μr) / μr
 
-    return Zc + Zr == 0 ? FT(-135) : 10 * (log10(Zc + Zr) - log_10_Z₀)
+    return max(FT(-150), 10 * (log10(Zc + Zr) - log_10_Z₀))
 end
 
 """
@@ -655,7 +654,7 @@ function effective_radius(
         (Br * C^μr)^(-(5 + 3 * νr) / (3 * μr)) *
         SF.gamma((5 + 3 * νr) / (3 * μr)) / μr
 
-    return M2_c + M2_r == 0 ? FT(0) : (M3_c + M3_r) / (M2_c + M2_r)
+    return M2_c + M2_r <= eps(FT) ? FT(0) : (M3_c + M3_r) / (M2_c + M2_r)
 end
 
 """
@@ -681,7 +680,7 @@ function effective_radius_Liu_Hallet_97(
 
     k = FT(0.8)
     r_vol =
-        ((N_liq + N_rai) == 0) ? FT(0) :
+        ((N_liq + N_rai) < eps(FT)) ? FT(0) :
         (
             (FT(3) * (q_liq + q_rai) * ρ_air) /
             (FT(4) * π * ρ_w * (N_liq + N_rai))

test/Project.toml

@@ -13,6 +13,7 @@ EnsembleKalmanProcesses = "aa8a2aa5-91d8-4396-bcef-d4f2ec43552d"
 EvoTrees = "f6006082-12f8-11e9-0c9c-0d5d367ab1e5"
 Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9"
+HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -21,8 +22,8 @@ MLJFlux = "094fc8d1-fd35-5302-93ea-dabda2abf845"
 MLJModels = "d491faf4-2d78-11e9-2867-c94bc002c0b7"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"

test/microphysics2M_tests.jl

@@ -150,8 +150,8 @@ function test_microphysics2M(FT)
     # 2M_microphysics - Seifert and Beheng 2006 double moment scheme tests
     TT.@testset "limiting lambda_r and x_r - Seifert and Beheng 2006" begin
         #setup
-        q_rai = [FT(0), FT(1e-4), FT(1e-2)]
-        N_rai = [FT(1e1), FT(1e3), FT(1e5)]
+        q_rai = [FT(0), FT(1e-3), FT(1e-4), FT(1e-2)]
+        N_rai = [FT(1e1), FT(1e1), FT(1e3), FT(1e5)]
         ρ = FT(1)
 
         (; xr_min, xr_max, λ_min, λ_max) = SB2006.pdf_r
@@ -162,7 +162,6 @@ function test_microphysics2M(FT)
                 λ = CM2.pdf_rain(SB2006.pdf_r, qr, ρ, Nr).λr
                 xr = CM2.pdf_rain(SB2006.pdf_r, qr, ρ, Nr).xr
 
-                #test
                 TT.@test λ_min <= λ <= λ_max
                 TT.@test xr_min <= xr <= xr_max
             end
@@ -322,9 +321,8 @@ function test_microphysics2M(FT)
             sc_br_rai =
                 CM2.rain_self_collection_and_breakup(SB, q_rai, ρ, N_rai)
 
-            λr =
-                CM2.pdf_rain(SB.pdf_r, q_rai, ρ, N_rai).λr *
-                FT(6 / π / 1000)^FT(1 / 3)
+            λr = CM2.pdf_rain(SB.pdf_r, q_rai, ρ, N_rai).Br
+
             dNrdt_sc =
                 -krr * N_rai * ρ * q_rai * (1 + κrr / λr)^-5 * sqrt(ρ0 / ρ)
 
@@ -394,12 +392,13 @@ function test_microphysics2M(FT)
             TT.@test vt_rai isa Tuple
             TT.@test vt_rai[1] ≈ vt0 rtol = 1e-6
             TT.@test vt_rai[2] ≈ vt1 rtol = 1e-6
+
             TT.@test CM2.rain_terminal_velocity(
                 SB,
                 SB2006Vel,
-                FT(0),
+                q_rai,
                 ρ,
-                N_rai,
+                FT(0),
             )[1] ≈ 0 atol = eps(FT)
             TT.@test CM2.rain_terminal_velocity(
                 SB,
@@ -432,9 +431,9 @@ function test_microphysics2M(FT)
             TT.@test CM2.rain_terminal_velocity(
                 SB,
                 Chen2022Vel,
-                FT(0),
+                q_rai,
                 ρ,
-                N_rai,
+                FT(0),
             )[1] ≈ 0 atol = eps(FT)
             TT.@test CM2.rain_terminal_velocity(
                 SB,
@@ -492,9 +491,9 @@ function test_microphysics2M(FT)
                 aps,
                 tps,
                 q,
-                FT(0),
+                q_rai,
                 ρ,
-                N_rai,
+                FT(0),
                 T,
             ).evap_rate_0 ≈ 0 atol = eps(FT)
             TT.@test CM2.rain_evaporation(
@@ -510,38 +509,30 @@ function test_microphysics2M(FT)
         end
     end
 
-    TT.@testset "2M_microphysics - Seifert and Beheng 2006 radar reflectivity" begin
+    TT.@testset "2M_microphysics - Seifert and Beheng 2006 effective radius and reflectivity" begin
         #setup
-        ρ_air = FT(1)
-        q_liq = FT(2.128e-4)
-        N_liq = FT(15053529)
-        q_rai = FT(1.573e-4)
-        N_rai = FT(510859)
+        ρₐ = FT(1)
 
-        for SB in [SB2006, SB2006_no_limiters]
-            #action
-            rr = CM2.radar_reflectivity(SB, q_liq, q_rai, N_liq, N_rai, ρ_air)
+        q_liq = [FT(2.128e-4), FT(2.128e-20), FT(1.6e-12), FT(0), FT(1.037e-25)]
+        N_liq = [FT(15053529), FT(3), FT(5512), FT(0), FT(5.225e-12)]
+        q_rai = [FT(1.573e-4), FT(1.573e-4), FT(1.9e-15), FT(0), FT(2.448e-27)]
+        N_rai = [FT(510859), FT(510859), FT(0), FT(0), FT(5.136e-18)]
 
-            # test
-            TT.@test rr ≈ FT(-13) atol = FT(0.5)
-        end
-    end
+        # reference values
+        rr = [FT(-12.561951), FT(-12.579899), FT(-150), FT(-150), FT(-150)]
+        reff = [FT(2.319383e-5), FT(6.91594e-5), FT(0), FT(0), FT(0)]
 
-    TT.@testset "2M_microphysics - Seifert and Beheng 2006 effective radius" begin
-        #setup
-        ρ_air = FT(1)
-        ρ_w = FT(1000)
-        q_liq = FT(2.128e-4)
-        N_liq = FT(15053529)
-        q_rai = FT(1.573e-4)
-        N_rai = FT(510859)
+        for (qₗ, Nₗ, qᵣ, Nᵣ, rₑ, Z) in zip(q_liq, N_liq, q_rai, N_rai, reff, rr)
+            for SB in [SB2006, SB2006_no_limiters]
 
-        for SB in [SB2006, SB2006_no_limiters]
-            #action
-            reff = CM2.effective_radius(SB, q_liq, q_rai, N_liq, N_rai, ρ_air)
+                #action
+                Z_val = CM2.radar_reflectivity(SB, qₗ, qᵣ, Nₗ, Nᵣ, ρₐ)
+                rₑ_val = CM2.effective_radius(SB, qₗ, qᵣ, Nₗ, Nᵣ, ρₐ)
 
-            #test
-            TT.@test reff ≈ FT(2.66e-05) atol = FT(8e-6)
+                #test
+                TT.@test rₑ_val ≈ rₑ atol = FT(1e-6)
+                TT.@test Z_val ≈ Z atol = FT(1e-4)
+            end
         end
     end
 
@@ -624,7 +615,7 @@ function test_microphysics2M(FT)
         TT.@test qD ≈ qᵣ atol = FT(1e-6)
         TT.@test qx ≈ qᵣ atol = FT(1e-6)
 
-        # The mass integrals don't work with limiters 
+        # The mass integrals don't work with limiters
         TT.@test qx_lim ≈ qᵣ atol = FT(1e-6)
         TT.@test qD_lim ≈ qx_lim atol = FT(1e-6)
     end
@@ -665,7 +656,7 @@ function test_microphysics2M(FT)
         TT.@test qx ≈ qₗ rtol = FT(1e-6)
         TT.@test Nx ≈ Nₗ rtol = FT(1e-6)
 
-        # mass of liquid droplets as a function of its diameter in mm and μg 
+        # mass of liquid droplets as a function of its diameter in mm and μg
         _m(D) = FT(π / 6) * ρₗ * FT(10)^(χ - 9) * D^3
 
         # integral bounds in millimiters