R/IFNscenario.R

In miquelcaceres/IFNdyn: IFN forest dynamics and scenario analysis

#' IFN scenario simulation engine
#'
#' Performs simulations of forest dynamics including management prescriptions (on a demand-basis) for a set of forest plots
#'
#' @param treeData A data frame with tree records in rows and columns 'ID', 'Species', 'DBH' (in cm), 'H (in m)' and 'N' (ha-1)
#' @param plotData A data frame with plot 'ID' in rows and columns 'SWHC' (mm), 'Rad' (MJ), 'Temp' (degrees C), 'Prec' (mm), 'PET' (mm), 'slope' (degrees), 'Managed' (boolean) and 'Estrato' (stratum code).
#'                 It can also include a column 'Area' (ha) with the surface corresponding to each forest plot. In this case, 'strata' will not be used.
#' @param plotDataDyn A data frame with dynamic plot data, with columns 'ID', 'Step', 'Temp' (degrees C), 'Prec' (mm) and 'PET' (mm)
#' @param prescription A data frame with plots in rows (row names should be equal to 'ID' in tree data) and prescription options in columns (see \code{\link{IFNmanagement}})
#' @param strata A data frame with surface areas per stratum (see \code{\link{estratosIFN3}}). Used to determine the surface corresponding to each plot.
#' @param demand A data frame with annual volume to be extracted by species (see \code{\link{defaultDemand}})
#' @param nb An object of class \code{listw} (see functions \code{\link{IFNknn}} and \code{\link{nb2listw}}).
#' @param numSteps Number of steps to simulate (1 step = 10 years)
#' @param demandOffset Initial demand offset (for simulations continuing previous ones).
#' @param stepOffset Offset of time steps to read dynamic plot data (for simulations continuing previous ones).
#' @param as.CO2 Boolean flag to indicate output as kg of CO2 instead of kg of dry weight (see \code{\link{IFNbiomass}}).
#' @param sequence Boolean flag to indicate that the sequence of tree lists by step should be returned.
#' If \code{sequence = TRUE}, an additional data frame is included in the output, called \code{treeDataSequence}.
#' This table includes a column 'Step' to indicate the simulation step.
#' @param dead.table Boolean flag to indicate that a table of trees dead during the simulation is also desired.
#' This table includes a column 'Step' to indicate at which step the trees died.
#' @param cut.table Boolean flag to indicate that a table of trees cut during the simulation is also desired.
#' This table includes a column 'Step' to indicate at which step the trees were cut.
#' @param verbose Boolean flag to indicate extra console output (reporting management interventions in each plot)
#' @param ... Additional options to be passed to function \code{\link{IFNdynamics}}
#'
#' @return A list with the following items:
#'  \itemize{
#'    \item{\code{plots}: number of IFN plots included in the simulation.}
#'    \item{\code{totalArea}: total area (in ha) that the set of plots represent.}
#'    \item{\code{plotDataSuppl}: A data frame with supplementary plot data, with columns:
#'      \itemize{
#'         \item{\code{Area}: Area (in ha) represented by each plot.}
#'         \item{\code{DominantSpecies}: Code corresponding to the dominant species of each plot.}
#'         \item{\code{DemandRow}: Demand species row index corresponding to each plot.}
#'         \item{\code{PrescriptionRow}: Prescription row index corresponding to each plot.}
#'         \item{\code{finalPreviousStage}: Previous stage of the final cuts (for plots with regular management).}
#'      }
#'    }
#'    \item{\code{steps}: number of 10-yr steps simulated.}
#'    \item{\code{extractedDemand}: a data frame with the volume (m3) extracted for each demand species and time step.
#'                The three last columns include: \code{AnnualDemand} the mean annual extraction (simulated), \code{Demand}  the annual demand (input), and \code{Coverage} the coverage of the demand (in percent), for each species.
#'                The last row ('Total') reports statistics across species.}
#'    \item{\code{managedArea}: a data frame with the area (ha) managed for each demand species and time step.
#'                The last row ('Total') reports statistics across species.}
#'    \item{\code{volume, density, rootBiomass, stemBiomass, branchBiomass}: simulation results by variable. Volume is in m3, density in individuals and biomass in Mg.
#'     Each of these elements is in turn a list with the following elements:
#'       \itemize{
#'         \item{\code{summary}: data frame with a summary of results for each time step (standing, variation, dead and extracted). Here values are given per hectare.}
#'         \item{\code{standingSpp}: volume, biomass or density of living trees, per species and time step.}
#'         \item{\code{standingPlot}: volume, biomass or density of living trees, per plot and time step.}
#'         \item{\code{standingDomSpp}: volume, biomass or density of living trees, per dominant species (i.e. pooling all plots with the same dominant species) and time step.}
#'         \item{\code{variationSpp}: variation (i.e. difference between the current time step and the previous one) in volume, biomass or density of living trees, per species and time step.}
#'         \item{\code{variationPlot}: variation (i.e. difference between the current time step and the previous one) in volume, biomass or density of living trees, per dominant species (i.e. pooling all plots with the same dominant species) and time step.}
#'         \item{\code{variationDomSpp}: variation (i.e. difference between the current time step and the previous one) in volume, biomass or density of living trees, per plot and time step.}
#'         \item{\code{standgrowthSpp}: stand growth in volume, biomass or density of living trees, per species and time step.}
#'         \item{\code{standgrowthPlot}: stand growth in volume, biomass or density of living trees, per plot and time step.}
#'         \item{\code{standgrowthDomSpp}: stand growth in volume, biomass or density of living trees, per dominant species (i.e. pooling all plots with the same dominant species) and time step.}
#'         \item{\code{deadSpp}: volume, biomass or density of dead trees (i.e. dead during the corresponding time step), per species and time step.}
#'         \item{\code{deadPlot}: volume, biomass or density of dead trees (i.e. dead during the corresponding time step), per plot and time step.}
#'         \item{\code{deadDomSpp}: volume, biomass or density of dead trees (i.e. dead during the corresponding time step), per dominant species (i.e. pooling all plots with the same dominant species) and time step.}
#'         \item{\code{extractedSpp}: volume, biomass or density of trees extracted (during the corresponding time step), per species and time step.}
#'         \item{\code{extractedPlot}: volume, biomass or density of trees extracted (during the corresponding time step), per plot and time step.}
#'         \item{\code{extractedDomSpp}: volume, biomass or density of trees extracted (during the corresponding time step), per dominant species (i.e. pooling all plots with the same dominant species) and time step.}
#'       }
#'     }
#'    \item{\code{treeDataNew}: tree data frame at the end of the simulated period.}
#'    \item{\code{treeDataSequence}: tree data frame for the whole sequence of steps (if \code{sequence = TRUE}).}
#'    \item{\code{treeDataDead}: tree data frame with dead trees for every step (if \code{dead.table = TRUE}).}
#'    \item{\code{treeDataCut}: tree data frame with cut trees for every step (if \code{cut.table = TRUE}).}
#'  }
#'
#' @examples
#'
#' \dontrun{
#' library(IFNread)
#' data(defaultDemand)
#' data(defaultPrescription)
#' data(estratosIFN3)
#' data(plotDataIFN)
#'
#' accessfilesIFN3 = list.files("D:/Recerca/Datasets/IFN/Sources/IFN3/MSAccess/",
#' pattern = "*[.]accdb", full.names = T)
#'
#' prov = 27 # Lugo
#'
#' treeData = readPiesMayoresIFN3(accessfilesIFN3[prov])
#'
#' IFNscenario(treeData, plotDataIFN, defaultPrescription,
#'             estratosIFN3[[prov]], defaultDemand[[prov]], numSteps = 3)
#' }
#'
IFNscenario<-function(treeData, plotData, prescription, strata, demand, nb = NULL, numSteps = 1,
                      plotDataDyn = NULL, demandOffset = NULL, stepOffset = 0,
                      as.CO2 = TRUE, sequence = FALSE, dead.table = FALSE, cut.table = FALSE,
                      verbose = FALSE, ...) {

  scen = NULL
  plotDataSuppl = NULL
  if(inherits(treeData, "data.frame")) {
    cat(paste("A. Data preparation...\n"))
    #Check variable type
    treeData$ID = as.character(treeData$ID)
    treeData$Species = as.character(treeData$Species)

    selDBH = is.na(treeData$DBH)
    if(sum(selDBH)>0) {
      cat(paste0("  Removing ", sum(selDBH), " records out of ", length(selDBH)," because of missing diameter.\n"))
      treeData = treeData[!selDBH,]
    }
    selDBH = (treeData$DBH==0)
    if(sum(selDBH)>0) {
      cat(paste0("  Removing ", sum(selDBH), " records out of ", length(selDBH)," because of zero diameter.\n"))
      treeData = treeData[!selDBH,]
    }
    selH = is.na(treeData$H)
    if(sum(selH)>0) {
      cat(paste0("  Removing ", sum(selH), " records out of ", length(selH)," because of missing height.\n"))
      treeData = treeData[!selH,]
    }
    selH = (treeData$H==0)
    if(sum(selH)>0) {
      cat(paste0("  Removing ", sum(selH), " records out of ", length(selH)," because of zero height.\n"))
      treeData = treeData[!selH,]
    }
    selSpecies = is.na(treeData$Species)
    if(sum(selSpecies)>0) {
      cat(paste0("  Removing ", sum(selSpecies), " records out of ",length(selSpecies)," because of missing species.\n"))
      treeData = treeData[!selSpecies,]
    }
    selSpecies = !(treeData$Species %in% biomass_species_match$ID_TAXON)
    if(sum(selSpecies)>0) {
      wc = unique(treeData$Species[selSpecies])
      cat(paste0("  Removing ", sum(selSpecies), " records out of ",length(selSpecies)," because of wrong species code(s): ",paste0(wc,collapse=","),".\n"))
      treeData = treeData[!selSpecies,]
    }

  }
  else if(inherits(treeData, "IFNscenario")) {
    cat(paste("A. Data loading from previous scenario...\n"))
    scen = treeData
    plotDataSuppl  = scen$plotDataSuppl
    treeData = scen$treeDataNew
    demandOffset = scen$demandOffset
  }

  #Check plotData
  selSWHC = is.na(plotData$SWHC)
  if(sum(selSWHC)>0) {
    cat(paste0("  Removing ", sum(selSWHC), " plots from plotData because of missing 'SWHC'.\n"))
    plotData = plotData[!selSWHC,]
  }
  selRad = is.na(plotData$Rad)
  if(sum(selRad)>0) {
    cat(paste0("  Removing ", sum(selRad), " plots from plotData because of missing 'Rad'.\n"))
    plotData = plotData[!selRad,]
  }
  selTemp = is.na(plotData$Temp)
  if(sum(selTemp)>0) {
    cat(paste0("  Removing ", sum(selTemp), " plots from plotData because of missing 'Temp'.\n"))
    plotData = plotData[!selTemp,]
  }
  selSlope = is.na(plotData$slope)
  if(sum(selSlope)>0) {
    cat(paste0("  Removing ", sum(selSlope), " plots from plotData because of missing 'slope'.\n"))
    plotData = plotData[!selSlope,]
  }
  selEstrato = is.na(plotData$Estrato)
  if(sum(selEstrato)>0) {
    cat(paste0("  Removing ", sum(selEstrato), " plots from plotData because of missing 'estrato'.\n"))
    plotData = plotData[!selEstrato,]
  }

  #Check plot codes
  plotIDs = unique(treeData$ID)
  sel = (plotIDs %in% row.names(plotData))
  nmis = sum(!sel)
  if(nmis>0) {
    cat(paste("  Removing tree data for", nmis, "plots not found in 'plotData'.\n"))
    plotIDs = plotIDs[sel]
    treeData = treeData[treeData$ID %in% plotIDs,]
  }
  cat(paste("  Subsetting", length(plotIDs), "plots from 'plotData'.\n"))
  plotData = plotData[row.names(plotData) %in% plotIDs,] # Subset records
  plotData = plotData[plotIDs,] # Sort plot data records
  if(!is.null(plotDataDyn)) {
    cat(paste("  Subsetting", length(plotIDs), "plots from 'plotDataDyn'.\n"))
    plotDataDyn = plotDataDyn[plotDataDyn$ID %in% plotIDs,] # Subset records
  }

  #If 'Area' is not defined in 'plotData', determine area (hectares) corresponding to each plot according to IFN strata definition
  nPar = table(plotData$Estrato)
  if(!("Area" %in% names(plotData))) {
    cat(paste("  Calculating plot area from strata.\n"))
    plotData$Area = NA
    for(i in 1:length(plotIDs)) {
      plotData[i, "Area"] = strata[plotData[i, "Estrato"],"Area"]/nPar[plotData[i, "Estrato"]]
    }
  } else {
    cat(paste("  Using plot area found in 'plotData'.\n"))
  }
  selArea = is.na(plotData$Area)
  if(sum(selArea)>0) {
    cat(paste("  Removing data for", sum(selArea), "plots with missing plot area.\n"))
    ids = row.names(plotData)[selArea]
    plotData = plotData[!selArea, ]
    treeData = treeData[!(treeData$ID %in% ids),]
    plotIDs = plotIDs[!(plotIDs %in% ids)]
  }

  #For all plots, determine dominant species (may be missing)
  BA = basalArea(treeData, bySpecies = T)
  plotData$DominantSpecies = colnames(BA)[apply(BA, 1, which.max)]
  if(!is.null(plotDataSuppl)) {
    cat(paste("  Replacing dominant species from input scenario.\n"))
    plotData[row.names(plotDataSuppl), "DominantSpecies"] = plotDataSuppl$DominantSpecies
  }

  #Determine prescription and demand rows for each forest plot
  if(!("PrescriptionRow" %in% names(plotData))) {
    cat(paste("  Initializing prescription row.\n"))
    plotData$PrescriptionRow = NA
    pr_s = strsplit(as.character(prescription$Espcodes), split=",")
    getPrescriptionRow<-function(sp) {
      for(i in 1:length(pr_s)) if(as.character(sp) %in% pr_s[[i]]) return(i)
      return(NA)
    }
    for(i in 1:length(plotIDs)) {
      plotData[i, "PrescriptionRow"] = getPrescriptionRow(plotData[i,"DominantSpecies"])
    }
  } else {
    cat(paste("  Using prescription row found in 'plotData'.\n"))
  }
  if(!is.null(plotDataSuppl)) {
    cat(paste("  Replacing prescription row from input scenario.\n"))
    plotData[row.names(plotDataSuppl), "PrescriptionRow"] = plotDataSuppl$PrescriptionRow
  }


  if(!("DemandRow" %in% names(plotData))) {
    cat(paste("  Initializing demand row.\n"))
    plotData$DemandRow = NA
    dem_s = strsplit(as.character(demand$Species), split=",")
    getDemandRow<-function(sp) {
      for(i in 1:length(dem_s)) if(as.character(sp) %in% dem_s[[i]]) return(i)
      return(NA)
    }
    for(i in 1:length(plotIDs)) {
      plotData[i, "DemandRow"] = getDemandRow(plotData[i,"DominantSpecies"])
    }
  } else {
    cat(paste("  Using demand row found in 'plotData'.\n"))
  }
  if(!is.null(plotDataSuppl)) {
    cat(paste("  Replacing demand row from input scenario.\n"))
    plotData[row.names(plotDataSuppl), "DemandRow"] = plotDataSuppl$DemandRow
  }
  if(!("finalPreviousStage" %in% names(plotData))) {
    cat(paste("  Initializing final previous stage.\n"))
    plotData$finalPreviousStage = NA
    for(i in 1:length(plotIDs)) {
      if(plotData[i, "Managed"]) plotData[i, "finalPreviousStage"] = prescription[plotData$PrescriptionRow[i], "finalPreviousStage"]
    }
  } else {
    cat(paste("  Using final previous stage found in 'plotData'.\n"))
  }
  if(!is.null(plotDataSuppl)) {
    cat(paste("  Replacing previous stage from input scenario.\n"))
    plotData[row.names(plotDataSuppl), "finalPreviousStage"] = plotDataSuppl$finalPreviousStage
  }


  #Initial demand offset for all species
  if(!is.null(demandOffset))  {
    cat("  Setting initial demand offset to input.\n")
    if(length(demandOffset)!= nrow(demand)) stop("Wrong demand offset length.")
    TenYrDemandOffset = demandOffset
  } else {
    cat("  Setting initial demand offset to zero.\n")
    TenYrDemandOffset = rep(0, nrow(demand))
  }
  manTurnOff = 0
  for(i in 1:length(plotIDs)) {
    if(is.na(plotData[i, "PrescriptionRow"]) || is.na(plotData[i, "DemandRow"]))  {
      if(plotData[i, "Managed"]) {
        plotData[i, "Managed"] = FALSE
        manTurnOff = manTurnOff+1
      }
    }
  }
  if(manTurnOff>0) cat(paste("  Management for", manTurnOff, "plots has been turned off because of missing prescription/demand.\n"))

  # print(table(plotData$DominantSpecies, plotData$Managed))
  # print(table(plotData$DominantSpecies, plotData$finalPreviousStage))
  #Divide treeData between managed and non-managed
  managedIDs = plotIDs[plotData$Managed]
  nonmanagedIDs = plotIDs[!plotData$Managed]
  cat(paste0("  ",length(managedIDs), " managed plots / ", length(nonmanagedIDs), " non-managed plots / ", nrow(treeData)," initial tree records.\n"))

  prescription$plantingSpecies = NA
  prescription$plantingDBH = 5
  prescription$plantingHeight = 5
  prescription$plantingDensity = 200

  #Init output matrices
  plotActions = matrix(NA, length(plotIDs), numSteps)
  rownames(plotActions) = plotIDs
  colnames(plotActions) = 1:numSteps


  extractedDemand = matrix(0, nrow(demand)+1, numSteps)
  rownames(extractedDemand) = c(demand$Name,"Total")
  colnames(extractedDemand) = 1:numSteps

  managedArea = matrix(0, nrow(demand)+1, numSteps)
  rownames(managedArea) = c(demand$Name,"Total")
  colnames(managedArea) = 1:numSteps

  domspp = as.character(sort(as.numeric(unique(plotData$DominantSpecies))))
  spp = as.character(sort(as.numeric(unique(c(treeData$Species, domspp))))) ## Add dominant species if not currently present


  #Standing volume
  standingVolumeSpp = matrix(0, length(spp), numSteps+1)
  rownames(standingVolumeSpp) = spp
  colnames(standingVolumeSpp) = 0:numSteps

  standingVolumePlot = matrix(0, length(plotIDs), numSteps+1)
  rownames(standingVolumePlot) = plotIDs
  colnames(standingVolumePlot) = 0:numSteps

  standingVolumeDomSpp = matrix(0, length(domspp), numSteps+1)
  rownames(standingVolumeDomSpp) = domspp
  colnames(standingVolumeDomSpp) = 0:numSteps


  #Dead volume
  deadVolumeSpp = matrix(0, length(spp), numSteps)
  rownames(deadVolumeSpp) = spp
  colnames(deadVolumeSpp) = 1:numSteps

  deadVolumePlot = matrix(0, length(plotIDs), numSteps)
  rownames(deadVolumePlot) = plotIDs
  colnames(deadVolumePlot) = 1:numSteps

  deadVolumeDomSpp = matrix(0, length(domspp), numSteps)
  rownames(deadVolumeDomSpp) = domspp
  colnames(deadVolumeDomSpp) = 1:numSteps


  #Extracted volume
  extractedVolumeSpp = matrix(0, length(spp), numSteps)
  rownames(extractedVolumeSpp) = spp
  colnames(extractedVolumeSpp) = 1:numSteps

  extractedVolumePlot = matrix(0, length(plotIDs), numSteps)
  rownames(extractedVolumePlot) = plotIDs
  colnames(extractedVolumePlot) = 1:numSteps

  extractedVolumeDomSpp = matrix(0, length(domspp), numSteps)
  rownames(extractedVolumeDomSpp) = domspp
  colnames(extractedVolumeDomSpp) = 1:numSteps

  #Standing density
  standingDensitySpp = matrix(0, length(spp), numSteps+1)
  rownames(standingDensitySpp) = spp
  colnames(standingDensitySpp) = 0:numSteps

  standingDensityPlot = matrix(0, length(plotIDs), numSteps+1)
  rownames(standingDensityPlot) = plotIDs
  colnames(standingDensityPlot) = 0:numSteps

  standingDensityDomSpp = matrix(0, length(domspp), numSteps+1)
  rownames(standingDensityDomSpp) = domspp
  colnames(standingDensityDomSpp) = 0:numSteps

  #Dead density
  deadDensitySpp = matrix(0, length(spp), numSteps)
  rownames(deadDensitySpp) = spp
  colnames(deadDensitySpp) = 1:numSteps

  deadDensityPlot = matrix(0, length(plotIDs), numSteps)
  rownames(deadDensityPlot) = plotIDs
  colnames(deadDensityPlot) = 1:numSteps

  deadDensityDomSpp = matrix(0, length(domspp), numSteps)
  rownames(deadDensityDomSpp) = domspp
  colnames(deadDensityDomSpp) = 1:numSteps


  #Extracted density
  extractedDensitySpp = matrix(0, length(spp), numSteps)
  rownames(extractedDensitySpp) = spp
  colnames(extractedDensitySpp) = 1:numSteps

  extractedDensityPlot = matrix(0, length(plotIDs), numSteps)
  rownames(extractedDensityPlot) = plotIDs
  colnames(extractedDensityPlot) = 1:numSteps

  extractedDensityDomSpp = matrix(0, length(domspp), numSteps)
  rownames(extractedDensityDomSpp) = domspp
  colnames(extractedDensityDomSpp) = 1:numSteps


  #Standing root biomass
  standingRootBiomassSpp = matrix(0, length(spp), numSteps+1)
  rownames(standingRootBiomassSpp) = spp
  colnames(standingRootBiomassSpp) = 0:numSteps

  standingRootBiomassPlot = matrix(0, length(plotIDs), numSteps+1)
  rownames(standingRootBiomassPlot) = plotIDs
  colnames(standingRootBiomassPlot) = 0:numSteps

  standingRootBiomassDomSpp = matrix(0, length(domspp), numSteps+1)
  rownames(standingRootBiomassDomSpp) = domspp
  colnames(standingRootBiomassDomSpp) = 0:numSteps

  #Standing stem biomass
  standingStemBiomassSpp = matrix(0, length(spp), numSteps+1)
  rownames(standingStemBiomassSpp) = spp
  colnames(standingStemBiomassSpp) = 0:numSteps

  standingStemBiomassPlot = matrix(0, length(plotIDs), numSteps+1)
  rownames(standingStemBiomassPlot) = plotIDs
  colnames(standingStemBiomassPlot) = 0:numSteps

  standingStemBiomassDomSpp = matrix(0, length(domspp), numSteps+1)
  rownames(standingStemBiomassDomSpp) = domspp
  colnames(standingStemBiomassDomSpp) = 0:numSteps

  #Standing branch biomass
  standingBranchBiomassSpp = matrix(0, length(spp), numSteps+1)
  rownames(standingBranchBiomassSpp) = spp
  colnames(standingBranchBiomassSpp) = 0:numSteps

  standingBranchBiomassPlot = matrix(0, length(plotIDs), numSteps+1)
  rownames(standingBranchBiomassPlot) = plotIDs
  colnames(standingBranchBiomassPlot) = 0:numSteps

  standingBranchBiomassDomSpp = matrix(0, length(domspp), numSteps+1)
  rownames(standingBranchBiomassDomSpp) = domspp
  colnames(standingBranchBiomassDomSpp) = 0:numSteps

  #Dead root biomass
  deadRootBiomassSpp = matrix(0, length(spp), numSteps)
  rownames(deadRootBiomassSpp) = spp
  colnames(deadRootBiomassSpp) = 1:numSteps

  deadRootBiomassPlot = matrix(0, length(plotIDs), numSteps)
  rownames(deadRootBiomassPlot) = plotIDs
  colnames(deadRootBiomassPlot) = 1:numSteps

  deadRootBiomassDomSpp = matrix(0, length(domspp), numSteps)
  rownames(deadRootBiomassDomSpp) = domspp
  colnames(deadRootBiomassDomSpp) = 1:numSteps

  #Dead stem biomass
  deadStemBiomassSpp = matrix(0, length(spp), numSteps)
  rownames(deadStemBiomassSpp) = spp
  colnames(deadStemBiomassSpp) = 1:numSteps

  deadStemBiomassPlot = matrix(0, length(plotIDs), numSteps)
  rownames(deadStemBiomassPlot) = plotIDs
  colnames(deadStemBiomassPlot) = 1:numSteps

  deadStemBiomassDomSpp = matrix(0, length(domspp), numSteps)
  rownames(deadStemBiomassDomSpp) = domspp
  colnames(deadStemBiomassDomSpp) = 1:numSteps


  #Dead branch biomass
  deadBranchBiomassSpp = matrix(0, length(spp), numSteps)
  rownames(deadBranchBiomassSpp) = spp
  colnames(deadBranchBiomassSpp) = 1:numSteps

  deadBranchBiomassPlot = matrix(0, length(plotIDs), numSteps)
  rownames(deadBranchBiomassPlot) = plotIDs
  colnames(deadBranchBiomassPlot) = 1:numSteps

  deadBranchBiomassDomSpp = matrix(0, length(domspp), numSteps)
  rownames(deadBranchBiomassDomSpp) = domspp
  colnames(deadBranchBiomassDomSpp) = 1:numSteps


  #Extracted root biomass
  extractedRootBiomassSpp = matrix(0, length(spp), numSteps)
  rownames(extractedRootBiomassSpp) = spp
  colnames(extractedRootBiomassSpp) = 1:numSteps

  extractedRootBiomassPlot = matrix(0, length(plotIDs), numSteps)
  rownames(extractedRootBiomassPlot) = plotIDs
  colnames(extractedRootBiomassPlot) = 1:numSteps

  extractedRootBiomassDomSpp = matrix(0, length(domspp), numSteps)
  rownames(extractedRootBiomassDomSpp) = domspp
  colnames(extractedRootBiomassDomSpp) = 1:numSteps

  #Extracted stem biomass
  extractedStemBiomassSpp = matrix(0, length(spp), numSteps)
  rownames(extractedStemBiomassSpp) = spp
  colnames(extractedStemBiomassSpp) = 1:numSteps

  extractedStemBiomassPlot = matrix(0, length(plotIDs), numSteps)
  rownames(extractedStemBiomassPlot) = plotIDs
  colnames(extractedStemBiomassPlot) = 1:numSteps

  extractedStemBiomassDomSpp = matrix(0, length(domspp), numSteps)
  rownames(extractedStemBiomassDomSpp) = domspp
  colnames(extractedStemBiomassDomSpp) = 1:numSteps

  #Extracted branch biomass
  extractedBranchBiomassSpp = matrix(0, length(spp), numSteps)
  rownames(extractedBranchBiomassSpp) = spp
  colnames(extractedBranchBiomassSpp) = 1:numSteps

  extractedBranchBiomassPlot = matrix(0, length(plotIDs), numSteps)
  rownames(extractedBranchBiomassPlot) = plotIDs
  colnames(extractedBranchBiomassPlot) = 1:numSteps

  extractedBranchBiomassDomSpp = matrix(0, length(domspp), numSteps)
  rownames(extractedBranchBiomassDomSpp) = domspp
  colnames(extractedBranchBiomassDomSpp) = 1:numSteps



  #Copying area data
  if(!("Area" %in% names(treeData))) treeData$Area = NA
  if(!("Provincia" %in% names(treeData))) treeData$Provincia = NA
  cat("  Copying area data to tree table.\n")
  for(i in 1:nrow(plotData)) {
    id = row.names(plotData)[i]
    treeData$Provincia[treeData$ID==id] = plotData[id, "Province"]
    treeData$Area[treeData$ID==id] = plotData[id, "Area"]
  }

  #Initial height
  cat("  Calculating initial height.\n")
  treeData$H = IFNheight(treeData, plotDataIFN, verbose = verbose)

  #Initial volume
  cat("  Calculating initial volume.\n")
  initialVolume = IFNvolume(treeData, verbose = verbose)
  volspp = tapply(initialVolume$VCC*treeData$Area, initialVolume$Species, sum, na.rm=T)
  standingVolumeSpp[names(volspp), 1] = as.numeric(volspp)
  volplot = tapply(initialVolume$VCC*treeData$Area, initialVolume$ID, sum, na.rm=T)
  standingVolumePlot[names(volplot), 1] = as.numeric(volplot)
  voldomspp = tapply(volplot, plotData[names(volplot), "DominantSpecies"], sum, na.rm=T)
  standingVolumeDomSpp[names(voldomspp), 1] = as.numeric(voldomspp)

  #Initial density
  cat("  Calculating initial density.\n")
  densspp = tapply(treeData$N*treeData$Area, treeData$Species, sum, na.rm=T)
  standingDensitySpp[names(volspp), 1] = as.numeric(densspp)
  densplot = tapply(treeData$N*treeData$Area, treeData$ID, sum, na.rm=T)
  standingDensityPlot[names(densplot), 1] = as.numeric(densplot)
  densdomspp = tapply(densplot, plotData[names(densplot), "DominantSpecies"], sum, na.rm=T)
  standingDensityDomSpp[names(densdomspp), 1] = as.numeric(densdomspp)


  #Initial biomass (Mg)
  cat("  Calculating initial biomass.\n")
  initialBiomass = IFNbiomass(treeData, as.CO2 = as.CO2, verbose = verbose)
  rootspp = 0.001*tapply(initialBiomass$Roots*treeData$Area, initialBiomass$Species, sum, na.rm=T)
  stemspp = 0.001*tapply(initialBiomass$Stem*treeData$Area, initialBiomass$Species, sum, na.rm=T)
  branchspp = 0.001*tapply(initialBiomass$Branches*treeData$Area, initialBiomass$Species, sum, na.rm=T)
  standingRootBiomassSpp[names(rootspp), 1] = as.numeric(rootspp)
  standingStemBiomassSpp[names(stemspp), 1] = as.numeric(stemspp)
  standingBranchBiomassSpp[names(branchspp), 1] = as.numeric(branchspp)
  rootplot = 0.001*tapply(initialBiomass$Roots*treeData$Area, initialBiomass$ID, sum, na.rm=T)
  stemplot = 0.001*tapply(initialBiomass$Stem*treeData$Area, initialBiomass$ID, sum, na.rm=T)
  branchplot = 0.001*tapply(initialBiomass$Branches*treeData$Area, initialBiomass$ID, sum, na.rm=T)
  standingRootBiomassPlot[names(rootplot), 1] = as.numeric(rootplot)
  standingStemBiomassPlot[names(stemplot), 1] = as.numeric(stemplot)
  standingBranchBiomassPlot[names(branchplot), 1] = as.numeric(branchplot)
  rootdomspp = tapply(rootplot, plotData[names(rootplot), "DominantSpecies"], sum, na.rm=T)
  stemdomspp = tapply(stemplot, plotData[names(stemplot), "DominantSpecies"], sum, na.rm=T)
  branchdomspp = tapply(branchplot, plotData[names(branchplot), "DominantSpecies"], sum, na.rm=T)
  standingRootBiomassDomSpp[names(rootdomspp), 1] = as.numeric(rootdomspp)
  standingStemBiomassDomSpp[names(stemdomspp), 1] = as.numeric(stemdomspp)
  standingBranchBiomassDomSpp[names(branchdomspp), 1] = as.numeric(branchdomspp)

  cat("\n") ## End of data preparation

  cat(paste("B. Simulations...\n"))

  # Store initial state as sequence
  if(sequence) {
    treeDataSeq = data.frame(Step = rep(0, nrow(treeData)),
                             ID = treeData$ID,
                             Species = treeData$Species,
                             DBH = treeData$DBH,
                             N = treeData$N,
                             H = treeData$H,
                             stringsAsFactors = F)
  }
  if(dead.table) treeDataDeadSeq = treeDataSeq[numeric(0),]
  if(cut.table) treeDataCutSeq = treeDataSeq[numeric(0),]


  for(step in 1:numSteps) { # START OF STEP LOOP
    cat(paste0(" *** Step ", step, "/", numSteps," ***\n"))

    #B.00 calculate NspNmatrix
    NspNmatrix = NULL
    if(!is.null(nb)) {
      cat(paste0(" Calculating Nsp/N matrix.\n"))
      Nmatrix = tapply(treeData[,"N"], list(treeData[,"ID"], treeData[,"Species"]),
                       FUN = sum, na.rm=T)
      NspNmatrix = sweep(Nmatrix,1, rowSums(Nmatrix, na.rm=T), "/")
      rm(Nmatrix)
    }

    #B.0 Replace dynamic plot data
    if(!is.null(plotDataDyn)) {
      cat(paste0(" Drawing step [",(step + stepOffset), "] data from plotDataDyn.\n"))
      plotDataDynStep = plotDataDyn[plotDataDyn$Step == (step + stepOffset),, drop=FALSE]
      if(nrow(plotDataDynStep)==0) stop(paste0("Step data [", (step + stepOffset), "] not defined in plotDataDyn."))
      row.names(plotDataDynStep) = plotDataDynStep$ID
    }

    #B. Process dynamics of non-managed plots
    treeDataNonManagedNew = NULL
    treeDataNonManagedDead = NULL
    if(length(nonmanagedIDs)>0) {
      cat(paste(" Processing", length(nonmanagedIDs),"non-managed plots:"))
      treeDataNonManaged = treeData[treeData$ID %in% nonmanagedIDs, , drop = FALSE]
      plotDataNonManaged = plotData[nonmanagedIDs,, drop = FALSE]
      if(!is.null(plotDataDyn)) {
        plotDataNonManaged$Temp = plotDataDynStep[nonmanagedIDs,"Temp"]
        plotDataNonManaged$Prec = plotDataDynStep[nonmanagedIDs,"Prec"]
        plotDataNonManaged$PET = plotDataDynStep[nonmanagedIDs,"PET"]
      }
      cat(paste0(" dynamics"))

      #B.1 Calculate dynamics
      dyn = IFNdynamics(treeDataNonManaged, plotDataNonManaged,
                        nb = nb, NspNmatrix = NspNmatrix,
                        numSteps = 1,
                        dead.table  = TRUE, ...)
      treeDataNonManagedNew = dyn$out[,c("ID", "Species", "Name", "DBH", "N", "H"), drop = FALSE]
      treeDataNonManagedDead = dyn$dead[, c("ID", "Species", "Name", "DBH", "N", "H"), drop = FALSE]

      #B.2 Add provincia/area
      treeDataNonManagedNew$Provincia = NA
      treeDataNonManagedNew$Area = NA
      for(i in 1:nrow(plotDataNonManaged)) {
        id = row.names(plotDataNonManaged)[i]
        treeDataNonManagedNew$Provincia[treeDataNonManagedNew$ID==id] = plotDataNonManaged[id, "Province"]
        treeDataNonManagedNew$Area[treeDataNonManagedNew$ID==id] = plotDataNonManaged[id, "Area"]
        treeDataNonManagedDead$Provincia[treeDataNonManagedDead$ID==id] = plotDataNonManaged[id, "Province"]
        treeDataNonManagedDead$Area[treeDataNonManagedDead$ID==id] = plotDataNonManaged[id, "Area"]
      }

      #B.3 Calculate standing/dead volume and density
      cat(paste(", volume"))
      nmvolume = IFNvolume(treeDataNonManagedNew)
      volspp = tapply(nmvolume$VCC*treeDataNonManagedNew$Area, nmvolume$Species, sum, na.rm=T)
      standingVolumeSpp[names(volspp), step+1] = standingVolumeSpp[names(volspp), step+1] + as.numeric(volspp)
      volplot = tapply(nmvolume$VCC*treeDataNonManagedNew$Area, nmvolume$ID, sum, na.rm=T)
      standingVolumePlot[names(volplot), step+1] = standingVolumePlot[names(volplot), step+1] + as.numeric(volplot)
      voldomspp = tapply(volplot, plotData[names(volplot), "DominantSpecies"], sum, na.rm=T)
      standingVolumeDomSpp[names(voldomspp), step+1] = standingVolumeDomSpp[names(voldomspp), step+1] + as.numeric(voldomspp)
      nmdeadvolume = IFNvolume(treeDataNonManagedDead)
      volspp = tapply(nmdeadvolume$VCC*treeDataNonManagedDead$Area, nmdeadvolume$Species, sum, na.rm=T)
      deadVolumeSpp[names(volspp), step] = deadVolumeSpp[names(volspp), step] + as.numeric(volspp)
      volplot = tapply(nmdeadvolume$VCC*treeDataNonManagedDead$Area, nmdeadvolume$ID, sum, na.rm=T)
      deadVolumePlot[names(volplot), step] = deadVolumePlot[names(volplot), step] + as.numeric(volplot)
      voldomspp = tapply(volplot, plotData[names(volplot), "DominantSpecies"], sum, na.rm=T)
      deadVolumeDomSpp[names(voldomspp), step] = deadVolumeDomSpp[names(voldomspp), step] + as.numeric(voldomspp)

      cat(paste(", density"))
      densspp = tapply(treeDataNonManagedNew$N*treeDataNonManagedNew$Area, nmvolume$Species, sum, na.rm=T)
      standingDensitySpp[names(densspp), step+1] =  standingDensitySpp[names(densspp), step+1] + as.numeric(densspp)
      densplot = tapply(treeDataNonManagedNew$N*treeDataNonManagedNew$Area, nmvolume$ID, sum, na.rm=T)
      standingDensityPlot[names(densplot), step+1] =  standingDensityPlot[names(densplot), step+1] + as.numeric(densplot)
      densdomspp = tapply(densplot, plotData[names(densplot), "DominantSpecies"], sum, na.rm=T)
      standingDensityDomSpp[names(densdomspp), step+1] =  standingDensityDomSpp[names(densdomspp), step+1] + as.numeric(densdomspp)
      densspp = tapply(treeDataNonManagedDead$N*treeDataNonManagedDead$Area, nmdeadvolume$Species, sum, na.rm=T)
      deadDensitySpp[names(densspp), step] = deadDensitySpp[names(densspp), step] +as.numeric(densspp)
      densplot = tapply(treeDataNonManagedDead$N*treeDataNonManagedDead$Area, nmdeadvolume$ID, sum, na.rm=T)
      deadDensityPlot[names(densplot), step] =  deadDensityPlot[names(densplot), step] + as.numeric(densplot)
      densdomspp = tapply(densplot, plotData[names(densplot), "DominantSpecies"], sum, na.rm=T)
      deadDensityDomSpp[names(densdomspp), step] =  deadDensityDomSpp[names(densdomspp), step] + as.numeric(densdomspp)


      #B.4 Calculate standing/dead biomass
      cat(paste(", biomass"))
      nmlivebiomass = IFNbiomass(treeDataNonManagedNew, as.CO2 = as.CO2)
      biomassspp = 0.001*tapply(nmlivebiomass$Roots*treeDataNonManagedNew$Area, nmlivebiomass$Species, sum, na.rm=T)
      standingRootBiomassSpp[names(biomassspp), step+1] = standingRootBiomassSpp[names(biomassspp), step+1] + as.numeric(biomassspp)
      biomassplot = 0.001*tapply(nmlivebiomass$Roots*treeDataNonManagedNew$Area, nmlivebiomass$ID, sum, na.rm=T)
      standingRootBiomassPlot[names(biomassplot), step+1] =  standingRootBiomassPlot[names(biomassplot), step+1] + as.numeric(biomassplot)
      biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
      standingRootBiomassDomSpp[names(biomassdomspp), step+1] =  standingRootBiomassDomSpp[names(biomassdomspp), step+1] + as.numeric(biomassdomspp)

      biomassspp = 0.001*tapply(nmlivebiomass$Stem*treeDataNonManagedNew$Area, nmlivebiomass$Species, sum, na.rm=T)
      standingStemBiomassSpp[names(biomassspp), step+1] = standingStemBiomassSpp[names(biomassspp), step+1] + as.numeric(biomassspp)
      biomassplot = 0.001*tapply(nmlivebiomass$Stem*treeDataNonManagedNew$Area, nmlivebiomass$ID, sum, na.rm=T)
      standingStemBiomassPlot[names(biomassplot), step+1] =  standingStemBiomassPlot[names(biomassplot), step+1] + as.numeric(biomassplot)
      biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
      standingStemBiomassDomSpp[names(biomassdomspp), step+1] =  standingStemBiomassDomSpp[names(biomassdomspp), step+1] + as.numeric(biomassdomspp)

      biomassspp = 0.001*tapply(nmlivebiomass$Branches*treeDataNonManagedNew$Area, nmlivebiomass$Species, sum, na.rm=T)
      standingBranchBiomassSpp[names(biomassspp), step+1] = standingBranchBiomassSpp[names(biomassspp), step+1] + as.numeric(biomassspp)
      biomassplot = 0.001*tapply(nmlivebiomass$Branches*treeDataNonManagedNew$Area, nmlivebiomass$ID, sum, na.rm=T)
      standingBranchBiomassPlot[names(biomassplot), step+1] =  standingBranchBiomassPlot[names(biomassplot), step+1] + as.numeric(biomassplot)
      biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
      standingBranchBiomassDomSpp[names(biomassdomspp), step+1] =  standingBranchBiomassDomSpp[names(biomassdomspp), step+1] + as.numeric(biomassdomspp)

      nmdeadbiomass = IFNbiomass(treeDataNonManagedDead, as.CO2 = as.CO2)
      biomassspp = 0.001*tapply(nmdeadbiomass$Roots*treeDataNonManagedDead$Area, nmdeadbiomass$Species, sum, na.rm=T)
      deadRootBiomassSpp[names(biomassspp), step] = deadRootBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
      biomassplot = 0.001*tapply(nmdeadbiomass$Roots*treeDataNonManagedDead$Area, nmdeadbiomass$ID, sum, na.rm=T)
      deadRootBiomassPlot[names(biomassplot), step] =  deadRootBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
      biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
      deadRootBiomassDomSpp[names(biomassdomspp), step] =  deadRootBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)

      biomassspp = 0.001*tapply(nmdeadbiomass$Stem*treeDataNonManagedDead$Area, nmdeadbiomass$Species, sum, na.rm=T)
      deadStemBiomassSpp[names(biomassspp), step] = deadStemBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
      biomassplot = 0.001*tapply(nmdeadbiomass$Stem*treeDataNonManagedDead$Area, nmdeadbiomass$ID, sum, na.rm=T)
      deadStemBiomassPlot[names(biomassplot), step] =  deadStemBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
      biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
      deadStemBiomassDomSpp[names(biomassdomspp), step] =  deadStemBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)

      biomassspp = 0.001*tapply(nmdeadbiomass$Branches*treeDataNonManagedDead$Area, nmdeadbiomass$Species, sum, na.rm=T)
      deadBranchBiomassSpp[names(biomassspp), step] = deadBranchBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
      biomassplot = 0.001*tapply(nmdeadbiomass$Branches*treeDataNonManagedDead$Area, nmdeadbiomass$ID, sum, na.rm=T)
      deadBranchBiomassPlot[names(biomassplot), step] =  deadBranchBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
      biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
      deadBranchBiomassDomSpp[names(biomassdomspp), step] =  deadBranchBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)

      cat(".\n") ## End of non-managed plots (B)
    }
    # print(head(treeDataNonManagedNew))
    # print(head(treeDataNonManagedDead))

    # stop()
    #C. Process dynamics of managed plots
    cat(paste(" Processing", length(managedIDs),"managed plots:\n"))
    plotDataManaged  = plotData[managedIDs,]
    treeDataManaged = treeData[treeData$ID %in% managedIDs, ]
    treeDataManagedNew = NULL
    treeDataManagedDead = NULL

    for(r in 1:nrow(demand)) {
      sel = (plotDataManaged$DemandRow == r)
      if(sum(sel)>0) {
        cat(paste0("  [",r,"/",nrow(demand),"] - ", sum(sel), " plots dominated by ", demand$Name[r], ":"))
        #C.0 Subset plots dominaged by species
        plotDataManagedSel = plotDataManaged[sel,]
        treeDataManagedSel = treeDataManaged[treeDataManaged$ID %in% row.names(plotDataManagedSel), ]
        prescriptionDataManagedSel = prescription[plotDataManagedSel$PrescriptionRow,-c(1,2)]
        # Copy final previous stage (from previous step) from plot data
        prescriptionDataManagedSel$finalPreviousStage = plotDataManaged$finalPreviousStage[sel]
        row.names(prescriptionDataManagedSel) = row.names(plotDataManagedSel)

        #C.2 Replace dynamic plot data
        if(!is.null(plotDataDyn)) {
          manSelIDs = row.names(plotDataManagedSel)
          plotDataManagedSel$Temp = plotDataDynStep[manSelIDs,"Temp"]
          plotDataManagedSel$Prec = plotDataDynStep[manSelIDs,"Prec"]
          plotDataManagedSel$PET = plotDataDynStep[manSelIDs,"PET"]
        }
        #C.3 Calculate Hart-Becking
        HB = hartBeckingIndex(treeDataManagedSel)
        plotDataManagedSel$HB  = HB[row.names(plotDataManagedSel)]
        #C.4 Set demand
        TenYrDemand = TenYrDemandOffset[r]+demand$AnnualDemand[r]*10


        managedAreaSel = 0
        Extracted = 0
        treeDataOut = NULL #Stores the result of management
        treeDataCut = NULL #Stores the result of management
        treated = 0
        if(!verbose) cat(paste0(" management"))
        regularPlots = prescriptionDataManagedSel$type=="regular"
        if(!verbose) cat(paste0(" [regular ", sum(regularPlots), " irregular ", sum(!regularPlots),"]"))
        # if(!verbose) cat(paste0(" management [Dem. ", round(TenYrDemand),"m3"))
        #C.5 Process first plots with regular models and previous stage > 0
        finalCutPlots = (regularPlots) & (prescriptionDataManagedSel$finalPreviousStage>0)
        if(sum(finalCutPlots)>0) {
          if(!verbose) cat(paste0(", final cuts ", sum(finalCutPlots)," "))
          plotDataFC = plotDataManagedSel[finalCutPlots, ]
          prescriptionDataFC = prescriptionDataManagedSel[finalCutPlots,]
          prescriptionDataFC$plantingSpecies = plotDataFC$DominantSpecies
          treeDataFC = treeDataManagedSel[treeDataManagedSel$ID %in% row.names(plotDataManagedSel)[finalCutPlots], ]
          if(nrow(treeDataFC)>0) {
            treeDataFCPost = IFNmanagement(treeDataFC, prescriptionDataFC, verbose = verbose)
            #store action
            action = treeDataFCPost$action
            plotActions[names(action), step] = action

            #add area to tree data cut
            treeDataFCPost$treeDataCut$Area = NA
            for(i in 1:nrow(plotDataFC)) {
              id = row.names(plotDataFC)[i]
              treeDataFCPost$treeDataCut$Area[treeDataFCPost$treeDataCut$ID==id] = plotDataFC[id, "Area"]
            }

            #Calculate extracted volume
            FCvolume = IFNvolume(treeDataFCPost$treeDataCut)
            VCCplot = tapply(FCvolume$VCC*treeDataFCPost$treeDataCut$Area, FCvolume$ID, "sum", na.rm=TRUE) #Remove NA (TODO check volume equations)

            #Update extractedVolumeSpp
            volspp = tapply(FCvolume$VCC*treeDataFCPost$treeDataCut$Area, FCvolume$Species, sum, na.rm=T)
            extractedVolumeSpp[names(volspp), step] = extractedVolumeSpp[names(volspp), step] + as.numeric(volspp)
            volplot = tapply(FCvolume$VCC*treeDataFCPost$treeDataCut$Area, FCvolume$ID, sum, na.rm=T)
            extractedVolumePlot[names(volplot), step] = extractedVolumePlot[names(volplot), step] + as.numeric(volplot)
            voldomspp = tapply(volplot, plotData[names(volplot), "DominantSpecies"], sum, na.rm=T)
            extractedVolumeDomSpp[names(voldomspp), step] = extractedVolumeDomSpp[names(voldomspp), step] + as.numeric(voldomspp)

            #Update extractedDensitySpp
            densspp= tapply(treeDataFCPost$treeDataCut$N*treeDataFCPost$treeDataCut$Area, FCvolume$Species, sum, na.rm=T)
            extractedDensitySpp[names(densspp), step] = extractedDensitySpp[names(densspp), step] + as.numeric(densspp)
            densplot= tapply(treeDataFCPost$treeDataCut$N*treeDataFCPost$treeDataCut$Area, FCvolume$ID, sum, na.rm=T)
            extractedDensityPlot[names(densplot), step] = extractedDensityPlot[names(densplot), step] + as.numeric(densplot)
            densdomspp = tapply(densplot, plotData[names(densplot), "DominantSpecies"], sum, na.rm=T)
            extractedDensityDomSpp[names(densdomspp), step] = extractedDensityDomSpp[names(densdomspp), step] + as.numeric(densdomspp)

            #Update extracted biomass
            FCbiomass = IFNbiomass(treeDataFCPost$treeDataCut, as.CO2 = as.CO2)
            biomassspp = 0.001*tapply(FCbiomass$Roots*treeDataFCPost$treeDataCut$Area, FCbiomass$Species, sum, na.rm=T)
            extractedRootBiomassSpp[names(biomassspp), step] = extractedRootBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
            biomassplot = 0.001*tapply(FCbiomass$Roots*treeDataFCPost$treeDataCut$Area, FCbiomass$ID, sum, na.rm=T)
            extractedRootBiomassPlot[names(biomassplot), step] = extractedRootBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
            biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
            extractedRootBiomassDomSpp[names(biomassdomspp), step] = extractedRootBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)

            biomassspp = 0.001*tapply(FCbiomass$Stem*treeDataFCPost$treeDataCut$Area, FCbiomass$Species, sum, na.rm=T)
            extractedStemBiomassSpp[names(biomassspp), step] = extractedStemBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
            biomassplot = 0.001*tapply(FCbiomass$Stem*treeDataFCPost$treeDataCut$Area, FCbiomass$ID, sum, na.rm=T)
            extractedStemBiomassPlot[names(biomassplot), step] = extractedStemBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
            biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
            extractedStemBiomassDomSpp[names(biomassdomspp), step] = extractedStemBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)

            biomassspp = 0.001*tapply(FCbiomass$Branches*treeDataFCPost$treeDataCut$Area, FCbiomass$Species, sum, na.rm=T)
            extractedBranchBiomassSpp[names(biomassspp), step] = extractedBranchBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
            biomassplot = 0.001*tapply(FCbiomass$Branches*treeDataFCPost$treeDataCut$Area, FCbiomass$ID, sum, na.rm=T)
            extractedBranchBiomassPlot[names(biomassplot), step] = extractedBranchBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
            biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
            extractedBranchBiomassDomSpp[names(biomassdomspp), step] = extractedBranchBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)

            #Add to extracted
            Extracted = Extracted + sum(VCCplot, na.rm=T)
            treated = treated + sum(finalCutPlots)
            #Add to managed area
            managedAreaSel = managedAreaSel + sum(plotDataFC$Area)

            #save final previous stage for next step
            prescriptionDataManagedSel$finalPreviousStage[finalCutPlots] = treeDataFCPost$nextPrescription$finalPreviousStage
            plotDataManagedSel$finalPreviousStage[finalCutPlots] = prescriptionDataManagedSel$finalPreviousStage[finalCutPlots]

            #Store result
            if(!is.null(treeDataOut)) {
              treeDataOut = rbind(treeDataOut, treeDataFCPost$treeDataOut)
            } else {
              treeDataOut = treeDataFCPost$treeDataOut
            }
            if(!is.null(treeDataCut)) {
              treeDataCut = rbind(treeDataCut, treeDataFCPost$treeDataCut)
            } else {
              treeDataCut = treeDataFCPost$treeDataCut
            }
          }
        }
        if(sum(!finalCutPlots)>0) {
          plotDataTh = plotDataManagedSel[!finalCutPlots, ]
          prescriptionDataTh = prescriptionDataManagedSel[!finalCutPlots,]
          treeDataTh = treeDataManagedSel[treeDataManagedSel$ID %in% row.names(plotDataManagedSel)[!finalCutPlots], ]
          plotActions[row.names(plotDataTh), step] = "none"
          #If demand has not been fulfilled continue
          if((Extracted<TenYrDemand)) {
            #Find difference between HB and thinningHB
            diffHB = plotDataTh$HB-prescriptionDataTh$thinningHB
            o = order(diffHB, decreasing=FALSE) #Sort from largest difference to more negative difference
            i = 1
            while((Extracted<TenYrDemand) && (i<=length(o))) {
              plotDataOne = plotDataTh[o[i], ]
              prescriptionDataOne = prescriptionDataTh[o[i],]
              prescriptionDataOne$plantingSpecies = plotDataOne$DominantSpecies
              treeDataOne = treeDataTh[treeDataTh$ID %in% row.names(plotDataOne), ]
              #update treeDataTh
              treeDataTh = treeDataTh[!(treeDataTh$ID %in% row.names(plotDataOne)), ]

              if(nrow(treeDataOne)>0) {
                treeDataOnePost = IFNmanagement(treeDataOne, prescriptionDataOne, verbose = verbose)

                #save next prescription
                prescriptionDataManagedSel$finalPreviousStage[!finalCutPlots][o[i]] = treeDataOnePost$nextPrescription$finalPreviousStage
                plotDataManagedSel$finalPreviousStage[!finalCutPlots][o[i]] = prescriptionDataManagedSel$finalPreviousStage[!finalCutPlots][o[i]]

                #store actions
                action = treeDataOnePost$action
                plotActions[names(action), step] = action

                #Calculate extracted volume
                if(nrow(treeDataOnePost$treeDataCut)>0) {

                  FCvolume = IFNvolume(treeDataOnePost$treeDataCut)
                  VCCarea = sum(FCvolume$VCC, na.rm=T)*as.numeric(plotDataOne$Area)
                  if(is.na(VCCarea)) {
                    stop("VCCarea = NA")
                  }
                  #Update extracted volume and density
                  volspp = tapply(FCvolume$VCC, FCvolume$Species, sum, na.rm=T)*as.numeric(plotDataOne$Area)
                  extractedVolumeSpp[names(volspp), step] = extractedVolumeSpp[names(volspp), step] + as.numeric(volspp)
                  volplot = tapply(FCvolume$VCC, FCvolume$ID, sum, na.rm=T)*as.numeric(plotDataOne$Area)
                  extractedVolumePlot[names(volplot), step] = extractedVolumePlot[names(volplot), step] + as.numeric(volplot)
                  voldomspp = tapply(volplot, plotData[names(volplot), "DominantSpecies"], sum, na.rm=T)
                  extractedVolumeDomSpp[names(voldomspp), step] = extractedVolumeDomSpp[names(voldomspp), step] + as.numeric(voldomspp)

                  densspp= tapply(treeDataOnePost$treeDataCut$N, FCvolume$Species, sum, na.rm=T)*as.numeric(plotDataOne$Area)
                  extractedDensitySpp[names(densspp), step] = extractedDensitySpp[names(densspp), step] + as.numeric(densspp)
                  densplot= tapply(treeDataOnePost$treeDataCut$N, FCvolume$ID, sum, na.rm=T)*as.numeric(plotDataOne$Area)
                  extractedDensityPlot[names(densplot), step] = extractedDensityPlot[names(densplot), step] + as.numeric(densplot)
                  densdomspp = tapply(densplot, plotData[names(densplot), "DominantSpecies"], sum, na.rm=T)
                  extractedDensityDomSpp[names(densdomspp), step] = extractedDensityDomSpp[names(densdomspp), step] + as.numeric(densdomspp)

                  #Update extracted biomass
                  FCbiomass = IFNbiomass(treeDataOnePost$treeDataCut, as.CO2 = as.CO2)
                  biomassspp = 0.001*tapply(FCbiomass$Roots, FCbiomass$Species, sum, na.rm=T)*as.numeric(plotDataOne$Area)
                  extractedRootBiomassSpp[names(biomassspp), step] = extractedRootBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
                  biomassplot = 0.001*tapply(FCbiomass$Roots, FCbiomass$ID, sum, na.rm=T)*as.numeric(plotDataOne$Area)
                  extractedRootBiomassPlot[names(biomassplot), step] = extractedRootBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
                  biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
                  extractedRootBiomassDomSpp[names(biomassdomspp), step] = extractedRootBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)

                  biomassspp = 0.001*tapply(FCbiomass$Stem, FCbiomass$Species, sum, na.rm=T)*as.numeric(plotDataOne$Area)
                  extractedStemBiomassSpp[names(biomassspp), step] = extractedStemBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
                  biomassplot = 0.001*tapply(FCbiomass$Stem, FCbiomass$ID, sum, na.rm=T)*as.numeric(plotDataOne$Area)
                  extractedStemBiomassPlot[names(biomassplot), step] = extractedStemBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
                  biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
                  extractedStemBiomassDomSpp[names(biomassdomspp), step] = extractedStemBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)

                  biomassspp = 0.001*tapply(FCbiomass$Branches, FCbiomass$Species, sum, na.rm=T)*as.numeric(plotDataOne$Area)
                  extractedBranchBiomassSpp[names(biomassspp), step] = extractedBranchBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
                  biomassplot = 0.001*tapply(FCbiomass$Branches, FCbiomass$ID, sum, na.rm=T)*as.numeric(plotDataOne$Area)
                  extractedBranchBiomassPlot[names(biomassplot), step] = extractedBranchBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
                  biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
                  extractedBranchBiomassDomSpp[names(biomassdomspp), step] = extractedBranchBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)

                  #Add to extracted
                  Extracted = Extracted + VCCarea
                  treated = treated+1
                  #Add to managed area
                  managedAreaSel = managedAreaSel + as.numeric(plotDataOne$Area)
                }
                #Store result
                if(!is.null(treeDataOut)) {
                  treeDataOut = rbind(treeDataOut, treeDataOnePost$treeDataOut)
                } else {
                  treeDataOut = treeDataOnePost$treeDataOut
                }
                if(!is.null(treeDataCut)) {
                  treeDataCut = rbind(treeDataCut, treeDataOnePost$treeDataCut)
                } else {
                  treeDataCut = treeDataOnePost$treeDataCut
                }
              }

              # print(treeDataOnePost$nextPrescription)
              #Increase plot counter
              i = i+1
            }
          }

          #Add the remaining (non-intervention plots) to treeDataOut
          if(nrow(treeDataTh)>0) {
            if(!is.null(treeDataOut)) {
              # print(head(treeDataTh))
              # print(head(treeDataOut))
              treeDataOut = rbind(treeDataOut, treeDataTh)
            } else {
              treeDataOut = treeDataTh
            }
          }
        }

        #Copy final previous stage to plotDataManaged
        plotDataManaged$finalPreviousStage[sel] = plotDataManagedSel$finalPreviousStage # copy to plotDataManaged (for next step)

        #Add to cut table (for output)
        if(cut.table) {
          # print("merging cut")
          if(!is.null(treeDataCut) && (nrow(treeDataCut)>0)) {
            treeDataCutStep = treeDataCut
            treeDataCutStep$Step = step
            if(is.null(treeDataCutSeq)) {
              treeDataCutSeq = treeDataCutStep[,c("Step","ID","Species", "DBH","N", "H")]
            } else {
              # print(names(treeDataCut))
              # print(names(treeDataCutSeq))
              treeDataCutSeq = rbind(treeDataCutSeq,treeDataCutStep[, names(treeDataCutSeq), drop=FALSE])
            }
          }
        }
        # Uodate demand offset (for next step)
        TenYrDemandOffset[r] = TenYrDemand - Extracted
        # if(!verbose) {
        #   cat(paste0(", Extr. ", round(Extracted),"m3, Diff. ", round(TenYrDemandOffset[r]),"m3, Treated ", treated,"]"))
        # } else {
        #   cat(paste0("\n  Management [Dem. ", round(TenYrDemand),"m3"))
        #   cat(paste0(", Extr. ", round(Extracted),"m3, Diff. ", round(TenYrDemandOffset[r]),"m3, Treated ", treated,"]"))
        # }
        extractedDemand[r,step] = Extracted
        extractedDemand[nrow(demand)+1,step] = extractedDemand[nrow(demand)+1,step]+Extracted
        managedArea[r,step] = managedAreaSel
        managedArea[nrow(demand)+1,step] = managedArea[nrow(demand)+1,step] + managedAreaSel

        # return(list(treeDataOut, plotDataManagedSel))
        #C.6 Run forest dynamics for all managed plots
        cat(paste(", dynamics"))
        treeDataManagedSelNew = NULL
        treeDataManagedSelDead = NULL
        idsIrregular = row.names(prescriptionDataManagedSel)[prescriptionDataManagedSel$type=="irregular"]
        tdoirr= treeDataOut[treeDataOut$ID %in% idsIrregular,, drop=FALSE]
        if(nrow(tdoirr)>0) {
          tdoirr$OIFFIN = as.character(1:nrow(tdoirr))
          ##Add cut records with OIFFIN = 0 to calculate BAL_ext
          if(!is.null(treeDataCut)) {
            tdci = treeDataCut[treeDataCut$ID %in% idsIrregular,, drop=FALSE]
            if(nrow(tdci)>0) {
              tdci$OIFFIN = "0"
              tdoirr = rbind(tdoirr,tdci)
            }
          }
          dyn = IFNdynamics(tdoirr, plotDataManagedSel[idsIrregular, , drop=FALSE],
                            nb = nb, NspNmatrix = NspNmatrix,
                            numSteps = 1,
                            ingrowth = T,
                            dead.table = T,
                            ...)
          treeDataManagedSelNew = dyn$out[,c("ID", "Species", "Name", "DBH", "N", "H"), drop = FALSE]
          treeDataManagedSelDead = dyn$dead[,c("ID", "Species", "Name", "DBH", "N", "H"), drop = FALSE]
        }
        idsRegular = row.names(prescriptionDataManagedSel)[prescriptionDataManagedSel$type!="irregular"]
        tdor= treeDataOut[treeDataOut$ID %in% idsRegular,, drop=FALSE]
        if(nrow(tdor)>0) {
          tdor$OIFFIN = as.character(1:nrow(tdor))
          ##Add cut records with OIFFIN = 0 to calculate BAL_ext
          if(!is.null(treeDataCut)) {
            tdci = treeDataCut[treeDataCut$ID %in% idsRegular,, drop=FALSE]
            if(nrow(tdci)>0) {
              tdci$OIFFIN = "0"
              tdor = rbind(tdor,tdci)
            }
          }
          dyn = IFNdynamics(tdor, plotDataManagedSel[idsRegular, , drop=FALSE],
                            nb = nb, NspNmatrix = NspNmatrix,
                            numSteps = 1,
                            ingrowth = F,
                            dead.table = T,
                            ...)
          if(!is.null(treeDataManagedSelNew)) treeDataManagedSelNew = rbind(treeDataManagedSelNew, dyn$out[,c("ID", "Species", "Name", "DBH", "N", "H"), drop = FALSE])
          else treeDataManagedSelNew = dyn$out[,c("ID", "Species", "Name", "DBH", "N", "H"), drop = FALSE]
          if(!is.null(treeDataManagedSelDead)) treeDataManagedSelDead = rbind(treeDataManagedSelDead, dyn$dead[,c("ID", "Species", "Name", "DBH", "N", "H"), drop = FALSE])
          else treeDataManagedSelDead = dyn$dead[,c("ID", "Species", "Name", "DBH", "N", "H"), drop = FALSE]
        }
        if(nrow(treeDataManagedSelNew)>0) {
          treeDataManagedSelNew$Provincia = plotDataManagedSel[treeDataManagedSelNew$ID, "Province"]
          treeDataManagedSelNew$Area = plotDataManagedSel[treeDataManagedSelNew$ID, "Area"]
        }
        if(nrow(treeDataManagedSelDead)>0) {
          treeDataManagedSelDead$Provincia = plotDataManagedSel[treeDataManagedSelDead$ID, "Province"]
          treeDataManagedSelDead$Area = plotDataManagedSel[treeDataManagedSelDead$ID, "Area"]
        }

        #C.7 Calculate standing/dead volume and density
        cat(paste(", volume"))
        mvolume = IFNvolume(treeDataManagedSelNew)
        volspp = tapply(mvolume$VCC*treeDataManagedSelNew$Area, mvolume$Species, sum, na.rm=T)
        standingVolumeSpp[names(volspp), step+1] = standingVolumeSpp[names(volspp), step+1] + as.numeric(volspp)
        volplot = tapply(mvolume$VCC*treeDataManagedSelNew$Area, mvolume$ID, sum, na.rm=T)
        standingVolumePlot[names(volplot), step+1] = standingVolumePlot[names(volplot), step+1] + as.numeric(volplot)
        voldomspp = tapply(volplot, plotData[names(volplot), "DominantSpecies"], sum, na.rm=T)
        standingVolumeDomSpp[names(voldomspp),step+ 1] = standingVolumeDomSpp[names(voldomspp), step+1] + as.numeric(voldomspp)

        mdeadvolume = IFNvolume(treeDataManagedSelDead)
        volspp = tapply(mdeadvolume$VCC*treeDataManagedSelDead$Area, mdeadvolume$Species, sum, na.rm=T)
        deadVolumeSpp[names(volspp), step] = deadVolumeSpp[names(volspp), step] + as.numeric(volspp)
        volplot = tapply(mdeadvolume$VCC*treeDataManagedSelDead$Area, mdeadvolume$ID, sum, na.rm=T)
        deadVolumePlot[names(volplot), step] = deadVolumePlot[names(volplot), step] + as.numeric(volplot)
        voldomspp = tapply(volplot, plotData[names(volplot), "DominantSpecies"], sum, na.rm=T)
        deadVolumeDomSpp[names(voldomspp), step] = deadVolumeDomSpp[names(voldomspp), step] + as.numeric(voldomspp)

        cat(paste(", density"))
        densspp = tapply(treeDataManagedSelNew$N*treeDataManagedSelNew$Area, mvolume$Species, sum, na.rm=T)
        standingDensitySpp[names(densspp), step+1] =  standingDensitySpp[names(densspp), step+1] + as.numeric(densspp)
        densplot =tapply(treeDataManagedSelNew$N*treeDataManagedSelNew$Area, mvolume$ID, sum, na.rm=T)
        standingDensityPlot[names(densplot), step+1] = standingVolumePlot[names(densplot), step+1] + as.numeric(densplot)
        densdomspp = tapply(densplot, plotData[names(densplot), "DominantSpecies"], sum, na.rm=T)
        standingDensityDomSpp[names(densdomspp),step+ 1] = standingDensityDomSpp[names(densdomspp), step+1] + as.numeric(densdomspp)

        densspp = tapply(treeDataManagedSelDead$N*treeDataManagedSelDead$Area, mdeadvolume$Species, sum, na.rm=T)
        deadDensitySpp[names(densspp), step] = as.numeric(densspp)
        densplot = tapply(treeDataManagedSelDead$N*treeDataManagedSelDead$Area, mdeadvolume$ID, sum, na.rm=T)
        deadDensityPlot[names(densplot), step] = as.numeric(densplot)
        densdomspp = tapply(densplot, plotData[names(densplot), "DominantSpecies"], sum, na.rm=T)
        deadDensityDomSpp[names(densdomspp), step] = deadDensityDomSpp[names(densdomspp), step] + as.numeric(densdomspp)

        #C.8 Calculate standing/dead biomass
        cat(paste(", biomass"))
        mlivebiomass = IFNbiomass(treeDataManagedSelNew, as.CO2 = as.CO2)
        biomassspp = 0.001*tapply(mlivebiomass$Roots*treeDataManagedSelNew$Area, mlivebiomass$Species, sum, na.rm=T)
        standingRootBiomassSpp[names(biomassspp), step+1] = standingRootBiomassSpp[names(biomassspp), step+1] + as.numeric(biomassspp)
        biomassplot = 0.001*tapply(mlivebiomass$Roots*treeDataManagedSelNew$Area, mlivebiomass$ID, sum, na.rm=T)
        standingRootBiomassPlot[names(biomassplot), step+1] = standingRootBiomassPlot[names(biomassplot), step+1] + as.numeric(biomassplot)
        biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
        standingRootBiomassDomSpp[names(biomassdomspp), step+1] = standingRootBiomassDomSpp[names(biomassdomspp), step+1] + as.numeric(biomassdomspp)

        biomassspp = 0.001*tapply(mlivebiomass$Stem*treeDataManagedSelNew$Area, mlivebiomass$Species, sum, na.rm=T)
        standingStemBiomassSpp[names(biomassspp), step+1] = standingStemBiomassSpp[names(biomassspp), step+1] + as.numeric(biomassspp)
        biomassplot = 0.001*tapply(mlivebiomass$Stem*treeDataManagedSelNew$Area, mlivebiomass$ID, sum, na.rm=T)
        standingStemBiomassPlot[names(biomassplot), step+1] = standingStemBiomassPlot[names(biomassplot), step+1] + as.numeric(biomassplot)
        biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
        standingStemBiomassDomSpp[names(biomassdomspp), step+1] = standingStemBiomassDomSpp[names(biomassdomspp), step+1] + as.numeric(biomassdomspp)

        biomassspp = 0.001*tapply(mlivebiomass$Branches*treeDataManagedSelNew$Area, mlivebiomass$Species, sum, na.rm=T)
        standingBranchBiomassSpp[names(biomassspp), step+1] = standingBranchBiomassSpp[names(biomassspp), step+1] + as.numeric(biomassspp)
        biomassplot = 0.001*tapply(mlivebiomass$Branches*treeDataManagedSelNew$Area, mlivebiomass$ID, sum, na.rm=T)
        standingBranchBiomassPlot[names(biomassplot), step+1] = standingBranchBiomassPlot[names(biomassplot), step+1] + as.numeric(biomassplot)
        biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
        standingBranchBiomassDomSpp[names(biomassdomspp), step+1] = standingBranchBiomassDomSpp[names(biomassdomspp), step+1] + as.numeric(biomassdomspp)

        mdeadbiomass = IFNbiomass(treeDataManagedSelDead, as.CO2 = as.CO2)
        biomassspp = 0.001*tapply(mdeadbiomass$Roots*treeDataManagedSelDead$Area, mdeadbiomass$Species, sum, na.rm=T)
        deadRootBiomassSpp[names(biomassspp), step] = deadRootBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
        biomassplot = 0.001*tapply(mdeadbiomass$Roots*treeDataManagedSelDead$Area, mdeadbiomass$ID, sum, na.rm=T)
        deadRootBiomassPlot[names(biomassplot), step] = deadRootBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
        biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
        deadRootBiomassDomSpp[names(biomassdomspp), step] = deadRootBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)

        biomassspp = 0.001*tapply(mdeadbiomass$Stem*treeDataManagedSelDead$Area, mdeadbiomass$Species, sum, na.rm=T)
        deadStemBiomassSpp[names(biomassspp), step] = deadStemBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
        biomassplot = 0.001*tapply(mdeadbiomass$Stem*treeDataManagedSelDead$Area, mdeadbiomass$ID, sum, na.rm=T)
        deadStemBiomassPlot[names(biomassplot), step] = deadStemBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
        biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
        deadStemBiomassDomSpp[names(biomassdomspp), step] = deadStemBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)

        biomassspp = 0.001*tapply(mdeadbiomass$Branches*treeDataManagedSelDead$Area, mdeadbiomass$Species, sum, na.rm=T)
        deadBranchBiomassSpp[names(biomassspp), step] = deadBranchBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
        biomassplot = 0.001*tapply(mdeadbiomass$Branches*treeDataManagedSelDead$Area, mdeadbiomass$ID, sum, na.rm=T)
        deadBranchBiomassPlot[names(biomassplot), step] = deadBranchBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
        biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
        deadBranchBiomassDomSpp[names(biomassdomspp), step] = deadBranchBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)

        cat(paste(".\n"))
        if(nrow(treeDataManagedSelDead)>0) {
          if(is.null(treeDataManagedDead)) {
            treeDataManagedDead = treeDataManagedSelDead
          } else {
            treeDataManagedDead = rbind(treeDataManagedDead,treeDataManagedSelDead)
          }
        }
        if(nrow(treeDataManagedSel)>0) {
          if(is.null(treeDataManagedNew)){
            treeDataManagedNew = treeDataManagedSelNew
          } else {
            treeDataManagedNew = rbind(treeDataManagedNew, treeDataManagedSelNew)
          }
        }

      } else {
        cat(paste0("  [",r,"/",nrow(demand),"] - ", sum(sel), " plots dominated by ", demand$Name[r], "\n"))
      }
    }

    #Copy final previous stage to plotData
    plotData[managedIDs, "finalPreviousStage"] = plotDataManaged$finalPreviousStage

    cat("\n") ## End of managed plots (C)

    # print(head(treeDataManagedNew))
    # print(head(treeDataManagedDead))

    ## Assemble tree data for next step

    treeData = rbind(treeDataNonManagedNew, treeDataManagedNew)
    treeDataDead = rbind(treeDataNonManagedDead, treeDataManagedDead)

    if(sequence) {
      treeDataStep = treeData
      treeDataStep$Step = step
      treeDataSeq = rbind(treeDataSeq, treeDataStep[,names(treeDataSeq)])
    }
    if(dead.table) {
      treeDataDead$Step = step
      if(is.null(treeDataDeadSeq)) {
        treeDataDeadSeq = treeDataDead[,c("Step","ID","Species", "DBH","N", "H"), drop=FALSE]
      } else {
        treeDataDeadSeq = rbind(treeDataDeadSeq,treeDataDead[,names(treeDataDeadSeq), drop=FALSE])
      }
    }

  } # END OF STEP LOOP




  # C. Assemble output
  cat(paste("C. Post-processing...\n"))
  cat("  Calculating coverage of demand.\n")
  extractedDemand = as.data.frame(extractedDemand)
  extractedDemand$AnnualMean = rowMeans(extractedDemand)/10
  extractedDemand$Demand = c(demand$AnnualDemand, sum(demand$AnnualDemand))
  extractedDemand$Coverage = 100*extractedDemand$AnnualMean/extractedDemand$Demand

  cat("  Calculating managed area.\n")
  managedArea = as.data.frame(managedArea)
  md = tapply(plotData$Area[plotData$Managed], plotData$DemandRow[plotData$Managed], sum, na.rm =T)
  managedArea$ManageableArea = 0
  managedArea$ManageableArea[as.numeric(names(md))] = md
  managedArea$ManageableArea[nrow(managedArea)] = sum(md)

  totalArea = sum(plotData$Area)


  #Calculate variation
  # variation_t = standing_t - standing_(t-1)
  cat("  Calculating variation of stands between steps.\n")
  variationVolumeSpp = standingVolumeSpp[,2:(numSteps+1), drop=FALSE] - standingVolumeSpp[,1:(numSteps), drop=FALSE]
  variationVolumePlot = standingVolumePlot[,2:(numSteps+1), drop=FALSE] - standingVolumePlot[,1:(numSteps), drop=FALSE]
  variationVolumeDomSpp = standingVolumeDomSpp[,2:(numSteps+1), drop=FALSE] - standingVolumeDomSpp[,1:(numSteps), drop=FALSE]

  variationDensitySpp = standingDensitySpp[,2:(numSteps+1), drop=FALSE] - standingDensitySpp[,1:(numSteps), drop=FALSE]
  variationDensityPlot = standingDensityPlot[,2:(numSteps+1), drop=FALSE] - standingDensityPlot[,1:(numSteps), drop=FALSE]
  variationDensityDomSpp = standingDensityDomSpp[,2:(numSteps+1), drop=FALSE] - standingDensityDomSpp[,1:(numSteps), drop=FALSE]

  variationRootBiomassSpp = standingRootBiomassSpp[,2:(numSteps+1), drop=FALSE] - standingRootBiomassSpp[,1:(numSteps), drop=FALSE]
  variationRootBiomassPlot = standingRootBiomassPlot[,2:(numSteps+1), drop=FALSE] - standingRootBiomassPlot[,1:(numSteps), drop=FALSE]
  variationRootBiomassDomSpp = standingRootBiomassDomSpp[,2:(numSteps+1), drop=FALSE] - standingRootBiomassDomSpp[,1:(numSteps), drop=FALSE]

  variationStemBiomassSpp = standingStemBiomassSpp[,2:(numSteps+1), drop=FALSE] - standingStemBiomassSpp[,1:(numSteps), drop=FALSE]
  variationStemBiomassPlot = standingStemBiomassPlot[,2:(numSteps+1), drop=FALSE] - standingStemBiomassPlot[,1:(numSteps), drop=FALSE]
  variationStemBiomassDomSpp = standingStemBiomassDomSpp[,2:(numSteps+1), drop=FALSE] - standingStemBiomassDomSpp[,1:(numSteps), drop=FALSE]

  variationBranchBiomassSpp = standingBranchBiomassSpp[,2:(numSteps+1), drop=FALSE] - standingBranchBiomassSpp[,1:(numSteps), drop=FALSE]
  variationBranchBiomassPlot = standingBranchBiomassPlot[,2:(numSteps+1), drop=FALSE] - standingBranchBiomassPlot[,1:(numSteps), drop=FALSE]
  variationBranchBiomassDomSpp = standingBranchBiomassDomSpp[,2:(numSteps+1), drop=FALSE] - standingBranchBiomassDomSpp[,1:(numSteps), drop=FALSE]

  #Calculate stand growth (tree growth + mortality + recruitment)
  cat("  Calculating stand growth between steps.\n")
  # Stand growth_t = standing_t - postmanagement_t = standing_t - (standing_(t-1) - extracted_t) = variation_t + extracted_t
  standgrowthVolumeSpp = variationVolumeSpp+extractedVolumeSpp
  standgrowthVolumePlot = variationVolumePlot+extractedVolumePlot
  standgrowthVolumeDomSpp = variationVolumeDomSpp+extractedVolumeDomSpp

  standgrowthDensitySpp = variationDensitySpp+extractedDensitySpp
  standgrowthDensityPlot = variationDensityPlot+extractedDensityPlot
  standgrowthDensityDomSpp = variationDensityDomSpp+extractedDensityDomSpp

  standgrowthRootBiomassSpp = variationRootBiomassSpp+extractedRootBiomassSpp
  standgrowthRootBiomassPlot = variationRootBiomassPlot+extractedRootBiomassPlot
  standgrowthRootBiomassDomSpp = variationRootBiomassDomSpp+extractedRootBiomassDomSpp

  standgrowthStemBiomassSpp = variationStemBiomassSpp+extractedStemBiomassSpp
  standgrowthStemBiomassPlot = variationStemBiomassPlot+extractedStemBiomassPlot
  standgrowthStemBiomassDomSpp = variationStemBiomassDomSpp+extractedStemBiomassDomSpp

  standgrowthBranchBiomassSpp = variationBranchBiomassSpp+extractedBranchBiomassSpp
  standgrowthBranchBiomassPlot = variationBranchBiomassPlot+extractedBranchBiomassPlot
  standgrowthBranchBiomassDomSpp = variationBranchBiomassDomSpp+extractedBranchBiomassDomSpp

  #Stats per hectare
  cat("  Calculating statistics per hectare.\n")
  volumeStats = data.frame(standing = colSums(standingVolumeSpp)/totalArea,
                           variation = c(NA,colSums(variationVolumeSpp)/totalArea),
                           standgrowth = c(NA,colSums(standgrowthVolumeSpp)/totalArea),
                           dead = c(NA,colSums(deadVolumeSpp)/totalArea),
                           extracted = c(NA,colSums(extractedVolumeSpp)/totalArea))
  densityStats = data.frame(standing = colSums(standingDensitySpp)/totalArea,
                           variation = c(NA,colSums(variationDensitySpp)/totalArea),
                           standgrowth = c(NA,colSums(standgrowthDensitySpp)/totalArea),
                           dead = c(NA,colSums(deadDensitySpp)/totalArea),
                           extracted = c(NA,colSums(extractedDensitySpp)/totalArea))
  rootBiomassStats = data.frame(standing = colSums(standingRootBiomassSpp)/totalArea,
                                variation = c(NA,colSums(variationRootBiomassSpp)/totalArea),
                                standgrowth = c(NA,colSums(standgrowthRootBiomassSpp)/totalArea),
                                dead = c(NA,colSums(deadRootBiomassSpp)/totalArea),
                                extracted = c(NA,colSums(extractedRootBiomassSpp)/totalArea))
  stemBiomassStats = data.frame(standing = colSums(standingStemBiomassSpp)/totalArea,
                                variation = c(NA,colSums(variationStemBiomassSpp)/totalArea),
                                standgrowth = c(NA,colSums(standgrowthStemBiomassSpp)/totalArea),
                                dead = c(NA,colSums(deadStemBiomassSpp)/totalArea),
                                extracted = c(NA,colSums(extractedStemBiomassSpp)/totalArea))
  branchBiomassStats = data.frame(standing = colSums(standingBranchBiomassSpp)/totalArea,
                                  variation = c(NA,colSums(variationBranchBiomassSpp)/totalArea),
                                  standgrowth = c(NA,colSums(standgrowthBranchBiomassSpp)/totalArea),
                                  dead = c(NA,colSums(deadBranchBiomassSpp)/totalArea),
                                  extracted = c(NA,colSums(extractedBranchBiomassSpp)/totalArea))

  cat("  Assembling output.\n")
  volume = list(summary = volumeStats,
                standingSpp = standingVolumeSpp,
                standingPlot = standingVolumePlot,
                standingDomSpp = standingVolumeDomSpp,
                variationSpp = variationVolumeSpp,
                variationPlot = variationVolumePlot,
                variationDomSpp = variationVolumeDomSpp,
                standgrowthSpp = standgrowthVolumeSpp,
                standgrowthPlot = standgrowthVolumePlot,
                standgrowthDomSpp = standgrowthVolumeDomSpp,
                deadSpp = deadVolumeSpp,
                deadPlot = deadVolumePlot,
                deadDomSpp = deadVolumeDomSpp,
                extractedSpp = extractedVolumeSpp,
                extractedPlot = extractedVolumePlot,
                extractedDomSpp = extractedVolumeDomSpp)
  density = list(summary = densityStats,
                 standingSpp = standingDensitySpp,
                 standingPlot = standingDensityPlot,
                 standingDomSpp = standingDensityDomSpp,
                 variationSpp = variationDensitySpp,
                 variationPlot = variationDensityPlot,
                 variationDomSpp = variationDensityDomSpp,
                 standgrowthSpp = standgrowthDensitySpp,
                 standgrowthPlot = standgrowthDensityPlot,
                 standgrowthDomSpp = standgrowthDensityDomSpp,
                 deadSpp = deadDensitySpp,
                 deadPlot = deadDensityPlot,
                 deadDomSpp = deadDensityDomSpp,
                 extractedSpp = extractedDensitySpp,
                 extractedPlot = extractedDensityPlot,
                 extractedDomSpp = extractedDensityDomSpp)
  rootBiomass = list(summary = rootBiomassStats,
                     standingSpp = standingRootBiomassSpp,
                     standingPlot = standingRootBiomassPlot,
                     standingDomSpp = standingRootBiomassDomSpp,
                     variationSpp = variationRootBiomassSpp,
                     variationPlot = variationRootBiomassPlot,
                     variationDomSpp = variationRootBiomassDomSpp,
                     standgrowthSpp = standgrowthRootBiomassSpp,
                     standgrowthPlot = standgrowthRootBiomassPlot,
                     standgrowthDomSpp = standgrowthRootBiomassDomSpp,
                     deadSpp = deadRootBiomassSpp,
                     deadPlot = deadRootBiomassPlot,
                     deadDomSpp = deadRootBiomassDomSpp,
                     extractedSpp = extractedRootBiomassSpp,
                     extractedPlot = extractedRootBiomassPlot,
                     extractedDomSpp = extractedRootBiomassDomSpp)
  stemBiomass = list(summary = stemBiomassStats,
                     standingSpp = standingStemBiomassSpp,
                     standingPlot = standingStemBiomassPlot,
                     standingDomSpp = standingStemBiomassDomSpp,
                     variationSpp = variationStemBiomassSpp,
                     variationPlot = variationStemBiomassPlot,
                     variationDomSpp = variationStemBiomassDomSpp,
                     standgrowthSpp = standgrowthStemBiomassSpp,
                     standgrowthPlot = standgrowthStemBiomassPlot,
                     standgrowthDomSpp = standgrowthStemBiomassDomSpp,
                     deadSpp = deadStemBiomassSpp,
                     deadPlot = deadStemBiomassPlot,
                     deadDomSpp = deadStemBiomassDomSpp,
                     extractedSpp = extractedStemBiomassSpp,
                     extractedPlot = extractedStemBiomassPlot,
                     extractedDomSpp = extractedStemBiomassDomSpp)
  branchBiomass = list(summary = branchBiomassStats,
                       standingSpp = standingBranchBiomassSpp,
                       standingPlot = standingBranchBiomassPlot,
                       standingDomSpp = standingBranchBiomassDomSpp,
                       variationSpp = variationBranchBiomassSpp,
                       variationPlot = variationBranchBiomassPlot,
                       variationDomSpp = variationBranchBiomassDomSpp,
                       standgrowthSpp = standgrowthBranchBiomassSpp,
                       standgrowthPlot = standgrowthBranchBiomassPlot,
                       standgrowthDomSpp = standgrowthBranchBiomassDomSpp,
                       deadSpp = deadBranchBiomassSpp,
                       deadPlot = deadBranchBiomassPlot,
                       deadDomSpp = deadBranchBiomassDomSpp,
                       extractedSpp = extractedBranchBiomassSpp,
                       extractedPlot = extractedBranchBiomassPlot,
                       extractedDomSpp = extractedBranchBiomassDomSpp)

  plotDataSuppl = plotData[,c("Area", "DominantSpecies", "DemandRow", "PrescriptionRow", "finalPreviousStage")]
  res = list(plots = nrow(plotData),
             steps = numSteps,
             stepOffset = stepOffset,
             totalArea = totalArea,
             plotDataSuppl = plotDataSuppl,
             plotActions = plotActions,
             extractedDemand=extractedDemand,
             demandOffset = TenYrDemandOffset,
             managedArea = managedArea,
             volume=volume,
             density=density,
             rootBiomass = rootBiomass,
             stemBiomass = stemBiomass,
             branchBiomass = branchBiomass,
             treeDataNew = treeData)
  if(sequence) {
    if((!is.null(treeDataSeq)) && (nrow(treeDataSeq)>0)) row.names(treeDataSeq) = 1:nrow(treeDataSeq)
    res[["treeDataSequence"]] = treeDataSeq
  }
  if(dead.table) {
    if((!is.null(treeDataDeadSeq)) && (nrow(treeDataDeadSeq)>0)) row.names(treeDataDeadSeq) = 1:nrow(treeDataDeadSeq)
    res[["treeDataDead"]] = treeDataDeadSeq
  }
  if(cut.table) {
    if((!is.null(treeDataCutSeq)) && (nrow(treeDataCutSeq)>0)) row.names(treeDataCutSeq) = 1:nrow(treeDataCutSeq)
    res[["treeDataCut"]] = treeDataCutSeq
  }
  class(res) <- "IFNscenario"
  return(res)
}