---
title: "Forest dynamics"
author: "Miquel De Caceres"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    toc: TRUE
vignette: >
  %\VignetteIndexEntry{Forest dynamics}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignettePackage{medfate}
  \usepackage[utf8]{inputenc}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(medfate)
```

## About this vignette

This document describes how to run the forest dynamics model of `medfate`. This document is meant to teach users to run the simulation model within R. Details of the model design and formulation can be found at https://emf-creaf.github.io/medfatebook/index.html. 

## Preparing model inputs
Any forest dynamics model needs information on climate, vegetation and soils of the forest stand to be simulated. Moreover, since models in `medfate` differentiate between species, information on species-specific model parameters is also needed. In this subsection we explain the different steps to prepare the data needed to run function `fordyn()`.

## Preparing model inputs
Model inputs are explained in greater detail in vignette '*Simulation inputs*'. Here we only review the different steps required to run function `fordyn()`.

### Soil, vegetation, meteorology and species data

Soil information needs to be entered as a `data frame` with soil layers in rows and physical attributes in columns. Soil physical attributes can be initialized to default values, for a given number of layers, using function `defaultSoilParams()`:

```{r}
spar = defaultSoilParams(2)
```

The soil input for water balance simulation is actually a list of class `soil` that is created using a function with the same name:
```{r}
examplesoil = soil(spar)
```

As explained in the package overview, models included in `medfate` were primarily designed to be ran on **forest inventory plots**. Here we use the example object provided with the package:
```{r}
data(exampleforestMED)
exampleforestMED
```

In the basic water balance, only mean temperature, precipitation and potential evapotranspiration is required, but radiation may also be necessary to simulate snow melt. 
```{r}
data(examplemeteo)
head(examplemeteo)
```

Finally, simulations in `medfate` require a data frame with species parameter values, which we load using defaults for Catalonia (NE Spain):

```{r}
data("SpParamsMED")
```

### Simulation control

Apart from data inputs, the behaviour of simulation models can be controlled using a set of global parameters. The default parameterization is obtained using function `defaultControl()`:

```{r}
control = defaultControl("Granier")
```
Here we will run simulations of forest dynamics using the basic water balance model (i.e. `transpirationMode = "Granier"`). The complexity of the soil water balance calculations can be changed by using `"Sperry"` as input to `defaultControl()`. However, when running `fordyn()` sub-daily output will never be stored (i.e. setting `subdailyResults = TRUE` is useless).

## Executing the forest dynamics model

In this vignette we will fake a three-year weather input by repeating the example weather data frame three times, while reducing precipitation of the second year to the half.
```{r}
meteo2001 = examplemeteo
meteo2002 = examplemeteo
meteo2002$Precipitation = meteo2002$Precipitation/2
meteo2003 = examplemeteo
row.names(meteo2002) = seq(as.Date("2002-01-01"), 
                           as.Date("2002-12-31"), by="day")
row.names(meteo2003) = seq(as.Date("2003-01-01"), 
                           as.Date("2003-12-31"), by="day")
meteo_01_03 = rbind(meteo2001, meteo2002, meteo2003)
```


Now we run the forest dynamics model using all inputs (note that no intermediate input object is needed, as in `spwb()` or  `growth()`):
```{r}
fd<-fordyn(exampleforestMED, examplesoil, SpParamsMED, meteo_01_03, control, 
           latitude = 41.82592, elevation = 100)
```
It is worth noting that, while `fordyn()` calls function `growth()` internally for each simulated year, the `verbose` option of the control parameters only affects function `fordyn()` (i.e. all console output from `growth()` is hidden). Recruitment and summaries are done only once a year at the level of function `fordyn()`.

## Inspecting model outputs

### Stand, species and cohort summaries

Among other outputs, function `fordyn()` calculates standard summary statistics that describe the structural and compositional state of the forest at each time step. For example, we can access stand-level statistics using:

```{r}
fd$StandSummary
```

where we can observe an increase in stand basal area and leaf area index during years 1 and 3, but a reduction during the second (drier) year. Species-level analogous statistics are shown using:
```{r}
fd$SpeciesSummary
```

### Tree/shrub tables

Another useful output of `fordyn()` are tables in long format with cohort structural information (i.e. DBH, height, density, etc) for each time step:

```{r}
fd$TreeTable
```

The same can be shown for dead trees, where we see larger mortality rates during the second (drier) year:
```{r}
fd$DeadTreeTable
```

### Accessing the output from function growth()

Since function `fordyn()` makes internal calls to function `growth()`, it stores the result in a vector called `GrowthResults`, which we can use to inspect intra-annual patterns of desired variables. For example, the following shows the leaf area for individuals of the three cohorts during the three consecutive years:

```{r, fig=TRUE, echo=TRUE, fig.width=7, fig.height=3}
plot(fd$GrowthResults[[1]], "LeafArea", bySpecies = T)
plot(fd$GrowthResults[[2]], "LeafArea", bySpecies = T)
plot(fd$GrowthResults[[3]], "LeafArea", bySpecies = T)
```
where we see that the second (drier) year, resulted in a decrease in leaf area, due to the sink limitation caused by drought, and the third (normal) year lead to a recovery of leaf area. Note also the sink limitation to sapwood area (hence diameter) growth during the second year:

```{r, fig=TRUE, echo=TRUE, fig.width=7, fig.height=3}
plot(fd$GrowthResults[[1]], "SapwoodArea", bySpecies = T)
plot(fd$GrowthResults[[2]], "SapwoodArea", bySpecies = T)
plot(fd$GrowthResults[[3]], "SapwoodArea", bySpecies = T)
```
