diff --git a/src/surface_roughness.jl b/src/surface_roughness.jl index 7422b03..c81f475 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=convert(FT,0.7), frac_z0m=convert(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=convert(FT,0.2), hs=convert(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 = 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 -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 924aa8d..b8fe36f 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(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; 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(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 .* constants.Pa2kPa + Delta = Delta_Pa .* convert(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 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]