Fix/type stability

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))

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
 

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]