R/plot.spwb.R

In medfate: Mediterranean Forest Simulation

#' Plots simulation results
#'
#' Function \code{plot} plots time series of the results of the soil plant water balance model (see \code{\link{spwb}}), 
#' plant water balance model (see \code{\link{pwb}}), the forest growth model (see \code{\link{growth}}) 
#' or the forest dynamics model (see \code{\link{fordyn}}).
#' 
#' @param x An object of class \code{spwb}, \code{pwb}, \code{growth} or \code{fordyn}.
#' @param type The information to be plotted (see details)
#' @param cohorts An integer, boolean or character vector to select the plant cohorts to be plotted. If \code{cohorts = "T"} (resp. \code{cohorts = "S"}) then all tree (resp. shrub) cohorts will be displayed.
#' @param bySpecies Allows aggregating output by species, before drawing plots (only has an effect with some values of \code{type}). Aggregation can involve a sum (as for plant lai or transpiration) or a LAI-weighted mean (as for plant stress or plant water potential), where LAI values are those of \code{LAIlive}.
#' @param dates A Date vector with a subset of dates to be plotted.
#' @param subdaily Whether subdaily results should be shown, only for simulations using \code{transpirationMode = "Sperry"} and having set \code{subdailyResults = TRUE} in the simulation control object. If \code{subdaily = TRUE}, then the valid strings for \code{type} are listed in \code{\link{plot.spwb_day}}.
#' @param xlim Range of values for x.
#' @param ylim Range of values for y.
#' @param xlab x-axis label.
#' @param ylab y-axis label.
#' @param summary.freq Frequency of summary statistics (see \code{\link{cut.Date}}).
#' @param ... Additional parameters for function \code{plot} (not used).
#' 
#' 
#' @details 
#' The following plots are currently available for \code{\link{spwb}} (most of them also for \code{\link{pwb}}):
#' \itemize{
#'   \item{\code{"PET_Precipitation"}:}{ Potential evapotranspiration and Precipitation.}
#'   \item{\code{"PET_NetRain"}:}{ Potential evapotranspiration and Net rainfall.}
#'   \item{\code{"Snow"}:}{ Snow precipitation and snowpack dynamics.}
#'   \item{\code{"Export"}:}{ Water exported through deep drainage and surface runoff.}
#'   \item{\code{"Evapotranspiration"}:}{ Plant transpiration and soil evaporation.}
#'   \item{\code{"SoilPsi"}:}{ Soil water potential.}
#'   \item{\code{"SoilRWC"}:}{ Soil relative water content (in percent of field capacity).}
#'   \item{\code{"SoilTheta"}:}{ Soil moisture water content (in percent volume).}
#'   \item{\code{"SoilVol"}:}{ Soil water volumetric content (in mm).}
#'   \item{\code{"PlantExtraction"}:}{ Water extracted by plants from each soil layer.}
#'   \item{\code{"HydraulicRedistribution"}: Water added to each soil layer coming from other soil layers, transported through the plant hydraulic network.}
#'   \item{\code{"WTD"}:}{ Water table depth.}
#'   \item{\code{"LAI"}:}{ Expanded and dead leaf area index of the whole stand.}
#'   \item{\code{"PlantLAI"}:}{ Plant cohort leaf area index (expanded leaves).}
#'   \item{\code{"PlantLAIlive"}:}{ Plant cohort leaf area index ("live" leaves).}
#'   \item{\code{"PlantStress"}:}{ Plant cohort average daily drought stress.}
#'   \item{\code{"PlantTranspiration"}:}{ Plant cohort transpiration.}
#'   \item{\code{"TranspirationPerLeaf"}:}{ Plant cohort transpiration per leaf area.}
#'   \item{\code{"PlantGrossPhotosynthesis"}:}{ Plant cohort photosynthesis.}
#'   \item{\code{"GrossPhotosynthesisPerLeaf"}:}{ Plant cohort photosynthesis per leaf area.}
#'   \item{\code{"StemRWC"}:  Average daily stem relative water content.}
#'   \item{\code{"LeafRWC"}:  Average daily leaf relative water content.}
#'   \item{\code{"LFMC"}:  Live fuel moisture content.}
#' }
#' The following plots are available for \code{\link{spwb}} and \code{\link{pwb}} \emph{only if} \code{transpirationMode = "Granier"}:
#'   \itemize{
#'     \item{\code{"PlantPsi"}:}{ Plant cohort water potential.}
#'     \item{\code{"FPAR"}:}{ Fraction of PAR at the canopy level of each plant cohort.}
#'     \item{\code{"AbsorbedSWRFraction"}:}{ Fraction of SWR absorbed by each plant cohort.}
#'   }
#' The following plots are available for \code{\link{spwb}} and \code{\link{pwb}} \emph{only if} \code{transpirationMode = "Sperry"}:
#'   \itemize{
#'     \item{\code{"SoilPlantConductance"}: Average instantaneous overall soil plant conductance (calculated as the derivative of the supply function).}
#'     \item{\code{"LeafPsiMin"}:  Midday leaf water potential.}
#'     \item{\code{"LeafPsiMax"}:  Pre-dawn leaf water potential.}
#'     \item{\code{"LeafPsiRange"}:  Range of leaf water potential.}
#'     \item{\code{"LeafPsiMin_SL"}: Minimum water potential of sunlit leaves.}
#'     \item{\code{"LeafPsiMax_SL"}: Maximum water potential of sunlit leaves.}
#'     \item{\code{"LeafPsiMin_SH"}: Minimum water potential of shade leaves.}
#'     \item{\code{"LeafPsiMax_SH"}: Maximum water potential of shade leaves.}
#'     \item{\code{"TempMin_SL"}: Minimum temperature of sunlit leaves.}
#'     \item{\code{"TempMax_SL"}: Maximum temperature of sunlit leaves.}
#'     \item{\code{"TempMin_SH"}: Minimum temperature of shade leaves.}
#'     \item{\code{"TempMax_SH"}: Maximum temperature of shade leaves.}
#'     \item{\code{"GSWMin_SL"}: Minimum stomatal conductance of sunlit leaves.}
#'     \item{\code{"GSWMax_SL"}: Maximum stomatal conductance of sunlit leaves.}
#'     \item{\code{"GSWMin_SH"}: Minimum stomatal conductance of shade leaves.}
#'     \item{\code{"GSWMax_SH"}: Maximum stomatal conductance of shade leaves.}
#'     \item{\code{"StemPsi"}:  Midday (upper) stem water potential.}
#'     \item{\code{"RootPsi"}:  Midday root crown water potential.}
#'     \item{\code{"PlantNetPhotosynthesis"}:  Plant cohort net photosynthesis.}
#'     \item{\code{"NetPhotosynthesisPerLeaf"}:  Plant cohort net photosynthesis per leaf area.}
#'     \item{\code{"PlantWUE"}:  Plant cohort daily water use efficiency.}
#'     \item{\code{"PlantAbsorbedSWR"}:  Plant cohort absorbed short wave radiation.}
#'     \item{\code{"AbsorbedSWRPerLeaf"}:  Plant cohort absorbed short wave radiation per leaf area.}
#'     \item{\code{"PlantNetLWR"}:  Plant cohort net long wave radiation.}
#'     \item{\code{"NetLWRPerLeaf"}:  Plant cohort net long wave radiation per leaf area.}
#'     \item{\code{"AirTemperature"}:  Minimum/maximum/mean daily temperatures above canopy.}
#'     \item{\code{"CanopyTemperature"}:  Minimum/maximum/mean daily temperatures inside canopy.}
#'     \item{\code{"SoilTemperature"}:  Minimum/maximum/mean daily temperatures inside the first soil layer.}
#'     \item{\code{"CanopyEnergyBalance"}:  Canopy energy balance components.}
#'     \item{\code{"SoilEnergyBalance"}:  Soil energy balance components.}
#'   }
#' In addition to the former, the following plots are available for objects \code{\link{growth}} or \code{\link{fordyn}}:
#'   \itemize{
#'     \item{\code{"CarbonBalance"}:  Stand-level carbon balance components.}
#'     \item{\code{"BiomassBalance"}:  Stand-level biomass balance components.}
#'     \item{\code{"GrossPhotosynthesis"}:  Gross photosynthesis rate per dry weight.}
#'     \item{\code{"MaintenanceRespiration"}:  Maintenance respiration cost per dry weight.}
#'     \item{\code{"PhotosynthesisMaintenanceRatio"}:  The ratio of gross photosynthesis over maintenance respiration.}
#'     \item{\code{"RootExudation"}:  Root exudation rate per dry weight.}
#'     \item{\code{"LabileCarbonBalance"}:  Labile carbon balance per dry weight.}
#'     \item{\code{"SugarLeaf"}:  Sugar concentration in leaves.}
#'     \item{\code{"StarchLeaf"}:  Starch concentration in leaves.}
#'     \item{\code{"SugarSapwood"}:  Sugar concentration in sapwood.}
#'     \item{\code{"StarchSapwood"}:  Starch concentration in sapwood.}
#'     \item{\code{"SugarTransport"}:  Phloem sugar transport rate.}
#'     \item{\code{"StructuralBiomassBalance"}:  Daily structural biomass balance (g dry · ind-2).}
#'     \item{\code{"LabileBiomassBalance"}:  Daily labile biomass balance (g dry · ind-2).}
#'     \item{\code{"PlantBiomassBalance"}:  Daily plant biomass balance, i.e. labile change + structural change (g dry · ind-2).}
#'     \item{\code{"MortalityBiomassLoss"}:  Biomass loss due to mortality (g dry · m-2).}    
#'     \item{\code{"PlantBiomassBalance"}:  Daily cohort biomass balance (including mortality) (g dry · m-2).}
#'     \item{\code{"LeafBiomass"}: Leaf structural dry biomass per individual.}
#'     \item{\code{"SapwoodBiomass"}: Sapwood dry biomass per individual.}
#'     \item{\code{"FineRootBiomass"}: Fine root dry biomass per individual.}
#'     \item{\code{"SapwoodArea"}:  Sapwood area per individual.}
#'     \item{\code{"LeafArea"}:  Leaf area per individual.}
#'     \item{\code{"FineRootArea"}: Fine root area per individual (only for \code{transpirationMode = "Sperry"} or \code{transpirationMode = "Cochard"}).}
#'     \item{\code{"DBH"}: Diameter at breast height (in cm) for an average individual of each plant cohort.}
#'     \item{\code{"Height"}: Height (in cm) for an average individual of each plant cohort.}
#'     \item{\code{"SAgrowth"}:  Sapwood area growth rate.}
#'     \item{\code{"LAgrowth"}:  Leaf area growth rate.}
#'     \item{\code{"FRAgrowth"}:  Fine root area growth rate (only for \code{transpirationMode = "Sperry"} or \code{transpirationMode = "Cochard"}).}
#'     \item{\code{"HuberValue"}:  Ratio of leaf area to sapwood area.}
#'     \item{\code{"RootAreaLeafArea"}:  Ratio of fine root area to leaf area (only for \code{transpirationMode = "Sperry"} or \code{transpirationMode = "Cochard"}).}
#'   }
#' Finally, the following plots are only available for \code{\link{fordyn}} simulation results:
#'   \itemize{
#'     \item{\code{"StandBasalArea"}: Stand basal area of living trees.}
#'     \item{\code{"StandDensity"}: Stand density of living trees.}
#'     \item{\code{"SpeciesBasalArea"}: Basal area of living trees by species.}
#'     \item{\code{"SpeciesDensity"}: Density of living trees by species.}
#'     \item{\code{"CohortBasalArea"}: Basal area of living trees by plant cohort.}
#'     \item{\code{"CohortDensity"}: Density of living trees by plant cohort.}
#'   }
#'
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' @return An ggplot object
#' 
#' @seealso  \code{\link{spwb}}, \code{\link{pwb}}, \code{\link{growth}}, \code{\link{fordyn}}, \code{\link{summary.spwb}}
#' 
#' @examples 
#' #Load example daily meteorological data
#' data(examplemeteo)
#' 
#' #Load example plot plant data
#' data(exampleforestMED)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Initialize soil with default soil params (2 layers)
#' examplesoil = soil(defaultSoilParams(2))
#' 
#' #Initialize control parameters
#' control = defaultControl("Granier")
#' 
#' #Initialize input
#' x = forest2spwbInput(exampleforestMED,examplesoil, SpParamsMED, control)
#' 
#' #Call simulation function
#' S1<-spwb(x, examplemeteo, latitude = 41.82592, elevation = 100)
#' 
#' #Plot results
#' plot(S1)
#' 
plot.spwb<-function(x, type="PET_Precipitation", cohorts = NULL, bySpecies = FALSE,
                    dates = NULL, subdaily = FALSE, 
                    xlim = NULL, ylim=NULL, xlab=NULL, ylab=NULL, 
                    summary.freq = NULL, ...) {
  
  if(subdaily) return(.plotsubdaily(x,type, cohorts, bySpecies, dates, 
                                    xlim, ylim, xlab, ylab))
  
  if(("spwb" %in% class(x)) || ("pwb" %in% class(x))) {
    input = x$spwbInput
  } else {
    input = x$growthInput
  }
  TYPES = .getDailySPWBPlotTypes(input$control$transpirationMode)  

  type = match.arg(type,TYPES)  
  if(is.null(xlab)) xlab = ""
  if(type %in% c("PET_Precipitation", "PET_NetRain", "Evapotranspiration", "Snow",
                 "WTD", "Export", "SoilVol")) {
    return(.plot_wb(WaterBalance = x$WaterBalance, Soil = x$Soil, input_soil = input$soil, 
                    type = type, dates = dates, 
                    xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                    summary.freq = summary.freq, ...))
  } 
  return(plot.pwb(x, type=type, cohorts = cohorts, bySpecies = bySpecies,
                  dates = dates, xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                  summary.freq = summary.freq, ...))
}

#' @rdname plot.spwb
plot.pwb<-function(x, type="PlantTranspiration", cohorts = NULL, bySpecies = FALSE,
                   dates = NULL, subdaily = FALSE,
                   xlim = NULL, ylim=NULL, xlab=NULL, ylab=NULL, 
                   summary.freq = NULL,...) {
  
  if(subdaily) return(.plotsubdaily(x,type, cohorts, bySpecies, dates, 
                                    xlim, ylim, xlab, ylab))
  
  if(("spwb" %in% class(x)) || ("pwb" %in% class(x))) {
    input = x$spwbInput
  } else {
    input = x$growthInput
  }
  nlayers = length(input$soil$W)
  transpMode = input$control$transpirationMode
  
  if(is.null(cohorts))  cohorts = row.names(input$cohorts)
  else if(length(cohorts)==1) {
    if(cohorts=="T") {
      cohorts = row.names(input$cohorts)
      cohorts = cohorts[substr(cohorts,1,1)=="T"]
    } else if(cohorts=="S") {
      cohorts = row.names(input$cohorts)
      cohorts = cohorts[substr(cohorts,1,1)=="S"]
    }
  }

  spnames = as.character(input$cohorts[cohorts,"Name"])
  
  PlantsLAIlive = x$Plants$LAIlive[,cohorts, drop=FALSE]
  PlantsLAI = x$Plants$LAI[,cohorts, drop=FALSE]
  
  TYPES = .getDailyPWBPlotTypes(transpMode)
  type = match.arg(type,TYPES)  
  if(is.null(xlab)) xlab = ""  
  
  if(type %in% c("SoilPsi", "SoilTheta", "SoilRWC", "PlantExtraction", "HydraulicRedistribution")) {
    return(.plot_soil(Soil = x$Soil, input_soil = input$soil, input_control = input$control,
                    type = type, dates = dates, 
                    xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                    summary.freq = summary.freq, ...))
  }
  else if(type %in% c("LAI", "GroundIrradiance")) {
    return(.plot_stand(Stand = x$Stand,
                      type = type, dates = dates, 
                      xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                      summary.freq = summary.freq, ...))
  }
  else if(type %in% c("FPAR", "AbsorbedSWRFraction")) {
    OM = x$Plants[[type]][,cohorts,drop=FALSE]
    return(.plot_plant_om(OM, PlantsLAIlive, spnames,
                          type, bySpecies = bySpecies, dates = dates, 
                          xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                          summary.freq = summary.freq, ...))
  } 
  else if(type %in% c("StemPLC", "StemRWC", "LeafRWC")) {
    OM = x$Plants[[type]][,cohorts,drop=FALSE]*100
    return(.plot_plant_om(OM, PlantsLAIlive, spnames,
                   type, bySpecies = bySpecies, dates = dates, 
                   xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                   summary.freq = summary.freq, ...))
  } 
  else if(type %in% c("PlantLAI","PlantLAIlive","PlantTranspiration","PlantNetPhotosynthesis", "PlantGrossPhotosynthesis",
                      "PlantAbsorbedSWR","PlantNetLWR")) {
    subtype = substr(type,6,nchar(type))
    OM = x$Plants[[subtype]][,cohorts,drop=FALSE]
    return(.plot_plant_om_sum(OM, spnames,
                          type, bySpecies = bySpecies, dates = dates, 
                          xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                          summary.freq = summary.freq, ...))
  } 
  else if(type %in% c("SoilPlantConductance","PlantPsi", "LeafPsiMin",
                      "LeafPsiMax", "StemPsi", "RootPsi", "PlantStress",
                      "PlantWaterBalance", "LFMC")) {
    if(type=="SoilPlantConductance") OM = x$Plants[["dEdP"]][,cohorts,drop=FALSE]
    else OM = x$Plants[[type]][,cohorts,drop=FALSE]
    return(.plot_plant_om(OM, PlantsLAIlive, spnames,
                   type, bySpecies = bySpecies, dates = dates, 
                   xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                   summary.freq = summary.freq, ...))
  } 
  else if(type %in% c("LeafPsiMin_SL", "LeafPsiMax_SL", "GSWMin_SL", "GSWMax_SL", "TempMin_SL", "TempMax_SL")) {
    subType = strsplit(type,"_")[[1]][1]
    OM = x$SunlitLeaves[[subType]][,cohorts,drop=FALSE]
    return(.plot_plant_om(OM, PlantsLAIlive, spnames,
                   type, bySpecies = bySpecies, dates = dates, 
                   xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                   summary.freq = summary.freq, ...))
  } 
  else if(type %in% c("LeafPsiMin_SH", "LeafPsiMax_SH", "GSWMin_SH", "GSWMax_SH", "TempMin_SH", "TempMax_SH")) {
    subType = strsplit(type,"_")[[1]][1]
    OM = x$ShadeLeaves[[subType]][,cohorts,drop=FALSE]
    return(.plot_plant_om(OM, PlantsLAIlive, spnames,
                   type, bySpecies = bySpecies, dates = dates, 
                   xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                   summary.freq = summary.freq, ...))
  } 
  else if(type %in% c("TranspirationPerLeaf","GrossPhotosynthesisPerLeaf", "NetPhotosynthesisPerLeaf",
                      "AbsorbedSWRPerLeaf", "NetLWRPerLeaf")) {
    subtype = substr(type, 1, nchar(type)-7)
    df = x$Plants[[subtype]][,cohorts,drop=FALSE]
    df = df/PlantsLAI
    df[PlantsLAI==0] = NA
    return(.plot_plant_om(df, PlantsLAIlive, spnames,
                   type, bySpecies = bySpecies, dates = dates, 
                   xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                   summary.freq = summary.freq, ...))
  } 
  else if(type == "LeafPsiRange") {
    OM1 = x$Plants$LeafPsiMax[,cohorts,drop=FALSE]
    OM2 = x$Plants$LeafPsiMin[,cohorts,drop=FALSE]
    if(bySpecies) {
      OM1 = .averageByLAISpecies(OM1, PlantsLAIlive, spnames)
      OM2 = .averageByLAISpecies(OM2, PlantsLAIlive, spnames)
    } 
    if(!is.null(dates)) {
      OM1 = OM1[row.names(OM1) %in% as.character(dates),,drop = FALSE]
      OM2 = OM2[row.names(OM2) %in% as.character(dates),,drop = FALSE]
    }
    if(!is.null(summary.freq)) {
      OM1 = .temporalSummary(OM1, summary.freq, mean, na.rm=TRUE)
      OM2 = .temporalSummary(OM2, summary.freq, mean, na.rm=TRUE)
    }
    return(.multiple_dynamics_range(as.matrix(OM1), as.matrix(OM2),  xlab = xlab, ylab = ylab, ylim = ylim))
  }
  else if(type=="PlantWUE") {
    OM = x$Plants$GrossPhotosynthesis/x$Plants$Transpiration
    OM[x$Plants$Transpiration==0] = 0
    OM = OM[,cohorts,drop=FALSE]
    return(.plot_plant_om(OM, PlantsLAIlive, spnames,
                          type, bySpecies = bySpecies, dates = dates, 
                          xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                          summary.freq = summary.freq, ...))
  } 
  else if(type %in% c("Temperature","TemperatureRange", "AirTemperature",
                      "CanopyTemperature", "SoilTemperature")) {
    return(.plot_temperature(x$Temperature, type,  
                                dates = dates, 
                                xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                                summary.freq = summary.freq, ...))
  }
  else if(type %in% c("CanopyEnergyBalance", "SoilEnergyBalance")) {
    return(.plot_energybalance(x$EnergyBalance, type,  
                               dates = dates, 
                               xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                               summary.freq = summary.freq, ...))
  } 
}

#' @rdname plot.spwb
plot.growth<-function(x, type="PET_Precipitation", cohorts = NULL, bySpecies = FALSE, 
                      dates = NULL, subdaily = FALSE, 
                      xlim = NULL, ylim=NULL, xlab=NULL, ylab=NULL, 
                      summary.freq = NULL, ...) {
  
  if(subdaily) return(.plotsubdaily(x,type, cohorts, bySpecies, dates, 
                                    xlim, ylim, xlab, ylab))
  
  # Get common elements
  input = x$growthInput
  nlayers = length(input$soil$W)
  transpMode = input$control$transpirationMode
  
  TYPES_GROWTH = .getDailyGROWTHPlotTypes(transpMode)
  TYPES_SWB = .getDailySPWBPlotTypes(transpMode)  

  type = match.arg(type,TYPES_GROWTH)  
  
  if(is.null(cohorts))  cohorts = row.names(input$cohorts)
  else if(length(cohorts)==1) {
    if(cohorts=="T") {
      cohorts = row.names(input$cohorts)
      cohorts = cohorts[substr(cohorts,1,1)=="T"]
    } else if(cohorts=="S") {
      cohorts = row.names(input$cohorts)
      cohorts = cohorts[substr(cohorts,1,1)=="S"]
    }
  }

  spnames = as.character(input$cohorts[cohorts,"Name"])
  PlantsLAIlive = x$Plants$LAIlive[,cohorts, drop=FALSE]
  PlantsLAI = x$Plants$LAI[,cohorts, drop=FALSE]
  
  if(type %in% TYPES_SWB) {
    return(plot.spwb(x,type, cohorts, bySpecies, dates, subdaily, xlim, ylim, xlab, ylab, 
              summary.freq, ...))
  } 
  else if(type == "BiomassBalance") {
    return(.plot_biomass(BiomassBalance= x$BiomassBalance,
                       type = type, dates = dates, 
                       xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                       summary.freq = summary.freq, ...))
  }
  else if(type == "CarbonBalance") {
    return(.plot_carbon(CarbonBalance= x$CarbonBalance,
                         type = type, dates = dates, 
                         xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                         summary.freq = summary.freq, ...))
  }
  else if(type %in% c("GrossPhotosynthesis", "MaintenanceRespiration",  "GrowthCosts", 
                      "LabileCarbonBalance", 
                      "SugarLeaf","StarchLeaf","SugarSapwood","StarchSapwood", "SugarTransport", "RootExudation")) {
      OM = x$LabileCarbonBalance[[type]][,cohorts,drop=FALSE]
  } 
  else if(type == "PhotosynthesisMaintenanceRatio") {
    OM1 = x$LabileCarbonBalance[["GrossPhotosynthesis"]][,cohorts,drop=FALSE]
    OM2 = x$LabileCarbonBalance[["MaintenanceRespiration"]][,cohorts,drop=FALSE]
    OM = OM1/OM2
  } 
  else if(type %in% c("StructuralBiomassBalance","LabileBiomassBalance", "PlantBiomassBalance",
                      "MortalityBiomassLoss","CohortBiomassBalance")) {
    OM = x$PlantBiomassBalance[[type]][,cohorts,drop=FALSE]
  }
  else if(type %in% c("SapwoodBiomass", "LeafBiomass", "FineRootBiomass", 
                      "SapwoodArea", "LeafArea", "FineRootArea", 
                      "DBH", "Height", "HuberValue", "RootAreaLeafArea")) {
    OM = x$PlantStructure[[type]][,cohorts,drop=FALSE]
  } 
  else if(type %in% c("SAgrowth", "LAgrowth", "FRAgrowth", "StarvationRate", "DessicationRate", "MortalityRate")) {
    OM = x$GrowthMortality[[type]][,cohorts,drop=FALSE]
  } 
  return(.plot_plant_om(OM, PlantsLAIlive, spnames,
                        type, bySpecies = bySpecies, dates = dates, 
                        xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                        summary.freq = summary.freq, ...))
}

#' @rdname plot.spwb
plot.fordyn<-function(x, type="StandBasalArea",  
                      cohorts = NULL, bySpecies = FALSE, dates = NULL, 
                      xlim = NULL, ylim=NULL, xlab = NULL, ylab=NULL, 
                      summary.freq = NULL,...) {
  
  input_control = x$GrowthResults[[1]]$growthInput$control
  TYPES_GROWTH = .getDailyGROWTHPlotTypes(input_control$transpirationMode)
  TYPES_GROWTH_UNIQUE = .getUniqueDailyGROWTHPlotTypes(input_control$transpirationMode)
  TYPES_FORDYN_UNIQUE = .getUniqueFORDYNPlotTypes(input_control$transpirationMode)  
  
  type = match.arg(type,c(TYPES_GROWTH, TYPES_FORDYN_UNIQUE))  
  
  if(type %in% TYPES_GROWTH) {
    
    if(is.null(cohorts))  cohorts = row.names(x$GrowthResults[[1]]$growthInput$cohorts)
    PlantsLAI = summary(x, freq = "days", output = "Plants$LAI")[,cohorts,drop=FALSE]
    PlantsLAIlive = summary(x, freq = "days", output = "Plants$LAIlive")[,cohorts,drop=FALSE]
    spnames = as.character(x$GrowthResults[[1]]$growthInput$cohorts[cohorts,"Name"])
    input_soil = x$GrowthResults[[1]]$growthInput$soil
    
    if(type %in% c("PET_Precipitation", "PET_NetRain", "Evapotranspiration" ,"Snow",
                   "WTD", "Export", "SoilVol")) {
      input_soil = x$GrowthResults[[1]]$growthInput$soil
      WaterBalance = summary(x, freq = "days", output = "WaterBalance")
      Soil = summary(x, freq = "days", output = "Soil")
      return(.plot_wb(WaterBalance = WaterBalance, Soil = Soil, input_soil = input_soil,
                      type=type, dates = dates, ylim = ylim, xlab = xlab, ylab = ylab,
                      summary.freq = summary.freq,...))
    }
    else if(type %in% c("SoilPsi", "SoilTheta", "SoilRWC", "PlantExtraction", "HydraulicRedistribution")) {
      Soil = summary(x, freq = "days", output = "Soil")
      return(.plot_soil(Soil = Soil, input_soil = input_soil, input_control = input_control,
                        type=type, dates = dates, 
                        xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab,
                        summary.freq = summary.freq, ...))
    }
    else if(type %in% c("LAI", "GroundIrradiance")) {
      Stand = summary(x, freq = "days", output = "Stand")
      return(.plot_stand(Stand = Stand,
                         type = type, dates = dates, 
                         xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab,
                         summary.freq = summary.freq,...))
    }
    else if(type %in% c("StemPLC", "StemRWC", "LeafRWC", "StemSympRWC", "LeafSympRWC")) {
      OM = summary(x, freq = "days", output = paste0("Plants$",type))[,cohorts,drop=FALSE]*100
      return(.plot_plant_om(OM, PlantsLAIlive, spnames,
                            type, bySpecies = bySpecies, dates = dates, 
                            xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                            summary.freq = summary.freq, ...))
    } 
    else if(type %in% c("PlantLAI","PlantLAIlive","PlantTranspiration","PlantNetPhotosynthesis", "PlantGrossPhotosynthesis",
                        "PlantAbsorbedSWR","PlantNetLWR")) {
      subtype = substr(type,6,nchar(type))
      OM = summary(x, freq = "days", output = paste0("Plants$",subtype))[,cohorts,drop=FALSE]
      return(.plot_plant_om_sum(OM, spnames,
                                type, bySpecies = bySpecies, dates = dates, 
                                xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                                summary.freq = summary.freq, ...))
    } 
    else if(type=="PlantWUE") {
      GP = summary(x, freq = "days", output = "Plants$GrossPhotosynthesis")[,cohorts,drop=FALSE]
      E = summary(x, freq = "days", output = "Plants$Transpiration")[,cohorts,drop=FALSE]
      OM = GP/E
      OM[E==0] = 0
      return(.plot_plant_om(OM, PlantsLAIlive, spnames,
                            type, bySpecies = bySpecies, dates = dates, 
                            xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                            summary.freq = summary.freq, ...))
    } 
    else if(type %in% c("SoilPlantConductance","PlantPsi", "LeafPsiMin",
                        "LeafPsiMax", "StemPsi", "RootPsi", "PlantStress",
                        "PlantWaterBalance", "FPAR", "AbsorbedSWRFraction")) {
      if(type=="SoilPlantConductance") OM = summary(x, freq = "days", output = "Plants$dEdP")[,cohorts,drop=FALSE]
      else OM = summary(x, freq = "days", output = paste0("Plants$",type))[,cohorts,drop=FALSE]
      return(.plot_plant_om(OM, PlantsLAIlive, spnames,
                            type, bySpecies = bySpecies, dates = dates, 
                            xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                            summary.freq = summary.freq, ...))
    } 
    else if(type == "LeafPsiRange") {
      OM1 = summary(x, freq = "days", output = "Plants$LeafPsiMax")[,cohorts,drop=FALSE]
      OM2 = summary(x, freq = "days", output = "Plants$LeafPsiMin")[,cohorts,drop=FALSE]
      if(bySpecies) {
        OM1 = .averageByLAISpecies(OM1, PlantsLAIlive, spnames)
        OM2 = .averageByLAISpecies(OM2, PlantsLAIlive, spnames)
      } 
      if(!is.null(dates)) {
        OM1 = OM1[row.names(OM1) %in% as.character(dates),,drop = FALSE]
        OM2 = OM2[row.names(OM2) %in% as.character(dates),,drop = FALSE]
      }
      if(!is.null(summary.freq)) {
        OM1 = .temporalSummary(OM1, summary.freq, mean, na.rm=TRUE)
        OM2 = .temporalSummary(OM2, summary.freq, mean, na.rm=TRUE)
      }
      return(.multiple_dynamics_range(as.matrix(OM1), as.matrix(OM2),  xlab = xlab, ylab = ylab, ylim = ylim))
    }
    else if(type %in% c("TranspirationPerLeaf","GrossPhotosynthesisPerLeaf", "NetPhotosynthesisPerLeaf",
                        "AbsorbedSWRPerLeaf", "NetLWRPerLeaf")) {
      subtype = substr(type, 1, nchar(type)-7)
      OM = summary(x, freq = "days", output = paste0("Plants$",subtype))[,cohorts,drop=FALSE]
      OM = OM/PlantsLAI
      OM[PlantsLAI==0] = NA
      return(.plot_plant_om(OM, PlantsLAIlive, spnames,
                            type, bySpecies = bySpecies, dates = dates, 
                            xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                            summary.freq = summary.freq, ...))
    } 
    else if(type %in% c("Temperature","TemperatureRange", "AirTemperature",
                        "CanopyTemperature", "SoilTemperature")) {
      OM = summary(x, freq = "days", output = "Temperature")
      return(.plot_temperature(OM, type,  
                               dates = dates, 
                               xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                               summary.freq = summary.freq, ...))
    }
    else if(type %in% c("CanopyEnergyBalance", "SoilEnergyBalance")) {
      OM = summary(x, freq = "days", output = "EnergyBalance")
      return(.plot_energybalance(OM, type,  
                                 dates = dates, 
                                 xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                                 summary.freq = summary.freq, ...))
    } 
    else if(type %in% TYPES_GROWTH_UNIQUE) {
      if(type=="BiomassBalance") {
        OM = summary(x, freq = "days", output = "BiomassBalance")
        return(.plot_biomass(OM, type,  
                             dates = dates, 
                             xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                             summary.freq = summary.freq, ...))
      }
      if(type=="CarbonBalance") {
        OM = summary(x, freq = "days", output = "CarbonBalance")
        return(.plot_carbon(OM, type,  
                             dates = dates, 
                             xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                             summary.freq = summary.freq, ...))
      }
      if(type %in% c("GrossPhotosynthesis", "MaintenanceRespiration",  "GrowthCosts", 
                     "LabileCarbonBalance", 
                     "SugarLeaf","StarchLeaf","SugarSapwood","StarchSapwood", "SugarTransport", "RootExudation")) {
        OM = summary(x, freq = "days", output = paste0("LabileCarbonBalance$",type))[,cohorts,drop=FALSE]
      } 
      if(type =="PhotosynthesisMaintenanceRatio") {
        OM1 = summary(x, freq = "days", output = "LabileCarbonBalance$GrossPhotosynthesis")[,cohorts,drop=FALSE]
        OM2 = summary(x, freq = "days", output = "LabileCarbonBalance$MaintenanceRespiration")[,cohorts,drop=FALSE]
        OM = OM1/OM2
      }
      else if(type %in% c("StructuralBiomassBalance", "LabileBiomassBalance", "PlantBiomassBalance", 
                          "MortalityBiomassLoss", "CohortBiomassBalance")) {
        OM = summary(x, freq = "days", output = paste0("PlantBiomassBalance$",type))[,cohorts,drop=FALSE]
      } 
      else if(type %in% c("LeafBiomass", "SapwoodBiomass", "FineRootBiomass", "SapwoodArea", "LeafArea", "FineRootArea","DBH", "Height",
                          "HuberValue", "RootAreaLeafArea")) {
        OM = summary(x, freq = "days", output = paste0("PlantStructure$",type))[,cohorts,drop=FALSE]
      } 
      else if(type %in% c("SAgrowth", "LAgrowth", "FRAgrowth", "StarvationRate", "MortalityRate", "DessicationRate")) {
        OM = summary(x, freq = "days", output = paste0("GrowthMortality$",type))[,cohorts,drop=FALSE]
      } 
      return(.plot_plant_om(OM, PlantsLAIlive, spnames,
                            type, bySpecies = bySpecies, dates = dates, 
                            xlim = xlim, ylim=ylim, xlab=xlab, ylab=ylab, 
                            summary.freq = summary.freq, ...))
    }
  }


  ## FORDYN PLOT

  if(is.null(ylab)) ylab = .getYLab(type)
  if(is.null(xlab)) xlab = "Step"

  if(type %in% c("NumTreeSpecies", "NumTreeCohorts", "NumShrubSpecies", "NumShrubCohorts"))  {
    out = x$StandSummary
    var = type
  }
  else if(type == "StandDensity")  {
    out = x$StandSummary
    var = "TreeDensityLive"
  }
  else if(type == "StandBasalArea")  {
    out = x$StandSummary
    var = "TreeBasalAreaLive"
  }
  else if(type == "StandShrubCover")  {
    out = x$StandSummary
    var = "ShrubCoverLive"
  }
  else if(type %in% c("DominantTreeHeight", "DominantTreeDiameter", "QuadraticMeanTreeDiameter", "HartBeckingIndex")) {
    out = x$StandSummary
    var = type
  }
  else if(type == "SpeciesDensity")  {
    out = x$SpeciesSummary
    var = "TreeDensityLive"
  }
  else if(type == "SpeciesBasalArea")  {
    out = x$SpeciesSummary
    var = "TreeBasalAreaLive"
  } 
  else if(type == "SpeciesShrubCover")  {
    out = x$SpeciesSummary
    var = "ShrubCoverLive"
  } 
  else if(type == "NumCohortsSpecies")  {
    out = x$SpeciesSummary
    var = "NumCohorts"
  } 
  else if(type == "CohortDensity")  {
    out = x$CohortSummary
    var = "TreeDensityLive"
  }
  else if(type == "CohortBasalArea")  {
    out = x$CohortSummary
    var = "TreeBasalAreaLive"
  }
  else if(type == "CohortShrubCover")  {
    out = x$CohortSummary
    var = "ShrubCoverLive"
  }
  
  df = data.frame(Step = out[["Step"]], y = out[[var]])
  if(type %in% c("SpeciesDensity", "SpeciesBasalArea", "SpeciesShrubCover", "NumCohortsSpecies")) df$group = as.character(out[["Species"]])
  else if(type %in% c("CohortBasalArea", "CohortDensity", "CohortShrubCover")) {
    df$group = paste0(as.character(out[["Cohort"]]), " (", as.character(out[["Species"]]),")")
  }
  df = df[!is.na(df$y),]
  g<-ggplot(df, aes(x=.data$Step, y=.data$y))
  if("group" %in% names(df)) {
    g <- g+ geom_line(aes(col=.data$group))+
      scale_color_discrete(name="")
  } else {
    g <- g+ geom_line()
  }
  if(!is.null(ylim)) g <- g+ylim(ylim)
  g<-g+theme_bw()+ylab(ylab)+xlab(xlab)
  return(g)
}