Fix error in reflected radiation; computed transmitted fraction differently (#1156)

AlexisRenchon · web-flow · commit f4b2405d134e · 2025-05-15T19:34:18.000Z
diff --git a/src/standalone/Vegetation/canopy_parameterizations.jl b/src/standalone/Vegetation/canopy_parameterizations.jl
@@ -236,15 +236,17 @@ function canopy_sw_rt_two_stream(
     α_soil::FT,
     frac_diff::FT,
 ) where {FT}
-
+    α_soil = max(eps(FT), α_soil) # this prevents division by zero, below.
+    cosθs = max(eps(FT), cosθs) # The insolations package returns θs > π/2 (nighttime), but this assumes cosθs >0
     # Compute μ̄, the average inverse diffuse optical length per LAI
-    μ̄ = 1 / (2G)
+    μ̄ = 1 / max(2G, eps(FT))
 
-    # Clip this to eps(FT) to prevent dividing by zero
-    ω = max(α_leaf + τ_leaf, 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))
 
     # Compute aₛ, the single scattering albedo
-    aₛ = 0.5 * ω * (1 - cosθs * log((abs(cosθs) + 1) / abs(cosθs)))
+    aₛ =
+        0.5 * ω * (1 - cosθs * log((abs(cosθs) + 1) / max(abs(cosθs), eps(FT))))
 
     # Compute β₀, the direct upscattering parameter
     β₀ = (1 / ω) * aₛ * (1 + μ̄ * K) / (μ̄ * K)
@@ -315,9 +317,12 @@ function canopy_sw_rt_two_stream(
     F_abs = 0
     i = 0
 
-    # Total light reflected form top of canopy
+    # Total light reflected from top of canopy
     F_refl = 0
 
+    # Total light transmitted through the canopy on the downward pass
+    F_trans = 0
+
     # Intialize vars to save computed fluxes from each layer for the next layer
     I_dir_up_prev = 0
     I_dir_dn_prev = 0
@@ -341,7 +346,7 @@ function canopy_sw_rt_two_stream(
             h₅ * exp(-h * L * Ω) +
             h₆ * exp(h * L * Ω)
 
-        # Add collimated radiation to downard flux
+        # Add collimated radiation to downward flux
         I_dir_dn += exp(-K * L * Ω)
 
         # Compute the diffuse fluxes into/out of the layer
@@ -357,8 +362,12 @@ function canopy_sw_rt_two_stream(
             I_dif_abs = I_dif_up - I_dif_up_prev - I_dif_dn + I_dif_dn_prev
         end
 
-        if i == 1
-            F_refl = (1 - frac_diff) * I_dir_up + (frac_diff) * I_dif_up
+        if i == 1 # prev = top of layer
+            F_refl =
+                (1 - frac_diff) * I_dir_up_prev + (frac_diff) * I_dif_up_prev
+        end
+        if i == n_layers # not prev = bottom of layer
+            F_trans = (1 - frac_diff) * I_dir_dn + (frac_diff) * I_dif_dn
         end
 
 
@@ -375,9 +384,13 @@ function canopy_sw_rt_two_stream(
         i += 1
     end
 
-    # Convert fractional absorption into absorption and return
-    # Ensure floating point precision is correct (it may be different for PAR)
-    F_trans = (1 - F_abs - F_refl) / (1 - α_soil)
+    # Reflected is reflected from the land surface. F_abs
+    # refers to radiation absorbed by the canopy on either pass.
+    # F_trans refers to the downward pass only. Note that
+    # (1-α_soil)*F_trans + F_abs = total absorbed fraction by land, so the following
+    # must hold
+    # @assert (1 - α_soil) * FT(F_trans) + FT(F_abs) + FT(F_refl) ≈ 1
+    # This is tested in test/standalone/Vegetation/test_two_stream.jl
     return (; abs = FT(F_abs), refl = FT(F_refl), trans = FT(F_trans))
 end
 
@@ -388,9 +401,12 @@ end
 Computes the vegetation extinction coefficient (`K`), as a function
 of the cosine of the sun zenith angle (`cosθs`),
 and the leaf angle distribution (`G`).
+
+In the two-stream scheme, values of K ~ 1/epsilon can lead to numerical issues.
+Here we clip it to 1e6.
 """
 function extinction_coeff(G::FT, cosθs::FT) where {FT}
-    K = G / max(cosθs, eps(FT))
+    K = min(G / max(cosθs, eps(FT)), FT(1e6))
     return K
 end
 
diff --git a/test/integrated/soil_canopy_lsm.jl b/test/integrated/soil_canopy_lsm.jl
@@ -6,6 +6,7 @@ using ClimaLand
 using ClimaLand.Soil
 using ClimaLand.Canopy
 using Dates
+using ClimaParams
 import ClimaLand.Parameters as LP
 using ClimaLand.Soil.Biogeochemistry
 using ClimaLand.Canopy.PlantHydraulics
diff --git a/test/standalone/Vegetation/test_two_stream.jl b/test/standalone/Vegetation/test_two_stream.jl
@@ -1,109 +1,153 @@
-# This script tests the output of the ClimaLand TwoStream implementation
+# This script tests the canopy absorbed output of the ClimaLand TwoStream implementation
 # against the T. Quaife pySellersTwoStream implementation by providing
 # the same setups to each model and checking that the outputs are equal.
 # The output of the python module for a variety of inputs is given in the
 # linked file which then gets read and compared to the Clima version, ensuring
 # the FAPAR of each model is within 1/2 of a percentage point.
 
+# It also tests that the sum of absorbed and reflected radiation is 1.
+
 using Test
 using ClimaLand
 import ClimaComms
 ClimaComms.@import_required_backends
 using ClimaLand.Canopy
 using DelimitedFiles
 using ClimaLand.Domains: Point
+import ClimaLand.Parameters as LP
+import ClimaParams
 
-include("../../Artifacts.jl")
+@testset "Comparison to pySellersTwoStream" begin
+    include("../../Artifacts.jl")
 
-# Read the test data from the ClimaArtifact
-datapath = twostr_test_data_path()
+    # Read the test data from the ClimaArtifact
+    datapath = twostr_test_data_path()
 
-data = joinpath(datapath, "twostr_test.csv")
-test_set = readdlm(data, ',')
+    data = joinpath(datapath, "twostr_test.csv")
+    test_set = readdlm(data, ',')
 
-# Floating point precision to use
-for FT in (Float32, Float64)
-    @testset "Two-Stream Model Correctness, FT = $FT" begin
-        # Read the conditions for each setup parameter from the test file
-        column_names = test_set[1, :]
-        cosθs = FT.(test_set[2:end, column_names .== "mu"])
-        LAI = FT.(test_set[2:end, column_names .== "LAI"])
-        a_soil = FT.(test_set[2:end, column_names .== "a_soil"])
-        n_layers = UInt64.(test_set[2:end, column_names .== "n_layers"])
-        PropDif = FT.(test_set[2:end, column_names .== "prop_diffuse"])
+    # Floating point precision to use
+    for FT in (Float32, Float64)
+        @testset "Two-Stream Model Correctness, FT = $FT" begin
+            # Read the conditions for each setup parameter from the test file
+            column_names = test_set[1, :]
+            cosθs = FT.(test_set[2:end, column_names .== "mu"])
+            LAI = FT.(test_set[2:end, column_names .== "LAI"])
+            a_soil = FT.(test_set[2:end, column_names .== "a_soil"])
+            n_layers = UInt64.(test_set[2:end, column_names .== "n_layers"])
+            PropDif = FT.(test_set[2:end, column_names .== "prop_diffuse"])
 
-        # setup spatially varying params as both float and spatially varying
-        domain = Point(; z_sfc = FT(0.0))
-        lds = FT.(test_set[2:end, column_names .== "ld"])
-        lds_field = map(x -> fill(x, domain.space.surface), lds)
-        α_PAR_leaf_scalars = FT.(test_set[2:end, column_names .== "rho"])
-        α_PAR_leaf_fields =
-            map(x -> fill(x, domain.space.surface), α_PAR_leaf_scalars)
-        τ_scalars = FT.(test_set[2:end, column_names .== "tau"])
-        τ_fields = map(x -> fill(x, domain.space.surface), τ_scalars)
-        # loop through once with params as floats, then with params as fields
-        Ω_cases = (FT(1), fill(FT(1), domain.space.surface))
-        α_PAR_leaf_cases = (α_PAR_leaf_scalars, α_PAR_leaf_fields)
-        τ_PAR_leaf_cases = (τ_scalars, τ_fields)
-        α_NIR_leaf_cases = (FT(0.4), fill(FT(0.4), domain.space.surface))
-        τ_NIR_leaf_cases = (FT(0.25), fill(FT(0.24), domain.space.surface))
-        lds_cases = (lds, lds_field)
-        zipped_params = zip(
-            Ω_cases,
-            α_PAR_leaf_cases,
-            τ_PAR_leaf_cases,
-            α_NIR_leaf_cases,
-            τ_NIR_leaf_cases,
-            lds_cases,
-        )
-        for (Ω, α_PAR_leaf, τ_PAR_leaf, α_NIR_leaf, τ_NIR_leaf, lds) in
-            zipped_params
-            # Read the result for each setup from the Python output
-            py_FAPAR = FT.(test_set[2:end, column_names .== "FAPAR"])
+            # setup spatially varying params as both float and spatially varying
+            domain = Point(; z_sfc = FT(0.0))
+            lds = FT.(test_set[2:end, column_names .== "ld"])
+            lds_field = map(x -> fill(x, domain.space.surface), lds)
+            α_PAR_leaf_scalars = FT.(test_set[2:end, column_names .== "rho"])
+            α_PAR_leaf_fields =
+                map(x -> fill(x, domain.space.surface), α_PAR_leaf_scalars)
+            τ_scalars = FT.(test_set[2:end, column_names .== "tau"])
+            τ_fields = map(x -> fill(x, domain.space.surface), τ_scalars)
+            # loop through once with params as floats, then with params as fields
+            Ω_cases = (FT(1), fill(FT(1), domain.space.surface))
+            α_PAR_leaf_cases = (α_PAR_leaf_scalars, α_PAR_leaf_fields)
+            τ_PAR_leaf_cases = (τ_scalars, τ_fields)
+            α_NIR_leaf_cases = (FT(0.4), fill(FT(0.4), domain.space.surface))
+            τ_NIR_leaf_cases = (FT(0.25), fill(FT(0.24), domain.space.surface))
+            lds_cases = (lds, lds_field)
+            zipped_params = zip(
+                Ω_cases,
+                α_PAR_leaf_cases,
+                τ_PAR_leaf_cases,
+                α_NIR_leaf_cases,
+                τ_NIR_leaf_cases,
+                lds_cases,
+            )
+            for (Ω, α_PAR_leaf, τ_PAR_leaf, α_NIR_leaf, τ_NIR_leaf, lds) in
+                zipped_params
+                # Read the result for each setup from the Python output
+                py_FAPAR = FT.(test_set[2:end, column_names .== "FAPAR"])
 
-            # Python code does not use clumping index, and λ_γ does not impact FAPAR
-            # Test over all rows in the stored output from the Python module
-            for i in 2:(size(test_set, 1) - 1)
+                # Python code does not use clumping index, and λ_γ does not impact FAPAR
+                # Test over all rows in the stored output from the Python module
+                for i in 2:(size(test_set, 1) - 1)
 
-                # Set the parameters based on the setup read from the file
-                RT_params = TwoStreamParameters(
-                    FT;
-                    Ω = Ω,
-                    G_Function = ConstantGFunction(FT.(lds[i])),
-                    α_PAR_leaf = α_PAR_leaf[i],
-                    τ_PAR_leaf = τ_PAR_leaf[i],
-                    α_NIR_leaf = α_NIR_leaf,
-                    τ_NIR_leaf = τ_NIR_leaf,
-                    n_layers = n_layers[i],
-                )
+                    # Set the parameters based on the setup read from the file
+                    RT_params = TwoStreamParameters(
+                        FT;
+                        Ω = Ω,
+                        G_Function = ConstantGFunction(FT.(lds[i])),
+                        α_PAR_leaf = α_PAR_leaf[i],
+                        τ_PAR_leaf = τ_PAR_leaf[i],
+                        α_NIR_leaf = α_NIR_leaf,
+                        τ_NIR_leaf = τ_NIR_leaf,
+                        n_layers = n_layers[i],
+                    )
 
-                # Initialize the TwoStream model
-                RT = TwoStreamModel(RT_params)
+                    # Initialize the TwoStream model
+                    RT = TwoStreamModel(RT_params)
 
-                # Compute the predicted FAPAR using the ClimaLand TwoStream implementation
-                G = compute_G(RT_params.G_Function, cosθs)
-                K = extinction_coeff.(G, cosθs[i])
-                output =
-                    canopy_sw_rt_two_stream.(
-                        G,
-                        RT_params.Ω,
-                        RT_params.n_layers,
-                        RT_params.α_PAR_leaf,
-                        RT_params.τ_PAR_leaf,
-                        LAI[i],
-                        K,
-                        cosθs[i],
-                        a_soil[i],
-                        PropDif[i],
+                    # Compute the predicted FAPAR using the ClimaLand TwoStream implementation
+                    G = compute_G(RT_params.G_Function, cosθs)
+                    K = extinction_coeff.(G, cosθs[i])
+                    output =
+                        canopy_sw_rt_two_stream.(
+                            G,
+                            RT_params.Ω,
+                            RT_params.n_layers,
+                            RT_params.α_PAR_leaf,
+                            RT_params.τ_PAR_leaf,
+                            LAI[i],
+                            K,
+                            cosθs[i],
+                            a_soil[i],
+                            PropDif[i],
+                        )
+                    FAPAR = output.abs
+                    # Check that the predictions are app. equivalent to the Python model
+                    # Create a field of the expect value because isapprox cannot be broadcast
+                    # over a field of floats. The domain is a point, so it makes no difference
+                    # to the error when FAPAR is a float
+                    expected_output = fill(py_FAPAR[i], domain.space.surface)
+                    @test isapprox(
+                        0,
+                        sum(FAPAR .- expected_output),
+                        atol = 0.005,
                     )
-                FAPAR = output.abs
-                # Check that the predictions are app. equivalent to the Python model
-                # Create a field of the expect value because isapprox cannot be broadcast
-                # over a field of floats. The domain is a point, so it makes no difference
-                # to the error when FAPAR is a float
-                expected_output = fill(py_FAPAR[i], domain.space.surface)
-                @test isapprox(0, sum(FAPAR .- expected_output), atol = 0.005)
+                end
             end
         end
     end
 end
+
+@testset "Test physicality" begin
+    N = 100000
+    θs = [rand(N - 1) * 2π..., π / 2]
+    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]
+    G = 0.5
+    K = ClimaLand.Canopy.extinction_coeff.(G, cosθs)
+    frac_diff = rand(N)
+    n_layers = UInt64(20)
+    Ω = rand(N)
+    LAI = round.(rand(N))
+    output =
+        ClimaLand.Canopy.canopy_sw_rt_two_stream.(
+            G,
+            Ω,
+            n_layers,
+            α_leaf,
+            τ_leaf,
+            LAI,
+            K,
+            cosθs,
+            α_soil,
+            frac_diff,
+        )
+    expected = zeros(N)
+    for i in 1:N
+        expected[i] =
+            output[i].trans * (1 - α_soil[i]) + output[i].abs + output[i].refl
+    end
+    @assert all(expected .≈ 1)
+end
