R/run_model.R

In aj-sykes92/soilc.ipcc: The IPCC (2019) Tier 2 Steady-State Method for Soil Carbon

#' @title \pkg{tidyverse} implementation of the IPCC (2019) steady-state model
#' @description Takes a tibble containing model input data and returns the same tibble with soil
#' carbon stocks estimated. For input data, see \code{toy_input} for example; run
#' \code{stoy_input} for variable interpretation and units.
#' @param data A tibble containing model input data.
#' @param runin_dur An integer value specifying the number of rows to enter into the run-in
#' calculation. Defaults to 10 years.
#' @param drop_prelim A logical indicating whether preliminary calculation results should be dropped
#' from the returned tibble. Defaults to \code{TRUE}.
#' @param drop_runin A logical indicating whether to drop the run-in row from the returned tibble.
#' Run \code{?run_in} for more details. Defaults to \code{TRUE}.
#' @param calculate_climfacs A logical indicating whether to calculate climate factors (\code{tfac}
#' and \code{wfac}) from the object passed to the \code{data} argument.
#' @param year A character string indicating the name of the variable containing the year input.
#' @param tfac A character string indicating the name of the variable containing the temperature
#' factor. Defaults to "tfac". Used to name the temperature factor variable if
#' \code{calculate_climfacs} is \code{TRUE}; used to find it if if \code{calculate_climfacs} is
#' \code{FALSE}.
#' @param wfac A character string indicating the name of the variable containing the water factor.
#' Defaults to "wfac". Used to name the water factor variable if \code{calculate_climfacs} is
#' \code{TRUE}; used to find it if \code{calculate_climfacs} is \code{FALSE}.
#' @param climdata A character string indicating the name of the list column containing climate
#' data. Defaults to "climdata"; unused if \code{calculate_climfacs} is \code{FALSE}.
#' @param temp A character string indicating the name of the (nested) variable containing monthly
#' temperature data. Defaults to "temp"; unused if \code{calculate_climfacs} is \code{FALSE}.
#' @param precip A character string indicating the name of the (nested) variable containing monthly
#' precipitation data. Defaults to "precip"; unused if \code{calculate_climfacs} is \code{FALSE}.
#' @param pet A character string indicating the name of the (nested) variable containing monthly
#' potential evapotranspiration data. Defaults to "pet"; unused if \code{calc_climfacs} is
#' \code{FALSE}.
#' @param c_input A character string naming the variable containing carbon inputs. Defaults to
#' "c_input".
#' @param lignin_frac A character string naming the variable containing lignin fraction. Defaults to
#' "lignin_frac".
#' @param n_frac A character string naming the variable containing nitrogen fraction. Defaults to
#' "n_frac".
#' @param sand_frac A character string naming the variable containing soil sand fraction. Defaults
#' to "sand_frac".
#' @param till_type A character string naming the variable containing tillage practice. Defaults to
#' "till_type".
#' @param params A named list containing model parameters. Defaults to the package default parameter
#' set \code{soilc_params}.
#' @import dplyr
#' @importFrom purrr map_dbl
#' @importFrom rlang :=
#' @examples
#' run_model(soilc.ipcc::toy_input)
#' @export
run_model <- function(
  data,
  runin_dur = 10,
  drop_prelim = TRUE,
  drop_runin = TRUE,
  calculate_climfacs = TRUE,
  year = "year",
  tfac = "tfac", wfac = "wfac",
  climdata = "climdata", temp = "temp", precip = "precip", pet = "pet",
  c_input = "c_input", lignin_frac = "lignin_frac", n_frac = "n_frac",
  sand_frac = "sand_frac", till_type = "till_type",
  params = soilc_params
) {

  # get climate factors if climdata is specfied
  if (isTRUE(calculate_climfacs)) {
    data <- data %>%
      mutate(
        {{ tfac }} := map_dbl(.data[[climdata]], ~tfac(.x[[temp]], params)),
        {{ wfac }} := map_dbl(.data[[climdata]], ~wfac(.x[[precip]], .x[[pet]], params))
      )
  }

  # run in model
  data <- data %>%
    run_in(
      dur = runin_dur,
      year = year,
      climdata = climdata
    )

  # run model
  data <- data %>%
    mutate(

      # model preliminaries + constants
      tillfac = tillfac(
        till_type = .data[[till_type]],
        params = params
      ),

      alpha = alpha(
        c_input = .data[[c_input]],
        lignin_frac = .data[[lignin_frac]],
        n_frac = .data[[n_frac]],
        sand_frac = .data[[sand_frac]],
        till_type = .data[[till_type]],
        params = params
      ),

      k_a = k_a(
        tfac = .data[[tfac]],
        wfac = .data[[wfac]],
        tillfac = .data$tillfac,
        sand_frac = .data[[sand_frac]],
        params = params
      ),

      k_s = k_s(
        tfac = .data[[tfac]],
        wfac = .data[[wfac]],
        tillfac = .data$tillfac,
        params = params
      ),

      k_p = k_p(
        tfac = .data[[tfac]],
        wfac = .data[[wfac]],
        params = params
      ),

      # model steady-state pools
      active_y_ss = active_y_ss(
        k_a = .data$k_a,
        alpha = .data$alpha
      ),

      slow_y_ss = slow_y_ss(
        c_input = .data[[c_input]],
        lignin_frac = .data[[lignin_frac]],
        active_y_ss = .data$active_y_ss,
        k_a = .data$k_a,
        k_s = .data$k_s,
        sand_frac = .data[[sand_frac]],
        params = params
      ),

      passive_y_ss = passive_y_ss(
        active_y_ss = .data$active_y_ss,
        slow_y_ss = .data$slow_y_ss,
        k_a = .data$k_a,
        k_s = .data$k_s,
        k_p = .data$k_p,
        params = params
      ),

      # dynamic soil c pools
      active_y = ss_to_dynamic(
        constant = .data$k_a,
        ss_pool = .data$active_y_ss
      ),

      slow_y = ss_to_dynamic(
        constant = .data$k_s,
        ss_pool = .data$slow_y_ss
      ),

      passive_y = ss_to_dynamic(
        constant = .data$k_p,
        ss_pool = .data$passive_y_ss
      ),

      # roundup and stock change
      total_y = .data$active_y + .data$slow_y + .data$passive_y,
      c_stock_change = soc_stock_change(total_y = .data$total_y)

    )

  # drop elements if specified
  if (isTRUE(drop_prelim)) data <- select(data, - (.data$tillfac:.data$passive_y_ss))
  if (isTRUE(drop_runin)) data <- slice(data, 2:nrow(data))

  return(data)
}