From dd921b77cbabfb2c985212b3b83b780c264d2ac1 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Thu, 24 Oct 2024 17:19:19 +0200
Subject: [PATCH] fix open water evaporation

---
 src/sbm.jl | 31 +++++++++++++++++++++++--------
 1 file changed, 23 insertions(+), 8 deletions(-)

diff --git a/src/sbm.jl b/src/sbm.jl
index 432408e40..3c32b9500 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -909,7 +909,18 @@ function update_until_recharge(sbm::SBM, config)
         ustoredepth = sum(@view sbm.ustorelayerdepth[i][1:sbm.nlayers[i]])
 
         runoff_river = min(1.0, sbm.riverfrac[i]) * avail_forinfilt
-        runoff_land = min(1.0, sbm.waterfrac[i]) * avail_forinfilt
+
+        # Calculate the initial capacity of the unsaturated store
+        ustorecapacity = sbm.soilwatercapacity[i] - sbm.satwaterdepth[i] - ustoredepth
+        # Add a flag if the soil is saturated
+        check_saturated = ustorecapacity < 0.1
+
+        if check_saturated
+            runoff_land = min(1.0, 1.0 - sbm.riverfrac[i]) * avail_forinfilt
+        else
+            runoff_land = min(1.0, sbm.waterfrac[i]) * avail_forinfilt
+        end
+
         if !isnothing(sbm.paddy) || !isnothing(sbm.nonpaddy)
             avail_forinfilt = avail_forinfilt + sbm.allocation.irri_alloc[i]
         end
@@ -928,10 +939,17 @@ function update_until_recharge(sbm::SBM, config)
         surface_water = sbm.waterlevel_land[i] * (1.0 - sbm.riverfrac[i])
 
         # Calculate open water evaporation
-        ae_openw_l = min(
-            surface_water,
-            sbm.potential_evaporation[i] * soilevap_fraction,
-        )
+        if check_saturated
+            ae_openw_l = min(
+                surface_water,
+                sbm.potential_evaporation[i] * soilevap_fraction,
+            )
+        else
+            ae_openw_l = min(
+                surface_water,
+                sbm.potential_evaporation[i] * sbm.waterfrac[i],
+            )
+        end
 
         # Update land waterlevel and scale to part of cell not covered by rivers
         if do_surface_water_infiltration
@@ -954,9 +972,6 @@ function update_until_recharge(sbm::SBM, config)
             evap_paddy_water = 0.0
         end
 
-        # Calculate the initial capacity of the unsaturated store
-        ustorecapacity = sbm.soilwatercapacity[i] - sbm.satwaterdepth[i] - ustoredepth
-
         # Calculate the infiltration flux into the soil column
         infiltsoilpath,
         infiltsoil,