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 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_mode write out all nitrogen state variables (default FALSE)
#' @param replace_input_file_names list with alternative names for output
#'        identifiers to replace the ones in inst/ext_files/metric_files.yml.
#'        e.g. list(npp="mnpp") would replace the expected output for npp with
#'        mnpp followed by the automatically detected file extension (.bin.json)
#' @param suppress_warnings suppress warnings - default: TRUE
#' @param external_variability use externally supplied variability for the
#'        reference period? experimental! (default: FALSE)
#' @param c2vr external variability array
#'
#' @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)
#'
#' @examples
#' \dontrun{
#' ecorisk_wrapper(
#'   path_ref = pnv_folder,
#'   path_scen = run_folder,
#'   read_saved_data = FALSE,
#'   nitrogen = TRUE,
#'   save_data = NULL,
#'   save_ecorisk = NULL,
#'   time_span_reference = c(1550:1579),
#'   time_span_scenario = c(1987:2016)
#' )
#' }
#'
#' @md
#' @export
ecorisk_wrapper <- function(path_ref,
                            path_scen,
                            read_saved_data = FALSE,
                            save_data = NULL,
                            save_ecorisk = NULL,
                            nitrogen = TRUE,
                            weighting = "equal",
                            time_span_reference,
                            time_span_scenario,
                            dimensions_only_local = FALSE,
                            overtime = FALSE,
                            window = 30,
                            debug_mode = FALSE,
                            replace_input_file_names = NULL,
                            external_variability = FALSE,
                            c2vr = NULL,
                            suppress_warnings = TRUE) {
  # check timespan consistency
  nyears <- length(time_span_reference)
  nyears_scen <- length(time_span_scenario)

  # prepare slices for overtime calculation
  slices <- (nyears_scen - window + 1)
  window_half <- round(window / 2.0)
  slice_years <- time_span_scenario[
    (window_half + 1):(nyears_scen - window_half + 1)
  ]

  if ((!nyears == window) || nyears_scen < window) {
    stop(
      "Timespan in reference is not equal to window size (",
      window,
      "), or scenario timespan is smaller than window size."
    )
  }

  if (!read_saved_data) {
    # translate output names (from metric_files.yml) and
    # folders to files_scenarios/reference lists
    metric_files <- system.file(
      "extdata",
      "metric_files.yml",
      package = "biospheremetrics"
    ) %>%
      yaml::read_yaml()

    file_extension <- get_major_file_ext(paste0(path_scen))
    files_scenario <- list()
    files_reference <- list()

    for (output in names(metric_files$metric$ecorisk_nitrogen$output)) {
      # Iterate over all outputs
      if (is.null(replace_input_file_names[[output]])) {
        for (file in metric_files$file_name[[output]]) {
          full_file_path_lu <- paste0(path_scen, file, ".", file_extension)
          if (file.exists(full_file_path_lu)) {
            files_scenario[[output]] <- full_file_path_lu
          }
          full_file_path_pnv <- paste0(path_ref, file, ".", file_extension)
          if (file.exists(full_file_path_pnv)) {
            files_reference[[output]] <- full_file_path_pnv
          }
        }
        if (is.null(files_scenario[[output]])) {
          stop(
            "None of the default file names for ", output,
            " were found in ", path_scen, " - please check or define manually",
            " using argument 'replace_input_file_names'. Stopping."
          )
        }
        if (is.null(files_reference[[output]])) {
          stop(
            "None of the default file names for ", output,
            " were found in ", path_ref, " - please check or define manually",
            " using argument 'replace_input_file_names'. Stopping."
          )
        }
      } else {
        files_scenario[[output]] <- paste0(
          path_scen,
          replace_input_file_names[[output]], ".", file_extension
        )
        files_reference[[output]] <- paste0(
          path_ref,
          replace_input_file_names[[output]], ".", file_extension
        )
      }
    }
  } # end if !read_saved_data
  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)) {
      message("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 EcoRisks
    returned_vars <- read_ecorisk_data(
      files_reference = files_reference,
      files_scenario = files_scenario,
      save_file = save_data,
      nitrogen = nitrogen,
      time_span_reference = time_span_reference,
      time_span_scenario = time_span_scenario,
      debug_mode = debug_mode,
      suppress_warnings = suppress_warnings
    )
    if (overtime == FALSE && slices > 1) {
      print("Overtime == FALSE, but too many time slices.
            I assume this is just a read all data operation. Exiting.")
      return(NULL)
    }

    # 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)
  cell_ids <- 0:(ncells - 1)

  ecorisk <- list(
    ecorisk_total = array(0,
      dim = c(ncells, slices),
      dimnames = list(cell = cell_ids, year = slice_years)
    ),
    vegetation_structure_change = array(0,
      dim = c(ncells, slices),
      dimnames = list(cell = cell_ids, year = slice_years)
    ),
    local_change = array(0,
      dim = c(ncells, slices),
      dimnames = list(cell = cell_ids, year = slice_years)
    ),
    global_importance = array(0,
      dim = c(ncells, slices),
      dimnames = list(cell = cell_ids, year = slice_years)
    ),
    ecosystem_balance = array(0,
      dim = c(ncells, slices),
      dimnames = list(cell = cell_ids, year = slice_years)
    ),
    c2vr = array(0,
      dim = c(4, ncells, slices),
      dimnames = list(part = 1:4, cell = cell_ids, year = slice_years)
    ),
    carbon_stocks = array(0,
      dim = c(ncells, slices),
      dimnames = list(cell = cell_ids, year = slice_years)
    ),
    carbon_fluxes = array(0,
      dim = c(ncells, slices),
      dimnames = list(cell = cell_ids, year = slice_years)
    ),
    carbon_total = array(0,
      dim = c(ncells, slices),
      dimnames = list(cell = cell_ids, year = slice_years)
    ),
    water_total = array(0,
      dim = c(ncells, slices),
      dimnames = list(cell = cell_ids, year = slice_years)
    ),
    water_fluxes = array(0,
      dim = c(ncells, slices),
      dimnames = list(cell = cell_ids, year = slice_years)
    ),
    nitrogen_stocks = array(0,
      dim = c(ncells, slices),
      dimnames = list(cell = cell_ids, year = slice_years)
    ),
    nitrogen_fluxes = array(0,
      dim = c(ncells, slices),
      dimnames = list(cell = cell_ids, year = slice_years)
    ),
    nitrogen_total = array(0,
      dim = c(ncells, slices),
      dimnames = list(cell = cell_ids, year = slice_years)
    ),
    lat = lat,
    lon = lon,
    cell_area = cell_area
  )
  for (y in slice_years) {
    year_range <- as.character((as.numeric(y) - window_half):
                                 (as.numeric(y) + window_half - 1))
    message("Calculating time slice ", y, " of ", slice_years[1], "-", slice_years[length(slice_years)])
    returned <- calc_ecorisk(
      fpc_ref = fpc_ref,
      fpc_scen = fpc_scen[, , year_range],
      bft_ref = bft_ref,
      bft_scen = bft_scen[, , year_range],
      cft_ref = cft_ref,
      cft_scen = cft_scen[, , year_range],
      state_ref = state_ref,
      state_scen = state_scen[, year_range, ],
      weighting = weighting,
      lat = lat,
      lon = lon,
      cell_area = cell_area,
      dimensions_only_local = dimensions_only_local,
      nitrogen = nitrogen,
      external_variability = external_variability,
      c2vr = c2vr
    )
    ecorisk$ecorisk_total[, as.character(y)] <- returned$ecorisk_total
    ecorisk$vegetation_structure_change[, as.character(y)] <- (
      returned$vegetation_structure_change
    )
    ecorisk$local_change[, as.character(y)] <- returned$local_change
    ecorisk$global_importance[, as.character(y)] <- returned$global_importance
    ecorisk$ecosystem_balance[, as.character(y)] <- returned$ecosystem_balance
    ecorisk$c2vr[, , as.character(y)] <- returned$c2vr
    ecorisk$carbon_stocks[, as.character(y)] <- returned$carbon_stocks
    ecorisk$carbon_fluxes[, as.character(y)] <- returned$carbon_fluxes
    ecorisk$carbon_total[, as.character(y)] <- returned$carbon_total
    ecorisk$water_total[, as.character(y)] <- returned$water_total
    ecorisk$water_fluxes[, as.character(y)] <- returned$water_fluxes
    if (nitrogen) {
      ecorisk$nitrogen_stocks[, as.character(y)] <- returned$nitrogen_stocks
      ecorisk$nitrogen_fluxes[, as.character(y)] <- returned$nitrogen_fluxes
      ecorisk$nitrogen_total[, as.character(y)] <- returned$nitrogen_total
    }
  }

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

#' Calculate the ecosystem change metric EcoRisk between 2 sets of states
#' This function is called by the wrapper function (ecorisk_wrapper),
#' unless you know what you are doing, don't use this function directly.
#'
#' 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)
#' @param external_variability include external change_to_variability_ratio?
#'        (default: FALSE)
#' @param c2vr list with external change_to_variability_ratios for each
#'        component (default: NULL)
#'
#' @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,
                         external_variability = FALSE,
                         c2vr = NULL) {
  if (external_variability && is.null(c2vr)) {
    stop("external_variability enabled, but not supplied (c2vr). Aborting.")
  }
  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 seq_len(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 ################
  # dimensions in the state vector
  # 1 "vegetation_carbon_pool"
  # 2 "soil_carbon_pool"
  # 3 "carbon_influx"
  # 4 "carbon_outflux"
  # 5 "soil_water_pool"
  # 6 "water_influx"
  # 7 "water_outflux"
  # 8 "other"
  # 9 "vegetation_nitrogen_pool"
  # 10 "soil_mineral_nitrogen_pool"
  # 11 "nitrogen_influx"
  # 12 "nitrogen_outflux"

  delta_var <- s_change_to_var_ratio(
    vegetation_structure_change,
    vegetation_structure_changesd
  )
  nitrogen_dimensions <- c(
    "vegetation_nitrogen_pool",
    "soil_mineral_nitrogen_pool",
    "nitrogen_influx",
    "nitrogen_outflux"
  )
  all_dimensions <- dimnames(state_scen)$class
  non_nitrogen_dimensions <- setdiff(all_dimensions, nitrogen_dimensions)
  if (nitrogen) {
    lc_raw <- calc_component(
      ref = state_ref,
      scen = state_scen,
      local = TRUE,
      cell_area = cell_area
    ) # local change

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

    eb_raw <- calc_ecosystem_balance(
      ref = state_ref,
      scen = state_scen
    ) # ecosystem balance
  } else {
    lc_raw <- calc_component(
      ref = state_ref[, , non_nitrogen_dimensions],
      scen = state_scen[, , non_nitrogen_dimensions],
      local = TRUE,
      cell_area = cell_area
    ) # local change

    gi_raw <- calc_component(
      ref = state_ref[, , non_nitrogen_dimensions],
      scen = state_scen[, , non_nitrogen_dimensions],
      local = FALSE,
      cell_area = cell_area
    ) # global importance

    eb_raw <- calc_ecosystem_balance(
      ref = state_ref[, , non_nitrogen_dimensions],
      scen = state_scen[, , non_nitrogen_dimensions]
    ) # ecosystem balance
  }
  if (dimensions_only_local == TRUE) {
    # carbon stocks (local change)
    cs <- calc_component(
      ref = state_ref[, , c("vegetation_carbon_pool", "soil_carbon_pool")],
      scen = state_scen[, , c("vegetation_carbon_pool", "soil_carbon_pool")],
      local = TRUE,
      cell_area = cell_area
    )$full
    # carbon fluxes (local change)
    cf <- calc_component(
      ref = state_ref[, , c("carbon_influx", "carbon_outflux")],
      scen = state_scen[, , c("carbon_influx", "carbon_outflux")],
      local = TRUE,
      cell_area = cell_area
    )$full
    # total carbon (local change)
    ct <- calc_component(
      ref = state_ref[, , c(
        "vegetation_carbon_pool",
        "soil_carbon_pool",
        "carbon_influx",
        "carbon_outflux"
      )],
      scen = state_scen[, , c(
        "vegetation_carbon_pool",
        "soil_carbon_pool",
        "carbon_influx",
        "carbon_outflux"
      )],
      local = TRUE,
      cell_area = cell_area
    )$full
    # water fluxes (local change)
    wf <- calc_component(
      ref = state_ref[, , c(
        "water_influx",
        "water_outflux"
      )],
      scen = state_scen[, , c(
        "water_influx",
        "water_outflux"
      )],
      local = TRUE,
      cell_area = cell_area
    )$full
    # total water (local change)
    wt <- calc_component(
      ref = state_ref[, , c(
        "water_influx",
        "water_outflux",
        "soil_water_pool"
      )],
      scen = state_scen[, , c(
        "water_influx",
        "water_outflux",
        "soil_water_pool"
      )],
      local = TRUE,
      cell_area = cell_area
    )$full

    # nitrogen stocks (local change)
    if (nitrogen) {
      ns <- calc_component(
        ref = state_ref[, , c(
          "vegetation_nitrogen_pool",
          "soil_mineral_nitrogen_pool"
        )],
        scen = state_scen[, , c(
          "vegetation_nitrogen_pool",
          "soil_mineral_nitrogen_pool"
        )],
        local = TRUE,
        cell_area = cell_area
      )$full
      # nitrogen fluxes (local change)
      nf <- calc_component(
        ref = state_ref[, , c("nitrogen_influx", "nitrogen_outflux")],
        scen = state_scen[, , c("nitrogen_influx", "nitrogen_outflux")],
        local = TRUE,
        cell_area = cell_area
      )$full
      # total nitrogen (local change)
      nt <- calc_component(
        ref = state_ref[, , c(
          "vegetation_nitrogen_pool",
          "soil_mineral_nitrogen_pool",
          "nitrogen_influx",
          "nitrogen_outflux"
        )],
        scen = state_scen[, , c(
          "vegetation_nitrogen_pool",
          "soil_mineral_nitrogen_pool",
          "nitrogen_influx",
          "nitrogen_outflux"
        )],
        local = TRUE,
        cell_area = cell_area
      )$full
    }
  } else { # local == FALSE
    cf <- (
      calc_component(
        ref = state_ref[, , c(
          "carbon_influx",
          "carbon_outflux"
        )],
        scen = state_scen[, , c(
          "carbon_influx",
          "carbon_outflux"
        )],
        local = TRUE,
        cell_area = cell_area
      )$full + # carbon fluxes
        calc_component(
          ref = state_ref[, , c(
            "carbon_influx",
            "carbon_outflux"
          )],
          scen = state_scen[, , c(
            "carbon_influx",
            "carbon_outflux"
          )],
          local = FALSE,
          cell_area = cell_area
        )$full +
        calc_ecosystem_balance(
          ref = state_ref[, , c(
            "carbon_influx",
            "carbon_outflux"
          )],
          scen = state_scen[, , c(
            "carbon_influx",
            "carbon_outflux"
          )]
        )$full
    ) / 3

    # carbon stocks
    cs <- (
      calc_component(
        ref = state_ref[, , c(
          "vegetation_carbon_pool",
          "soil_carbon_pool"
        )],
        scen = state_scen[, , c(
          "vegetation_carbon_pool",
          "soil_carbon_pool"
        )],
        local = TRUE,
        cell_area = cell_area
      )$full +
        calc_component(
          ref = state_ref[, , c(
            "vegetation_carbon_pool",
            "soil_carbon_pool"
          )],
          scen = state_scen[, , c(
            "vegetation_carbon_pool",
            "soil_carbon_pool"
          )],
          local = FALSE,
          cell_area = cell_area
        )$full +
        calc_ecosystem_balance(
          ref = state_ref[, , c(
            "vegetation_carbon_pool",
            "soil_carbon_pool"
          )],
          scen = state_scen[, , c(
            "vegetation_carbon_pool",
            "soil_carbon_pool"
          )]
        )$full
    ) / 3

    # carbon total
    ct <- (
      calc_component(
        ref = state_ref[, , c(
          "vegetation_carbon_pool",
          "soil_carbon_pool",
          "carbon_influx",
          "carbon_outflux"
        )],
        scen = state_scen[, , c(
          "vegetation_carbon_pool",
          "soil_carbon_pool",
          "carbon_influx",
          "carbon_outflux"
        )],
        local = TRUE,
        cell_area = cell_area
      )$full +
        calc_component(
          ref = state_ref[, , c(
            "vegetation_carbon_pool",
            "soil_carbon_pool",
            "carbon_influx",
            "carbon_outflux"
          )],
          scen = state_scen[, , c(
            "vegetation_carbon_pool",
            "soil_carbon_pool",
            "carbon_influx",
            "carbon_outflux"
          )],
          local = FALSE,
          cell_area = cell_area
        )$full +
        calc_ecosystem_balance(
          ref = state_ref[, , c(
            "vegetation_carbon_pool",
            "soil_carbon_pool",
            "carbon_influx",
            "carbon_outflux"
          )],
          scen = state_scen[, , c(
            "vegetation_carbon_pool",
            "soil_carbon_pool",
            "carbon_influx",
            "carbon_outflux"
          )]
        )$full
    ) / 3

    # water fluxes
    wf <- (
      calc_component(
        ref = state_ref[, , c(
          "water_influx",
          "water_outflux"
        )],
        scen = state_scen[, , c(
          "water_influx",
          "water_outflux"
        )],
        local = TRUE,
        cell_area = cell_area
      )$full +
        calc_component(
          ref = state_ref[, , c(
            "water_influx",
            "water_outflux"
          )],
          scen = state_scen[, , c(
            "water_influx",
            "water_outflux"
          )],
          local = FALSE,
          cell_area = cell_area
        )$full + calc_ecosystem_balance(
          ref = state_ref[, , c(
            "water_influx",
            "water_outflux"
          )],
          scen = state_scen[, , c(
            "water_influx",
            "water_outflux"
          )]
        )$full
    ) / 3

    # water total
    wt <- (
      calc_component(
        ref = state_ref[, , c(
          "water_influx",
          "water_outflux",
          "soil_water_pool"
        )],
        scen = state_scen[, , c(
          "water_influx",
          "water_outflux",
          "soil_water_pool"
        )],
        local = TRUE,
        cell_area = cell_area
      )$full +
        calc_component(
          ref = state_ref[, , c(
            "water_influx",
            "water_outflux",
            "soil_water_pool"
          )],
          scen = state_scen[, , c(
            "water_influx",
            "water_outflux",
            "soil_water_pool"
          )],
          local = FALSE,
          cell_area = cell_area
        )$full +
        calc_ecosystem_balance(
          ref = state_ref[, , c(
            "water_influx",
            "water_outflux",
            "soil_water_pool"
          )],
          scen = state_scen[, , c(
            "water_influx",
            "water_outflux",
            "soil_water_pool"
          )]
        )$full
    ) / 3

    if (nitrogen) {
      # nitrogen stocks (local change)
      ns <- (
        calc_component(
          ref = state_ref[, , c(
            "vegetation_nitrogen_pool",
            "soil_mineral_nitrogen_pool"
          )],
          scen = state_scen[, , c(
            "vegetation_nitrogen_pool",
            "soil_mineral_nitrogen_pool"
          )],
          local = TRUE,
          cell_area = cell_area
        )$full +
          calc_component(
            ref = state_ref[, , c(
              "vegetation_nitrogen_pool",
              "soil_mineral_nitrogen_pool"
            )],
            scen = state_scen[, , c(
              "vegetation_nitrogen_pool",
              "soil_mineral_nitrogen_pool"
            )],
            local = FALSE, cell_area = cell_area
          )$full +
          calc_ecosystem_balance(
            ref = state_ref[, , c(
              "vegetation_nitrogen_pool",
              "soil_mineral_nitrogen_pool"
            )],
            scen = state_scen[, , c(
              "vegetation_nitrogen_pool",
              "soil_mineral_nitrogen_pool"
            )]
          )$full
      ) / 3

      # nitrogen fluxes (local change)
      nf <- (
        calc_component(
          ref = state_ref[, , c(
            "nitrogen_influx",
            "nitrogen_outflux"
          )],
          scen = state_scen[, , c(
            "nitrogen_influx",
            "nitrogen_outflux"
          )],
          local = TRUE,
          cell_area = cell_area
        )$full +
          calc_component(
            ref = state_ref[, , c(
              "nitrogen_influx",
              "nitrogen_outflux"
            )],
            scen = state_scen[, , c(
              "nitrogen_influx",
              "nitrogen_outflux"
            )],
            local = FALSE,
            cell_area = cell_area
          )$full +
          calc_ecosystem_balance(
            ref = state_ref[, , c(
              "nitrogen_influx",
              "nitrogen_outflux"
            )],
            scen = state_scen[, , c(
              "nitrogen_influx",
              "nitrogen_outflux"
            )]
          )$full
      ) / 3

      nt <- (
        calc_component(
          ref = state_ref[, , c(
            "vegetation_nitrogen_pool",
            "soil_mineral_nitrogen_pool",
            "nitrogen_influx",
            "nitrogen_outflux"
          )],
          scen = state_scen[, , c(
            "vegetation_nitrogen_pool",
            "soil_mineral_nitrogen_pool",
            "nitrogen_influx",
            "nitrogen_outflux"
          )],
          local = TRUE,
          cell_area = cell_area
        )$full +
          calc_component(
            ref = state_ref[, , c(
              "vegetation_nitrogen_pool",
              "soil_mineral_nitrogen_pool",
              "nitrogen_influx",
              "nitrogen_outflux"
            )],
            scen = state_scen[, , c(
              "vegetation_nitrogen_pool",
              "soil_mineral_nitrogen_pool",
              "nitrogen_influx",
              "nitrogen_outflux"
            )],
            local = FALSE,
            cell_area = cell_area
          )$full +
          calc_ecosystem_balance(
            ref = state_ref[, , c(
              "vegetation_nitrogen_pool",
              "soil_mineral_nitrogen_pool",
              "nitrogen_influx",
              "nitrogen_outflux"
            )],
            scen = state_scen[, , c(
              "vegetation_nitrogen_pool",
              "soil_mineral_nitrogen_pool",
              "nitrogen_influx",
              "nitrogen_outflux"
            )]
          )$full
      ) / 3
    }
  }
  if (external_variability) {
    delta <- vegetation_structure_change * c2vr["vs", ] # vegetation_structure
    lc <- lc_raw$value * c2vr["lc", ]
    gi <- gi_raw$value * c2vr["gi", ]
    eb <- eb_raw$value * c2vr["eb", ]
  } else {
    delta <- vegetation_structure_change * delta_var # vegetation_structure
    lc <- lc_raw$value * lc_raw$var
    gi <- gi_raw$value * gi_raw$var
    eb <- eb_raw$value * eb_raw$var
    c2vr <- rbind(delta_var, lc_raw$var, gi_raw$var, eb_raw$var) # dim=(4,ncell)
    dimnames(c2vr) <- list(
      component = c("vs", "lc", "gi", "eb"),
      cell = 0:(ncells - 1)
    )
  }

  # 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,
      c2vr = c2vr,
      carbon_stocks = cs,
      carbon_fluxes = cf,
      carbon_total = ct,
      water_fluxes = wf,
      water_stocks = NA,
      water_total = wt,
      nitrogen_stocks = ns,
      nitrogen_fluxes = nf,
      nitrogen_total = nt
    )
  } else {
    ecorisk <- list(
      ecorisk_total = ecorisk_full,
      vegetation_structure_change = delta,
      local_change = lc,
      global_importance = gi,
      ecosystem_balance = eb,
      c2vr = c2vr,
      carbon_stocks = cs,
      carbon_fluxes = cf,
      carbon_total = ct,
      water_fluxes = wf,
      water_stocks = NA,
      water_total = wt,
      nitrogen_stocks = NA,
      nitrogen_fluxes = NA,
      nitrogen_total = NA
    )
  }
  ###
  return(ecorisk)
}

#' Read in output data from LPJmL to calculate the ecosystem change metric
#' EcoRisk. This function is called by the wrapper function (ecorisk_wrapper),
#' unless you know what you are doing, don't use this function directly.
#'
#' 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 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_mode write out all nitrogen state variables (default FALSE)
#' @param suppress_warnings suppress writing of Warnings, default: TRUE
#'