---
title: "Forest growth"
author: "Miquel De Caceres"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    toc: TRUE
vignette: >
  %\VignetteIndexEntry{Forest growth}
  %\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 growth model of `medfate`, implemented in function `growth()`. All the details of the model design and formulation can be found at https://emf-creaf.github.io/medfatebook/index.html. 

## 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 `growth()`.

### 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 growth simulations 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()`. 


### Growth input object

A last object is needed before calling simulation functions, called `growthInput`. This is analogous to `spwbInput` and consists in the compilation of aboveground and belowground parameters and the specification of additional parameter values for each plant cohort. The object can be obtained by calling function `growthInput()`, but if one has a `forest` object, it can be generated more directly using function `forest2growthInput()`:

```{r}
x = forest2growthInput(exampleforestMED, examplesoil, SpParamsMED, control)
```

All the input information for forest data and species parameter values can be inspected by printing different elements of the input object, which are: 
```{r}
names(x)
```
First, information about the cohort species is found in element `cohorts`  (i.e. code, species and name):
```{r}
x$cohorts
```

Element `above` contains the aboveground structure data that we already know, but with additional columns that describe the estimated initial amount of sapwood area and carbon reserves (if required):
```{r}
x$above
```

As with the soil input object, the `growthInput` object will be modified during simulations. In the case of `growth()`, these modifications concern `SA`, `fastCstorage`, `slowCstorage`, `LAI_expanded` and `LAI_dead` in element `above`, as well as canopy variables.

Elements starting with `params*` contain cohort-specific model parameters. Aboveground parameters related to plant transpiration can be seen in `paramsTransp`:
```{r}
x$paramsTransp
```

Belowground parameters can be seen in `below` and include root distribution as well as maximum root and rhizosphere conductances by soil layer:
```{r}
x$below
```

Parameters related to plant growth can be seen in `paramsGrowth` and `paramsAllometries`:
```{r}
x$paramsGrowth
x$paramsAllometries
```

Finally, elements starting with `internal*` contain state variables required to keep track of plant status. For example, the carbon storage levels can be seen in `internalCarbon`:
```{r}
x$internalCarbon
```

## Executing the growth model

In this vignette we will use one year of example meteorological data with the following precipitation and temperature seasonal trends:
```{r}
G1<-growth(x, examplemeteo, latitude = 41.82592, elevation = 100)

```

While water/carbon balances and growth of leaf/sapwood/root tissues is performed on a daily basis, structural variables are only updated once a year, here coinciding with the end of the simulation.

## Plotting results

In addition to all the plots that were available for result of simulations with function `spwb()`, one can display the following information out of the output of function `growth()`:

### Carbon balance components
To inspect components of the plant carbon balance we can first display gross photosynthesis expressed as assimilation relative to biomass:
```{r, fig=TRUE, echo=TRUE, fig.width=7, fig.height=3}
plot(G1, "GrossPhotosynthesis", bySpecies = T)
```
Then we can draw the respiration costs (maintenance and growth) in the same units:
```{r, fig=TRUE, echo=TRUE, fig.width=7, fig.height=3}
plot(G1, "MaintenanceRespiration", bySpecies = T)
plot(G1, "GrowthCosts", bySpecies = T)
```
Finally we can display the daily negative or positive balance of the plant, which determines changes in plant carbon pools:
```{r, fig=TRUE, echo=TRUE, fig.width=7, fig.height=3}
plot(G1, "CarbonBalance", bySpecies = T)
```

### Carbon storage
Carbon assimilation and respiration rates define the dynamics of stored carbon, which in medfate is divided into two organs (leaves and sapwood) and two chemical compounds (soluble sugars and starch):
```{r, fig=TRUE, echo=TRUE, fig.width=7, fig.height=3}
plot(G1, "SugarLeaf", bySpecies = T)
plot(G1, "StarchLeaf", bySpecies = T)
```

```{r, fig=TRUE, echo=TRUE, fig.width=7, fig.height=3}
plot(G1, "SugarSapwood", bySpecies = T)
plot(G1, "StarchSapwood", bySpecies = T)
```


### Leaf and sapwood area changes
Leaf and sapwood area dynamics arising from the interplay between growth and senescence of tissues can be inspected using:
```{r, fig=TRUE, echo=TRUE, fig.width=7, fig.height=3}
plot(G1, "LeafArea", bySpecies = T)
plot(G1, "SapwoodArea", bySpecies = T)
```
The corresponding daily growth rates can also be displayed:
```{r, fig=TRUE, echo=TRUE, fig.width=7, fig.height=3}
plot(G1, "LAgrowth", bySpecies = T)
plot(G1, "SAgrowth", bySpecies = T)
```
