ecorisk.R

# written by Fabian Stenzel, based on work by Sebastian Ostberg
# 2022-2023 - stenzel@pik-potsdam.de

################# EcoRisk calc functions  ###################

#' Wrapper for calculating the ecosystem change metric EcoRisk
#'
#' Function to read in data for ecorisk, and call the calculation function once,
#' if overtime is FALSE, or for each timeslice of length window years, if
#' overtime is TRUE
#'
#' @param path_ref folder of reference run
#' @param path_scen folder of scenario run
#' @param read_saved_data whether to read in previously saved data
#'        (default: FALSE)
#' @param save_data file to save read in data to (default NULL)
#' @param save_ecorisk file to save EcoRisk data to (default NULL)
#' @param nitrogen include nitrogen outputs for pools and fluxes into EcoRisk
#'        calculation (default FALSE)
#' @param weighting apply "old" (Ostberg-like), "new", or "equal" weighting of
#'        vegetation_structure_change weights (default "equal")
#' @param varnames data.frame with names of output files (outname) and time res.
#'        (timestep) -- can be specified to account for variable file names
#'        (default NULL -- standard names as below)
#' @param time_span_reference vector of years to use as scenario period
#' @param time_span_scenario vector of years to use as scenario period
#' @param dimensions_only_local flag whether to use only local change component
#'        for water/carbon/nitrogen fluxes and pools, or use an average of
#'        local change, global change and ecosystem balance (default FALSE)
#' @param overtime logical: calculate ecorisk as time-series? (default: FALSE)
#' @param window integer, number of years for window length (default: 30)
#' @param debug write out all nitrogen state variables (default FALSE)
#'
#' @return list data object containing arrays of ecorisk_total,
#'         vegetation_structure_change, local_change, global_importance,
#'         ecosystem_balance, carbon_stocks, carbon_fluxes, water_fluxes 
#'         (+ nitrogen_stocks and nitrogen_fluxes)
#'
#' @export
ecorisk_wrapper <- function(path_ref,
                            path_scen,
                            read_saved_data = FALSE,
                            save_data = NULL,
                            save_ecorisk = NULL,
                            nitrogen = TRUE,
                            weighting = "equal",
                            varnames = NULL,
                            time_span_reference,
                            time_span_scenario,
                            dimensions_only_local = FALSE,
                            overtime = FALSE,
                            window = 30,
                            debug = FALSE) {
  if (is.null(varnames)) {
    print("variable name list not provided, using standard list, which might
          not be applicable for this case ...")
    varnames <- data.frame(
      row.names = c(
        "grid", "fpc", "fpc_bft", "cftfrac", "firec", "rh_harvest", "npp",
        "evapinterc", "runoff", "transp", "soillitc", "vegc", "swcsum", "firef",
        "rh", "harvestc", "rharvestc", "pft_harvestc", "pft_rharvestc", "evap",
        "interc", "discharge", "soilc", "litc", "swc", "vegn", "soilnh4",
        "soilno3", "leaching", "n2o_denit", "n2o_nit", "n2o_denit", "n2_emis",
        "bnf", "n_volatilization"
      ),
      outname = c(
        "grid.bin", "fpc.bin", "fpc_bft.bin", "cftfrac.bin", "firec.bin",
        "rh_harvest.bin", "npp.bin", "evapinterc.bin", "runoff.bin",
        "transp.bin", "soillitc.bin", "vegc.bin", "swcsum.bin", "firef.bin",
        "rh.bin", "flux_harvest.bin", "flux_rharvest.bin",
        "pft_harvest.pft.bin", "pft_rharvest.pft.bin", "evap.bin", "interc.bin",
        "discharge.bin", "soilc.bin", "litc.bin", "swc.bin", "vegn.bin",
        "soilnh4.bin", "soilno3.bin", "leaching.bin", "n2o_denit.bin",
        "n2o_nit.bin", "n2o_denit.bin", "n2_emis.bin", "bnf.bin",
        "n_volatilization.bin"
      ),
      timestep = c(
        "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y",
        "Y", "Y", "Y", "Y", "Y", "Y", "Y", , "Y", "Y", "Y", "Y", , "Y", "Y",
        "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y"
      )
    )
  }

  nyears <- length(time_span_reference)
  nyears_scen <- length(time_span_scenario)
  if (nyears < 30 || nyears_scen < 30) {
    stop("Warning: timespan in reference or scenario is smaller than 30 years.")
  }
  # translate varnames and folders to files_scenarios/reference lists
  files_scenario <- list(
    grid = paste0(path_scen, varnames["grid", "outname"]),
    fpc = paste0(path_scen, varnames["fpc", "outname"]),
    fpc_bft = paste0(path_scen, varnames["fpc_bft", "outname"]),
    cftfrac = paste0(path_scen, varnames["cftfrac", "outname"]),
    firec = paste0(path_scen, varnames["firec", "outname"]),
    npp = paste0(path_scen, varnames["npp", "outname"]),
    runoff = paste0(path_scen, varnames["runoff", "outname"]),
    transp = paste0(path_scen, varnames["transp", "outname"]),
    vegc = paste0(path_scen, varnames["vegc", "outname"]),
    firef = paste0(path_scen, varnames["firef", "outname"]),
    rh = paste0(path_scen, varnames["rh", "outname"]),
    harvestc = paste0(path_scen, varnames["harvestc", "outname"]),
    rharvestc = paste0(path_scen, varnames["rharvestc", "outname"]),
    pft_harvestc = paste0(path_scen, varnames["pft_harvest", "outname"]),
    pft_rharvestc = paste0(path_scen, varnames["pft_rharvest", "outname"]),
    evap = paste0(path_scen, varnames["evap", "outname"]),
    interc = paste0(path_scen, varnames["interc", "outname"]),
    discharge = paste0(path_scen, varnames["discharge", "outname"]),
    soilc = paste0(path_scen, varnames["soilc", "outname"]),
    litc = paste0(path_scen, varnames["litc", "outname"]),
    swc = paste0(path_scen, varnames["swc", "outname"]),
    vegn = paste0(path_scen, varnames["vegn", "outname"]),
    soilnh4 = paste0(path_scen, varnames["soilnh4", "outname"]),
    soilno3 = paste0(path_scen, varnames["soilno3", "outname"]),
    leaching = paste0(path_scen, varnames["leaching", "outname"]),
    n2o_denit = paste0(path_scen, varnames["n2o_denit", "outname"]),
    n2o_nit = paste0(path_scen, varnames["n2o_nit", "outname"]),
    n2_emis = paste0(path_scen, varnames["n2_emis", "outname"]),
    bnf = paste0(path_scen, varnames["bnf", "outname"]),
    n_volatilization = paste0(path_scen, varnames["n_volatilization", "outname"])
  )
  files_reference <- list(
    grid = paste0(path_ref, varnames["grid", "outname"]),
    fpc = paste0(path_ref, varnames["fpc", "outname"]),
    fpc_bft = paste0(path_ref, varnames["fpc_bft", "outname"]),
    cftfrac = paste0(path_ref, varnames["cftfrac", "outname"]),
    firec = paste0(path_ref, varnames["firec", "outname"]),
    npp = paste0(path_ref, varnames["npp", "outname"]),
    runoff = paste0(path_ref, varnames["runoff", "outname"]),
    transp = paste0(path_ref, varnames["transp", "outname"]),
    vegc = paste0(path_ref, varnames["vegc", "outname"]),
    firef = paste0(path_ref, varnames["firef", "outname"]),
    rh = paste0(path_ref, varnames["rh", "outname"]),
    harvestc = paste0(path_ref, varnames["harvestc", "outname"]),
    rharvestc = paste0(path_ref, varnames["rharvestc", "outname"]),
    pft_harvestc = paste0(path_ref, varnames["pft_harvest", "outname"]),
    pft_rharvestc = paste0(path_ref, varnames["pft_rharvest", "outname"]),
    evap = paste0(path_ref, varnames["evap", "outname"]),
    interc = paste0(path_ref, varnames["interc", "outname"]),
    discharge = paste0(path_ref, varnames["discharge", "outname"]),
    soilc = paste0(path_ref, varnames["soilc", "outname"]),
    litc = paste0(path_ref, varnames["litc", "outname"]),
    swc = paste0(path_ref, varnames["swc", "outname"]),
    vegn = paste0(path_ref, varnames["vegn", "outname"]),
    soilnh4 = paste0(path_ref, varnames["soilnh4", "outname"]),
    soilno3 = paste0(path_ref, varnames["soilno3", "outname"]),
    leaching = paste0(path_ref, varnames["leaching", "outname"]),
    n2o_denit = paste0(path_ref, varnames["n2o_denit", "outname"]),
    n2o_nit = paste0(path_ref, varnames["n2o_nit", "outname"]),
    n2_emis = paste0(path_ref, varnames["n2_emis", "outname"]),
    bnf = paste0(path_ref, varnames["bnf", "outname"]),
    n_volatilization = paste0(path_ref, varnames["n_volatilization", "outname"])
  )

  if (overtime && window != nyears) stop("Overtime is enabled, but window \
                  length (", window, ") does not match the reference nyears.")

  if (read_saved_data) {
    if (!is.null(save_data)) {
      print(paste("Loading saved data from:", save_data))
      load(file = save_data)
    } else {
      stop("save_data is not specified as parameter, ",
           "nothing to load ... exiting")
    }
  } else {
    # first read in all lpjml output files required for computing EcoRisk
    returned_vars <- read_ecorisk_data(
      files_reference = files_reference,
      files_scenario = files_scenario,
      save_file = save_data,
      export = FALSE,
      nitrogen = nitrogen,
      time_span_reference = time_span_reference,
      time_span_scenario = time_span_scenario,
      debug = debug
    )
    # extract variables from return list object and give them proper names
    state_ref <- returned_vars$state_ref
    state_scen <- returned_vars$state_scen
    fpc_ref <- returned_vars$fpc_ref
    fpc_scen <- returned_vars$fpc_scen
    bft_ref <- returned_vars$bft_ref
    bft_scen <- returned_vars$bft_scen
    cft_ref <- returned_vars$cft_ref
    cft_scen <- returned_vars$cft_scen
    lat <- returned_vars$lat
    lon <- returned_vars$lon
    cell_area <- returned_vars$cell_area
    rm(returned_vars)
  }

  ncells <- length(cell_area)
  slices <- (nyears_scen - window + 1)
  ecorisk <- list(
    ecorisk_total = array(0, dim = c(ncells, slices)),
    vegetation_structure_change = array(0, dim = c(ncells, slices)),
    local_change = array(0, dim = c(ncells, slices)),
    global_importance = array(0, dim = c(ncells, slices)),
    ecosystem_balance = array(0, dim = c(ncells, slices)),
    carbon_stocks = array(0, dim = c(ncells, slices)),
    carbon_fluxes = array(0, dim = c(ncells, slices)),
    water_stocks = array(0, dim = c(ncells, slices)),
    water_fluxes = array(0, dim = c(ncells, slices)),
    nitrogen_stocks = array(0, dim = c(ncells, slices)),
    nitrogen_fluxes = array(0, dim = c(ncells, slices))
  )
  for (y in 1:slices) {
    print(paste0("Calculating time slice ", y, " of ", slices))
    returned <- calc_ecorisk(
      fpc_ref = fpc_ref,
      fpc_scen = fpc_scen[, , y:(y + window - 1)],
      bft_ref = bft_ref,
      bft_scen = bft_scen[, , y:(y + window - 1)],
      cft_ref = cft_ref,
      cft_scen = cft_scen[, , y:(y + window - 1)],
      state_ref = state_ref,
      state_scen = state_scen[, y:(y + window - 1), ],
      weighting = weighting,
      lat = lat,
      lon = lon,
      cell_area = cell_area,
      dimensions_only_local = dimensions_only_local,
      nitrogen = nitrogen
    )
    ecorisk$ecorisk_total[, y] <- returned$ecorisk_total
    ecorisk$vegetation_structure_change[, y] <- (
      returned$vegetation_structure_change
    )
    ecorisk$local_change[, y] <- returned$local_change
    ecorisk$global_importance[, y] <- returned$global_importance
    ecorisk$ecosystem_balance[, y] <- returned$ecosystem_balance
    ecorisk$carbon_stocks[, y] <- returned$carbon_stocks
    ecorisk$carbon_fluxes[, y] <- returned$carbon_fluxes
    ecorisk$water_stocks[, y] <- returned$water_stocks
    ecorisk$water_fluxes[, y] <- returned$water_fluxes
    if (nitrogen) {
      ecorisk$nitrogen_stocks[, y] <- returned$nitrogen_stocks
      ecorisk$nitrogen_fluxes[, y] <- returned$nitrogen_fluxes
    }
  }


  ############## export and save data if requested #############
  if (!(is.null(save_ecorisk))) {
    print(paste0("Saving EcoRisk data to: ", save_ecorisk))
    save(ecorisk, file = save_ecorisk)
  }
  #
  ###
  return(ecorisk)
}

#' Calculate the ecosystem change metric EcoRisk between 2 sets of states
#'
#' Function to calculate the ecosystem change metric EcoRisk, based on
#' gamma/vegetation_structure_change
#' work from Sykes (1999), Heyder (2011), and Ostberg (2015,2018).
#' This is a reformulated version in R, not producing 100% similar values
#' than the C/bash version from Ostberg et al. 2018, but similar the methodology
#'
#' @param fpc_ref reference run data for fpc
#' @param fpc_scen scenario run data for fpc
#' @param bft_ref reference run data for fpc_bft
#' @param bft_scen scenario run data for fpc_bft
#' @param cft_ref reference run data for cftfrac
#' @param cft_scen scenario run data for cftfrac
#' @param state_ref reference run data for state variables
#' @param state_scen scenario run data for state variables
#' @param weighting apply "old" (Ostberg-like), "new", or "equal" weighting of
#'        vegetation_structure_change weights (default "equal")
#' @param lat latitude array
#' @param lon longitude array
#' @param cell_area cellarea array
#' @param dimensions_only_local flag whether to use only local change component
#'        for water/carbon/nitrogen fluxes and pools, or use an average of
#'        local change, global change and ecosystem balance (default FALSE)
#' @param nitrogen include nitrogen outputs (default: TRUE)
#'
#' @return list data object containing arrays of ecorisk_total,
#'         vegetation_structure_change, local_change, global_importance,
#'         ecosystem_balance, carbon_stocks, carbon_fluxes, water_fluxes
#'         (+ nitrogen_stocks and nitrogen_fluxes)
#'
#' @export
calc_ecorisk <- function(fpc_ref,
                        fpc_scen,
                        bft_ref,
                        bft_scen,
                        cft_ref,
                        cft_scen,
                        state_ref,
                        state_scen,
                        weighting = "equal",
                        lat,
                        lon,
                        cell_area,
                        dimensions_only_local = FALSE,
                        nitrogen = TRUE) {
  di_ref <- dim(fpc_ref)
  di_scen <- dim(fpc_scen)
  ncells <- di_ref[1]
  nyears <- di_ref[3]
  if (di_ref[3] != di_scen[3]) {
    stop("Dimension year does not match between fpc_scen and fpc_ref.")
  }
  # calc vegetation_structure_change and variability of
  #   vegetation_structure_change within
  # reference period S(vegetation_structure_change,
  #   sigma_vegetation_structure_change)
  fpc_ref_mean <- apply(fpc_ref, c(1, 2), mean)
  bft_ref_mean <- apply(bft_ref, c(1, 2), mean)
  cft_ref_mean <- apply(cft_ref, c(1, 2), mean)


  sigma_vegetation_structure_change_ref_list <- array(
    0, dim = c(ncells, nyears)
  )
  # calculate for every year of the reference period,
  #   vegetation_structure_change between that year and the average reference
  #   period year
  # this gives the variability of vegetation_structure_change within the
  #   reference period
  for (y in 1:nyears) {
    sigma_vegetation_structure_change_ref_list[, y] <- calc_delta_v( # nolint
      fpc_ref = fpc_ref_mean,
      fpc_scen = fpc_ref[, , y],
      bft_ref = bft_ref_mean,
      bft_scen = bft_ref[, , y],
      cft_ref = cft_ref_mean,
      cft_scen = cft_ref[, , y],
      weighting = weighting
    )
  }

  # calculate the std deviation over the reference period for each gridcell
  vegetation_structure_changesd <- apply(
    sigma_vegetation_structure_change_ref_list,
    c(1),
    stats::sd
  )

  # calculate vegetation_structure_change between average reference and average
  #   scenario period
  vegetation_structure_change <- calc_delta_v(
    fpc_ref = fpc_ref_mean,
    fpc_scen = apply(fpc_scen, c(1, 2), mean),
    bft_ref = bft_ref_mean,
    bft_scen = apply(bft_scen, c(1, 2), mean),
    cft_ref = cft_ref_mean,
    cft_scen = apply(cft_scen, c(1, 2), mean),
    weighting = weighting
  )
  #
  ####
  ############## calc EcoRisk components ################
  # variable names for the state vector
  # 1:3 carbon fluxes
  # 4:6 water fluxes
  # 7:8 carbon pools/stocks,
  # 9:10 water pools
  # 11 additional variables for global/local difference, but not included in
  #   stocks/fluxes
  # 12:13 nitrogen pools/stocks
  # 14:16 nitrogen fluxes

  delta <- vegetation_structure_change * s_change_to_var_ratio(
    vegetation_structure_change,
    vegetation_structure_changesd
  ) # vegetation_structure_change

  lc <- calc_component(
    ref = state_ref,
    scen = state_scen,
    local = TRUE,
    cell_area = cell_area
  ) # local change

  gi <- calc_component(
    ref = state_ref,
    scen = state_scen,
    local = FALSE,
    cell_area = cell_area
  ) # global importance

  eb <- calc_ecosystem_balance(
    ref = state_ref,
    scen = state_scen
  ) # ecosystem balance

  if (dimensions_only_local == TRUE) {

    # carbon fluxes (local change)
    cf <- calc_component(
      ref = state_ref[, , 1:3],
      scen = state_scen[, , 1:3],
      local = TRUE,
      cell_area = cell_area)

    # carbon stocks (local change)
    cs <- calc_component(
      ref = state_ref[, , 7:8],
      scen = state_scen[, , 7:8],
      local = TRUE,
      cell_area = cell_area
    )

    # water fluxes (local change)
    wf <- calc_component(
      ref = state_ref[, , 4:6],
      scen = state_scen[, , 4:6],
      local = TRUE,
      cell_area = cell_area
    ) 

    # water pools (local change)
    ws <- calc_component(
      ref = state_ref[, , 9:10],
      scen = state_scen[, , 9:10],
      local = TRUE,
      cell_area = cell_area
    )

    # nitrogen stocks (local change)
    if (nitrogen) {
      ns <- calc_component(
        ref = state_ref[, , 12:13],
        scen = state_scen[, , 12:13],
        local = TRUE,
        cell_area = cell_area
      )

      # nitrogen fluxes (local change)
      nf <- calc_component(
        ref = state_ref[, , 14:16],
        scen = state_scen[, , 14:16],
        local = TRUE,
        cell_area = cell_area
      )
    }

  } else {
    cf <- (
      calc_component(
        ref = state_ref[, , 1:3],
        scen = state_scen[, , 1:3],
        local = TRUE,
        cell_area = cell_area
      ) + # carbon fluxes
      calc_component(
        ref = state_ref[, , 1:3],
        scen = state_scen[, , 1:3],
        local = FALSE,
        cell_area = cell_area
      ) +
      calc_ecosystem_balance(
        ref = state_ref[, , 1:3],
        scen = state_scen[, , 1:3]
      )
    ) / 3

    # carbon stocks
    cs <- (
      calc_component(
        ref = state_ref[, , 7:8],
        scen = state_scen[, , 7:8],
        local = TRUE,
        cell_area = cell_area
      ) +
      calc_component(
        ref = state_ref[, , 7:8],
        scen = state_scen[, , 7:8],
        local = FALSE,
        cell_area = cell_area
      ) +
      calc_ecosystem_balance(
        ref = state_ref[, , 7:8],
        scen = state_scen[, , 7:8]
      )
    ) / 3

    # water fluxes
    wf <- (
      calc_component(
        ref = state_ref[, , 4:6],
        scen = state_scen[, , 4:6],
        local = TRUE,
        cell_area = cell_area
      ) +
      calc_component(
        ref = state_ref[, , 4:6],
        scen = state_scen[, , 4:6],
        local = FALSE,
        cell_area = cell_area
      ) + calc_ecosystem_balance(
        ref = state_ref[, , 4:6],
        scen = state_scen[, , 4:6]
      )
    ) / 3

    # water pools
    ws <- (
      calc_component(
        ref = state_ref[, , 9:10],
        scen = state_scen[, , 9:10],
        local = TRUE,
        cell_area = cell_area
      ) +
      calc_component(
        ref = state_ref[, , 9:10],
        scen = state_scen[, , 9:10],
        local = FALSE,
        cell_area = cell_area
      ) +
      calc_ecosystem_balance(
        ref = state_ref[, , 9:10],
        scen = state_scen[, , 9:10]
      )
    ) / 3

    if (nitrogen) {

      # nitrogen stocks (local change)
      ns <- (
        calc_component(
          ref = state_ref[, , 12:13],
          scen = state_scen[, , 12:13],
          local = TRUE,
          cell_area = cell_area
        ) +
        calc_component(
          ref = state_ref[, , 12:13],
          scen = state_scen[, , 12:13],
          local = FALSE, cell_area = cell_area
          ) +
        calc_ecosystem_balance(
          ref = state_ref[, , 12:13],
          scen = state_scen[, , 12:13]
        )
      ) / 3

      # nitrogen fluxes (local change)
      nf <- (
        calc_component(
          ref = state_ref[, , 14:16],
          scen = state_scen[, , 14:16],
          local = TRUE,
          cell_area = cell_area
        ) +
        calc_component(
          ref = state_ref[, , 14:16],
          scen = state_scen[, , 14:16],
          local = FALSE,
          cell_area = cell_area
        ) +
        calc_ecosystem_balance(
          ref = state_ref[, , 14:16],
          scen = state_scen[, , 14:16]
        )
      ) / 3
    }
  }

  # calc total EcoRisk as the average of the 4 components
  ecorisk_full <- (delta + lc + gi + eb) / 4 # check for NAs

  if (nitrogen) {
    ecorisk <- list(
      ecorisk_total = ecorisk_full,
      vegetation_structure_change = delta,
      local_change = lc,
      global_importance = gi,
      ecosystem_balance = eb,
      carbon_stocks = cs,
      carbon_fluxes = cf,
      water_fluxes = wf,
      water_stocks = ws,
      nitrogen_stocks = ns,
      nitrogen_fluxes = nf
    )

  } else {
    ecorisk <- list(
      ecorisk_total = ecorisk_full,
      vegetation_structure_change = delta,
      local_change = lc,
      global_importance = gi,
      ecosystem_balance = eb,
      carbon_stocks = cs,
      carbon_fluxes = cf,
      water_fluxes = wf,
      water_stocks = ws
    )
  }
  ###
  return(ecorisk)
}

#' Read in output data from LPJmL to calculate the ecosystem change metric
#' EcoRisk
#'
#' Utility function to read in output data from LPJmL for calculation of EcoRisk
#'
#' @param files_reference folder of reference run
#' @param files_scenario folder of scenario run
#' @param save_file file to save read in data to (default NULL)
#' @param export flag whether to export réad in data to global environment
#'               (default FALSE)
#' @param time_span_reference vector of years to use as scenario period
#' @param time_span_scenario vector of years to use as scenario period
#' @param nitrogen include nitrogen outputs for pools and fluxes into EcoRisk
#'                 calculation (default FALSE)
#' @param debug write out all nitrogen state variables (default FALSE)
#'
#' @return list data object containing arrays of state_ref, mean_state_ref,
#'         state_scen, mean_state_scen, fpc_ref, fpc_scen, bft_ref, bft_scen,
#'         cft_ref, cft_scen, lat, lon, cell_area
#'
#' @export
read_ecorisk_data <- function(files_reference, # nolint
                              files_scenario,
                              save_file = NULL,
                              export = FALSE,
                              time_span_reference,
                              time_span_scenario,
                              nitrogen,
                              debug = FALSE) {
  file_type <- tools::file_ext(files_reference$grid)

  if (file_type %in% c("json", "clm")) {
    # read grid
    grid <- lpjmlkit::read_io(
      files_reference$grid,
    )
    # calculate cell area
    cell_area <- lpjmlkit::calc_cellarea(grid)
    lat <- grid$data[, , 2]
    lon <- grid$data[, , 1]

    ### read in lpjml output
    # for vegetation_structure_change (fpc,fpc_bft,cftfrac)
    print("Reading in fpc, fpc_bft, cftfrac")

    cft_scen <- aperm(lpjmlkit::read_io(
      files_scenario$cftfrac,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)), c(1, 3, 2))

    bft_scen <- aperm(lpjmlkit::read_io(
      files_scenario$fpc_bft,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)), c(1, 3, 2))

    fpc_scen <- aperm(lpjmlkit::read_io(
      files_scenario$fpc,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)), c(1, 3, 2))

    if (file.exists(files_reference$cftfrac)) {
      cft_ref <- aperm(lpjmlkit::read_io(
        files_reference$cftfrac,
        subset = list(year = as.character(time_span_reference))
      ) %>%
        lpjmlkit::transform(to = c("year_month_day")) %>%
        lpjmlkit::as_array(aggregate = list(month = sum)), c(1, 3, 2))
    } else {
      cft_ref <- cft_scen * 0
    }

    if (file.exists(files_reference$fpc_bft)) {
      bft_ref <- aperm(lpjmlkit::read_io(
        files_reference$fpc_bft,
        subset = list(year = as.character(time_span_reference))
      ) %>%
        lpjmlkit::transform(to = c("year_month_day")) %>%
        lpjmlkit::as_array(aggregate = list(month = sum)), c(1, 3, 2))
    } else {
      bft_ref <- bft_scen * 0
    }

    fpc_ref <- aperm(lpjmlkit::read_io(
      files_reference$fpc,
      subset = list(year = as.character(time_span_reference))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)), c(1, 3, 2))

    # cffiles = ( firec rh_harvest npp ) yearly carbon fluxes
    print("Reading in firec, rh_harvest, npp")

    rh_scen <- lpjmlkit::read_io(
      files_scenario$rh,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    if (file.exists(files_scenario$harvestc)) {
      harvest_scen <- lpjmlkit::read_io(
        files_scenario$harvestc,
        subset = list(year = as.character(time_span_scenario))
      ) %>%
        lpjmlkit::transform(to = c("year_month_day")) %>%
        lpjmlkit::as_array(aggregate = list(month = sum)) %>%
        drop()

    } else if (file.exists(files_scenario$pft_harvestc)) {
      harvest_scen <- lpjmlkit::read_io(
        files_scenario$pft_harvestc,
        subset = list(year = as.character(time_span_scenario))
      ) %>%
        lpjmlkit::transform(to = c("year_month_day")) %>%
        lpjmlkit::as_array(aggregate = list(month = sum, band = sum)) %>%
        drop()

    } else {
      stop("Missing harvestc output in scenario folder.")
    }

    if (file.exists(files_scenario$rharvestc)) {
      rharvest_scen <- lpjmlkit::read_io(
        files_scenario$rharvestc,
        subset = list(year = as.character(time_span_scenario))
      ) %>%
        lpjmlkit::transform(to = c("year_month_day")) %>%
        lpjmlkit::as_array(aggregate = list(month = sum)) %>%
        drop()

    } else if (file.exists(files_scenario$pft_rharvestc)) {
      rharvest_scen <- lpjmlkit::read_io(
        files_scenario$pft_rharvestc,
        subset = list(year = as.character(time_span_scenario))
      ) %>%
        lpjmlkit::transform(to = c("year_month_day")) %>%
        lpjmlkit::as_array(aggregate = list(month = sum, band = sum)) %>%
        drop()

    } else {
      stop("Missing rharvestc output in scenario folder.")
    }

    rh_harvest_scen <- rh_scen + harvest_scen + rharvest_scen

    firec_scen <- lpjmlkit::read_io(
      files_scenario$firec,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    npp_scen <- lpjmlkit::read_io(
      files_scenario$npp,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    rh_ref <- lpjmlkit::read_io(
      files_reference$rh,
      subset = list(year = as.character(time_span_reference))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    if (file.exists(files_reference$harvestc)) {
      harvest_ref <- lpjmlkit::read_io(
        files_reference$harvestc,
        subset = list(year = as.character(time_span_reference))
      ) %>%
        lpjmlkit::transform(to = c("year_month_day")) %>%
        lpjmlkit::as_array(aggregate = list(month = sum)) %>%
        drop()

    } else if (file.exists(files_reference$pft_harvestc)) {
      harvest_ref <- lpjmlkit::read_io(
        files_reference$pft_harvestc,
        subset = list(year = as.character(time_span_scenario))
      ) %>%
        lpjmlkit::transform(to = c("year_month_day")) %>%
        lpjmlkit::as_array(aggregate = list(month = sum, band = sum)) %>%
        drop()

    } else {
      print("No harvest output available for reference period, setting to 0.")
      harvest_ref <- harvest_scen * 0
    }

    if (file.exists(files_reference$rharvestc)) {
      rharvest_ref <- lpjmlkit::read_io(
        files_reference$rharvestc,
        subset = list(year = as.character(time_span_reference))
      ) %>%
        lpjmlkit::transform(to = c("year_month_day")) %>%
        lpjmlkit::as_array(aggregate = list(month = sum)) %>%
        drop()

    } else if (file.exists(files_reference$pft_rharvestc)) {
      rharvest_ref <- lpjmlkit::read_io(
        files_reference$pft_rharvestc,
        subset = list(year = as.character(time_span_scenario))
      ) %>%
        lpjmlkit::transform(to = c("year_month_day")) %>%
        lpjmlkit::as_array(aggregate = list(month = sum, band = sum)) %>%
        drop()

    } else {
      print("No rharvest output available for reference period, setting to 0.")
      rharvest_ref <- rharvest_scen * 0
    }

    rh_harvest_ref <- rh_ref + harvest_ref + rharvest_ref

    firec_ref <- lpjmlkit::read_io(
      files_reference$firec,
      subset = list(year = as.character(time_span_reference))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    npp_ref <- lpjmlkit::read_io(
      files_reference$npp,
      subset = list(year = as.character(time_span_reference))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    # wffiles = (evapinterc runoff transp) - yearly water fluxes
    print("Reading in evapinterc, runoff, transp")

    evap_ref <- lpjmlkit::read_io(
      files_reference$evap,
      subset = list(year = as.character(time_span_reference))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    interc_ref <- lpjmlkit::read_io(
      files_reference$interc,
      subset = list(year = as.character(time_span_reference))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    evapinterc_ref <- evap_ref + interc_ref

    runoff_ref <- lpjmlkit::read_io(
      files_reference$runoff,
      subset = list(year = as.character(time_span_reference))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    transp_ref <- lpjmlkit::read_io(
      files_reference$transp,
      subset = list(year = as.character(time_span_reference))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    evap_scen <- lpjmlkit::read_io(
      files_scenario$evap,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    interc_scen <- lpjmlkit::read_io(
      files_scenario$interc,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    evapinterc_scen <- evap_scen + interc_scen

    runoff_scen <- lpjmlkit::read_io(
      files_scenario$runoff,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    transp_scen <- lpjmlkit::read_io(
      files_scenario$transp,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    # csfiles = ( soillitc vegc_avg ) #carbon pools
    print("Reading in soillitc, vegc")

    soil_ref <- lpjmlkit::read_io(
      files_reference$soilc,
      subset = list(year = as.character(time_span_reference))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    litc_ref <- lpjmlkit::read_io(
      files_reference$litc,
      subset = list(year = as.character(time_span_reference))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    soillitc_ref <- soil_ref + litc_ref

    vegc_ref <- lpjmlkit::read_io(
      files_reference$vegc,
      subset = list(year = as.character(time_span_reference))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    soil_scen <- lpjmlkit::read_io(
      files_scenario$soilc,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    litc_scen <- lpjmlkit::read_io(
      files_scenario$litc,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    soillitc_scen <- soil_scen + litc_scen

    vegc_scen <- lpjmlkit::read_io(
      files_scenario$vegc,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    # water pools = (swcsum discharge)
    print("Reading in swcsum, discharge")
    swcsum_ref <- lpjmlkit::read_io(
      files_reference$swc,
      subset = list(year = as.character(time_span_reference))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum, band = sum)) %>%
      drop()

    swcsum_scen <- lpjmlkit::read_io(
      files_scenario$swc,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum, band = sum)) %>%
      drop()

    discharge_ref <- lpjmlkit::read_io(
      files_reference$discharge,
      subset = list(year = as.character(time_span_reference))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    discharge_scen <- lpjmlkit::read_io(
      files_scenario$discharge,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)) %>%
      drop()

    print("Reading in firef")