implement surface conductance

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

docs/src/walkthrough.md

@@ -295,24 +295,24 @@ 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,
 #        )
 #end;
 ##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}$ 

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``                             
 #'                                       

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) 

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 

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)  

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

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

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

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")]
+```
 
 
 

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

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
-