diff --git a/docs/src/walkthrough.md b/docs/src/walkthrough.md
index bd88b57..5107f81 100644
--- a/docs/src/walkthrough.md
+++ b/docs/src/walkthrough.md
@@ -345,7 +345,7 @@ We start with the simplest version, in which $G_b$ is calculated empirically bas
 the friction velocity ($u_*$) according to Thom 1972:
 
 ```@example doc
-aerodynamic_conductance!(thas)
+aerodynamic_conductance!(thas);
 thas[1:3, Cols(:datetime,Between(:zeta,:Ga_CO2))]
 ```
 
@@ -365,7 +365,7 @@ dimension ($D_l$, assumed to be 1cm here), and information on sensor and canopy
 
 ```@example doc
 aerodynamic_conductance!(thas;Gb_model=Val(:Su_2001),
-     LAI=thal.zh, zh=thal.zh, d=0.7*thal.zh, zr=thal.zr,Dl=thal.Dl)
+     LAI=thal.zh, zh=thal.zh, d=0.7*thal.zh, zr=thal.zr, Dl=thal.Dl);
 thas[1:3, Cols(:datetime,Between(:zeta,:Ga_CO2))]
 ```
 
@@ -383,8 +383,8 @@ Functin `add_Gb` calculates $G_b$ for other trace gases, provided that the respe
 number is known. 
 
 ```@example doc
-compute_Gb!(thas, Val(:Thom_1972)) # adds/modifies column Gb_h and Gb_CO2
-add_Gb!(thas, :Gb_O2 => 0.84, :Gb_CH4 => 0.99) # adds Gb_O2 and Gb_CH4
+compute_Gb!(thas, Val(:Thom_1972)); # adds/modifies column Gb_h and Gb_CO2
+add_Gb!(thas, :Gb_O2 => 0.84, :Gb_CH4 => 0.99); # adds Gb_O2 and Gb_CH4
 select(first(thas,3), r"Gb_")
 ```
 
@@ -403,10 +403,10 @@ close to the surface and weaker at greater heights:
 using Statistics
 wind_heights = 22:2:60.0
 d = 0.7 * thal.zh
-z0m = roughness_parameters(Val(:wind_profile), thas, thal.zh, thal.zr).z0m
+z0m = roughness_parameters(Val(:wind_profile), thas; zh=thal.zh, zr=thal.zr).z0m
 wp = map(wind_heights) do z
-  wind_profile(thas,z,d, z0m; zh=thal.zh, zr=thal.zr)
-end
+  wind_profile(z, thas,d, z0m; zh=thal.zh, zr=thal.zr)
+end;
 nothing # hide
 ```
 ```@setup doc
diff --git a/src/Bigleaf.jl b/src/Bigleaf.jl
index 8138108..b00c3c9 100644
--- a/src/Bigleaf.jl
+++ b/src/Bigleaf.jl
@@ -31,7 +31,8 @@ export setinvalid_range!, setinvalid_qualityflag!,
     setinvalid_nongrowingseason!, get_growingseason, setinvalid_afterprecip!
 export decoupling, surface_conductance, aerodynamic_conductance
 export compute_Gb, compute_Gb!, add_Gb, add_Gb!, 
-    Gb_Thom, Gb_Choudhury, Gb_Su, Gb_constant_kB1
+    Gb_Thom, Gb_Choudhury, Gb_Su, Gb_constant_kB1,
+    compute_Gb_quantities, compute_Gb_quantities!
 export wind_profile
 export Monin_Obukhov_length, Monin_Obukhov_length!, stability_parameter, 
     stability_parameter!, stability_correction, stability_correction!,
diff --git a/src/aerodynamic_conductance.jl b/src/aerodynamic_conductance.jl
index b6df2b9..d1609bf 100755
--- a/src/aerodynamic_conductance.jl
+++ b/src/aerodynamic_conductance.jl
@@ -83,32 +83,37 @@ function aerodynamic_conductance!(df; Gb_model = Val(:Thom_1972), Ram_model = Va
   zr=nothing,zh=nothing, d = isnothing(zh) ? nothing : 0.7*zh ,
   z0m=nothing,Dl=nothing,N=2,fc=nothing,LAI=nothing,Cd=0.2,hs=0.01,
   leafwidth=nothing,
-  wind_profile=false,
   stab_formulation=Val(:Dyer_1970),
   kB_h=nothing,constants=bigleaf_constants()
   )
   # add zeta
   if !isnothing(zr) && !isnothing(d) && !(stab_formulation isa Val{:no_stability_correction})
-    stability_parameter!(df::AbstractDataFrame; zr,d, constants)
+    stability_parameter!(df::AbstractDataFrame; z=zr, d, constants)
   else
     df[!,:zeta] .= missing
   end
   # adds columns psi_m and psi_h, (:no_stability_correction: just add 0s without using zeta)
-  stability_correction!(df; stab_formulation, constants)
+  stability_correction!(df; zeta=df.zeta, stab_formulation, constants)
   # pre-estimate z0m to use it in both Gb and Ga
   needs_windprofile = 
     Gb_model isa Union{Val{:Choudhury_1988}, Val{:Su_2001}} || 
     Ram_model isa Val{:wind_profile}
-  if isnothing(z0m) && needs_windprofile
-    z0m = roughness_parameters(Val(:wind_profile), df, zh, zr; psi_m = df.psi_m).z0m
+  if needs_windprofile
+    if isnothing(z0m) 
+      z0m = roughness_parameters(Val(:wind_profile), df; zh, zr, psi_m = df.psi_m).z0m
+    end
+    wind_zh = wind_profile(zh, df, d, z0m; zh, zr, stab_formulation, constants)
   end
+  #
   # calculate canopy boundary layer conductance (Gb)
   Gb_model isa Val{:Thom_1972} && compute_Gb!(df, Gb_model; constants)
   Gb_model isa Val{:constant_kB1} && compute_Gb!(df, Gb_model; kB_h, constants)
-  Gb_model isa Val{:Choudhury_1988} && compute_Gb!(df, Gb_model; 
-    leafwidth, LAI, zh, zr, d, z0m, stab_formulation, constants)
-  Gb_model isa Val{:Su_2001} && compute_Gb!(df, Gb_model; 
-    Dl, fc, N, Cd, hs, z0m, zh, zr, d, LAI, stab_formulation, constants)
+  Gb_model isa Val{:Choudhury_1988} && compute_Gb!(
+    df, Gb_model; leafwidth, LAI, wind_zh, constants)
+  Gb_model isa Val{:Su_2001} && compute_Gb!(
+    df, Gb_model; wind_zh, Dl, fc, N, Cd, hs, LAI, constants)
+  compute_Gb_quantities!(df)
+  #
   # calculate aerodynamic risistance for momentum (Ra_m)
   Ram_model isa Val{:wind_profile} && compute_Ram!(df, Ram_model; zr, d, z0m, constants) 
   Ram_model isa Val{:wind_zr} && compute_Ram!(df, Ram_model) 
@@ -161,7 +166,6 @@ Estimate bulk aerodynamic conductance.
 
 # Arguments
 - `ustar`          : Friction velocity (m s-1)
-- `wind`           : wind speed at measurement height (m s-1)
 - `df`             : DataFrame with above columns
 - `zr`             : Instrument (reference) height (m)
 - `d`              : Zero-plane displacement height (-), can be estimated using 
diff --git a/src/boundary_layer_conductance.jl b/src/boundary_layer_conductance.jl
index 9a1ba65..cbcfa80 100755
--- a/src/boundary_layer_conductance.jl
+++ b/src/boundary_layer_conductance.jl
@@ -4,7 +4,6 @@
 Estimate boundary layer conductance.
 
 # Arguments  
-For `Val(:Thom_1972)`
 - `df`       : DataFrame with required variables (depend on approach)
 - `approach` : one of
   - `Val(:Thom_1972)`: see [`Gb_Thom`](@ref)
@@ -12,42 +11,38 @@ For `Val(:Thom_1972)`
   - `Val(:Su_2001)`: see [`Gb_Su`](@ref)
   - `Val(:constant_kB1)`: see [`Gb_constant_kB1`](@ref)
  
-The different approaches required different variables present in `df` and
+The different approaches require different variables to be present in `df` and
 different keyword arguments.
 
 # Value
 updated DataFrame `df` with the following columns:
 - `Gb_h`: Boundary layer conductance for heat transfer (m s-1)
-- `Rb_h`: Boundary layer resistance for heat transfer (s m-1)
-- `kB_h`: kB-1 parameter for heat transfer
-- `Gb_CO2`: Boundary layer conductance for CO2 (m s-1). 
 
-To subsequently compute conductances for other species with different 
-Schmidt numbers see [`add_Gb!`](@ref).
+To subsequently derived quantities see
+- [`compute_Gb_quantities`](@ref) for Resistance, kB-1 constant, and CO2 conductance
+- [`add_Gb!`](@ref) for conductances of other species given their Schmidt numbers.
 
 # See also
-[`Gb_Thom`](@ref), `Gb_Choudhury`](@ref), [`Gb_Su`](@ref), `Gb_Choudhury`](@ref), 
+[`Gb_Thom`](@ref), [`Gb_Choudhury`](@ref), [`Gb_Su`](@ref), [`Gb_constant_kB1`](@ref), 
 [`aerodynamic_conductance!`](@ref)
 
 ```jldoctest; output = false
 using DataFrames
-df = DataFrame(Tair=25,pressure=100,wind=[3,4,5],
-  ustar=[0.5,0.6,0.65],H=[200,230,250]) 
-leafwidth=0.01;LAI=5;zh=25;zr=40
-compute_Gb!(df, Val(:Choudhury_1988); leafwidth, LAI, zh, zr)
-≈(df.Gb_h[1], 0.379, rtol=1e-3)
+df = DataFrame(wind=[3,4,5], ustar=[0.5,0.6,0.65]) 
+compute_Gb!(df, Val(:Thom_1972))
+≈(df.Gb_h[1], 0.102, rtol=1e-2)
 # output
 true
 ``` 
 """
 function compute_Gb!(df::AbstractDataFrame, approach::Val{:Thom_1972}; kwargs...)
   ft(ustar) = Gb_Thom(ustar; kwargs...)
-  transform!(df, :ustar => ByRow(ft) => AsTable)
+  transform!(df, :ustar => ByRow(ft) => :Gb_h)
 end
 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)
+  transform!(df, :ustar => ft => :Gb_h)
 end
 function compute_Gb!(df::AbstractDataFrame, approach::Val{:Choudhury_1988}; kwargs...)
   Gb_Choudhury!(df; kwargs...)
@@ -56,12 +51,41 @@ function compute_Gb!(df::AbstractDataFrame, approach::Val{:Su_2001}; kwargs...)
   Gb_Su!(df; kwargs...)
 end
 
+"""
+    compute_Gb_quantities(Gb_h)
+    compute_Gb_quantities!(df:::AbstractDataFrame)
+
+Based on boundary layer conductance for heat, compute derived quantities.
+
+# Arguments
+- `Gb_h` : Boundary layer conductance for heat transfer (m s-1)
+- `df`   : DataFrame with above columns
+- `constants=`[`bigleaf_constants`](@ref)`()`: entries `Sc_CO2` and `Pr`
+
+# Value
+NamedTuple with entries
+- `Gb_h`: Boundary layer conductance for heat transfer (m s-1)
+- `Rb_h`: Boundary layer resistance for heat transfer (s m-1)
+- `kB_h`: kB-1 parameter for heat transfer
+- `Gb_CO2`: Boundary layer conductance for CO2 (m s-1). 
+"""
+function compute_Gb_quantities(Gb_h, ustar; constants=bigleaf_constants())
+  Rb_h = 1/Gb_h
+  kB_h = Rb_h*constants[:k]*ustar
+  Gb_CO2 = Gb_h / (constants[:Sc_CO2]/constants[:Pr])^0.67
+  (;Rb_h, Gb_h, kB_h, Gb_CO2)
+end
+function compute_Gb_quantities!(df::AbstractDataFrame; constants=bigleaf_constants())
+  ft(args...) = compute_Gb_quantities.(args...; constants)
+  transform!(df, SA[:Gb_h, :ustar] => ft => AsTable)
+end
+
 
 """
     add_Gb(Gb_h::Union{Missing,Number}, Sc::Vararg{Pair,N}; constants)
     add_Gb!(df::AbstractDataFrame, Sc::Vararg{Pair,N}; Gb_h = df.Gb_h, kwargs...) 
 
-compute additional boundary layer conductance quantities for given Schmidt-numbers
+compute boundary layer conductance for additional quantities for given Schmidt-numbers
 
 # Arguments  
 - `Gb_h`     : Boundary layer conductance for heat transfer (m s-1)
@@ -115,6 +139,8 @@ function get_names_and_values(prefix::AbstractString, Sc::Vararg{Pair,N}) where
   Scn, Scv
 end
 
+
+
 """
     Gb_Thom(ustar; constants)
     compute_Gb!(df, Val{:Thom_1972})
@@ -130,7 +156,7 @@ for heat transfer based on a simple ustar (friction velocity) dependency.
 # Details
 The empirical equation for Rb suggested by Thom 1972 is:
  
-``Rb = 6.2 ustar^{-0.67}``
+``Rb = 6.2 {u^*}^{-0.67}``
  
 Gb (=1/Rb) for water vapor and heat are assumed to be equal in this package.
 
@@ -148,15 +174,12 @@ see [`compute_Gb!`](@ref)
 using DataFrames
 df = DataFrame(ustar = SA[0.1,missing,0.3])
 compute_Gb!(df, Val(:Thom_1972))
-propertynames(df) == [:ustar, :Rb_h, :Gb_h, :kB_h, :Gb_CO2]
+propertynames(df) == [:ustar, :Gb_h]
 ``` 
 """
 function Gb_Thom(ustar::Union{Missing,Number}; constants=bigleaf_constants())
   Rb_h = 6.2*ustar^-0.667
   Gb_h = 1/Rb_h
-  kB_h = Rb_h*constants[:k]*ustar
-  Gb_CO2 = Gb_h / (constants[:Sc_CO2]/constants[:Pr])^0.67
-  (;Rb_h, Gb_h, kB_h, Gb_CO2)
 end
 
 """
@@ -177,38 +200,22 @@ Rb_h computed by ``kB_h/(k * ustar)``, where k is the von Karman constant.
 function Gb_constant_kB1(ustar, kB_h; constants=bigleaf_constants())
   Rb_h = kB_h/(constants[:k] * ustar)
   Gb_h = 1/Rb_h
-  Gb_CO2 = Gb_h / (constants[:Sc_CO2]/constants[:Pr])^0.67
-  (;Rb_h, Gb_h, kB_h, Gb_CO2)
+  # Gb_CO2 = Gb_h / (constants[:Sc_CO2]/constants[:Pr])^0.67
+  # (;Rb_h, Gb_h, kB_h, Gb_CO2)
 end
   
 """
-    Gb_Choudhury(ustar; leafwidth, LAI, wind_zh, constants)
-      alpha   = 4.39 - 3.97*exp(-0.258*LAI)
-    Gb_Choudhury!(df::AbstractDataFrame; leafwidth, LAI, wind_zh=nothing, zh, zr, d=0.7*zh,
-      z0m = nothing, stab_formulation=Val(:Dyer_1970), constants)
+    Gb_Choudhury(; leafwidth, LAI, wind_zh, constants)
+    Gb_Choudhury!(df; leafwidth, LAI, wind_zh, constants)
 
 Estimate the canopy boundary layer conductance 
 for heat transfer according to Choudhury & Monteith 1988.
 
 # Arguments  
-- `Tair`             : Air temperature (degC)
-- `pressure`         : Atmospheric pressure (kPa)
-- `wind`             : Wind speed at sensor height (m s-1)
-- `ustar`            : Friction velocity (m s-1)
-- `H`                : Sensible heat flux (W m-2)
-- `df`               : DataFrame or matrix containing the above variables
+- `df`               : DataFrame where `Gb_h` is to be added/updated
 - `leafwidth`        : Leaf width (m)
 - `LAI`              : One-sided leaf area index
-- `wind_zh`          : Wind speed at canopy heihgt (m s-1), if not given then computed by
-                       [`wind_profile`](@ref) using the arguments below
-- `zh`               : Canopy height (m)
-- `zr`               : Instrument (reference) height (m)
-- `d`                : Zero-plane displacement height (-), can be calculated using 
-                      `roughness_parameters`
-- `z0m`              : Roughness length for momentum (m). If not provided, calculated 
-                        from `roughness_parameters`  within `wind_profile` using zh and zr
-- `stab_formulation` : Stability correction function used, 
-                       see [`stability_correction`](@ref).
+- `wind_zh`          : Wind speed at canopy heihgt (m s-1), see [`wind_profile`](@ref) 
 - `constants=`[`bigleaf_constants`](@ref)`()`
                         
 # Value
@@ -218,23 +225,18 @@ see [`compute_Gb!`](@ref)
 Boundary layer conductance according to Choudhury & Monteith 1988 is
 given by:
 
-``Gb_h = LAI((2a/\\alpha)*\\sqrt(u(zh)/w)*(1-exp(-\\alpha/2)))``
+``Gb_h = LAI \\left( 2a/\\alpha \\sqrt{u(z_h)/w} (1-e^{-\\alpha/2})\\right)``
 
 where ``\\alpha`` is modeled as an empirical relation to LAI (McNaughton & van den Hurk 1995):
 
-``\\alpha = 4.39 - 3.97*exp(-0.258*LAI)``
+``\\alpha = 4.39 - 3.97 e^{-0.258 \\, LAI}``
 
 ``w`` is leafwidth and ``u(zh)`` is the wind speed at the canopy surface.
 
 It can be approximated from measured wind speed at sensor height zr and a wind extinction 
-coefficient ``\\alpha``: ``u(zh) = u(zr) / (exp(\\alpha(zr/zh -1)))``.
+coefficient ``\\alpha``: ``u(z_h) = u(z_r) / e^{\\alpha(z_r/z_h -1)}``.
 However, here (if not explicitly given) it is estimated by [`wind_profile`](@ref)
 
-Gb (=1/Rb) for water vapor and heat are assumed to be equal in this package.
-Gb for other quantities x is calculated as (Hicks et al. 1987):
-``Gb_x = Gb / (Sc_x / Pr)^0.67``
-where Sc_x is the Schmidt number of quantity x, and Pr is the Prandtl number (0.71).
-         
 # References
 - Choudhury, B. J., Monteith J_L., 1988: A four-layer model for the heat
   budget of homogeneous land surfaces. Q. J. R. Meteorol. Soc. 114, 373-398.
@@ -245,36 +247,33 @@ where Sc_x is the Schmidt number of quantity x, and Pr is the Prandtl number (0.
   A preliminary multiple resistance routine for deriving dry deposition velocities
   from measured quantities. Water, Air, and Soil Pollution 36, 311-330.
 """
-function Gb_Choudhury(ustar; leafwidth, LAI, wind_zh, constants=bigleaf_constants())
+function Gb_Choudhury(; leafwidth, LAI, wind_zh, constants=bigleaf_constants())
   alpha   = 4.39 - 3.97*exp(-0.258*LAI)
   # (ismissing(wind_zh) || isnothing(wind_zh)) && return(
   #   (Rb_h=missing, Gb_h=missing, kB_h=missing, Gb_CO2=missing))
   wind_zh = max(0.01, wind_zh) ## avoid zero windspeed
   Gb_h = LAI*((0.02/alpha)*sqrt(wind_zh/leafwidth)*(1-exp(-alpha/2)))
-  # TODO facture out derving 3 others form Gb_h?
-  Rb_h = 1/Gb_h
-  kB_h = Rb_h*constants[:k]*ustar
-  Gb_CO2 = Gb_h / (constants[:Sc_CO2]/constants[:Pr])^0.67
-  (;Rb_h, Gb_h, kB_h, Gb_CO2)
+  # Rb_h = 1/Gb_h
+  # kB_h = Rb_h*constants[:k]*ustar
+  # Gb_CO2 = Gb_h / (constants[:Sc_CO2]/constants[:Pr])^0.67
+  # (;Rb_h, Gb_h, kB_h, Gb_CO2)
 end
-function Gb_Choudhury!(df::AbstractDataFrame; leafwidth, LAI, wind_zh=nothing, 
-  zh, zr, d=0.7*zh, z0m = nothing, stab_formulation=Val(:Dyer_1970),
-  constants=bigleaf_constants()
+function Gb_Choudhury!(
+  df::AbstractDataFrame; leafwidth, LAI, wind_zh, constants=bigleaf_constants()
   )
-  if isnothing(wind_zh)
-    wind_zh = wind_profile(df, zh, d, z0m; zh, zr, stab_formulation, constants)
-  end
+  # if isnothing(wind_zh)
+  #   wind_zh = wind_profile(zh, df, d, z0m; zh, zr, stab_formulation, constants)
+  # end
   # Broadcasting does not work over keyword arguments, need to pass as positional
-  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)
+  fwind(wind_zh, leafwidth, LAI; kwargs...) = Gb_Choudhury(;wind_zh, leafwidth, LAI, kwargs...)
+  ft() = fwind.(wind_zh, leafwidth, LAI; constants)
+  transform!(df, [] => ft => :Gb_h)
 end
 
 """
-    Gb_Su(Tair,pressure,ustar; zh, wind_zh, Dl, fc, N=2, Cd=0.2, hs=0.01, constants)
-    Gb_Su!(df; wind_zh=nothing, Dl, fc=nothing, N=2, Cd=0.2, hs=0.01,
-      z0m = nothing, zh, zr, d = 0.7*zh, LAI, stab_formulation=Val(:Dyer_1970), constants)
-
+    Gb_Su(Tair,pressure,ustar; wind_zh, Dl, fc, N=2, Cd=0.2, hs=0.01, constants)
+    Gb_Su!(df; wind_zh, Dl, fc=nothing, N=2, Cd=0.2, hs=0.01, LAI, constants)
+      
 Estimate Boundary Layer Conductance to heat transfer using the physically based 
 formulation according to Su et al. 2001. 
 
@@ -282,26 +281,13 @@ formulation according to Su et al. 2001.
 - `Tair`      : Air temperature (degC)
 - `pressure`  : Atmospheric pressure (kPa)
 - `ustar`     : Friction velocity (m s-1)
-- `H`         : Sensible heat flux (W m-2)
-- `wind`      : Wind speed at sensor height (m s-1)
 - `df`        : DataFrame or matrix containing the above variables
 - `Dl`        : Leaf characteristic dimension (m)
 - `fc`        : Fractional vegetation cover [0-1] (if not provided, calculated from LAI)
-- `LAI`       : One-sided leaf area index (-)
+- `LAI`       : One-sided leaf area index (-) - alternative to `fc`.
 - `N`         : Number of leaf sides participating in heat exchange (defaults to 2)
 - `Cd`        : Foliage drag coefficient (-)
 - `hs`        : Roughness height of the soil (m)
-- `LAI`       : One-sided leaf area index
-- `wind_zh`   : Wind speed at canopy heihgt (m s-1) (if not given then computed by
-                       [`wind_profile`](@ref) using the arguments below)
-- `zh`        : Canopy height (m)
-- `zr`        : Instrument (reference) height (m)
-- `d`         : Zero-plane displacement height (-), can be calculated using 
-         `roughness_parameters`
-- `z0m`       : Roughness length for momentum (m). If not provided, calculated 
-           from `roughness_parameters`  within `wind_profile` using zh and zr
-- `stab_formulation` : Stability correction function used, 
-                       see [`stability_correction`](@ref).
 - `constants=`[`bigleaf_constants`](@ref)`()`
 
 # Value
@@ -311,31 +297,27 @@ see [`compute_Gb!`](@ref)
 The formulation is based on the kB-1 model developed by Massman 1999. 
 Su et al. 2001 derived the following approximation:
  
-``kB-1 = (k Cd fc^2) / (4Ct ustar/u(zh)) + kBs-1(1 - fc)^2``
+``k_{B1} = (k C_d f_c^2) / (4C_t u^*/u(z_h)) + k_{Bs-1}(1 - f_c)^2``
 
-If fc (fractional vegetation cover) is missing, it is estimated from LAI:
-``fc = 1 - exp(-LAI/2)``
+If ``f_c`` (fractional vegetation cover) is missing, it is estimated from LAI:
+``f_c = 1 - e^{-LAI/2}``
 
 The wind speed at the top of the canopy is calculated using function
 [`wind_profile`](@ref).
 
 Ct is the heat transfer coefficient of the leaf (Massman 1999):
 
-``Ct = Pr^-2/3 Reh^-1/2 N``
+``C_t = P_r^{-2/3} R_{eh}^{-1/2}  N``
 
-where Pr is the Prandtl number (set to 0.71), and Reh is the Reynolds number for leaves:
+where ``P_r`` is the Prandtl number (set to 0.71), 
+and ``R_{eh}`` is the Reynolds number for leaves:
 
-``Reh = Dl wind(zh) / v``
+``R_{eh} = D_l \\, wind(z_h) / v``
  
-kBs-1, the kB-1 value for bare soil surface, is calculated according 
+k_{Bs-1}, the k_{B-1} value for bare soil surface, is calculated according 
 to Su et al. 2001:
 
-``kBs^-1 = 2.46(Re)^0.25 - ln(7.4)``
-
-Gb (=1/Rb) for water vapor and heat are assumed to be equal in this package.
-Gb for other quantities x is calculated as (Hicks et al. 1987):
-  ``Gb_x = Gb / (Sc_x / Pr)^0.67``
-where Sc_x is the Schmidt number of quantity x, and Pr is the Prandtl number (0.71).
+``k_{Bs-1} = 2.46(Re)^{0.25} - ln(7.4)``
 
 # References
 - Su, Z., Schmugge, T., Kustas, W. & Massman, W., 2001: An evaluation of
@@ -350,11 +332,13 @@ where Sc_x is the Schmidt number of quantity x, and Pr is the Prandtl number (0.
 ```jldoctest; output = false
 using DataFrames
 df = DataFrame(Tair=25,pressure=100,wind=[3,4,5],ustar=[0.5,0.6,0.65],H=[200,230,250]) 
-compute_Gb!(df,Val(:Su_2001), zh=25,zr=40,Dl=0.01,LAI=5)
+zh = 25; zr = 40
+wind_zh = wind_profile(zh, df, 0.7*zh; zh, zr)
+compute_Gb!(df,Val(:Su_2001); wind_zh, Dl=0.01, LAI=5)
 # the same meteorological conditions, but larger leaves
-compute_Gb!(df,Val(:Su_2001), zh=25,zr=40,Dl=0.1,LAI=5)
+compute_Gb!(df,Val(:Su_2001); wind_zh, Dl=0.1,LAI=5)
 # same conditions, large leaves, and sparse canopy cover (LAI = 1.5)
-compute_Gb!(df,Val(:Su_2001), zh=25,zr=40,Dl=0.1,LAI=1.5)
+compute_Gb!(df,Val(:Su_2001); wind_zh, Dl=0.1,LAI=1.5)
 ≈(df.Gb_h[1], 0.0638, rtol=1e-3)
 # output
 true
@@ -373,16 +357,10 @@ function Gb_Su(Tair,pressure,ustar; wind_zh, Dl, fc, N=2, Cd=0.2, hs=0.01,
   kB_h = (constants[:k]*Cd)/(4*Ct*ustar/wind_zh)*fc^2 + kBs*(1 - fc)^2
   Rb_h = kB_h/(constants[:k]*ustar)
   Gb_h = 1/Rb_h
-  Gb_CO2 = Gb_h / (constants[:Sc_CO2]/constants[:Pr])^0.67
-  (;Rb_h, Gb_h, kB_h, Gb_CO2)
 end
-function Gb_Su!(df::AbstractDataFrame; wind_zh=nothing, Dl, fc=nothing, 
-  N=2, Cd=0.2, hs=0.01, z0m = nothing, zh, zr, d = 0.7*zh, LAI,
-  stab_formulation=Val(:Dyer_1970), constants=bigleaf_constants()
+function Gb_Su!(df::AbstractDataFrame; wind_zh, Dl, fc=nothing, 
+  N=2, Cd=0.2, hs=0.01, LAI, constants=bigleaf_constants()
   )
-  if isnothing(wind_zh)
-    wind_zh = wind_profile(df, zh, d, z0m; zh, zr, stab_formulation, constants)
-  end
   if isnothing(fc)
     isnothing(LAI) && error("one of 'fc' or 'LAI' must be provided")
     fc = length(LAI) == 1 ? (1-exp(-LAI/2)) : @. (1-exp(-LAI/2))
@@ -393,5 +371,5 @@ function Gb_Su!(df::AbstractDataFrame; wind_zh=nothing, Dl, fc=nothing,
   # Broadcasting does not work over keyword arguments, need to pass as positional
   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)
+  transform!(df, inputcols => ft => :Gb_h)
 end
diff --git a/src/stability_correction.jl b/src/stability_correction.jl
index 52160b7..3806f0c 100755
--- a/src/stability_correction.jl
+++ b/src/stability_correction.jl
@@ -56,37 +56,35 @@ 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::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
+# function Monin_Obukhov_length(df::DFTable; kwargs...)
+#     tmp = Monin_Obukhov_length.(df.Tair, df.pressure, df.ustar, df.H; kwargs...)
+# end
 
 """
-    stability_parameter(zr,d,MOL)
-    stability_parameter!(df::AbstractDataFrame; zr,d, MOL=nothing, constants)
-    stability_parameter(df::AbstractDataFrame; zr,d, MOL=nothing, constants)
+    stability_parameter(z,d,MOL)
+    stability_parameter(z,d,Tair, pressure, ustar, H; constants)
+    stability_parameter!(df::AbstractDataFrame; z,d, MOL=nothing, constants)
 
 calculates stability parameter "zeta", a parameter characterizing stratification in the 
 lower atmosphere.
 
 # Arguments             
-- `zr`  : Instrument (reference) height (m)
+- `z`   : height (m)
 - `d`   : Zero-plane displacement height (m)
 - `MOL` : Monin-Obukhov-length L (m)
 - `df`  : DataFrame containting the variables required by [`Monin_Obukhov_length`](@ref)
 optional
 - `constants=`[`bigleaf_constants`](@ref)`()`
 
-If `MOL=nothing` in the DataFrame variants, it is computed by 
-[`Monin_Obukhov_length`](@ref) using `df` and `constants`.
+In the second variant and if `MOL=nothing` in the DataFrame variants, MOL is computed by 
+[`Monin_Obukhov_length`](@ref).
 
 # Details
 The stability parameter ``\\zeta`` is given by:
 
-``\\zeta = (zr - d) / L``
+``\\zeta = (z - d) / L``
 
-where L is the Monin-Obukhov length (m), calculated from the ction
+where L is the Monin-Obukhov length (m), calculated by
 [`Monin_Obukhov_length`](@ref). The displacement height can 
 be estimated from the function [`roughness_parameters`](@ref).
          
@@ -97,43 +95,54 @@ The mutating variant modifies or adds column [`:zeta`].
 ```jldoctest; output = false
 using DataFrames
 df = DataFrame(Tair=25, pressure=100, ustar=0.2:0.1:1.0, H=40:20:200)
-zeta = stability_parameter(df;zr=40,d=15)
+z=40;d=15
+zeta = stability_parameter.(z,d, df.Tair, df.pressure, df.ustar, df.H)
 all(zeta .< 0)
 # output
 true
 ``` 
 """ 
-function stability_parameter(zr,d,MOL)
-  zeta = (zr - d) / MOL
+function stability_parameter(z,d,MOL)
+  zeta = (z - d) / MOL
 end
-function stability_parameter!(df::AbstractDataFrame; zr,d, MOL=nothing,
+function stability_parameter(z,d,Tair, pressure, ustar, H; 
   constants=bigleaf_constants())
-  df[!,:zeta] .= stability_parameter(df; zr,d,MOL,constants)
-  df
+  MOL = Monin_Obukhov_length(Tair, pressure, ustar, H; constants);
+  stability_parameter(z,d,MOL)
 end
-function stability_parameter(df::DFTable; zr,d, MOL=nothing,
+function stability_parameter!(df::AbstractDataFrame; z,d, MOL=nothing,
   constants=bigleaf_constants())
-  if isnothing(MOL) MOL = Monin_Obukhov_length(df; constants); end
-  stability_parameter.(zr,d,MOL)
+  if isnothing(MOL)
+    ft(args...) = stability_parameter.(z,d,args...; constants)
+    transform!(df, SA[:Tair, :pressure, :ustar, :H] => ft => :zeta)
+  else
+    ft2() = stability_parameter.(z,d,MOL)
+    transform!(df, [] => ft2 => :zeta)
+  end
 end
+# function stability_parameter(df::DFTable; z,d, MOL=nothing,
+#   constants=bigleaf_constants())
+#   if isnothing(MOL) MOL = Monin_Obukhov_length(df; constants); end
+#   stability_parameter.(z,d,MOL)
+# end
 
 """
     stability_correction(zeta; 
       stab_formulation=Val(:Dyer_1970))
-    stability_correction(Tair,pressure,ustar,H, z,d; constants,
+    stability_correction(z,d, Tair,pressure,ustar,H; constants,
       stab_formulation=Val(:Dyer_1970))
-    stability_correction!(df, z, d; 
+    stability_correction!(df; zeta, z, d; 
       stab_formulation=Val(:Dyer_1970), constants = bigleaf_constants())
     
 Integrated Stability Correction Functions for Heat and Momentum
 
 # Arguments
 - `zeta`             : Stability parameter zeta (-)
-- `stab_formulation` : Formulation for the stability function. Either 
-            `Val(:Dyer_1970)`, or `Val(:Businger_1971)` or `Val(:no_stability_correction)`
 - `Tair`,`pressure`,`ustar`,`H` : see [`Monin_Obukhov_length`](@ref)
 - `z`,`d`            : see [`stability_parameter`](@ref)
 - `df`  : DataFrame containting the variables required by [`Monin_Obukhov_length`](@ref)
+- `stab_formulation` : Formulation for the stability function. Either 
+            `Val(:Dyer_1970)`, or `Val(:Businger_1971)` or `Val(:no_stability_correction)`
 
 In the second and third form computes `zeta` by [`stability_parameter`](@ref) and
 [`Monin_Obukhov_length`](@ref) and requires respective arguments.
@@ -219,8 +228,7 @@ function get_stability_coefs_unstable(::Val{:Dyer_1970}, zeta)
   (;y_h, y_m)
 end
 
-#TODO z or zr here?
-function stability_correction(Tair::Union{Missing,Number},pressure,ustar,H, z,d; 
+function stability_correction(z,d, Tair::Union{Missing,Number},pressure,ustar,H; 
   stab_formulation=Val(:Dyer_1970), constants = bigleaf_constants())
   stab_formulation isa Val{:no_stability_correction} && return(
     (psi_h = zero(z), psi_m = zero(z)))
@@ -228,16 +236,15 @@ function stability_correction(Tair::Union{Missing,Number},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, 
+# 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=nothing, z=nothing, d=nothing, 
   stab_formulation=Val(:Dyer_1970), constants = bigleaf_constants())
   # cannot dispatch on keyword argument, hence need if-clause
   if stab_formulation isa Val{:no_stability_correction}
@@ -245,12 +252,30 @@ function stability_correction!(df; zeta=df.zeta,
     df[!,:psi_m] .= 0.0
     return(df)
   end
-  ft() = stability_correction.(zeta; stab_formulation)
-  transform!(df, [] => ft => AsTable)
+  if !isnothing(zeta)
+    ft() = stability_correction.(zeta; stab_formulation)
+    transform!(df, [] => ft => AsTable)
+  else
+    function ft_met(args...)
+      stability_correction.(z, d, args...; stab_formulation, constants)
+    end
+    transform!(df, SA[:Tair,:pressure,:ustar,:H] => ft_met => AsTable)
+  end
 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] => ft => AsTable)
+function stability_correction(df::DFTable; z, d, 
+  stab_formulation=Val(:Dyer_1970), constants = bigleaf_constants()
+  )
+  # cannot dispatch on keyword argument, hence need if-clause
+  if stab_formulation isa Val{:no_stability_correction}
+    z = zero(first(skipmissing(df.ustar)))
+    rows = map(Tables.rows(df)) do row
+      (psi_h = z, psi_m = z)
+    end
+    Tables.columns(rows)
+  end
+  Tables.columns(stability_correction.(
+    z, d, df.Tair, df.pressure, df.ustar, df.H; stab_formulation, constants))
 end
+
+
diff --git a/src/surface_roughness.jl b/src/surface_roughness.jl
index bebdd10..13c9082 100755
--- a/src/surface_roughness.jl
+++ b/src/surface_roughness.jl
@@ -37,12 +37,17 @@ end
 
 
 """
-    roughness_parameters(::Val{:canopy_height},     zh; frac_d=0.7, frac_z0m=0.1)
+    roughness_parameters(::Val{:canopy_height}    , zh; frac_d=0.7, frac_z0m=0.1)
     roughness_parameters(::Val{:canopy_height_LAI}, zh, LAI; cd=0.2, hs=0.01)
-    roughness_parameters(::Val{:wind_profile},      df, zh, zr;
-      d = 0.7*zh, psi_m = nothing, constants=bigleaf_constants())
+    roughness_parameters(::Val{:wind_profile}     , ustar, wind, psi_m; 
+      zh, zr, d = 0.7*zh, constants)
 
-A simple approximation of the two roughness parameters displacement height (d)
+    roughness_parameters(method::Val{:wind_profile}, utar, wind, Tair, pressure, H; 
+      zh, zr, d = 0.7*zh, stab_formulation=Val(:Dyer_1970), constants)
+    roughness_parameters(method::Val{:wind_profile}, df::DFTable; ...)
+    
+
+A approximations of the two roughness parameters displacement height (d)
 and roughness length for momentum (z0m).
              
 # Arguments              
@@ -59,16 +64,17 @@ By canopy height and LAI
 - `hs`        : roughness length of the soil surface (m). 
 
 By wind profile
-- `df`        : DataFrame or matrix containing all required variables
-  - `wind`      : Wind speed at height zr (m s-1)
-  - `ustar`     : Friction velocity (m s-1)
-  - further variables required by [`stability_correction`](@ref) if `psi_m` is not given
+- `wind`      : Wind speed at height zr (m s-1)
+- `ustar`     : Friction velocity (m s-1)
 - `zr`        : Instrument (reference) height (m)
-- `d`         : Zero-plane displacement height (m)
+- `d`         : Zero-plane displacement height (-)
 - `psi_m`     : value of the stability function for heat, see [`stability_correction`](@ref)
-  Pass `psi_m = 0.0` to neglect stability correction.
-- `z0m`       : Roughness length for momentum (m)
-                 
+
+Another variant estimates of `psi_m` by [`stability_correction`](@ref), which
+requires further input arguments.
+For convenience, these arguments can be provided using a DataFrame, however, the result
+then is not type stable.
+
 # Details
 The two main roughness parameters, the displacement height (d)
 and the roughness length for momentum (z0m) can be estimated from simple
@@ -126,10 +132,10 @@ roughness_parameters(Val(:canopy_height_LAI),zh,2)
    
 # fix d to 0.7*zh and estimate z0m from the wind profile
 df = DataFrame(Tair=[25,25,25],pressure=100,wind=[3,4,5],ustar=[0.5,0.6,0.65],H=200)
-roughness_parameters(Val(:wind_profile),df,zh,40;d=0.7*zh)
+roughness_parameters(Val(:wind_profile),df;zh,zr=40,d=0.7*zh)
 
 # assume d = 0.8*zh
-rp = roughness_parameters(Val(:wind_profile),df,zh,40;d=0.8*zh)
+rp = roughness_parameters(Val(:wind_profile),df;zh,zr=40,d=0.8*zh)
 ≈(rp.z0m, 0.55, rtol=0.1)
 # output
 true
@@ -151,16 +157,10 @@ function roughness_parameters(::Val{:canopy_height_LAI}, zh, LAI;
   (;d, z0m, z0m_se)
 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::DFTable, zh, zr;
-  d = 0.7*zh, psi_m = nothing, stab_formulation=Val(:Dyer_1970), 
-  constants=bigleaf_constants()
+function roughness_parameters(::Val{:wind_profile}, ustar::AbstractVector, wind, psi_m; 
+  zh, zr, d = 0.7*zh, constants=bigleaf_constants()
   )
-  if isnothing(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 = allowmissing(@. (zr - d) * exp(-constants[:k]*wind / ustar - psi_m))
   #z0m_all[(z0m_all .> zh)] .= missing # problems with missings
   replace!(x -> !ismissing(x) && x > zh ? missing : x, z0m_all)
   nval = sum(.!ismissing.(z0m_all))
@@ -168,6 +168,35 @@ function roughness_parameters(::Val{:wind_profile}, df::DFTable, zh, zr;
   z0m_se = constants[:se_median] * (std(skipmissing(z0m_all)) / sqrt(nval))
   (;d, z0m, z0m_se)
 end
+
+function roughness_parameters(method::Val{:wind_profile}, 
+  ustar::AbstractVector, wind, Tair, pressure, H; zh, zr, d = 0.7*zh,
+  stab_formulation=Val(:Dyer_1970), constants=bigleaf_constants(), kwargs...
+  )
+  # psi_m = Tables.Columns(stability_correction.(
+  #   zr, d, Tair, pressure, ustar, H; stab_formulation, constants)).psi_m
+  psi_m = getindex.(stability_correction.(
+    zr, d, Tair, pressure, ustar, H; stab_formulation, constants), :psi_m)
+  roughness_parameters(method, ustar, wind, psi_m; zh, zr, d, constants, kwargs...)
+end
+
+function roughness_parameters(method::Val{:wind_profile}, df::DFTable; 
+  psi_m = nothing, stab_formulation=Val(:Dyer_1970), kwargs...
+  )
+  !isnothing(psi_m) && return(roughness_parameters(
+    method, df.ustar, df.wind, psi_m; kwargs...))
+  if stab_formulation isa Val{:no_stability_correction}
+    psi_m = 0.0
+    roughness_parameters(method, df.ustar, df.wind, psi_m; kwargs...)
+  else
+    roughness_parameters(method, df.ustar, df.wind, df.Tair, df.pressure, df.H; 
+      stab_formulation, kwargs...)
+  end
+end
+
+
+
+
   
 """                                                                                                                       
     wind_profile(z::Number, ustar, d, z0m, psi_m = zero(z); constants)
@@ -237,52 +266,43 @@ using DataFrames
 heights = 18:2:40  # heights above ground for which to calculate wind speed
 df = DataFrame(Tair=25,pressure=100,wind=[3,4,5],ustar=[0.5,0.6,0.65],H=[200,230,250]) 
 zr=40;zh=25;d=16
-psi_m = stability_correction!(copy(df, copycols=false), zr, d).psi_m
-z0m = roughness_parameters(Val(:wind_profile), df, zh, zr; psi_m).z0m
+z0m = roughness_parameters(Val(:wind_profile), df; zh, zr).z0m
 ws = map(heights) do z
-  wind_profile(df,z,d,z0m,psi_m)
+  wind_profile(z,df,d,z0m; zr, zh)
 end
 using Plots # plot wind profiles for the three rows in df
-plot(heights, first.(ws), xlab = "height (m)", ylab = "wind speed (m/s)")
-plot!(heights, getindex.(ws, 2))
-plot!(heights, getindex.(ws, 3))
+plot(first.(ws), heights, ylab = "height (m)", xlab = "wind speed (m/s)", legend=:topleft)
+plot!(getindex.(ws, 2), heights)
+plot!(getindex.(ws, 3), heights)
 nothing
 # output
 ``` 
 """  
-function wind_profile(z::Number, ustar, d, z0m, psi_m = zero(z); constants=bigleaf_constants())
+function wind_profile(z::Number, ustar, d, z0m, psi_m; constants=bigleaf_constants())
   wind_heights = max(0,(ustar / constants[:k]) * (log(max(0,(z - d)) / z0m) - psi_m))
 end
 
-function wind_profile(stab_formulation, z, ustar, Tair,pressure,H, d, z0m;
-  constants=bigleaf_constants())
-  # in order to comput psi_m need Monin_Obukhov_length with additional vars
-  MOL = Monin_Obukhov_length(Tair,pressure,ustar,H; constants)
-  zeta  = stability_parameter(z,d,MOL)
-  psi_m = stability_correction(zeta; stab_formulation).psi_m
+function wind_profile(z, ustar::Union{Missing,Number}, d, z0m, Tair,pressure,H,
+  stab_formulation=Val(:Dyer_1970), constants=bigleaf_constants())
+  psi_m = stability_correction(z,d, Tair,pressure,ustar,H; stab_formulation, constants).psi_m
   wind_profile(z, ustar, d, z0m, psi_m)
 end
 
-function wind_profile(df::DFTable, z, d, z0m = nothing; 
-  zh = nothing, zr = nothing, 
+function wind_profile(z, df::DFTable, d, z0m = nothing; 
+  zh = nothing, zr = nothing, psi_m = nothing,
   stab_formulation = Val(:Dyer_1970), constants = bigleaf_constants()
   )
-  # TODO providingn z or zr to stabiity correction?
   if isnothing(z0m)
     (isnothing(zh) || isnothing(zr)) && error(
       "wind_profile: If 'z0m' is not given, must specify both 'zh' and 'zr' optional " *
       "arguments to estimate 'z0m' from 'df.wind' and 'df.ustar' measured at height zr.")
     # uses psi_m at zr
-    z0m = roughness_parameters(Val(:wind_profile), df, zh, zr).z0m
+    z0m = roughness_parameters(Val(:wind_profile), df; zh, zr).z0m
+  end
+  if isnothing(psi_m)
+    psi_m = stability_correction(df; z, d, stab_formulation, constants).psi_m
   end
-  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::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)
 end
 
diff --git a/test/aerodynamic_conductance.jl b/test/aerodynamic_conductance.jl
index dd18205..98f7b95 100644
--- a/test/aerodynamic_conductance.jl
+++ b/test/aerodynamic_conductance.jl
@@ -29,12 +29,12 @@ end
     zr, zh, d = thal.zr, thal.zh, 0.7*thal.zh
     ustar, wind = tha48[24, [:ustar, :wind]]
     df = copy(tha48)
-    stability_parameter!(df; zr, d)
-    stability_correction!(df)
-    z0m = roughness_parameters(Val(:wind_profile), df, zh, zr; psi_m = df.psi_m).z0m
+    stability_correction!(df; z=zr, d)
+    z0m = roughness_parameters(Val(:wind_profile), df; zh, zr, psi_m = df.psi_m).z0m
     Ram_r = @inferred compute_Ram(Val(:wind_zr), ustar, wind)
     @test Ram_r ≈ 6.31 rtol = 1/100 # TODO check with R
-    Ram_p =  @inferred 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]]
@@ -62,12 +62,13 @@ end
     # missing due to missing zr
     @test all(ismissing.(df.psi_h))        
     @test all(ismissing.(df.psi_m))        
-    @inferred aerodynamic_conductance!(df; Gb_model=Val(:Thom_1972), zr = thal.zr, zh = thal.zh
-        ,stab_formulation = Val(:no_stability_correction))
+    @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]
+        SA[:Gb_h, :Rb_h, :kB_h, :Gb_CO2, :Ra_m, :Ga_m, :Ga_h, :Ra_h, :Ga_CO2]
     # 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))
@@ -77,10 +78,11 @@ end
 @testset "aerodynamic_conductance! Gb_Thom" begin
     df = copy(tha48)
     # test Ga_m by wind_profile
-    @inferred aerodynamic_conductance!(df; Gb_model=Val(:Thom_1972), Ram_model=Val(: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]
+        SA[:Gb_h, :Rb_h, :kB_h, :Gb_CO2, :Ra_m, :Ga_m, :Ga_h, :Ra_h, :Ga_CO2]
 end
 
 @testset "aerodynamic_conductance! constant_kB1" begin
@@ -91,7 +93,7 @@ end
     @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]
+        SA[:Gb_h, :Rb_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
 
@@ -99,10 +101,11 @@ end
     leafwidth=0.1
     df = copy(tha48)
     # test Ga_m by wind_profile
-    @inferred aerodynamic_conductance!(df; Gb_model=Val(:Choudhury_1988), Ram_model=Val(: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]
+        SA[:Gb_h, :Rb_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))
@@ -117,7 +120,7 @@ end
     @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]
+        SA[:Gb_h, :Rb_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))
diff --git a/test/boundary_layer_conductance.jl b/test/boundary_layer_conductance.jl
index 70f8478..46b7ec6 100644
--- a/test/boundary_layer_conductance.jl
+++ b/test/boundary_layer_conductance.jl
@@ -22,76 +22,73 @@
     @test df2.Gb_N20[1] == Gb.Gb_N20
 end
 
-@testset "compute_Gb Gb_constant_kB1" begin
+@testset "compute_Gb Gb_constant_kB1 and compute_Gb_quantities" begin
     kB_h = 1.18
-    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))
+    ustar = 0.1
+    Gb = @inferred Gb_constant_kB1(ustar, kB_h)
+    @test Gb ≈ 0.0347 rtol = 1e-2
+    Gbq = compute_Gb_quantities(Gb, ustar)
+    @test keys(Gbq) == (:Rb_h, :Gb_h, :kB_h, :Gb_CO2)
+    @test all(isapprox.(values(Gbq), values(
+        (Rb_h = 28.8, Gb_h = 0.0347, kB_h = 1.18, Gb_CO2 = 0.0264)), rtol = 1e-2))
     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]))))
+    @test ismissing(Gb)
+    Gbq = compute_Gb_quantities(Gb, ustar)
+    @test keys(Gbq) == (:Rb_h, :Gb_h, :kB_h, :Gb_CO2)
+    @test all(ismissing.(values(getindex.(Ref(Gbq), SA[:Rb_h, :Gb_h, :Gb_CO2]))))
     #
     # DataFrame variant
     dfo = DataFrame(ustar = SA[0.1,missing,0.3])
     df = copy(dfo)
     @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))
+    @test propertynames(df) == [:ustar, :Gb_h]
+    compute_Gb_quantities!(df)
+    @test propertynames(df) == [:ustar, :Gb_h, :Rb_h, :kB_h, :Gb_CO2]
+    @test all(isapprox.(values(df[1,2:end]), values(
+        (Gb_h = 0.0347, Rb_h = 28.8, kB_h = 1.18, Gb_CO2 = 0.0264)), rtol = 1e-2))
 end
 
 @testset "compute_Gb Gb_Thom" begin
-    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)
-    @test keys(Gb) == (:Rb_h, :Gb_h, :kB_h, :Gb_CO2)
-    @test all(ismissing.(values(Gb)))
+    ustar = 0.1
+    Gb1 = @inferred Gb_Thom(ustar)
+    @test Gb1 ≈ 0.0347 rtol = 1e-2
+    Gbm = Gb_Thom(missing)
+    @test ismissing(Gbm)
     #
     # DataFrame variant
-    dfo = DataFrame(ustar = SA[0.1,missing,0.3])
+    dfo = DataFrame(ustar = SA[ustar,missing,0.3])
     df = copy(dfo)
     @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))
+    @test propertynames(df) == [:ustar, :Gb_h]
+    @test df.Gb_h[1] == Gb1
 end
 
-@testset "compute_Gb Gb_Choudhury" begin
+@testset "compute_Gb Gb_Choudhury and Gb_Su" begin
     zh, zr, LAI = thal.zh, thal.zr, thal.LAI
     leafwidth=0.1
-    wind_zh = 1.2
-    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 = @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
+    wind_zh = 2.1662688
+    Gb_Choud = @inferred Gb_Choudhury(; leafwidth, LAI, wind_zh)
+    @test Gb_Choud ≈ 0.157 rtol=1e-2 # from R
+    Gbm = @inferred Gb_Choudhury(; leafwidth, LAI, wind_zh=missing)
+    @test ismissing(Gbm)
+    #
+    Tair, pressure, ustar = values(tha48[1, SA[:Tair, :pressure, :ustar]])
+    Dl=0.01; fc = (1-exp(-LAI/2)) 
+    Gb_S = @inferred Gb_Su(Tair, pressure, ustar; wind_zh, Dl, fc)
+    @test Gb_S ≈ 0.185 rtol=1e-2 # from R
+    Gbm = @inferred Gb_Su(missing, pressure, ustar; wind_zh, Dl, fc)
+    @test ismissing(Gbm)
     #
     # DataFrame variant
     df = copy(tha48)
-    @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))
-end
-
-@testset "compute_Gb Gb_Su" begin
-    zh, zr, LAI = thal.zh, thal.zr, thal.LAI
-    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 = @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 = @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)))
+    wind_zh = wind_profile(zh, df, 0.7*zh; zh, zr)
+    @inferred compute_Gb!(df, Val(:Choudhury_1988); leafwidth, LAI, wind_zh)
+    @test last(propertynames(df)) == :Gb_h
+    @test df.Gb_h[1] ≈ Gb_Choud rtol=1e-6
     #
-    # DataFrame variant
     df = copy(tha48)
-    @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
+    @inferred compute_Gb!(df, Val(:Su_2001); wind_zh, Dl, LAI)
+    @test last(propertynames(df)) == :Gb_h
+    @test df.Gb_h[1] ≈ Gb_S rtol=1e-6
 end
 
diff --git a/test/stability_correction.jl b/test/stability_correction.jl
index 53dcecf..3bff533 100644
--- a/test/stability_correction.jl
+++ b/test/stability_correction.jl
@@ -3,46 +3,44 @@
     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)
     dfd = disallowmissing(df)
     dfm = allowmissing(df)
     @inferred Monin_Obukhov_length!(df)
     @test df.MOL[24] == MOL24
-    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 = @inferred stability_parameter(zr,d,MOL24)
+    df24 = df[24,:]
+    MOL24 = Monin_Obukhov_length(df24.Tair, df24.pressure, df24.ustar, df24.H)
+    z=40;d=15
+    zeta = @inferred stability_parameter(z,d,MOL24)
     @test ≈(zeta , -0.240, rtol = 1/100)
+    zeta2 = @inferred stability_parameter(z,d,df24.Tair, df24.pressure, df24.ustar, df24.H)
+    @test zeta2 == zeta
     #
-    @inferred stability_parameter!(df;zr,d)
+    @inferred stability_parameter!(df;z,d)
     @test df.zeta[24] == zeta
     zeta_scalar = df.zeta
     # 
-    # non-mutatingvariant
     # also test zr as a vector
     df = copy(tha48)
-    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
+    df[!,:zi] .= z
+    df.zi[1] = z/2
+    zetas = stability_parameter!(df; z=df.zi, d).zeta
     @test zetas[2:24] == zeta_scalar[2:24]
     @test zetas[1] != zeta_scalar[1]
+    #
+    # specifyngn MOL directly
+    df4 = copy(tha48)
+    Monin_Obukhov_length!(df4) 
+    zetas2 = stability_parameter!(df4; z, d, MOL = df4.MOL./2).zeta
+    @test all(zetas2[2:end] .== zetas[2:end] .* 2)
 end
 
-@testset "stability_correction" begin
+@testset "stability_correction zeta" begin
     zeta = -2:0.5:0.5
     @inferred stability_correction(first(zeta))
     df2 = DataFrame(stability_correction.(zeta))
@@ -61,51 +59,50 @@ end
     @test resm == (psi_h = 0.0, psi_m = 0.0)
 end
 
+@testset "stability_correction metvars" begin
+    z=40;d=15
+    df = DataFrame(Tair=25, pressure=100, ustar=0.2:0.1:1.0, H=40:20:200)
+    df1 = df[1,:]
+    zeta1 = stability_parameter(z,d,df1.Tair, df1.pressure, df1.ustar, df1.H)
+    res1 = @inferred stability_correction(z,d, df1.Tair, df1.pressure, df1.ustar, df1.H) 
+    @test res1 == stability_correction(zeta1)
+end
+
 @testset "stability_correction DataFrame variant" begin
-    zr=40;d=15
+    z=40;d=15
     dfo = DataFrame(Tair=25, pressure=100, ustar=0.2:0.1:1.0, H=40:20:200)
     df = copy(dfo)
+    res = stability_correction(df; z, d)
     # for type stability, use columntable(df)
-    tmp = @inferred stability_correction(columntable(df), zr, d)
-    @inferred stability_correction!(df, zr, d)
+    stability_correction!(df; z, d)
     propertynames(df)[(end-1):end] == SA[:psi_h, :psi_m]
+    @test DataFrame(res) == df[!, (end-1):end]
     #
-    dfm = copy(dfo)
-    @inferred stability_correction!(dfm, zr, d; stab_formulation = Val(:no_stability_correction))
+    dfm = allowmissing(dfo)
+    dfm.ustar[1] = missing
+    res = stability_correction(
+        dfm; z, d, stab_formulation = Val(:no_stability_correction))
+    stability_correction!(
+        dfm; z, 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 DataFrame(res) == dfm[!, (end-1):end]
     #
-    # test computing zeta first
+    # test specifying zeta instead of z and d
     df2 = copy(dfo)
-    stability_parameter!(df2; zr, d) # adds zeta
-    stability_correction!(df2)
+    stability_parameter!(df2; z, d) # adds zeta
+    stability_correction!(df2, zeta=df2.zeta)
     @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)
+    df3[!,:zi] .= z
+    df3.zi[1] = z/2
+    res = stability_correction(df3; z=df3.zi, d)
+    @inferred stability_correction!(df3; z=df3.zi, 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 = @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 = @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
-    @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)))
+    @test DataFrame(res) == df3[!, (end-1):end]
 end
diff --git a/test/surface_roughness.jl b/test/surface_roughness.jl
index 268902b..6f8fcc6 100644
--- a/test/surface_roughness.jl
+++ b/test/surface_roughness.jl
@@ -20,32 +20,45 @@ end
     @test all(isapproxm.(values(rp), (21.77, 1.419, missing), rtol=1e-3))
     #
     df = copy(tha48)
-    dfd = disallowmissing(df)
-    #df.wind[1] = missing
+    #df.wind[1] = missing # obscures comparing to R
     d=0.7*zh
-    psi_m = (@inferred stability_correction(columntable(dfd), zr, d)).psi_m
-    psi_m = stability_correction(df, zr, d).psi_m
+    psi_m = stability_correction!(df; z=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)
+    # not type-stable if some columns allow missings in Julia 1.6
+    #rp = @inferred roughness_parameters(Val(:wind_profile), df.ustar, df.wind, psi_m; zh, zr)
+    rp = roughness_parameters(Val(:wind_profile), df.ustar, df.wind, psi_m; zh, zr)
     @test keys(rp) == keys_exp
     #@test all(isapproxm.(values(rp), (18.55, 1.879, 0.3561), rtol=1e-3))
     #from R:
     @test all(isapproxm.(values(rp), (18.55, 1.879402, 0.356108), rtol=1e-3))
     #
     # no stability correction
-    rp0 = @inferred roughness_parameters(Val(:wind_profile), columntable(dfd), zh, zr; 
+    # broadcast across Missings is not type-stable 
+    # rp0 = @inferred roughness_parameters(Val(:wind_profile), df.ustar, df.wind, 
+    #     df.Tair, df.pressure, df.H; zh, zr, 
+    #     stab_formulation = Val(:no_stability_correction))
+    dfd = disallowmissing(df)
+    rp0 = @inferred roughness_parameters(Val(:wind_profile), dfd.ustar, dfd.wind, 
+        dfd.Tair, dfd.pressure, dfd.H; 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))
+    # providing DataFrame is not type stable
+    rp0b = roughness_parameters(Val(:wind_profile), df; zh, zr, 
+        stab_formulation = Val(:no_stability_correction))
+    @test rp0b == rp0         
     #
     # estimate psi
     #@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
+    rp_psiauto = roughness_parameters(Val(:wind_profile), df; zh, zr)
     @test rp_psiauto == rp
+    #
+    # DataFrame with only columns ustar and wind
+    rp0c = roughness_parameters(Val(:wind_profile), df[!,Cols(:ustar, :wind)]; zh, zr, 
+        stab_formulation = Val(:no_stability_correction))
+    @test rp0c == rp0         
 end
 
 @testset "wind_profile" begin
@@ -53,33 +66,33 @@ end
     z = 30
     d=0.7*thal.zh
     z0m= 2.65
-    u30 = @inferred wind_profile(z, ustar, d, z0m)
+    u30 = @inferred wind_profile(z, ustar, d, z0m, 0.0)
     @test ≈(u30, 1.93, rtol = 1/100 ) # from R
     #
-    u30c = @inferred wind_profile(Val(:Dyer_1970), z, ustar, Tair,pressure, H, d, z0m)
+    u30c = @inferred wind_profile(z, ustar, d, z0m, Tair,pressure, H)
     @test ≈(u30c, 2.31, rtol = 1/100 ) # from R
     #
     z0m=1.9 #2.14 #2.65
-    u30 = @inferred wind_profile(z, ustar, d, z0m) # used below
-    u30c = @inferred wind_profile(Val(:Dyer_1970), z, ustar, Tair,pressure, H, d, z0m)
+    u30 = @inferred wind_profile(z, ustar, d, z0m, 0.0) # used below
+    u30c = @inferred wind_profile(z, ustar, d, z0m, Tair,pressure, H)
     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))    
+    windz = @inferred wind_profile(z, columntable(dfd), d, z0m; stab_formulation = Val(:no_stability_correction))    
+    windz = wind_profile(z, df, d, z0m; stab_formulation = Val(:no_stability_correction))    
     @test length(windz) == 48
     @test windz[1] == u30
-    windzc = @inferred wind_profile(columntable(dfd), z, d, z0m; stab_formulation = Val(:Dyer_1970))    
+    windzc = @inferred wind_profile(z, columntable(dfd), d, z0m; stab_formulation = Val(:Dyer_1970))    
     @test windzc[1] == u30c
     #plot(windz)
     #plot!(windz2)
-    psi_m = stability_correction(df, z, d).psi_m
-    windzc2 = @inferred wind_profile(columntable(dfd), z, d, z0m, psi_m)    
+    psi_m = stability_correction(df; z, d).psi_m
+    windzc2 = @inferred wind_profile(z, columntable(dfd), 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 = @inferred wind_profile(columntable(dfd), z, d; zh=thal.zh, zr=thal.zr, 
+    @test_throws Exception wind_profile(z, df, d)    
+    windzc3 = @inferred wind_profile(z, columntable(dfd), d; zh=thal.zh, zr=thal.zr, 
         stab_formulation = Val(:Dyer_1970))    
     # may have used slightly different estimated z0m
     #windzc3 - windzc