#' Run the LWF-Brook90 hydrologic model
#'
#' Sets up the input objects for the LWF-Brook90 hydrologic model, starts the
#' model, and returns the selected results.
#'
#' @param options_b90 Named list of model control options. Use
#' \code{\link{set_optionsLWFB90}} to generate a list with default model
#' control options.
#' @param param_b90 Named list of model input parameters. Use
#' \code{\link{set_paramLWFB90}} to generate a list with default model
#' parameters.
#' @param climate Data.frame with daily climatic data, or a function that
#' returns a suitable data.frame. See details for the required variables.
#' @param precip Data.frame with columns 'dates' and 'prec' to supply
#' precipitation data separately from climate data. Can be used to provide
#' sub-day resolution precipitation data to LWFBrook90. For each day in dates,
#' 1 (daily resolution) to 240 values of precipitation can be provided, with
#' the number of values per day defined in \code{options_b90$prec_interval}.
#' @param soil Data.frame containing the hydraulic properties of the soil
#' layers. See section 'Soil parameters'
#' @param output_fun A function or a list of functions of the form
#' \code{f(x,...)}, where \code{x} is the object regularly returned by
#' \code{run_LWFB90}. During function evaluation, \code{x} contains model
#' input and selected output objects, irrespective of \code{rtrn_input} and
#' \code{rtrn_output}. Can be used to aggregate output on-the-fly, and is
#' especially useful if the function is evaluated within a large multi-run
#' application, for which the output might overload the memory. (see
#' \code{\link{run_multi_LWFB90}} and \code{\link{run_multisite_LWFB90}}).
#' @param rtrn_input Logical: append \code{param_b90}, \code{options_b90}, and
#' daily plant properties (\code{standprop_daily}, as derived from parameters)
#' to the result?
#' @param rtrn_output Logical: return the simulation results select via
#' \code{output}?
#' @param chk_input Logical wether to check \code{param_b90},
#' \code{options_b90}, \code{climate}, \code{precip}, and \code{soil} for
#' completeness and consistency.
#' @param run Logical: run LWF-Brook90 or only return model input objects?
#' Useful to inspect the effects of options and parameters on model input.
#' Default is TRUE.
#' @param timelimit Integer to set elapsed time limits (seconds) for running LWF-Brook90.
#' @param verbose Logical: print messages to the console? Default is FALSE.
#' @param ... Additional arguments passed to \code{output_fun} and/or
#' \code{climate}, if the latter is a function.
#'
#' @return A list containing the selected model output (if \code{rtrn_output ==
#' TRUE}), the model input (if \code{rtrn_input == TRUE}, except for
#' \code{climate}), and the return values of \code{output_fun} if specified.
#'
#'@section Climate input data: The \code{climate} data.frame (or function) must
#' contain (return) the following variables in columns named 'dates' (Date),
#' 'tmax' (deg C), 'tmin' (deg C), 'tmean' (deg C), 'windspeed' (m/s), 'prec'
#' (mm) , 'vappres' (kPa), and either 'globrad' ( MJ/(m²d) ) or 'sunhours' (h).
#' When using \code{sunhours}, please set \code{options_b90$fornetrad =
#' 'sunhours'}.
#'
#'@section Soil input parameters: Each row of \code{soil} represents one layer,
#' containing the layers' boundaries and soil hydraulic parameters. The column
#' names for the upper and lower layer boundaries are 'upper' and 'lower' (m,
#' negative downwards). When using \code{options_b90$imodel = 'MvG'}, the
#' hydraulic parameters are 'ths', 'thr', 'alpha' (1/m), 'npar', 'ksat' (mm/d)
#' and 'tort'. With \code{options_b90$imodel = 'CH'}, the parameters are
#' 'thsat', 'thetaf', 'psif' (kPa), 'bexp', 'kf' (mm/d), and 'wetinf'. For both
#' parameterizations, the volume fraction of stones has to be named 'gravel'.
#' If the \code{soil} argument is not provided, list items \code{soil_nodes}
#' and \code{soil_materials} of \code{param_b90} are used for the simulation.
#' These have to be set up in advance, see \code{\link{soil_to_param}}.
#'
#' @section Outputs:
#' \tabular{llcl}{
#' \strong{Name} \tab \strong{Description} \tab \strong{Unit} \cr
#' yr \tab year \tab - \cr
#' mo \tab month \tab - \cr
#' da \tab day of month \tab - \cr
#' doy \tab day of year \tab - \cr
#' aa \tab average available energy above canopy \tab W/m2 \cr
#' adef \tab available water deficit in root zone \tab mm \cr
#' asubs \tab average available energy below canopy \tab W/m2 \cr
#' awat \tab total available soil water in layers with roots between -6.18 kPa and \code{param_b90$psicr} \tab mm \cr
#' balerr \tab error in water balance (daily value, output at the day's last precipitation interval) \tab mm \cr
#' byfl \tab total bypass flow \tab mm/d \cr
#' cc \tab cold content of snowpack (positive) \tab MJ m-2 \cr
#' dsfl \tab downslope flow \tab mm/d \cr
#' evap \tab evapotranspiration \tab mm/d \cr
#' flow \tab total streamflow \tab mm/d \cr
#' gwat \tab groundwater storage below soil layers \tab mm \cr
#' gwfl \tab groundwater flow \tab mm/d \cr
#' intr \tab intercepted rain \tab mm \cr
#' ints \tab intercepted snow \tab mm \cr
#' irvp \tab evaporation of intercepted rain \tab mm/d \cr
#' isvp \tab evaporation of intercepted snow \tab mm/d \cr
#' lngnet \tab net longwave radiation \tab W/m2 \cr
#' nits \tab total number of iterations \tab - \cr
#' pint \tab potential interception for a canopy always wet \tab mm/d \cr
#' pslvp \tab potential soil evaporation \tab mm/d \cr
#' ptran \tab potential transpiration \tab mm/d \cr
#' relawat \tab relative available soil water in layers with roots \tab - \cr
#' rfal \tab rainfall \tab mm/d \cr
#' rint \tab rain interception catch rate \tab mm/d \cr
#' rnet \tab rainfall to soil surface \tab mm/d \cr
#' rsno \tab rain on snow \tab mm/d \cr
#' rthr \tab rain throughfall rate \tab mm/d \cr
#' sthr \tab snow throughfall rate \tab mm/d \cr
#' safrac \tab source area fraction \tab - \cr
#' seep \tab seepage loss \tab mm/d \cr
#' sfal \tab snowfall \tab mm/d \cr
#' sint \tab snow interception catch rate \tab mm/d \cr
#' slfl \tab input to soil surface \tab mm/d \cr
#' slvp \tab evaporation rate from soil \tab mm/d \cr
#' slrad \tab average solar radiation on slope over daytime \tab W/m2 \cr
#' solnet \tab net solar radiation on slope over daytime \tab W/m2 \cr
#' smlt \tab snowmelt \tab mm/d \cr
#' snow \tab snowpack water equivalent \tab mm \cr
#' snowlq \tab liquid water content of snow on the ground \tab mm \cr
#' snvp \tab evaporation from snowpack \tab mm/d \cr
#' srfl \tab source area flow \tab mm/d \cr
#' stres \tab tran / ptran (daily value, output at the day's last precipitation interval) \tab - \cr
#' swat \tab total soil water in all layers\tab mm \cr
#' tran \tab transpiration \tab mm/d \cr
#' vrfln \tab vertical matrix drainage from lowest layer \tab mm/d \cr
#'}
#'
#' @section Layer outputs:
#' \tabular{llcl}{
#' \strong{Name} \tab \strong{Description} \tab \strong{Unit} \cr
#' yr \tab year \tab - \cr
#' mo \tab month \tab - \cr
#' da \tab day of month \tab - \cr
#' doy \tab day of year \tab - \cr
#' nl \tab index of soil layer \tab \cr
#' swati \tab soil water volume in layer \tab mm \cr
#' theta \tab water content of soil layer, mm water / mm soil matrix \tab - \cr
#' wetnes \tab wetness of soil layer, fraction of saturation \tab - \cr
#' psimi \tab matric soil water potential for soil layer \tab kPa \cr
#' infl \tab infiltration to soil water in soil layer \tab mm/d \cr
#' byfl \tab bypass flow from soil layer \tab mm/d \cr
#' tran \tab transpiration from soil layer \tab mm/d \cr
#' vrfl \tab vertical matrix drainage from soil layer \tab mm/d \cr
#' dsfl \tab downslope drainage from layer \tab mm/d \cr
#' ntfl \tab net flow into soil layer \tab mm/d \cr
#'}
#'
#' @export
#' @example inst/examples/run_LWFB90-help.R
run_LWFB90 <- function(options_b90,
param_b90,
climate,
precip = NULL,
soil = NULL,
output_fun = NULL,
rtrn_input = TRUE,
rtrn_output = TRUE,
chk_input = TRUE,
run = TRUE,
timelimit = Inf,
verbose = FALSE,
...) {
xfunargs <- list(...)
if (is.function(climate)) {
if (verbose == TRUE) {
message("Applying climate input function")
}
climfunargsnms <- match.arg(methods::formalArgs(climate), names(xfunargs),several.ok = TRUE)
climate <- do.call(climate,xfunargs[climfunargsnms])
}
# input checks -----
if (chk_input) {
chk_options()
chk_param()
chk_clim()
chk_soil()
}
# Simulation period ----
climyears <- unique(as.integer(format(climate$dates,"%Y")))
simyears <- seq(from = as.integer(format(options_b90$startdate,"%Y")),
to = as.integer(format(options_b90$enddate,"%Y")),
by = 1)
## Number of Simulation days
param_b90$ndays <- as.integer(difftime(options_b90$enddate,options_b90$startdate)) + 1
## Vegetation-Period: calculate budburst and leaffall days of year ----
budburst_leaffall <- calc_vegperiod(dates = climate$dates, tavg = climate$tmean,
out_yrs = simyears,
budburst_method = options_b90$budburst_method,
leaffall_method = options_b90$leaffall_method,
budburstdoy.fixed = param_b90$budburstdoy,
leaffalldoy.fixed = param_b90$leaffalldoy,
species = param_b90$budburst_species,
est.prev = ifelse(length(climyears) <= 5,
length(climyears) - 1, 5),
check.data = FALSE)
param_b90$budburstdoy <- budburst_leaffall$start
param_b90$leaffalldoy <- budburst_leaffall$end
# Vegetation ----
# Prepare vegetation parameters
if (tolower(options_b90$standprop_input) == "table") {
if (verbose == TRUE) {message("Creating long term stand dynamics from table 'standprop_table'...")}
param_b90 <- standprop_yearly_to_param(param_b90$standprop_table,
param_b90,
out_yrs = simyears)
} else {
if (verbose == TRUE) {message("Creating stand properties from parameters...")}
# derive age from age_ini for simyears
param_b90$age <- seq(from = param_b90$age_ini + 1,
by = 1, length.out = length(simyears))
}
## Create daily standproperties from parameters ----
standprop_daily <- make_standprop(options_b90, param_b90, out_yrs = simyears)
### constrain to simulation period
standprop_daily <- standprop_daily[which(standprop_daily$dates >= options_b90$startdate
& standprop_daily$dates <= options_b90$enddate),]
if (verbose == TRUE) {
message("Standproperties created succesfully")
}
# Prepare climate ----
# constrain data to simulation period
climate <- climate[which(climate$dates >= options_b90$startdate
& climate$dates <= options_b90$enddate),]
if (!is.null(precip)){
precip <- precip[which(precip$dates >= options_b90$startdate
& precip$dates <= options_b90$enddate),]
}
## Precipitation correction (Richter) ----
if (options_b90$correct_prec == TRUE) {
if (!is.null(precip)) {
warning("Correction of precipitation not possible for sub-daily precipitation data! Doing nothing.")
}
climate$prec <- with(climate, correct_prec(mo, tmean, prec,
station.exposure = param_b90$prec_corr_statexp))
}
## Calculate global radiation from sunshine duration ----
if (options_b90$fornetrad == "sunhours") {
climate$globrad <- calc_globrad(climate$dates, climate$sunhours,
param_b90$coords_y)
}
## Make soilnodes & soil materials ----
# soil provided as argument -> create soil nodes and materials and add them to param_b90
if (!is.null(soil)) {
param_b90[c("soil_nodes","soil_materials" )] <- soil_to_param(soil, options_b90$imodel)
}
if (options_b90$imodel == "MvG") {
param_b90$soil_materials <- param_b90$soil_materials[,c("mat","ths","thr","alpha","npar","ksat","tort","gravel")]
} else {
param_b90$soil_materials <- param_b90$soil_materials[,c("mat","thsat","thetaf","psif","bexp","kf","wetinf","gravel")]
}
### add initial water potential from parameters
param_b90$soil_nodes$psiini <- param_b90$psiini
## Make Roots -----
if (options_b90$root_method != "soilvar") {
param_b90$soil_nodes$rootden <- make_rootden(soilnodes = c(max(param_b90$soil_nodes$upper),
param_b90$soil_nodes$lower),
maxrootdepth = param_b90$maxrootdepth,
method = options_b90$root_method,
beta = param_b90$betaroot,
rootdat = param_b90$rootden_table)
} else {
if (!is.null(soil)) {
param_b90$soil_nodes$rootden <- soil$rootden
} else {
stopifnot(!is.null(param_b90$soil_nodes$rootden))
}
}
# Execute LWF-Brook90 ----
if (run) {
if (verbose == TRUE) {
message("Running model..." )
}
start <- Sys.time()
simout <- r_lwfbrook90(
siteparam = data.frame(simyears[1],
as.integer(format(options_b90$startdate, "%j")),
param_b90$coords_y, param_b90$snowini, param_b90$gwatini,
options_b90$prec_interval, param_b90$snowlqini,
param_b90$snowccini, param_b90$water_table_depth),
climveg = cbind(climate[, c("yr", "mo", "da","globrad","tmax","tmin",
"vappres","windspeed","prec","mesfl")],
standprop_daily[, c("densef", "height", "lai", "sai", "age")]),
precdat = precip[,c("yr", "mo", "da","ii","prec", "mesfl")],
param = param_to_rlwfbrook90(param_b90, options_b90$imodel),
pdur = param_b90$pdur,
soil_materials = param_b90$soil_materials,
soil_nodes = param_b90$soil_nodes[,c("layer","midpoint", "thick", "mat", "psiini", "rootden")],
output_log = verbose,
chk_input = chk_input,
timelimit = timelimit
)
finishing_time <- Sys.time()
simtime <- finishing_time - start
units(simtime) <- "secs"
if (verbose == TRUE) {
message(paste("Simulation successful! Duration:", round(simtime,2), "seconds"))
}
## initialize return value ----
simres <- list(simulation_duration = simtime,
finishing_time = finishing_time)
## append model input ----
# might be needed for access from output_fun. Will be removed again later, if not required
simres$model_input <- list(options_b90 = options_b90,
param_b90 = param_b90,
standprop_daily = standprop_daily)
## append simulation results ----
simres[names(simout)[-1]] <- simout[-1] # append without error-code
## apply functions on simulation output ----
if (!is.null(output_fun)) {
if (verbose == TRUE) {
message("Applying custom functions on simulation output...")
}
if (!is.list(output_fun)){
output_fun <- list(output_fun)
}
outfunargs <- list(x = simres, ...)
outfunargsnms <- lapply(output_fun, FUN = function(x,argsnms) {
match.arg(methods::formalArgs(x),
argsnms,
several.ok = TRUE)},
argsnms = names(outfunargs))
# TODO: simres is copied for use in each output_fun.
# Better to name -object directly in the call to output_fun,
# instead of addressing the whole list x?
simres$output_fun <- tryCatch( {
Map(do.call, output_fun, lapply(outfunargsnms, function(x,args) args[x], args = outfunargs))
},
warning = function(wrn){return(wrn)},
error = function(err){return(err)})
}
## remove the basic results again if they are not required ----
if (!rtrn_output) {
simres <- simres[-which(names(simres) %in% c("output", "layer_output"))]
}
## remove the model_input if not required ----
if (!rtrn_input) {
simres <- simres[-which(names(simres) == "model_input")]
}
} else {
# 'dry' run = FALSE -> always return model input
return(list(options_b90 = options_b90,
param_b90 = param_b90,
standprop_daily = standprop_daily))
}
if (verbose == TRUE) {
message("Finished!")
}
return(simres)
}
## check-functions for input ---------------------------------------------------
chk_options <- function(){
eval.parent(quote({ # manipulate the calling environment
names(options_b90) <- tolower(names(options_b90))
stopifnot(all(names(options_b90) %in% c("startdate","enddate","fornetrad","prec_interval",
"correct_prec","budburst_method",
"leaffall_method", "standprop_input", "standprop_interp",
"use_growthperiod","lai_method","imodel", "root_method")))
options_b90$fornetrad <- match.arg(options_b90$fornetrad, choices = c("globrad","sunhours"))
options_b90$standprop_input <- match.arg(options_b90$standprop_input, choices = c("parameters", "table"))
options_b90$lai_method <- match.arg(options_b90$lai_method, choices = c("b90", "linear", "Coupmodel"))
options_b90$root_method <- match.arg(options_b90$root_method, choices = c("betamodel", "table", "linear", "constant", "soilvar"))
options_b90$imodel <- match.arg(options_b90$imodel, choices = c("MvG", "CH"))
if (!inherits(options_b90$startdate, "Date")) {
stop("Please provide 'options_b90$startdate' as Date-object")}
if (!inherits(options_b90$enddate, "Date")) {
stop("Please provide 'options_b90$enddate' as Date-object")}
if (!(options_b90$startdate < options_b90$enddate)) {
stop("Check options_b90: 'startdate > enddate ")}
if (options_b90$budburst_method %in% c("const", "const.", "constant")) options_b90$budburst_method <- "fixed"
if (options_b90$leaffall_method %in% c("const", "const.", "constant")) options_b90$leaffall_method <- "fixed"
}))
}
chk_param <- function() {
eval.parent(quote({ # manipulate the calling environment
names(param_b90) <- tolower(names(param_b90))
nms <- c("maxlai","sai","sai_ini","height","height_ini","densef",
"densef_ini","age_ini","winlaifrac","budburst_species","budburstdoy",
"leaffalldoy","shp_budburst","shp_leaffall","shp_optdoy","emergedur","leaffalldur",
"lai_doy","lai_frac","alb","albsn","ksnvp","fxylem",
"mxkpl","lwidth","psicr","nooutf","lpc","cs",
"czs","czr","hs","hr","rhotp","nn",
"maxrlen","initrlen","initrdep","rrad","rgrorate","rgroper",
"maxrootdepth","betaroot","radex","glmax","glmin",
"rm","r5","cvpd","tl","t1","t2",
"th","frintlai","frintsai","fsintlai","fsintsai","cintrl",
"cintrs","cintsl","cintss","infexp","bypar","qfpar",
"qffc","imperv","drain","gsc","gsp","ilayer",
"qlayer","z0s","rstemp","melfac","ccfac","laimlt",
"saimlt","grdmlt","maxlqf","snoden","obsheight","rssa",
"rssb","dtimax","dswmax","dpsimax",
"wndrat","fetch","z0w","zw","zminh","coords_x",
"coords_y","c1","c2","c3","pdur","eslope",
"aspect","dslope","slopelen","intrainini","intsnowini","gwatini",
"snowini","psiini")
if ( any(!nms %in% names(param_b90) )) {
stop(paste("param_b90-list is incomplete. Missing list items:",
paste(nms[which(!nms %in% names(param_b90))], collapse = ", ")))
}
}))
}
chk_clim <- function() {
eval.parent(quote({
# climate name checks
names(climate) <- tolower(names(climate))
if (!all(c("dates", "tmax", "tmin", options_b90$fornetrad, "vappres", "windspeed") %in% tolower(names(climate)))) {
if (is.null(precip)) {}
stop("Please check column names of 'climate'. Must contain: \n
'dates', 'tmax', 'tmin', 'vappres', 'windspeed', and 'globrad' (or 'sunhours') ")
}
stopifnot(inherits(climate$dates, "Date"))
if (min(climate$dates) > options_b90$startdate | max(climate$dates) < options_b90$enddate){
stop("climate not covering requested simulation period completely.")
}
# Climate
# if (options_b90$fornetrad == "globrad" & !any( names(climate) == "globrad" )) {
# if (any( names(climate) == "sunhour" )) {
# options_b90$fornetrad <- "sunhour"
# warning("Global radiation missing, will be calculated from sunshine duration!")
# } else {
# stop("Please either provide globrad or sunhour with your climate data!")
# }
# } else {
# if (any( names(climate) == "globrad" )) {
# options_b90$fornetrad <- "globrad"
# stop("Please either provide 'globrad' or 'sunhours' with your climate!")}
# }
if (any(!c("yr", "mo", "da") %in% names(climate))) {
climate$yr <- data.table::year(climate$dates)
climate$mo <- data.table::month(climate$dates)
climate$da <- data.table::mday(climate$dates)
}
if (!any( names(climate) == "mesfl") ) {
climate$mesfl <- 0
}
if (is.null(precip) ){
stopifnot("prec" %in% names(climate))
} else {
names(precip) <- tolower(names(precip))
stopifnot(all(c("dates","prec") %in% tolower(names(precip))))
if (!any( names(precip) == "mesfl") ) {
precip$mesfl <- 0
}
if (nrow(climate)*options_b90$prec_interval != nrow(precip)) {
stop("Climate and Precipitation data provided do not fit to the precipitation
interval defined in options_b90$prec_interval.")
} else {
if (options_b90$prec_interval == 1) {
climate$prec <- precip$prec
precip <- NULL
} else {
precip$ii <- rep(1:options_b90$prec_interval,nrow(climate))
precip$yr <- data.table::year(precip$dates)
precip$mo <- data.table::month(precip$dates)
precip$da <- data.table::mday(precip$dates)
climate$prec <- -999
}
}
}
}))
}
chk_soil <- function(){
eval.parent(quote({
# soil names
if (is.null(soil)) {
if (is.null(param_b90$soil_nodes) | is.null(param_b90$soil_materials)) {
stop("Please provide soil parameters as items 'soil_nodes' and 'soil_materials' via 'param_b90',
when not using 'soil'-argument in run_LWFB90.")
}
names(param_b90$soil_nodes) <- tolower(names(param_b90$soil_nodes))
names(param_b90$soil_materials) <- tolower(names(param_b90$soil_materials))
stopifnot(all(c( "upper", "lower", "mat") %in% names(param_b90$soil_nodes)))
if (anyNA(param_b90$soil_nodes[,c( "upper", "lower", "mat")])) {
stop("No NAs allowed in param_b90$soil_nodes")
}
if (options_b90$imodel == "MvG" ) {
stopifnot(all(c("mat","ths","thr","alpha","npar","ksat","tort","gravel") %in% names(param_b90$soil_materials)))
if (anyNA(param_b90$soil_materials[,c("mat","ths","thr","alpha","npar","ksat","tort","gravel")])) {
stop("No NAs allowed in param_b90$soil_materials!")
}
} else {
stopifnot(all(c("mat","thsat","thetaf","psif","bexp","kf","wetinf","gravel") %in% names(param_b90$soil_materials)))
if (anyNA(param_b90$soil_materials[,c("mat","thsat","thetaf","psif","bexp","kf","wetinf","gravel")])) {
stop("No NAs allowed in param_b90$soil_materials!")
}
}
if (options_b90$root_method == "soilvar" & is.null(param_b90$soil_nodes$rootden)) {
stop("Please provide column 'rootden' in param_b90$soil_nodes when using options_b90$root_method = 'soilvar'.")
}
} else {
names(soil) <- tolower(names(soil))
if (options_b90$imodel == "MvG") {
stopifnot(all(c("upper","lower", "ths","thr","alpha","npar","ksat","tort","gravel") %in% names(soil)))
if (anyNA(soil[,c("upper","lower","ths","thr","alpha","npar","ksat","tort","gravel")])) {
stop("No NAs allowed in 'soil'!")
}
} else {
stopifnot(all(c("upper","lower", "thsat","thetaf","psif","bexp","kf","wetinf","gravel") %in% names(soil)))
if (anyNA(soil[,c("upper","lower","thsat","thetaf","psif","bexp","kf","wetinf","gravel")])) {
stop("No NAs allowed in 'soil'!")
}
}
if (options_b90$root_method == "soilvar" & is.null(soil$rootden)) {
stop("Please provide column 'rootden' in 'soil'-data.frame when using options_b90$root_method = 'soilvar'.")
}
}
}))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.