Merge pull request #407 from CliMA/cc/1MomentAdjustment

1M radarrefl adjustment

src/Microphysics1M.jl

@@ -91,7 +91,7 @@ function lambda(
     # mass(size)
     (; r0, m0, me, Δm, χm) = mass
 
-    return q > FT(0) ?
+    return q > eps(FT) ?
            (
         χm * m0 * n0 * SF.gamma(me + Δm + FT(1)) / ρ / q / r0^(me + Δm)
     )^FT(1 / (me + Δm + 1)) : FT(0)
@@ -113,11 +113,15 @@ function radar_reflectivity(
     ρ::FT,
 ) where {FT}
 
-    n0 = get_n0(pdf)
-    λ = lambda(pdf, mass, q, ρ)
-    Z₀ = FT(1e-18)
+    # change units for accuracy
+    n0 = get_n0(pdf) * FT(1e-12)
+    λ = lambda(pdf, mass, q, ρ) * FT(1e-3)
 
-    return 10 * log10((720 * n0 / λ^7) / Z₀)
+    Z = (λ == FT(0)) ? FT(0) : (720 * n0 / λ^7)
+    log_10_Z₀ = FT(-18)
+    log_Z = FT(10) * (log10(Z) - log_10_Z₀ - FT(9))
+
+    return max(FT(-430), log_Z)
 end
 
 """
@@ -140,7 +144,7 @@ function terminal_velocity(
     ρ::FT,
     q::FT,
 ) where {FT}
-    if q > FT(0)
+    if q > eps(FT)
         # terminal_velocity(size)
         (; χv, ve, Δv) = vel
         v0 = get_v0(vel, ρ)
@@ -164,7 +168,7 @@ function terminal_velocity(
     q::FT,
 ) where {FT}
     fall_w = FT(0)
-    if q > FT(0)
+    if q > eps(FT)
         # coefficients from Table B1 from Chen et. al. 2022
         aiu, bi, ciu = CO.Chen2022_vel_coeffs_small(vel, ρ)
         # size distribution parameter
@@ -183,7 +187,7 @@ function terminal_velocity(
     q::FT,
 ) where {FT}
     fall_w = FT(0)
-    if q > FT(0)
+    if q > eps(FT)
 
         (; r0, m0, me, Δm, χm) = mass
         (; a0, ae, Δa, χa) = area
@@ -322,7 +326,7 @@ function accretion(
 ) where {FT}
 
     accr_rate = FT(0)
-    if (q_clo > FT(0) && q_pre > FT(0))
+    if (q_clo > eps(FT) && q_pre > eps(FT))
 
         n0::FT = get_n0(precip.pdf, q_pre, ρ)
         v0::FT = get_v0(vel, ρ)
@@ -479,7 +483,7 @@ function evaporation_sublimation(
     evap_subl_rate = FT(0)
     S = TD.supersaturation(tps, q, ρ, T, TD.Liquid())
 
-    if (q_rai > FT(0) && S < FT(0))
+    if (q_rai > eps(FT) && S < FT(0))
 
         (; ν_air, D_vapor) = aps
         G = CO.G_func(aps, tps, T, TD.Liquid())
@@ -515,7 +519,7 @@ function evaporation_sublimation(
     T::FT,
 ) where {FT}
     evap_subl_rate = FT(0)
-    if q_sno > FT(0)
+    if q_sno > eps(FT)
         (; ν_air, D_vapor) = aps
 
         S = TD.supersaturation(tps, q, ρ, T, TD.Ice())

test/microphysics1M_tests.jl

@@ -45,10 +45,28 @@ function test_microphysics1M(FT)
         ρ_air, q_tot, ρ_air_ground = FT(1.2), FT(20 * 1e-3), FT(1.22)
 
         for q_rai in q_rain_range
-            TT.@test CM1.terminal_velocity(rain, blk1mvel.rain, ρ_air, q_rai) ≈
-                     terminal_velocity_empir(q_rai, q_tot, ρ_air, ρ_air_ground) atol =
-                0.2 * terminal_velocity_empir(q_rai, q_tot, ρ_air, ρ_air_ground)
-
+            if q_rai > eps(FT)
+                TT.@test CM1.terminal_velocity(
+                    rain,
+                    blk1mvel.rain,
+                    ρ_air,
+                    q_rai,
+                ) ≈ terminal_velocity_empir(
+                    q_rai,
+                    q_tot,
+                    ρ_air,
+                    ρ_air_ground,
+                ) atol =
+                    0.2 *
+                    terminal_velocity_empir(q_rai, q_tot, ρ_air, ρ_air_ground)
+            else
+                TT.@test CM1.terminal_velocity(
+                    rain,
+                    blk1mvel.rain,
+                    ρ_air,
+                    q_rai,
+                ) ≈ FT(0) atol = FT(0.2)
+            end
         end
     end
 
@@ -228,16 +246,28 @@ function test_microphysics1M(FT)
         ρ_air, q_liq, q_tot = FT(1.2), FT(5e-4), FT(20e-3)
 
         for q_rai in q_rain_range
-            TT.@test CM1.accretion(
-                liquid,
-                rain,
-                blk1mvel.rain,
-                ce,
-                q_liq,
-                q_rai,
-                ρ_air,
-            ) ≈ accretion_empir(q_rai, q_liq, q_tot) atol =
-                (0.1 * accretion_empir(q_rai, q_liq, q_tot))
+            if q_rai > eps(FT)
+                TT.@test CM1.accretion(
+                    liquid,
+                    rain,
+                    blk1mvel.rain,
+                    ce,
+                    q_liq,
+                    q_rai,
+                    ρ_air,
+                ) ≈ accretion_empir(q_rai, q_liq, q_tot) atol =
+                    (0.1 * accretion_empir(q_rai, q_liq, q_tot))
+            else
+                TT.@test CM1.accretion(
+                    liquid,
+                    rain,
+                    blk1mvel.rain,
+                    ce,
+                    q_liq,
+                    q_rai,
+                    ρ_air,
+                ) ≈ FT(0) atol = FT(0.2)
+            end
         end
     end
 
@@ -352,17 +382,30 @@ function test_microphysics1M(FT)
         ρ = p / R / T
 
         for q_rai in q_rain_range
-            TT.@test CM1.evaporation_sublimation(
-                rain,
-                blk1mvel.rain,
-                aps,
-                tps,
-                q,
-                q_rai,
-                ρ,
-                T,
-            ) ≈ rain_evap_empir(tps, q_rai, q, T, p, ρ) atol =
-                -0.5 * rain_evap_empir(tps, q_rai, q, T, p, ρ)
+            if q_rai > eps(FT)
+                TT.@test CM1.evaporation_sublimation(
+                    rain,
+                    blk1mvel.rain,
+                    aps,
+                    tps,
+                    q,
+                    q_rai,
+                    ρ,
+                    T,
+                ) ≈ rain_evap_empir(tps, q_rai, q, T, p, ρ) atol =
+                    -0.5 * rain_evap_empir(tps, q_rai, q, T, p, ρ)
+            else
+                TT.@test CM1.evaporation_sublimation(
+                    rain,
+                    blk1mvel.rain,
+                    aps,
+                    tps,
+                    q,
+                    q_rai,
+                    ρ,
+                    T,
+                ) ≈ FT(0) atol = FT(0.2)
+            end
         end
 
         # no condensational growth for rain