Bug in two stream scheme (#1193)

Author: kmdeck

src/standalone/Vegetation/canopy_parameterizations.jl

@@ -235,8 +235,9 @@ function canopy_sw_rt_two_stream(
     # Compute μ̄, the average inverse diffuse optical length per LAI
     μ̄ = 1 / max(2G, eps(FT))
 
-    # Clip this to eps(FT) to prevent dividing by zero; the sum must also be < 1
-    ω = min(max(α_leaf + τ_leaf, eps(FT)), 1 - eps(FT))
+    # Clip this to prevent dividing by zero; the sum must also be < 1. Note that using eps(FT)
+    # as the threshold leads to numerical errors, so we use 1e-4
+    ω = min(max(α_leaf + τ_leaf, FT(1e-4)), 1 - FT(1e-4))
 
     # Compute aₛ, the single scattering albedo
     aₛ =

test/standalone/Vegetation/test_two_stream.jl

@@ -125,7 +125,8 @@ end
     cosθs = cos.(θs)
     α_leaf = [rand(N - 4)..., 0.0, 0.0, 1.0, 1.0]
     τ_leaf = (1.0 .- α_leaf) .* rand(N)
-    α_soil = [rand(N - 6)..., 0.0, 1.0, 0.2, 0.2, 0.2, 0.2]
+    τ_leaf[end - 2] = 0.0
+    α_soil = [rand(N - 6)..., 0.0, 1.0, 0.2, 0.99, 0.2, 0.2]
     G_Function = ConstantGFunction(0.5)
     K = ClimaLand.Canopy.extinction_coeff.(G_Function, cosθs)
     frac_diff = rand(N)