diff --git a/ext/CreateParametersExt.jl b/ext/CreateParametersExt.jl index 1e7a4935e3..1a460ff607 100644 --- a/ext/CreateParametersExt.jl +++ b/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,6 +550,7 @@ end α_NIR_leaf = 0.4, τ_NIR_leaf = 0.25, Ω = 1, + LS_c = 0.05, n_layers = UInt64(20), kwargs... ) @@ -556,7 +558,11 @@ end 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,6 +605,7 @@ end α_PAR_leaf = 0.1, α_NIR_leaf = 0.4, Ω = 1, + LS_c = 0.05, kwargs... ) function BeerLambertParameters(toml_dict; @@ -604,13 +613,18 @@ end α_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..., ) diff --git a/src/standalone/Vegetation/Canopy.jl b/src/standalone/Vegetation/Canopy.jl index 1e4baf9194..4f5aff3456 100644 --- a/src/standalone/Vegetation/Canopy.jl +++ b/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] diff --git a/src/standalone/Vegetation/canopy_parameterizations.jl b/src/standalone/Vegetation/canopy_parameterizations.jl index ed44a8bf4e..9b319610fd 100644 --- a/src/standalone/Vegetation/canopy_parameterizations.jl +++ b/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,6 +433,9 @@ 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}}, @@ -438,7 +443,7 @@ function extinction_coeff( ) 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 diff --git a/src/standalone/Vegetation/radiation.jl b/src/standalone/Vegetation/radiation.jl index 0eeb926a7e..1d013b831a 100644 --- a/src/standalone/Vegetation/radiation.jl +++ b/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 diff --git a/test/standalone/Vegetation/radiation_roundoff.jl b/test/standalone/Vegetation/radiation_roundoff.jl new file mode 100644 index 0000000000..3abf866df5 --- /dev/null +++ b/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