From 165e7ecdec3bd31735e13199c3157d44bac4dbd6 Mon Sep 17 00:00:00 2001 From: Olivier Bonte Date: Fri, 1 Aug 2025 10:03:14 +0200 Subject: [PATCH 1/3] Add improved type stability (conversion and roughness) --- src/surface_roughness.jl | 17 +++++++++-------- src/unit_conversions.jl | 16 ++++++++-------- 2 files changed, 17 insertions(+), 16 deletions(-) diff --git a/src/surface_roughness.jl b/src/surface_roughness.jl index 7422b03..7ac6c6f 100755 --- a/src/surface_roughness.jl +++ b/src/surface_roughness.jl @@ -146,25 +146,26 @@ rp = roughness_parameters(Roughness_wind_profile(),df;zh,zr=40,d=0.8*zh) ≈(rp.z0m, 0.55, rtol=0.1) # output true -``` -""" -function roughness_parameters(::RoughnessCanopyHeight, zh; frac_d=0.7, frac_z0m=0.1) +``` +""" +function roughness_parameters(::RoughnessCanopyHeight, zh::FT; + frac_d=FT(0.7), frac_z0m=FT(0.1)) where FT d = frac_d*zh z0m = frac_z0m*zh z0m_se = missing (;d, z0m, z0m_se) end -function roughness_parameters(::RoughnessCanopyHeightLAI, zh, LAI; - cd=0.2, hs=0.01) +function roughness_parameters(::RoughnessCanopyHeightLAI, zh::FT, LAI; + cd=FT(0.2), hs=FT(0.01)) where FT X = cd * LAI - d = 1.1 * zh * log(1 + X^(1/4)) - z0m = ifelse(0 <= X <= 0.2, hs + 0.3 * X^(1/2), 0.3 * zh * (1 - d/zh)) + d = FT(1.1) * zh * log(1 + X^FT(1/4)) + z0m = ifelse(0 <= X <= FT(0.2), hs + FT(0.3) * X^(1/2), FT(0.3) * zh * (1 - d/zh)) z0m_se = missing (;d, z0m, z0m_se) end -function roughness_parameters(::Roughness_wind_profile, ustar::AbstractVector{FT}, wind, psi_m; +function roughness_parameters(::Roughness_wind_profile, ustar::AbstractVector{FT}, wind, psi_m; zh, zr, d = 0.7*zh, constants=BigleafConstants() ) where FT FT == Missing && return((d = d, z0m = missing, z0m_se = missing)) diff --git a/src/unit_conversions.jl b/src/unit_conversions.jl index bb7a4f9..8e301a4 100644 --- a/src/unit_conversions.jl +++ b/src/unit_conversions.jl @@ -60,17 +60,17 @@ function Esat_slope(Tair::Number; Esat_formula=Sonntag1990(), constants=BigleafC Delta = Esat_from_Tair_deriv(Tair; Esat_formula, constants) Esat, Delta end, -function Esat_from_Tair(Tair; Esat_formula=Sonntag1990(), constants=BigleafConstants()) - a,b,c = get_EsatCoef(Esat_formula) - Esat = a * exp((b * Tair) / (c + Tair)) * constants.Pa2kPa +function Esat_from_Tair(Tair::FT; Esat_formula=Sonntag1990(), constants=BigleafConstants()) where FT + a,b,c = map(FT, get_EsatCoef(Esat_formula)) + Esat = a * exp((b * Tair) / (c + Tair)) * FT(constants.Pa2kPa) end, -function Esat_from_Tair_deriv(Tair; Esat_formula=Sonntag1990(), constants=BigleafConstants()) +function Esat_from_Tair_deriv(Tair::FT; Esat_formula=Sonntag1990(), constants=BigleafConstants()) where FT # slope of the saturation vapor pressure curve #Delta = eval(D(expression(a * exp((b * Tair) / (c + Tair))),name="Tair")) - a,b,c = get_EsatCoef(Esat_formula) + a,b,c = map(FT, get_EsatCoef(Esat_formula)) #Delta_Pa = @. a*(b / (Tair + c) + (-Tair*b) / ((Tair + c)^2))*exp((Tair*b) / (Tair + c)) Delta_Pa = @. a * (exp((b * Tair)/(c + Tair)) * (b/(c + Tair) - (b * Tair)/(c + Tair)^2)) - Delta = Delta_Pa .* constants.Pa2kPa + Delta = Delta_Pa .* FT(constants.Pa2kPa) end get_EsatCoef(::Sonntag1990) = (a=611.2,b=17.62,c=243.12) @@ -145,12 +145,12 @@ true """ function ms_to_mol(G_ms,Tair,pressure; constants=BigleafConstants()) Tair = Tair + oftype(Tair,constants.Kelvin) - pressure = pressure * constants.kPa2Pa + pressure = pressure * oftype(Tair, constants.kPa2Pa) G_mol = G_ms * pressure / (oftype(Tair, constants.Rgas) * Tair) end, function mol_to_ms(G_mol,Tair,pressure; constants=BigleafConstants()) Tair = Tair + oftype(Tair,constants.Kelvin) - pressure = pressure * constants.kPa2Pa + pressure = pressure * oftype(Tair, constants.kPa2Pa) G_ms = G_mol * (oftype(Tair, constants.Rgas) * Tair) / (pressure) end From 06ab098eeb13d71a7deadba519810d150cf2691d Mon Sep 17 00:00:00 2001 From: Olivier Bonte Date: Fri, 1 Aug 2025 10:57:41 +0200 Subject: [PATCH 2/3] Tair to Float Esat_from_Tair test --- test/unit_conversions_test.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/unit_conversions_test.jl b/test/unit_conversions_test.jl index e507759..04d272f 100644 --- a/test/unit_conversions_test.jl +++ b/test/unit_conversions_test.jl @@ -1,7 +1,7 @@ using Bigleaf, Test @testset "Esat_from_Tair" begin - Tair = 15 + Tair = 15.0 Esat_formula=Sonntag1990() constants=BigleafConstants() eSat = Esat_from_Tair(Tair; Esat_formula, constants) @@ -25,7 +25,7 @@ end constants=BigleafConstants() Esat = Esat_from_Tair.(Tair) delta = Esat_from_Tair_deriv.(Tair) - delta2 = diff(Esat)/step + delta2 = diff(Esat)/step @test all(isapprox.(delta[2:end] - delta2, 0, atol=1e-3)) Esat3, delta3 = Esat_slope(Tair[1]) @test Esat3 ≈ Esat[1] From a801ea70e6b0c548e92d5b9618036feea447c7f5 Mon Sep 17 00:00:00 2001 From: Olivier Bonte Date: Fri, 1 Aug 2025 13:50:09 +0200 Subject: [PATCH 3/3] Use convert function to change type --- src/surface_roughness.jl | 8 ++++---- src/unit_conversions.jl | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/surface_roughness.jl b/src/surface_roughness.jl index 7ac6c6f..c81f475 100755 --- a/src/surface_roughness.jl +++ b/src/surface_roughness.jl @@ -149,7 +149,7 @@ true ``` """ function roughness_parameters(::RoughnessCanopyHeight, zh::FT; - frac_d=FT(0.7), frac_z0m=FT(0.1)) where FT + frac_d=convert(FT,0.7), frac_z0m=convert(FT,0.1)) where FT d = frac_d*zh z0m = frac_z0m*zh z0m_se = missing @@ -157,10 +157,10 @@ function roughness_parameters(::RoughnessCanopyHeight, zh::FT; end function roughness_parameters(::RoughnessCanopyHeightLAI, zh::FT, LAI; - cd=FT(0.2), hs=FT(0.01)) where FT + cd=convert(FT,0.2), hs=convert(FT,0.01)) where FT X = cd * LAI - d = FT(1.1) * zh * log(1 + X^FT(1/4)) - z0m = ifelse(0 <= X <= FT(0.2), hs + FT(0.3) * X^(1/2), FT(0.3) * zh * (1 - d/zh)) + d = convert(FT,1.1) * zh * log(1 + X^convert(FT,1/4)) + z0m = ifelse(0 <= X <= convert(FT,0.2), hs + convert(FT,0.3) * X^(1/2), convert(FT,0.3) * zh * (1 - d/zh)) z0m_se = missing (;d, z0m, z0m_se) end diff --git a/src/unit_conversions.jl b/src/unit_conversions.jl index 152bd89..b8fe36f 100644 --- a/src/unit_conversions.jl +++ b/src/unit_conversions.jl @@ -61,16 +61,16 @@ function Esat_slope(Tair::Number; Esat_formula=Sonntag1990(), constants=BigleafC Esat, Delta end, function Esat_from_Tair(Tair::FT; Esat_formula=Sonntag1990(), constants=BigleafConstants()) where FT - a,b,c = map(FT, get_EsatCoef(Esat_formula)) + a,b,c = map(x -> convert(FT,x), get_EsatCoef(Esat_formula)) Esat = a * exp((b * Tair) / (c + Tair)) * FT(constants.Pa2kPa) end, function Esat_from_Tair_deriv(Tair::FT; Esat_formula=Sonntag1990(), constants=BigleafConstants()) where FT # slope of the saturation vapor pressure curve #Delta = eval(D(expression(a * exp((b * Tair) / (c + Tair))),name="Tair")) - a,b,c = map(FT, get_EsatCoef(Esat_formula)) + a,b,c = map(x -> convert(FT,x), get_EsatCoef(Esat_formula)) #Delta_Pa = @. a*(b / (Tair + c) + (-Tair*b) / ((Tair + c)^2))*exp((Tair*b) / (Tair + c)) Delta_Pa = a * (exp((b * Tair)/(c + Tair)) * (b/(c + Tair) - (b * Tair)/(c + Tair)^2)) - Delta = Delta_Pa .* FT(constants.Pa2kPa) + Delta = Delta_Pa .* convert(FT, constants.Pa2kPa) end get_EsatCoef(::Sonntag1990) = (a=611.2,b=17.62,c=243.12)