diff --git a/DESCRIPTION b/DESCRIPTION index e0c95db96..6038742d2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: rsofun Title: The P-Model and BiomeE Modelling Framework -Version: 5.1.0 +Version: 5.1.0.9000 Authors@R: c( person( family = "Stocker", @@ -82,6 +82,7 @@ Suggests: ggplot2, knitr, sensitivity, + rpmodel, gridExtra, patchwork VignetteBuilder: knitr diff --git a/NAMESPACE b/NAMESPACE index 65bc96355..599002e41 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -8,6 +8,7 @@ export(cost_rmse_biomee) export(cost_rmse_pmodel) export(run_biomee_f_bysite) export(run_pmodel_f_bysite) +export(run_pmodel_onestep_f_bysite) export(runread_biomee_f) export(runread_pmodel_f) import(BayesianTools) diff --git a/NEWS.md b/NEWS.md index 28fb50335..c02240591 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +# rsofun (development version) + +* New `run_pmodel_onestep_f_bysite()` adds single-step leaf-level simulations of the P-model. + # rsofun 5.1.0 ## New features diff --git a/R/calib_sofun.R b/R/calib_sofun.R index fcd572fdd..bf066b011 100644 --- a/R/calib_sofun.R +++ b/R/calib_sofun.R @@ -95,18 +95,18 @@ calib_sofun <- function( lower <- upper <- out_optim <- NULL # check input variables - if(missing(obs) | missing(drivers) | missing(settings)){ + if (missing(obs) | missing(drivers) | missing(settings)){ stop("missing input arguments, please check all parameters") } # check data structure - if(is.data.frame(obs)){ + if (is.data.frame(obs)){ if (nrow(obs) == 0){ warning("no validation data available, returning NA parameters") return(lapply(settings$par, function(x) NA)) } - }else{ + } else { stop("obs must be a (nested) data.frame") } @@ -131,9 +131,10 @@ calib_sofun <- function( drivers = drivers, ... ) - if(optim_out){ + + if (optim_out){ out_optim <- list(par = out$par, mod = out) - }else{ + } else { out_optim <- list(par = out$par) } @@ -163,14 +164,15 @@ calib_sofun <- function( # setup the bayes run, no message forwarding is provided # so wrap the function in a do.call setup <- BayesianTools::createBayesianSetup( - likelihood = function( - random_par) { - do.call("cost", - list( - par = random_par, - obs = obs, - drivers = drivers - )) + likelihood = function(random_par){ + do.call( + "cost", + list( + par = random_par, + obs = obs, + drivers = drivers + ) + ) }, prior = priors, names = names(settings$par) @@ -189,9 +191,10 @@ calib_sofun <- function( # drop last value bt_par <- BayesianTools::MAP(out)$parametersMAP bt_par <- bt_par[1:(length(bt_par))] - if(optim_out){ + + if (optim_out){ out_optim <- list(par = bt_par, mod = out) - }else{ + } else { out_optim <- list(par = bt_par) } diff --git a/R/data.R b/R/data.R index 4c9545c63..d26a04360 100644 --- a/R/data.R +++ b/R/data.R @@ -26,7 +26,7 @@ #' \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{co2}{Atmospheric CO\eqn{_2} concentration in ppm.} #' \item{ccov}{Cloud coverage in \%. This is only used when either PPFD or #' net radiation are not prescribed.} #' } @@ -214,13 +214,13 @@ #' }} #' \item{forcing}{Forcing data.frame used as input #' \describe{ -#' \item{ppfd}{Photosynthetic photon flux density (mol s-1 m-2)} +#' \item{ppfd}{Photosynthetic photon flux density (mol m-2 s-1)} #' \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{pair}{Atmospheric pressure (Pa)} +#' \item{co2}{Atmospheric CO\eqn{_2} concentration in ppm.} #' }} #' \item{params_tile}{Tile-level model parameters, into a single row data.frame, including the following data: #' \describe{ diff --git a/R/init_dates_dataframe.R b/R/init_dates_dataframe.R deleted file mode 100644 index 84a4baf5c..000000000 --- a/R/init_dates_dataframe.R +++ /dev/null @@ -1,101 +0,0 @@ -#' Initialises a tibble with dates -#' -#' 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. -#' @param yrend An integer defining the end year of dates -#' covered by the dataframe. -#' @param startmoy An integer defining the start month-of-year of dates -#' covered by the dataframe. Defaults to 1. -#' @param startdoy An integer defining the start day-of-year of -#' dates covered by the dataframe. Defaults to 1. -#' @param freq A character string specifying the time steps of dates -#' (in rows). Defaults to \code{"days"}. Any of \code{"days", "months", "years"}. If -#' \code{freq = "months"} the 15\eqn{^{th}} day of the months is used as date, -#' and if \code{freq = "years"} the 1\eqn{^{st}} of January of each year is returned. -#' @param endmoy An integer defining the end month-of-year of dates covered -#' by the dataframe. Defaults to 12. -#' @param enddom An integer defining the end day-of-year of dates -#' covered by the dataframe. Defaults to 31. -#' @param noleap Whether leap years are ignored, that is, whether the 29\eqn{^{th}} -#' of February is removed. Defaults to \code{FALSE}. -#' -#' @return A tibble with dates. -#' - -init_dates_dataframe <- function( - yrstart, - yrend, - startmoy=1, - startdoy=1, - freq="days", - endmoy=12, - enddom=31, - noleap=FALSE ){ - - if (freq=="days"){ - - start_date <- as.Date( - sprintf("%04d-%02d-01", - yrstart, startmoy)) + (startdoy - 1) - - end_date <- as.Date( - sprintf("%04d-%02d-%02d", - yrend, endmoy, enddom)) - - } else if (freq=="months"){ - - start_date <- as.Date( - sprintf("%04d-%02d-15", - yrstart, startmoy)) - - end_date <- as.Date( - sprintf("%04d-%02d-15", - yrend, endmoy)) - - } else if (freq=="years"){ - - start_date <- as.Date( - sprintf("%04d-%02d-01", - yrstart, 1)) - - end_date <- as.Date( - sprintf("%04d-%02d-01", - yrend, 7)) - } - - # define date range - date_range <- data.frame( - date = seq.Date( - from = start_date, - to = end_date, - by = freq - )) - - # convert to decimal date - 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) - ) - } - 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") - ) - - } - - return(date_range) -} diff --git a/R/run_biomee_f_bysite.R b/R/run_biomee_f_bysite.R index ca4efed5a..952a41585 100644 --- a/R/run_biomee_f_bysite.R +++ b/R/run_biomee_f_bysite.R @@ -271,13 +271,13 @@ run_biomee_f_bysite <- function( n_annual_trans = as.integer(n_annual_trans) ) - out <- build_out(biomeeout, init_lu$name, sitename, params_siml$do_daily_diagnostics) + out <- build_out_biomee(biomeeout, init_lu$name, sitename, params_siml$do_daily_diagnostics) return(out) } # Build R output -build_out <- function(biomeeout, lu_names, sitename, do_daily_diagnostics){ +build_out_biomee <- function(biomeeout, lu_names, sitename, do_daily_diagnostics){ if (do_daily_diagnostics) { # If simulation is very long, output gets massive. # E.g., In a 3000 years-simulation 'biomeeout' is 11.5 GB. diff --git a/R/run_pmodel_f_bysite.R b/R/run_pmodel_f_bysite.R index a74bc0bf5..4bfe7de0f 100644 --- a/R/run_pmodel_f_bysite.R +++ b/R/run_pmodel_f_bysite.R @@ -1,6 +1,6 @@ -#' Run P-model (R wrapper) +#' Run P-model (time series) #' -#' Run P-model Fortran model on single site. +#' Run P-model on a single site for a forcing time series. #' #' @param sitename Site name. #' @param params_siml Simulation parameters. @@ -29,27 +29,31 @@ #' \item{\code{le}}{Latent heat flux (in J m\eqn{^{-2}} d\eqn{^{-1}}).} #' \item{\code{pet}}{Potential evapotranspiration (PET), calculated by SPLASH following Priestly-Taylor (in mm d\eqn{^{-1}}).} #' \item{\code{vcmax}}{Maximum rate of RuBisCO carboxylation -#' (Vcmax) (in mol C m\eqn{^{-2}} d\eqn{^{-1}}).} +#' (Vcmax) (in mol C m\eqn{^{-2}} s\eqn{^{-1}}).} #' \item{\code{jmax}}{Maximum rate of electron transport for RuBP regeneration #' (in mol CO\eqn{_2} m\eqn{^{-2}} s\eqn{^{-1}}).} #' \item{\code{vcmax25}}{Maximum rate of carboxylation (Vcmax), -#' normalised to 25\eqn{^o}C (in mol C m\eqn{^{-2}} d\eqn{^{-1}}).} +#' normalised to 25\eqn{^o}C (in mol C m\eqn{^{-2}} s\eqn{^{-1}}).} #' \item{\code{jmax25}}{Maximum rate of electron transport, normalised to #' 25\eqn{^o}C (in mol C m\eqn{^{-2}} s\eqn{^{-1}}).} #' \item{\code{gs_accl}}{Acclimated stomatal conductance (in -#' mol C m\eqn{^{-2}} d\eqn{^{-1}} Pa\eqn{^{-1}}).} +#' mol C (mol photons)\eqn{^{-1}} Pa\eqn{^{-1}}. (Multiply by +#' ppfd (mol photons m\eqn{^{-2}} d\eqn{^{-1}}) and fapar +#' to express per unit ground area and time.)} #' \item{\code{wscal}}{Relative soil water content, between 0 (permanent wilting #' point, PWP) and 1 (field capacity, FC).} #' \item{\code{chi}}{Ratio of leaf-internal to ambient CO\eqn{_{2}}, ci:ca (unitless).} -#' \item{\code{iwue}}{Intrinsic water use efficiency (iWUE) (in Pa).} -#' \item{\code{rd}}{Dark respiration (Rd) in gC m\eqn{^{-2}} d\eqn{^{-1}}.} +#' \item{\code{iwue}}{Intrinsic water use efficiency (iWUE) (unitless, +#' multiply with patm (Pa) to get iWUE in Pa).} +#' \item{\code{rd}}{Dark respiration (Rd) in gC m\eqn{^{-2}} s\eqn{^{-1}}. +#' (Multiply by 1/12 (mol C / gC) to convert to mol C m\eqn{^{-2}} s\eqn{^{-1}}.)} #' \item{\code{tsoil}}{Soil temperature, in \eqn{^{o}}C.} #' \item{\code{netrad}}{Net radiation, in W m\eqn{^{-2}}. WARNING: this is currently ignored as a model forcing. Instead, net radiation is internally calculated by SPLASH.} #' \item{\code{wcont}}{Soil water content, in mm.} #' \item{\code{snow}}{Snow water equivalents, in mm.} #' \item{\code{cond}}{Water input by condensation, in mm d\eqn{^{-1}}} -#' } -#' +#' } +#' #' @details Depending on the input model parameters, it's possible to run the #' different P-model setups presented in Stocker et al. 2020 GMD. The P-model #' version implemented in this package allows more flexibility than the one @@ -121,7 +125,7 @@ run_pmodel_f_bysite <- function( ){ # predefine variables for CRAN check compliance - ccov <- fsun <- . <- NULL + ccov <- fsun <- NULL # base state, always execute the call continue <- TRUE @@ -217,6 +221,12 @@ run_pmodel_f_bysite <- function( 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 + correspond to full years.") + continue <- FALSE + } # simulation parameters to check check_param <- c( "spinup", @@ -231,10 +241,10 @@ run_pmodel_f_bysite <- function( "lgr4" ) - parameter_integrity <- lapply(check_param, function(check_var){ - if (any(is.nanull(params_siml[check_var]))){ + parameter_integrity <- lapply(check_param, function(check_par){ + if (any(is.nanull(params_siml[check_par]))){ warning(sprintf("Error: Missing value in %s for %s", - check_var, sitename)) + check_par, sitename)) return(FALSE) } else { return(TRUE) @@ -244,19 +254,12 @@ 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 - correspond to full years.") - continue <- FALSE - } - + # model parameters to check - 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') + 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') ) != 9){ warning(" Returning a dummy data frame. Incorrect model parameters.") continue <- FALSE @@ -273,29 +276,25 @@ run_pmodel_f_bysite <- function( 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)){ + if (! (in_ppfd & in_netrad)){ # fsun must be available when one of ppfd or netrad is missing - if(any(is.na(forcing$fsun))) continue <- FALSE + if (any(is.na(forcing$fsun))) continue <- FALSE } } - if(continue){ - - + if (continue){ ## C wrapper call - out <- .Call( - + pmodelout <- .Call( 'pmodel_f_C', - + secs_per_tstep = as.integer(secs_per_tstep), + in_ppfd = as.logical(in_ppfd), + in_netrad = as.logical(in_netrad), ## Simulation parameters spinup = as.logical(params_siml$spinup), spinupyears = as.integer(params_siml$spinupyears), recycle = as.integer(params_siml$recycle), firstyeartrend = as.integer(params_siml$firstyeartrend), nyeartrend = as.integer(params_siml$nyeartrend), - secs_per_tstep = as.integer(secs_per_tstep), - in_ppfd = as.logical(in_ppfd), - in_netrad = as.logical(in_netrad), outdt = as.integer(params_siml$outdt), ltre = as.logical(params_siml$ltre), ltne = as.logical(params_siml$ltne), @@ -320,67 +319,169 @@ run_pmodel_f_bysite <- function( as.numeric(params_modl$kc_jmax)), forcing = as.matrix(forcing) ) + } else { + pmodelout <- array(dim = c(1,19), data = NA_real_) + } + + out <- build_out_pmodel(pmodelout, params_siml$firstyeartrend, params_siml$nyeartrend) + + return(out) +} + + +# Build R output +build_out_pmodel <- function(pmodelout, firstyeartrend, nyeartrend){ + # predefine variables for CRAN check compliance + . <- NULL + + # Prepare output to be a nice looking tidy data frame (tibble) + ddf <- init_dates_dataframe( + yrstart = firstyeartrend, + yrend = firstyeartrend + nyeartrend - 1, + noleap = TRUE) + + out <- pmodelout %>% + as.matrix() %>% + as.data.frame() %>% + stats::setNames( + c("fapar", + "gpp", + "aet", + "le", + "pet", + "vcmax", + "jmax", + "vcmax25", + "jmax25", + "gs_accl", + "wscal", + "chi", + "iwue", + "rd", + "tsoil", + "netrad", + "wcont", + "snow", + "cond") + ) %>% + as_tibble(.name_repair = "check_unique") %>% + dplyr::bind_cols(ddf, .) + + if (all(is.na(pmodelout))){ + # return single row output + out <- out[1,] + out$date <- as.Date("2000-01-01") + out$year_dec <- 2000.000 + } else { + # return full output + } + + return(out) +} + +#' Initialises a tibble with dates +#' +#' 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. +#' @param yrend An integer defining the end year of dates +#' covered by the dataframe. +#' @param startmoy An integer defining the start month-of-year of dates +#' covered by the dataframe. Defaults to 1. +#' @param startdoy An integer defining the start day-of-year of +#' dates covered by the dataframe. Defaults to 1. +#' @param freq A character string specifying the time steps of dates +#' (in rows). Defaults to \code{"days"}. Any of \code{"days", "months", "years"}. If +#' \code{freq = "months"} the 15\eqn{^{th}} day of the months is used as date, +#' and if \code{freq = "years"} the 1\eqn{^{st}} of January of each year is returned. +#' @param endmoy An integer defining the end month-of-year of dates covered +#' by the dataframe. Defaults to 12. +#' @param enddom An integer defining the end day-of-year of dates +#' covered by the dataframe. Defaults to 31. +#' @param noleap Whether leap years are ignored, that is, whether the 29\eqn{^{th}} +#' of February is removed. Defaults to \code{FALSE}. +#' +#' @return A tibble with dates. +#' + +init_dates_dataframe <- function( + yrstart, + yrend, + startmoy=1, + startdoy=1, + freq="days", + endmoy=12, + enddom=31, + noleap=FALSE ){ + + if (freq=="days"){ - # Prepare output to be a nice looking tidy data frame (tibble) - ddf <- init_dates_dataframe( - yrstart = params_siml$firstyeartrend, - yrend = params_siml$firstyeartrend + params_siml$nyeartrend - 1, - noleap = TRUE) + start_date <- as.Date( + sprintf("%04d-%02d-01", + yrstart, startmoy)) + (startdoy - 1) - out <- out %>% - as.matrix() %>% - as.data.frame() %>% - stats::setNames( - c("fapar", - "gpp", - "aet", - "le", - "pet", - "vcmax", - "jmax", - "vcmax25", - "jmax25", - "gs_accl", - "wscal", - "chi", - "iwue", - "rd", - "tsoil", - "netrad", - "wcont", - "snow", - "cond") - ) %>% - as_tibble(.name_repair = "check_unique") %>% - dplyr::bind_cols(ddf,.) + end_date <- as.Date( + sprintf("%04d-%02d-%02d", + yrend, endmoy, enddom)) - } else { - out <- tibble(date = as.Date("2000-01-01"), - fapar = NA, - gpp = NA, - transp = NA, - latenth = NA, - pet = NA, - vcmax = NA, - jmax = NA, - vcmax25 = NA, - jmax25 = NA, - gs_accl = NA, - wscal = NA, - chi = NA, - iwue = NA, - rd = NA, - tsoil = NA, - netrad = NA, - wcont = NA, - snow = NA, - cond = NA) + } else if (freq=="months"){ + + start_date <- as.Date( + sprintf("%04d-%02d-15", + yrstart, startmoy)) + + end_date <- as.Date( + sprintf("%04d-%02d-15", + yrend, endmoy)) + + } else if (freq=="years"){ + + start_date <- as.Date( + sprintf("%04d-%02d-01", + yrstart, 1)) + + end_date <- as.Date( + sprintf("%04d-%02d-01", + yrend, 7)) } - return(out) + # define date range + date_range <- data.frame( + date = seq.Date( + from = start_date, + to = end_date, + by = freq + )) + + # convert to decimal date + 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) + ) + } + 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") + ) + + } + return(date_range) } + .onUnload <- function(libpath) { library.dynam.unload("rsofun", libpath) } diff --git a/R/run_pmodel_onestep_f_bysite.R b/R/run_pmodel_onestep_f_bysite.R new file mode 100644 index 000000000..5004a41d5 --- /dev/null +++ b/R/run_pmodel_onestep_f_bysite.R @@ -0,0 +1,190 @@ +#' Run P-model (single time step) +#' +#' Run P-model on a single site for a single time step. This does not include the simulation of ecosystem-level quantities, +#' water limitation, nor a simulation of water fluxes. Instead, this corresponds to a leaf-level representation of the +#' acclimation of photosynthesis. +#' +#' @param lc4 Locigical specifying whether P-model simulation is for C4 (as opposed to C3). Defaults to \code{FALSE}. +#' @param forcing A data frame of forcing climate data, used as input (single row). +#' @param params_modl A named list of free (calibratable) model parameters. See \code{\link{runread_pmodel_f}} +#' @param makecheck A logical specifying whether checks are performed +#' to verify forcings and model parameters. \code{TRUE} by default. +#' +#' For further specifications of above inputs and examples see \code{\link{p_model_drivers}} or \code{\link{p_model_drivers_vcmax25}} + +#' @import dplyr +#' +#' @returns Model output is provided as a tidy dataframe, with columns: +#' \describe{ +#' \item{\code{vcmax}}{Maximum rate of RuBisCO carboxylation +#' (Vcmax) (in mol C m\eqn{^{-2}} s\eqn{^{-1}}).} +#' \item{\code{jmax}}{Maximum rate of electron transport for RuBP regeneration +#' (in mol CO\eqn{_2} m\eqn{^{-2}} s\eqn{^{-1}}).} +#' \item{\code{vcmax25}}{Maximum rate of carboxylation (Vcmax), +#' normalised to 25\eqn{^o}C (in mol C m\eqn{^{-2}} s\eqn{^{-1}}).} +#' \item{\code{jmax25}}{Maximum rate of electron transport, normalised to +#' 25\eqn{^o}C (in mol C m\eqn{^{-2}} s\eqn{^{-1}}).} +#' \item{\code{gs_accl}}{Acclimated stomatal conductance (in +#' mol C (mol photons)\eqn{^{-1}} Pa\eqn{^{-1}}. (Multiply by +#' ppfd (mol photons m\eqn{^{-2}} d\eqn{^{-1}}) and fapar +#' to express per unit ground area and time.)} +#' \item{\code{wscal}}{Relative soil water content, between 0 (permanent wilting +#' point, PWP) and 1 (field capacity, FC).} +#' \item{\code{chi}}{Ratio of leaf-internal to ambient CO\eqn{_{2}}, ci:ca (unitless).} +#' \item{\code{iwue}}{Intrinsic water use efficiency (iWUE) (unitless, +#' multiply with patm (Pa) to get iWUE in Pa).} +#' \item{\code{rd}}{Dark respiration (Rd) in gC m\eqn{^{-2}} s\eqn{^{-1}}. +#' (Multiply by 1/12 (mol C / gC) to convert to mol C m\eqn{^{-2}} s\eqn{^{-1}}.)} +#' } +#' +#' @details TBC +#' +#' @export +#' @useDynLib rsofun +#' +#' @examples +#' # 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, +#' beta_unitcostratio = 146.0, +#' rd_to_vcmax = 0.014, # from Atkin et al. 2015 for C3 herbaceous +#' kc_jmax = 0.41 +#' ) +#' +#' # Run the Fortran P-model +#' run_pmodel_onestep_f_bysite( +#' lc4 = FALSE, +#' forcing = data.frame( +#' temp = 20, # temperature, deg C +#' vpd = 1000, # Pa, +#' ppfd = 300/10^6, # mol/m2/s +#' co2 = 400, # ppm, +#' patm = 101325 # Pa +#' ), +#' 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, +#' beta_unitcostratio = 146.0, +#' rd_to_vcmax = 0.014, # from Atkin et al. 2015 for C3 herbaceous +#' kc_jmax = 0.41 +#' ), +#' makecheck = TRUE +#' ) + +run_pmodel_onestep_f_bysite <- function( + lc4, + forcing, + params_modl, + makecheck = TRUE +){ + + # base state, always execute the call + continue <- TRUE + + # re-define units and naming of forcing dataframe + # keep the order of columns - it's critical for Fortran (reading by column number) + forcing_features <- c( + 'temp', + 'vpd', + 'ppfd', + 'co2', + 'patm' + ) + forcing <- forcing %>% + dplyr::select(all_of(forcing_features)) + + # validate input + if (makecheck){ + + is.nanull <- function(x) ifelse(any(is.null(x), is.na(x)), TRUE, FALSE) + + # list variable to check for + check_vars <- c( + "temp", + "vpd", + 'ppfd', + "co2", + "patm" + ) + + # create a loop to loop over a list of variables + # to check validity + data_integrity <- lapply(check_vars, function(check_var){ + if (any(is.nanull(forcing[check_var]))){ + sitename <- "undefined site" + warning(sprintf("Error: Missing value in %s for %s", + check_var, sitename)) + return(FALSE) + } else { + return(TRUE) + } + }) + + # only run simulation if all checked variables are valid + # suppress warning on coercion of list to single logical + if (suppressWarnings(!all(data_integrity))){ + continue <- FALSE + } + + + # model parameters to check + if ( sum( names(params_modl) %in% c('kphio', 'kphio_par_a', 'kphio_par_b', + 'beta_unitcostratio', 'rd_to_vcmax', 'kc_jmax') + ) != 6){ + warning(" Returning a dummy data frame. Incorrect model parameters.") + continue <- FALSE + } + } + + if (continue){ + ## C wrapper call + pmodelonestep_out <- .Call( + 'pmodel_onestep_f_C', + ## Simulation parameters + lc4 = as.logical(lc4), + par = c(as.numeric(params_modl$kphio), # model parameters as vector in order + as.numeric(params_modl$kphio_par_a), + as.numeric(params_modl$kphio_par_b), + as.numeric(params_modl$beta_unitcostratio), + as.numeric(params_modl$rd_to_vcmax), + as.numeric(params_modl$kc_jmax)), + forcing = as.matrix(forcing) + ) + pmodelonestep_out <- t(pmodelonestep_out) + + } else { + pmodelonestep_out <- array(dim = c(1,9)) + } + + out <- build_out_pmodel_onestep(pmodelonestep_out) + + return(out) +} + +# Build R output +build_out_pmodel_onestep <- function(pmodelonestep_out){ + out <- pmodelonestep_out %>% + as.matrix() %>% + as.data.frame() %>% + stats::setNames( + c("vcmax", + "jmax", + "vcmax25", + "jmax25", + "gs_accl", + "wscal", + "chi", + "iwue", + "rd") + ) %>% + as_tibble(.name_repair = "check_unique") + + return(out) +} + +.onUnload <- function(libpath) { + library.dynam.unload("rsofun", libpath) +} diff --git a/cran-comments.md b/cran-comments.md index 6561170d6..721bcb749 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,7 +1,11 @@ -## Resubmission +## R CMD check results -- Fix parallel make errors +0 errors ✔ | 0 warnings ✔ | 0 notes ✔ -## R CMD check +## revdepcheck results + +We checked 0 reverse dependencies, comparing R CMD check results across CRAN and dev versions of this package. + + * We saw 0 new problems + * We failed to check 0 packages -Status: Ok diff --git a/man/biomee_gs_leuning_drivers.Rd b/man/biomee_gs_leuning_drivers.Rd index f2e6f8d84..983916b30 100644 --- a/man/biomee_gs_leuning_drivers.Rd +++ b/man/biomee_gs_leuning_drivers.Rd @@ -35,13 +35,13 @@ This data structure can be freely used for documenting the dataset, but must inc }} \item{forcing}{Forcing data.frame used as input \describe{ - \item{ppfd}{Photosynthetic photon flux density (mol s-1 m-2)} + \item{ppfd}{Photosynthetic photon flux density (mol m-2 s-1)} \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{pair}{Atmospheric pressure (Pa)} + \item{co2}{Atmospheric CO\eqn{_2} concentration in ppm.} }} \item{params_tile}{Tile-level model parameters, into a single row data.frame, including the following data: \describe{ diff --git a/man/init_dates_dataframe.Rd b/man/init_dates_dataframe.Rd index 881b91128..ecb6f4e31 100644 --- a/man/init_dates_dataframe.Rd +++ b/man/init_dates_dataframe.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/init_dates_dataframe.R +% Please edit documentation in R/run_pmodel_f_bysite.R \name{init_dates_dataframe} \alias{init_dates_dataframe} \title{Initialises a tibble with dates} diff --git a/man/p_model_drivers.Rd b/man/p_model_drivers.Rd index 11fb28db1..0bd6bd3cc 100644 --- a/man/p_model_drivers.Rd +++ b/man/p_model_drivers.Rd @@ -26,7 +26,7 @@ A tibble of driver data: \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{co2}{Atmospheric CO\eqn{_2} concentration in ppm.} \item{ccov}{Cloud coverage in \%. This is only used when either PPFD or net radiation are not prescribed.} } diff --git a/man/run_pmodel_f_bysite.Rd b/man/run_pmodel_f_bysite.Rd index 57ae93ec5..62c172a36 100644 --- a/man/run_pmodel_f_bysite.Rd +++ b/man/run_pmodel_f_bysite.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/run_pmodel_f_bysite.R \name{run_pmodel_f_bysite} \alias{run_pmodel_f_bysite} -\title{Run P-model (R wrapper)} +\title{Run P-model (time series)} \usage{ run_pmodel_f_bysite( sitename, @@ -47,20 +47,24 @@ Model output is provided as a tidy dataframe, with columns: \item{\code{le}}{Latent heat flux (in J m\eqn{^{-2}} d\eqn{^{-1}}).} \item{\code{pet}}{Potential evapotranspiration (PET), calculated by SPLASH following Priestly-Taylor (in mm d\eqn{^{-1}}).} \item{\code{vcmax}}{Maximum rate of RuBisCO carboxylation - (Vcmax) (in mol C m\eqn{^{-2}} d\eqn{^{-1}}).} + (Vcmax) (in mol C m\eqn{^{-2}} s\eqn{^{-1}}).} \item{\code{jmax}}{Maximum rate of electron transport for RuBP regeneration (in mol CO\eqn{_2} m\eqn{^{-2}} s\eqn{^{-1}}).} \item{\code{vcmax25}}{Maximum rate of carboxylation (Vcmax), - normalised to 25\eqn{^o}C (in mol C m\eqn{^{-2}} d\eqn{^{-1}}).} + normalised to 25\eqn{^o}C (in mol C m\eqn{^{-2}} s\eqn{^{-1}}).} \item{\code{jmax25}}{Maximum rate of electron transport, normalised to 25\eqn{^o}C (in mol C m\eqn{^{-2}} s\eqn{^{-1}}).} \item{\code{gs_accl}}{Acclimated stomatal conductance (in - mol C m\eqn{^{-2}} d\eqn{^{-1}} Pa\eqn{^{-1}}).} + mol C (mol photons)\eqn{^{-1}} Pa\eqn{^{-1}}. (Multiply by + ppfd (mol photons m\eqn{^{-2}} d\eqn{^{-1}}) and fapar + to express per unit ground area and time.)} \item{\code{wscal}}{Relative soil water content, between 0 (permanent wilting point, PWP) and 1 (field capacity, FC).} \item{\code{chi}}{Ratio of leaf-internal to ambient CO\eqn{_{2}}, ci:ca (unitless).} - \item{\code{iwue}}{Intrinsic water use efficiency (iWUE) (in Pa).} - \item{\code{rd}}{Dark respiration (Rd) in gC m\eqn{^{-2}} d\eqn{^{-1}}.} + \item{\code{iwue}}{Intrinsic water use efficiency (iWUE) (unitless, + multiply with patm (Pa) to get iWUE in Pa).} + \item{\code{rd}}{Dark respiration (Rd) in gC m\eqn{^{-2}} s\eqn{^{-1}}. + (Multiply by 1/12 (mol C / gC) to convert to mol C m\eqn{^{-2}} s\eqn{^{-1}}.)} \item{\code{tsoil}}{Soil temperature, in \eqn{^{o}}C.} \item{\code{netrad}}{Net radiation, in W m\eqn{^{-2}}. WARNING: this is currently ignored as a model forcing. Instead, net radiation is internally calculated by SPLASH.} \item{\code{wcont}}{Soil water content, in mm.} @@ -69,7 +73,7 @@ Model output is provided as a tidy dataframe, with columns: } } \description{ -Run P-model Fortran model on single site. +Run P-model on a single site for a forcing time series. } \details{ Depending on the input model parameters, it's possible to run the diff --git a/man/run_pmodel_onestep_f_bysite.Rd b/man/run_pmodel_onestep_f_bysite.Rd new file mode 100644 index 000000000..802f8014a --- /dev/null +++ b/man/run_pmodel_onestep_f_bysite.Rd @@ -0,0 +1,84 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/run_pmodel_onestep_f_bysite.R +\name{run_pmodel_onestep_f_bysite} +\alias{run_pmodel_onestep_f_bysite} +\title{Run P-model (single time step)} +\usage{ +run_pmodel_onestep_f_bysite(lc4, forcing, params_modl, makecheck = TRUE) +} +\arguments{ +\item{lc4}{Locigical specifying whether P-model simulation is for C4 (as opposed to C3). Defaults to \code{FALSE}.} + +\item{forcing}{A data frame of forcing climate data, used as input (single row).} + +\item{params_modl}{A named list of free (calibratable) model parameters. See \code{\link{runread_pmodel_f}}} + +\item{makecheck}{A logical specifying whether checks are performed + to verify forcings and model parameters. \code{TRUE} by default. + +For further specifications of above inputs and examples see \code{\link{p_model_drivers}} or \code{\link{p_model_drivers_vcmax25}}} +} +\value{ +Model output is provided as a tidy dataframe, with columns: +\describe{ + \item{\code{vcmax}}{Maximum rate of RuBisCO carboxylation + (Vcmax) (in mol C m\eqn{^{-2}} s\eqn{^{-1}}).} + \item{\code{jmax}}{Maximum rate of electron transport for RuBP regeneration + (in mol CO\eqn{_2} m\eqn{^{-2}} s\eqn{^{-1}}).} + \item{\code{vcmax25}}{Maximum rate of carboxylation (Vcmax), + normalised to 25\eqn{^o}C (in mol C m\eqn{^{-2}} s\eqn{^{-1}}).} + \item{\code{jmax25}}{Maximum rate of electron transport, normalised to + 25\eqn{^o}C (in mol C m\eqn{^{-2}} s\eqn{^{-1}}).} + \item{\code{gs_accl}}{Acclimated stomatal conductance (in + mol C (mol photons)\eqn{^{-1}} Pa\eqn{^{-1}}. (Multiply by + ppfd (mol photons m\eqn{^{-2}} d\eqn{^{-1}}) and fapar + to express per unit ground area and time.)} + \item{\code{wscal}}{Relative soil water content, between 0 (permanent wilting + point, PWP) and 1 (field capacity, FC).} + \item{\code{chi}}{Ratio of leaf-internal to ambient CO\eqn{_{2}}, ci:ca (unitless).} + \item{\code{iwue}}{Intrinsic water use efficiency (iWUE) (unitless, + multiply with patm (Pa) to get iWUE in Pa).} + \item{\code{rd}}{Dark respiration (Rd) in gC m\eqn{^{-2}} s\eqn{^{-1}}. + (Multiply by 1/12 (mol C / gC) to convert to mol C m\eqn{^{-2}} s\eqn{^{-1}}.)} + } +} +\description{ +Run P-model on a single site for a single time step. This does not include the simulation of ecosystem-level quantities, +water limitation, nor a simulation of water fluxes. Instead, this corresponds to a leaf-level representation of the +acclimation of photosynthesis. +} +\details{ +TBC +} +\examples{ +# 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, + beta_unitcostratio = 146.0, + rd_to_vcmax = 0.014, # from Atkin et al. 2015 for C3 herbaceous + kc_jmax = 0.41 +) + +# Run the Fortran P-model +run_pmodel_onestep_f_bysite( + lc4 = FALSE, + forcing = data.frame( + temp = 20, # temperature, deg C + vpd = 1000, # Pa, + ppfd = 300/10^6, # mol/m2/s + co2 = 400, # ppm, + patm = 101325 # Pa + ), + 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, + beta_unitcostratio = 146.0, + rd_to_vcmax = 0.014, # from Atkin et al. 2015 for C3 herbaceous + kc_jmax = 0.41 + ), + makecheck = TRUE +) +} diff --git a/src/cohort_biomee.mod.f90 b/src/cohort_biomee.mod.f90 index c183f9b42..288fcda47 100644 --- a/src/cohort_biomee.mod.f90 +++ b/src/cohort_biomee.mod.f90 @@ -33,7 +33,7 @@ module md_cohort integer :: species = 1 ! vegetation species !===== Population structure - real :: density = 0.0 ! density of vegetation, tree/m2 + real :: density = 0.0 ! density of vegetation, tree m-2 real :: deathrate = 0.0 ! Deathrate (0 to 1) of the original cohort. Only set for cohorts in the killed list. real :: age = 0.0 ! age of cohort, years real :: topyear = 0.0 ! number of years the cohort is in top layer diff --git a/src/gpp_biomee.mod.f90 b/src/gpp_biomee.mod.f90 index ba536526e..709c4a1e8 100644 --- a/src/gpp_biomee.mod.f90 +++ b/src/gpp_biomee.mod.f90 @@ -30,7 +30,7 @@ module md_gpp_biomee real :: kc_jmax = 0.41 ! Jmax cost ratio ! these should be species-specific, temporary solution to put them here - real :: kphio = 0.05 ! quantum yield efficiency at optimal temperature, phi_0 (Stocker et al., 2020 GMD Eq. 10) + real :: kphio = 0.05 ! quantum yield efficiency at optimal temperature, phi_0 (Stocker et al., 2020 GMD Eq. 10 and Eq. 20) real :: kphio_par_a = 0.0 ! shape parameter of temperature-dependency of quantum yield efficiency real :: kphio_par_b = 25.0 ! optimal temperature of quantum yield efficiency (deg C) diff --git a/src/gpp_pmodel.mod.f90 b/src/gpp_pmodel.mod.f90 index daf11bdae..d8ca469c6 100644 --- a/src/gpp_pmodel.mod.f90 +++ b/src/gpp_pmodel.mod.f90 @@ -183,11 +183,11 @@ subroutine gpp( tile, tile_fluxes, co2, climate, grid, init, in_ppfd) ! but not too dangerous... !---------------------------------------------------------------- if( in_ppfd ) then - ! Take input daily PPFD (in mol/m^2) + ! Take input daily PPFD (dppfd in mol m-2 s-1) tile_fluxes(lu)%plant(pft)%dgpp = tile(lu)%plant(pft)%fpc_grid * tile(lu)%canopy%fapar & * climate%dppfd * myinterface%params_siml%secs_per_tstep * out_pmodel%lue * soilmstress else - ! Take daily PPFD generated by SPLASH (in mol/m^2/d) + ! Take daily PPFD generated by SPLASH (ppfd_splash in mol m-2 d-1) tile_fluxes(lu)%plant(pft)%dgpp = tile(lu)%plant(pft)%fpc_grid * tile(lu)%canopy%fapar & * tile_fluxes(lu)%canopy%ppfd_splash * out_pmodel%lue * soilmstress end if @@ -231,7 +231,7 @@ end subroutine gpp ! ! arguments ! real, intent(in) :: fapar ! fraction of absorbed photosynthetically active radiation (unitless) ! real, intent(in) :: fpc_grid ! foliar projective cover, used for dividing grid cell area (unitless) - ! real, intent(in) :: dppfd ! daily total photon flux density (mol m-2) + ! real, intent(in) :: dppfd ! daily total photon flux density (mol m-2 d-1) ! real, intent(in) :: lue ! light use efficiency (g CO2 mol-1) ! real, intent(in) :: kphio_temp ! air temperature (deg C) ! real, intent(in) :: soilmstress ! soil moisture stress factor (unitless) @@ -351,7 +351,7 @@ end subroutine gpp ! real, intent(in) :: meanmppfd ! ! function return variable - ! real :: my_vcmax ! canopy-level Vcmax [gCO2/m2-ground/s] + ! real :: my_vcmax ! canopy-level Vcmax (gCO2/m2-ground/s) ! ! Calculate leafy-scale Rubisco-N as a function of LAI and current LUE ! my_vcmax = fapar * meanmppfd * vcmax_unitiabs diff --git a/src/photosynth_pmodel.mod.f90 b/src/photosynth_pmodel.mod.f90 index 795637aa1..c2a5c2425 100644 --- a/src/photosynth_pmodel.mod.f90 +++ b/src/photosynth_pmodel.mod.f90 @@ -27,10 +27,10 @@ module md_photosynth real :: iwue ! intrinsic water use efficiency = A / gs = ca - ci = ca ( 1 - chi ) , unitless real :: lue ! light use efficiency (mol CO2 / mol photon) ! real :: assim ! leaf-level assimilation rate (mol CO2 m-2 s-1) - real :: gs_setpoint ! stomatal conductance to CO2 (mol C Pa-1 m-2 s-1) + real :: gs_setpoint ! stomatal conductance to CO2 (mol CO2 Pa-1 (mol photons)-1) ! real :: gs_unitfapar ! stomatal conductance to CO2 per unit fapar (mol C Pa-1 m-2 s-1) - ! real :: gs_unitiabs ! stomatal conductance to CO2 per unit absorbed light (mol C Pa-1 m-2 s-1) - ! real :: gpp ! gross primary productivity (g CO2 m-2 d-1) + ! real :: gs_unitiabs ! stomatal conductance to CO2 per unit absorbed light (mol CO2 Pa-1 (mol photons)-1) + ! real :: gpp ! gross primary productivity (g CO2 m-2 s-1) ! real :: vcmax ! canopy-level maximum carboxylation capacity per unit ground area (mol CO2 m-2 s-1) real :: jmax25 ! canopy-level maximum rate of electron transport, normalized to 25 deg C (mol m-2 s-1) real :: vcmax25 ! canopy-level Vcmax25 (Vcmax normalized to 25 deg C) (mol CO2 m-2 s-1) @@ -38,9 +38,9 @@ module md_photosynth ! real :: vcmax_unitiabs ! Vcmax per unit absorbed light (mol CO2 m-2 s-1 mol-1) ! real :: ftemp_inst_vcmax ! Instantaneous temperature response factor of Vcmax (unitless) ! real :: ftemp_inst_rd ! Instantaneous temperature response factor of Rd (unitless) - ! real :: rd ! Dark respiration (mol CO2 m-2 s-1) - ! real :: rd_unitfapar ! Dark respiration per unit fAPAR (mol CO2 m-2 s-1) - ! real :: rd_unitiabs ! Dark respiration per unit absorbed light (mol CO2 m-2 s-1) + ! real :: rd ! Dark respiration (g C m-2 s-1) + ! real :: rd_unitfapar ! Dark respiration per unit fAPAR (g C m-2 s-1) + ! real :: rd_unitiabs ! Dark respiration per unit absorbed light (g C (mol photons)-1) real :: actnv ! Canopy-level total metabolic leaf N per unit ground area (g N m-2) real :: actnv_unitfapar ! Metabolic leaf N per unit fAPAR (g N m-2) real :: actnv_unitiabs ! Metabolic leaf N per unit absorbed light (g N m-2 mol-1) @@ -78,9 +78,9 @@ function pmodel( kphio, beta, kc_jmax, ppfd, co2, tc, vpd, patm, c4, method_optc ! function calc_dgpp(). !------------------------------------------------------------------ ! arguments - real, intent(in) :: kphio ! apparent quantum yield efficiency - real, intent(in) :: beta ! parameter for the unit cost ratio (corresponding to beta in Prentice et al., 2014) - real, intent(in) :: kc_jmax ! parameter Jmax cost ratio (corresponding to c in Prentice et al., 2014) + real, intent(in) :: kphio ! apparent quantum yield efficiency (mol mol-1) + real, intent(in) :: beta ! parameter for the unit cost ratio (-) (corresponding to beta in Prentice et al., 2014) + real, intent(in) :: kc_jmax ! parameter Jmax cost ratio (-) (corresponding to c in Prentice et al., 2014) ! real, intent(in) :: fapar ! fraction of absorbed photosynthetically active radiation (unitless) real, intent(in) :: ppfd ! photosynthetic photon flux density (mol m-2 s-1), relevant for acclimated response real, intent(in) :: co2 ! atmospheric CO2 concentration (ppm), relevant for acclimated response @@ -99,9 +99,9 @@ function pmodel( kphio, beta, kc_jmax, ppfd, co2, tc, vpd, patm, c4, method_optc real :: kmm ! Michaelis-Menten coefficient (Pa) real :: gammastar ! photorespiratory compensation point - Gamma-star (Pa) real :: ca ! ambient CO2 partial pressure, (Pa) - real :: gs_setpoint ! stomatal conductance to CO2 (mol CO2 Pa-1 m-2 s-1) + real :: gs_setpoint ! stomatal conductance to CO2 (mol CO2 Pa-1 (mol photons)-1) ! real :: gs_unitfapar ! stomatal conductance to CO2 (mol CO2 Pa-1 m-2 s-1) - ! real :: gs_unitiabs ! stomatal conductance to CO2 (mol CO2 Pa-1 m-2 s-1) + ! real :: gs_unitiabs ! stomatal conductance to CO2 (mol CO2 Pa-1 (mol photons)-1) real :: ci ! leaf-internal partial pressure, (Pa) real :: chi ! = ci/ca, leaf-internal to ambient CO2 partial pressure, ci/ca (unitless) real :: ns ! viscosity of H2O at ambient temperatures (Pa s) @@ -110,7 +110,7 @@ function pmodel( kphio, beta, kc_jmax, ppfd, co2, tc, vpd, patm, c4, method_optc real :: mprime ! factor in light use model with Jmax limitation real :: iwue ! intrinsic water use efficiency = A / gs = ca - ci = ca ( 1 - chi ) , unitless real :: lue ! light use efficiency (mol CO2 / mol photon) - ! real :: gpp ! gross primary productivity (g CO2 m-2 d-1) + ! real :: gpp ! gross primary productivity (g CO2 m-2 s-1) real :: jmax ! canopy-level maximum rate of electron transport (XXX) real :: jmax25 ! canopy-level maximum rate of electron transport (XXX) real :: vcmax ! canopy-level maximum carboxylation capacity per unit ground area (mol CO2 m-2 s-1) @@ -122,9 +122,9 @@ function pmodel( kphio, beta, kc_jmax, ppfd, co2, tc, vpd, patm, c4, method_optc ! real :: vcmax25_unitiabs ! Vcmax25 per unit absorbed light (mol CO2 m-2 s-1 mol-1) real :: ftemp_inst_vcmax ! Instantaneous temperature response factor of Vcmax (unitless) real :: ftemp_inst_jmax ! Instantaneous temperature response factor of Jmax (unitless) - ! real :: rd ! Dark respiration (mol CO2 m-2 s-1) - ! real :: rd_unitfapar ! Dark respiration per unit fAPAR (mol CO2 m-2 s-1) - ! real :: rd_unitiabs ! Dark respiration per unit absorbed light (mol CO2 m-2 s-1) + ! real :: rd ! Dark respiration (g C m-2 s-1) + ! real :: rd_unitfapar ! Dark respiration per unit fAPAR (g C m-2 s-1) + ! real :: rd_unitiabs ! Dark respiration per unit absorbed light (g C (mol photons)-1) real :: actnv ! Canopy-level total metabolic leaf N per unit ground area (g N m-2) ! real :: actnv_unitfapar ! Metabolic leaf N per unit fAPAR (g N m-2) ! real :: actnv_unitiabs ! Metabolic leaf N per unit absorbed light (g N m-2 mol-1) @@ -239,7 +239,8 @@ function pmodel( kphio, beta, kc_jmax, ppfd, co2, tc, vpd, patm, c4, method_optc lue = kphio * mprime * c_molmass ! in g CO2 m-2 s-1 / (mol light m-2 s-1) ! Vcmax after accounting for Jmax limitation - vcmax = kphio * ppfd * out_optchi%mjoc * mprime / out_optchi%mj + vcmax = kphio * ppfd * out_optchi%mjoc * mprime / out_optchi%mj ! mol m-2 s-1 + ! (-) * (mol m-2 s-1) * (-) * (-) / (-) ! xxx test ! print*,'out_optchi%mjoc : ', out_optchi%mjoc @@ -327,7 +328,8 @@ function pmodel( kphio, beta, kc_jmax, ppfd, co2, tc, vpd, patm, c4, method_optc if (fact_jmaxlim >= 1 .or. fact_jmaxlim <= 0) then jmax = dummy else - jmax = 4.0 * kphio * ppfd / sqrt( (1.0/fact_jmaxlim)**2 - 1.0 ) + jmax = 4.0 * kphio * ppfd / sqrt( (1.0/fact_jmaxlim)**2 - 1.0 ) ! mol m-2 s-1 + ! (-) * (-) * (mol m-2 s-1) / (-) end if ! for normalization using temperature response from Duursma et al., 2015, implemented in plantecophys R package ftemp_inst_jmax = calc_ftemp_inst_jmax( tc, tc, tcref = 25.0 ) @@ -838,9 +840,9 @@ function calc_kphio_temp( dtemp, c4, kphio, kphio_par_a, kphio_par_b ) result( k ! arguments real, intent(in) :: dtemp ! (leaf) temperature in degrees celsius logical, intent(in) :: c4 - real, intent(in) :: kphio - real, intent(in) :: kphio_par_a - real, intent(in) :: kphio_par_b + real, intent(in) :: kphio ! (mol mol-1) + real, intent(in) :: kphio_par_a ! ((deg C)-2) + real, intent(in) :: kphio_par_b ! (deg C) ! function return variable real :: kphio_temp diff --git a/src/plant_pmodel.mod.f90 b/src/plant_pmodel.mod.f90 index b8136cde7..6d937c851 100644 --- a/src/plant_pmodel.mod.f90 +++ b/src/plant_pmodel.mod.f90 @@ -59,8 +59,8 @@ module md_plant_pmodel real :: lma ! leaf mass per area (gC m-2) real :: sla ! specific leaf area (m2 gC-1) real :: nmass ! leaf N per unit leaf mass, g N / g-dry mass - real :: r_cton_leaf ! leaf C:N ratio [gC/gN] - real :: r_ntoc_leaf ! leaf N:C ratio [gN/gC] + real :: r_cton_leaf ! leaf C:N ratio (gC/gN) + real :: r_ntoc_leaf ! leaf N:C ratio (gN/gC) end type plant_type @@ -71,22 +71,22 @@ module md_plant_pmodel type plant_fluxes_type ! daily updated variables - real :: dgpp ! daily gross primary production [gC/m2/d] - real :: drd ! daily dark respiration [gC/m2/d] + real :: dgpp ! daily gross primary production (gC m-2 d-1) + real :: drd ! daily dark respiration (gC m-2 d-1) real :: assim ! daily assimilation (mol CO2 m-2 s-1) - real :: dtransp ! daily transpiration [mm] - real :: dlatenth ! daily latent heat flux [J m-2 d-1] + real :: dtransp ! daily transpiration (mm) + real :: dlatenth ! daily latent heat flux (J m-2 d-1) real :: vcmax25 ! acclimated Vcmax, normalised to 25 deg C (mol CO2 m-2 s-1) real :: jmax25 ! acclimated Jmax, normalised to 25 deg C (mol CO2 m-2 s-1) real :: vcmax ! daily varying Vcmax (mol CO2 m-2 s-1) real :: jmax ! daily varying Jmax (mol CO2 m-2 s-1) - real :: gs_accl ! acclimated stomatal conductance (xxx) + real :: gs_accl ! acclimated stomatal conductance (mol CO2 Pa-1 (mol photons)-1) real :: chi ! ci:ca ratio (unitless) real :: iwue ! intrinsic water use efficiency (A/gs = ca*(1-chi)) ! ! annual variables - ! real :: agpp ! annual total gross primary production [gC/m2/yr] + ! real :: agpp ! annual total gross primary production (gC m-2 yr-1) ! real :: avcmax25_mean ! annual Vcmax, normalised to 25 deg C, GPP-weighted mean ! real :: avcmax25_max ! annual Vcmax, normalised to 25 deg C, annual maximum diff --git a/src/pmodel.mod.f90 b/src/pmodel.mod.f90 index 3869fb2d6..543f2d746 100644 --- a/src/pmodel.mod.f90 +++ b/src/pmodel.mod.f90 @@ -7,9 +7,116 @@ module pmodel_mod implicit none private - public :: pmodel_f + public :: pmodel_f, pmodel_onestep_f contains + subroutine pmodel_onestep_f( & + lc4, & + par, & + forcing, & + output & + ) bind(C, name = "pmodel_onestep_f_") + + !//////////////////////////////////////////////////////////////// + ! Main subroutine to handle I/O with C and R. + ! Receives simulation parameters, site parameters. Forcing is for + ! one step. + !---------------------------------------------------------------- + use md_params_core, only: eps, c_molmass + use md_photosynth + + implicit none + + ! arguments + integer(kind=c_int), intent(in) :: lc4 ! logical type is not supported in the C interface (LTO) + real(kind=c_double), dimension(6), intent(in) :: par ! free (calibratable) model parameters + real(kind=c_double), dimension(1,5), intent(in) :: forcing ! 1-row array containing forcing data (rows: time steps; columns: 1=air temperature, 2=vpd, 3=ppfd, 4=co2, 5=patm) + real(kind=c_double), dimension(9), intent(out) :: output + + ! local variables + real :: kphio, kphio_par_a, kphio_par_b, beta_unitcostratio, rd_to_vcmax, kc_jmax, & + temp, vpd, ppfd, co2, patm, kphio_temp, vcmax, jmax, rd + logical :: c4 + type(outtype_pmodel) :: out_pmodel ! list of P-model output variables + + ! convert integer to logcial + c4 = lc4 /= 0 + + !---------------------------------------------------------------- + ! GET CALIBRATABLE MODEL PARAMETERS (so far a small list) + !---------------------------------------------------------------- + kphio = real(par(1)) ! (mol mol-1) + kphio_par_a = real(par(2)) ! ((deg C)-2) + kphio_par_b = real(par(3)) ! (deg C) + beta_unitcostratio = real(par(4)) ! (-) + rd_to_vcmax = real(par(5)) ! (-/-), ratio of rd to Vcmax + kc_jmax = real(par(6)) ! (-), jmax cost ratio, a.k.a c* + + !---------------------------------------------------------------- + ! GET FORCING + !---------------------------------------------------------------- + temp = real(forcing(1,1)) ! (deg C) + vpd = real(forcing(1,2)) ! (Pa) + ppfd = real(forcing(1,3)) ! (mol m-2 s-1) (while netrad is W m-2) + co2 = real(forcing(1,4)) ! (ppm) + patm = real(forcing(1,5)) ! (Pa) + + !---------------------------------------------------------------- + ! Low-temperature effect on quantum yield efficiency + !---------------------------------------------------------------- + ! take the instananeously varying temperature for governing quantum yield variations + if (abs(kphio_par_a) < eps) then + + kphio_temp = kphio + + else + + kphio_temp = calc_kphio_temp( & + temp, & + c4, & + kphio, & + kphio_par_a, & + kphio_par_b & + ) + + end if + + !================================================================ + ! P-model call to get acclimated quantities as a function of the + ! damped climate forcing. + !---------------------------------------------------------------- + out_pmodel = pmodel( & + kphio = kphio_temp, & + beta = beta_unitcostratio, & + kc_jmax = kc_jmax, & + ppfd = ppfd, & + co2 = co2, & + tc = temp, & + vpd = vpd, & + patm = patm, & + c4 = c4, & + method_optci = "prentice14", & + method_jmaxlim = "wang17" & + ) + + ! quantities with instantaneous temperature response + vcmax = calc_ftemp_inst_vcmax( temp, temp, tcref = 25.0 ) * out_pmodel%vcmax25 + jmax = calc_ftemp_inst_jmax( temp, temp, tcref = 25.0 ) * out_pmodel%jmax25 + rd = out_pmodel%vcmax25 * rd_to_vcmax * calc_ftemp_inst_rd( temp ) * c_molmass + + !---------------------------------------------------------------- + ! Populate Fortran output array which is passed back to C/R + !---------------------------------------------------------------- + output(1) = dble(vcmax) + output(2) = dble(jmax) + output(3) = dble(out_pmodel%vcmax25) + output(4) = dble(out_pmodel%jmax25) + output(5) = dble(out_pmodel%gs_setpoint) + output(7) = dble(out_pmodel%chi) + output(8) = dble(out_pmodel%iwue) + output(9) = dble(rd) + + end subroutine pmodel_onestep_f subroutine pmodel_f( & spinup, & @@ -124,16 +231,16 @@ subroutine pmodel_f( & !---------------------------------------------------------------- ! GET CALIBRATABLE MODEL PARAMETERS (so far a small list) !---------------------------------------------------------------- - myinterface%params_calib%kphio = real(par(1)) - myinterface%params_calib%kphio_par_a = real(par(2)) - myinterface%params_calib%kphio_par_b = real(par(3)) + myinterface%params_calib%kphio = real(par(1)) ! (mol mol-1) + myinterface%params_calib%kphio_par_a = real(par(2)) ! ((deg C)-2) + myinterface%params_calib%kphio_par_b = real(par(3)) ! (deg C) myinterface%params_calib%soilm_thetastar = real(par(4)) myinterface%params_calib%soilm_betao = real(par(5)) - myinterface%params_calib%beta_unitcostratio = real(par(6)) - myinterface%params_calib%rd_to_vcmax = real(par(7)) + myinterface%params_calib%beta_unitcostratio = real(par(6)) ! (-) + myinterface%params_calib%rd_to_vcmax = real(par(7)) ! (-/-), ratio of rd to Vcmax myinterface%params_calib%tau_acclim = real(par(8)) - myinterface%params_calib%kc_jmax = real(par(9)) - + myinterface%params_calib%kc_jmax = real(par(9)) ! (-), jmax cost ratio, a.k.a c* + !---------------------------------------------------------------- ! GET VEGETATION COVER (fractional projective cover by PFT) !---------------------------------------------------------------- diff --git a/src/soil_biomee.mod.f90 b/src/soil_biomee.mod.f90 index 049cc526e..e95d3f4fc 100644 --- a/src/soil_biomee.mod.f90 +++ b/src/soil_biomee.mod.f90 @@ -35,7 +35,7 @@ subroutine water_supply_layer( vegn ) type(cohort_type),pointer :: cc real :: fWup(MAX_LEVELS) ! fraction to the actual soil water real :: freewater(MAX_LEVELS) - real :: totWsup(MAX_LEVELS) ! potential water uptake, mol s-1 m-2 + real :: totWsup(MAX_LEVELS) ! potential water uptake, mol m-2 s-1 real :: thetaS(MAX_LEVELS) ! soil moisture index (0~1) real :: dpsiSR(MAX_LEVELS) ! pressure difference between soil water and root water, Pa integer :: i diff --git a/src/tile_pmodel.mod.f90 b/src/tile_pmodel.mod.f90 index ce7359e78..df8f910f6 100644 --- a/src/tile_pmodel.mod.f90 +++ b/src/tile_pmodel.mod.f90 @@ -15,10 +15,10 @@ module md_tile_pmodel ! physical soil state variables with memory from year to year (~pools) !---------------------------------------------------------------- type psoilphystype - real :: temp ! soil temperature [deg C] - real :: wcont ! liquid soil water mass [mm = kg/m2] + real :: temp ! soil temperature (deg C) + real :: wcont ! liquid soil water mass (mm = kg m-2) real :: wscal ! relative soil water content, between 0 (PWP) and 1 (FC) - real :: snow ! snow depth in liquid-water-equivalents [mm = kg/m2] + real :: snow ! snow depth in liquid-water-equivalents (mm = kg m-2) ! real :: rlmalpha ! rolling mean of annual mean alpha (AET/PET) end type psoilphystype @@ -88,8 +88,8 @@ module md_tile_pmodel real :: dro ! daily runoff (mm d-1) real :: dfleach ! daily fraction of soil water going to runoff (used for calculating leaching) real :: dwbal ! daily water balance as precipitation and snow melt minus runoff and evapotranspiration (mm d-1) - real :: econ ! water-to-energy conversion factor (m^3/J) - real :: drn ! daily total net radiation (J/m2/d) + real :: econ ! water-to-energy conversion factor (m3 J-1) + real :: drn ! daily total net radiation (J m-2 d-1) real :: drnn ! nighttime total net radiation (J m-1 d-1) real :: rnl ! net longwave radiation (W m-2) real :: dcn ! daily total condensation (mm d-1) @@ -116,21 +116,21 @@ module md_tile_pmodel ! real :: rho_water ! density of water (g m-3) ! carbon - real :: dgpp - real :: drd - real :: assim ! leaf-level assimilation rate + real :: dgpp ! daily gross primary production (gC m-2 d-1) + real :: drd ! daily dark respiration (gC m-2 d-1) + real :: assim ! leaf-level assimilation rate (mol CO2 m-2 s-1) real :: vcmax25 ! acclimated Vcmax, normalised to 25 deg C (mol CO2 m-2 s-1) real :: jmax25 ! acclimated Jmax, normalised to 25 deg C (mol CO2 m-2 s-1) real :: vcmax ! daily varying Vcmax (mol CO2 m-2 s-1) real :: jmax ! daily varying Jmax (mol CO2 m-2 s-1) - real :: gs_accl ! acclimated stomatal conductance (xxx) + real :: gs_accl ! acclimated stomatal conductance (mol CO2 Pa-1 (mol photons)-1) real :: chi ! ci:ca ratio (unitless) real :: iwue ! intrinsic water use efficiency (A/gs = ca*(1-chi)) ! radiation - real :: ppfd_splash - real :: dra ! daily top-of-atmosphere solar radiation (J/m^2/d) + real :: ppfd_splash ! daily photosynthetic photon flux density (mol m-2 d-1, Note that this differs from ppfd input 'dppfd', which is mol m-2 s-1.) + real :: dra ! daily top-of-atmosphere solar radiation (J m-2 d-1) ! ! annual ! !---------------------------------------------------------------- @@ -139,7 +139,7 @@ module md_tile_pmodel ! real :: avcmax25_mean ! annual Vcmax, normalised to 25 deg C, GPP-weighted mean ! real :: avcmax25_max ! annual Vcmax, normalised to 25 deg C, annual maximum - ! real, dimension(ndayyear) :: dra ! daily TOA solar irradiation (J/m2) + ! real, dimension(ndayyear) :: dra ! daily TOA solar irradiation (J m-2) ! real, dimension(ndayyear) :: dppfd_splash ! daily total PPFD (mol m-2 d-1) ! real, dimension(nmonth) :: mppfd_splash ! monthly total PPFD (mol m-2 month-1) ! real, dimension(nmonth) :: meanmppfd_splash ! monthly mean PPFD, averaged over daylight seconds (mol m-2 s-1) diff --git a/src/vegetation_processes_biomee.mod.f90 b/src/vegetation_processes_biomee.mod.f90 index aa2da8016..79b3f602f 100755 --- a/src/vegetation_processes_biomee.mod.f90 +++ b/src/vegetation_processes_biomee.mod.f90 @@ -833,8 +833,8 @@ subroutine update_plant_pools( cc, vegn, dL, dR, dStem) type(orgpool), intent(in) :: dL, dR, dStem ! leaf and fine root pool tendencies ! local variables - type(orgpool) :: loss_coarse, loss_fine, dtot, dAleaf_pool - real :: dAleaf ! leaf area decrease due to dBL + type(orgpool) :: loss_coarse, loss_fine, dtot !, dAleaf_pool + ! real :: dAleaf ! leaf area decrease due to dBL associate ( sp => cc%sp() ) ! Retranslocation to NSC and NSN diff --git a/src/waterbal_splash.mod.f90 b/src/waterbal_splash.mod.f90 index a83125081..1851d8a76 100644 --- a/src/waterbal_splash.mod.f90 +++ b/src/waterbal_splash.mod.f90 @@ -141,8 +141,8 @@ subroutine solar( tile_fluxes, grid, climate, doy ) ! 'evap', adopted from the evap() function in GePiSaT (Python version). ! This subroutine ('get_solar') is called before the daily loop. ! Output: - ! - daily extraterrestrial solar radiation (dra), J/m^2 - ! - daily PPFD (dppfd), mol/m^2 + ! - daily extraterrestrial solar radiation (dra), J m-2 d-1 + ! - daily PPFD (ppfd_splash), mol m-2 d-1 !------------------------------------------------------------------------- ! arguments type(tile_fluxes_type), dimension(nlu), intent(inout) :: tile_fluxes @@ -190,7 +190,7 @@ subroutine solar( tile_fluxes, grid, climate, doy ) grid%dayl = 24.0 * 60 * 60 * hs / 180.0 ! hs is in degrees (pi = 180 deg) !--------------------------------------------------------- - ! 7. Calculate daily extraterrestrial solar radiation (dra), J/m^2/d + ! 7. Calculate daily extraterrestrial solar radiation (dra), J m-2 d-1 !--------------------------------------------------------- ! Eq. 1.10.3, Duffy & Beckman (1993) tile_fluxes(:)%canopy%dra = ( secs_per_day / pi ) * kGsc * dr * ( radians(ru) * hs + rv * dgsin(hs) ) @@ -201,7 +201,7 @@ subroutine solar( tile_fluxes, grid, climate, doy ) tau = calc_tau( climate%dfsun, grid%elv ) !--------------------------------------------------------- - ! 9. Calculate daily PPFD (dppfd), mol/m^2 + ! 9. Calculate daily PPFD (ppfd_splash), mol m-2 d-1 !--------------------------------------------------------- ! Eq. 57, SPLASH 2.0 Documentation tile_fluxes(:)%canopy%ppfd_splash = (1.0e-6) * kfFEC * ( 1.0 - kalb_vis ) * tau * tile_fluxes(:)%canopy%dra diff --git a/src/wrappersc.c b/src/wrappersc.c index b2698d945..ea2fdd2cd 100644 --- a/src/wrappersc.c +++ b/src/wrappersc.c @@ -8,7 +8,7 @@ // See https://www.stat.berkeley.edu/~spector/s243/calling.pdf for more details ///////////////////////////////////////////////////////////// -// P-model +// P-model for time series ///////////////////////////////////////////////////////////// void F77_NAME(pmodel_f)( @@ -40,14 +40,14 @@ void F77_NAME(pmodel_f)( // C wrapper function for P-model extern SEXP pmodel_f_C( + SEXP secs_per_tstep, + SEXP in_ppfd, + SEXP in_netrad, SEXP spinup, SEXP spinupyears, SEXP recycle, SEXP firstyeartrend, SEXP nyeartrend, - SEXP secs_per_tstep, - SEXP in_ppfd, - SEXP in_netrad, SEXP outdt, SEXP ltre, SEXP ltne, @@ -105,11 +105,53 @@ extern SEXP pmodel_f_C( return output; } + +///////////////////////////////////////////////////////////// +// P-model for one step +///////////////////////////////////////////////////////////// + +void F77_NAME(pmodel_onestep_f)( + int *lc4, // LOGICAL + double *par, + double *forcing, + double *output + ); + +// C wrapper function for P-model +extern SEXP pmodel_onestep_f_C( + SEXP lc4, + SEXP par, + SEXP forcing + ){ + + // Specify output + // 2nd agument to allocVector is number of rows, 3rd is number of columns + SEXP output = PROTECT( allocVector(REALSXP, 9) ); + + // Fortran subroutine call + F77_CALL(pmodel_onestep_f)( + INTEGER(lc4), + REAL(par), + REAL(forcing), + REAL(output) + ); + + // // Output as list + // SEXP out_full = PROTECT( allocVector(VECSXP, 1) ); + // SET_VECTOR_ELT(out_full, 0, output); + + UNPROTECT(1); + + return output; +} + + + ///////////////////////////////////////////////////////////// // biomee ///////////////////////////////////////////////////////////// void F77_NAME(biomee_f)( - double *params_sim, + double *params_siml, double *site_info, double *params_tile, int *n_params_species, @@ -228,6 +270,7 @@ extern SEXP biomee_f_C( ///////////////////////////////////////////////////////////// static const R_CallMethodDef CallEntries[] = { {"pmodel_f_C", (DL_FUNC) &pmodel_f_C, 23}, // Specify number of arguments to C wrapper as the last number here + {"pmodel_onestep_f_C", (DL_FUNC) &pmodel_onestep_f_C, 3}, // Specify number of arguments to C wrapper as the last number here {"biomee_f_C", (DL_FUNC) &biomee_f_C, 12}, // Number of arguments of the C wrapper function for biomee (the SEXP variables, not the output) { NULL, NULL, 0 } }; @@ -238,5 +281,6 @@ void R_init_rsofun(DllInfo *dll) R_useDynamicSymbols(dll, FALSE); R_RegisterCCallable("rsofun", "pmodel_f_C", (DL_FUNC) &pmodel_f_C); + R_RegisterCCallable("rsofun", "pmodel_onestep_f_C", (DL_FUNC) &pmodel_onestep_f_C); R_RegisterCCallable("rsofun", "biomee_f_C", (DL_FUNC) &biomee_f_C); } diff --git a/tests/testthat/test-model-runs.R b/tests/testthat/test-model-runs.R index 351325b7d..bc631f59b 100644 --- a/tests/testthat/test-model-runs.R +++ b/tests/testthat/test-model-runs.R @@ -195,6 +195,97 @@ test_that("p-model run check Vcmax25", { expect_type(df_output_p, "list") }) +test_that("p-model onestep output check (run_pmodel_onestep_f_bysite())", { + skip_on_cran() + + # Define simulation inputs + inputs <- list( + # for forcing: + temp = 20, # temperature, deg C + vpd = 1000, # Pa, + ppfd = 300/10^6, # mol/m2/s + co2 = 400, # ppm, + patm = 101325, # Pa + fapar = 1, # - + # for params_modl + 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, + beta_unitcostratio = 146.0, + rd_to_vcmax = 0.014, # from Atkin et al. 2015 for C3 herbaceous + kc_jmax = 0.41 + ) + + # compute reference value: + library(rpmodel) + resR <- rpmodel( + tc=inputs$temp, vpd=inputs$vpd, co2=inputs$co2, + patm=inputs$patm, kphio=inputs$kphio, beta=inputs$beta_unitcostratio, + ppfd=inputs$ppfd, # rpmodel docs state that units of ppfd define output units of: lue,gpp,vcmax,rd + # this also affects: vcmax25,gs + # NOTE: unused inputs: kphio_par_a, kphio_par_b, rd_to_vcmax, kc_jmax + fapar = inputs$fapar, # fraction , + c4 = FALSE, + method_jmaxlim = "wang17", + do_ftemp_kphio = FALSE, # corresponding to setup ORG + do_soilmstress = FALSE, # corresponding to setup ORG + verbose = TRUE + ) |> + tidyr::as_tibble() + + # compute value with run_pmodel_onestep_f_bysite: + resF <- run_pmodel_onestep_f_bysite( + lc4 = FALSE, + forcing = data.frame(temp = inputs$temp, vpd = inputs$vpd, ppfd = inputs$ppfd, + co2 = inputs$co2, patm = inputs$patm), + params_modl = list( + kphio = inputs$kphio, + kphio_par_a = inputs$kphio_par_a, + kphio_par_b = inputs$kphio_par_b, + beta_unitcostratio = inputs$beta_unitcostratio, + rd_to_vcmax = inputs$rd_to_vcmax, + kc_jmax = inputs$kc_jmax + ), + makecheck = TRUE + ) + # testthat::expect_equal(resF, resR) + # NOTE: this fails because of different units. + + # Fix units for comparison + resR_units_fixed <- resR |> + dplyr::mutate( + gs = gs, # NOTE: keep units as-is: (mol C m-2 Pa-1) (computed as A/(ca-ci)) + iwue = iwue / inputs$patm, # NOTE: was initially: (Pa) (computed as (ca-ci)/1.6) + # NOTE: is now: (unitless) + rd = rd * 12.0107) |> # NOTE: was initially: 2.82e-7 mol C m-2 s-1 (computed as 0.015*Vcmax*(fr/fv)) + # NOTE: is now: (3.39e-6 g C m-2 s-1) # still slightly different from 3.16e-6, but remains within tolerance + dplyr::select(-ns_star, -xi, -mj, -mc, -ci, + -gpp, -ca, -gammastar, -kmm) |> + dplyr::select(vcmax, jmax, vcmax25, jmax25, chi, gs, iwue, rd) + + resF_units_fixed <- resF |> + dplyr::mutate( + gs = gs_accl * inputs$ppfd*inputs$fapar, + # NOTE: was initially: mol C / mol photons Pa-1 (computed as lue/molmass / (ca-ci+0.1)) + # NOTE: is now: mol C m-2 Pa-1 s-1 (computed as lue/molmass*iabs / (ca-ci+0.1)) + # NOTE: is now:gs with 4.75e-7 still slightly different from 4.79e-7, but remains within tolerance + iwue = iwue, # NOTE: keep units as-is: (-) (computed as (ca-ci)/1.6/patm) + rd = rd) |> # NOTE: keep units as-is: 3.16e-6 g C m-2 s-1 (computed as rd_to_vcmax*vcmax25*calc_ftemp_inst_rd(Temp)*c_molmass + dplyr::select(-wscal, -gs_accl) |> + dplyr::select(vcmax, jmax, vcmax25, jmax25, chi, gs, iwue, rd) + + + # Now comparison must pass + testthat::expect_equal(resR_units_fixed, resF_units_fixed, tolerance = 1e-5) + # If this test fails it means that the output of the model is incompatible with + # the package {rpmodel}. + # It could either mean that: + # - the model was accidentally altered and should be fixed to deliver the expected output + # - the package {rpmodel} has been altered + +}) + + test_that("biomeE output check (gs leuning)", { skip_on_cran() diff --git a/vignettes/files/sensitivity_analysis.Rmd__par_calib.RDS b/vignettes/files/sensitivity_analysis.Rmd__par_calib.RDS index 59bb325e8..743eea74e 100644 Binary files a/vignettes/files/sensitivity_analysis.Rmd__par_calib.RDS and b/vignettes/files/sensitivity_analysis.Rmd__par_calib.RDS differ diff --git a/vignettes/pmodel_use.Rmd b/vignettes/pmodel_use.Rmd index 27d7a6be6..861fcdc5f 100644 --- a/vignettes/pmodel_use.Rmd +++ b/vignettes/pmodel_use.Rmd @@ -237,3 +237,68 @@ ggplot(data = df_gpp_plot) + ``` For details on the optimization settings we refer to the manuals of [GenSA](https://cran.r-project.org/package=GenSA) and [BayesianTools](https://github.com/florianhartig/BayesianTools). + +## One-step P-model call + +For certain applications, the core theory of the P-model, predicting the acclimation of photosynthesis at the leaf-level, may be all that's needed for simulations. {rsofun} also provides a stripped-down version of "one-step" P-model call. That is, without considering forcing (and output) time series, delayed acclimation, ecosystem-level up-scaled quantities, and the coupling with the water balance and related water limitation on photosynthesis. Applications of the stripped-down version may include simulations of leaf-level photosynthetic traits. The function `run_pmodel_onestep_f_bysite()` provides this functionality and further improves computational efficiency over a call to `run_pmodel_f_bysite()`. This is particularly useful when calibrating the model to photosynthetic traits data. + +```{r} +run_pmodel_onestep_f_bysite( + lc4 = FALSE, + forcing = data.frame( + temp = 20, # temperature, deg C + vpd = 1000, # Pa, + ppfd = 300/10^6, # mol/m2/s + co2 = 400, # ppm, + patm = 101325 # Pa + ), + 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, + beta_unitcostratio = 146.0, + rd_to_vcmax = 0.014, # from Atkin et al. 2015 for C3 herbaceous + kc_jmax = 0.41 + ), + makecheck = TRUE +) +# +# +# make sure to check ?run_pmodel_onestep_f_bysite for units of output quantities. +``` + +`run_pmodel_onestep_f_bysite()` uses the same low-level code module as `run_pmodel_f_bysite()` (`photosynth_pmodel.mod.f90`) and therefore provides fully consistent outputs with the latter. Using `run_pmodel_onestep_f_bysite()` instead of `rpmodel::rpmodel()` guarantees consistency of the implementation across time-series and one-step model applications within the same package. Note that unit transformations as shown below is needed to obtain identical outputs. + +```{r} +library(rpmodel) +out_rpmodel <- rpmodel( + tc = 20, # temperature, deg C + vpd = 1000, # Pa, + co2 = 400, # ppm, + fapar = 1, # fraction, + ppfd = 300/10^6, # ~~mol/m2/d~~ or mol/m2/s + # rpmodel docs state that units of ppfd define + # output units of: lue, gpp, vcmax, rd + # (as well as vcmax25, gs) + patm = 101325, # Pa + kphio = 0.04998, # quantum yield efficiency as calibrated + # for setup ORG by Stocker et al. 2020 GMD, + beta = 146.0, # unit cost ratio a/b, + c4 = FALSE, + method_jmaxlim = "wang17", + do_ftemp_kphio = FALSE, # corresponding to setup ORG + do_soilmstress = FALSE, # corresponding to setup ORG + verbose = TRUE + ) + +# out_rpmodel +tidyr::as_tibble(out_rpmodel) |> + dplyr::select(vcmax, jmax, vcmax25, jmax25, gs, chi, iwue, rd) |> + # bring units to same as output of run_pmodel_onestep_f_bysite() + dplyr::mutate(gs = gs/1/(300/10^6), # divide by fapar (-) and by ppfd (mol/m2/s) + iwue = iwue/101325, # divide by patm (Pa) + rd = rd*12 # mutliply by c_molmass (gC/mol) + ) +``` + + diff --git a/vignettes/sensitivity_analysis.Rmd b/vignettes/sensitivity_analysis.Rmd index 00ee187d2..4bf0793b0 100644 --- a/vignettes/sensitivity_analysis.Rmd +++ b/vignettes/sensitivity_analysis.Rmd @@ -282,7 +282,7 @@ par_calib <- calib_sofun( par_fixed = par_fixed, targets = "gpp" ) -saveRDS(par_calib, "files/sensitivity_analysis.Rmd__par_calib.RDS") +saveRDS(par_calib, "files/sensitivity_analysis.Rmd__par_calib.RDS", compress = "xz") # This code takes 15 minutes to run ```