absorbed radiation is not zero if LAI is zero

ext/CreateParametersExt.jl

@@ -539,6 +539,7 @@ end
         α_NIR_leaf = 0.4,
         τ_NIR_leaf = 0.25,
         Ω = 1,
+        LS_c = 0.05,
         n_layers = UInt64(20),
         kwargs...
     )
@@ -549,14 +550,19 @@ end
         α_NIR_leaf = 0.4,
         τ_NIR_leaf = 0.25,
         Ω = 1,
+        LS_c = 0.05,
         n_layers = UInt64(20),
         kwargs...
     )
 
 Floating-point and toml dict based constructor supplying default values
 for the TwoStreamParameters struct. Additional parameter values can be directly set via kwargs.
 """
-TwoStreamParameters(::Type{FT}; kwargs...) where {FT <: AbstractFloat} =
+TwoStreamParameters(
+    ::Type{FT};
+    LS_c = FT(0.05),
+    kwargs...,
+) where {FT <: AbstractFloat} =
     TwoStreamParameters(CP.create_toml_dict(FT); kwargs...)
 
 function TwoStreamParameters(
@@ -567,6 +573,7 @@ function TwoStreamParameters(
     α_NIR_leaf = 0.4,
     τ_NIR_leaf = 0.25,
     Ω = 1,
+    LS_c = 0.05,
     n_layers = UInt64(20),
     kwargs...,
 )
@@ -585,6 +592,7 @@ function TwoStreamParameters(
         α_NIR_leaf,
         τ_NIR_leaf,
         Ω,
+        LS_c,
         n_layers,
         parameters...,
         kwargs...,
@@ -597,20 +605,26 @@ end
         α_PAR_leaf = 0.1,
         α_NIR_leaf = 0.4,
         Ω = 1,
+        LS_c = 0.05,
         kwargs...
     )
     function BeerLambertParameters(toml_dict;
         ld = (_) -> 0.5,
         α_PAR_leaf = 0.1,
         α_NIR_leaf = 0.4,
         Ω = 1,
+        LS_c = 0.05,
         kwargs...
     )
 
 Floating-point and toml dict based constructor supplying default values
 for the BeerLambertParameters struct. Additional parameter values can be directly set via kwargs.
 """
-BeerLambertParameters(::Type{FT}; kwargs...) where {FT <: AbstractFloat} =
+BeerLambertParameters(
+    ::Type{FT};
+    LS_c = FT(0.05),
+    kwargs...,
+) where {FT <: AbstractFloat} =
     BeerLambertParameters(CP.create_toml_dict(FT); kwargs...)
 
 function BeerLambertParameters(
@@ -619,6 +633,7 @@ function BeerLambertParameters(
     α_PAR_leaf = 0.1,
     α_NIR_leaf = 0.4,
     Ω = 1,
+    LS_c = 0.05,
     kwargs...,
 )
     name_map = (;
@@ -634,6 +649,7 @@ function BeerLambertParameters(
         α_PAR_leaf,
         α_NIR_leaf,
         Ω,
+        LS_c,
         parameters...,
         kwargs...,
     )

src/standalone/Vegetation/Canopy.jl

@@ -435,7 +435,7 @@ function ClimaLand.make_update_aux(
         ρ_l = FT(LP.ρ_cloud_liq(earth_param_set))
         R = FT(LP.gas_constant(earth_param_set))
         thermo_params = earth_param_set.thermo_params
-        (; G_Function, Ω, λ_γ_PAR, λ_γ_NIR) =
+        (; G_Function, Ω, λ_γ_PAR, λ_γ_NIR, LS_c) =
             canopy.radiative_transfer.parameters
         energy_per_photon_PAR = planck_h * c / λ_γ_PAR
         energy_per_photon_NIR = planck_h * c / λ_γ_NIR
@@ -563,7 +563,7 @@ function ClimaLand.make_update_aux(
             c_co2_air,
             R,
         )
-        @. GPP = compute_GPP(An, K, LAI, Ω)
+        @. GPP = compute_GPP(An, K, LAI, Ω, LS_c)
         @. gs = medlyn_conductance(g0, Drel, medlyn_factor, An, c_co2_air)
         # update autotrophic respiration
         h_canopy = hydraulics.compartment_surfaces[end]

src/standalone/Vegetation/canopy_parameterizations.jl

@@ -188,12 +188,15 @@ function plant_absorbed_pfd(
     α_soil::FT,
 ) where {FT}
     RTP = RT.parameters
-    AR = SW_IN * (1 - α_leaf) * (1 - exp(-K * LAI * RTP.Ω)) * (1 - α_soil)
-    TR = SW_IN * exp(-K * LAI * RTP.Ω)
+    KLAIΩ = LAI > RTP.LS_c ? K * LAI * RTP.Ω : FT(0)
+    AR = SW_IN * (1 - α_leaf) * (1 - exp(-KLAIΩ)) * (1 - α_soil)
+    TR = SW_IN * exp(-KLAIΩ)
     RR = SW_IN - AR - TR * (1 - α_soil)
     return (; abs = AR, refl = RR, trans = TR)
+
 end
 
+
 """
     plant_absorbed_pfd(
         RT::TwoStreamModel{FT},
@@ -229,8 +232,8 @@ function plant_absorbed_pfd(
     α_soil::FT,
     frac_diff::FT,
 ) where {FT}
-
-    (; G_Function, Ω, n_layers) = RT.parameters
+    RTP = RT.parameters
+    (; G_Function, Ω, n_layers) = RTP
 
     # Get the current leaf angular distribution value based on the solar zenith angle
     G = compute_G(G_Function, θs)
@@ -241,6 +244,7 @@ function plant_absorbed_pfd(
     ω = α_leaf + τ_leaf
 
     # Compute aₛ, the single scattering albedo
+    # if θs = π/2 this seems OK
     aₛ = 0.5 * ω * (1 - cos(θs) * log((abs(cos(θs)) + 1) / abs(cos(θs))))
 
     # Compute β₀, the direct upscattering parameter
@@ -264,8 +268,10 @@ function plant_absorbed_pfd(
     u₂ = b - c * α_soil
     u₃ = f + c * α_soil
 
-    s₁ = exp(-h * LAI * Ω)
-    s₂ = exp(-K * LAI * Ω)
+    hLAIΩ = LAI > RTP.LS_c ? h * LAI * Ω : FT(0)
+    KLAIΩ = LAI > RTP.LS_c ? K * LAI * Ω : FT(0)
+    s₁ = exp(-hLAIΩ)
+    s₂ = exp(-KLAIΩ)
 
     p₁ = b + μ̄ * h
     p₂ = b - μ̄ * h
@@ -306,7 +312,7 @@ function plant_absorbed_pfd(
     h₁₀ = -s₁ * (u₂ - μ̄ * h) / d₂
 
     # Compute the LAI per layer for this canopy
-    Lₗ = LAI / n_layers
+    LAI_per_layer = LAI / n_layers
 
     # Initialize the fraction absorbed value and layer counter
     F_abs = 0
@@ -321,28 +327,24 @@ function plant_absorbed_pfd(
     I_dif_up_prev = 0
     I_dif_dn_prev = 0
 
+
     # Compute F_abs in each canopy layer
     while i <= n_layers
 
         # Compute cumulative LAI at this layer
-        L = i * Lₗ
-
+        L = i * LAI_per_layer
+        KΩL = LAI > RTP.LS_c ? K * L * Ω : FT(0)
+        hΩL = LAI > RTP.LS_c ? h * L * Ω : FT(0)
         # Compute the direct fluxes into/out of the layer 
-        I_dir_up =
-            h₁ * exp(-K * L * Ω) / σ +
-            h₂ * exp(-h * L * Ω) +
-            h₃ * exp(h * L * Ω)
-        I_dir_dn =
-            h₄ * exp(-K * L * Ω) / σ +
-            h₅ * exp(-h * L * Ω) +
-            h₆ * exp(h * L * Ω)
+        I_dir_up = h₁ * exp(-KΩL) / σ + h₂ * exp(-hΩL) + h₃ * exp(hΩL)
+        I_dir_dn = h₄ * exp(-KΩL) / σ + h₅ * exp(-hΩL) + h₆ * exp(hΩL)
 
         # Add collimated radiation to downard flux
-        I_dir_dn += exp(-K * L * Ω)
+        I_dir_dn += exp(-KΩL)
 
         # Compute the diffuse fluxes into/out of the layer
-        I_dif_up = h₇ * exp(-h * L * Ω) + h₈ * exp(h * L * Ω)
-        I_dif_dn = h₉ * exp(-h * L * Ω) + h₁₀ * exp(h * L * Ω)
+        I_dif_up = h₇ * exp(-hΩL) + h₈ * exp(hΩL)
+        I_dif_dn = h₉ * exp(-hΩL) + h₁₀ * exp(hΩL)
 
         # Energy balance giving radiation absorbed in the layer 
         if i == 0
@@ -370,7 +372,6 @@ function plant_absorbed_pfd(
         # Move on to the next layer
         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)
@@ -379,6 +380,7 @@ function plant_absorbed_pfd(
         refl = FT(SW_IN * F_refl),
         trans = FT(SW_IN * F_trans),
     )
+
 end
 
 """
@@ -431,14 +433,17 @@ end
 
 Computes the vegetation extinction coefficient (`K`), as a function
 of the sun zenith angle (`θs`), and the leaf angle distribution (`ld`).
+
+In order to prevent errors as zenith angle → π/2, we clip K
+to a maximum value of 100.
 """
 function extinction_coeff(
     G::Union{ConstantGFunction{FT}, CLMGFunction{FT}},
     θs::FT,
 ) where {FT}
     ld = FT(compute_G(G, θs))
     K = ld / max(cos(θs), eps(FT))
-    return K
+    return min(K, FT(1e2))
 end
 
 # 2. Photosynthesis, Farquhar model
@@ -738,7 +743,8 @@ Computes the total canopy photosynthesis (`GPP`) as a function of
 the total net carbon assimilation (`An`), the extinction coefficient (`K`),
 leaf area index (`LAI`) and the clumping index (`Ω`).
 """
-function compute_GPP(An::FT, K::FT, LAI::FT, Ω::FT) where {FT}
+function compute_GPP(An::FT, K::FT, LAI::FT, Ω::FT, LS_c::FT) where {FT}
+    KLAIΩ = K * LAI * Ω#LAI > LS_c ? K * LAI * Ω : FT(0)
     GPP = An * (1 - exp(-K * LAI * Ω)) / (K * Ω)
     return GPP
 end

src/standalone/Vegetation/radiation.jl

@@ -64,6 +64,8 @@ Base.@kwdef struct BeerLambertParameters{
     λ_γ_NIR::FT
     "Leaf angle distribution function"
     G_Function::G
+    "Critical value of LAI+SAI below which land is treated as unvegetated"
+    LS_c::FT
 end
 
 Base.eltype(::BeerLambertParameters{FT}) where {FT} = FT
@@ -112,23 +114,9 @@ Base.@kwdef struct TwoStreamParameters{
     n_layers::UInt64
     "Leaf angle distribution function"
     G_Function::G
+    "Critical value of LAI+SAI below which land is treated as unvegetated"
+    LS_c::FT
 end
-"""
-    function TwoStreamParameters{FT, G}(;
-        ld = ConstantGFunction(FT(0.5)),
-        α_PAR_leaf = FT(0.3),
-        τ_PAR_leaf = FT(0.2),
-        α_NIR_leaf = FT(0.4),
-        τ_NIR_leaf = FT(0.25),
-        ϵ_canopy = FT(0.98),
-        Ω = FT(1),
-        λ_γ_PAR = FT(5e-7),
-        λ_γ_NIR = FT(1.65e-6),
-        n_layers = UInt64(20),
-    ) where {FT}
-
-A constructor supplying default values for the TwoStreamParameters struct.
-"""
 
 Base.eltype(::TwoStreamParameters{FT}) where {FT} = FT
 

test/standalone/Vegetation/radiation_roundoff.jl

@@ -0,0 +1,89 @@
+using Test
+using ClimaLand
+
+
+for FT in (Float32, Float64)
+    @testset "Roundoff bug, FT = $FT" begin
+        Ω = FT(0.69)
+        χl = FT(0.1)
+        G_Function = ClimaLand.Canopy.CLMGFunction(χl)
+        α_PAR_leaf = FT(0.1)
+        λ_γ_PAR = FT(5e-7)
+        λ_γ_NIR = FT(1.65e-6)
+        τ_PAR_leaf = FT(0.05)
+        α_NIR_leaf = FT(0.45)
+        τ_NIR_leaf = FT(0.25)
+        ϵ_canopy = FT(0.97)
+        ts_params = ClimaLand.Canopy.TwoStreamParameters(
+            α_PAR_leaf,
+            τ_PAR_leaf,
+            α_NIR_leaf,
+            τ_NIR_leaf,
+            ϵ_canopy,
+            Ω,
+            λ_γ_PAR,
+            λ_γ_NIR,
+            UInt64(20),
+            G_Function,
+            FT(0.05),
+        )
+
+        BL_params = ClimaLand.Canopy.BeerLambertParameters(
+            α_PAR_leaf,
+            α_NIR_leaf,
+            ϵ_canopy,
+            Ω,
+            λ_γ_PAR,
+            λ_γ_NIR,
+            G_Function,
+            FT(0.05),
+        )
+
+        RTmodels = [
+            ClimaLand.Canopy.TwoStreamModel{FT, typeof(ts_params)}(ts_params),
+            ClimaLand.Canopy.BeerLambertModel{FT, typeof(BL_params)}(BL_params),
+        ]
+
+        PAR = FT(1.0)
+
+        LAI = eps(FT)
+        θs = FT(π / 2)
+        K = ClimaLand.Canopy.extinction_coeff(G_Function, θs)
+        α_soil_PAR = FT(0.2)
+        frac_diff = FT(0.0)
+
+        RT = RTmodels[1]
+        canopy_par_ts = @. ClimaLand.Canopy.plant_absorbed_pfd(
+            RT,
+            PAR,
+            RT.parameters.α_PAR_leaf,
+            RT.parameters.τ_PAR_leaf,
+            LAI,
+            K,
+            θs,
+            α_soil_PAR,
+            frac_diff,
+        )[1]
+        # INC - OUT = CANOPY_ABS + (1-α_soil)*CANOPY_TRANS
+        # OUT = REFL
+        @test canopy_par_ts.abs < eps(FT)
+        @test canopy_par_ts.refl ≈ FT(α_soil_PAR)
+        @test canopy_par_ts.trans ≈ FT(1)
+
+        RT = RTmodels[2]
+        canopy_par_bl = @. ClimaLand.Canopy.plant_absorbed_pfd(
+            RT,
+            PAR,
+            RT.parameters.α_PAR_leaf,
+            LAI,
+            K,
+            α_soil_PAR,
+        )[1]
+        # INC - OUT = CANOPY_ABS + (1-α_soil)*CANOPY_TRANS
+        # OUT = REFL
+        @test canopy_par_bl.abs < eps(FT)
+        @test canopy_par_bl.refl ≈ FT(α_soil_PAR)
+        @test canopy_par_bl.trans ≈ FT(1)
+
+    end
+end