diff --git a/Project.toml b/Project.toml
index c5775ae..47a6ff0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,6 +19,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 TimeZones = "f269a46b-ccf7-5d73-abea-4c690281aa53"
 
 [compat]
@@ -35,6 +36,7 @@ RecursiveArrayTools = "2.20.0"
 StaticArrays = "1.2.13"
 StatsBase = "0.33.12"
 Suppressor = "0.2"
+Tables = "1.6.0"
 TimeZones = "1.6.1"
 julia = "1.6"
 
diff --git a/inst/explore_broadcast_missing.jl b/inst/explore_broadcast_missing.jl
new file mode 100644
index 0000000..2a59b09
--- /dev/null
+++ b/inst/explore_broadcast_missing.jl
@@ -0,0 +1,14 @@
+x = [1.0, 2.0, missing]
+y = [3.0, 3.0, 3.0]
+f(x,y) = x+y
+function g(x,y)
+    z = f.(x,y)
+    if z isa Union
+        return typeof(z)
+        convert(typeof(z).b.a, z)
+    else
+        z
+    end
+end
+z = g(x,y)
+@code_warntype g(x,y)
\ No newline at end of file
diff --git a/src/Bigleaf.jl b/src/Bigleaf.jl
index 52f2403..8138108 100644
--- a/src/Bigleaf.jl
+++ b/src/Bigleaf.jl
@@ -4,6 +4,7 @@ module Bigleaf
 using Optim
 using FillArrays
 using DataFrames
+using Tables: columntable
 using Dates, TimeZones
 using Pipe
 using AstroLib
@@ -37,6 +38,9 @@ export Monin_Obukhov_length, Monin_Obukhov_length!, stability_parameter,
     roughness_parameters, Reynolds_Number
 export aerodynamic_conductance!, add_Ga, add_Ga!, compute_Ram, compute_Ram!, roughness_z0h
 
+# shorthand type: either DataFrame or column table (Tables.columntable)
+DFTable = Union{NamedTuple,AbstractDataFrame}
+
 include("util.jl")    
 include("bigleaf_constants.jl")
 include("unit_conversions.jl")
diff --git a/src/bigleaf_constants.jl b/src/bigleaf_constants.jl
index 2d1189f..6cc31ac 100644
--- a/src/bigleaf_constants.jl
+++ b/src/bigleaf_constants.jl
@@ -39,6 +39,7 @@ Default values can be overridden by the named arguments:
 - frac2percent Conversion between fraction and percent
 """
 function bigleaf_constants(;
+  # make sure to provide Float64 values (rather than Int) so that return is type-stable
   ## Physical constants
   cp         = 1004.834,        # specific heat of air for constant pressure (J K-1 kg-1)
   Rgas       = 8.31451,         # universal gas constant (J mol-1 K-1)
@@ -49,7 +50,7 @@ function bigleaf_constants(;
   eps        = 0.622,           # ratio of the molecular weight of water vapor to dry air (=Mw/Md)
   g          = 9.81,            # gravitational acceleration (m s-2)
   solar_constant = 1366.1,      # solar constant, i.e. solar radation at earth distance from the sun (W m-2)
-  pressure0  = 101325,          # reference atmospheric pressure at sea level (Pa)
+  pressure0  = 101325.0,        # reference atmospheric pressure at sea level (Pa)
   Tair0      = 273.15,          # reference air temperature (K)
   k          = 0.41,            # von Karman constant
   Cmol       = 0.012011,        # molar mass of carbon (kg mol-1)
@@ -61,17 +62,17 @@ function bigleaf_constants(;
   ## Conversion constants
   Kelvin       = 273.15,         # conversion degree Celsius to Kelvin
   DwDc         = 1.6,            # Ratio of the molecular diffusivities for water vapor and CO2
-  days2seconds = 86400,          # seconds per day
-  kPa2Pa       = 1000,           # conversion kilopascal (kPa) to pascal (Pa)
+  days2seconds = 86400.0,        # seconds per day
+  kPa2Pa       = 1000.0,         # conversion kilopascal (kPa) to pascal (Pa)
   Pa2kPa       = 0.001,          # conversion pascal (Pa) to kilopascal (kPa)
-  umol2mol     = 1e-06,          # conversion micromole (umol) to mole (mol)
-  mol2umol     = 1e06,           # conversion mole (mol) to micromole (umol)
-  kg2g         = 1000,           # conversion kilogram (kg) to gram (g)
+  umol2mol     = 1.0e-6,         # conversion micromole (umol) to mole (mol)
+  mol2umol     = 1.0e6,          # conversion mole (mol) to micromole (umol)
+  kg2g         = 1000.0,         # conversion kilogram (kg) to gram (g)
   g2kg         = 0.001,          # conversion gram (g) to kilogram (kg)
-  kJ2J         = 1000,           # conversion kilojoule (kJ) to joule (J)
+  kJ2J         = 1000.0,         # conversion kilojoule (kJ) to joule (J)
   J2kJ         = 0.001,          # conversion joule (J) to kilojoule (kJ)
   se_median    = 1.253,          # conversion standard error (SE) of the mean to SE of the median (http://influentialpoints.com/Training/standard_error_of_median.htm)
-  frac2percent = 100             # conversion between fraction and percent
+  frac2percent = 100.0           # conversion between fraction and percent
 )
   Dict(
 :cp => cp, :Rgas => Rgas, :Rv => Rv, :Rd => Rd, :Md => Md, :Mw => Mw, :eps => eps, :g => g,
diff --git a/src/boundary_layer_conductance.jl b/src/boundary_layer_conductance.jl
index 51802a8..9a1ba65 100755
--- a/src/boundary_layer_conductance.jl
+++ b/src/boundary_layer_conductance.jl
@@ -41,10 +41,13 @@ true
 ``` 
 """
 function compute_Gb!(df::AbstractDataFrame, approach::Val{:Thom_1972}; kwargs...)
-  compute_Gb_!(df, approach, :ustar; kwargs...) # inputcols
+  ft(ustar) = Gb_Thom(ustar; kwargs...)
+  transform!(df, :ustar => ByRow(ft) => AsTable)
 end
-function compute_Gb!(df::AbstractDataFrame, approach::Val{:constant_kB1}; kwargs...)
-  compute_Gb_!(df, approach, :ustar; kwargs...) # inputcols
+function compute_Gb!(df::AbstractDataFrame, approach::Val{:constant_kB1}; kB_h, kwargs...)
+  # do not use ByRow because kb_H can be a vector
+  ft(ustar) = Gb_constant_kB1.(ustar, kB_h; kwargs...)
+  transform!(df, :ustar => ft => AsTable)
 end
 function compute_Gb!(df::AbstractDataFrame, approach::Val{:Choudhury_1988}; kwargs...)
   Gb_Choudhury!(df; kwargs...)
@@ -53,14 +56,6 @@ function compute_Gb!(df::AbstractDataFrame, approach::Val{:Su_2001}; kwargs...)
   Gb_Su!(df; kwargs...)
 end
 
-function compute_Gb_!(df::AbstractDataFrame, approach, inputcols; 
-  Sc::AbstractVector=SA[], kwargs...
-  )
-  fGb(args...) = compute_Gb(approach, args...; kwargs...)
-  transform!(df, :ustar => ByRow(fGb) => SA[:Rb_h, :Gb_h, :kB_h, :Gb_CO2])
-end
-compute_Gb(::Val{:Thom_1972}, args...; kwargs...) = Gb_Thom(args...; kwargs...)
-compute_Gb(::Val{:constant_kB1}, args...; kB_h, kwargs...) = Gb_constant_kB1(args..., kB_h; kwargs...)
 
 """
     add_Gb(Gb_h::Union{Missing,Number}, Sc::Vararg{Pair,N}; constants)
@@ -270,8 +265,8 @@ function Gb_Choudhury!(df::AbstractDataFrame; leafwidth, LAI, wind_zh=nothing,
     wind_zh = wind_profile(df, zh, d, z0m; zh, zr, stab_formulation, constants)
   end
   # Broadcasting does not work over keyword arguments, need to pass as positional
-  fwind(ustar, wind_zh; kwargs...) = Gb_Choudhury(ustar;wind_zh, kwargs...)
-  ft(ustar) = fwind.(ustar, wind_zh; leafwidth, LAI, constants)
+  fwind(ustar, wind_zh, leafwidth, LAI; kwargs...) = Gb_Choudhury(ustar;wind_zh, leafwidth, LAI, kwargs...)
+  ft(ustar) = fwind.(ustar, wind_zh, leafwidth, LAI; constants)
   transform!(df, :ustar => ft => AsTable)
 end
 
@@ -365,7 +360,7 @@ compute_Gb!(df,Val(:Su_2001), zh=25,zr=40,Dl=0.1,LAI=1.5)
 true
 ``` 
 """
-function Gb_Su(Tair,pressure,ustar; zh, wind_zh, Dl, fc, N=2, Cd=0.2, hs=0.01,
+function Gb_Su(Tair,pressure,ustar; wind_zh, Dl, fc, N=2, Cd=0.2, hs=0.01,
   constants=bigleaf_constants()
   )
   wind_zh = max(0.01,wind_zh) ## avoid zero windspeed
@@ -390,13 +385,13 @@ function Gb_Su!(df::AbstractDataFrame; wind_zh=nothing, Dl, fc=nothing,
   end
   if isnothing(fc)
     isnothing(LAI) && error("one of 'fc' or 'LAI' must be provided")
-    fc = (1-exp(-LAI/2)) 
+    fc = length(LAI) == 1 ? (1-exp(-LAI/2)) : @. (1-exp(-LAI/2))
   end
   isnothing(Dl) && error(
     "need to provide keyword argument Dl with :Su_2001 method")
   inputcols = SA[:Tair,:pressure,:ustar]
   # Broadcasting does not work over keyword arguments, need to pass as positional
-  fwind(wind_zh, args...; kwargs...) = Gb_Su(args...; wind_zh, kwargs...)
-  ft(args...) = fwind.(wind_zh, args...; zh, Dl, fc, N=2, Cd=0.2, hs=0.01, constants)
+  fwind(wind_zh, Dl, fc, N, Cd, hs, args...; kwargs...) = Gb_Su(args...; wind_zh, Dl, fc, N, Cd, hs, kwargs...)
+  ft(args...) = fwind.(wind_zh, Dl, fc, N, Cd, hs, args...; constants)
   transform!(df, inputcols => ft => AsTable)
 end
diff --git a/src/stability_correction.jl b/src/stability_correction.jl
index 580cdcc..52160b7 100755
--- a/src/stability_correction.jl
+++ b/src/stability_correction.jl
@@ -56,9 +56,10 @@ function Monin_Obukhov_length!(df;constants=bigleaf_constants())
   ft(args...) = Monin_Obukhov_length(args...; constants)
   transform!(df, SA[:Tair, :pressure, :ustar, :H] => ByRow(ft) => :MOL)
 end
-function Monin_Obukhov_length(df;constants=bigleaf_constants())
-  Monin_Obukhov_length!(copy(df, copycols = false); 
-    constants).MOL
+function Monin_Obukhov_length(df::DFTable; kwargs...)
+    tmp = Monin_Obukhov_length.(df.Tair, df.pressure, df.ustar, df.H; kwargs...)
+    # eltype_agg = Union{eltype(df.Tair), eltype(df.pressure), eltype(df.ustar), eltype(df.H)}
+    # convert(Vector{eltype_agg}, tmp)
 end
 
 """
@@ -110,7 +111,7 @@ function stability_parameter!(df::AbstractDataFrame; zr,d, MOL=nothing,
   df[!,:zeta] .= stability_parameter(df; zr,d,MOL,constants)
   df
 end
-function stability_parameter(df::AbstractDataFrame; zr,d, MOL=nothing,
+function stability_parameter(df::DFTable; zr,d, MOL=nothing,
   constants=bigleaf_constants())
   if isnothing(MOL) MOL = Monin_Obukhov_length(df; constants); end
   stability_parameter.(zr,d,MOL)
@@ -219,7 +220,7 @@ function get_stability_coefs_unstable(::Val{:Dyer_1970}, zeta)
 end
 
 #TODO z or zr here?
-function stability_correction(Tair,pressure,ustar,H, z,d; 
+function stability_correction(Tair::Union{Missing,Number},pressure,ustar,H, z,d; 
   stab_formulation=Val(:Dyer_1970), constants = bigleaf_constants())
   stab_formulation isa Val{:no_stability_correction} && return(
     (psi_h = zero(z), psi_m = zero(z)))
@@ -227,6 +228,14 @@ function stability_correction(Tair,pressure,ustar,H, z,d;
   zeta  = stability_parameter(z,d,MOL)
   stability_correction(zeta; stab_formulation)
 end
+function stability_correction(Tair,pressure,ustar,H, z,d; kwargs...) 
+  Tables.columns(
+    stability_correction.(Tair,pressure,ustar,H, z,d; kwargs...))
+end
+function stability_correction(df::DFTable, z,d; kwargs...)   
+  tmp1 = stability_correction.(df.Tair, df.pressure, df.ustar, df.H, z, d; kwargs...)
+  tmp2 = Tables.columns(tmp1)
+end
 
 function stability_correction!(df; zeta=df.zeta, 
   stab_formulation=Val(:Dyer_1970), constants = bigleaf_constants())
@@ -242,6 +251,6 @@ end
 
 function stability_correction!(df, z, d; 
   stab_formulation=Val(:Dyer_1970), constants = bigleaf_constants())
-  ft(args...) = stability_correction(args..., z, d; stab_formulation, constants)
-  transform!(df, SA[:Tair,:pressure,:ustar,:H] => ByRow(ft) => AsTable)
+  ft(args...) = stability_correction.(args..., z, d; stab_formulation, constants)
+  transform!(df, SA[:Tair,:pressure,:ustar,:H] => ft => AsTable)
 end
diff --git a/src/surface_roughness.jl b/src/surface_roughness.jl
index f36b16f..bebdd10 100755
--- a/src/surface_roughness.jl
+++ b/src/surface_roughness.jl
@@ -142,7 +142,8 @@ function roughness_parameters(::Val{:canopy_height}, zh; frac_d=0.7, frac_z0m=0.
   (;d, z0m, z0m_se)
 end
 
-function roughness_parameters(::Val{:canopy_height_LAI}, zh, LAI; cd=0.2, hs=0.01)
+function roughness_parameters(::Val{:canopy_height_LAI}, zh, LAI; 
+  cd=0.2, hs=0.01)
   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))
@@ -152,13 +153,12 @@ end
 
 # docu: supply psi_m = 0 for no stability correction, default method
 # if psi_m is given df only needs wind and ustar
-function roughness_parameters(::Val{:wind_profile}, df, zh, zr;
+function roughness_parameters(::Val{:wind_profile}, df::DFTable, zh, zr;
   d = 0.7*zh, psi_m = nothing, stab_formulation=Val(:Dyer_1970), 
   constants=bigleaf_constants()
   )
   if isnothing(psi_m)
-    psi_m = stability_correction!(copy(df, copycols=false), zr, d; 
-      stab_formulation, constants).psi_m
+    psi_m = stability_correction(df, zr, d; stab_formulation, constants).psi_m
   end
   z0m_all = allowmissing(@. (zr - d) * exp(-constants[:k]*df.wind / df.ustar - psi_m))
   #z0m_all[(z0m_all .> zh)] .= missing # problems with missings
@@ -263,7 +263,7 @@ function wind_profile(stab_formulation, z, ustar, Tair,pressure,H, d, z0m;
   wind_profile(z, ustar, d, z0m, psi_m)
 end
 
-function wind_profile(df::AbstractDataFrame, z, d, z0m = nothing; 
+function wind_profile(df::DFTable, z, d, z0m = nothing; 
   zh = nothing, zr = nothing, 
   stab_formulation = Val(:Dyer_1970), constants = bigleaf_constants()
   )
@@ -275,13 +275,12 @@ function wind_profile(df::AbstractDataFrame, z, d, z0m = nothing;
     # uses psi_m at zr
     z0m = roughness_parameters(Val(:wind_profile), df, zh, zr).z0m
   end
-  psi_m = stability_correction!(
-    copy(df, copycols=false), z, d; stab_formulation, constants).psi_m
+  psi_m = stability_correction(df, z, d; stab_formulation, constants).psi_m
   #wind_profile(df, z, d, z0m, psi_m; constants)
   wind_profile(df, z, d, z0m, psi_m; constants)
 end
 
-function wind_profile(df::AbstractDataFrame, z, d, z0m, psi_m::AbstractVector;
+function wind_profile(df::DFTable, z, d, z0m, psi_m::AbstractVector;
   constants = bigleaf_constants())
   # when psi_m is given, df is not used any more, but keep for consistency
   wind_profile.(z, df.ustar, d, z0m, psi_m; constants)
diff --git a/test/Project.toml b/test/Project.toml
index 749d32a..595a8a4 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -7,5 +7,6 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimeZones = "f269a46b-ccf7-5d73-abea-4c690281aa53"
diff --git a/test/aerodynamic_conductance.jl b/test/aerodynamic_conductance.jl
index 5b25c33..dd18205 100644
--- a/test/aerodynamic_conductance.jl
+++ b/test/aerodynamic_conductance.jl
@@ -1,8 +1,10 @@
 @testset "add_Ga" begin
     Gb_h = 0.0347
     Ga_m = 0.3
-    @test add_Ga(Gb_h, Ga_m) == NamedTuple()
+    @test @inferred add_Ga(Gb_h, Ga_m) == NamedTuple()
     #
+    # here cannot @inferred the names in the return
+    # TODO
     Ga = add_Ga(Gb_h, Ga_m, :N20 => 2, :CH4 => 4)
     @test keys(Ga) == (:Ga_N20, :Ga_CH4)
     @test ≈(Ga.Ga_N20, 0.0164, rtol = 1/100)
@@ -15,9 +17,9 @@
     # DataFrame variant
     dfo = DataFrame(Gb_h=SA[Gb_h, missing, 0.055], Ga_m = SA[0.3,0.3,0.3])
     df = copy(dfo)
-    df2 = add_Gb!(df)
+    df2 = @inferred add_Gb!(df)
     @test isequal(df2, dfo)
-    df2 = add_Ga!(df, :N20 => 2, :CH4 => 4)
+    df2 = @inferred add_Ga!(df, :N20 => 2, :CH4 => 4)
     @test df2 === df # mutating
     @test propertynames(df2)[end-1:end] == [:Ga_N20, :Ga_CH4]
     @test df2.Ga_N20[1] == Ga.Ga_N20
@@ -30,73 +32,100 @@ end
     stability_parameter!(df; zr, d)
     stability_correction!(df)
     z0m = roughness_parameters(Val(:wind_profile), df, zh, zr; psi_m = df.psi_m).z0m
-    Ram_r = compute_Ram(Val(:wind_zr), ustar, wind)
+    Ram_r = @inferred compute_Ram(Val(:wind_zr), ustar, wind)
     @test Ram_r ≈ 6.31 rtol = 1/100 # TODO check with R
-    Ram_p = compute_Ram(Val(:wind_profile), ustar; zr, d, z0m, psi_h = df.psi_h[24])
+    Ram_p =  @inferred compute_Ram(Val(:wind_profile), ustar; zr, d, z0m, psi_h = df.psi_h[24])
     @test Ram_p ≈ 5.41 rtol = 1/100 # TODO check with R
-
+    #
     #Gb_h, Gb_CO2 = compute_Gb!(df, Val(:Thom_1972))[24, [:Gb_h, :Gb_CO2]]
-    compute_Ram!(df, Val(:wind_zr))
+    @inferred compute_Ram!(df, Val(:wind_zr))
     #@test all(propertynames(df)[(end+1-length(Ram)):end] .== keys(Ram))
     @test propertynames(df)[end] == :Ra_m
     @test df.Ra_m[24] == Ram_r
-    compute_Ram!(df, Val(:wind_profile); zr, d, z0m)
+    @inferred compute_Ram!(df, Val(:wind_profile); zr, d, z0m)
     @test df.Ra_m[24] == Ram_p
+    Ram_skalar = df.Ra_m
+    #
+    # test zr as a vecttor
+    df[!,:zri] .= zr
+    df.zri[1] = zr/2
+    @inferred compute_Ram!(df, Val(:wind_profile); zr=df.zri, d, z0m)
+    @test df.Ra_m[2:end] == Ram_skalar[2:end]
+    @test df.Ra_m[1] != Ram_skalar[1]
 end
 
 @testset "aerodynamic_conductance! only ustar and wind" begin
     # using only ustar and wind
     dfo = DataFrame(ustar = SA[0.1,missing,0.3], wind = [3.4, 2.8, 3.3])
     df = copy(dfo)
-    aerodynamic_conductance!(df; Gb_model=Val(:Thom_1972))
+    @inferred aerodynamic_conductance!(df; Gb_model=Val(:Thom_1972))
     # missing due to missing zr
     @test all(ismissing.(df.psi_h))        
     @test all(ismissing.(df.psi_m))        
-    aerodynamic_conductance!(df; Gb_model=Val(:Thom_1972), zr = thal.zr, zh = thal.zh
+    @inferred aerodynamic_conductance!(df; Gb_model=Val(:Thom_1972), zr = thal.zr, zh = thal.zh
         ,stab_formulation = Val(:no_stability_correction))
     @test all(iszero.(df.psi_h))        
     @test all(iszero.(df.psi_m))        
     @test propertynames(df)[(end-8):end] == 
         SA[:Rb_h, :Gb_h, :kB_h, :Gb_CO2, :Ra_m, :Ga_m, :Ga_h, :Ra_h, :Ga_CO2]
-    # TODO compare with R        
+    # compare with R        
+    @test all(isapproxm.(df.Gb_h, (0.0347, missing, 0.0723), rtol = 1e-3))
+    @test all(isapproxm.(df.Ga_m, (0.00294, missing, 0.0273), rtol = 1e-3))
+    @test all(isapproxm.(df.Ga_h, (0.00271, missing, 0.0198), rtol = 1e-3))
 end
 
 @testset "aerodynamic_conductance! Gb_Thom" begin
     df = copy(tha48)
-    aerodynamic_conductance!(df; Gb_model=Val(:Thom_1972), zh = thal.zh, zr = thal.zr)
+    # test Ga_m by wind_profile
+    @inferred aerodynamic_conductance!(df; Gb_model=Val(:Thom_1972), Ram_model=Val(:wind_profile), 
+        zh = thal.zh, zr = thal.zr)
     @test propertynames(df)[(end-8):end] == 
         SA[:Rb_h, :Gb_h, :kB_h, :Gb_CO2, :Ra_m, :Ga_m, :Ga_h, :Ra_h, :Ga_CO2]
 end
 
 @testset "aerodynamic_conductance! constant_kB1" begin
-    kB_h = 1.18
-    # DataFrame variant
+    # specify kB1 as a vecttor
     df = copy(tha48)
-    aerodynamic_conductance!(df; Gb_model=Val(:constant_kB1), kB_h, zh = thal.zh, zr = thal.zr)
+    df[!,:kB_hi] .= 1.18
+    df.kB_hi[1] = 1.18/2
+    @inferred aerodynamic_conductance!(df; Gb_model=Val(:constant_kB1), kB_h = df.kB_hi, 
+        zh = thal.zh, zr = thal.zr)
     @test propertynames(df)[(end-8):end] == 
         SA[:Rb_h, :Gb_h, :kB_h, :Gb_CO2, :Ra_m, :Ga_m, :Ga_h, :Ra_h, :Ga_CO2]
+    @test all(isapproxm.(df.Gb_h[1:3], (0.375, 0.17, 0.167), rtol = 1e-2))
 end
 
 @testset "aerodynamic_conductance! Gb_Choudhury" begin
     leafwidth=0.1
     df = copy(tha48)
-    aerodynamic_conductance!(df; Gb_model=Val(:Choudhury_1988),
+    # test Ga_m by wind_profile
+    @inferred aerodynamic_conductance!(df; Gb_model=Val(:Choudhury_1988), Ram_model=Val(:wind_profile), 
         leafwidth, LAI=thal.LAI, zh = thal.zh, zr = thal.zr)
     @test propertynames(df)[(end-8):end] == 
         SA[:Rb_h, :Gb_h, :kB_h, :Gb_CO2, :Ra_m, :Ga_m, :Ga_h, :Ra_h, :Ga_CO2]
+    # compare with R        
+    @test all(isapproxm.(df.Gb_h[1:3], (0.157, 0.15, 0.151), rtol = 1e-2))
+    @test all(isapproxm.(df.Ga_m[1:3], (0.0709, 0.0649, 0.0603), rtol = 1e-2))
 end
 
 @testset "aerodynamic_conductance! Gb_Su" begin
     Dl=0.01
     df = copy(tha48)
-    aerodynamic_conductance!(df; Gb_model=Val(:Su_2001), Dl, LAI=thal.LAI, zh=thal.zh, zr=thal.zr)
+    df[!,:LAI] .= thal.LAI
+    df.LAI[1] = thal.LAI/2
+    # test providing changing LAI
+    @inferred aerodynamic_conductance!(df; Gb_model=Val(:Su_2001), 
+        Dl, LAI=df.LAI, zh=thal.zh, zr=thal.zr);
     @test propertynames(df)[(end-8):end] == 
         SA[:Rb_h, :Gb_h, :kB_h, :Gb_CO2, :Ra_m, :Ga_m, :Ga_h, :Ra_h, :Ga_CO2]
+    # compare with R
+    @test all(isapproxm.(df.Gb_h[1:3], (0.202, 0.177, 0.167), rtol = 1e-2))
+    @test all(isapproxm.(df.Ga_m[1:3], (0.0693, 0.0538, 0.0507), rtol = 1e-2))
 end
 
 @testset "roughness_z0h" begin
     z0m, kB_h = 1.2, 3.4
-    z0h = roughness_z0h(z0m, kB_h)
+    z0h = @inferred roughness_z0h(z0m, kB_h)
     @test z0h ≈ 0.0400 rtol = 1e-2
 end
 
diff --git a/test/boundary_layer_conductance.jl b/test/boundary_layer_conductance.jl
index 9bc0756..70f8478 100644
--- a/test/boundary_layer_conductance.jl
+++ b/test/boundary_layer_conductance.jl
@@ -14,9 +14,9 @@
     # DataFrame variant
     dfo = DataFrame(Gb_h=SA[Gb_h, missing, 0.055])
     df = copy(dfo)
-    df2 = add_Gb!(df)
+    df2 = @inferred add_Gb!(df)
     @test isequal(df2, dfo)
-    df2 = add_Gb!(df, :N20 => 2, :CH4 => 4)
+    df2 = @inferred add_Gb!(df, :N20 => 2, :CH4 => 4)
     @test df2 === df # mutating
     @test propertynames(df2)[end-1:end] == [:Gb_N20, :Gb_CH4]
     @test df2.Gb_N20[1] == Gb.Gb_N20
@@ -24,10 +24,10 @@ end
 
 @testset "compute_Gb Gb_constant_kB1" begin
     kB_h = 1.18
-    Gb = Gb_constant_kB1(0.1, kB_h)
+    Gb = @inferred Gb_constant_kB1(0.1, kB_h)
     @test keys(Gb) == (:Rb_h, :Gb_h, :kB_h, :Gb_CO2)
     @test all(isapprox.(values(Gb), values((Rb_h = 28.8, Gb_h = 0.0347, kB_h = 1.18, Gb_CO2 = 0.0264)), rtol = 1e-2))
-    Gb = Gb_constant_kB1(missing, kB_h)
+    Gb = @inferred Gb_constant_kB1(missing, kB_h)
     @test keys(Gb) == (:Rb_h, :Gb_h, :kB_h, :Gb_CO2)
     @test Gb.kB_h == kB_h
     @test all(ismissing.(values(getindex.(Ref(Gb), SA[:Rb_h, :Gb_h, :Gb_CO2]))))
@@ -35,13 +35,13 @@ end
     # DataFrame variant
     dfo = DataFrame(ustar = SA[0.1,missing,0.3])
     df = copy(dfo)
-    compute_Gb!(df, Val(:constant_kB1); kB_h)
+    @inferred compute_Gb!(df, Val(:constant_kB1); kB_h)
     @test propertynames(df) == [:ustar, :Rb_h, :Gb_h, :kB_h, :Gb_CO2]
     @test all(isapprox.(values(df[1,2:end]), values((Rb_h = 28.8, Gb_h = 0.0347, kB_h = 1.18, Gb_CO2 = 0.0264)), rtol = 1e-2))
 end
 
 @testset "compute_Gb Gb_Thom" begin
-    Gb = Gb_Thom(0.1)
+    Gb = @inferred Gb_Thom(0.1)
     @test keys(Gb) == (:Rb_h, :Gb_h, :kB_h, :Gb_CO2)
     @test all(isapprox.(values(Gb), values((Rb_h = 28.8, Gb_h = 0.0347, kB_h = 1.18, Gb_CO2 = 0.0264)), rtol = 1e-2))
     Gb = Gb_Thom(missing)
@@ -51,7 +51,7 @@ end
     # DataFrame variant
     dfo = DataFrame(ustar = SA[0.1,missing,0.3])
     df = copy(dfo)
-    compute_Gb!(df, Val(:Thom_1972))
+    @inferred compute_Gb!(df, Val(:Thom_1972))
     @test propertynames(df) == [:ustar, :Rb_h, :Gb_h, :kB_h, :Gb_CO2]
     @test all(isapprox.(values(df[1,2:end]), values((Rb_h = 28.8, Gb_h = 0.0347, kB_h = 1.18, Gb_CO2 = 0.0264)), rtol = 1e-2))
 end
@@ -60,16 +60,16 @@ end
     zh, zr, LAI = thal.zh, thal.zr, thal.LAI
     leafwidth=0.1
     wind_zh = 1.2
-    Gb = Gb_Choudhury(tha48.ustar[24]; leafwidth, LAI, wind_zh)
+    Gb = @inferred Gb_Choudhury(tha48.ustar[24]; leafwidth, LAI, wind_zh)
     @test keys(Gb) == (:Rb_h, :Gb_h, :kB_h, :Gb_CO2)
     @test Gb.Rb_h ≈ 8.533 rtol=1e-3 # regression to first implementation
-    Gbm = Gb_Choudhury(missing; leafwidth, LAI, wind_zh)
+    Gbm = @inferred Gb_Choudhury(missing; leafwidth, LAI, wind_zh)
     @test keys(Gb) == (:Rb_h, :Gb_h, :kB_h, :Gb_CO2)
     @test Gbm.Rb_h == Gb.Rb_h
     #
     # DataFrame variant
     df = copy(tha48)
-    compute_Gb!(df, Val(:Choudhury_1988); leafwidth, LAI, zh, zr)
+    @inferred compute_Gb!(df, Val(:Choudhury_1988); leafwidth, LAI, zh, zr)
     @test propertynames(df)[(end-3):end] == SA[:Rb_h, :Gb_h, :kB_h, :Gb_CO2]
     all(isapprox.(values(df[24,SA[:Rb_h, :Gb_h, :kB_h, :Gb_CO2]]), 
         (7.534, 0.1327, 2.255, 0.1008), rtol=1e-3))
@@ -80,16 +80,16 @@ end
     Dl=0.01; leafwidth=0.1; fc = (1-exp(-LAI/2)) 
     wind_zh = 1.2
     Tair, pressure, ustar = values(tha48[24, SA[:Tair, :pressure, :ustar]])
-    Gb = Gb_Su(Tair, pressure, ustar; zh, wind_zh, Dl, fc)
+    Gb = @inferred Gb_Su(Tair, pressure, ustar; wind_zh, Dl, fc)
     @test keys(Gb) == (:Rb_h, :Gb_h, :kB_h, :Gb_CO2)
     @test Gb.Rb_h ≈ 1.221 rtol=1e-3 # regression to first implementation
-    Gb = Gb_Su(missing, pressure, ustar; zh, wind_zh, Dl, fc)
+    Gb = @inferred Gb_Su(missing, pressure, ustar; wind_zh, Dl, fc)
     @test keys(Gb) == (:Rb_h, :Gb_h, :kB_h, :Gb_CO2)
     @test all(ismissing.(values(Gb)))
     #
     # DataFrame variant
     df = copy(tha48)
-    compute_Gb!(df, Val(:Su_2001); Dl, LAI, zh, zr)
+    @inferred compute_Gb!(df, Val(:Su_2001); Dl, LAI, zh, zr)
     @test propertynames(df)[(end-3):end] == SA[:Rb_h, :Gb_h, :kB_h, :Gb_CO2]
     all(isapprox.(values(df[24,SA[:Rb_h, :Gb_h, :kB_h, :Gb_CO2]]), 
         (1.767, 0.5659, 0.5289, 0.4299), rtol=1e-3)) # regression to first implementation
diff --git a/test/meteorological_variables.jl b/test/meteorological_variables.jl
index f454725..7d0c713 100644
--- a/test/meteorological_variables.jl
+++ b/test/meteorological_variables.jl
@@ -1,35 +1,35 @@
 @testset "air_density" begin
-  ad = air_density(25.0,100.0) # Tair, pressure
+  ad = @inferred air_density(25.0,100.0) # Tair, pressure
   # regression test
   @test ≈(ad, 1.168, atol =0.001)
 end
 
 @testset "pressure_from_elevation" begin
-  pressure = pressure_from_elevation(500.0, 25.0) # elev, Tair
+  pressure = @inferred pressure_from_elevation(500.0, 25.0) # elev, Tair
   # regression test
   @test ≈(pressure, 95.681, atol =0.001)
   #
-  pressure = pressure_from_elevation(500.0, 25.0, 1.5) # elev, Tair, VPD
+  pressure = @inferred pressure_from_elevation(500.0, 25.0, 1.5) # elev, Tair, VPD
   # regression test
   @test ≈(pressure, 95.717, atol =0.001)
 end
 
 @testset "virtual_temp" begin
   Tair,pressure,VPD = 25,100,1.5
-  vt = virtual_temp(Tair,pressure,VPD)  
+  vt = @inferred virtual_temp(Tair,pressure,VPD)  
   # regression test
   @test ≈(vt, 26.9, atol =0.1)
 end
 
 @testset "kinematic_viscosity" begin
   Tair,pressure = 25,100
-  vis = kinematic_viscosity(Tair,pressure)
+  vis = @inferred kinematic_viscosity(Tair,pressure)
   # regression test
   @test ≈(vis, 1.58e-5, atol =1e-7)
 end
 
 @testset "air_density" begin
-  ad = air_density(25.0,100.0) # Tair, pressure
+  ad = @inferred air_density(25.0,100.0) # Tair, pressure
   # regression test
   @test ≈(ad, 1.168, atol =0.001)
 end
@@ -41,11 +41,11 @@ end
   VPD = 1.5
   Tair = 25.0
   accuracy = 1e-2
-  ea = VPD_to_e(VPD,Tair)
-  Td = dew_point_from_e(ea,accuracy = accuracy)                
+  ea = @inferred VPD_to_e(VPD,Tair)
+  Td = @inferred dew_point_from_e(ea,accuracy = accuracy)                
   @test ≈(ea, Esat_from_Tair(Td), atol = accuracy)
   @test ≈(dew_point(Tair, VPD; accuracy = accuracy), Td, atol = accuracy)
-  Td1 = @test_logs (:warn,r"set to 1")  dew_point(Tair, VPD; accuracy = 1.2)
+  Td1 = @test_logs (:warn,r"set to 1")  @inferred dew_point(Tair, VPD; accuracy = 1.2)
   @test ≈(Td1, Td, atol = 1.0)
 end
 
@@ -55,12 +55,12 @@ end
   pressure = 100
   accuracy = 1e-2
   #
-  gamma  = psychrometric_constant(Tair,pressure)
+  gamma  = @inferred psychrometric_constant(Tair,pressure)
   ea     = VPD_to_e(VPD,Tair)
-  Tw = wetbulb_temp_from_e_Tair_gamma(ea, Tair, gamma; accuracy = accuracy)
+  Tw = @inferred wetbulb_temp_from_e_Tair_gamma(ea, Tair, gamma; accuracy = accuracy)
   #@test ≈(ea, Esat_from_Tair(Tw) - gamma* (Tair - Tw), atol = accuracy)
   @test ≈(ea, Esat_from_Tair(Tw) - 0.93*gamma* (Tair - Tw), atol = accuracy)
   @test ≈(wetbulb_temp(Tair, pressure, VPD; accuracy = accuracy), Tw, atol = accuracy)
-  Tw1 = @test_logs (:warn,r"set to 1") wetbulb_temp(Tair, pressure, VPD; accuracy = 1.2)
+  Tw1 = @test_logs (:warn,r"set to 1") @inferred wetbulb_temp(Tair, pressure, VPD; accuracy = 1.2)
   @test ≈(Tw1, Tw, atol = 1.0)
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index feb745c..ce64676 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,6 +2,7 @@ using Bigleaf
 using Test, StableRNGs
 using Pipe: @pipe 
 using DataFrames, Dates, TimeZones
+using Tables: columntable
 using StaticArrays
 using Statistics, StatsBase
 
diff --git a/test/stability_correction.jl b/test/stability_correction.jl
index 87fd35b..53dcecf 100644
--- a/test/stability_correction.jl
+++ b/test/stability_correction.jl
@@ -1,31 +1,50 @@
 @testset "Monin_Obukhov_length" begin
     datetime, ustar, Tair, pressure, H = values(tha48[24,1:5])
-    MOL24 = Monin_Obukhov_length(Tair, pressure, ustar, H)
+    MOL24 = @inferred Monin_Obukhov_length(Tair, pressure, ustar, H)
+    MOL24 = @inferred Monin_Obukhov_length(Tair, pressure, ustar, H; constants = bigleaf_constants())
     @test ≈(MOL24, -104.3, rtol = 1/1000)
     df = copy(tha48)
-    Monin_Obukhov_length!(df)
+    dfd = disallowmissing(df)
+    dfm = allowmissing(df)
+    @inferred Monin_Obukhov_length!(df)
     @test df.MOL[24] == MOL24
-    MOL = Monin_Obukhov_length(df)
+    MOL = Monin_Obukhov_length(columntable(df); constants = bigleaf_constants()) 
+    @test MOL == df.MOL
+    # with non-missings types can be inferred, but with missings its not inferrable
+    #@code_warntype Monin_Obukhov_length(columntable(dfd); constants = bigleaf_constants()) 
+    MOL = @inferred Monin_Obukhov_length(columntable(dfd); constants = bigleaf_constants()) 
     @test MOL == df.MOL
 end
 
 @testset "stability_parameter" begin
     df = copy(tha48)
+    dfd = allowmissing(df)
     MOL24 = first(Monin_Obukhov_length(df[SA[24],:]))
     zr=40;d=15
-    zeta = stability_parameter(zr,d,MOL24)
+    zeta = @inferred stability_parameter(zr,d,MOL24)
     @test ≈(zeta , -0.240, rtol = 1/100)
     #
-    stability_parameter!(df;zr,d)
+    @inferred stability_parameter!(df;zr,d)
     @test df.zeta[24] == zeta
+    zeta_scalar = df.zeta
+    # 
+    # non-mutatingvariant
+    # also test zr as a vector
     df = copy(tha48)
-    zetas = stability_parameter(df;zr,d)
-    @test df == tha48 # did not modify original df
-    @test zetas[24] == zeta
+    df[!,:zri] .= zr
+    df.zri[1] = zr/2
+    df_bak = copy(df)
+    # again type-stable only with disallowmissing
+    #@code_warntype stability_parameter(columntable(disallowmissing(df));zr=df.zri,d) 
+    zetas = @inferred stability_parameter(columntable(disallowmissing(df));zr=df.zri,d) 
+    @test df == df_bak # did not modify original df
+    @test zetas[2:24] == zeta_scalar[2:24]
+    @test zetas[1] != zeta_scalar[1]
 end
 
 @testset "stability_correction" begin
     zeta = -2:0.5:0.5
+    @inferred stability_correction(first(zeta))
     df2 = DataFrame(stability_correction.(zeta))
     @test all(isapprox.(df2.psi_h, SA[2.431,2.197,1.881,1.386,0,-2.5], rtol=1e-3))
     @test all(isapprox.(df2.psi_m, SA[2.275, 2.061, 1.772, 1.317, 0, -2.5], rtol=1e-3))
@@ -33,11 +52,11 @@ end
     @test all(isapprox.(df2.psi_h, SA[2.085, 1.862, 1.564, 1.106,0, -3.9], rtol=1e-3))
     @test all(isapprox.(df2.psi_m, SA[2.418, 2.200, 1.904, 1.435,0, -3], rtol=1e-3))
     #
-    resm = stability_correction(missing)
+    resm = @inferred stability_correction(missing)
     @test keys(resm) == (:psi_h, :psi_m)
     @test all(ismissing.(values(resm)))
     #
-    resm = stability_correction(first(zeta); 
+    resm = @inferred stability_correction(first(zeta); 
         stab_formulation = Val(:no_stability_correction))
     @test resm == (psi_h = 0.0, psi_m = 0.0)
 end
@@ -46,36 +65,47 @@ end
     zr=40;d=15
     dfo = DataFrame(Tair=25, pressure=100, ustar=0.2:0.1:1.0, H=40:20:200)
     df = copy(dfo)
-    stability_correction!(df, zr, d)
+    # for type stability, use columntable(df)
+    tmp = @inferred stability_correction(columntable(df), zr, d)
+    @inferred stability_correction!(df, zr, d)
     propertynames(df)[(end-1):end] == SA[:psi_h, :psi_m]
     #
     dfm = copy(dfo)
-    stability_correction!(dfm, zr, d; stab_formulation = Val(:no_stability_correction))
+    @inferred stability_correction!(dfm, zr, d; stab_formulation = Val(:no_stability_correction))
     propertynames(dfm)[(end-1):end] == SA[:psi_h, :psi_m]
     @test all(iszero.(dfm.psi_h))
     @test all(iszero.(dfm.psi_m))
     #
+    # test computing zeta first
     df2 = copy(dfo)
     stability_parameter!(df2; zr, d) # adds zeta
     stability_correction!(df2)
     @test df2.psi_h == df.psi_h
     @test df2.psi_m == df.psi_m
+    #
+    # test specifying zr as a vector
+    df3 = copy(dfo)
+    df3[!,:zri] .= zr
+    df3.zri[1] = zr/2
+    @inferred stability_correction!(df3, df3.zri, d)
+    @test df3.psi_h[2:end] == df.psi_h[2:end]
+    @test df3.psi_h[1] != df.psi_h[1]
 end
 
 @testset "stability_correction from raw" begin
     datetime, ustar, Tair, pressure, H = values(tha48[24,1:5])
     z=40.0;d=15.0
-    resm = stability_correction(Tair,pressure,ustar,H, z,d)
+    resm = @inferred stability_correction(Tair,pressure,ustar,H, z,d)
     @test resm.psi_h ≈ 0.940 rtol=1/1000
     @test resm.psi_m ≈ 0.902 rtol=1/1000
     #
-    resm0 = stability_correction(Tair,pressure,ustar,H, z,d; 
+    resm0 = @inferred stability_correction(Tair,pressure,ustar,H, z,d; 
         stab_formulation = Val(:no_stability_correction))
     @test resm0 == (psi_h = 0.0, psi_m = 0.0)
     #
     df = copy(tha48)
     df.ustar[3] = missing
-    stability_correction!(df,z,d)
+    @inferred stability_correction!(df,z,d)
     @test all(ismissing.((df.psi_h[3], df.psi_m[3])))
     @test all(isapprox.((df.psi_h[24], df.psi_m[24]),values(resm)))
-end
\ No newline at end of file
+end
diff --git a/test/surface_roughness.jl b/test/surface_roughness.jl
index 264dbfc..268902b 100644
--- a/test/surface_roughness.jl
+++ b/test/surface_roughness.jl
@@ -1,29 +1,32 @@
 @testset "Reynolds_Number" begin
     Tair,pressure,ustar,z0m = 25,100,0.5,0.5
-    R = Reynolds_Number(Tair,pressure,ustar,z0m)                             
+    R = @inferred Reynolds_Number(Tair,pressure,ustar,z0m)                             
     @test ≈(R, 15870, rtol=1e-3) 
 end
 
-
 @testset "roughness_parameters" begin
     zh = thal.zh
     zr = thal.zr
     LAI = thal.LAI
     keys_exp = (:d, :z0m, :z0m_se)
-    rp = roughness_parameters(Val(:canopy_height), zh)
+    rp = @inferred roughness_parameters(Val(:canopy_height), zh)
     #round.(values(rp); sigdigits = 4)
     @test keys(rp) == keys_exp
     @test all(isapproxm.(values(rp), (18.55, 2.65, missing), rtol=1e-3))
     #
-    rp = roughness_parameters(Val(:canopy_height_LAI), zh, LAI)
+    rp = @inferred roughness_parameters(Val(:canopy_height_LAI), zh, LAI)
     #round.(values(rp); sigdigits = 4)
     @test keys(rp) == keys_exp
     @test all(isapproxm.(values(rp), (21.77, 1.419, missing), rtol=1e-3))
     #
     df = copy(tha48)
+    dfd = disallowmissing(df)
+    #df.wind[1] = missing
     d=0.7*zh
-    psi_m = stability_correction!(copy(df, copycols=false), zr, d).psi_m
-    rp = roughness_parameters(Val(:wind_profile), df, zh, zr; psi_m)
+    psi_m = (@inferred stability_correction(columntable(dfd), zr, d)).psi_m
+    psi_m = stability_correction(df, zr, d).psi_m
+    # note: must use columntable for type stability - but needs compilation timede
+    rp = @inferred roughness_parameters(Val(:wind_profile), columntable(dfd), zh, zr; psi_m)
     #round.(values(rp); sigdigits = 4)
     @test keys(rp) == keys_exp
     #@test all(isapproxm.(values(rp), (18.55, 1.879, 0.3561), rtol=1e-3))
@@ -31,13 +34,16 @@ end
     @test all(isapproxm.(values(rp), (18.55, 1.879402, 0.356108), rtol=1e-3))
     #
     # no stability correction
-    rp0 = roughness_parameters(Val(:wind_profile), df, zh, zr; psi_m = 0.0)
+    rp0 = @inferred roughness_parameters(Val(:wind_profile), columntable(dfd), zh, zr; 
+        stab_formulation = Val(:no_stability_correction))
     @test keys(rp0) == keys_exp
     # same magnitude as with stability correction
     @test all(isapproxm.(values(rp0), values(rp), rtol=0.5))
     #
     # estimate psi
-    rp_psiauto = roughness_parameters(Val(:wind_profile), df, zh, zr)
+    #@code_warntype stability_correction(columntable(df), zr, 0.7*zh)
+    #@code_warntype roughness_parameters(Val(:wind_profile), columntable(df), zh, zr)
+    rp_psiauto = @inferred roughness_parameters(Val(:wind_profile), columntable(dfd), zh, zr)
     @test propertynames(df) == propertynames(tha48) # not changed
     @test rp_psiauto == rp
 end
@@ -47,33 +53,37 @@ end
     z = 30
     d=0.7*thal.zh
     z0m= 2.65
-    u30 = wind_profile(z, ustar, d, z0m)
-    @test ≈(u30, 1.93, rtol = 1/100 )
+    u30 = @inferred wind_profile(z, ustar, d, z0m)
+    @test ≈(u30, 1.93, rtol = 1/100 ) # from R
     #
-    u30c = wind_profile(Val(:Dyer_1970), z, ustar, Tair,pressure, H, d, z0m)
-    @test ≈(u30c, 2.31, rtol = 1/100 )
+    u30c = @inferred wind_profile(Val(:Dyer_1970), z, ustar, Tair,pressure, H, d, z0m)
+    @test ≈(u30c, 2.31, rtol = 1/100 ) # from R
     #
     z0m=1.9 #2.14 #2.65
-    u30 = wind_profile(z, ustar, d, z0m) # used below
-    u30c = wind_profile(Val(:Dyer_1970), z, ustar, Tair,pressure, H, d, z0m)
+    u30 = @inferred wind_profile(z, ustar, d, z0m) # used below
+    u30c = @inferred wind_profile(Val(:Dyer_1970), z, ustar, Tair,pressure, H, d, z0m)
     df = copy(tha48)
+    dfd = disallowmissing(df)
+    windz = @inferred wind_profile(columntable(dfd), z, d, z0m; stab_formulation = Val(:no_stability_correction))    
     windz = wind_profile(df, z, d, z0m; stab_formulation = Val(:no_stability_correction))    
     @test length(windz) == 48
     @test windz[1] == u30
-    windzc = wind_profile(df, z, d, z0m; stab_formulation = Val(:Dyer_1970))    
+    windzc = @inferred wind_profile(columntable(dfd), z, d, z0m; stab_formulation = Val(:Dyer_1970))    
     @test windzc[1] == u30c
     #plot(windz)
     #plot!(windz2)
-    psi_m = stability_correction!(copy(df, copycols=false), z, d).psi_m
-    windzc2 = wind_profile(df, z, d, z0m, psi_m)    
+    psi_m = stability_correction(df, z, d).psi_m
+    windzc2 = @inferred wind_profile(columntable(dfd), z, d, z0m, psi_m)    
     @test windzc2 == windzc
     #
     # estimate z0m
     # need to give zh and zr in addition to many variables in df
     @test_throws Exception wind_profile(df, z, d)    
-    windzc3 = wind_profile(df, z, d; zh=thal.zh, zr=thal.zr)    
+    windzc3 = @inferred wind_profile(columntable(dfd), z, d; zh=thal.zh, zr=thal.zr, 
+        stab_formulation = Val(:Dyer_1970))    
     # may have used slightly different estimated z0m
     #windzc3 - windzc
     @test all(isapprox.(windzc3, windzc, atol=0.1))
+    @test windzc3[1] ≈ 2.764203 rtol=1e-3 # from R
 end