From c901b71d3903a0029ebe4917621281e167d05213 Mon Sep 17 00:00:00 2001
From: Sajjad Azimi <azimi@caltech.edu>
Date: Thu, 4 Jul 2024 11:44:56 -0700
Subject: [PATCH] bugfix; evaporation should return zero (instead of NaN) when
 xr = 0

---
 NEWS.md                      |  1 +
 Project.toml                 |  2 +-
 src/Microphysics2M.jl        |  5 ++++-
 test/microphysics2M_tests.jl | 13 +++++++++++++
 4 files changed, 19 insertions(+), 2 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index f4d2ffe93..1311d1ebc 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,6 +5,7 @@ main
 ------
 <!--- # Add changes since the most recent release here --->
 
+- Bugfix; fixed the evaporation scheme of SB2006 to prevent returning NaN values when limiters are not applied. ([#420](https://github.com/CliMA/CloudMicrophysics.jl/pull/415))
 
 v0.20.0
 ------
diff --git a/Project.toml b/Project.toml
index 839a6979d..c2b0aa48d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "CloudMicrophysics"
 uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
 authors = ["Climate Modeling Alliance"]
-version = "0.22.0"
+version = "0.22.1"
 
 [deps]
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
diff --git a/src/Microphysics2M.jl b/src/Microphysics2M.jl
index a1aa5f187..4ef299097 100644
--- a/src/Microphysics2M.jl
+++ b/src/Microphysics2M.jl
@@ -560,7 +560,10 @@ 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
+        # When xr = 0 evap_rate_0 becomes NaN. We replace NaN with 0 which is the limit of 
+        # evap_rate_0 for xr -> 0.
+        evap_rate_0 =
+            N_rai < eps(FT) || xr / x_star < eps(FT) ? FT(0) : evap_rate_0
         evap_rate_1 = q_rai < eps(FT) ? FT(0) : evap_rate_1
     end
 
diff --git a/test/microphysics2M_tests.jl b/test/microphysics2M_tests.jl
index 40df43d42..ad59bc214 100644
--- a/test/microphysics2M_tests.jl
+++ b/test/microphysics2M_tests.jl
@@ -553,6 +553,19 @@ function test_microphysics2M(FT)
                 T,
             ).evap_rate_1 ≈ 0 atol = eps(FT)
         end
+
+        # test limit case: xr = 0 for SB with no limiters
+        TT.@test CM2.rain_evaporation(
+            SB2006_no_limiters,
+            aps,
+            tps,
+            q,
+            FT(0),
+            ρ,
+            N_rai,
+            T,
+        ).evap_rate_0 ≈ 0 atol = eps(FT)
+
     end
 
     TT.@testset "2M_microphysics - Seifert and Beheng 2006 effective radius and reflectivity" begin