diff --git a/README.md b/README.md
index 907046f..57945cd 100644
--- a/README.md
+++ b/README.md
@@ -46,6 +46,8 @@ file an [issue](https://github.com/bgctw/Bigleaf.jl/issues) if you need a specif
 
 At the current state we ported
 - Meteorological variables
+- Boundary layer and Aerodynamic conductance
+- Surface conductance
 - Evapotranspiration
 - Potential radiation
 - Unit conversions
diff --git a/docs/src/walkthrough.md b/docs/src/walkthrough.md
index ac95705..d5a22ca 100644
--- a/docs/src/walkthrough.md
+++ b/docs/src/walkthrough.md
@@ -295,15 +295,15 @@ The following figure compares them at absole scale and as difference to the
 #Tair = 0:0.25:12
 ##Tair = [10.0,20.0]
 #eform_def = Val(:Sonntag_1990)
-#Esat_def = Esat_from_Tair.(Tair; formula = eform_def)
+#Esat_def = Esat_from_Tair.(Tair; Esat_formula = eform_def)
 #eforms = (Val(:Sonntag_1990), Val(:Alduchov_1996), Val(:Allen_1998))
 #eform = eforms[2]
 #string.(eforms)
 #df = mapreduce(vcat, eforms) do eform 
-#    Esat = Esat_from_Tair.(Tair; formula = eform)
+#    Esat = Esat_from_Tair.(Tair; Esat_formula = eform)
 #    local dff # make sure to not override previous results
 #    dff = DataFrame(
-#        formula = eform, Tair = Tair, 
+#        Esat_formula = eform, Tair = Tair, 
 #        Esat = Esat,
 #        dEsat = Esat - Esat_def,
 #        )
@@ -311,8 +311,8 @@ The following figure compares them at absole scale and as difference to the
 ##using Chain
 #using Pipe
 #using Plots, StatsPlots
-#dfw = @pipe df |> select(_, 1,2, :Esat) |> unstack(_, :formula, 3)
-#dfws = @pipe df |> select(_, 1,2, :dEsat) |> unstack(_, :formula, 3)
+#dfw = @pipe df |> select(_, 1,2, :Esat) |> unstack(_, :Esat_formula, 3)
+#dfws = @pipe df |> select(_, 1,2, :dEsat) |> unstack(_, :Esat_formula, 3)
 #@df dfw plot(:Tair, cols(2:4), legend = :topleft, xlab="Tair (degC)", #ylab="Esat (kPa)")
 #savefig("Esat_abs.svg")
 #@df dfws plot(:Tair, cols(2:4), legend = :topleft, xlab="Tair (degC)", #ylab="Esat -ESat_Sonntag_1990 (kPa)")
@@ -388,6 +388,32 @@ add_Gb!(thas, :Gb_O2 => 0.84, :Gb_CH4 => 0.99); # adds Gb_O2 and Gb_CH4
 select(first(thas,3), r"Gb_")
 ```
 
+## Surface conductance
+
+Knowledge of aerodynamic conductance $G_a$ 
+allows us to calculate the bulk surface conductance ($G_s$) of the site 
+(In this case by inverting the Penman-Monteith equation). Gs represents the combined 
+conductance of the vegetation and the soil to water vapor transfer (and as such it is not 
+a purely physiological quantity). Calculating $G_s$ in `Bigleaf.jl` is simple:
+
+```@example doc
+surface_conductance!(thas, Val(:PenmanMonteith));
+thas[1:3,Cols(:datetime, r"Gs")]
+```
+
+The two columns only differ in the unit of $G_s$. 
+One in m s$^{-1}$ and one in mol m$^{-2}$ s$^{-1}$. 
+In this function we have ignored the ground heat flux ($G$) and the storage fluxes ($S$).
+By default they are assumed zero.
+In our example we do not have information on the storage fluxes, but we have measurements 
+on the ground heat flux, which we should add to the function call:
+
+```@example doc
+surface_conductance!(thas, Val(:PenmanMonteith); G=thas.G);
+thas[1:3,Cols(:datetime, r"Gs")]
+```
+
+
 ## Wind profile
 
 The 'big-leaf' framework assumes that wind speed is zero at height d + $z_{0m}$ 
diff --git a/inst/fromR/decoupling.jl b/inst/fromR/decoupling.jl
index 618c39e..6ddf1e5 100755
--- a/inst/fromR/decoupling.jl
+++ b/inst/fromR/decoupling.jl
@@ -14,7 +14,7 @@
 #' - approach    Approach used to calculate omega. Either `"Jarvis&McNaughton_1986"` (default)
 #'                    or `"Martin_1989"`.
 #' - LAI         Leaf area index (m2 m-2), only used if `approach = "Martin_1989"`.
-#' - Esat_formula  Optional: formula to be used for the calculation of esat and the slope of esat.
+#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
 #'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
 #'                      See [`Esat_slope`](@ref). 
 #' - constants   Kelvin - conversion degree Celsius to Kelvin 
@@ -121,7 +121,7 @@ end
 #'                  cp - specific heat of air for constant pressure (J K-1 kg-1)
 #'                  
 #' # Details
- the following formula is used (Martin, 1989):
+ the following Esat_formula is used (Martin, 1989):
 #' 
 #'          ``Gr = 4 \sigma Tair^3 LAI / cp``                             
 #'                                       
diff --git a/inst/fromR/evapotranspiration.jl b/inst/fromR/evapotranspiration.jl
index 9758871..ec2b1da 100755
--- a/inst/fromR/evapotranspiration.jl
+++ b/inst/fromR/evapotranspiration.jl
@@ -21,7 +21,7 @@
 #'                  only used if `approach = Val(:PenmanMonteith)`.
 #' - missing_G_as_NA  if `true`, missing G are treated as `missing`s, otherwise set to 0. 
 #' - missing_S_as_NA  if `true`, missing S are treated as `missing`s, otherwise set to 0. 
-#' - Esat_formula  Optional: formula to be used for the calculation of esat and the slope of esat.
+#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
 #'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
 #'                      See [`Esat_slope`](@ref). 
 #' - constants cp - specific heat of air for constant pressure (J K-1 kg-1) 
@@ -165,7 +165,7 @@ end
 #' - S         Sum of all storage fluxes (W m-2); optional
 #' - missing_G_as_NA  if `true`, missing G are treated as `missing`s, otherwise set to 0. 
 #' - missing_S_as_NA  if `true`, missing S are treated as `missing`s, otherwise set to 0. 
-#' - Esat_formula  Optional: formula to be used for the calculation of esat and the slope of esat.
+#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
 #'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
 #'                      See [`Esat_slope`](@ref). 
 #' - constants cp - specific heat of air for constant pressure (J K-1 kg-1) 
@@ -202,7 +202,7 @@ end
 #' - S         Sum of all storage fluxes (W m-2); optional
 #' - missing_G_as_NA  if `true`, missing G are treated as `missing`s, otherwise set to 0. 
 #' - missing_S_as_NA  if `true`, missing S are treated as `missing`s, otherwise set to 0.
-#' - Esat_formula  Optional: formula to be used for the calculation of esat and the slope of esat.
+#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
 #'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
 #'                      See [`Esat_slope`](@ref). 
 #' - constants cp - specific heat of air for constant pressure (J K-1 kg-1) 
diff --git a/inst/fromR/meteorological_variables.jl b/inst/fromR/meteorological_variables.jl
index 924002c..aaea875 100755
--- a/inst/fromR/meteorological_variables.jl
+++ b/inst/fromR/meteorological_variables.jl
@@ -111,7 +111,7 @@ end
 #'              corresponding slope of the saturation vapor pressure curve.
 #' 
 #' - Tair      Air temperature (deg C)
-#' - formula   Formula to be used. Either `"Sonntag_1990"` (Default), 
+#' - Esat_formula   Esat_formula to be used. Either `"Sonntag_1990"` (Default), 
 #'                  `"Alduchov_1996"`, or `"Allen_1998"`.
 #' - constants Pa2kPa - conversion pascal (Pa) to kilopascal (kPa)
 #' 
@@ -160,20 +160,20 @@ end
 #' @importFrom stats D                  
 """
 """
-function Esat_slope(Tair,formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
+function Esat_slope(Tair,Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                        constants=bigleaf_constants())
   
-  formula = match_arg(formula)
+  Esat_formula = match_arg(Esat_formula)
   
-  if (formula == "Sonntag_1990")
+  if (Esat_formula == "Sonntag_1990")
     a = 611.2
     b = 17.62
     c = 243.12
-elseif (formula == "Alduchov_1996")
+elseif (Esat_formula == "Alduchov_1996")
     a = 610.94
     b = 17.625
     c = 243.04
-elseif (formula == "Allen_1998")
+elseif (Esat_formula == "Allen_1998")
     a = 610.8
     b = 17.27
     c = 237.3
@@ -239,7 +239,7 @@ end
 #' - Tair  Air temperature (deg C)
 #' 
 #' # Details
- The following formula is used:
+ The following Esat_formula is used:
 #' 
 #'   ``\\lambda = (2.501 - 0.00237*Tair)10^6``
 #' 
@@ -280,7 +280,7 @@ end
 #' - Tair         Air temperature (degC)
 #' - gamma        Psychrometric constant (kPa K-1)
 #' - accuracy     Accuracy of the result (degC)
-#' - Esat_formula Optional: formula to be used for the calculation of esat and the slope of esat.
+#' - Esat_formula Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
 #'                     One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
 #'                     See [`Esat_slope`](@ref). 
 #' - constants    Pa2kPa - conversion pascal (Pa) to kilopascal (kPa)
@@ -308,7 +308,7 @@ end
 #' - pressure  Atmospheric pressure (kPa)
 #' - VPD       Vapor pressure deficit (kPa)
 #' - accuracy  Accuracy of the result (deg C)
-#' - Esat_formula  Optional: formula to be used for the calculation of esat and the slope of esat.
+#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
 #'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
 #'                      See [`Esat_slope`](@ref). 
 #' - constants cp - specific heat of air for constant pressure (J K-1 kg-1) 
@@ -381,7 +381,7 @@ end
 #' 
 #' - ea           Air vapor pressure (kPa)
 #' - accuracy     Accuracy of the result (degC)
-#' - Esat_formula Optional: formula to be used for the calculation of esat and the slope of esat.
+#' - Esat_formula Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
 #'                     One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
 #'                     See [`Esat_slope`](@ref). 
 #' - constants    Pa2kPa - conversion pascal (Pa) to kilopascal (kPa)
@@ -409,7 +409,7 @@ end
 #' - Tair     Air temperature (degC)
 #' - VPD      Vapor pressure deficit (kPa)
 #' - accuracy Accuracy of the result (deg C)
-#' - Esat_formula  Optional: formula to be used for the calculation of esat and the slope of esat.
+#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
 #'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
 #'                      See [`Esat_slope`](@ref). 
 #' - constants Pa2kPa - conversion pascal (Pa) to kilopascal (kPa)
@@ -472,7 +472,7 @@ end
 #' - Tair      Air temperature (deg C)
 #' - pressure  Atmospheric pressure (kPa)
 #' - VPD       Vapor pressure deficit (kPa)
-#' - Esat_formula  Optional: formula to be used for the calculation of esat and the slope of esat. 
+#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat. 
 #'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
 #'                      See [`Esat_slope`](@ref). 
 #' - constants Kelvin - conversion degree Celsius to Kelvin 
diff --git a/inst/fromR/surface_conditions.jl b/inst/fromR/surface_conditions.jl
index ec08643..0a992ac 100755
--- a/inst/fromR/surface_conditions.jl
+++ b/inst/fromR/surface_conditions.jl
@@ -19,7 +19,7 @@
 #' - Ca               Atmospheric CO2 concentration (mol mol-1). Required if `calc_surface_CO2 = true`.
 #' - NEE              Net ecosystem exchange (umol m-2 s-1). Required if `calc_surface_CO2 = true`.
 #' - Ga_CO2           Aerodynamic conductance for CO2 (m s-1). Required if `calc_surface_CO2 = true`.
-#' - Esat_formula     Optional: formula to be used for the calculation of esat and the slope of esat.
+#' - Esat_formula     Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
 #'                         One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
 #'                         See [`Esat_slope`](@ref).           
 #' - constants        cp - specific heat of air for constant pressure (J K-1 kg-1)  
diff --git a/inst/fromR/surface_conductance.jl b/inst/fromR/surface_conductance.jl
index 7e95d7c..b51d50d 100755
--- a/inst/fromR/surface_conductance.jl
+++ b/inst/fromR/surface_conductance.jl
@@ -21,9 +21,9 @@
 #' - missing_S_as_NA  if `true`, missing S are treated as `missing`s, otherwise they are set to 0. 
 #'                          Only used if `formulation = Val(:PenmanMonteith)`.
 #' - formulation Formulation used. Either `Val(:PenmanMonteith)` (the default) 
-#'                    using the inverted Penman-Monteith equation, or `"Flux-Gradient"`,
+#'                    using the inverted Penman-Monteith equation, or `Val(:FluxGradient)`,
 #'                    for a simple flux-gradient approach requiring ET, pressure, and VPD only.
-#' - Esat_formula  Optional: formula to be used for the calculation of esat and the slope of esat.
+#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
 #'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`. 
 #'                      Only used if `formulation = Val(:PenmanMonteith)`. See [`Esat_slope`](@ref).
 #' - constants   cp - specific heat of air for constant pressure (J K-1 kg-1) 
@@ -48,7 +48,7 @@
 #'  By default, any missing data in G and S are set to 0. If `missing_S_as_NA = true`
 #'  or `missing_S_as_NA = true`, Gs will give `missing` for these timesteps.
 #'  
-#'  If `formulation="Flux-Gradient"`, Gs (in mol m-2 s-1) is calculated from VPD and ET only:
+#'  If `formulation=Val(:FluxGradient)`, Gs (in mol m-2 s-1) is calculated from VPD and ET only:
 #'  
 #'     ``Gs = ET/pressure * VPD``
 #'  
@@ -84,7 +84,7 @@
 #' 
 #' # calculate Gs based on a simple gradient approach
 #' Gs_gradient = surface_conductance(DE_Tha_Jun_2014_2,Tair="Tair",pressure="pressure",
-#'                                    VPD="VPD",formulation="Flux-Gradient")
+#'                                    VPD="VPD",formulation=Val(:FluxGradient))
 #' summary(Gs_gradient)
 #' 
 #' # calculate Gs from the the inverted PM equation (now Rn, and Ga are needed),
@@ -120,13 +120,13 @@
 """
 function surface_conductance(data,Tair="Tair",pressure="pressure",Rn="Rn",G=nothing,S=nothing,
                                 VPD="VPD",LE="LE",Ga="Ga_h",missing_G_as_NA=false,missing_S_as_NA=false,
-                                formulation=c(Val(:PenmanMonteith),"Flux-Gradient"),
+                                formulation=c(Val(:PenmanMonteith),Val(:FluxGradient)),
                                 Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                                 constants=bigleaf_constants())
   
   formulation = match_arg(formulation)
   
-  if (formulation == "Flux-Gradient")
+  if (formulation == Val(:FluxGradient))
   
     check_input(data,list(Tair,pressure,VPD,LE))
     
diff --git a/inst/fromR/unit_conversions.jl b/inst/fromR/unit_conversions.jl
index 7f38a59..b65a156 100755
--- a/inst/fromR/unit_conversions.jl
+++ b/inst/fromR/unit_conversions.jl
@@ -125,7 +125,7 @@ end
 #' - q         Specific humidity (kg kg-1)
 #' - VPD       Vapor pressure deficit (kPa)
 #' - rH        Relative humidity (-)
-#' - Esat_formula  Optional: formula to be used for the calculation of esat and the slope of esat.
+#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
 #'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
 #'                      See [`Esat_slope`](@ref).
 #' - constants eps - ratio of the molecular weight of water vapor to dry air (-) 
diff --git a/inst/tha.jl b/inst/tha.jl
index df3e419..f2efd1f 100644
--- a/inst/tha.jl
+++ b/inst/tha.jl
@@ -38,13 +38,20 @@ thas = subset(thaf, :valid)
 
 
 function tmpf()
-    
+    Dl=0.01
+    aerodynamic_conductance!(thas; Gb_model=Val(:Su_2001), 
+        Dl, LAI=thal.LAI, zh=thal.zh, zr=thal.zr);
 end
 
 # tha48 and thal see runtests.jl
 
 
 show(thaf.wind[1:48])
+show(thaf.VPD[1:48])
+show(thaf.LE[1:48])
+show(thaf.Rn[1:48])
+show(thaf.G[1:48])
+
 
 dfGPPd = @pipe tha |> 
     transform(_, :datetime => ByRow(yearmonthday) => :ymd, copycols = false) |>
diff --git a/inst/walkthrough_todo.jmd b/inst/walkthrough_todo.jmd
index 6526dae..b599d98 100644
--- a/inst/walkthrough_todo.jmd
+++ b/inst/walkthrough_todo.jmd
@@ -236,26 +236,27 @@ Both surface temperature and VPD are in most cases higher than the ones measured
 
 ## Surface conductance
 
-Knowledge on $G_a$ allows us to calculate the bulk surface conductance ($G_s$) of the site (In this case by inverting the Penman-Monteith equation). Gs represents the combined conductance of the vegetation and the soil to water vapor transfer (and as such it is not a purely physiological quantity). Calculating $G_s$ in `Bigleaf.jl` is simple:
+Knowledge on $G_a$ allows us to calculate the bulk surface conductance ($G_s$) of the site 
+(In this case by inverting the Penman-Monteith equation). Gs represents the combined 
+conductance of the vegetation and the soil to water vapor transfer (and as such it is not 
+a purely physiological quantity). Calculating $G_s$ in `Bigleaf.jl` is simple:
 
-```julia
-summary(surface_conductance(thaf))
-```
-The function output is another DataFrame with two columns which only differ in the unit of $G_s$ (i.e. a hopeless attempt to make both physicists and physiologists happy). One in m s$^{-1}$ and one in mol m$^{-2}$ s$^{-1}$. In this function we have ignored the ground heat flux ($G$) and the storage fluxes ($S$), and the function politely reminds us of this omission by printing the first two lines of the output (it also tells us what it does, it assumes they are 0 in each time step). In this case we do not have information on the storage fluxes, but we have measurements on the ground heat flux, which we should add to the function call:
-
-```julia
-Gs = surface_conductance(thaf,G="G")
-summary(Gs)
-thaf = cbind(thaf,Gs)
+```@example doc
+surface_conductance!(thas, Val(:PenmanMonteith));
+thas[1:3,Cols(:datetime, r"Gs")]
 ```
 
-Again, we have added the two output columns to our DataFrame `thaf`.
-
-
-
-
-
+The two columns only differ in the unit of $G_s$. 
+One in m s$^{-1}$ and one in mol m$^{-2}$ s$^{-1}$. 
+In this function we have ignored the ground heat flux ($G$) and the storage fluxes ($S$).
+By default they are assumed zero.
+In our example we do not have information on the storage fluxes, but we have measurements 
+on the ground heat flux, which we should add to the function call:
 
+```@example doc
+surface_conductance!(thas, Val(:PenmanMonteith); G=thas.G);
+thas[1:3,Cols(:datetime, r"Gs")]
+```
 
 
 
diff --git a/src/Bigleaf.jl b/src/Bigleaf.jl
index b00c3c9..691dce8 100644
--- a/src/Bigleaf.jl
+++ b/src/Bigleaf.jl
@@ -38,6 +38,7 @@ export Monin_Obukhov_length, Monin_Obukhov_length!, stability_parameter,
     stability_parameter!, stability_correction, stability_correction!,
     roughness_parameters, Reynolds_Number
 export aerodynamic_conductance!, add_Ga, add_Ga!, compute_Ram, compute_Ram!, roughness_z0h
+export surface_conductance, surface_conductance!
 
 # shorthand type: either DataFrame or column table (Tables.columntable)
 DFTable = Union{NamedTuple,AbstractDataFrame}
@@ -54,5 +55,6 @@ include("stability_correction.jl")
 include("surface_roughness.jl")
 include("boundary_layer_conductance.jl")
 include("aerodynamic_conductance.jl")
+include("surface_conductance.jl")
 
 end
diff --git a/src/evapotranspiration.jl b/src/evapotranspiration.jl
index bd26dc6..9ac2382 100755
--- a/src/evapotranspiration.jl
+++ b/src/evapotranspiration.jl
@@ -139,7 +139,7 @@ function potential_ET(Tair, pressure, Rn, G, S, ::Val{:PriestleyTaylor};
   constants=bigleaf_constants())
   #
   gamma  = psychrometric_constant(Tair,pressure;constants)
-  Delta  = Esat_from_Tair_deriv(Tair; formula = Esat_formula,constants)
+  Delta  = Esat_from_Tair_deriv(Tair; Esat_formula = Esat_formula,constants)
   LE_pot = (alpha * Delta * (Rn - G - S)) / (Delta + gamma)
   ET_pot = LE_to_ET(LE_pot,Tair)
   (ET_pot = ET_pot, LE_pot = LE_pot)
@@ -155,7 +155,7 @@ function potential_ET(Tair, pressure, Rn, VPD, Ga, G, S, ::Val{:PenmanMonteith};
   constants=bigleaf_constants())
   #
   gamma  = psychrometric_constant(Tair,pressure;constants)
-  Delta  = Esat_from_Tair_deriv(Tair; formula = Esat_formula,constants)
+  Delta  = Esat_from_Tair_deriv(Tair; Esat_formula = Esat_formula,constants)
   Gs_pot = mol_to_ms(Gs_pot,Tair,pressure;constants)
   rho    = air_density(Tair,pressure;constants)
   LE_pot = (Delta * (Rn - G - S) + rho * constants[:cp] * VPD * Ga) / 
@@ -321,7 +321,7 @@ function equilibrium_imposed_ET(Tair,pressure,VPD,Gs, Rn, G, S;
   # 
   rho    = air_density(Tair, pressure; constants)
   gamma  = psychrometric_constant(Tair, pressure; constants)
-  Delta  = Esat_from_Tair_deriv(Tair; formula = Esat_formula, constants)
+  Delta  = Esat_from_Tair_deriv(Tair; Esat_formula = Esat_formula, constants)
   LE_eq  = (Delta * (Rn - G - S)) / (gamma + Delta)
   LE_imp = (rho * constants[:cp] * Gs * VPD) / gamma
   #
@@ -363,12 +363,3 @@ function decoupling()
 end
 
 
-"""
-    TODO; implement surface_conductance.
-
-This stub is there to satisfy links im Help-pages.
-"""
-function surface_conductance()
-  error("not yet implemented.")
-end
-
diff --git a/src/meteorological_variables.jl b/src/meteorological_variables.jl
index c1883fa..0e1a7a5 100644
--- a/src/meteorological_variables.jl
+++ b/src/meteorological_variables.jl
@@ -203,7 +203,7 @@ function wetbulb_temp_from_e_Tair_gamma(ea, Tair, gamma; accuracy=1e-03,
     @warn ("'accuracy' is set to 1 degC")
     accuracy = one(accuracy)
   end
-  fopt(Tw) = abs(ea - (Esat_from_Tair(Tw; formula = Esat_formula,constants) - 
+  fopt(Tw) = abs(ea - (Esat_from_Tair(Tw; Esat_formula = Esat_formula,constants) - 
     0.93*gamma*(Tair - Tw)))
   resopt = optimize(fopt, -100, 100, abs_tol = accuracy)
   roundmult(resopt.minimizer, accuracy)
@@ -261,7 +261,7 @@ function dew_point_from_e(ea;accuracy=1e-03,Esat_formula=Val(:Sonntag_1990), con
     @warn ("'accuracy' is set to 1 degC")
     accuracy = one(accuracy)
   end
-  fopt(Td) = abs(ea - Esat_from_Tair(Td;formula = Esat_formula,constants))
+  fopt(Td) = abs(ea - Esat_from_Tair(Td;Esat_formula = Esat_formula,constants))
   resopt = optimize(fopt, -100, 100, abs_tol = accuracy)
   roundmult(resopt.minimizer, accuracy)
 end
diff --git a/src/surface_conductance.jl b/src/surface_conductance.jl
new file mode 100755
index 0000000..b5ebd95
--- /dev/null
+++ b/src/surface_conductance.jl
@@ -0,0 +1,111 @@
+"""
+      surface_conductance(::Val{:FluxGradient},   Tair,pressure,VPD,LE; constants)
+      surface_conductance(::Val{:PenmanMonteith}, Tair,pressure,VPD,LE,Rn,Ga_h; 
+            G = 0.0, S = 0.0, Esat_formula=Val(:Sonntag_1990), constants)
+
+      surface_conductance!(df, method::Val{:FluxGradient}; kwargs...)
+      surface_conductance!(df, method::Val{:PenmanMonteith}; G = 0.0, S = 0.0, kwargs...)
+
+Calculate surface conductance to water vapor from the inverted Penman-Monteith equation 
+or from a simple flux-gradient approach.
+
+# Arguments
+- `Tair`      : Air temperature (deg C)
+- `pressure`  : Atmospheric pressure (kPa)
+- `Rn`        : Net radiation (W m-2)
+- `df`        : DataFrame with above variables
+- `constants=`[`bigleaf_constants`](@ref)`()`: Dictionary with physical constants
+
+additional for PenmanMonteith
+- `LE`        : Latent heat flux (W m-2)
+- `VPD`       : Vapor pressure deficit (kPa)
+- `Ga_h`      : Aerodynamic conductance towater vapor (m s-1), assumed equal to that of heat
+- `G=0.0`     : Ground heat flux (W m-2), defaults to zero
+- `S=0.0`     : Sum of all storage fluxes (W m-2), defaults to zero
+- `Esat_formula=Val(:Sonntag_1990)`: formula used in [`Esat_from_Tair`](@ref)
+
+# Details
+For `Val(:PenmanMonteith)`, surface conductance (Gs) in m s-1 
+is calculated from the inverted Penman-Monteith equation:
+
+``G_s = ( LE \\. G_a \\, \\gamma ) / ( \\Delta \\, A + \\rho \\, c_p \\, G_a \\, VPD - 
+LE \\, (\\Delta + \\gamma ))``
+
+Where ``\\gamma`` is the psychrometric constant (kPa K-1), ``\\Delta`` is the slope of the 
+saturation vapor pressure curve (kPa K^-1), and ``\\rho`` is air density (kg m-3).
+Available energy (A) is defined as ``A = R_n - G - S``. 
+
+Ground heat flux and total storage flux can be provided as scalars or vectors of
+the lenght of the DataFrame in the DataFrame variant.
+While the bigleaf R package by default converts any missings in `G` and `S` to 0,
+in `Bigleaf.jl` the caller must take care, e.g. by using `G = coalesce(myGvector, 0.0)`.
+
+For `Val(:FluxGradient)`, Gs (in mol m-2 s-1) is calculated from VPD and ET only:
+
+   ``Gs = ET/p \\, VPD``
+
+where ET is in mol m-2 s-1 and p is pressure. 
+Note that this formulation assumes fully coupled conditions (i.e. Ga = inf).
+This formulation is equivalent to the inverted form of Eq.6 in McNaughton & Black 1973:
+
+   ``Gs = LE \\, \\gamma / (\\rho \\, c_p \\, VPD)``
+
+which gives Gs in m s-1. Note that Gs > Gc (canopy conductance) under conditions 
+when a significant fraction of ET comes from interception or soil evaporation. 
+
+If `pressure` is not available, it can be approximated by elevation using the 
+function [`pressure_from_elevation`](@ref)
+ 
+# Value
+NamedTuple with entries: 
+ - Gs_ms: Surface conductance in m s-1
+ - Gs_mol: Surface conductance in mol m-2 s-1
+
+# References
+- Monteith, J., 1965: Evaporation and environment. In Fogg, G. E. (Ed.),
+  The state and movement of water in living organisms (pp.205-234).
+  19th Symp. Soc. Exp. Biol., Cambridge University Press, Cambridge
+- McNaughton, K_G., Black, T_A., 1973: A study of evapotranspiration
+  from a Douglas Fir forest using the energy balance approach. Water
+  Resources Research 9, 1579-1590.
+
+# Examples
+```jldoctest; output = false
+Tair,pressure,VPD,LE,Rn,Ga_h,G = (14.8, 97.7, 1.08, 183.0, 778.0, 0.116, 15.6)
+Gs = surface_conductance(Val(:PenmanMonteith), Tair,pressure,VPD,LE,Rn,Ga_h;G)
+isapprox(Gs.Gs_mol, 0.28, atol=0.1)
+# output
+true
+``` 
+"""
+function surface_conductance(::Val{:FluxGradient}, Tair,pressure,VPD,LE; constants=bigleaf_constants())
+    Gs_mol = (LE_to_ET(LE,Tair)/constants[:Mw]) * pressure / VPD
+    Gs_ms  = mol_to_ms(Gs_mol,Tair,pressure)
+    (;Gs_ms,Gs_mol)   
+end     
+function surface_conductance!(df::AbstractDataFrame, method::Val{:FluxGradient}; kwargs...)
+      fr = (args...) -> surface_conductance(method, args...; kwargs...) 
+      transform!(df, SA[:Tair,:pressure,:VPD,:LE] => ByRow(fr) => AsTable)
+end
+
+function surface_conductance(::Val{:PenmanMonteith}, Tair,pressure,VPD,LE,Rn,Ga_h; 
+      G = 0.0, S = 0.0,
+      #Ga="Ga_h",missing_G_as_NA=false,missing_S_as_NA=false,
+      Esat_formula=Val(:Sonntag_1990), constants=bigleaf_constants()
+      )
+      Delta = Esat_from_Tair_deriv(Tair; Esat_formula, constants)
+      gamma = psychrometric_constant(Tair, pressure; constants)
+      rho   = air_density(Tair, pressure; constants)
+      Gs_ms  = ( LE * Ga_h * gamma ) / 
+      (Delta * (Rn-G-S) + rho * constants[:cp] * Ga_h * VPD - LE * (Delta + gamma))
+      Gs_mol = ms_to_mol(Gs_ms,Tair,pressure)
+      (;Gs_ms,Gs_mol)    
+end
+function surface_conductance!(df::AbstractDataFrame, method::Val{:PenmanMonteith}; 
+      G = 0.0, S = 0.0, kwargs...)
+      # need to convert G and S to positional arguments for proper broadcast
+      fpos = (G, S, args...) -> surface_conductance(method, args...; G,S, kwargs...) 
+      ft = (args...) -> fpos.(G, S, args...)
+      transform!(df, SA[:Tair,:pressure,:VPD,:LE,:Rn,:Ga_h] => ft => AsTable)
+end
+    
diff --git a/src/unit_conversions.jl b/src/unit_conversions.jl
index a0cc66f..b3fc29e 100644
--- a/src/unit_conversions.jl
+++ b/src/unit_conversions.jl
@@ -1,13 +1,13 @@
 """
-    Esat_slope(Tair; formula, constants) 
-    Esat_from_Tair(Tair; formula, constants) 
-    Esat_from_Tair_deriv(Tair; formula, constants) 
+    Esat_slope(Tair; Esat_formula, constants) 
+    Esat_from_Tair(Tair; Esat_formula, constants) 
+    Esat_from_Tair_deriv(Tair; Esat_formula, constants) 
 
 Saturation Vapor Pressure (Esat) and Slope of the Esat Curve
 
 # Arguemtns
 - `Tair`:      Air temperature (deg C)
-- `formula=Val(:Sonntag_1990)`:   Formula to be used. Either 
+- `Esat_formula=Val(:Sonntag_1990)`:   Esat_formula to be used. Either 
    `Val(:Sonntag_1990)` (Default), `Val(:Alduchov_1996`, or `Val(:Allen_1998)`.
 - `constants=bigleaf_constants()`: Dictionary with entry  :Pa2kPa 
 
@@ -15,7 +15,7 @@ Saturation Vapor Pressure (Esat) and Slope of the Esat Curve
 Esat (kPa) is calculated using the Magnus equation:
 
 ```Esat = a * exp((b * Tair) / (c + Tair)) / 1000```
-where the coefficients a, b, c take different values depending on the formula use
+where the coefficients a, b, c take different values depending on the Esat_formula use
 The default values are from Sonntag 1990 (a=611.2, b=17.62, c=243.12). This versi
 of the Magnus equation is recommended by the WMO (WMO 2008; p1.4-29). Alternativel
 parameter values determined by Alduchov & Eskridge 1996 or Allen et al. 1998 can b
@@ -50,19 +50,19 @@ Esat_from_Tair_deriv(20.0)    # its derivative to temperature in kPa K-1
 Esat_slope(20.0)              # both as a tuple
 ```
 """
-function Esat_slope(Tair::Number; formula=Val(:Sonntag_1990), constants=bigleaf_constants()) 
-  Esat = Esat_from_Tair(Tair; formula, constants)
-  Delta = Esat_from_Tair_deriv(Tair; formula, constants)
+function Esat_slope(Tair::Number; Esat_formula=Val(:Sonntag_1990), constants=bigleaf_constants()) 
+  Esat = Esat_from_Tair(Tair; Esat_formula, constants)
+  Delta = Esat_from_Tair_deriv(Tair; Esat_formula, constants)
   Esat, Delta
 end,
-function Esat_from_Tair(Tair; formula=Val(:Sonntag_1990), constants=bigleaf_constants()) 
-  a,b,c = get_EsatCoef(formula)
+function Esat_from_Tair(Tair; Esat_formula=Val(:Sonntag_1990), constants=bigleaf_constants()) 
+  a,b,c = get_EsatCoef(Esat_formula)
   Esat = a * exp((b * Tair) / (c + Tair)) * constants[:Pa2kPa]
 end,
-function Esat_from_Tair_deriv(Tair; formula=Val(:Sonntag_1990), constants=bigleaf_constants()) 
+function Esat_from_Tair_deriv(Tair; Esat_formula=Val(:Sonntag_1990), constants=bigleaf_constants()) 
   # slope of the saturation vapor pressure curve
   #Delta = eval(D(expression(a * exp((b * Tair) / (c + Tair))),name="Tair"))
-  a,b,c = get_EsatCoef(formula)
+  a,b,c = get_EsatCoef(Esat_formula)
   #Delta_Pa = @. a*(b / (Tair + c) + (-Tair*b) / ((Tair + c)^2))*exp((Tair*b) / (Tair + c))
   Delta_Pa = @. a * (exp((b * Tair)/(c + Tair)) * (b/(c + Tair) - (b * Tair)/(c + Tair)^2))
   Delta = Delta_Pa .* constants[:Pa2kPa]
@@ -182,7 +182,7 @@ Foken, T, 2008: Micrometeorology_ Springer, Berlin, Germany.
 """
 function VPD_to_rH(VPD,Tair; Esat_formula=Val(:Sonntag_1990),
                       constants=bigleaf_constants())
-  esat = Esat_from_Tair(Tair; formula = Esat_formula,constants)
+  esat = Esat_from_Tair(Tair; Esat_formula = Esat_formula,constants)
   rH   = 1 - VPD/esat
 end,
 function rH_to_VPD(rH,Tair; Esat_formula=Val(:Sonntag_1990),
@@ -190,12 +190,12 @@ function rH_to_VPD(rH,Tair; Esat_formula=Val(:Sonntag_1990),
   if !ismissing(rH) && rH > 1 
     @warn("Expected relative humidity (rH) between 0 and 1, but was ", rH)
   end
-  esat = Esat_from_Tair(Tair; formula = Esat_formula,constants)
+  esat = Esat_from_Tair(Tair; Esat_formula = Esat_formula,constants)
   VPD  = esat - rH*esat
 end,
 function e_to_rH(e,Tair; Esat_formula=Val(:Sonntag_1990),
                     constants=bigleaf_constants())
-  esat = Esat_from_Tair(Tair; formula = Esat_formula,constants)
+  esat = Esat_from_Tair(Tair; Esat_formula = Esat_formula,constants)
   if !ismissing(e) && (e > esat + sqrt(eps()))
     @warn("Provided vapour pressure that was higher than saturation_
              Returning rH=1 for those cases.")
@@ -204,12 +204,12 @@ function e_to_rH(e,Tair; Esat_formula=Val(:Sonntag_1990),
 end,
 function VPD_to_e(VPD,Tair; Esat_formula=Val(:Sonntag_1990),
                      constants=bigleaf_constants())
-  esat = Esat_from_Tair(Tair; formula = Esat_formula,constants)
+  esat = Esat_from_Tair(Tair; Esat_formula = Esat_formula,constants)
   e    = esat - VPD
 end,
 function e_to_VPD(e,Tair; Esat_formula=Val(:Sonntag_1990),
                      constants=bigleaf_constants())
-  esat = Esat_from_Tair(Tair; formula = Esat_formula,constants)
+  esat = Esat_from_Tair(Tair; Esat_formula = Esat_formula,constants)
   VPD  = esat - e
 end,
 function e_to_q(e,pressure; constants=bigleaf_constants())
@@ -220,13 +220,13 @@ function q_to_e(q,pressure; constants=bigleaf_constants())
 end,
 function q_to_VPD(q,Tair,pressure; Esat_formula=Val(:Sonntag_1990),
                      constants=bigleaf_constants())
-  esat = Esat_from_Tair(Tair; formula = Esat_formula,constants)
+  esat = Esat_from_Tair(Tair; Esat_formula = Esat_formula,constants)
   e    = q_to_e(q,pressure; constants)
   VPD  = esat - e
 end,
 function VPD_to_q(VPD,Tair,pressure; Esat_formula=Val(:Sonntag_1990),
                      constants=bigleaf_constants())
-  esat = Esat_from_Tair(Tair; formula = Esat_formula,constants)
+  esat = Esat_from_Tair(Tair; Esat_formula = Esat_formula,constants)
   e    = esat - VPD
   q    = e_to_q(e,pressure;constants)
 end
diff --git a/test/aerodynamic_conductance.jl b/test/aerodynamic_conductance.jl
index 98f7b95..d7f6beb 100644
--- a/test/aerodynamic_conductance.jl
+++ b/test/aerodynamic_conductance.jl
@@ -76,7 +76,7 @@ end
 end
 
 @testset "aerodynamic_conductance! Gb_Thom" begin
-    df = copy(tha48)
+    df = tha48[:,Not(Cols(:Gb_h, :Ga_h))]
     # test Ga_m by wind_profile
     @inferred aerodynamic_conductance!(
         df; Gb_model=Val(:Thom_1972), Ram_model=Val(:wind_profile), 
@@ -87,7 +87,7 @@ end
 
 @testset "aerodynamic_conductance! constant_kB1" begin
     # specify kB1 as a vecttor
-    df = copy(tha48)
+    df = tha48[:,Not(Cols(:Gb_h, :Ga_h))]
     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, 
@@ -99,7 +99,7 @@ end
 
 @testset "aerodynamic_conductance! Gb_Choudhury" begin
     leafwidth=0.1
-    df = copy(tha48)
+    df = tha48[:,Not(Cols(:Gb_h, :Ga_h))]
     # test Ga_m by wind_profile
     @inferred aerodynamic_conductance!(
         df; Gb_model=Val(:Choudhury_1988), Ram_model=Val(:wind_profile), 
@@ -113,7 +113,7 @@ end
 
 @testset "aerodynamic_conductance! Gb_Su" begin
     Dl=0.01
-    df = copy(tha48)
+    df = tha48[:,Not(Cols(:Gb_h, :Ga_h))]
     df[!,:LAI] .= thal.LAI
     df.LAI[1] = thal.LAI/2
     # test providing changing LAI
diff --git a/test/boundary_layer_conductance.jl b/test/boundary_layer_conductance.jl
index d7ae722..74d2db7 100644
--- a/test/boundary_layer_conductance.jl
+++ b/test/boundary_layer_conductance.jl
@@ -80,7 +80,7 @@ end
     @test ismissing(Gbm)
     #
     # DataFrame variant
-    df = copy(tha48)
+    df = tha48[:,Not(:Gb_h)]
     # sum(skipmissing(...)) not inferrable in Julia 1.6
     # @descend_code_warntype roughness_parameters(
     #     Val(:wind_profile), df.ustar, df.wind, df.Tair, df.pressure, df.H; 
@@ -104,7 +104,7 @@ end
     @test last(propertynames(df)) == :Gb_h
     @test df.Gb_h[1] ≈ Gb_Choud rtol=1e-6
     #
-    df = copy(tha48)
+    df = tha48[:,Not(:Gb_h)]
     @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
diff --git a/test/runtests.jl b/test/runtests.jl
index ce64676..d522fb6 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -16,9 +16,27 @@ tha48 = DataFrame(
     Tair = [11.88000011444092, 11.67000007629395, 11.1899995803833, 10.80000019073486, 10.67000007629395, 10.13000011444092, 9.779999732971191, 9.5, 9.09000015258789, 8.800000190734863, 8.6899995803833, 8.899999618530273, 9.430000305175781, 9.710000038146973, 10.55000019073486, 11.19999980926514, 11.72999954223633, 12.68000030517578, 13.1899995803833, 13.46000003814697, 14.1899995803833, 14.73999977111816, 14.65999984741211, 14.8100004196167, 15.02999973297119, 14.98999977111816, 14.77999973297119, 15.35000038146973, 15.5, 15.61999988555908, 15.47000026702881, 16.20000076293945, 16.13999938964844, 15.56999969482422, 15.14999961853027, 15.18000030517578, 14.81999969482422, 14.60000038146973, 14.10000038146973, 13.23999977111816, 12.65999984741211, 12.46000003814697, 12.09000015258789, 11.85999965667725, 11.72000026702881, 11.64000034332275, 11.52999973297119, 11.43000030517578],
     pressure = [97.63999938964844, 97.62999725341797, 97.61000061035156, 97.61000061035156, 97.61000061035156, 97.61000061035156, 97.61000061035156, 97.61000061035156, 97.62999725341797, 97.66000366210938, 97.66999816894531, 97.69000244140625, 97.69000244140625, 97.69999694824219, 97.70999908447266, 97.69999694824219, 97.7300033569336, 97.72000122070312, 97.70999908447266, 97.69999694824219, 97.69999694824219, 97.69999694824219, 97.69999694824219, 97.70999908447266, 97.70999908447266, 97.72000122070312, 97.70999908447266, 97.70999908447266, 97.69999694824219, 97.68000030517578, 97.68000030517578, 97.66999816894531, 97.68000030517578, 97.66999816894531, 97.66999816894531, 97.66000366210938, 97.6500015258789, 97.63999938964844, 97.63999938964844, 97.6500015258789, 97.66000366210938, 97.66999816894531, 97.69000244140625, 97.68000030517578, 97.68000030517578, 97.69000244140625, 97.69000244140625, 97.69000244140625],
     H = Union{Missing, Float64}[-68.18000030517578, -48.54000091552734, -59.09999847412109, -60.11000061035156, -59.40000152587891, -48.91999816894531, -53.52999877929688, -40.75, -29.28000068664551, -33.36999893188477, -9.050000190734863, 2.400000095367432, 38.29999923706055, 38.45000076293945, 85.08000183105469, 108.7799987792969, 176.6799926757812, 210.6900024414062, 230.0800018310547, 227.4199981689453, 321.6499938964844, 358.8999938964844, 329.9299926757812, 317.3099975585938, 375.1900024414062, 336.489990234375, 260.7099914550781, 360.6400146484375, 264.2200012207031, 185.3200073242188, 123.7699966430664, 236.0200042724609, 207.5200042724609, 128.3399963378906, 66.05999755859375, 79.6500015258789, 44.93999862670898, 16.22999954223633, -22.35000038146973, -29.10000038146973, -47.79999923706055, -63.61999893188477, -61.38000106811523, -65.11000061035156, -59.77000045776367, -60.52000045776367, -49.02999877929688, -53.45000076293945],
+    LE = Union{Missing, Float64}[9.9399995803833, 5.269999980926514, 3.980000019073486, missing, missing, missing, missing, missing, 5.190000057220459, 10.11999988555908, 26.27000045776367, 27.3700008392334, 23.52000045776367, 6.639999866485596, 52.2599983215332, 54.02000045776367, 92.27999877929688, 100.9199981689453, 111.7200012207031, 223.1000061035156, 200.7400054931641, 185.0500030517578, 170.9799957275391, 183.4900054931641, 187.6900024414062, 223.0299987792969, 161.6600036621094, 160.9700012207031, 153.7599945068359, 77.20999908447266, 118.4300003051758, 31.86000061035156, 116.4100036621094, 84.98999786376953, 62.2400016784668, 68.12999725341797, 29.88999938964844, 41.65000152587891, 28.77000045776367, 16.53000068664551, 16.04000091552734, 10.67000007629395, 2.880000114440918, 5.369999885559082, 1.820000052452087, 1.580000042915344, 4.139999866485596, 4.880000114440918],
+    Rn = [-86.48999786376953, -84.19999694824219, -81.80999755859375, -77.9000015258789, -79.08999633789062, -79.08999633789062, -80.13999938964844, -73.36000061035156, -66.63999938964844, -37.93999862670898, -25.25, 78.0, 113.2399978637695, 121.9499969482422, 230.6199951171875, 302.1700134277344, 394.5400085449219, 453.8800048828125, 518.2100219726562, 609.9000244140625, 693.4099731445312, 729.1400146484375, 767.3599853515625, 778.1699829101562, 778.5599975585938, 778.239990234375, 606.7899780273438, 724.239990234375, 524.8499755859375, 425.1499938964844, 399.75, 515.6400146484375, 371.9100036621094, 304.1300048828125, 161.6799926757812, 176.8000030517578, 108.5299987792969, 35.54999923706055, -25.29999923706055, -72.19999694824219, -91.4800033569336, -92.30000305175781, -92.12000274658203, -90.70999908447266, -90.12000274658203, -89.77999877929688, -87.5, -86.76000213623047],
+    G = [-4.934999942779541, -5.085000038146973, -5.135000228881836, -5.210000038146973, -5.28000020980835, -5.34499979019165, -5.5, -5.804999828338623, -6.239999771118164, -6.550000190734863, -6.385000228881836, -6.039999961853027, -5.034999847412109, -4.070000171661377, -3.019999980926514, -1.475000023841858, 0.5550000071525574, 1.534999966621399, 2.630000114440918, 7.550000190734863, 22.06500053405762, 17.09499931335449, 18.93000030517578, 15.5649995803833, 16.90500068664551, 12.55000019073486, 30.14999961853027, 23.71500015258789, 14.81999969482422, 15.06999969482422, 9.925000190734863, 6.525000095367432, 4.864999771118164, 5.394999980926514, 4.920000076293945, 3.289999961853027, 2.345000028610229, 1.559999942779541, 0.6499999761581421, -0.4550000131130219, -1.669999957084656, -2.795000076293945, -3.674999952316284, -4.295000076293945, -4.775000095367432, -5.190000057220459, -5.364999771118164, -5.440000057220459],
+    VPD = [0.5745999813079834, 0.5633999824523925, 0.513700008392334, 0.4560999870300293, 0.4184000015258789, 0.3608999967575073, 0.3033999919891358, 0.2822999954223633, 0.2641999959945679, 0.2621000051498413, 0.2773000001907349, 0.3075999975204468, 0.3339999914169312, 0.3270999908447266, 0.373799991607666, 0.4267000198364258, 0.4580999851226807, 0.5138000011444092, 0.5421999931335449, 0.6954999923706054, 0.8619000434875488, 1.010499954223633, 1.014999961853027, 1.075800037384033, 1.090100002288818, 1.094499969482422, 0.9982000350952148, 1.085700035095215, 1.094299983978271, 1.09379997253418, 1.096599960327148, 1.164799976348877, 1.048600006103516, 0.9133999824523926, 0.8229999542236328, 0.8409000396728515, 0.7889999866485595, 0.7486999988555908, 0.6928999900817872, 0.601800012588501, 0.5798999786376953, 0.5870999813079834, 0.5499000072479248, 0.5432000160217285, 0.5427000045776367, 0.5219999790191651, 0.5199999809265137, 0.511299991607666],
     wind = Union{Missing, Float64}[4.210000038146973, 4.460000038146973, 4.539999961853027, 4.079999923706055, 3.950000047683716, 4.019999980926514, 3.269999980926514, 3.480000019073486, 3.240000009536743, 3.420000076293945, 3.440000057220459, 2.970000028610229, 3.349999904632568, 2.660000085830688, 2.779999971389771, 2.329999923706055, 2.220000028610229, 2.160000085830688, 2.519999980926514, 2.980000019073486, 2.359999895095825, 2.420000076293945, 2.920000076293945, 3.359999895095825, 2.759999990463257, 3.279999971389771, 3.410000085830688, 3.480000019073486, 3.039999961853027, 3.950000047683716, 2.700000047683716, 2.829999923706055, 2.730000019073486, 3.640000104904175, 2.75, 3.299999952316284, 2.630000114440918, 2.049999952316284, 1.879999995231628, 2.130000114440918, 2.400000095367432, 2.619999885559082, 2.640000104904175, 2.589999914169312, 2.589999914169312, 2.819999933242798, 2.779999971389771, 2.660000085830688],
+    #
+    # from Su_2001 of entire tha-dataset (see tmpf below)
+    Gb_h = Union{Missing, Float64}[0.1822661462615249, 0.1750029796115819, 0.16439640538330097, 0.15115204867322993, 0.17536657313260987, 0.1625005163806663, 0.18840537219804904, 0.18879242029415302, 0.17645201287966256, 0.1538975543389697, 0.21088990351067222, 0.2159981919177606, 0.26124601608024695, 0.2773512635093451, 0.31590536172758477, 0.40377826077805495, 0.4767947347202397, 0.56854692623082, 0.42750700682547943, 0.4686532838500691, 0.4713224944023968, 0.5367616432054471, 0.4390590733192645, 0.4365982669032558, 0.45028995760447604, 0.44342264360628164, 0.3870369356275578, 0.4091200789143718, 0.38965757992457295, 0.3741648297018856, 0.3598397227673043, 0.41828517646506513, 0.37765531002288344, 0.346427027859903, 0.2926818717409249, 0.3033011874663078, 0.2735927127688745, 0.23189051804555683, 0.13765984163452955, 0.1230228546581407, 0.10194295317288511, 0.11062842016540558, 0.11745785563008848, 0.11461165845075606, 0.107937932508708, 0.11269011080949713, 0.11690251453131727, 0.1020485594839973],
+    Ga_h = Union{Missing, Float64}[0.050190555728311555, 0.04116957600962841, 0.03877814728450016, 0.03736361775914192, 0.047872524670215494, 0.03975837173989957, 0.06052647712603172, 0.05504593321238934, 0.04766632363024636, 0.03725385088675965, 0.05886451561135729, 0.06231168613359922, 0.06166420574946217, 0.05146825750242876, 0.07874566459614012, 0.07414026084897143, 0.09405110045117922, 0.09616429374872118, 0.1150966610946552, 0.09097297096976756, 0.13839881444493318, 0.1317378026984531, 0.13140678312765053, 0.11633916710154416, 0.1454360006411824, 0.12128624346147493, 0.12212101975187602, 0.14200569128979465, 0.1322214598437316, 0.0871694750533624, 0.08780626502144728, 0.11251759597304214, 0.12401547627348973, 0.07893328170976276, 0.08836491703204545, 0.08219874492644172, 0.07681770620097074, 0.06927204867424476, 0.042503890501350394, 0.036116991360095725, 0.03139917116540786, 0.03548953025999413, 0.03731866634395528, 0.03750740234735876, 0.034189028583981027, 0.03392950934600214, 0.03332805495857941, 0.030478849592838437],
 );
 
+function tmpf() 
+    Dl=0.01
+    df = aerodynamic_conductance!(copy(df, copycols=false); Gb_model=Val(:Su_2001), 
+        Dl, LAI=thal.LAI, zh=thal.zh, zr=thal.zr);
+    @show df.Gb_h[1:48];
+    @show df.Ga_h[1:48];
+end
+
+
+
 thal = (
     LAI = 7.6,   # leaf area index
     zh  = 26.5,  # average vegetation height (m)
@@ -60,6 +78,9 @@ thal = (
     @testset "aerodynamic_conductance" begin
         include("aerodynamic_conductance.jl")
     end
+    @testset "surface_conductance" begin
+        include("surface_conductance.jl")
+    end
     @testset "evapotranspiration" begin
         include("evapotranspiration.jl")
     end
diff --git a/test/stability_correction.jl b/test/stability_correction.jl
index 3bff533..1830a90 100644
--- a/test/stability_correction.jl
+++ b/test/stability_correction.jl
@@ -1,12 +1,10 @@
 @testset "Monin_Obukhov_length" begin
-    datetime, ustar, Tair, pressure, H = values(tha48[24,1:5])
+    datetime, ustar, Tair, pressure, H = tha48[24,Cols(:datetime, :ustar, :Tair, :pressure, :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)
-    dfd = disallowmissing(df)
-    dfm = allowmissing(df)
     @inferred Monin_Obukhov_length!(df)
     @test df.MOL[24] == MOL24
 end
diff --git a/test/surface_conductance.jl b/test/surface_conductance.jl
new file mode 100644
index 0000000..b3ecaba
--- /dev/null
+++ b/test/surface_conductance.jl
@@ -0,0 +1,27 @@
+@testset "FluxGradient" begin
+    ir = 24
+    Tair,pressure,VPD,LE = tha48[ir, Cols(:Tair,:pressure,:VPD,:LE)]
+    Gs24 = surface_conductance(Val(:FluxGradient), Tair,pressure,VPD,LE)
+    @test keys(Gs24) == (:Gs_ms, :Gs_mol)
+    # compare with R
+    @test all(isapprox.(values(Gs24), (0.00919248, 0.3751482), rtol=1e-3))
+    #
+    df = copy(tha48)
+    surface_conductance!(df, Val(:FluxGradient))
+    @test propertynames(df)[(end-1):end] == SA[:Gs_ms, :Gs_mol]
+    @test df[ir, Cols(:Gs_ms, :Gs_mol)] == Gs24
+end
+
+@testset "PenmanMonteith" begin
+    ir = 24
+    Tair,pressure,VPD,LE,Rn,Ga_h,G = tha48[ir, Cols(:Tair,:pressure,:VPD,:LE, :Rn,:Ga_h,:G)]
+    Gs24 = surface_conductance(Val(:PenmanMonteith), Tair,pressure,VPD,LE,Rn,Ga_h;G)
+    @test keys(Gs24) == (:Gs_ms, :Gs_mol)
+    # compare with R
+    @test all(isapprox.(values(Gs24), (0.006846274, 0.2793988), rtol=1e-3))
+    #
+    df = copy(tha48)
+    surface_conductance!(df, Val(:PenmanMonteith); G=df.G)
+    @test propertynames(df)[(end-1):end] == SA[:Gs_ms, :Gs_mol]
+    @test df[ir, Cols(:Gs_ms, :Gs_mol)] == Gs24
+end
\ No newline at end of file
diff --git a/test/surface_roughness.jl b/test/surface_roughness.jl
index 01be270..d40fb3d 100644
--- a/test/surface_roughness.jl
+++ b/test/surface_roughness.jl
@@ -37,7 +37,7 @@ end
     # 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)
+    dfd = disallowmissing(df[!,Not(:LE)])
     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))
@@ -76,7 +76,7 @@ end
     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)
+    dfd = disallowmissing(df[!,Not(:LE)])
     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
diff --git a/test/unit_conversions.jl b/test/unit_conversions.jl
index aa806f9..ae41fe6 100644
--- a/test/unit_conversions.jl
+++ b/test/unit_conversions.jl
@@ -1,12 +1,12 @@
 @testset "Esat_from_Tair" begin
   Tair = 15
-  formula=Val(:Sonntag_1990)
+  Esat_formula=Val(:Sonntag_1990)
   constants=bigleaf_constants()
-  eSat = Esat_from_Tair(Tair; formula, constants)
+  eSat = Esat_from_Tair(Tair; Esat_formula, constants)
   # regression test
   @test isapprox(eSat, 1.702, atol=1e-3)
-  eSat2 = Esat_from_Tair(Tair; formula = Val(:Alduchov_1996), constants)
-  eSat3 = Esat_from_Tair(Tair; formula = Val(:Allen_1998), constants)
+  eSat2 = Esat_from_Tair(Tair; Esat_formula = Val(:Alduchov_1996), constants)
+  eSat3 = Esat_from_Tair(Tair; Esat_formula = Val(:Allen_1998), constants)
   # different results with different methods
   @test !(eSat2 ≈ eSat)
   @test !(eSat3 ≈ eSat)
@@ -19,7 +19,7 @@ end
 @testset "Esat_from_Tair_slope" begin
     step = 0.1
     Tair = range(20.0, 22.0, step = step)
-    formula=Val(:Sonntag_1990)
+    Esat_formula=Val(:Sonntag_1990)
     constants=bigleaf_constants()
     Esat = Esat_from_Tair.(Tair)
     delta = Esat_from_Tair_deriv.(Tair)