diff --git a/.Rbuildignore b/.Rbuildignore index e6d99f24..24ec969f 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -8,6 +8,7 @@ ^benchmarking$ ^Rscripts$ ^Rold$ +^.idea$ # R-markdown files in the root folder. ^.*.ods$ diff --git a/COPYING b/COPYING index a57d3c7c..98736771 100644 --- a/COPYING +++ b/COPYING @@ -1,4 +1,4 @@ -Copyright 2023, Geocomputation and Earth Observation, University of Bern +Copyright 2024, Geocomputation and Earth Observation, University of Bern This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. diff --git a/DESCRIPTION b/DESCRIPTION index d9203148..c0e63ddd 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,10 +1,10 @@ Package: rsofun Title: The P-Model and BiomeE Modelling Framework -Version: 4.4.1 +Version: 5.0.0 Authors@R: c( person( - family = "Benjamin", - given = "Stocker", + family = "Stocker", + given = "Benjamin", email = "benjamin.stocker@gmail.com", comment = c(ORCID = "0000-0003-2697-9096"), role = c("aut", "cre")), @@ -26,12 +26,24 @@ Authors@R: c( email = "lauramarqueslopez@gmail.com", comment = c(ORCID = "0000-0002-3593-5557"), role = c("ctb")), + person( + family = "Marcadella", + given = "Mayeul", + email = "mayeul.marcadella@gmail.com", + comment = c(ORCID = "0000-0001-8555-3808"), + role = c("ctb")), person( family = "Weng", given = "Ensheng", email = "wengensheng@gmail.com", comment = c(ORCID = "0000-0002-1858-4847"), role = c("ctb")), + person( + family = "Bernhard", + given = "Fabian", + email = "fabian.bernhard@alumni.epfl.ch", + role = c("aut"), + comment = c(ORCID = "0000-0003-0338-0961")), person(given = "Geocomputation and Earth Observation, University of Bern", role = c("cph", "fnd")) ) diff --git a/NAMESPACE b/NAMESPACE index 088db2a3..35e645ca 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -5,7 +5,6 @@ export(cost_likelihood_biomee) export(cost_likelihood_pmodel) export(cost_rmse_biomee) export(cost_rmse_pmodel) -export(init_dates_dataframe) export(run_biomee_f_bysite) export(run_pmodel_f_bysite) export(runread_biomee_f) diff --git a/NEWS.md b/NEWS.md index 97dcd97b..35492369 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,12 @@ +# rsofun v5.0.0 + +* new BiomeE forcing data matching that of P-model + * `prec` is now called `rain` + * `rh` is now provided as `vpd` + * See `biomee_gs_leuning_drivers` for an example +* fix Fortran modules leading to segmentation faults using BiomeE model +* improved documentation + # rsofun v4.4.1 * bugfix Fortran modules and derived types diff --git a/R/calib_sofun.R b/R/calib_sofun.R index 67a6f904..fcd572fd 100644 --- a/R/calib_sofun.R +++ b/R/calib_sofun.R @@ -92,8 +92,7 @@ calib_sofun <- function( ... ){ # predefine variables for CRAN check compliance - cost <- lower <- upper <- pars <- out <- out_optim <- priors <- setup <- - bt_par <- bt_settings <- NULL + lower <- upper <- out_optim <- NULL # check input variables if(missing(obs) | missing(drivers) | missing(settings)){ @@ -166,12 +165,6 @@ calib_sofun <- function( setup <- BayesianTools::createBayesianSetup( likelihood = function( random_par) { - # cost( - # par = random_par, - # obs = obs, - # drivers = drivers, - # ... - # ) do.call("cost", list( par = random_par, diff --git a/R/cost_likelihood_biomee.R b/R/cost_likelihood_biomee.R index 34249f13..2882c99d 100644 --- a/R/cost_likelihood_biomee.R +++ b/R/cost_likelihood_biomee.R @@ -43,10 +43,10 @@ #' } cost_likelihood_biomee <- function( - par, - obs, - drivers, - targets + par, + obs, + drivers, + targets ){ # predefine variables for CRAN check compliance @@ -57,7 +57,7 @@ cost_likelihood_biomee <- function( drivers$params_species[[1]]$LAI_light[] <- par[2] drivers$params_tile[[1]]$tf_base <- par[3] drivers$params_tile[[1]]$par_mort <- par[4] - + # run model df <- runread_biomee_f( drivers, @@ -67,7 +67,7 @@ cost_likelihood_biomee <- function( # did we spin up spin_up <- drivers$params_siml[[1]]$spinup - + # drop spinup years if activated # see below if (spin_up){ @@ -103,7 +103,7 @@ cost_likelihood_biomee <- function( }) |> unlist() |> sum() # sum log-likelihoods - + # trap boundary conditions if(is.nan(ll) || is.na(ll) | ll == 0){ ll <- -Inf diff --git a/R/cost_likelihood_pmodel.R b/R/cost_likelihood_pmodel.R index 5aac9304..74b16a52 100644 --- a/R/cost_likelihood_pmodel.R +++ b/R/cost_likelihood_pmodel.R @@ -212,7 +212,7 @@ cost_likelihood_pmodel <- function( }) |> unlist() |> sum() - + # trap boundary conditions if(is.nan(ll) | is.na(ll) | ll == 0){ll <- -Inf} diff --git a/R/cost_rmse_pmodel.R b/R/cost_rmse_pmodel.R index 7236bfb0..a9c6dfd9 100644 --- a/R/cost_rmse_pmodel.R +++ b/R/cost_rmse_pmodel.R @@ -73,7 +73,7 @@ cost_rmse_pmodel <- function( targets, par_fixed = NULL, # non-calibrated model parameters target_weights = NULL, # if using several targets, how are the individual - # RMSE weighted? named vector + # # RMSE weighted? named vector parallel = FALSE, ncores = 2 ){ @@ -144,10 +144,10 @@ cost_rmse_pmodel <- function( }else{ # Join P-model output and flux observations df_flux <- df |> - dplyr::filter(sitename %in% flux_sites) |> - dplyr::left_join( - obs_flux, - by = c('sitename', 'date')) # observations with missing date are ignored + dplyr::filter(sitename %in% flux_sites) |> + dplyr::left_join( + obs_flux, + by = c('sitename', 'date')) # observations with missing date are ignored } }else{ df_flux <- data.frame() @@ -170,7 +170,7 @@ cost_rmse_pmodel <- function( df_trait <- df |> dplyr::filter(sitename %in% trait_sites) |> dplyr::group_by(sitename) |> - # get growing season average traits + # get growing season average traits dplyr::summarise(across(ends_with("_mod") & !starts_with('gpp'), ~ sum(.x * gpp_mod/sum(gpp_mod)), .names = "{.col}")) |> @@ -192,7 +192,7 @@ cost_rmse_pmodel <- function( } if(target %in% colnames(df_trait)){ error <- c(error, - (df_trait[[target]] - df_trait[[paste0(target, '_mod')]])^2) + (df_trait[[target]] - df_trait[[paste0(target, '_mod')]])^2) } sqrt(mean(error, na.rm = TRUE)) }) |> @@ -204,6 +204,6 @@ cost_rmse_pmodel <- function( }else{ cost <- mean(rmse, na.rm = TRUE) } - + return(cost) } diff --git a/R/data.R b/R/data.R index fc5c593b..b627d7e6 100644 --- a/R/data.R +++ b/R/data.R @@ -1,21 +1,23 @@ #' rsofun P-model driver data #' -#' Small tests dataset to validate if compiled code -#' and optimization routines can run +#' Small dataset representing the driver to run the P-model at the FR-Pue site. +#' It can also be used together with daily GPP flux time series data from CH-LAE +#' (\code{\link{p_model_validation}}) to optimize model parameters. +#' To optimize model parameters to leaf traits data use the datasets \code{\link{p_model_drivers_vcmax25}} and \code{\link{p_model_validation_vcmax25}}. #' #' @format A tibble of driver data: #' \describe{ #' \item{sitename}{A character string containing the site name.} -#' \item{forcing}{A tibble of a time series of forcing climate data, including -#' the following variables: +#' \item{forcing}{A tibble of a time series of forcing climate data, including +#' the following data: #' \describe{ #' \item{date}{Date of the observation in YYYY-MM-DD format.} #' \item{temp}{Daytime average air temperature in \eqn{^\circ}C.} #' \item{vpd}{Daytime average vapour pressure deficit in Pa.} -#' \item{ppfd}{Photosynthetic photon flux density (PPFD) in +#' \item{ppfd}{Photosynthetic photon flux density (PPFD) in #' mol m\eqn{^{-2}} s\eqn{^{-1}}. If all values are NA, it indicates that #' PPFD should be calculated by the SPLASH model.} -#' \item{netrad}{Net radiation in W m\eqn{^{-2}}. This is currently +#' \item{netrad}{Net radiation in W m\eqn{^{-2}}. This is currently #' ignored as a model forcing.} #' \item{patm}{Atmospheric pressure in Pa.} #' \item{snow}{Snow in water equivalents mm s\eqn{^{-1}}.} @@ -29,7 +31,8 @@ #' net radiation are not prescribed.} #' } #' } -#' \item{params_siml}{A tibble of simulation parameters. +#' \item{params_siml}{A tibble of simulation parameters, including +#' the following data: #' \describe{ #' \item{spinup}{A logical value indicating whether this simulation does spin-up.} #' \item{spinupyears}{Number of spin-up years.} @@ -46,6 +49,7 @@ #' } #' } #' \item{site_info}{A tibble containing site meta information. +#' This data structure can be freely used for documenting the dataset, but must include at least the following data: #' \describe{ #' \item{lon}{Longitude of the site location in degrees east.} #' \item{lat}{Latitude of the site location in degrees north.} @@ -74,10 +78,10 @@ #' surfaces for global land areas. International Journal of Climatology 37 (12): 4302-4315. "p_model_drivers" -#' SOFUN p-model GPP validation data +#' rsofun P-model GPP validation data #' -#' Small tests dataset to validate -#' calibration routines for a time series of fluxes. +#' Small example dataset of target observations (daily GPP flux data) to optimize +#' model parameters with the function \code{\link{calib_sofun}} #' #' @format A tibble of validation data: #' \describe{ @@ -101,60 +105,12 @@ #' rsofun P-model driver data (for leaf traits) #' -#' Small tests dataset to validate if compiled code -#' and optimization routines can run for leaf traits data -#' -#' @format A tibble of driver data: -#' \describe{ -#' \item{sitename}{A character string containing the site name.} -#' \item{forcing}{A tibble of a time series of forcing climate data, including -#' the following variables: -#' \describe{ -#' \item{date}{Date of the observation in YYYY-MM-DD format.} -#' \item{temp}{Daytime average air temperature in \eqn{^\circ}C.} -#' \item{vpd}{Daytime average vapour pressure deficit in Pa.} -#' \item{ppfd}{Photosynthetic photon flux density (PPFD) in -#' mol m\eqn{^{-2}} s\eqn{^{-1}}. If all values are NA, it indicates that -#' PPFD should be calculated by the SPLASH model.} -#' \item{netrad}{Net radiation in W m\eqn{^{-2}}. This is currently -#' ignored as a model forcing.} -#' \item{patm}{Atmospheric pressure in Pa.} -#' \item{snow}{Snow in water equivalents mm s\eqn{^{-1}}.} -#' \item{rain}{Rain as precipitation in liquid form in mm s\eqn{^{-1}}.} -#' \item{tmin}{Daily minimum air temperature in \eqn{^\circ}C.} -#' \item{tmax}{Daily maximum air temperature in \eqn{^\circ}C.} -#' \item{fapar}{Fraction of photosynthetic active radiation (fAPAR), taking -#' values between 0 and 1.} -#' \item{co2}{Atmospheric CO\eqn{_2} concentration.} -#' \item{ccov}{Cloud coverage in \%. This is only used when either PPFD or -#' net radiation are not prescribed.} -#' } -#' } -#' \item{params_siml}{A tibble of simulation parameters. -#' \describe{ -#' \item{spinup}{A logical value indicating whether this simulation does spin-up.} -#' \item{spinupyears}{Number of spin-up years.} -#' \item{recycle}{Length of standard recycling period, in years.} -#' \item{outdt}{An integer indicating the output periodicity.} -#' \item{ltre}{A logical value, \code{TRUE} if evergreen tree.} -#' \item{ltne}{A logical value, \code{TRUE} if evergreen tree and N-fixing.} -#' \item{ltrd}{A logical value, \code{TRUE} if deciduous tree.} -#' \item{ltnd}{A logical value, \code{TRUE} if deciduous tree and N-fixing.} -#' \item{lgr3}{A logical value, \code{TRUE} if grass with C3 photosynthetic pathway.} -#' \item{lgn3}{A logical value, \code{TRUE} if grass with C3 photosynthetic -#' pathway and N-fixing.} -#' \item{lgr4}{A logical value, \code{TRUE} if grass with C4 photosynthetic pathway.} -#' } -#' } -#' \item{site_info}{A tibble containing site meta information. -#' \describe{ -#' \item{lon}{Longitude of the site location in degrees east.} -#' \item{lat}{Latitude of the site location in degrees north.} -#' \item{elv}{Elevation of the site location, in meters above sea level.} -#' \item{whc}{A numeric value for the rooting zone water holding capacity (in mm)} -#' } -#' } -#' } +#' Small dataset representing the driver to run the P-model at four separate sites. +#' It can also be used together with leaf traits data from these four sites +#' (\code{\link{p_model_validation_vcmax25}}) to optimize model parameters. +#' To optimize model parameters to GPP flux data use the datasets \code{\link{p_model_drivers}} and \code{\link{p_model_validation}}. +#' +#' @format See \code{\link{p_model_drivers}} #' #' @source Atkin, O. K., Bloomfield, K. J., Reich, P. B., Tjoelker, M. G., Asner, G. P., Bonal, D., et al. (2015). #' Global variability in leaf respiration in relation to climate, plant functional types and leaf traits. @@ -179,10 +135,10 @@ #' Atmospheric carbon dioxide variations at Mauna Loa Observatory, Hawaii, Tellus, vol. 28, 538-551 "p_model_drivers_vcmax25" -#' SOFUN p-model Vcmax25 validation data +#' rsofun P-model Vcmax25 validation data #' -#' Small tests dataset to validate -#' calibration routines for leaf traits. +#' Small example dataset of target observations (leaf trait data) to optimize +#' model parameters with the function \code{\link{calib_sofun}} #' #' @format A tibble of validation data: #' \describe{ @@ -207,42 +163,180 @@ #' New Phytol. 206 (2), 614–636. doi:10.1111/nph.13253 "p_model_validation_vcmax25" -#' rsofun BiomeE driver data +#' rsofun P-model output data #' -#' Small tests dataset to validate if compiled code -#' and optimization routines can run using the -#' p-model specifications -#' -#' @format A tibble of driver data: -#' \describe{ -#' \item{sitename}{site name} -#' \item{params_siml}{model parameters} -#' \item{site_info}{site information} -#' \item{soil_texture}{soil texture data} -#' \item{forcing}{forcing data} -#' } -"biomee_p_model_drivers" +#' Example output dataset from a p-model run using \code{\link{p_model_drivers}} +"p_model_output" -#' rsofun BiomeE driver data +#' rsofun P-model output data (using vcmax25 drivers) #' -#' Small tests dataset to validate if compiled code -#' and optimization routines can run using the -#' Leuning specifications +#' Example output dataset from a p-model run using \code{\link{p_model_drivers_vcmax25}} +"p_model_output_vcmax25" + +#' rsofun BiomeE driver data (Leuning photosynthesis model) +#' +#' Small dataset representing the driver to run the BiomeE-model at the CH-LAE site +#' using the Leuning photosynthesis specification (and half-hourly time step) +#' It can also be used together with leaf trait data from CH-LAE (\code{\link{biomee_validation}}) +#' to optimize model parameters. #' -#' @format A tibble of driver data: +#' @format A tibble of driver data. #' \describe{ -#' \item{sitename}{site name} -#' \item{params_siml}{model parameters} -#' \item{site_info}{site information} -#' \item{soil_texture}{soil texture data} -#' \item{forcing}{forcing data} +#' \item{sitename}{Site name} +#' \item{params_siml}{Simulation parameters as a data.frame, including +#' the following data: +#' \describe{ +#' \item{spinup}{Flag indicating whether this simulation does spin-up.} +#' \item{spinupyears}{Number of spin-up years.} +#' \item{recycle}{Length of standard recycling period (years).} +#' \item{firstyeartrend}{First transient year.} +#' \item{nyeartrend}{Number of transient years.} +#' \item{steps_per_day}{Time resolution (day-1).} +#' \item{do_U_shaped_mortality}{Flag indicating whether U-shaped +#' mortality is used.} +#' \item{update_annualLAImax}{Flag indicating whether updating +#' LAImax according to mineral N in soil.} +#' \item{do_closedN_run}{Flag indicating whether doing N closed +#' runs to recover N balance enforcing 0.2 kg N m-2 in the inorganic N pool.} +#' \item{code_method_photosynth}{String specifying the method of photosynthesis +#' used in the model, either "pmodel" or "gs_leuning".document()} +#' \item{code_method_mortality}{String indicating the type of mortality in the +#' model. One of the following: "dbh" is size-dependent mortality, "const_selfthin" +#' is constant self thinning (in development), "cstarvation" is carbon starvation, and +#' "growthrate" is growth rate dependent mortality.} +#' }} +#' \item{site_info}{Site meta info in a data.frame. +#' This data structure can be freely used for documenting the dataset, but must include at least the following data: +#' \describe{ +#' \item{lon}{Longitude of the site location.} +#' \item{lat}{Latitude of the site location.} +#' \item{elv}{Elevation of the site location, in meters.} +#' }} +#' \item{forcing}{Forcing data.frame used as input +#' \describe{ +#' \item{ppfd}{Photosynthetic photon flux density (mol s-1 m-2)} +#' \item{tair}{Air temperature (deg C)} +#' \item{vpd}{Vapor pressure deficit (Pa)} +#' \item{rain}{Precipitation (kgH2O m-2 s-1 == mm s-1)} +#' \item{wind}{Wind velocity (m s-1)} +#' \item{pair}{Atmospheric pressure (pa)} +#' \item{co2}{CO2 atmospheric concentration (ppm)} +#' }} +#' \item{params_tile}{Tile-level model parameters, into a single row data.frame, including +#' the following data: +#' \describe{ +#' \item{soiltype}{Integer indicating the type of soil: Sand = 1, LoamySand = 2, +#' SandyLoam = 3, SiltLoam = 4, FrittedClay = 5, Loam = 6, Clay = 7.} +#' \item{FLDCAP}{Field capacity (vol/vol). Water remaining in a soil after it +#' has been thoroughly saturated and allowed to drain freely.} +#' \item{WILTPT}{Wilting point (vol/vol). Water content of a soil at which +#' plants wilt and fail to recover.} +#' \item{K1}{Fast soil C decomposition rate (year\eqn{^{-1}}).} +#' \item{K2}{Slow soil C decomposition rate (year\eqn{^{-1}}).} +#' \item{K_nitrogen}{Mineral Nitrogen turnover rate (year\eqn{^{-1}}).} +#' \item{MLmixRatio}{Ratio of C and N returned to litters from microbes.} +#' \item{etaN}{N loss rate through runoff (organic and mineral) (year\eqn{^{-1}}).} +#' \item{LMAmin}{Minimum LMA, leaf mass per unit area, kg C m\eqn{^{-2}}.} +#' \item{fsc_fine}{Fraction of fast turnover carbon in fine biomass.} +#' \item{fsc_wood}{Fraction of fast turnover carbon in wood biomass.} +#' \item{GR_factor}{Growth respiration factor.} +#' \item{l_fract}{Fraction of the carbon retained after leaf drop.} +#' \item{retransN}{Retranslocation coefficient of nitrogen.} +#' \item{f_initialBSW}{Coefficient for setting up initial sapwood.} +#' \item{f_N_add}{Re-fill of N for sapwood.} +#' \item{tf_base}{Calibratable scalar for respiration, used to increase LUE +#' levels.} +#' \item{par_mort}{Canopy mortality parameter.} +#' \item{par_mort_under}{Parameter for understory mortality.} +#' }} +#' \item{params_species}{A data.frame containing species-specific model parameters, +#' with one species per row, including the following data: +#' \describe{ +#' \item{lifeform}{Integer set to 0 for grasses and 1 for trees.} +#' \item{phenotype}{Integer set to 0 for deciduous and 1 for evergreen.} +#' \item{pt}{Integer indicating the type of plant according to photosynthesis: +#' 0 for C3; 1 for C4} +#' \item{alpha_FR}{Fine root turnonver rate (year\eqn{^{-1}}).} +#' \item{rho_FR}{Material density of fine roots (kg C m\eqn{^{-3}}).} +#' \item{root_r}{Radius of the fine roots, in m.} +#' \item{root_zeta}{e-folding parameter of root vertical distribution, in m.} +#' \item{Kw_root}{Fine root water conductivity (mol m\eqn{^{-2}} +#' s\eqn{^{-1}} MPa\eqn{^{-1}}).} +#' \item{leaf_size}{Characteristic leaf size.} +#' \item{Vmax}{Max RuBisCo rate, in mol m\eqn{^{-2}} s\eqn{^{-1}}.} +#' \item{Vannual}{Annual productivity per unit area at full sun (kg C +#' m\eqn{^{-2}} year\eqn{^{-2}}).} +#' \item{wet_leaf_dreg}{Wet leaf photosynthesis down-regulation.} +#' \item{m_cond}{Factor of stomatal conductance.} +#' \item{alpha_phot}{Photosynthesis efficiency.} +#' \item{gamma_L}{Leaf respiration coefficient, in year\eqn{^{-1}}.} +#' \item{gamma_LN}{Leaf respiration coefficient per unit N.} +#' \item{gamma_SW}{Sapwood respiration rate, in kg C m\eqn{^{-2}} year\eqn{^{-1}}.} +#' \item{gamma_FR}{Fine root respiration rate, kg C kg C\eqn{^{-1}} +#' year\eqn{^{-1}}.} +#' \item{tc_crit}{Critical temperature triggerng offset of phenology, in Kelvin.} +#' \item{tc_crit_on}{Critical temperature triggerng onset of phenology, in Kelvin.} +#' \item{gdd_crit}{Critical value of GDD5 for turning ON growth season.} +#' \item{betaON}{Critical soil moisture for phenology onset.} +#' \item{betaOFF}{Critical soil moisture for phenology offset.} +#' \item{seedlingsize}{Initial size of seedlings, in kg C per individual.} +#' \item{LNbase}{Basal leaf N per unit area, in kg N m\eqn{^{-2}}.} +#' \item{lAImax}{Maximum crown LAI (leaf area index).} +#' \item{Nfixrate0}{Reference N fixation rate (kg N kg C\eqn{^{-1}} root).} +#' \item{NfixCost0}{Carbon cost of N fixation (kg C kg N\eqn{^{-1}}).} +#' \item{phiCSA}{Ratio of sapwood area to leaf area.} +#' \item{mortrate_d_c}{Canopy tree mortality rate (year\eqn{^{-1}}).} +#' \item{mortrate_d_u}{Understory tree mortality rate (year\eqn{^{-1}}).} +#' \item{maturalage}{Age at which trees can reproduce (years).} +#' \item{v_seed}{Fraction of G_SF to G_F.} +#' \item{fNSmax}{Multiplier for NSNmax as sum of potential bl and br.} +#' \item{LMA}{Leaf mass per unit area (kg C m\eqn{^{-2}}).} +#' \item{rho_wood}{Wood density (kg C m\eqn{^{-3}}).} +#' \item{alphaBM}{Coefficient for allometry (biomass = alphaBM * DBH ** thetaBM).} +#' \item{thetaBM}{Coefficient for allometry (biomass = alphaBM * DBH ** thetaBM).} +#' \item{kphio}{Quantum yield efficiency \eqn{\varphi_0}, +#' in mol mol\eqn{^{-1}}.} +#' \item{phiRL}{Ratio of fine root to leaf area.} +#' \item{LAI_light}{Maximum LAI limited by light.} +#' }} +#' \item{init_cohort}{A data.frame of initial cohort specifications, including +#' the following data: +#' \describe{ +#' \item{init_cohort_species}{Index of a species described in param_species.} +#' \item{init_cohort_nindivs}{Initial individual density, in individuals per +#' m\eqn{^{2}}.} +#' \item{init_cohort_bsw}{Initial biomass of sapwood, in kg C per individual.} +#' \item{init_cohort_bHW}{Initial biomass of heartwood, in kg C per tree.} +#' \item{init_cohort_nsc}{Initial non-structural biomass.} +#' }} +#' \item{init_soil}{A data.frame of initial soil pools, including +#' the following data: +#' \describe{ +#' \item{init_fast_soil_C}{Initial fast soil carbon, in kg C m\eqn{^{-2}}.} +#' \item{init_slow_soil_C}{Initial slow soil carbon, in kg C m\eqn{^{-2}}.} +#' \item{init_Nmineral}{Mineral nitrogen pool, in kg N m\eqn{^{-2}}.} +#' \item{N_input}{Annual nitrogen input to soil N pool, in kg N m\eqn{^{-2}} +#' year\eqn{^{-1}}.} +#' }} #' } "biomee_gs_leuning_drivers" +#' rsofun BiomeE driver data (P-model photosynthesis model) +#' +#' Small dataset representing the driver to run the BiomeE-model at the CH-LAE site +#' using the P-model photosynthesis specification (and daily time step). +#' It can also be used together with leaf trait data from CH-LAE (\code{\link{biomee_validation}}) +#' to optimize model parameters. +#' +#' @format See \code{\link{biomee_gs_leuning_drivers}} +#' +#' @inherit biomee_gs_leuning_drivers source +"biomee_p_model_drivers" + #' rsofun BiomeE targets validation data #' -#' Small tests dataset to validate -#' calibration routines +#' Small example dataset of target observations (leaf trait data) at the CH-LAE site +#' to optimize model parameters with the function \code{\link{calib_sofun}} #' #' @format A tibble of validation data: #' \describe{ @@ -255,13 +349,12 @@ #' Dataset. https://doi.org/10.18140/FLX/1440134 "biomee_validation" -#' rsofun BiomeE (p-model) output data +#' rsofun BiomeE (P-model) output data #' -#' Example output dataset using BiomeE (p-model) +#' Example output dataset from a BiomeE-model run (p-model) "biomee_p_model_output" -#' rsofun BiomeE (GS leuning) output data +#' rsofun BiomeE (gs_leuning) output data #' -#' Example output dataset using BiomeE (GS leuning) +#' Example output dataset from a BiomeE-model run (gs_leuning) "biomee_gs_leuning_output" - diff --git a/R/init_dates_dataframe.R b/R/init_dates_dataframe.R index 1192075a..932b79a1 100644 --- a/R/init_dates_dataframe.R +++ b/R/init_dates_dataframe.R @@ -3,6 +3,8 @@ #' Creates a tibble with rows for each date from \code{'yrstart'} to \code{'yrend'} #' in \code{'yyyy-mm-dd'} format. Intervals of dates are specified by argument #'\code{'freq'}. +#' ddf <- init_dates_dataframe(2000, 2003, startmoy=1, startdoy=1, +#' freq="days", endmoy=12, enddom=31, noleap=FALSE) #' #' @param yrstart An integer defining the start year #' of dates covered by the dataframe. @@ -24,21 +26,17 @@ #' of February is removed. Defaults to \code{FALSE}. #' #' @return A tibble with dates. -#' @export #' -#' @examples -#' ddf <- init_dates_dataframe( 2000, 2003, startmoy=1, startdoy=1, -#' freq="days", endmoy=12, enddom=31, noleap=FALSE ) init_dates_dataframe <- function( - yrstart, - yrend, - startmoy=1, - startdoy=1, - freq="days", - endmoy=12, - enddom=31, - noleap=FALSE ){ + yrstart, + yrend, + startmoy=1, + startdoy=1, + freq="days", + endmoy=12, + enddom=31, + noleap=FALSE ){ if (freq=="days"){ @@ -49,7 +47,7 @@ init_dates_dataframe <- function( end_date <- as.Date( sprintf("%04d-%02d-%02d", yrend, endmoy, enddom)) - + } else if (freq=="months"){ start_date <- as.Date( @@ -59,7 +57,7 @@ init_dates_dataframe <- function( end_date <- as.Date( sprintf("%04d-%02d-15", yrend, endmoy)) - + } else if (freq=="years"){ start_date <- as.Date( @@ -70,7 +68,7 @@ init_dates_dataframe <- function( sprintf("%04d-%02d-01", yrend, 7)) } - + # define date range date_range <- data.frame( date = seq.Date( @@ -81,14 +79,14 @@ init_dates_dataframe <- function( # convert to decimal date date_range$year_dec <- numeric_year(date_range$date) - + # leap year filter if (noleap) { date_range <- dplyr::filter(date_range, - !(format(date, "%m-%d") == "02-29") - ) + !(format(date, "%m-%d") == "02-29") + ) - } - + } + return(date_range) } diff --git a/R/run_biomee_f_bysite.R b/R/run_biomee_f_bysite.R index c10932b5..c48faa1d 100644 --- a/R/run_biomee_f_bysite.R +++ b/R/run_biomee_f_bysite.R @@ -4,137 +4,19 @@ #' #' @param sitename Site name. #' @param params_siml Simulation parameters. -#' \describe{ -#' \item{spinup}{Flag indicating whether this simulation does spin-up.} -#' \item{spinupyears}{Number of spin-up years.} -#' \item{recycle}{Length of standard recycling period (years).} -#' \item{firstyeartrend}{First transient year.} -#' \item{nyeartrend}{Number of transient years.} -#' \item{steps_per_day}{Time resolution (day-1).} -#' \item{do_U_shaped_mortality}{Flag indicating whether U-shaped -#' mortality is used.} -#' \item{update_annualLAImax}{Flag indicating whether updating -#' LAImax according to mineral N in soil.} -#' \item{do_closedN_run}{Flag indicating whether doing N closed -#' runs to recover N balance enforcing 0.2 kg N m-2 in the inorganic N pool.} -#' \item{code_method_photosynth}{String specifying the method of photosynthesis -#' used in the model, either "pmodel" or "gs_leuning".document()} -#' \item{code_method_mortality}{String indicating the type of mortality in the -#' model. One of the following: "dbh" is size-dependent mortality, "const_selfthin" -#' is constant self thinning (in development), "cstarvation" is carbon starvation, and -#' "growthrate" is growth rate dependent mortality.} -#' } +#' See examples \code{\link{biomee_gs_leuning_drivers}} or \code{\link{biomee_p_model_drivers}} #' @param site_info Site meta info in a data.frame. -#' \describe{ -#' \item{sitename}{Name of the site.} -#' \item{lon}{Longitud of the site location.} -#' \item{lat}{Latitude of the site location.} -#' \item{elv}{Elevation of the site location, in meters.} -#' } +#' See examples \code{\link{biomee_gs_leuning_drivers}} or \code{\link{biomee_p_model_drivers}} #' @param forcing Forcing data.frame used as input. -#' \describe{ -#' \item{ppfd}{Photosynthetic photon flux densisty(mol s-1 m-2)} -#' \item{tair}{Air temperature (deg C)} -#' \item{vpd}{Vapor pressure defficit (Pa)} -#' \item{rain}{Precipitation (kgH2O m-2 s-1 == mm s-1)} -#' \item{wind}{Wind velocity (m s-1)} -#' \item{pair}{Atmoshperic pressure (pa)} -#' \item{co2}{CO2 athmospheric concentration (ppm)} -#' } -#' @param params_tile Tile-level model parameters, into a single row data.frame -#' with columns: -#' \describe{ -#' \item{soiltype}{Integer indicating the type of soil: Sand = 1, LoamySand = 2, -#' SandyLoam = 3, SiltLoam = 4, FrittedClay = 5, Loam = 6, Clay = 7.} -#' \item{FLDCAP}{Field capacity (vol/vol). Water remaining in a soil after it -#' has been thoroughly saturated and allowed to drain freely.} -#' \item{WILTPT}{Wilting point (vol/vol). Water content of a soil at which -#' plants wilt and fail to recover.} -#' \item{K1}{Fast soil C decomposition rate (year\eqn{^{-1}}).} -#' \item{K2}{Slow soil C decomposition rate (year\eqn{^{-1}}).} -#' \item{K_nitrogen}{Mineral Nitrogen turnover rate (year\eqn{^{-1}}).} -#' \item{MLmixRatio}{Ratio of C and N returned to litters from microbes.} -#' \item{etaN}{N loss rate through runoff (organic and mineral) (year\eqn{^{-1}}).} -#' \item{LMAmin}{Minimum LMA, leaf mass per unit area, kg C m\eqn{^{-2}}.} -#' \item{fsc_fine}{Fraction of fast turnover carbon in fine biomass.} -#' \item{fsc_wood}{Fraction of fast turnover carbon in wood biomass.} -#' \item{GR_factor}{Growth respiration factor.} -#' \item{l_fract}{Fraction of the carbon retained after leaf drop.} -#' \item{retransN}{Retranslocation coefficient of nitrogen.} -#' \item{f_initialBSW}{Coefficient for setting up initial sapwood.} -#' \item{f_N_add}{Re-fill of N for sapwood.} -#' \item{tf_base}{Calibratable scalar for respiration, used to increase LUE -#' levels.} -#' \item{par_mort}{Canopy mortality parameter.} -#' \item{par_mort_under}{Parameter for understory mortality.} -#' } +#' See examples \code{\link{biomee_gs_leuning_drivers}} or \code{\link{biomee_p_model_drivers}} +#' @param params_tile Tile-level model parameters, into a single row data.frame. +#' See examples \code{\link{biomee_gs_leuning_drivers}} or \code{\link{biomee_p_model_drivers}} #' @param params_species A data.frame containing species-specific model parameters, -#' with one species per row. The columns of this data.frame are: -#' \describe{ -#' \item{lifeform}{Integer set to 0 for grasses and 1 for trees.} -#' \item{phenotype}{Integer set to 0 for deciduous and 1 for evergreen.} -#' \item{pt}{Integer indicating the type of plant according to photosynthesis: -#' 0 for C3; 1 for C4} -#' \item{alpha_FR}{Fine root turnonver rate (year\eqn{^{-1}}).} -#' \item{rho_FR}{Material density of fine roots (kg C m\eqn{^{-3}}).} -#' \item{root_r}{Radious of the fine roots, in m.} -#' \item{root_zeta}{e-folding parameter of root vertical distribution, in m.} -#' \item{Kw_root}{Fine root water conductivity (mol m\eqn{^{-2}} -#' s\eqn{^{-1}} MPa\eqn{^{-1}}).} -#' \item{leaf_size}{Characteristic leaf size.} -#' \item{Vmax}{Max RuBisCo rate, in mol m\eqn{^{-2}} s\eqn{^{-1}}.} -#' \item{Vannual}{Annual productivity per unit area at full sun (kg C -#' m\eqn{^{-2}} year\eqn{^{-2}}).} -#' \item{wet_leaf_dreg}{Wet leaf photosynthesis down-regulation.} -#' \item{m_cond}{Factor of stomatal conductance.} -#' \item{alpha_phot}{Photosynthesis efficiency.} -#' \item{gamma_L}{Leaf respiration coefficient, in year\eqn{^{-1}}.} -#' \item{gamma_LN}{Leaf respiration coefficient per unit N.} -#' \item{gamma_SW}{Sapwood respiration rate, in kg C m\eqn{^{-2}} year\eqn{^{-1}}.} -#' \item{gamma_FR}{Fine root respiration rate, kg C kg C\eqn{^{-1}} -#' year\eqn{^{-1}}.} -#' \item{tc_crit}{Critical temperature triggerng offset of phenology, in Kelvin.} -#' \item{tc_crit_on}{Critical temperature triggerng onset of phenology, in Kelvin.} -#' \item{gdd_crit}{Critical value of GDD5 for turning ON growth season.} -#' \item{betaON}{Critical soil moisture for phenology onset.} -#' \item{betaOFF}{Critical soil moisture for phenology offset.} -#' \item{seedlingsize}{Initial size of seedlings, in kg C per individual.} -#' \item{LNbase}{Basal leaf N per unit area, in kg N m\eqn{^{-2}}.} -#' \item{lAImax}{Maximum crown LAI (leaf area index).} -#' \item{Nfixrate0}{Reference N fixation rate (kg N kg C\eqn{^{-1}} root).} -#' \item{NfixCost0}{Carbon cost of N fixation (kg C kg N\eqn{^{-1}}).} -#' \item{phiCSA}{Ratio of sapwood area to leaf area.} -#' \item{mortrate_d_c}{Canopy tree mortality rate (year\eqn{^{-1}}).} -#' \item{mortrate_d_u}{Understory tree mortality rate (year\eqn{^{-1}}).} -#' \item{maturalage}{Age at which trees can reproduce (years).} -#' \item{v_seed}{Fraction of G_SF to G_F.} -#' \item{fNSmax}{Multiplier for NSNmax as sum of potential bl and br.} -#' \item{LMA}{Leaf mass per unit area (kg C m\eqn{^{-2}}).} -#' \item{rho_wood}{Wood density (kg C m\eqn{^{-3}}).} -#' \item{alphaBM}{Coefficient for allometry (biomass = alphaBM * DBH ** thetaBM).} -#' \item{thetaBM}{Coefficient for allometry (biomass = alphaBM * DBH ** thetaBM).} -#' \item{kphio}{Quantum yield efficiency \eqn{\varphi_0}, -#' in mol mol\eqn{^{-1}}.} -#' \item{phiRL}{Ratio of fine root to leaf area.} -#' \item{LAI_light}{Maximum LAI limited by light.} -#' } +#' with one species per row. See examples \code{\link{biomee_gs_leuning_drivers}} or \code{\link{biomee_p_model_drivers}} #' @param init_cohort A data.frame of initial cohort specifications. -#' \describe{ -#' \item{init_cohort_species}{Index of a species described in param_species.} -#' \item{init_cohort_nindivs}{Initial individual density, in individuals per -#' m\eqn{^{2}}.} -#' \item{init_cohort_bsw}{Initial biomass of sapwood, in kg C per individual.} -#' \item{init_cohort_bHW}{Initial biomass of heartwood, in kg C per tree.} -#' \item{init_cohort_nsc}{Initial non-structural biomass.} -#' } +#' See examples \code{\link{biomee_gs_leuning_drivers}} or \code{\link{biomee_p_model_drivers}} #' @param init_soil A data.frame of initial soil pools. -#' \describe{ -#' \item{init_fast_soil_C}{Initial fast soil carbon, in kg C m\eqn{^{-2}}.} -#' \item{init_slow_soil_C}{Initial slow soil carbon, in kg C m\eqn{^{-2}}.} -#' \item{init_Nmineral}{Mineral nitrogen pool, in kg N m\eqn{^{-2}}.} -#' \item{N_input}{Annual nitrogen input to soil N pool, in kg N m\eqn{^{-2}} -#' year\eqn{^{-1}}.} -#' } +#' See examples \code{\link{biomee_gs_leuning_drivers}} or \code{\link{biomee_p_model_drivers}} #' @param makecheck Flag specifying whether checks are performed to verify model inputs and parameters. #' #' @export @@ -365,20 +247,20 @@ #' } run_biomee_f_bysite <- function( - sitename, - params_siml, - site_info, - forcing, - params_tile, - params_species, - init_cohort, - init_soil, - makecheck = TRUE - ){ + sitename, + params_siml, + site_info, + forcing, + params_tile, + params_species, + init_cohort, + init_soil, + makecheck = TRUE +){ # predefine variables for CRAN check compliance type <- NULL - + forcing_features <- c( 'ppfd', 'temp', @@ -394,14 +276,14 @@ run_biomee_f_bysite <- function( select( any_of(forcing_features) ) - + runyears <- ifelse( params_siml$spinup, (params_siml$spinupyears + params_siml$nyeartrend), params_siml$nyeartrend) n_daily <- params_siml$nyeartrend * 365 - + # Types of photosynthesis model if (params_siml$method_photosynth == "gs_leuning"){ code_method_photosynth <- 1 @@ -417,17 +299,17 @@ run_biomee_f_bysite <- function( stop( "run_biomee_f_bysite: time step must be daily for P-model photosynthesis setup." - ) - } + ) + } } else { stop( paste("run_biomee_f_bysite: params_siml$method_photosynth not recognised:", params_siml$method_photosynth)) } - -# Types of mortality formulations - if (params_siml$method_mortality == "cstarvation"){ + + # Types of mortality formulations + if (params_siml$method_mortality == "cstarvation"){ code_method_mortality <- 1 } else if (params_siml$method_mortality == "growthrate"){ code_method_mortality <- 2 @@ -442,10 +324,10 @@ run_biomee_f_bysite <- function( paste("run_biomee_f_bysite: params_siml$method_mortality not recognised:", params_siml$method_mortality)) } - + # base state, always execute the call continue <- TRUE - + # validate input if (makecheck){ # Add input and parameter checks here if applicable. @@ -461,19 +343,19 @@ run_biomee_f_bysite <- function( return(TRUE) } }) - + # only return true if all checked variables are TRUE # suppress warning on coercion of list to single logical continue <- suppressWarnings(all(as.vector(data_integrity))) } - + if (continue) { - + ## C wrapper call biomeeout <- .Call( - + 'biomee_f_C', - + ## Simulation parameters spinup = as.logical(params_siml$spinup), spinupyears = as.integer(params_siml$spinupyears), @@ -490,7 +372,7 @@ run_biomee_f_bysite <- function( longitude = as.numeric(site_info$lon), latitude = as.numeric(site_info$lat), altitude = as.numeric(site_info$elv), - + ## Tile-level parameters soiltype = as.integer(params_tile$soiltype), FLDCAP = as.numeric(params_tile$FLDCAP), @@ -511,15 +393,15 @@ run_biomee_f_bysite <- function( tf_base = as.numeric(params_tile$tf_base), par_mort = as.numeric(params_tile$par_mort), par_mort_under = as.numeric(params_tile$par_mort_under), - + ## Species-specific parameters n_params_species = as.integer(nrow(params_species)), params_species = as.matrix(params_species), - + ## initial cohort sizes n_init_cohort = as.integer(nrow(init_cohort)), init_cohort = as.matrix(init_cohort), - + ## initial soil pools init_fast_soil_C = as.numeric(init_soil$init_fast_soil_C), init_slow_soil_C = as.numeric(init_soil$init_slow_soil_C), @@ -544,9 +426,9 @@ run_biomee_f_bysite <- function( format( utils::object.size(biomeeout), units = "GB" - ) ) ) + ) if (size_of_object_gb >= 5){ warning( @@ -605,66 +487,66 @@ run_biomee_f_bysite <- function( # annual tile output_annual_tile <- as.data.frame(biomeeout[[2]], stringAsFactor = FALSE) colnames(output_annual_tile) <- c( - "year", - "CAI", - "LAI", - "Density", - "DBH", - "Density12", - "DBH12", - "QMD12", - "NPP", - "GPP", - "Rauto", - "Rh", - "rain", - "SoilWater", - "Transp", - "Evap", - "Runoff", - "plantC", - "soilC", - "plantN", - "soilN", - "totN", - "NSC", - "SeedC", - "leafC", - "rootC", - "SapwoodC", - "WoodC", - "NSN", - "SeedN", - "leafN", - "rootN", - "SapwoodN", - "WoodN", - "McrbC", - "fastSOM", - "SlowSOM", - "McrbN", - "fastSoilN", - "slowSoilN", - "mineralN", - "N_fxed", - "N_uptk", - "N_yrMin", - "N_P2S", - "N_loss", - "totseedC", - "totseedN", - "Seedling_C", - "Seedling_N", - "MaxAge", - "MaxVolume", - "MaxDBH", - "NPPL", - "NPPW", - "n_deadtrees", - "c_deadtrees", - "m_turnover", - "c_turnover_time" - ) + "year", + "CAI", + "LAI", + "Density", + "DBH", + "Density12", + "DBH12", + "QMD12", + "NPP", + "GPP", + "Rauto", + "Rh", + "rain", + "SoilWater", + "Transp", + "Evap", + "Runoff", + "plantC", + "soilC", + "plantN", + "soilN", + "totN", + "NSC", + "SeedC", + "leafC", + "rootC", + "SapwoodC", + "WoodC", + "NSN", + "SeedN", + "leafN", + "rootN", + "SapwoodN", + "WoodN", + "McrbC", + "fastSOM", + "SlowSOM", + "McrbN", + "fastSoilN", + "slowSoilN", + "mineralN", + "N_fxed", + "N_uptk", + "N_yrMin", + "N_P2S", + "N_loss", + "totseedC", + "totseedN", + "Seedling_C", + "Seedling_N", + "MaxAge", + "MaxVolume", + "MaxDBH", + "NPPL", + "NPPW", + "n_deadtrees", + "c_deadtrees", + "m_turnover", + "c_turnover_time" + ) #--- annual cohorts ---- annual_values <- c( @@ -720,7 +602,7 @@ run_biomee_f_bysite <- function( # drop rows (cohorts) with no values output_annual_cohorts$year[output_annual_cohorts$year == -9999 | - output_annual_cohorts$year == 0] <- NA + output_annual_cohorts$year == 0] <- NA output_annual_cohorts <- output_annual_cohorts[!is.na(output_annual_cohorts$year),] @@ -735,9 +617,9 @@ run_biomee_f_bysite <- function( } else { out <- NA } - + return(out) - + } .onUnload <- function(libpath) { diff --git a/R/run_pmodel_f_bysite.R b/R/run_pmodel_f_bysite.R index 8ebdbbdd..9e4367c8 100644 --- a/R/run_pmodel_f_bysite.R +++ b/R/run_pmodel_f_bysite.R @@ -20,7 +20,7 @@ #' } #' @param site_info A list of site meta info. Required: #' \describe{ -#' \item{lon}{Longitud of the site location.} +#' \item{lon}{Longitude of the site location.} #' \item{lat}{Latitude of the site location.} #' \item{elv}{Elevation of the site location, in meters.} #' \item{whc}{A numeric value for the total root zone water holding capacity (in mm), used @@ -155,20 +155,20 @@ #' ) run_pmodel_f_bysite <- function( - sitename, - params_siml, - site_info, - forcing, - params_modl, - makecheck = TRUE, - verbose = TRUE - ){ + sitename, + params_siml, + site_info, + forcing, + params_modl, + makecheck = TRUE, + verbose = TRUE +){ # predefine variables for CRAN check compliance ccov <- temp <- rain <- vpd <- ppfd <- netrad <- - fsun <- snow <- co2 <- fapar <- patm <- - nyeartrend_forcing <- firstyeartrend_forcing <- - tmin <- tmax <- . <- NULL + fsun <- snow <- co2 <- fapar <- patm <- + nyeartrend_forcing <- firstyeartrend_forcing <- + tmin <- tmax <- . <- NULL # base state, always execute the call continue <- TRUE @@ -211,7 +211,7 @@ run_pmodel_f_bysite <- function( patm, tmin, tmax - ) + ) # validate input if (makecheck){ @@ -235,7 +235,7 @@ run_pmodel_f_bysite <- function( data_integrity <- lapply(check_vars, function(check_var){ if (any(is.nanull(forcing[check_var]))){ warning(sprintf("Error: Missing value in %s for %s", - check_var, sitename)) + check_var, sitename)) return(FALSE) } else { return(TRUE) @@ -263,7 +263,7 @@ run_pmodel_f_bysite <- function( parameter_integrity <- lapply(check_param, function(check_var){ if (any(is.nanull(params_siml[check_var]))){ warning(sprintf("Error: Missing value in %s for %s", - check_var, sitename)) + check_var, sitename)) return(FALSE) } else { return(TRUE) @@ -273,7 +273,7 @@ run_pmodel_f_bysite <- function( if (suppressWarnings(!all(parameter_integrity))){ continue <- FALSE } - + if (nrow(forcing) %% ndayyear != 0){ # something weird more fundamentally -> don't run the model warning(" Returning a dummy data frame. Forcing data does not @@ -283,9 +283,9 @@ run_pmodel_f_bysite <- function( # Check model parameters if( sum( names(params_modl) %in% c('kphio', 'kphio_par_a', 'kphio_par_b', - 'soilm_thetastar', 'soilm_betao', - 'beta_unitcostratio', 'rd_to_vcmax', - 'tau_acclim', 'kc_jmax') + 'soilm_thetastar', 'soilm_betao', + 'beta_unitcostratio', 'rd_to_vcmax', + 'tau_acclim', 'kc_jmax') ) != 9){ warning(" Returning a dummy data frame. Incorrect model parameters.") continue <- FALSE @@ -293,14 +293,14 @@ run_pmodel_f_bysite <- function( } if (continue){ - + # determine whether to read PPFD from forcing or to calculate internally in_ppfd <- ifelse(any(is.na(forcing$ppfd)), FALSE, TRUE) - + # determine whether to read PPFD from forcing or to calculate internally # in_netrad <- ifelse(any(is.na(forcing$netrad)), FALSE, TRUE) in_netrad <- FALSE # net radiation is currently ignored as a model forcing, but is internally simulated by SPLASH. - + # Check if fsun is available if(! (in_ppfd & in_netrad)){ # fsun must be available when one of ppfd or netrad is missing @@ -315,7 +315,7 @@ run_pmodel_f_bysite <- function( # number of rows in matrix (pre-allocation of memory) n <- as.integer(nrow(forcing)) - + # Model parameters as vector in order par <- c( as.numeric(params_modl$kphio), @@ -327,11 +327,11 @@ run_pmodel_f_bysite <- function( as.numeric(params_modl$rd_to_vcmax), as.numeric(params_modl$tau_acclim), as.numeric(params_modl$kc_jmax) - ) - + ) + ## C wrapper call out <- .Call( - + 'pmodel_f_C', ## Simulation parameters @@ -358,14 +358,14 @@ run_pmodel_f_bysite <- function( n = n, par = par, forcing = forcing - ) + ) # Prepare output to be a nice looking tidy data frame (tibble) ddf <- init_dates_dataframe( yrstart = firstyeartrend_forcing, yrend = firstyeartrend_forcing + nyeartrend_forcing - 1, noleap = TRUE) - + out <- out %>% as.matrix() %>% as.data.frame() %>% @@ -389,10 +389,10 @@ run_pmodel_f_bysite <- function( "wcont", "snow", "cond") - ) %>% + ) %>% as_tibble(.name_repair = "check_unique") %>% dplyr::bind_cols(ddf,.) - + } else { out <- tibble(date = as.Date("2000-01-01"), fapar = NA, @@ -415,9 +415,9 @@ run_pmodel_f_bysite <- function( snow = NA, cond = NA) } - + return(out) - + } .onUnload <- function(libpath) { diff --git a/R/runread_biomee_f.R b/R/runread_biomee_f.R index 6e182836..d118f49b 100644 --- a/R/runread_biomee_f.R +++ b/R/runread_biomee_f.R @@ -21,17 +21,20 @@ #' \donttest{ #' # Example BiomeE model run #' -#' mod_output <- runread_biomee_f( +#' runread_biomee_f( #' drivers = biomee_gs_leuning_drivers #' ) +#' runread_biomee_f( +#' drivers = biomee_p_model_drivers +#' ) #' } runread_biomee_f <- function( - drivers, - makecheck = TRUE, - parallel = FALSE, - ncores = 2 - ){ + drivers, + makecheck = TRUE, + parallel = FALSE, + ncores = 2 +){ # predefine variables for CRAN check compliance forcing <- init_cohort <- init_soil <- data <- @@ -44,32 +47,32 @@ runread_biomee_f <- function( multidplyr::cluster_assign(makecheck = FALSE) %>% multidplyr::cluster_library( packages = c("dplyr", "purrr", "rsofun") - ) + ) ## distribute to to cores, making sure all data from a specific site is sent to the same core df_out <- drivers %>% dplyr::group_by( id = row_number() ) %>% tidyr::nest('input' = c("sitename", - "params_siml", - "site_info", - "forcing", - "params_tile", - "params_species", - "init_cohort", - "init_soil")) %>% + "params_siml", + "site_info", + "forcing", + "params_tile", + "params_species", + "init_cohort", + "init_soil")) %>% multidplyr::partition(cl) %>% dplyr::mutate('data' = purrr::map(input, - ~run_biomee_f_bysite( - sitename = .x$sitename[[1]], - params_siml = .x$params_siml[[1]], - site_info = .x$site_info[[1]], - forcing = .x$forcing[[1]], - params_tile = .x$params_tile[[1]], - params_species = .x$params_species[[1]], - init_cohort = .x$init_cohort[[1]], - init_soil = .x$init_soil[[1]], - makecheck = makecheck ) - + ~run_biomee_f_bysite( + sitename = .x$sitename[[1]], + params_siml = .x$params_siml[[1]], + site_info = .x$site_info[[1]], + forcing = .x$forcing[[1]], + params_tile = .x$params_tile[[1]], + params_species = .x$params_species[[1]], + init_cohort = .x$init_cohort[[1]], + init_soil = .x$init_soil[[1]], + makecheck = makecheck ) + )) %>% dplyr::collect() %>% dplyr::ungroup() %>% @@ -87,7 +90,7 @@ runread_biomee_f <- function( "params_species", "init_cohort", "init_soil" - ) %>% + ) %>% dplyr::mutate(data = purrr::pmap( ., run_biomee_f_bysite, diff --git a/R/runread_pmodel_f.R b/R/runread_pmodel_f.R index 5120e706..f5a0f9a8 100644 --- a/R/runread_pmodel_f.R +++ b/R/runread_pmodel_f.R @@ -94,13 +94,16 @@ #' output <- rsofun::runread_pmodel_f( #' drivers = rsofun::p_model_drivers, #' par = params_modl) +#' output_vcmax25 <- rsofun::runread_pmodel_f( +#' drivers = rsofun::p_model_drivers_vcmax25, +#' par = params_modl) runread_pmodel_f <- function( - drivers, - par, - makecheck = TRUE, - parallel = FALSE, - ncores = 1){ + drivers, + par, + makecheck = TRUE, + parallel = FALSE, + ncores = 1){ # predefine variables for CRAN check compliance sitename <- params_siml <- site_info <- @@ -174,9 +177,9 @@ runread_pmodel_f <- function( df_out <- drivers %>% dplyr::mutate( data = purrr::pmap(., - run_pmodel_f_bysite, - params_modl = par, - makecheck = makecheck + run_pmodel_f_bysite, + params_modl = par, + makecheck = makecheck ) ) |> dplyr::select(sitename, site_info, data) diff --git a/R/zzz.R b/R/zzz.R index 8660c217..0b009768 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -6,9 +6,9 @@ is.nanull <- function(x) ifelse(any(is.null(x), is.na(x)), TRUE, FALSE) numeric_year <- function(x){ y <- as.numeric(format(x, format="%Y")) doy <- as.numeric(format(x, format="%j")) - 1 - - ifelse(y %% 4 == 0, - round(y + doy/366, 3), - round(y + doy/365, 3) - ) + + ifelse(y %% 4 == 0, + round(y + doy/366, 3), + round(y + doy/365, 3) + ) } diff --git a/README.md b/README.md index 8fe43c4a..2883535a 100644 --- a/README.md +++ b/README.md @@ -14,14 +14,12 @@ An R Simulating Optimal FUNctioning (RSOFUN) framework for site-scale simulation ### Stable release -**rsofun is not currently available on CRAN.** - - +``` ### Development release diff --git a/analysis/01-sensitivity-analysis.R b/analysis/01-sensitivity-analysis.R index 822e30f9..ebefbdbb 100644 --- a/analysis/01-sensitivity-analysis.R +++ b/analysis/01-sensitivity-analysis.R @@ -13,8 +13,7 @@ set.seed(432) # Define log-likelihood function ll_pmodel <- function( - par_v # a vector of all calibratable parameters - # including errors + par_v # a vector of all calibratable parameters including errors ){ rsofun::cost_likelihood_pmodel( # likelihood cost function from package par_v, @@ -118,5 +117,5 @@ gg <- morrisOut.df |> ) + coord_flip() # make horizontal -# ggsave("./analysis/paper_results_files/morris.pdf", plot = gg, width = 5, height = 3) -# ggsave("./analysis/paper_results_files/morris.png", plot = gg, width = 5, height = 3) +ggsave("./analysis/paper_results_files/morris.pdf", plot = gg, width = 5, height = 3) +ggsave("./analysis/paper_results_files/morris.png", plot = gg, width = 5, height = 3) diff --git a/analysis/02-bayesian-calibration.R b/analysis/02-bayesian-calibration.R index f6851a68..eec4c298 100644 --- a/analysis/02-bayesian-calibration.R +++ b/analysis/02-bayesian-calibration.R @@ -44,7 +44,6 @@ plot_prior_posterior_density <- function( priorMat <- getSetup(x)$prior$sampler(10000) # nPriorDraws = 10000 # Parameter names - nPar <- ncol(posteriorMat) parNames <- colnames(posteriorMat) # rename columns priorMat colnames(priorMat) <- parNames @@ -81,15 +80,15 @@ settings_calib <- list( control = list( sampler = "DEzs", settings = list( - burnin = 1500, + burnin = 3000, iterations = 12000, nrChains = 3, # number of independent chains startValue = 3 # number of internal chains to be sampled )), par = list( - kphio = list(lower = 0.03, upper = 0.15, init = 0.05), - soilm_betao = list(lower = 0, upper = 1, init = 0.2), - kc_jmax = list(lower = 0.2, upper = 0.8, init = 0.41), + kphio = list(lower = 0.02, upper = 0.15, init = 0.05), + kphio_par_b = list(lower = 10, upper = 30, init = 20), + kc_jmax = list(lower = 0.1, upper = 0.8, init = 0.40), err_gpp = list(lower = 0.1, upper = 3, init = 0.8) ) ) @@ -102,9 +101,9 @@ par_calib <- calib_sofun( obs = rsofun::p_model_validation, settings = settings_calib, par_fixed = list( - kphio_par_a = -0.0025, # define model parameter values from - kphio_par_b = 20, # Stocker et al. 2020 + kphio_par_a = -0.0025, # define model parameter values from Stocker et al. 2020 soilm_thetastar = 0.6*240, + soilm_betao = 0.2, beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, tau_acclim = 30.0), @@ -116,11 +115,18 @@ toc() # Stop measuring time # Plot prior and posterior distributions gg <- plot_prior_posterior_density(par_calib$mod) +get_runtime <- function(par_calib) {# function(settings_calib){ + total_time_secs <- sum(unlist(lapply( + par_calib$mod, + function(curr_chain){curr_chain$settings$runtime[["elapsed"]]}))) + return(sprintf("Total runtime: %.0f secs", total_time_secs)) +} + # Plot MCMC diagnostics -# plot(par_calib$mod) -# summary(par_calib$mod) # Gives Gelman Rubin multivariate of 1.019 -# summary(par_calib$par) -# print(get_runtime(par_calib)) +plot(par_calib$mod) +summary(par_calib$mod) # Gives Gelman Rubin multivariate of 1.019 +summary(par_calib$par) +print(get_runtime(par_calib)) # Output ## Define MCMC postprocessing @@ -128,17 +134,17 @@ get_settings_str <- function(par_calib) {# function(settings_calib){ stopifnot(is(par_calib$mod, "mcmcSamplerList")) # explore what's in a mcmcSamplerList: - # summary(par_calib$mod) - # plot(par_calib$mod) + # summary(par_calib$mod) + # plot(par_calib$mod) individual_chains <- par_calib$mod nrChains <- length(individual_chains) # number of chains - # plot(individual_chains[[1]]) # chain 1 - # plot(individual_chains[[2]]) # chain 2 - # plot(individual_chains[[3]]) # chain 3 - # class(individual_chains[[1]]$setup); individual_chains[[1]]$setup # Bayesian Setup - # individual_chains[[1]]$chain - # individual_chains[[1]]$X - # individual_chains[[1]]$Z + # plot(individual_chains[[1]]) # chain 1 + # plot(individual_chains[[2]]) # chain 2 + # plot(individual_chains[[3]]) # chain 3 + # class(individual_chains[[1]]$setup); individual_chains[[1]]$setup # Bayesian Setup + # individual_chains[[1]]$chain + # individual_chains[[1]]$X + # individual_chains[[1]]$Z nrInternalChains <- lapply(individual_chains, function(curr_chain){curr_chain$settings$nrChains}) |> unlist() |> unique() |> paste0(collapse = "-") nrIterations <- lapply(individual_chains, function(curr_chain){curr_chain$settings$iterations})|> unlist() |> unique() |> paste0(collapse = "-") nrBurnin <- lapply(individual_chains, function(curr_chain){curr_chain$settings$burnin}) |> unlist() |> unique() |> paste0(collapse = "-") @@ -149,17 +155,10 @@ get_settings_str <- function(par_calib) {# function(settings_calib){ "Sampler-%s-%siterations_ofwhich%sburnin_chains(%sx%s)", sampler_name, nrIterations, nrBurnin, nrChains, nrInternalChains)) } -get_runtime <- function(par_calib) {# function(settings_calib){ - total_time_secs <- sum(unlist(lapply( - par_calib$mod, - function(curr_chain){curr_chain$settings$runtime[["elapsed"]]}))) - return(sprintf("Total runtime: %.0f secs", total_time_secs)) -} -dir.create("./analysis/paper_results_files2") settings_string <- get_settings_str(par_calib) -saveRDS(par_calib, file = paste0("./analysis/paper_results_files2/",settings_string,"_prior_posterior.RDS")) -ggsave(paste0("./analysis/paper_results_files2/",settings_string,"_prior_posterior.pdf"), plot = gg, width = 6, height = 5) -ggsave(paste0("./analysis/paper_results_files2/",settings_string,"_prior_posterior.png"), plot = gg, width = 6, height = 5) +saveRDS(par_calib, file = paste0("./analysis/paper_results_files/",settings_string,"_prior_posterior.RDS")) +ggsave(paste0("./analysis/paper_results_files/",settings_string,"_prior_posterior.pdf"), plot = gg, width = 6, height = 5) +ggsave(paste0("./analysis/paper_results_files/",settings_string,"_prior_posterior.png"), plot = gg, width = 6, height = 5) diff --git a/analysis/03-uncertainty-estimation.R b/analysis/03-uncertainty-estimation.R index 2087e48a..f346799e 100644 --- a/analysis/03-uncertainty-estimation.R +++ b/analysis/03-uncertainty-estimation.R @@ -39,9 +39,9 @@ run_pmodel <- function(sample_par){ par = list( # copied from par_fixed above kphio = sample_par$kphio, kphio_par_a = -0.0025, - kphio_par_b = 20, + kphio_par_b = sample_par$kphio_par_b, soilm_thetastar = 0.6*240, - soilm_betao = sample_par$soilm_betao, + soilm_betao = 0.2, beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, tau_acclim = 30.0, @@ -87,21 +87,30 @@ plot_gpp_error <- ggplot(data = data_to_plot) + aes(ymin = gpp_pred_q05, ymax = gpp_pred_q95, x = date, - fill = "Model uncertainty")) + + fill = "Model uncertainty" + ), + alpha = 0.5) + geom_ribbon( aes(ymin = gpp_q05, ymax = gpp_q95, x = date, - fill = "Parameter uncertainty")) + - geom_line( - aes(x = date, - y = gpp_q50, - color = "Predictions")) + + fill = "Parameter uncertainty" + ), + alpha = 0.5) + # Include observations in the plot - geom_point(shape = 3, + geom_point( aes(x = date, y = gpp_obs, - color = "Observations")) + + color = "Observations" + ), + ) + + geom_line( + aes(x = date, + y = gpp_q50, + color = "Predictions" + ), + alpha = 0.6 + ) + theme_classic() + theme(panel.grid.major.y = element_line(), legend.position = "bottom") + @@ -110,27 +119,21 @@ plot_gpp_error <- ggplot(data = data_to_plot) + y = expression(paste("GPP (g C m"^-2, "s"^-1, ")")) ) -plot_gpp_error <- plot_gpp_error + +plot_gpp_error <- plot_gpp_error + scale_color_manual(name = "", breaks = c("Observations", - "Predictions", - "Non-calibrated predictions"), + "Predictions"), values = c(t_col("black", 10), - t_col("#E69F00", 10), - t_col("#56B4E9", 10))) + - scale_fill_manual(name = "", - breaks = c("Model uncertainty", - "Parameter uncertainty"), - values = c(t_col("#E69F00", 60), - t_col("#009E73", 40))) + t_col("grey40", 10)))# + +# scale_fill_manual(name = "", +# breaks = c("Model uncertainty", +# "Parameter uncertainty"), +# values = c(t_col("#E69F00", 60), +# t_col("#009E73", 40))) plot_gpp_error -dir.create("./analysis/paper_results_files2") settings_string <- get_settings_str(par_calib) -ggsave(paste0("./analysis/paper_results_files2/",settings_string,"_gpp_predictions_observations.pdf"), plot = plot_gpp_error, width = 6, height = 5) -ggsave(paste0("./analysis/paper_results_files2/",settings_string,"_gpp_predictions_observations.png"), plot = plot_gpp_error, width = 6, height = 5) - - - +ggsave(paste0("./analysis/paper_results_files/",settings_string,"_gpp_predictions_observations.pdf"), plot = plot_gpp_error, width = 6, height = 5) +ggsave(paste0("./analysis/paper_results_files/",settings_string,"_gpp_predictions_observations.png"), plot = plot_gpp_error, width = 6, height = 5) plot_gpp_error2 <- ggplot(data = data_to_plot) + @@ -175,7 +178,6 @@ plot_gpp_error2 <- plot_gpp_error2 + values = c(t_col("#E69F00", 60), t_col("#009E73", 40))) plot_gpp_error2 -dir.create("./analysis/paper_results_files2") -ggsave(paste0("./analysis/paper_results_files2/",settings_string,"_gpp_predictions_observations_2.pdf"), plot = plot_gpp_error2, width = 6, height = 5) -ggsave(paste0("./analysis/paper_results_files2/",settings_string,"_gpp_predictions_observations_2.png"), plot = plot_gpp_error2, width = 6, height = 5) +ggsave(paste0("./analysis/paper_results_files/",settings_string,"_gpp_predictions_observations_2.pdf"), plot = plot_gpp_error2, width = 6, height = 5) +ggsave(paste0("./analysis/paper_results_files/",settings_string,"_gpp_predictions_observations_2.png"), plot = plot_gpp_error2, width = 6, height = 5) diff --git a/cran-comments.md b/cran-comments.md index 4883d78f..6e3635a9 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,43 +1,3 @@ -Dear CRAN team, +## R CMD check -This is the re-submission of the {rsofun} package v4.4.1. - -In order to correct the outcome of the automatic CRAN checks, we have removed the --fc-prototypes-external flag for flang compilers. Our local and GitHub checks -pass, and we hope that the previous changes on the Fortran code alone correct -the failure of our previous release. - -Thank you for your time and support. - ---- - -The {rsofun} package provides the implementation of a modelling framework for site-scale simulations of ecosystem processes, with low level routines in Fortran 90. It contains the following models: -- P-model for leaf-level acclimation of photosynthesis from Stocker et al. (2019). -- SPLASH for bioclimatic variables, including the surface radiation budget and the soil water balance from Davis et al. (2017). -- BiomeE for comprehensive simulations of ecosystem carbon and water cycling, tree growth, and tree cohort-explicit forest dynamics following the Perfect Plasticity Approximation, from Weng et al. (2015). - -This package is an extension of {rpmodel} in the sense that it expands the P-model implementation and provides functions for multiple-site simulations and model parameter calibration. - -The full documentation can be found at the github repository link: https://geco-bern.github.io/rsofun - -Code coverage sits at ~76%, with remaining uncovered code pertaining to minor input data format checks of the main functions. The underlying P-model implementation is based on the {rpmodel} and the parameter calibration routines use packages {GenSA} and {BayesianTools}. - -I hope this package is useful for other earth system scientists and the larger CRAN community. Kind regards, Josefa Arán. - ---- - -I have read and agree to the CRAN policies enumerated here: https://cran.r-project.org/web/packages/policies.html - -## Local, github actions and r-hub checks - -- Pop!_OS 22.04 install on R 4.3 - -- Ubuntu 22.04, MacOS and Windows on github actions (devel / release) - -- rhub::check_on_cran() with only notes for latex elements - -- codecove.io code coverage at ~76% - -## Github actions R CMD check results - -0 errors | 0 warnings | 0 notes +Status: Ok diff --git a/data-raw/generate_biomee_drivers.R b/data-raw/generate_biomee_drivers.R index 2582d8bb..9060e406 100644 --- a/data-raw/generate_biomee_drivers.R +++ b/data-raw/generate_biomee_drivers.R @@ -19,16 +19,12 @@ sitename <- "CH-Lae" # Take only year 2004 to 2014, corresponding to subset of data for site CH-Lae siteinfo <- tibble( - sitename = "CH-Lae", lon = 8.365, lat = 47.47808, elv = 700, year_start = 2004, year_end = 2014, - classid = NA, c4 = FALSE, - whc = NA, - koeppen_code = NA, igbp_land_use = "Mixed Forests", plant_functional_type = "Broadleaf trees") @@ -166,14 +162,14 @@ rh_to_vpd <- function(temp, # Air temperature (deg C) rh # Relative humidity (< 1) ) { esat <- 611.0 * exp( (17.27 * temp)/(temp + 237.3) ) - + return(esat * (1.0 - rh)) # VPD (Pa) } rad_to_ppfd <- function(rad # Downwelling radiation (W m-2) ) { kcFCE <- 2.04 # from flux to energy conversion, umol/J (Meek et al., 1984) - + return(rad * kcFCE * 1.0e-6) # PPFD (mol m-2 s-1) } @@ -197,11 +193,11 @@ build_forcing <- function(forcing_data, hourly) { vpd = rh_to_vpd(temp, rh / 100.0), ppfd = rad_to_ppfd(rad)) %>% select(date,hod,temp,rain,vpd,ppfd,patm,wind,co2,) - return(forcing) + return(forcing) } build_driver <- function(params_siml, forcing) { - drivers <- tibble( + drivers <- tibble( sitename, site_info = list(tibble(siteinfo)), params_siml = list(tibble(params_siml)), diff --git a/data-raw/generate_biomee_output.R b/data-raw/generate_biomee_output.R index f53556f4..10dc434c 100644 --- a/data-raw/generate_biomee_output.R +++ b/data-raw/generate_biomee_output.R @@ -7,7 +7,7 @@ generate_output <- function(drivers) { drivers, makecheck = TRUE, parallel = FALSE) - + return(out) } diff --git a/data-raw/generate_pmodel_output.R b/data-raw/generate_pmodel_output.R new file mode 100644 index 00000000..9ec53d53 --- /dev/null +++ b/data-raw/generate_pmodel_output.R @@ -0,0 +1,33 @@ +#!/usr/bin/env Rscript + +library(rsofun) + +# Define model parameter values from previous work +params_modl <- list( + kphio = 0.04998, # setup ORG in Stocker et al. 2020 GMD + kphio_par_a = 0.0, # disable temperature-dependence of kphio + kphio_par_b = 1.0, + soilm_thetastar = 0.6 * 240, # old setup with soil moisture stress + soilm_betao = 0.0, + beta_unitcostratio = 146.0, + rd_to_vcmax = 0.014, # from Atkin et al. 2015 for C3 herbaceous + tau_acclim = 30.0, + kc_jmax = 0.41 +) + +# Run the model for these parameters and the example drivers +p_model_output <- rsofun::runread_pmodel_f( + drivers = rsofun::p_model_drivers, + par = params_modl) + +p_model_output_vcmax25 <- rsofun::runread_pmodel_f( + drivers = rsofun::p_model_drivers_vcmax25, + par = params_modl) + +save(p_model_output, + file ="data/p_model_output.rda", + compress = "xz") + +save(p_model_output_vcmax25, + file ="data/p_model_output_vcmax25.rda", + compress = "xz") diff --git a/data/biomee_gs_leuning_drivers.rda b/data/biomee_gs_leuning_drivers.rda index 9ccaa532..5df8b5be 100644 Binary files a/data/biomee_gs_leuning_drivers.rda and b/data/biomee_gs_leuning_drivers.rda differ diff --git a/data/biomee_p_model_drivers.rda b/data/biomee_p_model_drivers.rda index 61dbaf62..18d2f22e 100644 Binary files a/data/biomee_p_model_drivers.rda and b/data/biomee_p_model_drivers.rda differ diff --git a/data/p_model_output.rda b/data/p_model_output.rda new file mode 100644 index 00000000..20230dd2 Binary files /dev/null and b/data/p_model_output.rda differ diff --git a/data/p_model_output_vcmax25.rda b/data/p_model_output_vcmax25.rda new file mode 100644 index 00000000..62f294d7 Binary files /dev/null and b/data/p_model_output_vcmax25.rda differ diff --git a/man/biomee_gs_leuning_drivers.Rd b/man/biomee_gs_leuning_drivers.Rd index 65c1135a..c43ee693 100644 --- a/man/biomee_gs_leuning_drivers.Rd +++ b/man/biomee_gs_leuning_drivers.Rd @@ -3,23 +3,155 @@ \docType{data} \name{biomee_gs_leuning_drivers} \alias{biomee_gs_leuning_drivers} -\title{rsofun BiomeE driver data} +\title{rsofun BiomeE driver data (Leuning photosynthesis model)} \format{ -A tibble of driver data: +A tibble of driver data. \describe{ - \item{sitename}{site name} - \item{params_siml}{model parameters} - \item{site_info}{site information} - \item{soil_texture}{soil texture data} - \item{forcing}{forcing data} + \item{sitename}{Site name} + \item{params_siml}{Simulation parameters as a data.frame, including + the following data: + \describe{ + \item{spinup}{Flag indicating whether this simulation does spin-up.} + \item{spinupyears}{Number of spin-up years.} + \item{recycle}{Length of standard recycling period (years).} + \item{firstyeartrend}{First transient year.} + \item{nyeartrend}{Number of transient years.} + \item{steps_per_day}{Time resolution (day-1).} + \item{do_U_shaped_mortality}{Flag indicating whether U-shaped + mortality is used.} + \item{update_annualLAImax}{Flag indicating whether updating + LAImax according to mineral N in soil.} + \item{do_closedN_run}{Flag indicating whether doing N closed + runs to recover N balance enforcing 0.2 kg N m-2 in the inorganic N pool.} + \item{code_method_photosynth}{String specifying the method of photosynthesis + used in the model, either "pmodel" or "gs_leuning".document()} + \item{code_method_mortality}{String indicating the type of mortality in the + model. One of the following: "dbh" is size-dependent mortality, "const_selfthin" + is constant self thinning (in development), "cstarvation" is carbon starvation, and + "growthrate" is growth rate dependent mortality.} + }} + \item{site_info}{Site meta info in a data.frame. +This data structure can be freely used for documenting the dataset, but must include at least the following data: + \describe{ + \item{lon}{Longitude of the site location.} + \item{lat}{Latitude of the site location.} + \item{elv}{Elevation of the site location, in meters.} + }} + \item{forcing}{Forcing data.frame used as input + \describe{ + \item{ppfd}{Photosynthetic photon flux density (mol s-1 m-2)} + \item{tair}{Air temperature (deg C)} + \item{vpd}{Vapor pressure deficit (Pa)} + \item{rain}{Precipitation (kgH2O m-2 s-1 == mm s-1)} + \item{wind}{Wind velocity (m s-1)} + \item{pair}{Atmospheric pressure (pa)} + \item{co2}{CO2 atmospheric concentration (ppm)} + }} + \item{params_tile}{Tile-level model parameters, into a single row data.frame, including + the following data: + \describe{ + \item{soiltype}{Integer indicating the type of soil: Sand = 1, LoamySand = 2, + SandyLoam = 3, SiltLoam = 4, FrittedClay = 5, Loam = 6, Clay = 7.} + \item{FLDCAP}{Field capacity (vol/vol). Water remaining in a soil after it + has been thoroughly saturated and allowed to drain freely.} + \item{WILTPT}{Wilting point (vol/vol). Water content of a soil at which + plants wilt and fail to recover.} + \item{K1}{Fast soil C decomposition rate (year\eqn{^{-1}}).} + \item{K2}{Slow soil C decomposition rate (year\eqn{^{-1}}).} + \item{K_nitrogen}{Mineral Nitrogen turnover rate (year\eqn{^{-1}}).} + \item{MLmixRatio}{Ratio of C and N returned to litters from microbes.} + \item{etaN}{N loss rate through runoff (organic and mineral) (year\eqn{^{-1}}).} + \item{LMAmin}{Minimum LMA, leaf mass per unit area, kg C m\eqn{^{-2}}.} + \item{fsc_fine}{Fraction of fast turnover carbon in fine biomass.} + \item{fsc_wood}{Fraction of fast turnover carbon in wood biomass.} + \item{GR_factor}{Growth respiration factor.} + \item{l_fract}{Fraction of the carbon retained after leaf drop.} + \item{retransN}{Retranslocation coefficient of nitrogen.} + \item{f_initialBSW}{Coefficient for setting up initial sapwood.} + \item{f_N_add}{Re-fill of N for sapwood.} + \item{tf_base}{Calibratable scalar for respiration, used to increase LUE + levels.} + \item{par_mort}{Canopy mortality parameter.} + \item{par_mort_under}{Parameter for understory mortality.} + }} + \item{params_species}{A data.frame containing species-specific model parameters, + with one species per row, including the following data: + \describe{ + \item{lifeform}{Integer set to 0 for grasses and 1 for trees.} + \item{phenotype}{Integer set to 0 for deciduous and 1 for evergreen.} + \item{pt}{Integer indicating the type of plant according to photosynthesis: + 0 for C3; 1 for C4} + \item{alpha_FR}{Fine root turnonver rate (year\eqn{^{-1}}).} + \item{rho_FR}{Material density of fine roots (kg C m\eqn{^{-3}}).} + \item{root_r}{Radius of the fine roots, in m.} + \item{root_zeta}{e-folding parameter of root vertical distribution, in m.} + \item{Kw_root}{Fine root water conductivity (mol m\eqn{^{-2}} + s\eqn{^{-1}} MPa\eqn{^{-1}}).} + \item{leaf_size}{Characteristic leaf size.} + \item{Vmax}{Max RuBisCo rate, in mol m\eqn{^{-2}} s\eqn{^{-1}}.} + \item{Vannual}{Annual productivity per unit area at full sun (kg C + m\eqn{^{-2}} year\eqn{^{-2}}).} + \item{wet_leaf_dreg}{Wet leaf photosynthesis down-regulation.} + \item{m_cond}{Factor of stomatal conductance.} + \item{alpha_phot}{Photosynthesis efficiency.} + \item{gamma_L}{Leaf respiration coefficient, in year\eqn{^{-1}}.} + \item{gamma_LN}{Leaf respiration coefficient per unit N.} + \item{gamma_SW}{Sapwood respiration rate, in kg C m\eqn{^{-2}} year\eqn{^{-1}}.} + \item{gamma_FR}{Fine root respiration rate, kg C kg C\eqn{^{-1}} + year\eqn{^{-1}}.} + \item{tc_crit}{Critical temperature triggerng offset of phenology, in Kelvin.} + \item{tc_crit_on}{Critical temperature triggerng onset of phenology, in Kelvin.} + \item{gdd_crit}{Critical value of GDD5 for turning ON growth season.} + \item{betaON}{Critical soil moisture for phenology onset.} + \item{betaOFF}{Critical soil moisture for phenology offset.} + \item{seedlingsize}{Initial size of seedlings, in kg C per individual.} + \item{LNbase}{Basal leaf N per unit area, in kg N m\eqn{^{-2}}.} + \item{lAImax}{Maximum crown LAI (leaf area index).} + \item{Nfixrate0}{Reference N fixation rate (kg N kg C\eqn{^{-1}} root).} + \item{NfixCost0}{Carbon cost of N fixation (kg C kg N\eqn{^{-1}}).} + \item{phiCSA}{Ratio of sapwood area to leaf area.} + \item{mortrate_d_c}{Canopy tree mortality rate (year\eqn{^{-1}}).} + \item{mortrate_d_u}{Understory tree mortality rate (year\eqn{^{-1}}).} + \item{maturalage}{Age at which trees can reproduce (years).} + \item{v_seed}{Fraction of G_SF to G_F.} + \item{fNSmax}{Multiplier for NSNmax as sum of potential bl and br.} + \item{LMA}{Leaf mass per unit area (kg C m\eqn{^{-2}}).} + \item{rho_wood}{Wood density (kg C m\eqn{^{-3}}).} + \item{alphaBM}{Coefficient for allometry (biomass = alphaBM * DBH ** thetaBM).} + \item{thetaBM}{Coefficient for allometry (biomass = alphaBM * DBH ** thetaBM).} + \item{kphio}{Quantum yield efficiency \eqn{\varphi_0}, + in mol mol\eqn{^{-1}}.} + \item{phiRL}{Ratio of fine root to leaf area.} + \item{LAI_light}{Maximum LAI limited by light.} + }} + \item{init_cohort}{A data.frame of initial cohort specifications, including + the following data: + \describe{ + \item{init_cohort_species}{Index of a species described in param_species.} + \item{init_cohort_nindivs}{Initial individual density, in individuals per + m\eqn{^{2}}.} + \item{init_cohort_bsw}{Initial biomass of sapwood, in kg C per individual.} + \item{init_cohort_bHW}{Initial biomass of heartwood, in kg C per tree.} + \item{init_cohort_nsc}{Initial non-structural biomass.} + }} + \item{init_soil}{A data.frame of initial soil pools, including + the following data: + \describe{ + \item{init_fast_soil_C}{Initial fast soil carbon, in kg C m\eqn{^{-2}}.} + \item{init_slow_soil_C}{Initial slow soil carbon, in kg C m\eqn{^{-2}}.} + \item{init_Nmineral}{Mineral nitrogen pool, in kg N m\eqn{^{-2}}.} + \item{N_input}{Annual nitrogen input to soil N pool, in kg N m\eqn{^{-2}} + year\eqn{^{-1}}.} + }} } } \usage{ biomee_gs_leuning_drivers } \description{ -Small tests dataset to validate if compiled code -and optimization routines can run using the -Leuning specifications +Small dataset representing the driver to run the BiomeE-model at the CH-LAE site +using the Leuning photosynthesis specification (and half-hourly time step) +It can also be used together with leaf trait data from CH-LAE (\code{\link{biomee_validation}}) +to optimize model parameters. } \keyword{datasets} diff --git a/man/biomee_gs_leuning_output.Rd b/man/biomee_gs_leuning_output.Rd index af550ad0..8bd4ff89 100644 --- a/man/biomee_gs_leuning_output.Rd +++ b/man/biomee_gs_leuning_output.Rd @@ -3,7 +3,7 @@ \docType{data} \name{biomee_gs_leuning_output} \alias{biomee_gs_leuning_output} -\title{rsofun BiomeE (GS leuning) output data} +\title{rsofun BiomeE (gs_leuning) output data} \format{ An object of class \code{tbl_df} (inherits from \code{tbl}, \code{data.frame}) with 1 rows and 2 columns. } @@ -11,6 +11,6 @@ An object of class \code{tbl_df} (inherits from \code{tbl}, \code{data.frame}) w biomee_gs_leuning_output } \description{ -Example output dataset using BiomeE (GS leuning) +Example output dataset from a BiomeE-model run (gs_leuning) } \keyword{datasets} diff --git a/man/biomee_p_model_drivers.Rd b/man/biomee_p_model_drivers.Rd index 0044ad1a..aa144df1 100644 --- a/man/biomee_p_model_drivers.Rd +++ b/man/biomee_p_model_drivers.Rd @@ -3,23 +3,17 @@ \docType{data} \name{biomee_p_model_drivers} \alias{biomee_p_model_drivers} -\title{rsofun BiomeE driver data} +\title{rsofun BiomeE driver data (P-model photosynthesis model)} \format{ -A tibble of driver data: -\describe{ - \item{sitename}{site name} - \item{params_siml}{model parameters} - \item{site_info}{site information} - \item{soil_texture}{soil texture data} - \item{forcing}{forcing data} -} +See \code{\link{biomee_gs_leuning_drivers}} } \usage{ biomee_p_model_drivers } \description{ -Small tests dataset to validate if compiled code -and optimization routines can run using the -p-model specifications +Small dataset representing the driver to run the BiomeE-model at the CH-LAE site +using the P-model photosynthesis specification (and daily time step). +It can also be used together with leaf trait data from CH-LAE (\code{\link{biomee_validation}}) +to optimize model parameters. } \keyword{datasets} diff --git a/man/biomee_p_model_output.Rd b/man/biomee_p_model_output.Rd index 3f1696fa..6a4ef5b2 100644 --- a/man/biomee_p_model_output.Rd +++ b/man/biomee_p_model_output.Rd @@ -3,7 +3,7 @@ \docType{data} \name{biomee_p_model_output} \alias{biomee_p_model_output} -\title{rsofun BiomeE (p-model) output data} +\title{rsofun BiomeE (P-model) output data} \format{ An object of class \code{tbl_df} (inherits from \code{tbl}, \code{data.frame}) with 1 rows and 2 columns. } @@ -11,6 +11,6 @@ An object of class \code{tbl_df} (inherits from \code{tbl}, \code{data.frame}) w biomee_p_model_output } \description{ -Example output dataset using BiomeE (p-model) +Example output dataset from a BiomeE-model run (p-model) } \keyword{datasets} diff --git a/man/biomee_validation.Rd b/man/biomee_validation.Rd index 94ad3d8d..0ef500b6 100644 --- a/man/biomee_validation.Rd +++ b/man/biomee_validation.Rd @@ -20,7 +20,7 @@ Dataset. https://doi.org/10.18140/FLX/1440134 biomee_validation } \description{ -Small tests dataset to validate -calibration routines +Small example dataset of target observations (leaf trait data) at the CH-LAE site +to optimize model parameters with the function \code{\link{calib_sofun}} } \keyword{datasets} diff --git a/man/init_dates_dataframe.Rd b/man/init_dates_dataframe.Rd index 837c81a7..881b9112 100644 --- a/man/init_dates_dataframe.Rd +++ b/man/init_dates_dataframe.Rd @@ -48,9 +48,7 @@ A tibble with dates. \description{ Creates a tibble with rows for each date from \code{'yrstart'} to \code{'yrend'} in \code{'yyyy-mm-dd'} format. Intervals of dates are specified by argument -\code{'freq'}. -} -\examples{ - ddf <- init_dates_dataframe( 2000, 2003, startmoy=1, startdoy=1, - freq="days", endmoy=12, enddom=31, noleap=FALSE ) +\code{'freq'}. + ddf <- init_dates_dataframe(2000, 2003, startmoy=1, startdoy=1, + freq="days", endmoy=12, enddom=31, noleap=FALSE) } diff --git a/man/p_model_drivers.Rd b/man/p_model_drivers.Rd index 60f695e9..abda6ae0 100644 --- a/man/p_model_drivers.Rd +++ b/man/p_model_drivers.Rd @@ -8,16 +8,16 @@ A tibble of driver data: \describe{ \item{sitename}{A character string containing the site name.} - \item{forcing}{A tibble of a time series of forcing climate data, including - the following variables: + \item{forcing}{A tibble of a time series of forcing climate data, including + the following data: \describe{ \item{date}{Date of the observation in YYYY-MM-DD format.} \item{temp}{Daytime average air temperature in \eqn{^\circ}C.} \item{vpd}{Daytime average vapour pressure deficit in Pa.} - \item{ppfd}{Photosynthetic photon flux density (PPFD) in + \item{ppfd}{Photosynthetic photon flux density (PPFD) in mol m\eqn{^{-2}} s\eqn{^{-1}}. If all values are NA, it indicates that PPFD should be calculated by the SPLASH model.} - \item{netrad}{Net radiation in W m\eqn{^{-2}}. This is currently + \item{netrad}{Net radiation in W m\eqn{^{-2}}. This is currently ignored as a model forcing.} \item{patm}{Atmospheric pressure in Pa.} \item{snow}{Snow in water equivalents mm s\eqn{^{-1}}.} @@ -31,7 +31,8 @@ A tibble of driver data: net radiation are not prescribed.} } } - \item{params_siml}{A tibble of simulation parameters. + \item{params_siml}{A tibble of simulation parameters, including + the following data: \describe{ \item{spinup}{A logical value indicating whether this simulation does spin-up.} \item{spinupyears}{Number of spin-up years.} @@ -48,6 +49,7 @@ A tibble of driver data: } } \item{site_info}{A tibble containing site meta information. +This data structure can be freely used for documenting the dataset, but must include at least the following data: \describe{ \item{lon}{Longitude of the site location in degrees east.} \item{lat}{Latitude of the site location in degrees north.} @@ -80,7 +82,9 @@ surfaces for global land areas. International Journal of Climatology 37 (12): 43 p_model_drivers } \description{ -Small tests dataset to validate if compiled code -and optimization routines can run +Small dataset representing the driver to run the P-model at the FR-Pue site. +It can also be used together with daily GPP flux time series data from CH-LAE +(\code{\link{p_model_validation}}) to optimize model parameters. +To optimize model parameters to leaf traits data use the datasets \code{\link{p_model_drivers_vcmax25}} and \code{\link{p_model_validation_vcmax25}}. } \keyword{datasets} diff --git a/man/p_model_drivers_vcmax25.Rd b/man/p_model_drivers_vcmax25.Rd index eb5278f8..6b857d40 100644 --- a/man/p_model_drivers_vcmax25.Rd +++ b/man/p_model_drivers_vcmax25.Rd @@ -5,57 +5,7 @@ \alias{p_model_drivers_vcmax25} \title{rsofun P-model driver data (for leaf traits)} \format{ -A tibble of driver data: -\describe{ - \item{sitename}{A character string containing the site name.} - \item{forcing}{A tibble of a time series of forcing climate data, including - the following variables: - \describe{ - \item{date}{Date of the observation in YYYY-MM-DD format.} - \item{temp}{Daytime average air temperature in \eqn{^\circ}C.} - \item{vpd}{Daytime average vapour pressure deficit in Pa.} - \item{ppfd}{Photosynthetic photon flux density (PPFD) in - mol m\eqn{^{-2}} s\eqn{^{-1}}. If all values are NA, it indicates that - PPFD should be calculated by the SPLASH model.} - \item{netrad}{Net radiation in W m\eqn{^{-2}}. This is currently - ignored as a model forcing.} - \item{patm}{Atmospheric pressure in Pa.} - \item{snow}{Snow in water equivalents mm s\eqn{^{-1}}.} - \item{rain}{Rain as precipitation in liquid form in mm s\eqn{^{-1}}.} - \item{tmin}{Daily minimum air temperature in \eqn{^\circ}C.} - \item{tmax}{Daily maximum air temperature in \eqn{^\circ}C.} - \item{fapar}{Fraction of photosynthetic active radiation (fAPAR), taking - values between 0 and 1.} - \item{co2}{Atmospheric CO\eqn{_2} concentration.} - \item{ccov}{Cloud coverage in \%. This is only used when either PPFD or - net radiation are not prescribed.} - } - } - \item{params_siml}{A tibble of simulation parameters. - \describe{ - \item{spinup}{A logical value indicating whether this simulation does spin-up.} - \item{spinupyears}{Number of spin-up years.} - \item{recycle}{Length of standard recycling period, in years.} - \item{outdt}{An integer indicating the output periodicity.} - \item{ltre}{A logical value, \code{TRUE} if evergreen tree.} - \item{ltne}{A logical value, \code{TRUE} if evergreen tree and N-fixing.} - \item{ltrd}{A logical value, \code{TRUE} if deciduous tree.} - \item{ltnd}{A logical value, \code{TRUE} if deciduous tree and N-fixing.} - \item{lgr3}{A logical value, \code{TRUE} if grass with C3 photosynthetic pathway.} - \item{lgn3}{A logical value, \code{TRUE} if grass with C3 photosynthetic - pathway and N-fixing.} - \item{lgr4}{A logical value, \code{TRUE} if grass with C4 photosynthetic pathway.} - } - } - \item{site_info}{A tibble containing site meta information. - \describe{ - \item{lon}{Longitude of the site location in degrees east.} - \item{lat}{Latitude of the site location in degrees north.} - \item{elv}{Elevation of the site location, in meters above sea level.} - \item{whc}{A numeric value for the rooting zone water holding capacity (in mm)} - } - } -} +See \code{\link{p_model_drivers}} } \source{ Atkin, O. K., Bloomfield, K. J., Reich, P. B., Tjoelker, M. G., Asner, G. P., Bonal, D., et al. (2015). @@ -84,7 +34,9 @@ Atmospheric carbon dioxide variations at Mauna Loa Observatory, Hawaii, Tellus, p_model_drivers_vcmax25 } \description{ -Small tests dataset to validate if compiled code -and optimization routines can run for leaf traits data +Small dataset representing the driver to run the P-model at four separate sites. +It can also be used together with leaf traits data from these four sites +(\code{\link{p_model_validation_vcmax25}}) to optimize model parameters. +To optimize model parameters to GPP flux data use the datasets \code{\link{p_model_drivers}} and \code{\link{p_model_validation}}. } \keyword{datasets} diff --git a/man/p_model_output.Rd b/man/p_model_output.Rd new file mode 100644 index 00000000..050d52c1 --- /dev/null +++ b/man/p_model_output.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data.R +\docType{data} +\name{p_model_output} +\alias{p_model_output} +\title{rsofun P-model output data} +\format{ +An object of class \code{tbl_df} (inherits from \code{tbl}, \code{data.frame}) with 1 rows and 3 columns. +} +\usage{ +p_model_output +} +\description{ +Example output dataset from a p-model run using \code{\link{p_model_drivers}} +} +\keyword{datasets} diff --git a/man/p_model_output_vcmax25.Rd b/man/p_model_output_vcmax25.Rd new file mode 100644 index 00000000..c205ef18 --- /dev/null +++ b/man/p_model_output_vcmax25.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data.R +\docType{data} +\name{p_model_output_vcmax25} +\alias{p_model_output_vcmax25} +\title{rsofun P-model output data (using vcmax25 drivers)} +\format{ +An object of class \code{tbl_df} (inherits from \code{tbl}, \code{data.frame}) with 4 rows and 3 columns. +} +\usage{ +p_model_output_vcmax25 +} +\description{ +Example output dataset from a p-model run using \code{\link{p_model_drivers_vcmax25}} +} +\keyword{datasets} diff --git a/man/p_model_validation.Rd b/man/p_model_validation.Rd index c742bff1..48113152 100644 --- a/man/p_model_validation.Rd +++ b/man/p_model_validation.Rd @@ -3,7 +3,7 @@ \docType{data} \name{p_model_validation} \alias{p_model_validation} -\title{SOFUN p-model GPP validation data} +\title{rsofun P-model GPP validation data} \format{ A tibble of validation data: \describe{ @@ -27,8 +27,8 @@ Sci Data 7, 225 (2020). https://doi.org/10.1038/s41597-020-0534-3 p_model_validation } \description{ -Small tests dataset to validate -calibration routines for a time series of fluxes. +Small example dataset of target observations (daily GPP flux data) to optimize +model parameters with the function \code{\link{calib_sofun}} } \examples{ require(ggplot2); require(tidyr) diff --git a/man/p_model_validation_vcmax25.Rd b/man/p_model_validation_vcmax25.Rd index 2a4009e5..77bbbef6 100644 --- a/man/p_model_validation_vcmax25.Rd +++ b/man/p_model_validation_vcmax25.Rd @@ -3,7 +3,7 @@ \docType{data} \name{p_model_validation_vcmax25} \alias{p_model_validation_vcmax25} -\title{SOFUN p-model Vcmax25 validation data} +\title{rsofun P-model Vcmax25 validation data} \format{ A tibble of validation data: \describe{ @@ -30,8 +30,8 @@ New Phytol. 206 (2), 614–636. doi:10.1111/nph.13253 p_model_validation_vcmax25 } \description{ -Small tests dataset to validate -calibration routines for leaf traits. +Small example dataset of target observations (leaf trait data) to optimize +model parameters with the function \code{\link{calib_sofun}} } \examples{ require(ggplot2); require(tidyr) diff --git a/man/run_biomee_f_bysite.Rd b/man/run_biomee_f_bysite.Rd index 2339bb4c..8062f92d 100644 --- a/man/run_biomee_f_bysite.Rd +++ b/man/run_biomee_f_bysite.Rd @@ -20,143 +20,25 @@ run_biomee_f_bysite( \item{sitename}{Site name.} \item{params_siml}{Simulation parameters. -\describe{ - \item{spinup}{Flag indicating whether this simulation does spin-up.} - \item{spinupyears}{Number of spin-up years.} - \item{recycle}{Length of standard recycling period (years).} - \item{firstyeartrend}{First transient year.} - \item{nyeartrend}{Number of transient years.} - \item{steps_per_day}{Time resolution (day-1).} - \item{do_U_shaped_mortality}{Flag indicating whether U-shaped - mortality is used.} - \item{update_annualLAImax}{Flag indicating whether updating - LAImax according to mineral N in soil.} - \item{do_closedN_run}{Flag indicating whether doing N closed - runs to recover N balance enforcing 0.2 kg N m-2 in the inorganic N pool.} - \item{code_method_photosynth}{String specifying the method of photosynthesis - used in the model, either "pmodel" or "gs_leuning".document()} - \item{code_method_mortality}{String indicating the type of mortality in the - model. One of the following: "dbh" is size-dependent mortality, "const_selfthin" - is constant self thinning (in development), "cstarvation" is carbon starvation, and - "growthrate" is growth rate dependent mortality.} -}} +See examples \code{\link{biomee_gs_leuning_drivers}} or \code{\link{biomee_p_model_drivers}}} \item{site_info}{Site meta info in a data.frame. -\describe{ - \item{sitename}{Name of the site.} - \item{lon}{Longitud of the site location.} - \item{lat}{Latitude of the site location.} - \item{elv}{Elevation of the site location, in meters.} -}} +See examples \code{\link{biomee_gs_leuning_drivers}} or \code{\link{biomee_p_model_drivers}}} \item{forcing}{Forcing data.frame used as input. -\describe{ - \item{ppfd}{Photosynthetic photon flux densisty(mol s-1 m-2)} - \item{tair}{Air temperature (deg C)} - \item{vpd}{Vapor pressure defficit (Pa)} - \item{rain}{Precipitation (kgH2O m-2 s-1 == mm s-1)} - \item{wind}{Wind velocity (m s-1)} - \item{pair}{Atmoshperic pressure (pa)} - \item{co2}{CO2 athmospheric concentration (ppm)} -}} +See examples \code{\link{biomee_gs_leuning_drivers}} or \code{\link{biomee_p_model_drivers}}} -\item{params_tile}{Tile-level model parameters, into a single row data.frame - with columns: -\describe{ - \item{soiltype}{Integer indicating the type of soil: Sand = 1, LoamySand = 2, - SandyLoam = 3, SiltLoam = 4, FrittedClay = 5, Loam = 6, Clay = 7.} - \item{FLDCAP}{Field capacity (vol/vol). Water remaining in a soil after it - has been thoroughly saturated and allowed to drain freely.} - \item{WILTPT}{Wilting point (vol/vol). Water content of a soil at which - plants wilt and fail to recover.} - \item{K1}{Fast soil C decomposition rate (year\eqn{^{-1}}).} - \item{K2}{Slow soil C decomposition rate (year\eqn{^{-1}}).} - \item{K_nitrogen}{Mineral Nitrogen turnover rate (year\eqn{^{-1}}).} - \item{MLmixRatio}{Ratio of C and N returned to litters from microbes.} - \item{etaN}{N loss rate through runoff (organic and mineral) (year\eqn{^{-1}}).} - \item{LMAmin}{Minimum LMA, leaf mass per unit area, kg C m\eqn{^{-2}}.} - \item{fsc_fine}{Fraction of fast turnover carbon in fine biomass.} - \item{fsc_wood}{Fraction of fast turnover carbon in wood biomass.} - \item{GR_factor}{Growth respiration factor.} - \item{l_fract}{Fraction of the carbon retained after leaf drop.} - \item{retransN}{Retranslocation coefficient of nitrogen.} - \item{f_initialBSW}{Coefficient for setting up initial sapwood.} - \item{f_N_add}{Re-fill of N for sapwood.} - \item{tf_base}{Calibratable scalar for respiration, used to increase LUE - levels.} - \item{par_mort}{Canopy mortality parameter.} - \item{par_mort_under}{Parameter for understory mortality.} -}} +\item{params_tile}{Tile-level model parameters, into a single row data.frame. +See examples \code{\link{biomee_gs_leuning_drivers}} or \code{\link{biomee_p_model_drivers}}} \item{params_species}{A data.frame containing species-specific model parameters, - with one species per row. The columns of this data.frame are: -\describe{ - \item{lifeform}{Integer set to 0 for grasses and 1 for trees.} - \item{phenotype}{Integer set to 0 for deciduous and 1 for evergreen.} - \item{pt}{Integer indicating the type of plant according to photosynthesis: - 0 for C3; 1 for C4} - \item{alpha_FR}{Fine root turnonver rate (year\eqn{^{-1}}).} - \item{rho_FR}{Material density of fine roots (kg C m\eqn{^{-3}}).} - \item{root_r}{Radious of the fine roots, in m.} - \item{root_zeta}{e-folding parameter of root vertical distribution, in m.} - \item{Kw_root}{Fine root water conductivity (mol m\eqn{^{-2}} - s\eqn{^{-1}} MPa\eqn{^{-1}}).} - \item{leaf_size}{Characteristic leaf size.} - \item{Vmax}{Max RuBisCo rate, in mol m\eqn{^{-2}} s\eqn{^{-1}}.} - \item{Vannual}{Annual productivity per unit area at full sun (kg C - m\eqn{^{-2}} year\eqn{^{-2}}).} - \item{wet_leaf_dreg}{Wet leaf photosynthesis down-regulation.} - \item{m_cond}{Factor of stomatal conductance.} - \item{alpha_phot}{Photosynthesis efficiency.} - \item{gamma_L}{Leaf respiration coefficient, in year\eqn{^{-1}}.} - \item{gamma_LN}{Leaf respiration coefficient per unit N.} - \item{gamma_SW}{Sapwood respiration rate, in kg C m\eqn{^{-2}} year\eqn{^{-1}}.} - \item{gamma_FR}{Fine root respiration rate, kg C kg C\eqn{^{-1}} - year\eqn{^{-1}}.} - \item{tc_crit}{Critical temperature triggerng offset of phenology, in Kelvin.} - \item{tc_crit_on}{Critical temperature triggerng onset of phenology, in Kelvin.} - \item{gdd_crit}{Critical value of GDD5 for turning ON growth season.} - \item{betaON}{Critical soil moisture for phenology onset.} - \item{betaOFF}{Critical soil moisture for phenology offset.} - \item{seedlingsize}{Initial size of seedlings, in kg C per individual.} - \item{LNbase}{Basal leaf N per unit area, in kg N m\eqn{^{-2}}.} - \item{lAImax}{Maximum crown LAI (leaf area index).} - \item{Nfixrate0}{Reference N fixation rate (kg N kg C\eqn{^{-1}} root).} - \item{NfixCost0}{Carbon cost of N fixation (kg C kg N\eqn{^{-1}}).} - \item{phiCSA}{Ratio of sapwood area to leaf area.} - \item{mortrate_d_c}{Canopy tree mortality rate (year\eqn{^{-1}}).} - \item{mortrate_d_u}{Understory tree mortality rate (year\eqn{^{-1}}).} - \item{maturalage}{Age at which trees can reproduce (years).} - \item{v_seed}{Fraction of G_SF to G_F.} - \item{fNSmax}{Multiplier for NSNmax as sum of potential bl and br.} - \item{LMA}{Leaf mass per unit area (kg C m\eqn{^{-2}}).} - \item{rho_wood}{Wood density (kg C m\eqn{^{-3}}).} - \item{alphaBM}{Coefficient for allometry (biomass = alphaBM * DBH ** thetaBM).} - \item{thetaBM}{Coefficient for allometry (biomass = alphaBM * DBH ** thetaBM).} - \item{kphio}{Quantum yield efficiency \eqn{\varphi_0}, - in mol mol\eqn{^{-1}}.} - \item{phiRL}{Ratio of fine root to leaf area.} - \item{LAI_light}{Maximum LAI limited by light.} -}} +with one species per row. See examples \code{\link{biomee_gs_leuning_drivers}} or \code{\link{biomee_p_model_drivers}}} \item{init_cohort}{A data.frame of initial cohort specifications. -\describe{ - \item{init_cohort_species}{Index of a species described in param_species.} - \item{init_cohort_nindivs}{Initial individual density, in individuals per - m\eqn{^{2}}.} - \item{init_cohort_bsw}{Initial biomass of sapwood, in kg C per individual.} - \item{init_cohort_bHW}{Initial biomass of heartwood, in kg C per tree.} - \item{init_cohort_nsc}{Initial non-structural biomass.} -}} +See examples \code{\link{biomee_gs_leuning_drivers}} or \code{\link{biomee_p_model_drivers}}} \item{init_soil}{A data.frame of initial soil pools. -\describe{ - \item{init_fast_soil_C}{Initial fast soil carbon, in kg C m\eqn{^{-2}}.} - \item{init_slow_soil_C}{Initial slow soil carbon, in kg C m\eqn{^{-2}}.} - \item{init_Nmineral}{Mineral nitrogen pool, in kg N m\eqn{^{-2}}.} - \item{N_input}{Annual nitrogen input to soil N pool, in kg N m\eqn{^{-2}} - year\eqn{^{-1}}.} -}} +See examples \code{\link{biomee_gs_leuning_drivers}} or \code{\link{biomee_p_model_drivers}}} \item{makecheck}{Flag specifying whether checks are performed to verify model inputs and parameters.} } diff --git a/man/run_pmodel_f_bysite.Rd b/man/run_pmodel_f_bysite.Rd index 59c8d557..59decebd 100644 --- a/man/run_pmodel_f_bysite.Rd +++ b/man/run_pmodel_f_bysite.Rd @@ -35,7 +35,7 @@ run_pmodel_f_bysite( \item{site_info}{A list of site meta info. Required: \describe{ - \item{lon}{Longitud of the site location.} + \item{lon}{Longitude of the site location.} \item{lat}{Latitude of the site location.} \item{elv}{Elevation of the site location, in meters.} \item{whc}{A numeric value for the total root zone water holding capacity (in mm), used diff --git a/man/runread_biomee_f.Rd b/man/runread_biomee_f.Rd index 9551388a..29de7ba4 100644 --- a/man/runread_biomee_f.Rd +++ b/man/runread_biomee_f.Rd @@ -32,8 +32,11 @@ Runs BiomeE model for multiple sites. \donttest{ # Example BiomeE model run -mod_output <- runread_biomee_f( +runread_biomee_f( drivers = biomee_gs_leuning_drivers ) +runread_biomee_f( + drivers = biomee_p_model_drivers +) } } diff --git a/man/runread_pmodel_f.Rd b/man/runread_pmodel_f.Rd index 6c50b38e..7eda5c30 100644 --- a/man/runread_pmodel_f.Rd +++ b/man/runread_pmodel_f.Rd @@ -107,4 +107,7 @@ params_modl <- list( output <- rsofun::runread_pmodel_f( drivers = rsofun::p_model_drivers, par = params_modl) +output_vcmax25 <- rsofun::runread_pmodel_f( + drivers = rsofun::p_model_drivers_vcmax25, + par = params_modl) } diff --git a/src/datatypes.mod.f90 b/src/datatypes.mod.f90 index 63d65099..58796d62 100644 --- a/src/datatypes.mod.f90 +++ b/src/datatypes.mod.f90 @@ -156,12 +156,10 @@ module datatypes type :: vegn_tile_type integer :: n_cohorts = 0.0 - integer :: n_years = 0.0 - integer :: n_canopycc = 0.0 !===== For reseting vegetation integer :: n_initialCC = 0 - type(cohort_type), pointer :: initialCC(:)=>NULL() + type(cohort_type), pointer :: initialCC(:) => NULL() !===== Cohorts nested inside tile type(cohort_type), pointer :: cohorts(:) => NULL() diff --git a/src/params_siml_biomee.mod.f90 b/src/params_siml_biomee.mod.f90 index 9781d2e2..be43019e 100644 --- a/src/params_siml_biomee.mod.f90 +++ b/src/params_siml_biomee.mod.f90 @@ -16,7 +16,6 @@ module md_params_siml_biomee type(steering_parameters) :: steering - ! integer :: model_run_years logical :: do_U_shaped_mortality logical :: update_annualLAImax logical :: do_closedN_run diff --git a/tests/testthat/test-calibration-pmodel.R b/tests/testthat/test-calibration-pmodel.R index c2546abf..c1327655 100644 --- a/tests/testthat/test-calibration-pmodel.R +++ b/tests/testthat/test-calibration-pmodel.R @@ -178,9 +178,9 @@ test_that("test Vcmax25 calibration routine p-model (GenSA, rmse)", { test_that("test joint calibration routine p-model (BT, likelihood maximization)", { skip_on_cran() drivers <- rbind(gpp = rsofun::p_model_drivers, - vcmax25 = rsofun::p_model_drivers_vcmax25) + vcmax25 = rsofun::p_model_drivers_vcmax25) obs <- rbind(gpp = rsofun::p_model_validation, - vcmax25 = rsofun::p_model_validation_vcmax25) + vcmax25 = rsofun::p_model_validation_vcmax25) params_fix <- list( # kphio = 0.04998, # setup ORG in Stocker et al. 2020 GMD kphio_par_a = 0.01, # set to zero to disable temperature-dependence of kphio, setup ORG in Stocker et al. 2020 GMD diff --git a/tests/testthat/test-model-runs.R b/tests/testthat/test-model-runs.R index 0489420d..ca48078a 100644 --- a/tests/testthat/test-model-runs.R +++ b/tests/testthat/test-model-runs.R @@ -41,7 +41,7 @@ test_that("p-model run check GPP", { makecheck = FALSE, parallel = FALSE ) - + # test for correctly returned values expect_type(df_output, "list") @@ -116,16 +116,16 @@ test_that("p-model run check Vcmax25", { test_that("biomeE output check (gs leuning)", { skip_on_cran() - + out <- runread_biomee_f( biomee_gs_leuning_drivers, makecheck = TRUE, parallel = FALSE) - + output_annual_tile <- out$data[[1]]$output_annual_tile - + expect_true(all.equal(colMeans(output_annual_tile), colMeans(biomee_gs_leuning_output$data[[1]]$output_annual_tile), tolerance = 1e-4)) - + # If this test fails it means that the output of the model is out of sync with the data in the data directory. # It could either mean that: # - the model was accidentally altered and should be fixed to deliver the expected output @@ -136,16 +136,16 @@ test_that("biomeE output check (gs leuning)", { test_that("biomee output check (p-model)", { skip_on_cran() - + out <- runread_biomee_f( biomee_p_model_drivers, makecheck = TRUE, parallel = FALSE) - + output_annual_tile <- out$data[[1]]$output_annual_tile - + expect_true(all.equal(colMeans(output_annual_tile), colMeans(biomee_p_model_output$data[[1]]$output_annual_tile), tolerance = 1e-4)) - + # If this test fails it means that the output of the model is out of sync with the data in the data directory. # It could either mean that: # - the model was accidentally altered and should be fixed to deliver the expected output diff --git a/tests/testthat/test-quantitative-validation.R b/tests/testthat/test-quantitative-validation.R index 376b6072..b9321162 100644 --- a/tests/testthat/test-quantitative-validation.R +++ b/tests/testthat/test-quantitative-validation.R @@ -29,7 +29,7 @@ test_that("p-model quantitative check", { # normal tolerance ~ 0.305 tolerance <- mean(abs(output - gpp), na.rm = TRUE)/ - mean(abs(gpp), na.rm = TRUE) + mean(abs(gpp), na.rm = TRUE) # test for correctly returned values expect_equal(tolerance, 0.4201191, tolerance = 0.04) diff --git a/vignettes/biomee_use.Rmd b/vignettes/biomee_use.Rmd index d883a3d3..06c6363a 100644 --- a/vignettes/biomee_use.Rmd +++ b/vignettes/biomee_use.Rmd @@ -19,11 +19,15 @@ library(dplyr) library(ggplot2) ``` -The `rsofun` package and framework includes two main models. The `pmodel` and `biomee` (which in part relies on pmodel component). Here we give a short example on how to run the `biomee` model on the included demo datasets to familiarize yourself with both the data structure and the outputs. +The `rsofun` package and framework includes two distinct simulation models. +The `p-model` and `biomee` (which in part relies on p-model component). +Below we give a short example on how to run the `biomee` model on the included demo datasets to familiarize yourself with both the data structure and the outputs. ## Demo data -The package includes two demo datasets to run. These files can be directly loaded into your workspace by typing: +The package includes two demo datasets to run simulations at the site CH-LAE with either the P-model or Leuning photosynthesis models. +See below under 'Two model approaches' for a brief description of the differences. +These demo files can be directly loaded into your workspace by typing: ```{r eval = FALSE} library(rsofun) @@ -33,7 +37,7 @@ biomee_p_model_drivers biomee_validation ``` -These are real data from the Swiss CH-Lae fluxnet site. We can use these data to run the model, together with observations of GPP we can also parameterize `biomee` parameters. +These are real data from the Swiss CH-LAE fluxnet site. We can use these data to run the model, together with observations of GPP we can also parameterize `biomee` parameters. # Two model approaches @@ -66,9 +70,13 @@ biomee_gs_leuning_output_annual_tile <- out$data[[1]]$output_annual_tile biomee_gs_leuning_output_annual_cohorts <- out$data[[1]]$output_annual_cohorts ``` -```{r echo = FALSE} -biomee_gs_leuning_output_annual_tile <- biomee_gs_leuning_output$data[[1]]$output_annual_tile -biomee_gs_leuning_output_annual_cohorts <- biomee_gs_leuning_output$data[[1]]$output_annual_cohorts +```{r simulate_biomee_gs_leuning_run, include = FALSE} +# TODO: get rid of this and always fully run the vignettes +# fake output since model isn't run +# saveRDS(out, "files/biomee_use.Rmd__biomee_gs_leuning_output___out.RDS") +out <- readRDS("files/biomee_use.Rmd__biomee_gs_leuning_output___out.RDS") +biomee_gs_leuning_output_annual_tile <- out$data[[1]]$output_annual_tile +biomee_gs_leuning_output_annual_cohorts <- out$data[[1]]$output_annual_cohorts ``` ### Plotting output @@ -121,9 +129,13 @@ biomee_p_model_output_annual_tile <- out$data[[1]]$output_annual_tile biomee_p_model_output_annual_cohorts <- out$data[[1]]$output_annual_cohorts ``` -```{r echo = FALSE} -biomee_p_model_output_annual_tile <- biomee_p_model_output$data[[1]]$output_annual_tile -biomee_p_model_output_annual_cohorts <- biomee_p_model_output$data[[1]]$output_annual_cohorts +```{r simulate_biomee_p_model_run, include = FALSE} +# TODO: get rid of this and always fully run the vignettes +# fake output since model isn't run +# saveRDS(out, "files/biomee_use.Rmd__biomee_p_model_output___out.RDS") +out <- readRDS("files/biomee_use.Rmd__biomee_p_model_output___out.RDS") +biomee_p_model_output_annual_tile <- out$data[[1]]$output_annual_tile +biomee_p_model_output_annual_cohorts <- out$data[[1]]$output_annual_cohorts ``` ### Plotting output @@ -180,14 +192,17 @@ pars <- calib_sofun( ) ``` -```{r echo = FALSE} -# Take values from the chunk before, which was run locally +```{r include = FALSE} +# TODO: get rid of this and always fully run the vignettes +# fake output since calibration isn't run +# saveRDS(pars, "files/biomee_use.Rmd__biomee_p_model_output___pars.RDS") +# pars <- readRDS("files/biomee_use.Rmd__biomee_p_model_output___pars.RDS") pars <- list( - par = c(phiRL = 2.0492210, - LAI_light = 4.5462360, - tf_base = 0.5006033, - par_mort = 1.9317278) -) + par = c(phiRL = 1.5597893, #2.0492210, + LAI_light = 4.9657286,#4.5462360, + tf_base = 0.9885088, #0.5006033, + par_mort = 0.1665635 #1.9317278) +)) ``` Using the calibrated parameter values, we can run again the BiomeE simulations. @@ -200,41 +215,40 @@ drivers$params_tile[[1]]$tf_base <- pars$par[3] drivers$params_tile[[1]]$par_mort <- pars$par[4] # run the model with new parameter values -biomee_p_model_output_calib <- runread_biomee_f( +calibrated_out <- runread_biomee_f( drivers, makecheck = TRUE, parallel = FALSE ) # split out the annual data -biomee_p_model_output_calib <- biomee_p_model_output_calib$data[[1]]$output_annual_tile -``` - -```{r echo = FALSE, eval = FALSE} -# Save output -save(biomee_p_model_output_calib, file = "files/biomee_p_model_output_calib.rda") +biomee_p_model_calibratedOutput_annual_tile <- calibrated_out$data[[1]]$output_annual_tile ``` -```{r echo = FALSE} -load("files/biomee_p_model_output_calib.rda") +```{r simulate_biomee_p_model_calibration, include = FALSE} +# TODO: get rid of this and always fully run the vignettes +# fake output since calibrated simulation isn't run +# saveRDS(calibrated_out, "files/biomee_use.Rmd__biomee_p_model_output___calibrated_out.RDS") +calibrated_out <- readRDS("files/biomee_use.Rmd__biomee_p_model_output___calibrated_out.RDS") +biomee_p_model_calibratedOutput_annual_tile <- calibrated_out$data[[1]]$output_annual_tile ``` Finally, we can visually compare the two model runs. We plot the original model run in black and the run using the calibrated parameters in grey. The dashed line represents the validation data, i.e. the observed GPP and Biomass. ```{r fig.width=7} # unnest model output for our single site ggplot() + - geom_line(data = biomee_p_model_output$data[[1]]$output_annual_tile, + geom_line(data = biomee_p_model_output_annual_tile, aes(x = year, y = GPP)) + - geom_line(data = biomee_p_model_output_calib, + geom_line(data = biomee_p_model_calibratedOutput_annual_tile, aes(x = year, y = GPP), color = "grey50") + theme_classic() + labs(x = "Year", y = "GPP") ggplot() + - geom_line(data = biomee_p_model_output$data[[1]]$output_annual_tile, + geom_line(data = biomee_p_model_output_annual_tile, aes(x = year, y = plantC)) + - geom_line(data = biomee_p_model_output_calib, + geom_line(data = biomee_p_model_calibratedOutput_annual_tile, aes(x = year, y = plantC), color = "grey50") + theme_classic() + diff --git a/vignettes/files/biomee_gs_leuning_output.rda b/vignettes/files/biomee_gs_leuning_output.rda deleted file mode 100644 index 400e9f99..00000000 Binary files a/vignettes/files/biomee_gs_leuning_output.rda and /dev/null differ diff --git a/vignettes/files/biomee_p_model_output.rda b/vignettes/files/biomee_p_model_output.rda deleted file mode 100644 index f6006193..00000000 Binary files a/vignettes/files/biomee_p_model_output.rda and /dev/null differ diff --git a/vignettes/files/biomee_p_model_output_calib.rda b/vignettes/files/biomee_p_model_output_calib.rda deleted file mode 100644 index cf66f4a2..00000000 Binary files a/vignettes/files/biomee_p_model_output_calib.rda and /dev/null differ diff --git a/vignettes/files/biomee_use.Rmd__biomee_gs_leuning_output___out.RDS b/vignettes/files/biomee_use.Rmd__biomee_gs_leuning_output___out.RDS new file mode 100644 index 00000000..95f5fae6 Binary files /dev/null and b/vignettes/files/biomee_use.Rmd__biomee_gs_leuning_output___out.RDS differ diff --git a/vignettes/files/biomee_use.Rmd__biomee_p_model_output___calibrated_out.RDS b/vignettes/files/biomee_use.Rmd__biomee_p_model_output___calibrated_out.RDS new file mode 100644 index 00000000..5d2017f6 Binary files /dev/null and b/vignettes/files/biomee_use.Rmd__biomee_p_model_output___calibrated_out.RDS differ diff --git a/vignettes/files/biomee_use.Rmd__biomee_p_model_output___out.RDS b/vignettes/files/biomee_use.Rmd__biomee_p_model_output___out.RDS new file mode 100644 index 00000000..77d53f25 Binary files /dev/null and b/vignettes/files/biomee_use.Rmd__biomee_p_model_output___out.RDS differ diff --git a/vignettes/files/morrisOut.rda b/vignettes/files/morrisOut.rda deleted file mode 100644 index 5fd2aad7..00000000 Binary files a/vignettes/files/morrisOut.rda and /dev/null differ diff --git a/vignettes/files/par_calib.rda b/vignettes/files/par_calib.rda deleted file mode 100644 index 939febe5..00000000 Binary files a/vignettes/files/par_calib.rda and /dev/null differ diff --git a/vignettes/files/pmodel_runs.rda b/vignettes/files/pmodel_runs.rda deleted file mode 100644 index c7d9fa08..00000000 Binary files a/vignettes/files/pmodel_runs.rda and /dev/null differ diff --git a/vignettes/files/sensitivity_analysis.Rmd__morrisOut.RDS b/vignettes/files/sensitivity_analysis.Rmd__morrisOut.RDS new file mode 100644 index 00000000..1e08006d Binary files /dev/null and b/vignettes/files/sensitivity_analysis.Rmd__morrisOut.RDS differ diff --git a/vignettes/files/sensitivity_analysis.Rmd__par_calib.RDS b/vignettes/files/sensitivity_analysis.Rmd__par_calib.RDS new file mode 100644 index 00000000..bf0c2f09 Binary files /dev/null and b/vignettes/files/sensitivity_analysis.Rmd__par_calib.RDS differ diff --git a/vignettes/files/sensitivity_analysis.Rmd__pmodel_runs.RDS b/vignettes/files/sensitivity_analysis.Rmd__pmodel_runs.RDS new file mode 100644 index 00000000..36404bd4 Binary files /dev/null and b/vignettes/files/sensitivity_analysis.Rmd__pmodel_runs.RDS differ diff --git a/vignettes/pmodel_use.Rmd b/vignettes/pmodel_use.Rmd index a3220ca8..393c6b9b 100644 --- a/vignettes/pmodel_use.Rmd +++ b/vignettes/pmodel_use.Rmd @@ -20,17 +20,15 @@ knitr::opts_chunk$set( library(rsofun) library(dplyr) library(ggplot2) - -# fake variable as optimization isn't run -pars <- list() -pars$par["kphio"] <- 0.04478049 ``` -The `rsofun` package and framework includes two main models. The `pmodel` and `biomee` (which in part relies on P-model components). Here we give a short example on how to run the `pmodel` on the included demo datasets to familiarize yourself with both the data structure and the outputs. +The `rsofun` package and framework includes two distinct simulation models. +The `p-model` and `biomee` (which in part relies on p-model component). +Here we give a short example on how to run the `p-model` on the included demo datasets to familiarize yourself with both the data structure and the outputs. ## Demo data -The package includes two demo datasets to run and validate pmodel output using GPP observations. These files can be directly loaded into your workspace by typing: +The package includes two demo datasets to run and validate p-model output using GPP observations. These files can be directly loaded into your workspace by typing: ```{r} library(rsofun) @@ -41,7 +39,7 @@ p_model_drivers p_model_validation ``` -These are real data from the French FR-Pue fluxnet site. Information about data structure, variable names, and their meaning and units can be found in the reference pages of `p_model_drivers` and `p_model_validation`. We can use these data to run the model, together with observations of GPP we can also calibrate `pmodel` parameters. +These are real data from the French FR-Pue fluxnet site. Information about data structure, variable names, and their meaning and units can be found in the reference pages of `p_model_drivers` and `p_model_validation`. We can use these data to run the model, together with observations of GPP we can also calibrate `p-model` parameters. Another two datasets are provided as an example to validate the model against leaf traits data, rather than fluxes. Measurements of Vcmax25 (aggregated over species) for a subset of 4 sites from the GlobResp database (Atkin et al., 2015) are given in `p_model_validation_vcmax25` and the corresponding forcing for the P-model is given in `p_model_drivers_vcmax25`. Since leaf traits are only measured once per site, the forcing used is a single year of average climate (the average measurements between 2001 and 2015 on each day of the year). @@ -163,6 +161,11 @@ pars <- calib_sofun( ) ``` +```{r simulate_calibration_run, include = FALSE} +# fake variable as optimization isn't run +pars <- list() +pars$par["kphio"] <- 0.03580962 +``` When successful the optimized parameters can be used to run subsequent modelling efforts, in this case slightly improving the model fit over a more global parameter set. ```{r} diff --git a/vignettes/sensitivity_analysis.Rmd b/vignettes/sensitivity_analysis.Rmd index 7a0b5146..9e584de8 100644 --- a/vignettes/sensitivity_analysis.Rmd +++ b/vignettes/sensitivity_analysis.Rmd @@ -160,14 +160,11 @@ morrisOut <- sensitivity::morris( scale = TRUE) ``` -```{r eval = FALSE, echo = FALSE} -# Save Morris sensitivity output because it takes very long to compute -save(morrisOut, file = "files/morrisOut.rda") -``` - -```{r eval = TRUE, echo = FALSE} -# Load Morris sensitivity output -load("files/morrisOut.rda") +```{r simulate_morrisOut, include = FALSE} +# TODO: get rid of this and always fully run the vignettes +# fake output since Morris sensitivity isn't run +# saveRDS(morrisOut, "files/sensitivity_analysis.Rmd__morrisOut.RDS") +morrisOut <- readRDS("files/sensitivity_analysis.Rmd__morrisOut.RDS") ``` The analysis evaluates the variability of the target function, i.e. the @@ -238,6 +235,7 @@ It is always important to check the convergence of the MCMC algorithm used for t According to the previous sensitivity analysis, calibrating the error parameter for GPP and the quantum yield efficiency parameters will have a high impact on the model fit. Let's run the calibration: ```{r eval = FALSE, echo = TRUE} +# Calibrates kphio, kphio_par_b, kc_jmax - top 3 model params set.seed(2023) # Define calibration settings @@ -254,76 +252,37 @@ settings_calib <- list( )), par = list( kphio = list(lower = 0.03, upper = 0.15, init = 0.05), - kphio_par_a = list(lower = -0.004, upper = -0.001, init = -0.0025), + #kphio_par_a = list(lower = -0.004, upper = -0.001, init = -0.0025), kphio_par_b = list(lower = 10, upper = 30, init =25), + kc_jmax = list(lower = 0.2, upper = 0.8, init = 0.41), err_gpp = list(lower = 0.1, upper = 3, init = 0.8) ) ) - -# Calibrate kphio-related parameters and err_gpp -par_calib <- calib_sofun( - drivers = p_model_drivers, - obs = p_model_validation, - settings = settings_calib, - par_fixed = list( +par_fixed <- list( soilm_thetastar = 0.6*240, soilm_betao = 0.2, beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, - tau_acclim = 30.0, - kc_jmax = 0.41), - targets = "gpp" -) - -# This code takes 15 minutes to run -``` - -```{r eval = FALSE, echo = FALSE} -# Calibrates kphio, betao, kc_jmax - top 3 model params -set.seed(2023) - -# Define calibration settings -settings_calib <- list( - method = "BayesianTools", - metric = rsofun::cost_likelihood_pmodel, - control = list( - sampler = "DEzs", - settings = list( - burnin = 1500, - iterations = 6000, - nrChains = 3 # number of chains to be sampled - )), - par = list( - kphio = list(lower = 0.03, upper = 0.15, init = 0.05), - soilm_betao = list(lower = 0, upper = 1, init = 0.2), - kc_jmax = list(lower = 0.2, upper = 0.8, init = 0.41), - err_gpp = list(lower = 0.1, upper = 3, init = 0.8) - ) -) + kphio_par_a = -0.0025, + tau_acclim = 30.0) +# Calibrate kphio-related parameters and err_gpp par_calib <- calib_sofun( drivers = p_model_drivers, obs = p_model_validation, settings = settings_calib, - par_fixed = list( - kphio_par_a = -0.0025, - kphio_par_b = 20, - soilm_thetastar = 0.6*240, - beta_unitcostratio = 146.0, - rd_to_vcmax = 0.014, - tau_acclim = 30.0), + par_fixed = par_fixed, targets = "gpp" ) -``` -```{r eval = FALSE, echo = FALSE} -# Save calibration output because it takes very long to compute -save(par_calib, file = "files/par_calib.rda") +# This code takes 15 minutes to run ``` -```{r eval = TRUE, echo = FALSE} -# Load calibration output -load("files/par_calib.rda") +```{r simulate_par_calib, include = FALSE} +# TODO: get rid of this and always fully run the vignettes +# fake output since calibration isn't run +# saveRDS(par_calib, "files/sensitivity_analysis.Rmd__par_calib.RDS") +par_calib <- readRDS("files/sensitivity_analysis.Rmd__par_calib.RDS") ``` `BayesianTools` makes it easy to produce the trace plot of the MCMC chains and the posterior density plot for the parameters. Trace plots show the time series of the sampled chains, which should reach a stationary state. One can also choose a burnin visually, to discard the early iterations and keep only the samples from the stationary distribution to which they converge. We set \code{burnin = 3000} above from previous runs, and those iterations are not shown by the following trace plot. The samples after the burnin period should be used for inference. @@ -363,7 +322,7 @@ plot_acf_mcmc <- function(chains, par_names){ } } } -plot_acf_mcmc(par_calib$mod, c("kphio", "kphio_par_a", "kphio_par_b", "err_gpp")) +plot_acf_mcmc(par_calib$mod, c("kphio", "kc_jmax", "kphio_par_b", "err_gpp")) ``` Looking at the correlation between chains for different parameters is also helpful because parameter correlation may slow down convergence, or the chains may oscillate in the multivariate posterior space. In this calibration we expect parameter samples to be somewhat correlated, especially `kphio_par_a` and `kphio_par_b` because they specify the shape of the temperature dependence of the quantum yield efficiency, $\varphi_o(T)$. We can also see that `err_gpp` is correlated with `kphio` (to which the P-model is very sensitive), since the error represents how good the model fits the observed GPP. @@ -399,9 +358,7 @@ incorporates the Gaussian model error. # Evaluation of the uncertainty coming from the model parameters' uncertainty # Sample parameter values from the posterior distribution -samples_par <- getSample(par_calib$mod, - thin = 30, # get every 30th sample - whichParameters = 1:4) |> +samples_par <- getSample(par_calib$mod, numSamples = 600) |> as.data.frame() |> dplyr::mutate(mcmc_id = 1:n()) |> tidyr::nest(.by = mcmc_id, .key = "pars") @@ -409,19 +366,9 @@ samples_par <- getSample(par_calib$mod, run_pmodel <- function(sample_par){ # Function that runs the P-model for a sample of parameters # and also adds the new observation error - out <- runread_pmodel_f( drivers = p_model_drivers, - par = list( # copied from par_fixed above - kphio = sample_par$kphio, - kphio_par_a = sample_par$kphio_par_a, - kphio_par_b = sample_par$kphio_par_b, - soilm_thetastar = 0.6*240, - soilm_betao = 0.2, - beta_unitcostratio = 146.0, - rd_to_vcmax = 0.014, - tau_acclim = 30.0, - kc_jmax = 0.41) # value from posterior + par = c(par_fixed, sample_par) # err_gpp is also here, but simply ignored ) # return modelled GPP and prediction for a new GPP observation @@ -450,14 +397,11 @@ pmodel_runs <- samples_par |> ) ``` -```{r eval = FALSE, echo = FALSE} -save(pmodel_runs, file = "files/pmodel_runs.rda") -``` - -```{r echo = FALSE, eval = TRUE} -load("files/pmodel_runs.rda") - -pmodel_runs <- pmodel_runs %>% dplyr::rename(gpp_q50 = gpp) # Temporary fix until rda is updated +```{r simulate_pmodel_runs, include = FALSE} +# TODO: get rid of this and always fully run the vignettes +# fake output since calibration isn't run +# saveRDS(pmodel_runs, "files/sensitivity_analysis.Rmd__pmodel_runs.RDS") +pmodel_runs <- readRDS("files/sensitivity_analysis.Rmd__pmodel_runs.RDS") ``` Below we plot the first year of observed GPP (in black) against the predicted @@ -471,40 +415,111 @@ captured by the predictive interval. ```{r fig.width=7, fig.height=5} # Plot the credible intervals computed above # for the first year only -plot_gpp_error <- ggplot(data = pmodel_runs |> - dplyr::slice(1:365)) + # Plot only first year +data_to_plot <- pmodel_runs |> + # Plot only first year + dplyr::slice(1:365) |> + dplyr::left_join( + # Merge GPP validation data (first year) + p_model_validation$data[[1]][1:365, ] |> + dplyr::rename(gpp_obs = gpp), + by = "date") +# TODO: clean below +# plot_gpp_error <- ggplot(data = data_to_plot) + +# # geom_ribbon( +# # aes(ymin = gpp_pred_q05, +# # ymax = gpp_pred_q95, +# # x = date), +# # fill = 'grey40', alpha = 0.2) + +# geom_ribbon( +# aes(ymin = gpp_q05, +# ymax = gpp_q95, +# x = date), +# fill = 'blue', alpha = 0.5) + +# geom_ribbon( +# aes(ymin = gpp_pred_q05, +# ymax = gpp_pred_q95, +# x = date), +# fill = 'grey40', alpha = 0.2) + +# geom_line( +# aes(x = date, +# y = gpp_q50), +# colour = "grey40", +# alpha = 0.8) + +# # # Include observations in the plot +# geom_point(shape = 3, +# aes(x = date, +# y = gpp_obs),# color = "Observations"), +# alpha = 0.8) + +# theme_classic() + +# theme(panel.grid.major.y = element_line(), +# legend.position = "bottom") + +# labs( +# x = 'Date', +# y = expression(paste("GPP (g C m"^-2, "s"^-1, ")")) +# ) +# plot_gpp_error + +plot_gpp_error <- ggplot(data = data_to_plot) + + geom_ribbon( + aes(ymin = gpp_pred_q05, + ymax = gpp_pred_q95, + x = date, + fill = "Model uncertainty")) + geom_ribbon( aes(ymin = gpp_q05, ymax = gpp_q95, - x = date), - fill = 'blue', alpha = 0.5) + - geom_ribbon( - aes(ymin = gpp_pred_q05, - ymax = gpp_pred_q95, - x = date), - fill = 'grey40', alpha = 0.2) + + x = date, + fill = "Parameter uncertainty")) + geom_line( aes(x = date, - y = gpp_q50), - colour = "grey40", - alpha = 0.8) + + y = gpp_q50, + color = "Predictions")) + + # Include observations in the plot + geom_point(shape = 3, + aes(x = date, + y = gpp_obs, + color = "Observations")) + theme_classic() + - theme(panel.grid.major.y = element_line()) + + theme(panel.grid.major.y = element_line(), + legend.position = "bottom") + labs( x = 'Date', y = expression(paste("GPP (g C m"^-2, "s"^-1, ")")) ) -# Define GPP validation data (first year) -validation_data <- p_model_validation$data[[1]][1:365, ] +t_col <- function(color, percent = 50, name = NULL) { + # color = color name + # percent = % transparency + # name = an optional name for the color + + ## Get RGB values for named color + rgb.val <- col2rgb(color) + + ## Make new color using input color as base and alpha set by transparency + t.col <- rgb(rgb.val[1], rgb.val[2], rgb.val[3], + max = 255, + alpha = (100 - percent) * 255 / 100, + names = name) + + ## Save the color + invisible(t.col) +} -# Include observations in the plot -plot_gpp_error + - geom_line( - data = validation_data, - aes(x = date, - y = gpp), - alpha = 0.8) +plot_gpp_error <- plot_gpp_error + + scale_color_manual(name = "", + breaks = c("Observations", + "Predictions", + "Non-calibrated predictions"), + values = c(t_col("black", 10), + t_col("#E69F00", 10), + t_col("#56B4E9", 10))) + + scale_fill_manual(name = "", + breaks = c("Model uncertainty", + "Parameter uncertainty"), + values = c(t_col("#E69F00", 60), + t_col("#009E73", 40))) +plot_gpp_error +#dir.create("./analysis/paper_results_files2") ``` @@ -515,22 +530,13 @@ plot_gpp_error + ```{r fig.width=7, fig.height=5, echo = FALSE, eval = FALSE} # # Plot observed and predicted GPP, with a 95% confidence interval using err_gpp -plot_gpp_error <- ggplot(data = runread_pmodel_f( - drivers = p_model_drivers, - par = list( - kphio = par_calib$par[1], # values from posterior - kphio_par_a = par_calib$par[2], - kphio_par_b = par_calib$par[3], - soilm_thetastar = 0.6*240, # copied from par_fixed above - soilm_betao = 0.2, - beta_unitcostratio = 146.0, - rd_to_vcmax = 0.014, - tau_acclim = 30.0, - kc_jmax = 0.41) - ) |> +data_to_plot_MAP <- runread_pmodel_f(drivers = p_model_drivers, + par = c(par_fixed, par_calib$par)) |> dplyr::select(sitename, data) |> tidyr::unnest(data) |> - dplyr::slice(1:365)) + # Plot only first year + dplyr::slice(1:365) + +plot_gpp_error <- ggplot(data = data_to_plot_MAP) + # Plot only first year geom_ribbon( aes(ymin = gpp - 2*par_calib$par[4], ymax = gpp + 2*par_calib$par[4], @@ -556,7 +562,7 @@ validation_data <- p_model_validation$data[[1]][1:365, ] # Include observations in the plot plot_gpp_error + - geom_line( + geom_point( data = validation_data, aes( date,