baseline_model.Rmd

---
title: "QE framework result summary"
author: "Mingkai Jiang"
date: "05-06-2018"
output:
  html_document:
    toc: true
    toc_float: true
    theme: united
---

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

# 1. Introduction

This R markdown file introduces the Quasi-equilibrium (QE) framework baseline model. The basline QE model is a carbon-nitrogen (C-N) coupled model that is based on the Generic Decomposition And Yield (GDAY) model, originally developed by Comins and McMurtrie (1993).

The baseline model in this study compares the effect of wood NC stoichiometric flexiblity, with the 1st model has variable wood NC ratio and 2nd model has fixed wood NC ratio. 

Here we are evaluating the alternative assumption for wood stoichiometric flexiblity on plant's CO2 response. Time is also a factor considered here, with VL (very-long), L (long), and M (medium) terms considered as the final model equilibrium point, passive SOM pool equilibrated, passive and slow SOM pools equilibrated, respectively. 

<br/><br/>

# 2. Model structure

<p align="center">
<img src="Model_framework.png" width = "650" height = "650">
</p>

The entire QE framework is built based on the model structure of the Generic Decomposition And Yield (G'DAY) model, as shown in the figure above. Boxes with dotted boundaries are the M term pools, box filled with diamonds is the L term pool. 

<br/><br/>

# 3. Procedures

Below I use step-by-step guidance to generate the results.

``` {r, message=FALSE}
    source("QE_Analyses/prepare_R.R")

```

Next, we run analytical run scripts to generate the pdf files and store the rest dataframes into data list:

``` {r, message = FALSE,warning=F}
    source("QE_Analyses/Run_analytical_solutions.R")

```

Here I am preparing the necessary color and plot configurations. 

``` {r, message = FALSE, warning=FALSE, results="hide"}
    # library
    require(ggplot2)
    require(knitr)

    # two nice color palette for color blind
    # The palette with grey:
    cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

    # The palette with black:
    cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
    
    save.plot <- T

```


Now, I am plotting the graphic illustration of the QE framework, based on Comins and McMurtrie's old approach.

``` {r, echo = FALSE, fig.width = 7, fig.align="center", message = FALSE, warning = FALSE, results="hide"}
    bDF <- mcm_model()

    ### split the data.list into two data frames
    cdf <- bDF$cDF
    edf <- bDF$eDF

    ### plot
    p<-ggplot() + 
        geom_line(data=cdf, aes(x=nc, y=P350, col="C400"), size=0.5, linetype="dotted")+   
        geom_line(data=cdf, aes(x=nc, y=VL, col="VL"), size=0.5, linetype="dotted") + 
        geom_line(data=cdf, aes(x=nc, y=L, col="L"), size=0.5, linetype="dotted") +  
        geom_line(data=cdf, aes(x=nc, y=M, col="M"), size=0.5, linetype="dotted") +            
        geom_line(data=cdf, aes(x=nc, y=P700, col="C800"), size=0.5, linetype="dotted") +     
        geom_point(data=edf, aes(x=nc_VL, y=NPP_VL, fill="A"), 
                color="black", shape=21, size=5) + 
        geom_point(data=edf, aes(x=nc_L_700, y=NPP_L_700, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf, aes(x=nc_VL_700, y=NPP_VL_700, fill="E"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf, aes(x=nc_M_700, y=NPP_M_700, fill="C"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf, aes(x=nc_I, y=NPP_I, fill="B"), 
                   shape=21, color="black", size=5) +
        ylim(1.0, 4.0) + 
        xlim(0.001, 0.04) +
        labs(x="Leaf N:C Ratio", 
             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=12),
              axis.title=element_text(size=14),
              legend.text=element_text(size=12),
              legend.title=element_text(size=14),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = cbbPalette[7],
                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4], "E" = cbPalette[7])) 
        ### only save this to pdf and then add arrows and recall it
        pdf("Plots/QE_illustration.pdf", width=10, height=10)
        plot(p)
        dev.off()

```

<p align="center">
<img src="Transient_behavior.png" width = "650" height = "650">
</p>


<br/><br/>

Now, I am plotting QE framework for each individual model assumption.

### Run 1: Baseline N model with variable wood NC ratio
``` {r, echo = FALSE, fig.width = 7, fig.align="center", message = FALSE, warning = FALSE, results="hide"}
    ### split the data.list into two data frames
    cdf.r1 <- r1$cDF
    edf.r1 <- r1$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r1 <- cdf.r1[cdf.r1$CO2 == 400, ]
    cdf2.r1 <- cdf.r1[cdf.r1$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r1 <- edf.r1[edf.r1$CO2 == 400, ]
    edf2.r1 <- edf.r1[edf.r1$CO2 == 800, ]
    ### plot
    p<-ggplot() + 
        geom_line(data=cdf1.r1, aes(x=nc, y=NPP_photo, col="C400")) +   
        geom_line(data=cdf1.r1, aes(x=nc, y=NPP_VL, col="VL")) + 
        geom_line(data=cdf1.r1, aes(x=nc, y=NPP_L, col="L")) +  
        geom_line(data=cdf1.r1, aes(x=nc, y=NPP_M, col="M")) +            
        geom_line(data=cdf2.r1, aes(x=nc, y=NPP_photo, col="C800")) +     
        geom_point(data=edf1.r1, aes(x=nc_VL, y=NPP_VL, fill="A"), 
                color="black", shape=21, size=5) + 
        geom_point(data=edf2.r1, aes(x=nc_L, y=NPP_L, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r1, aes(x=nc_M, y=NPP_M, fill="C"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r1, aes(x=nc_VL, y=NPP_VL, fill="E"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf1.r1, aes(x=nc_VL, y=NPP_I, fill="B"), 
                   shape=21, color="black", size=5) +
        ylim(1.25, 2.25) + 
        xlim(0.001, 0.03) +
        labs(x="Leaf N:C Ratio", 
             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = cbbPalette[7],
                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4], "E" = cbPalette[7])) 

    plot(p)
    if (save.plot) {
        pdf("Plots/Run01.pdf", width=8, height=6)
        plot(p)
        dev.off()
    } 
```

Here, the QE framework indicates that, the VL and L line constraints are very much similar. 

<br/>
<br/>

### Run 2: Baseline N model with fixed wood NC ratio
```{r, echo=FALSE, fig.width = 7, fig.align="center", message = FALSE, warning = FALSE, results="hide"}

    ### split the data.list into two data frames
    cdf.r2 <- r2$cDF
    edf.r2 <- r2$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r2 <- cdf.r2[cdf.r2$CO2 == 400, ]
    cdf2.r2 <- cdf.r2[cdf.r2$CO2 == 800, ]
    
    ### extract CO2 conc. from eq DF
    edf1.r2 <- edf.r2[edf.r2$CO2 == 400, ]
    edf2.r2 <- edf.r2[edf.r2$CO2 == 800, ]
    
    ### plot
    p<-ggplot() + 
        geom_line(data=cdf1.r2, aes(x=nc, y=NPP_photo, col="C400")) +   
        geom_line(data=cdf1.r2, aes(x=nc, y=NPP_VL, col="VL")) + 
        geom_line(data=cdf1.r2, aes(x=nc, y=NPP_L, col="L")) +  
        geom_line(data=cdf1.r2, aes(x=nc, y=NPP_M, col="M")) +            
        geom_line(data=cdf2.r2, aes(x=nc, y=NPP_photo, col="C800")) +     
        geom_point(data=edf1.r2, aes(x=nc_VL, y=NPP_VL, fill="A"), 
                color="black", shape=21, size=5) + 
        geom_point(data=edf2.r2, aes(x=nc_L, y=NPP_L, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r2, aes(x=nc_M, y=NPP_M, fill="C"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r2, aes(x=nc_VL, y=NPP_VL, fill="E"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf1.r2, aes(x=nc_VL, y=NPP_I, fill="B"), 
                   shape=21, color="black", size=5) +
        ylim(1.25, 2.25) + 
        xlim(0.001, 0.03) +
        labs(x="Leaf N:C Ratio", 
             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = cbbPalette[7],
                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4], "E" = cbPalette[7])) 
    plot(p)
    if (save.plot) {
        pdf("Plots/Run02.pdf", width=8, height=6)
        plot(p)
        dev.off()
    } 
```

Comparing variable (Run 1) and fixed (Run 2) wood NC ratio treatment, two key features emerge:

        - 1. The equilibrium points shift towards lower N:C ratio and NPP regions under fixed wood NC assumption.
        - 2. VL and L line are more separable under fixed wood NC assumption.

Variable wood stoichiometry means that there's less nutrient stored in the wood storage pool, hence more N is recycling at fast timescale for higher NPP production at higher N:C ratio (higher N:C means C is assimilated at higher N availability). 

Percent change in shoot N:C ratio and NPP at I, M, L and VL timepoints are summarized below:

``` {r, echo=FALSE}
    ### compute % change in NC and NPP under eCO2 condition
    chang1 <- (edf2.r1[,2:8] - edf1.r1[,2:8]) / edf1.r1[,2:8] * 100
    chang1$NPP_I <- (edf1.r1[,"NPP_I"] - edf1.r1[,"NPP_VL"]) / edf1.r1[,"NPP_VL"] * 100

    chang2 <- (edf2.r2[,2:8] - edf1.r2[,2:8]) / edf1.r2[,2:8] * 100
    chang2$NPP_I <- (edf1.r2[,"NPP_I"] - edf1.r2[,"NPP_VL"]) / edf1.r2[,"NPP_VL"] * 100
    
    ### organize the df for output
    changeDF <- round(rbind(chang1, chang2),1)
    changeDF[1,"Run"] <- "variable NC" 
    changeDF[2,"Run"] <- "fixed NC" 
    changeDF <- changeDF[,c("Run", "NPP_I", "NPP_M", "NPP_L", "NPP_VL", 
                            "nc_M", "nc_L", "nc_VL")]
    
    ### table output
    kable(changeDF, row.names=F, format = "markdown")

```

Comparing variable to fixed wood stoichiometric treatment, we have less CO2 fertilization effect at the I term but more CO2 fertilization effect at M, L and VL terms. NC ratios for M, L and VL terms decreased more under fixed than variable wood NC assumption. 

Under fixed wood NC scenario, more N is stored in woody pool at the M term, leaving less N available for M-term NPP, thus the significant less M term CO2 fertilization effect. 

<br/>
<br/>  

### Run 3: Explicit mineral N pool, variable wood NC ratio
``` {r, echo = FALSE, fig.width = 7, fig.align="center", message = FALSE, warning = FALSE, results="hide"}
    ### split the data.list into two data frames
    cdf.r3 <- r3$cDF
    edf.r3 <- r3$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r3 <- cdf.r3[cdf.r3$CO2 == 400, ]
    cdf2.r3 <- cdf.r3[cdf.r3$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r3 <- edf.r3[edf.r3$CO2 == 400, ]
    edf2.r3 <- edf.r3[edf.r3$CO2 == 800, ]

    ### plot
    p<-ggplot() + 
        geom_line(data=cdf1.r3, aes(x=nc, y=NPP_photo, col="C400")) +   
        geom_line(data=cdf1.r3, aes(x=nc, y=NPP_VL, col="VL")) + 
        geom_line(data=cdf1.r3, aes(x=nc, y=NPP_L, col="L")) +  
        geom_line(data=cdf1.r3, aes(x=nc, y=NPP_M, col="M")) +            
        geom_line(data=cdf2.r3, aes(x=nc, y=NPP_photo, col="C800")) +     
        geom_point(data=edf1.r3, aes(x=nc_VL, y=NPP_VL, fill="A"), 
                color="black", shape=21, size=5) + 
        geom_point(data=edf2.r3, aes(x=nc_L, y=NPP_L, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r3, aes(x=nc_M, y=NPP_M, fill="C"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r3, aes(x=nc_VL, y=NPP_VL, fill="E"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf1.r3, aes(x=nc_VL, y=NPP_I, fill="B"), 
                   shape=21, color="black", size=5) +
        ylim(1.25, 2.25) + 
        xlim(0.001, 0.03) +
        labs(x="Leaf N:C Ratio", 
             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = cbbPalette[7],
                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4], "E" = cbPalette[7])) 

    plot(p)
    if (save.plot) {
        pdf("Plots/Run03.pdf", width=8, height=6)
        plot(p)
        dev.off()
    } 
```

Given that N uptake rate is assumed to be 1, here this QE framework should be exactly the same as the Run 1 simulation (or at least be visually the same). Visually, there is minor difference comparing this graph to Run 1 graph, but I think this small difference can be neglected. 

<br/>
<br/>

### Run 4: Explicit mineral N pool, fixed wood NC ratio
``` {r, echo = FALSE, fig.width = 7, fig.align="center", message = FALSE, warning = FALSE, results="hide"}
    ### split the data.list into two data frames
    cdf.r4 <- r4$cDF
    edf.r4 <- r4$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r4 <- cdf.r4[cdf.r4$CO2 == 400, ]
    cdf2.r4 <- cdf.r4[cdf.r4$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r4 <- edf.r4[edf.r4$CO2 == 400, ]
    edf2.r4 <- edf.r4[edf.r4$CO2 == 800, ]

    ### plot
    p<-ggplot() + 
        geom_line(data=cdf1.r4, aes(x=nc, y=NPP_photo, col="C400")) +   
        geom_line(data=cdf1.r4, aes(x=nc, y=NPP_VL, col="VL")) + 
        geom_line(data=cdf1.r4, aes(x=nc, y=NPP_L, col="L")) +  
        geom_line(data=cdf1.r4, aes(x=nc, y=NPP_M, col="M")) +            
        geom_line(data=cdf2.r4, aes(x=nc, y=NPP_photo, col="C800")) +     
        geom_point(data=edf1.r4, aes(x=nc_VL, y=NPP_VL, fill="A"), 
                color="black", shape=21, size=5) + 
        geom_point(data=edf2.r4, aes(x=nc_L, y=NPP_L, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r4, aes(x=nc_M, y=NPP_M, fill="C"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r4, aes(x=nc_VL, y=NPP_VL, fill="E"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf1.r4, aes(x=nc_VL, y=NPP_I, fill="B"), 
                   shape=21, color="black", size=5) +
        ylim(1.25, 2.25) + 
        xlim(0.001, 0.03) +
        labs(x="Leaf N:C Ratio", 
             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = cbbPalette[7],
                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4], "E" = cbPalette[7])) 

    plot(p)
    if (save.plot) {
        pdf("Plots/Run04.pdf", width=8, height=6)
        plot(p)
        dev.off()
    } 
```

Given that N uptake rate is assumed to be 1, here this QE framework should be exactly the same as the Run 2 simulation (or at least be visually the same).

Here (run 3 and 4), minral N pool is made explicit. By modifying mineral N uptake rate, we should be able to see different QE points as the uptake of N will change. 

<br/>
<br/>

### Run 5: N uptake as a saturating function of root biomass, variable wood NC (GDAY approach)
``` {r, echo = FALSE, fig.width = 7, fig.align="center", message = FALSE, warning = FALSE, results="hide"}
    ### split the data.list into two data frames
    cdf.r5 <- r5$cDF
    edf.r5 <- r5$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r5 <- cdf.r5[cdf.r5$CO2 == 400, ]
    cdf2.r5 <- cdf.r5[cdf.r5$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r5 <- edf.r5[edf.r5$CO2 == 400, ]
    edf2.r5 <- edf.r5[edf.r5$CO2 == 800, ]

    ### plot
    p<-ggplot() + 
        geom_line(data=cdf1.r5, aes(x=nc, y=NPP_photo, col="C400")) +   
        geom_line(data=cdf1.r5, aes(x=nc, y=NPP_VL, col="VL")) + 
        geom_line(data=cdf1.r5, aes(x=nc, y=NPP_L, col="L")) +  
        geom_line(data=cdf1.r5, aes(x=nc, y=NPP_M, col="M")) +            
        geom_line(data=cdf2.r5, aes(x=nc, y=NPP_photo, col="C800")) +     
        geom_point(data=edf1.r5, aes(x=nc_VL, y=NPP_VL, fill="A"), 
                color="black", shape=21, size=5) + 
        geom_point(data=edf2.r5, aes(x=nc_L, y=NPP_L, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r5, aes(x=nc_M, y=NPP_M, fill="C"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r5, aes(x=nc_VL, y=NPP_VL, fill="E"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf1.r5, aes(x=nc_VL, y=NPP_I, fill="B"), 
                   shape=21, color="black", size=5) +
        ylim(1.0, 2.2) + 
        xlim(0.001, 0.03) +
        labs(x="Leaf N:C Ratio", 
             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = cbbPalette[7],
                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4], "E" = cbPalette[7])) 

    plot(p)
    if (save.plot) {
        pdf("Plots/Run05.pdf", width=8, height=6)
        plot(p)
        dev.off()
    } 
```

Here, uptake rate of N is represented mathematically as a saturating function of root biomass and linearly depend upon mineral N pool. We also assumed variable wood NC ratio to make it comparable to Run 1. 

<br/>
<br/>

### Run 6: N uptake as a saturating function of mineral N, variable wood NC ratio (O-CN approach)
``` {r, echo = FALSE, fig.width = 7, fig.align="center", message = FALSE, warning = FALSE, results="hide"}
    ### split the data.list into two data frames
    cdf.r6 <- r6$cDF
    edf.r6 <- r6$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r6 <- cdf.r6[cdf.r6$CO2 == 400, ]
    cdf2.r6 <- cdf.r6[cdf.r6$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r6 <- edf.r6[edf.r6$CO2 == 400, ]
    edf2.r6 <- edf.r6[edf.r6$CO2 == 800, ]

    ### plot
    p<-ggplot() + 
        geom_line(data=cdf1.r6, aes(x=nc, y=NPP_photo, col="C400")) +   
        geom_line(data=cdf1.r6, aes(x=nc, y=NPP_L, col="L")) +  
        geom_line(data=cdf1.r6, aes(x=nc, y=NPP_M, col="M")) +            
        geom_line(data=cdf2.r6, aes(x=nc, y=NPP_photo, col="C800")) +     
        geom_point(data=edf1.r6, aes(x=nc_VL, y=NPP_VL, fill="A"), 
                   color="black", shape=21, size=5) + 
        geom_point(data=edf2.r6, aes(x=nc_L, y=NPP_L, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r6, aes(x=nc_M, y=NPP_M, fill="C"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r6, aes(x=nc_VL, y=NPP_VL, fill="E"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf1.r6, aes(x=nc_VL, y=NPP_I, fill="B"), 
                   shape=21, color="black", size=1) +
        ylim(1.0, 2.2) + 
        xlim(0.001, 0.03) +
        labs(x="Leaf N:C Ratio", 
             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], 
                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4], "E" = cbPalette[7])) 

    plot(p)
    if (save.plot) {
        pdf("Plots/Run06.pdf", width=8, height=6)
        plot(p)
        dev.off()
    } 
```

Here, mineral N uptake rate is a saturating function of mineral N pool and also linearly depend upon root biomass. The mathematical formular is similar to those represented in O-CN model.

There is no VL constraint curve because NPP cancels out for all possible leaf NC values (not shown here but it's in the equation). This means that, the VL-term CO2 fertilization effect is the same as the initial response, which means it is a full CO2 fertilization response. Over shorter timescales, we the NPP responses are similar, as the table below demonstrated. 

``` {r, echo=FALSE}
    ### compute % change in NC and NPP under eCO2 condition
    chang5 <- (edf2.r5[,2:8] - edf1.r5[,2:8]) / edf1.r5[,2:8] * 100
    chang5$NPP_I <- (edf1.r5[,"NPP_I"] - edf1.r5[,"NPP_VL"]) / edf1.r5[,"NPP_VL"] * 100

    chang6 <- (edf2.r6[,2:8] - edf1.r6[,2:8]) / edf1.r6[,2:8] * 100
    chang6$NPP_I <- (edf1.r6[,"NPP_I"] - edf1.r6[,"NPP_VL"]) / edf1.r6[,"NPP_VL"] * 100
    
    ### organize the df for output
    changeDF <- round(rbind(chang5, chang6),1)
    changeDF[1,"Run"] <- "GDAY" 
    changeDF[2,"Run"] <- "OCN" 
    changeDF <- changeDF[,c("Run", "NPP_I", "NPP_M", "NPP_L", "NPP_VL", 
                            "nc_M", "nc_L", "nc_VL")]
    
    ### table output
    kable(changeDF, row.names=F, format = "markdown")

```

<br/>
<br/>

### Run 7: Priming effect, fixed wood NC ratio
``` {r, echo = FALSE, fig.width = 7, fig.align="center", message = FALSE, warning = FALSE, results="hide"}
    ### split the data.list into two data frames
    cdf.r7 <- r7$cDF
    edf.r7 <- r7$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r7 <- cdf.r7[cdf.r7$CO2 == 400, ]
    cdf2.r7 <- cdf.r7[cdf.r7$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r7 <- edf.r7[edf.r7$CO2 == 400, ]
    edf2.r7 <- edf.r7[edf.r7$CO2 == 800, ]

    ### plot
    p<-ggplot() + 
        geom_line(data=cdf1.r7, aes(x=nc, y=NPP_photo, col="C400")) +   
        geom_line(data=cdf1.r7, aes(x=nc, y=NPP_VL, col="VL")) + 
        geom_line(data=cdf1.r7, aes(x=nc, y=NPP_L, col="L")) +  
        geom_line(data=cdf1.r7, aes(x=nc, y=NPP_M, col="M")) +            
        geom_line(data=cdf2.r7, aes(x=nc, y=NPP_photo, col="C800")) +     
        geom_point(data=edf1.r7, aes(x=nc_VL, y=NPP_VL, fill="A"), 
                color="black", shape=21, size=5) + 
        geom_point(data=edf2.r7, aes(x=nc_L, y=NPP_L, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r7, aes(x=nc_M, y=NPP_M, fill="C"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r7, aes(x=nc_VL, y=NPP_VL, fill="E"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf1.r7, aes(x=nc_VL, y=NPP_I, fill="B"), 
                   shape=21, color="black", size=5) +
        ylim(1.25, 2.25) + 
        xlim(0.001, 0.03) +
        labs(x="Leaf N:C Ratio", 
             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = cbbPalette[7],
                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4], "E" = cbPalette[7])) 
    plot(p)
    if (save.plot) {
        pdf("Plots/Run07.pdf", width=8, height=6)
        plot(p)
        dev.off()
    } 

```

Here we see that priming is not really doing a lot. The treatment for priming effect is that, we have certain C allocated to root exudated and help decomposition of slow SOM, hence more N is recycled back at the M term. 

<br/>
<br/>

### Run 8: Priming effect, variable wood NC ratio
``` {r, echo = FALSE, fig.width = 7, fig.align="center", message = FALSE, warning = FALSE, results="hide"}
    ### split the data.list into two data frames
    cdf.r8 <- r8$cDF
    edf.r8 <- r8$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r8 <- cdf.r8[cdf.r8$CO2 == 400, ]
    cdf2.r8 <- cdf.r8[cdf.r8$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r8 <- edf.r8[edf.r8$CO2 == 400, ]
    edf2.r8 <- edf.r8[edf.r8$CO2 == 800, ]

    ### plot
    p<-ggplot() + 
        geom_line(data=cdf1.r8, aes(x=nc, y=NPP_photo, col="C400")) +   
        geom_line(data=cdf1.r8, aes(x=nc, y=NPP_VL, col="VL")) + 
        geom_line(data=cdf1.r8, aes(x=nc, y=NPP_L, col="L")) +  
        geom_line(data=cdf1.r8, aes(x=nc, y=NPP_M, col="M")) +            
        geom_line(data=cdf2.r8, aes(x=nc, y=NPP_photo, col="C800")) +     
        geom_point(data=edf1.r8, aes(x=nc_VL, y=NPP_VL, fill="A"), 
                color="black", shape=21, size=5) + 
        geom_point(data=edf2.r8, aes(x=nc_L, y=NPP_L, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r8, aes(x=nc_M, y=NPP_M, fill="C"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r8, aes(x=nc_VL, y=NPP_VL, fill="E"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf1.r8, aes(x=nc_VL, y=NPP_I, fill="B"), 
                   shape=21, color="black", size=5) +
        ylim(1.25, 2.25) + 
        xlim(0.001, 0.03) +
        labs(x="Leaf N:C Ratio", 
             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = cbbPalette[7],
                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4], "E" = cbPalette[7])) 

    plot(p)
    if (save.plot) {
        pdf("Plots/Run08.pdf", width=8, height=6)
        plot(p)
        dev.off()
    } 
```

Under variable wood NC assumption, we see that the M, L and VL line overlap with each other. This means that, the slow pool decomposion is hugely stimulated by root exudates, making this as large a pool as the passive N release pool. Previously (under variable wood NC but without priming effect) the combination of slow SOM and wood pools should be very small, because little N is stored in wood pool and the decomposition of slow SOM is not as fast as the decomposition of slow SOM stimualted by root exudates. 

``` {r, echo=FALSE}

    ### split the data.list into two data frames
    cdf.r70 <- r70$cDF
    edf.r70 <- r70$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r70 <- cdf.r70[cdf.r70$CO2 == 400, ]
    cdf2.r70 <- cdf.r70[cdf.r70$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r70 <- edf.r70[edf.r70$CO2 == 400, ]
    edf2.r70 <- edf.r70[edf.r70$CO2 == 800, ]
    
    ### split the data.list into two data frames
    cdf.r80 <- r80$cDF
    edf.r80 <- r80$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r80 <- cdf.r80[cdf.r80$CO2 == 400, ]
    cdf2.r80 <- cdf.r80[cdf.r80$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r80 <- edf.r80[edf.r80$CO2 == 400, ]
    edf2.r80 <- edf.r80[edf.r80$CO2 == 800, ]
    
    
    
    ### compute % change in NC and NPP under eCO2 condition
    chang7 <- (edf2.r7[,2:8] - edf1.r7[,2:8]) / edf1.r7[,2:8] * 100
    chang7$NPP_I <- (edf1.r7[,"NPP_I"] - edf1.r7[,"NPP_VL"]) / edf1.r7[,"NPP_VL"] * 100

    chang8 <- (edf2.r8[,2:8] - edf1.r8[,2:8]) / edf1.r8[,2:8] * 100
    chang8$NPP_I <- (edf1.r8[,"NPP_I"] - edf1.r8[,"NPP_VL"]) / edf1.r8[,"NPP_VL"] * 100
    
    chang70 <- (edf2.r70[,2:8] - edf1.r70[,2:8]) / edf1.r70[,2:8] * 100
    chang70$NPP_I <- (edf1.r70[,"NPP_I"] - edf1.r70[,"NPP_VL"]) / edf1.r70[,"NPP_VL"] * 100

    chang80 <- (edf2.r80[,2:8] - edf1.r80[,2:8]) / edf1.r80[,2:8] * 100
    chang80$NPP_I <- (edf1.r80[,"NPP_I"] - edf1.r80[,"NPP_VL"]) / edf1.r80[,"NPP_VL"] * 100
    
    ### organize the df for output
    changeDF <- round(rbind(chang2, chang1, chang7, chang8, chang70, chang80),1)
    changeDF[1,"Run"] <- "Baseline, fixed NC" 
    changeDF[2,"Run"] <- "Baseline, variable NC" 
    changeDF[3,"Run"] <- "Priming, fixed NC" 
    changeDF[4,"Run"] <- "Priming, variable NC" 
    changeDF[5,"Run"] <- "Exudation, fixed NC" 
    changeDF[6,"Run"] <- "Exudation, variable NC" 
    changeDF <- changeDF[,c("Run", "NPP_I", "NPP_M", "NPP_L", "NPP_VL", 
                            "nc_M", "nc_L", "nc_VL")]
    
    subDF <- changeDF[c(2,4,6),]
    
    write.csv(subDF, "Tables/Exudation_and_priming_effect.csv", row.names=F)
    
    ### table output
    kable(changeDF, row.names=F, format = "markdown")
    
    
    ### Make a bar plot of M-term response
    subDF[1,"Run"] <- "Baseline"
    subDF[2,"Run"] <- "Priming"
    subDF[3,"Run"] <- "Exudation"
    
    subDF[1,"F1"] <- "M"
    subDF[2,"F1"] <- "M"
    subDF[3,"F1"] <- "M"
    
    subDF[1,"F2"] <- "VL"
    subDF[2,"F2"] <- "VL"
    subDF[3,"F2"] <- "VL"

    p<-ggplot(data=subDF) + 
        geom_bar(aes(x=Run, y=NPP_VL, fill=F2), stat="identity") +
        geom_bar(aes(x=Run, y=NPP_M, fill=F1), stat="identity") +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              panel.grid.major=element_line(color="grey")) + 
        ylab("% change") + xlab("") +
        ylim(0, 15)+
        scale_fill_manual(name="Timescale", values = c("VL"="black",
    "M"="grey"))
    plot(p)
    
    if (save.plot) {
        pdf("Plots/exudation_effect.pdf", width=10, height=6)
        plot(p)
        dev.off()
    } 
    

```

Comparing the CO2 effect across models (baseline vs. priming effect), we see that priming effect only stimulated little extra M term NPP under eCO2 as compared to the baseline model, under fixed wood NC assumption. However, the fixed model comparison tells a completely different story. M term response hugely differ, with priming stimulated an almost double M-term NPP response to eCO2 as compared to the baseline model. 

<br/>
<br/>


### Run 9: CLM potential NPP, Fixed wood N:C ratio

CLM has a potential NPP approach, with the following set of equations:
$$P_{dem} = NPP_{pot}n_p $$

where *P<sub>dem</sub>* is the plant N demand, and *n<sub>p</sub>* is the plant N:C ratio.

The immobilization N demand is calculated as (*I<sub>dem</sub>*):
$$ I_{dem} = fC_{lit}t_{lit}(n_s - n_p)$$

where *f* is the fraction of litter C becomes soil C, *C<sub>lit</sub>* is the total litter C, *t<sub>lit</sub>* is the turnover time of the litter pool, and *<sub>ns</sub>* is the soil NC ratio. *C<sub>lit</sub>* can also be expressed as:
$$C_{lit} = NPP_{pot}/t_{lit} $$


Actual plant uptake is expressed as:
$$P_{act} = min(N_{min}P_{dem}/(I_{dem} + P_{dem}), P_{dem}) $$

and actually immobilisation is expressed as:
$$I_{act} = min(N_{min}I_{dem}/(I_{dem} + P_{dem}), I_{dem}) $$

Actual NPP is expressed as:
$$NPP_{act} = NPP_{pot}P_{act}/P_{dem}$$

Now, for VL-term constraint, we have
$$N_{in} = N_{loss} = l_nN_{min}$$

Under potential NPP, we can calculate NPP by re-arranging the equations, as:
$$NPP = N_{in} / [l_n/(1-l_n)n_p] $$

For an actual NPP, we need to consider immobilization demand, expressed as:
$$N_{min} = l_nN_{min} + P_{act} + I_{act} = l_nN_{min} + NPP_{act}n_p + NPP_{act}f(n_s-n_p) $$

Re-arrange, we get:

$$NPP_{act} = N_{in}(1-l_n) / l_n[n_p + f(n_s-n_p)] $$

Similar equations can be derived for L and M term constraints. 

Visionally, with potential NPP and actual NPP, we have the QE framework looks like this:

``` {r, echo = FALSE, fig.width = 7, fig.align="center", message = FALSE, warning = FALSE, results="hide"}
    ### split the data.list into two data frames
    cdf.r9 <- r9$cDF
    edf.r9 <- r9$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r9 <- cdf.r9[cdf.r9$CO2 == 400, ]
    cdf2.r9 <- cdf.r9[cdf.r9$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r9 <- edf.r9[edf.r9$CO2 == 400, ]
    edf2.r9 <- edf.r9[edf.r9$CO2 == 800, ]

    ### plot
    p<-ggplot() + 
        geom_line(data=cdf1.r9, aes(x=nc, y=NPP_photo, col="C400")) +   
        geom_line(data=cdf1.r9, aes(x=nc, y=NPP_VL, col="VL")) + 
        geom_line(data=cdf1.r9, aes(x=nc, y=NPP_L, col="L")) +  
        geom_line(data=cdf1.r9, aes(x=nc, y=NPP_M, col="M")) +            
        geom_line(data=cdf2.r9, aes(x=nc, y=NPP_photo, col="C800")) +     
        geom_point(data=edf1.r9, aes(x=nc_act_VL, y=NPP_I, fill="B"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf1.r9, aes(x=nc_act_VL, y=NPP_act_VL, fill="A"), 
                color="black", shape=21, size=5) + 
        geom_point(data=edf2.r9, aes(x=nc_act_L, y=NPP_act_L, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r9, aes(x=nc_act_M, y=NPP_act_M, fill="C"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r9, aes(x=nc_act_VL, y=NPP_act_VL, fill="E"), 
                   shape=21, color="black", size=5) +
        ylim(0.5, 2.0) + 
        xlim(0.001, 0.01) +
        labs(x="Leaf N:C Ratio", 
             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = cbbPalette[7],
                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4], "E" = cbPalette[7])) 

    plot(p)
    if (save.plot) {
        pdf("Plots/Run09.pdf", width=8, height=6)
        plot(p)
        dev.off()
    } 
```

The unfilled dots are the potential equilibrated NPP, whereas the filled dots are the actually equilibrated NPP at various timesteps. 

Here we can see that, at equilibrium, the actual NPP is not intersecting with the photosynthetic constraint line, and is much lower than the potential NPP, due to the N immobilization demand. We have effectively removed actual NPP dependence on potential NPP. 


### Run 10: CLM potential NPP, variable wood N:C ratio
``` {r, echo = FALSE, fig.width = 7, fig.align="center", message = FALSE, warning = FALSE, results="hide"}
    ### split the data.list into two data frames
    cdf.r10 <- r10$cDF
    edf.r10 <- r10$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r10 <- cdf.r10[cdf.r10$CO2 == 400, ]
    cdf2.r10 <- cdf.r10[cdf.r10$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r10 <- edf.r10[edf.r10$CO2 == 400, ]
    edf2.r10 <- edf.r10[edf.r10$CO2 == 800, ]

    ### plot
    p<-ggplot() + 
        geom_line(data=cdf1.r10, aes(x=nc, y=NPP_photo, col="C400")) +   
        geom_line(data=cdf1.r10, aes(x=nc, y=NPP_VL, col="VL")) + 
        geom_line(data=cdf1.r10, aes(x=nc, y=NPP_L, col="L")) +  
        geom_line(data=cdf1.r10, aes(x=nc, y=NPP_M, col="M")) +            
        geom_line(data=cdf2.r10, aes(x=nc, y=NPP_photo, col="C800")) +     
        geom_point(data=edf1.r10, aes(x=nc_act_VL, y=NPP_I, fill="B"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf1.r10, aes(x=nc_act_VL, y=NPP_act_VL, fill="A"), 
                color="black", shape=21, size=5) + 
        geom_point(data=edf2.r10, aes(x=nc_act_L, y=NPP_act_L, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r10, aes(x=nc_act_M, y=NPP_act_M, fill="C"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r10, aes(x=nc_act_VL, y=NPP_act_VL, fill="E"), 
                   shape=21, color="black", size=5) +
        ylim(0.5, 2.0) + 
        xlim(0.001, 0.01) +
        labs(x="Leaf N:C Ratio", 
             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = cbbPalette[7],
                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4], "E" = cbPalette[7])) 

    plot(p)
    if (save.plot) {
        pdf("Plots/Run10.pdf", width=8, height=6)
        plot(p)
        dev.off()
    } 
```

<br/><br/>

### CLM N uptake 

Based on equations from Ghimire et al. (2016), CLM representation of N uptake is the same as O-CN. So here I can ignore the testing of this function.

<br/><br/>

### Run 11. CLM model C cost of N uptake, fixed wood N:C ratio
This treatment is based on Fisher et al. (2010) and Shi et al. (2015).

There are 4 pathways through which N is absorbed:

    - Diffusion - as part of advection
    - Retranslocation
    - Active uptake
    - BNF

Retranslocation or resorption requires C to synthesize the enzymes and regulatory elements that degrade and remobilize leaf nutrients, and to drive the translocation stream in which the nutrients are suspended. 

Active uptake is an ion-specific enzyme-catalyzed process analogous to Michaelis-Menten kinetics. Energy demand for ion uptake can consume a substantial amount of total respiratory C to move N against concentration gradients. 

BNF is performed by bacteria living in symbiosis within root modules on certain types of plants.

Here, I am calculating C cost of N uptake based on passive uptake, retranslocation, active uptake and fixation. 

``` {r, echo = FALSE, fig.width = 7, fig.align="center", message = FALSE, warning = FALSE, results="hide"}
#    ### split the data.list into two data frames
#    cdf.r11 <- r11$cDF
#    edf.r11 <- r11$eDF
#
#    ### extract CO2 conc. from constraint DF
#    cdf1.r11 <- cdf.r11[cdf.r11$CO2 == 400, ]
#    cdf2.r11 <- cdf.r11[cdf.r11$CO2 == 800, ]
#
#    ### extract CO2 conc. from eq DF
#    edf1.r11 <- edf.r11[edf.r11$CO2 == 400, ]
#    edf2.r11 <- edf.r11[edf.r11$CO2 == 800, ]
#
#    ### plot
#    p<-ggplot() + 
#        geom_line(data=cdf1.r11, aes(x=nc, y=NPP_photo, col="C400")) +   
#        geom_line(data=cdf1.r11, aes(x=nc, y=NPP_VL, col="VL")) + 
#        geom_line(data=cdf1.r11, aes(x=nc, y=NPP_L, col="L")) +  
#        geom_line(data=cdf1.r11, aes(x=nc, y=NPP_M, col="M")) +            
#        geom_line(data=cdf2.r11, aes(x=nc, y=NPP_photo, col="C800")) +     
#        geom_point(data=edf1.r11, aes(x=nc_VL, y=NPP_I, fill="B"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf1.r11, aes(x=nc_VL, y=NPP_VL, fill="A"), 
#                color="black", shape=21, size=5) + 
#        geom_point(data=edf2.r11, aes(x=nc_L, y=NPP_L, fill="D"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf2.r11, aes(x=nc_M, y=NPP_M, fill="C"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf2.r11, aes(x=nc_VL, y=NPP_VL, fill="E"), 
#                   shape=21, color="black", size=5) +
#        ylim(0.0, 2.0) + 
#        xlim(0.001, 0.01) +
#        labs(x="Leaf N:C Ratio", 
#             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
#        theme_linedraw() +
#        theme(panel.grid.minor=element_blank(),
#              axis.text=element_text(size=14),
#              axis.title=element_text(size=16),
#              legend.text=element_text(size=14),
#              legend.title=element_text(size=16),
#              panel.grid.major=element_line(color="grey")) +
#              scale_colour_manual(name="Constraint line", 
#                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = #cbbPalette[7],
#                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
#              scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = #cbPalette[2], "C" = cbPalette[3],
#                                                            "D" = cbPalette[4], "E" = #cbPalette[7])) 
#
#    plot(p)
#    if (save.plot) {
#        pdf("Plots/Run11.pdf", width=8, height=6)
#        plot(p)
#        dev.off()
#    } 
```


<br/><br/>

### Run 12. CLM model C cost of N uptake, variable wood N:C ratio
``` {r, echo = FALSE, fig.width = 7, fig.align="center", message = FALSE, warning = FALSE, results="hide"}
#    ### split the data.list into two data frames
#    cdf.r12 <- r12$cDF
#    edf.r12 <- r12$eDF
#
#    ### extract CO2 conc. from constraint DF
#    cdf1.r12 <- cdf.r12[cdf.r12$CO2 == 400, ]
#    cdf2.r12 <- cdf.r12[cdf.r12$CO2 == 800, ]
#
#    ### extract CO2 conc. from eq DF
#    edf1.r12 <- edf.r12[edf.r12$CO2 == 400, ]
#    edf2.r12 <- edf.r12[edf.r12$CO2 == 800, ]
#
#    ### plot
#    p<-ggplot() + 
#        geom_line(data=cdf1.r12, aes(x=nc, y=NPP_photo, col="C400")) +   
#        geom_line(data=cdf1.r12, aes(x=nc, y=NPP_VL, col="VL")) + 
#        geom_line(data=cdf1.r12, aes(x=nc, y=NPP_L, col="L")) +  
#        geom_line(data=cdf1.r12, aes(x=nc, y=NPP_M, col="M")) +            
#        geom_line(data=cdf2.r12, aes(x=nc, y=NPP_photo, col="C800")) +     
#        geom_point(data=edf1.r12, aes(x=nc_VL, y=NPP_I, fill="B"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf1.r12, aes(x=nc_VL, y=NPP_VL, fill="A"), 
#                color="black", shape=21, size=5) + 
#        geom_point(data=edf2.r12, aes(x=nc_L, y=NPP_L, fill="D"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf2.r12, aes(x=nc_M, y=NPP_M, fill="C"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf2.r12, aes(x=nc_VL, y=NPP_VL, fill="E"), 
#                   shape=21, color="black", size=5) +
#        ylim(0.0, 2.0) + 
#        xlim(0.001, 0.01) +
#        labs(x="Leaf N:C Ratio", 
#             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
#        theme_linedraw() +
#        theme(panel.grid.minor=element_blank(),
#              axis.text=element_text(size=14),
#              axis.title=element_text(size=16),
#              legend.text=element_text(size=14),
#              legend.title=element_text(size=16),
#              panel.grid.major=element_line(color="grey")) +
#              scale_colour_manual(name="Constraint line", 
#                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = #cbbPalette[7],
#                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
#              scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = #cbPalette[2], "C" = cbPalette[3],
#                                                            "D" = cbPalette[4], "E" = #cbPalette[7])) 
#
#    plot(p)
#    if (save.plot) {
#        pdf("Plots/Run12.pdf", width=8, height=6)
#        plot(p)
#        dev.off()
#    } 
```

# 4. Summary table:
``` {r, echo=FALSE}
#    source("QE_Analyses/make_summary_table.R")
#    sumDF <- make_summary_table(destDir = "Tables")
#    
#    ### table output
#    kable(sumDF, row.names=F, format = "markdown")

```

The table above summarizes scenarios tested in this study. Next, I am going to make a barplot to visually compare the results. 

``` {r, echo=FALSE, warning=F, message=F}
#    ### Firstly, remove unwanted scenarios
#    sumDF2 <- sumDF[grepl("variable wood", sumDF$model),]
#    sumDF3 <-  rbind(sumDF2, sumDF[sumDF$model == "Baseline fixed wood", ])     
#    sumDF3$model <- gsub(", variable wood", "", sumDF3$model)
#    
#    ### get % change df and abs df
#    pDF <- data.frame(sumDF3$model, sumDF3$I, sumDF3$M, sumDF3$L, sumDF3$VL)
#    colnames(pDF) <- c("Model", "I", "M", "L", "VL")
#    pcDF <- melt(pDF, id="Model")
#    nDF <- sumDF3[,1:3]
#    nppDF <- melt(nDF, id="model")
#    
#    ### plot
#    p<-ggplot() + 
#        geom_bar(data=pcDF, aes(x=Model, y=value, fill=variable), 
#                 stat="identity", position=position_dodge()) +
#        theme_linedraw() +
#        theme(panel.grid.minor=element_blank(),
#              axis.text=element_text(size=14),
#              axis.title=element_text(size=16),
#              legend.text=element_text(size=14),
#              legend.title=element_text(size=16),
#              axis.text.x = element_text(angle = 90, hjust=1),
#              panel.grid.major=element_line(color="grey")) + 
#        ylab("% change") + xlab("") + 
#        scale_y_log10(breaks=c(1,3,10,30,100)) + 
#        scale_fill_manual(name="Timescale", values = c("I" = cbPalette[2], "M" = cbPalette[3],
#                                                      "L" = cbPalette[4], "VL" = cbPalette[7])) #+
#        coord_flip()
#    plot(p)
#    if (save.plot) {
#        pdf("Plots/summary_plot.pdf", width=10, height=6)
#        plot(p)
#        dev.off()
#    } 
#    
#
```


# 5. Sensitivity analyses:

### 5.1 N uptake coefficient
``` {r, echo=FALSE, warning=F, message=F}

    sensDF <- matrix(nrow=5*4, ncol=3)
    colnames(sensDF) <- c("u", "value", "Time")
    sensDF <- as.data.frame(sensDF)
    sensDF$u <- rep(c(0.2, 0.5, 1, 2, 5), each = 4)
    sensDF$Time <- rep(c("I", "M", "L", "VL"), by = 5)

    sensDF[sensDF$u == 0.2 & sensDF$Time == "I", "value"] <- (r32$eDF$NPP_I[2] - r32$eDF$NPP_VL[1])/r32$eDF$NPP_VL[1]
    sensDF[sensDF$u == 0.2 & sensDF$Time == "M", "value"] <- (r32$eDF$NPP_M[2] - r32$eDF$NPP_VL[1])/r32$eDF$NPP_VL[1]
    sensDF[sensDF$u == 0.2 & sensDF$Time == "L", "value"] <- (r32$eDF$NPP_L[2] - r32$eDF$NPP_VL[1])/r32$eDF$NPP_VL[1]
    sensDF[sensDF$u == 0.2 & sensDF$Time == "VL", "value"] <- (r32$eDF$NPP_VL[2] - r32$eDF$NPP_VL[1])/r32$eDF$NPP_VL[1]

    sensDF[sensDF$u == 0.5 & sensDF$Time == "I", "value"] <- (r31$eDF$NPP_I[2] - r31$eDF$NPP_VL[1])/r31$eDF$NPP_VL[1]
    sensDF[sensDF$u == 0.5 & sensDF$Time == "M", "value"] <- (r31$eDF$NPP_M[2] - r31$eDF$NPP_VL[1])/r31$eDF$NPP_VL[1]
    sensDF[sensDF$u == 0.5 & sensDF$Time == "L", "value"] <- (r31$eDF$NPP_L[2] - r31$eDF$NPP_VL[1])/r31$eDF$NPP_VL[1]
    sensDF[sensDF$u == 0.5 & sensDF$Time == "VL", "value"] <- (r31$eDF$NPP_VL[2] - r31$eDF$NPP_VL[1])/r31$eDF$NPP_VL[1]
    
    sensDF[sensDF$u == 1.0 & sensDF$Time == "I", "value"] <- (r3$eDF$NPP_I[2] - r3$eDF$NPP_VL[1])/r3$eDF$NPP_VL[1]
    sensDF[sensDF$u == 1.0 & sensDF$Time == "M", "value"] <- (r3$eDF$NPP_M[2] - r3$eDF$NPP_VL[1])/r3$eDF$NPP_VL[1]
    sensDF[sensDF$u == 1.0 & sensDF$Time == "L", "value"] <- (r3$eDF$NPP_L[2] - r3$eDF$NPP_VL[1])/r3$eDF$NPP_VL[1]
    sensDF[sensDF$u == 1.0 & sensDF$Time == "VL", "value"] <- (r3$eDF$NPP_VL[2] - r3$eDF$NPP_VL[1])/r3$eDF$NPP_VL[1]
    
    sensDF[sensDF$u == 2 & sensDF$Time == "I", "value"] <- (r33$eDF$NPP_I[2] - r33$eDF$NPP_VL[1])/r33$eDF$NPP_VL[1]
    sensDF[sensDF$u == 2 & sensDF$Time == "M", "value"] <- (r33$eDF$NPP_M[2] - r33$eDF$NPP_VL[1])/r33$eDF$NPP_VL[1]
    sensDF[sensDF$u == 2 & sensDF$Time == "L", "value"] <- (r33$eDF$NPP_L[2] - r33$eDF$NPP_VL[1])/r33$eDF$NPP_VL[1]
    sensDF[sensDF$u == 2 & sensDF$Time == "VL", "value"] <- (r33$eDF$NPP_VL[2] - r33$eDF$NPP_VL[1])/r33$eDF$NPP_VL[1]
    
    sensDF[sensDF$u == 5 & sensDF$Time == "I", "value"] <- (r34$eDF$NPP_I[2] - r34$eDF$NPP_VL[1])/r34$eDF$NPP_VL[1]
    sensDF[sensDF$u == 5 & sensDF$Time == "M", "value"] <- (r34$eDF$NPP_M[2] - r34$eDF$NPP_VL[1])/r34$eDF$NPP_VL[1]
    sensDF[sensDF$u == 5 & sensDF$Time == "L", "value"] <- (r34$eDF$NPP_L[2] - r34$eDF$NPP_VL[1])/r34$eDF$NPP_VL[1]
    sensDF[sensDF$u == 5 & sensDF$Time == "VL", "value"] <- (r34$eDF$NPP_VL[2] - r34$eDF$NPP_VL[1])/r34$eDF$NPP_VL[1]
    
    sensDF$u <- as.character(sensDF$u)
    sensDF$value <- sensDF$value * 100
    
    ### plot
    p<-ggplot() + 
        geom_bar(data=sensDF, aes(x=u, y=value, fill=Time), 
                 stat="identity", position=position_dodge()) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              axis.text.x = element_text(angle = 90, hjust=1),
              panel.grid.major=element_line(color="grey")) + 
        ylab("% change") + xlab("") + 
        scale_fill_manual(name="Timescale", values = c("I" = cbPalette[2], "M" = cbPalette[3],
                                                      "L" = cbPalette[4], "VL" = cbPalette[7])) 
    plot(p)
    if (save.plot) {
        pdf("Plots/sensitivity_n_uptake_coefficient.pdf", width=10, height=6)
        plot(p)
        dev.off()
    } 
    
    
    ### save a table on equilibrated NPP
    u1 <- r32$eDF
    u2 <- r31$eDF
    u3 <- r3$eDF
    u4 <- r33$eDF
    u5 <- r34$eDF
    u1$u <- 0.2
    u2$u <- 0.5
    u3$u <- 1.0
    u4$u <- 2.0
    u5$u <- 5.0
    outDF <- rbind(u1, u2, u3, u4, u5)
    write.csv(outDF, "Tables/varying_u_summary_table.csv", row.names=F)
    

```

### 5.2 Varying ar and sr for GDAY
``` {r, echo=FALSE, warning=F, message=F}
    sensDF <- matrix(nrow=5*4, ncol=4)
    colnames(sensDF) <- c("ar", "sr", "value", "Time")
    sensDF <- as.data.frame(sensDF)
    sensDF$ar <- rep(c(0.2, 0.4, 0.6, 0.2, 0.2), each = 4)
    sensDF$sr <- rep(c(1.5, 1.5, 1.5, 5, 10), each = 4)

    sensDF$Time <- rep(c("I", "M", "L", "VL"), by = 5)

    sensDF[sensDF$ar == 0.2 & sensDF$sr == 1.5 & sensDF$Time == "I", "value"] <- (r5$eDF$NPP_I[2] - r5$eDF$NPP_VL[1])/r5$eDF$NPP_VL[1]
    sensDF[sensDF$ar == 0.2 & sensDF$sr == 1.5 & sensDF$Time == "M", "value"] <- (r5$eDF$NPP_M[2] - r5$eDF$NPP_VL[1])/r5$eDF$NPP_VL[1]
    sensDF[sensDF$ar == 0.2 & sensDF$sr == 1.5 & sensDF$Time == "L", "value"] <- (r5$eDF$NPP_L[2] - r5$eDF$NPP_VL[1])/r5$eDF$NPP_VL[1]
    sensDF[sensDF$ar == 0.2 & sensDF$sr == 1.5 & sensDF$Time == "VL", "value"] <- (r5$eDF$NPP_VL[2] - r5$eDF$NPP_VL[1])/r5$eDF$NPP_VL[1]

    sensDF[sensDF$ar == 0.4 & sensDF$sr == 1.5 & sensDF$Time == "I", "value"] <- (r51$eDF$NPP_I[2] - r51$eDF$NPP_VL[1])/r51$eDF$NPP_VL[1]
    sensDF[sensDF$ar == 0.4 & sensDF$sr == 1.5 & sensDF$Time == "M", "value"] <- (r51$eDF$NPP_M[2] - r51$eDF$NPP_VL[1])/r51$eDF$NPP_VL[1]
    sensDF[sensDF$ar == 0.4 & sensDF$sr == 1.5 & sensDF$Time == "L", "value"] <- (r51$eDF$NPP_L[2] - r51$eDF$NPP_VL[1])/r51$eDF$NPP_VL[1]
    sensDF[sensDF$ar == 0.4 & sensDF$sr == 1.5 & sensDF$Time == "VL", "value"] <- (r51$eDF$NPP_VL[2] - r51$eDF$NPP_VL[1])/r51$eDF$NPP_VL[1]
    
    sensDF[sensDF$ar == 0.6 & sensDF$sr == 1.5 & sensDF$Time == "I", "value"] <- (r52$eDF$NPP_I[2] - r52$eDF$NPP_VL[1])/r52$eDF$NPP_VL[1]
    sensDF[sensDF$ar == 0.6 & sensDF$sr == 1.5 & sensDF$Time == "M", "value"] <- (r52$eDF$NPP_M[2] - r52$eDF$NPP_VL[1])/r52$eDF$NPP_VL[1]
    sensDF[sensDF$ar == 0.6 & sensDF$sr == 1.5 & sensDF$Time == "L", "value"] <- (r52$eDF$NPP_L[2] - r52$eDF$NPP_VL[1])/r52$eDF$NPP_VL[1]
    sensDF[sensDF$ar == 0.6 & sensDF$sr == 1.5 & sensDF$Time == "VL", "value"] <- (r52$eDF$NPP_VL[2] - r52$eDF$NPP_VL[1])/r52$eDF$NPP_VL[1]
    
    sensDF[sensDF$ar == 0.2 & sensDF$sr == 5 & sensDF$Time == "I", "value"] <- (r53$eDF$NPP_I[2] - r53$eDF$NPP_VL[1])/r53$eDF$NPP_VL[1]
    sensDF[sensDF$ar == 0.2 & sensDF$sr == 5 & sensDF$Time == "M", "value"] <- (r53$eDF$NPP_M[2] - r53$eDF$NPP_VL[1])/r53$eDF$NPP_VL[1]
    sensDF[sensDF$ar == 0.2 & sensDF$sr == 5 & sensDF$Time == "L", "value"] <- (r53$eDF$NPP_L[2] - r53$eDF$NPP_VL[1])/r53$eDF$NPP_VL[1]
    sensDF[sensDF$ar == 0.2 & sensDF$sr == 5 & sensDF$Time == "VL", "value"] <- (r53$eDF$NPP_VL[2] - r53$eDF$NPP_VL[1])/r53$eDF$NPP_VL[1]
    
    sensDF[sensDF$ar == 0.2 & sensDF$sr == 10 & sensDF$Time == "I", "value"] <- (r54$eDF$NPP_I[2] - r54$eDF$NPP_VL[1])/r54$eDF$NPP_VL[1]
    sensDF[sensDF$ar == 0.2 & sensDF$sr == 10 & sensDF$Time == "M", "value"] <- (r54$eDF$NPP_M[2] - r54$eDF$NPP_VL[1])/r54$eDF$NPP_VL[1]
    sensDF[sensDF$ar == 0.2 & sensDF$sr == 10 & sensDF$Time == "L", "value"] <- (r54$eDF$NPP_L[2] - r54$eDF$NPP_VL[1])/r54$eDF$NPP_VL[1]
    sensDF[sensDF$ar == 0.2 & sensDF$sr == 10 & sensDF$Time == "VL", "value"] <- (r54$eDF$NPP_VL[2] - r54$eDF$NPP_VL[1])/r54$eDF$NPP_VL[1]
    
    sensDF$ar <- as.character(sensDF$ar)
    sensDF$sr <- as.character(sensDF$sr)
    sensDF$value <- sensDF$value * 100
    
    plotDF1 <- subset(sensDF, sr == "1.5")
    plotDF2 <- subset(sensDF, ar == "0.2")

    ### plot
    p1<-ggplot() + 
        geom_bar(data=plotDF1, aes(x=ar, y=value, fill=Time), 
                 stat="identity", position=position_dodge()) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              axis.text.x = element_text(angle = 90, hjust=1),
              panel.grid.major=element_line(color="grey")) + 
        ylab("% change") + xlab("") + 
        scale_fill_manual(name="Timescale", values = c("I" = cbPalette[2], "M" = cbPalette[3],
                                                      "L" = cbPalette[4], "VL" = cbPalette[7])) 
    
    p2<-ggplot() + 
        geom_bar(data=plotDF2, aes(x=sr, y=value, fill=Time), 
                 stat="identity", position=position_dodge()) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              axis.text.x = element_text(angle = 90, hjust=1),
              panel.grid.major=element_line(color="grey")) + 
        ylab("% change") + xlab("") + 
        scale_fill_manual(name="Timescale", values = c("I" = cbPalette[2], "M" = cbPalette[3],
                                                      "L" = cbPalette[4], "VL" = cbPalette[7])) 
    plot(p1)
    plot(p2)
    if (save.plot) {
        pdf("Plots/sensitivity_ar_GDAY.pdf", width=10, height=6)
        plot(p1)
        dev.off()
        
        pdf("Plots/sensitivity_sr_GDAY.pdf", width=10, height=6)
        plot(p2)
        dev.off()
    } 
    
    
    ### save a table on equilibrated NPP
    u1 <- r32$eDF
    u2 <- r31$eDF
    u3 <- r3$eDF
    u4 <- r33$eDF
    u5 <- r34$eDF
    u1$u <- 0.2
    u2$u <- 0.5
    u3$u <- 1.0
    u4$u <- 2.0
    u5$u <- 5.0
    outDF <- rbind(u1, u2, u3, u4, u5)
    write.csv(outDF, "Tables/varying_u_summary_table.csv", row.names=F)
    

```

### 5.3 Priming effect coefficients
``` {r, echo=FALSE, warning=F, message=F}

    sensDF <- matrix(nrow=9, ncol=3)
    colnames(sensDF) <- c("variable", "variable_value", "M")
    sensDF <- as.data.frame(sensDF)
    sensDF$variable <- c("a1", "a1", "a1", "a0", "a0", "a0",
                         "fcue", "fcue", "fcue")
    sensDF$variable_value <- c(0.2, 0.5, 1, 0.01, 0.05, 0.1, 0.3, 0.5, 0.7)
    
    sensDF[sensDF$variable == "a1" & sensDF$variable_value == 0.2, "M"] <- (r81$eDF$NPP_M[2] - r81$eDF$NPP_VL[1])/r81$eDF$NPP_VL[1]
    sensDF[sensDF$variable == "a1" & sensDF$variable_value == 0.5, "M"] <- (r82$eDF$NPP_M[2] - r82$eDF$NPP_VL[1])/r82$eDF$NPP_VL[1]
    sensDF[sensDF$variable == "a1" & sensDF$variable_value == 1.0, "M"] <- (r8$eDF$NPP_M[2] - r8$eDF$NPP_VL[1])/r8$eDF$NPP_VL[1]

    sensDF[sensDF$variable == "a0" & sensDF$variable_value == 0.01, "M"] <- (r8$eDF$NPP_M[2] - r8$eDF$NPP_VL[1])/r8$eDF$NPP_VL[1]
    sensDF[sensDF$variable == "a0" & sensDF$variable_value == 0.05, "M"] <- (r83$eDF$NPP_M[2] - r83$eDF$NPP_VL[1])/r83$eDF$NPP_VL[1]
    sensDF[sensDF$variable == "a0" & sensDF$variable_value == 0.1, "M"] <- (r84$eDF$NPP_M[2] - r84$eDF$NPP_VL[1])/r4$eDF$NPP_VL[1]
    
    sensDF[sensDF$variable == "fcue" & sensDF$variable_value == 0.3, "M"] <- (r8$eDF$NPP_M[2] - r8$eDF$NPP_VL[1])/r8$eDF$NPP_VL[1]
    sensDF[sensDF$variable == "fcue" & sensDF$variable_value == 0.5, "M"] <- (r85$eDF$NPP_M[2] - r85$eDF$NPP_VL[1])/r85$eDF$NPP_VL[1]
    sensDF[sensDF$variable == "fcue" & sensDF$variable_value == 0.7, "M"] <- (r86$eDF$NPP_M[2] - r86$eDF$NPP_VL[1])/r86$eDF$NPP_VL[1]
    
    
    sensDF$variable_value <- as.character(sensDF$variable_value)
    sensDF$M <- sensDF$M * 100
    
    ### plot
    p<-ggplot() + 
        geom_bar(data=sensDF, aes(x=variable, y=M, fill=variable_value), 
                 stat="identity", position=position_dodge()) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              axis.text.x = element_text(angle = 90, hjust=1),
              panel.grid.major=element_line(color="grey")) + 
        ylab("% change") + xlab("") + 
        scale_fill_manual(name="coefficient values", values = c(cbPalette)) 
    
    plot(p)
    if (save.plot) {
        pdf("Plots/sensitivity_priming_effect_coefficient.pdf", width=10, height=6)
        plot(p)
        dev.off()
    } 
    
    sensDF <- read.csv("Tables/Priming_sensitivity.csv")
    
    p<-ggplot() + 
        geom_bar(data=sensDF, aes(x=variable, y=newks, fill=as.factor(variable_value)), 
                 stat="identity", position=position_dodge()) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              axis.text.x = element_text(angle = 90, hjust=1),
              panel.grid.major=element_line(color="grey"),
              legend.position="none") + 
        ylab(expression(k[slow])) + xlab("") + 
        scale_fill_manual(name="coefficient values", 
                          values = c(cbPalette[1:5],cbPalette[1:5],cbPalette[1:5]))+
        geom_hline(yintercept=0.067)+
        scale_x_discrete(labels=c("arhizo" = expression(a[rhizo]), "fcue" = expression(f[cue]), "km" = expression(k[m])))
    
    plot(p)
    if (save.plot) {
        pdf("Plots/sensitivity_priming_effect_coefficient2.pdf", width=10, height=6)
        plot(p)
        dev.off()
    } 

```

### 5.4 Priming effect on NUE

Here NUE is defined as: 
$$NUE = 1 / (n_fa_f + n_wa_w + n_ra_r) $$

#### 5.4.1 Fixed wood NC assumption

``` {r, echo=FALSE, warning=F, message=F}

    af <- 0.2
    ar <- 0.2
    aw <- 0.6
    nwood <- 0.005
    
    ### split the data.list into two data frames
    cdf.r7 <- r7$cDF
    edf.r7 <- r7$eDF
    
    cdf.r7$NUE <- 1 / (cdf.r7$nc * af + cdf.r7$nc * 0.7 * ar + nwood * aw)

    edf.r7$NUE_VL <- 1 / (edf.r7$nc_VL * af + edf.r7$nc_VL * 0.7 * ar + nwood * aw)
    edf.r7$NUE_L <- 1 / (edf.r7$nc_L * af + edf.r7$nc_L * 0.7 * ar + nwood * aw)
    edf.r7$NUE_M <- 1 / (edf.r7$nc_M* af + edf.r7$nc_M * 0.7 * ar + nwood * aw) 

    ### extract CO2 conc. from constraint DF
    cdf1.r7 <- cdf.r7[cdf.r7$CO2 == 400, ]
    cdf2.r7 <- cdf.r7[cdf.r7$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r7 <- edf.r7[edf.r7$CO2 == 400, ]
    edf2.r7 <- edf.r7[edf.r7$CO2 == 800, ]
    
    ### plot
    p<-ggplot() + 
        geom_line(data=cdf1.r7, aes(x=nc, y=NUE, col="NUE trajectory")) +   
        geom_point(data=edf1.r7, aes(x=nc_VL, y=NUE_VL, fill="A"), 
                color="black", shape=21, size=5) + 
        geom_point(data=edf2.r7, aes(x=nc_L, y=NUE_L, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r7, aes(x=nc_M, y=NUE_M, fill="C"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r7, aes(x=nc_VL, y=NUE_VL, fill="E"), 
                   shape=21, color="black", size=5) +
        ylim(100, 250) + 
        xlim(0.001, 0.03) +
        labs(x="Leaf N:C Ratio", 
             y="NUE") +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=16),
              axis.title=element_text(size=18),
              legend.text=element_text(size=16),
              legend.title=element_text(size=18),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("NUE trajectory" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = cbbPalette[7],
                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4], "E" = cbPalette[7])) 

    plot(p)
    if (save.plot) {
        pdf("Plots/sensitivity_priming_NUE_fixed_wood.pdf", width=10, height=6)
        plot(p)
        dev.off()
    } 
    

```

#### 5.4.2 Fixed wood NC assumption

``` {r, echo=FALSE, warning=F, message=F}

    af <- 0.2
    ar <- 0.2
    aw <- 0.6
    nwood <- 0.005
    
    ### split the data.list into two data frames
    cdf.r8 <- r8$cDF
    edf.r8 <- r8$eDF
    
    cdf.r8$NUE <- 1 / (cdf.r8$nc * af + cdf.r8$nc * 0.7 * ar + cdf.r8$nc * nwood * aw)
    
    edf.r8$NUE_VL <- 1 / (edf.r8$nc_VL * af + edf.r8$nc_VL * 0.7 * ar + edf.r8$nc_M * nwood * aw)
    edf.r8$NUE_L <- 1 / (edf.r8$nc_L * af + edf.r8$nc_L * 0.7 * ar + edf.r8$nc_M * nwood * aw)
    edf.r8$NUE_M <- 1 / (edf.r8$nc_M* af + edf.r8$nc_M * 0.7 * ar + edf.r8$nc_M * nwood * aw) 

    ### extract CO2 conc. from constraint DF
    cdf1.r8 <- cdf.r8[cdf.r8$CO2 == 400, ]
    cdf2.r8 <- cdf.r8[cdf.r8$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r8 <- edf.r8[edf.r8$CO2 == 400, ]
    edf2.r8 <- edf.r8[edf.r8$CO2 == 800, ]
    
    ### plot
    p<-ggplot() + 
        geom_line(data=cdf1.r8, aes(x=nc, y=NUE, col="NUE trajectory")) +   
        geom_point(data=edf1.r8, aes(x=nc_VL, y=NUE_VL, fill="A"), 
                color="black", shape=21, size=5) + 
        geom_point(data=edf2.r8, aes(x=nc_L, y=NUE_L, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r8, aes(x=nc_M, y=NUE_M, fill="C"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r8, aes(x=nc_VL, y=NUE_VL, fill="E"), 
                   shape=21, color="black", size=5) +
        ylim(100, 250) + 
        xlim(0.001, 0.03) +
        labs(x="Leaf N:C Ratio", 
             y="NUE") +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=16),
              axis.title=element_text(size=18),
              legend.text=element_text(size=16),
              legend.title=element_text(size=18),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("NUE trajectory" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = cbbPalette[7],
                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4], "E" = cbPalette[7])) 

    plot(p)
    if (save.plot) {
        pdf("Plots/sensitivity_priming_NUE_variable_wood.pdf", width=10, height=6)
        plot(p)
        dev.off()
    } 
    

```

### 5.5 NUE of baseline model

#### 5.5.1 Variable wood NC assumption

``` {r, echo=FALSE, warning=F, message=F}

    af <- 0.2
    ar <- 0.2
    aw <- 0.6
    nwood <- 0.005
    
    ### split the data.list into two data frames
    cdf.r1 <- r1$cDF
    edf.r1 <- r1$eDF
    
    cdf.r1$NUE <- 1 / (cdf.r1$nc * af + cdf.r1$nc * 0.7 * ar + nwood * aw)

    edf.r1$NUE_VL <- 1 / (edf.r1$nc_VL * af + edf.r1$nc_VL * 0.7 * ar + nwood * aw)
    edf.r1$NUE_L <- 1 / (edf.r1$nc_L * af + edf.r1$nc_L * 0.7 * ar + nwood * aw)
    edf.r1$NUE_M <- 1 / (edf.r1$nc_M* af + edf.r1$nc_M * 0.7 * ar + nwood * aw) 

    ### extract CO2 conc. from constraint DF
    cdf1.r1 <- cdf.r1[cdf.r1$CO2 == 400, ]
    cdf2.r1 <- cdf.r1[cdf.r1$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r1 <- edf.r1[edf.r1$CO2 == 400, ]
    edf2.r1 <- edf.r1[edf.r1$CO2 == 800, ]
    
    ### plot
    p<-ggplot() + 
        geom_line(data=cdf1.r1, aes(x=nc, y=NUE, col="NUE trajectory")) +   
        geom_point(data=edf1.r1, aes(x=nc_VL, y=NUE_VL, fill="A"), 
                color="black", shape=21, size=5) + 
        geom_point(data=edf2.r1, aes(x=nc_L, y=NUE_L, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r1, aes(x=nc_M, y=NUE_M, fill="C"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r1, aes(x=nc_VL, y=NUE_VL, fill="E"), 
                   shape=21, color="black", size=5) +
        ylim(100, 250) + 
        xlim(0.001, 0.03) +
        labs(x="Leaf N:C Ratio", 
             y="NUE") +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=16),
              axis.title=element_text(size=18),
              legend.text=element_text(size=16),
              legend.title=element_text(size=18),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("NUE trajectory" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = cbbPalette[7],
                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4], "E" = cbPalette[7])) 

    plot(p)
    if (save.plot) {
        pdf("Plots/sensitivity_baseline_NUE_variable_wood.pdf", width=10, height=6)
        plot(p)
        dev.off()
    } 
    

```

#### 5.5.2 Fixed wood NC assumption

``` {r, echo=FALSE, warning=F, message=F, results="hide"}

    af <- 0.2
    ar <- 0.2
    aw <- 0.6
    nwood <- 0.005
    
    ### split the data.list into two data frames
    cdf.r2 <- r2$cDF
    edf.r2 <- r2$eDF
    
    cdf.r2$NUE <- 1 / (cdf.r2$nc * af + cdf.r2$nc * 0.7 * ar + cdf.r2$nc * nwood * aw)
    
    edf.r2$NUE_VL <- 1 / (edf.r2$nc_VL * af + edf.r2$nc_VL * 0.7 * ar + edf.r2$nc_M * nwood * aw)
    edf.r2$NUE_L <- 1 / (edf.r2$nc_L * af + edf.r2$nc_L * 0.7 * ar + edf.r2$nc_M * nwood * aw)
    edf.r2$NUE_M <- 1 / (edf.r2$nc_M* af + edf.r2$nc_M * 0.7 * ar + edf.r2$nc_M * nwood * aw) 

    ### extract CO2 conc. from constraint DF
    cdf1.r2 <- cdf.r2[cdf.r2$CO2 == 400, ]
    cdf2.r2 <- cdf.r2[cdf.r2$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r2 <- edf.r2[edf.r2$CO2 == 400, ]
    edf2.r2 <- edf.r2[edf.r2$CO2 == 800, ]
    
    ### plot
    p<-ggplot() + 
        geom_line(data=cdf1.r2, aes(x=nc, y=NUE, col="NUE trajectory")) +   
        geom_point(data=edf1.r2, aes(x=nc_VL, y=NUE_VL, fill="A"), 
                color="black", shape=21, size=5) + 
        geom_point(data=edf2.r2, aes(x=nc_L, y=NUE_L, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r2, aes(x=nc_M, y=NUE_M, fill="C"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf2.r2, aes(x=nc_VL, y=NUE_VL, fill="E"), 
                   shape=21, color="black", size=5) +
        ylim(100, 800) + 
        xlim(0.001, 0.03) +
        labs(x="Leaf N:C Ratio", 
             y="NUE") +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=16),
              axis.title=element_text(size=18),
              legend.text=element_text(size=16),
              legend.title=element_text(size=18),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("NUE trajectory" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = cbbPalette[7],
                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4], "E" = cbPalette[7])) 

    plot(p)
    if (save.plot) {
        pdf("Plots/sensitivity_baseline_NUE_variable_wood.pdf", width=10, height=6)
        plot(p)
        dev.off()
    } 
    

```


### 5.6. Compare CLM with CLM-FUN
#### 5.6.1 variable wood NC
``` {r, echo=FALSE, warning=F, message=F, results="hide"}
#
#    #### Run 9
#    ### split the data.list into two data frames
#    cdf.r9 <- r9$cDF
#    edf.r9 <- r9$eDF
#
#    ### extract CO2 conc. from constraint DF
#    cdf1.r9 <- cdf.r9[cdf.r9$CO2 == 400, ]
#    cdf2.r9 <- cdf.r9[cdf.r9$CO2 == 800, ]
#
#    ### extract CO2 conc. from eq DF
#    edf1.r9 <- edf.r9[edf.r9$CO2 == 400, ]
#    edf2.r9 <- edf.r9[edf.r9$CO2 == 800, ]
#    
#    #### Run 11
#    ### split the data.list into two data frames
#    cdf.r11 <- r11$cDF
#    edf.r11 <- r11$eDF
#
#    ### extract CO2 conc. from constraint DF
#    cdf1.r11 <- cdf.r11[cdf.r11$CO2 == 400, ]
#    cdf2.r11 <- cdf.r11[cdf.r11$CO2 == 800, ]
#
#    ### extract CO2 conc. from eq DF
#    edf1.r11 <- edf.r11[edf.r11$CO2 == 400, ]
#    edf2.r11 <- edf.r11[edf.r11$CO2 == 800, ]
#    
#    ### plot
#    p1<-ggplot() + 
#        geom_line(data=cdf1.r11, aes(x=nc, y=NPP_photo, col="C400")) +   
#        #geom_line(data=cdf1.r9, aes(x=nc, y=NPP_photo, col="C400"), linetype="dotted") +   
#        geom_line(data=cdf1.r11, aes(x=nc, y=NPP_VL, col="VL")) + 
#        geom_line(data=cdf1.r9, aes(x=nc, y=NPP_VL, col="VL"), linetype="dotted") + 
#        #geom_line(data=cdf1.r11, aes(x=nc, y=NPP_L, col="L")) +  
#        #geom_line(data=cdf1.r9, aes(x=nc, y=NPP_L, col="L"), linetype="dotted") +  
#        #geom_line(data=cdf1.r11, aes(x=nc, y=NPP_M, col="M")) +    
#        #geom_line(data=cdf1.r9, aes(x=nc, y=NPP_M, col="M"), linetype="dotted") +            
#        geom_line(data=cdf2.r11, aes(x=nc, y=NPP_photo, col="C800")) +     
#        #geom_line(data=cdf2.r9, aes(x=nc, y=NPP_photo, col="C800"), linetype="dotted") +     
#        geom_point(data=edf1.r11, aes(x=nc_VL, y=NPP_I, fill="B"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf1.r9, aes(x=nc_act_VL, y=NPP_I), 
#                   shape=19, color=cbPalette[2], size=5) +
#        geom_point(data=edf1.r11, aes(x=nc_VL, y=NPP_VL, fill="A"), 
#                color="black", shape=21, size=5) + 
#        geom_point(data=edf1.r9, aes(x=nc_act_VL, y=NPP_act_VL), 
#                color=cbPalette[1], shape=19, size=5) + 
#        #geom_point(data=edf2.r11, aes(x=nc_L, y=NPP_L, fill="D"), 
#        #           shape=21, color="black", size=5) +
#        #geom_point(data=edf2.r9, aes(x=nc_act_L, y=NPP_act_L), 
#        #           shape=19, color=cbPalette[4], size=5) +
#        #geom_point(data=edf2.r11, aes(x=nc_M, y=NPP_M, fill="C"), 
#        #           shape=21, color="black", size=5) +
#        #geom_point(data=edf2.r9, aes(x=nc_act_M, y=NPP_act_M), 
#        #           shape=19, color=cbPalette[3], size=5) +
#        geom_point(data=edf2.r11, aes(x=nc_VL, y=NPP_VL, fill="E"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf2.r9, aes(x=nc_act_VL, y=NPP_act_VL), 
#                   shape=19, color=cbPalette[7], size=5) +
#        ylim(0.0, 2.0) + 
#        xlim(0.001, 0.005) +
#        labs(x="Leaf N:C Ratio", 
#             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
#        theme_linedraw() +
#        theme(panel.grid.minor=element_blank(),
#              axis.text=element_text(size=14),
#              axis.title=element_text(size=16),
#              legend.text=element_text(size=14),
#              legend.title=element_text(size=16),
#              panel.grid.major=element_line(color="grey")) +
#              scale_colour_manual(name="Constraint line", 
#                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = #cbbPalette[7],
#                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
#              scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = #cbPalette[2], "C" = cbPalette[3],
#                                                            "D" = cbPalette[4], "E" = #cbPalette[7])) 
#    
#    
#    p2<-ggplot() + 
#        geom_line(data=cdf1.r11, aes(x=nc, y=NPP_photo, col="C400")) +   
#        #geom_line(data=cdf1.r9, aes(x=nc, y=NPP_photo, col="C400"), linetype="dotted") +   
#        #geom_line(data=cdf1.r11, aes(x=nc, y=NPP_VL, col="VL")) + 
#        #geom_line(data=cdf1.r9, aes(x=nc, y=NPP_VL, col="VL"), linetype="dotted") + 
#        geom_line(data=cdf1.r11, aes(x=nc, y=NPP_L, col="L")) +  
#        geom_line(data=cdf1.r9, aes(x=nc, y=NPP_L, col="L"), linetype="dotted") +  
#        #geom_line(data=cdf1.r11, aes(x=nc, y=NPP_M, col="M")) +    
#        #geom_line(data=cdf1.r9, aes(x=nc, y=NPP_M, col="M"), linetype="dotted") +            
#        geom_line(data=cdf2.r11, aes(x=nc, y=NPP_photo, col="C800")) +     
#        #geom_line(data=cdf2.r9, aes(x=nc, y=NPP_photo, col="C800"), linetype="dotted") +     
#        geom_point(data=edf1.r11, aes(x=nc_VL, y=NPP_I, fill="B"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf1.r9, aes(x=nc_act_VL, y=NPP_I), 
#                   shape=19, color=cbPalette[2], size=5) +
#        geom_point(data=edf1.r11, aes(x=nc_VL, y=NPP_VL, fill="A"), 
#                color="black", shape=21, size=5) + 
#        geom_point(data=edf1.r9, aes(x=nc_act_VL, y=NPP_act_VL), 
#                color=cbPalette[1], shape=19, size=5) + 
#        geom_point(data=edf2.r11, aes(x=nc_L, y=NPP_L, fill="D"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf2.r9, aes(x=nc_act_L, y=NPP_act_L), 
#                   shape=19, color=cbPalette[4], size=5) +
#        #geom_point(data=edf2.r11, aes(x=nc_M, y=NPP_M, fill="C"), 
#        #           shape=21, color="black", size=5) +
#        #geom_point(data=edf2.r9, aes(x=nc_act_M, y=NPP_act_M), 
#        #           shape=19, color=cbPalette[3], size=5) +
#        #geom_point(data=edf2.r11, aes(x=nc_VL, y=NPP_VL, fill="E"), 
#        #           shape=21, color="black", size=5) +
#        #geom_point(data=edf2.r9, aes(x=nc_act_VL, y=NPP_act_VL), 
#        #           shape=19, color=cbPalette[7], size=5) +
#        ylim(0.0, 2.0) + 
#        xlim(0.001, 0.005) +
#        labs(x="Leaf N:C Ratio", 
#             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
#        theme_linedraw() +
#        theme(panel.grid.minor=element_blank(),
#              axis.text=element_text(size=14),
#              axis.title=element_text(size=16),
#              legend.text=element_text(size=14),
#              legend.title=element_text(size=16),
#              panel.grid.major=element_line(color="grey")) +
#              scale_colour_manual(name="Constraint line", 
#                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = #cbbPalette[7],
#                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
#              scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = #cbPalette[2], "C" = cbPalette[3],
#                                                            "D" = cbPalette[4], "E" = #cbPalette[7])) 
#
#    
#    p3<-ggplot() + 
#        geom_line(data=cdf1.r11, aes(x=nc, y=NPP_photo, col="C400")) +   
#        #geom_line(data=cdf1.r9, aes(x=nc, y=NPP_photo, col="C400"), linetype="dotted") +   
#        #geom_line(data=cdf1.r11, aes(x=nc, y=NPP_VL, col="VL")) + 
#        #geom_line(data=cdf1.r9, aes(x=nc, y=NPP_VL, col="VL"), linetype="dotted") + 
#        #geom_line(data=cdf1.r11, aes(x=nc, y=NPP_L, col="L")) +  
#        #geom_line(data=cdf1.r9, aes(x=nc, y=NPP_L, col="L"), linetype="dotted") +  
#        geom_line(data=cdf1.r11, aes(x=nc, y=NPP_M, col="M")) +    
#        geom_line(data=cdf1.r9, aes(x=nc, y=NPP_M, col="M"), linetype="dotted") +            
#        geom_line(data=cdf2.r11, aes(x=nc, y=NPP_photo, col="C800")) +     
#        #geom_line(data=cdf2.r9, aes(x=nc, y=NPP_photo, col="C800"), linetype="dotted") +     
#        geom_point(data=edf1.r11, aes(x=nc_VL, y=NPP_I, fill="B"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf1.r9, aes(x=nc_act_VL, y=NPP_I), 
#                   shape=19, color=cbPalette[2], size=5) +
#        geom_point(data=edf1.r11, aes(x=nc_VL, y=NPP_VL, fill="A"), 
#                color="black", shape=21, size=5) + 
#        geom_point(data=edf1.r9, aes(x=nc_act_VL, y=NPP_act_VL), 
#                color=cbPalette[1], shape=19, size=5) + 
#        #geom_point(data=edf2.r11, aes(x=nc_L, y=NPP_L, fill="D"), 
#        #           shape=21, color="black", size=5) +
#        #geom_point(data=edf2.r9, aes(x=nc_act_L, y=NPP_act_L), 
#        #           shape=19, color=cbPalette[4], size=5) +
#        geom_point(data=edf2.r11, aes(x=nc_M, y=NPP_M, fill="C"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf2.r9, aes(x=nc_act_M, y=NPP_act_M), 
#                   shape=19, color=cbPalette[3], size=5) +
#        #geom_point(data=edf2.r11, aes(x=nc_VL, y=NPP_VL, fill="E"), 
#        #           shape=21, color="black", size=5) +
#        #geom_point(data=edf2.r9, aes(x=nc_act_VL, y=NPP_act_VL), 
#        #           shape=19, color=cbPalette[7], size=5) +
#        ylim(0.0, 2.0) + 
#        xlim(0.001, 0.005) +
#        labs(x="Leaf N:C Ratio", 
#             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
#        theme_linedraw() +
#        theme(panel.grid.minor=element_blank(),
#              axis.text=element_text(size=14),
#              axis.title=element_text(size=16),
#              legend.text=element_text(size=14),
#              legend.title=element_text(size=16),
#              panel.grid.major=element_line(color="grey")) +
#              scale_colour_manual(name="Constraint line", 
#                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = #cbbPalette[7],
#                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
#              scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = #cbPalette[2], "C" = cbPalette[3],
#                                                            "D" = cbPalette[4], "E" = #cbPalette[7]))
#    
#    plot(p1)
#    
#    plot(p2)
#    
#    plot(p3)
#        
#    if (save.plot) {
#        pdf("Plots/CLM_vs_CLM_FUN_fixed_wood_1.pdf", width=8, height=6)
#        plot(p1)
#        dev.off()
#        
#        pdf("Plots/CLM_vs_CLM_FUN_fixed_wood_2.pdf", width=8, height=6)
#        plot(p2)
#        dev.off()
#        
#        pdf("Plots/CLM_vs_CLM_FUN_fixed_wood_3.pdf", width=8, height=6)
#        plot(p3)
#        dev.off()
#    } 
#    

```



#### 5.6.2 variable wood NC
``` {r, echo=FALSE, warning=F, message=F, results="hide"}
#
#    #### Run 10
#    ### split the data.list into two data frames
#    cdf.r10 <- r10$cDF
#    edf.r10 <- r10$eDF
#
#    ### extract CO2 conc. from constraint DF
#    cdf1.r10 <- cdf.r10[cdf.r10$CO2 == 400, ]
#    cdf2.r10 <- cdf.r10[cdf.r10$CO2 == 800, ]
#
#    ### extract CO2 conc. from eq DF
#    edf1.r10 <- edf.r10[edf.r10$CO2 == 400, ]
#    edf2.r10 <- edf.r10[edf.r10$CO2 == 800, ]
#    
#    #### Run 11
#    ### split the data.list into two data frames
#    cdf.r12 <- r12$cDF
#    edf.r12 <- r12$eDF
#
#    ### extract CO2 conc. from constraint DF
#    cdf1.r12 <- cdf.r12[cdf.r12$CO2 == 400, ]
#    cdf2.r12 <- cdf.r12[cdf.r12$CO2 == 800, ]
#
#    ### extract CO2 conc. from eq DF
#    edf1.r12 <- edf.r12[edf.r12$CO2 == 400, ]
#    edf2.r12 <- edf.r12[edf.r12$CO2 == 800, ]
#    
#    ### plot
#    p1<-ggplot() + 
#        geom_line(data=cdf1.r12, aes(x=nc, y=NPP_photo, col="C400")) +   
#        #geom_line(data=cdf1.r10, aes(x=nc, y=NPP_photo, col="C400"), linetype="dotted") +   
#        geom_line(data=cdf1.r12, aes(x=nc, y=NPP_VL, col="VL")) + 
#        geom_line(data=cdf1.r10, aes(x=nc, y=NPP_VL, col="VL"), linetype="dotted") + 
#        #geom_line(data=cdf1.r12, aes(x=nc, y=NPP_L, col="L")) +  
#        #geom_line(data=cdf1.r10, aes(x=nc, y=NPP_L, col="L"), linetype="dotted") +  
#        #geom_line(data=cdf1.r12, aes(x=nc, y=NPP_M, col="M")) +    
#        #geom_line(data=cdf1.r10, aes(x=nc, y=NPP_M, col="M"), linetype="dotted") +            
#        geom_line(data=cdf2.r12, aes(x=nc, y=NPP_photo, col="C800")) +     
#        #geom_line(data=cdf2.r10, aes(x=nc, y=NPP_photo, col="C800"), linetype="dotted") +     
#        geom_point(data=edf1.r12, aes(x=nc_VL, y=NPP_I, fill="B"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf1.r10, aes(x=nc_act_VL, y=NPP_I), 
#                   shape=19, color=cbPalette[2], size=5) +
#        geom_point(data=edf1.r12, aes(x=nc_VL, y=NPP_VL, fill="A"), 
#                color="black", shape=21, size=5) + 
#        geom_point(data=edf1.r10, aes(x=nc_act_VL, y=NPP_act_VL), 
#                color=cbPalette[1], shape=19, size=5) + 
#        #geom_point(data=edf2.r12, aes(x=nc_L, y=NPP_L, fill="D"), 
#        #           shape=21, color="black", size=5) +
#        #geom_point(data=edf2.r10, aes(x=nc_act_L, y=NPP_act_L), 
#        #           shape=19, color=cbPalette[4], size=5) +
#        #geom_point(data=edf2.r12, aes(x=nc_M, y=NPP_M, fill="C"), 
#        #           shape=21, color="black", size=5) +
#        #geom_point(data=edf2.r10, aes(x=nc_act_M, y=NPP_act_M), 
#        #           shape=19, color=cbPalette[3], size=5) +
#        geom_point(data=edf2.r12, aes(x=nc_VL, y=NPP_VL, fill="E"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf2.r10, aes(x=nc_act_VL, y=NPP_act_VL), 
#                   shape=19, color=cbPalette[7], size=5) +
#        ylim(0.0, 2.0) + 
#        xlim(0.001, 0.01) +
#        labs(x="Leaf N:C Ratio", 
#             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
#        theme_linedraw() +
#        theme(panel.grid.minor=element_blank(),
#              axis.text=element_text(size=14),
#              axis.title=element_text(size=16),
#              legend.text=element_text(size=14),
#              legend.title=element_text(size=16),
#              panel.grid.major=element_line(color="grey")) +
#              scale_colour_manual(name="Constraint line", 
#                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = #cbbPalette[7],
#                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
#              scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = #cbPalette[2], "C" = cbPalette[3],
#                                                            "D" = cbPalette[4], "E" = #cbPalette[7])) 
#    
#    
#    p2<-ggplot() + 
#        geom_line(data=cdf1.r12, aes(x=nc, y=NPP_photo, col="C400")) +   
#        #geom_line(data=cdf1.r10, aes(x=nc, y=NPP_photo, col="C400"), linetype="dotted") +   
#        #geom_line(data=cdf1.r12, aes(x=nc, y=NPP_VL, col="VL")) + 
#        #geom_line(data=cdf1.r10, aes(x=nc, y=NPP_VL, col="VL"), linetype="dotted") + 
#        geom_line(data=cdf1.r12, aes(x=nc, y=NPP_L, col="L")) +  
#        geom_line(data=cdf1.r10, aes(x=nc, y=NPP_L, col="L"), linetype="dotted") +  
#        #geom_line(data=cdf1.r12, aes(x=nc, y=NPP_M, col="M")) +    
#        #geom_line(data=cdf1.r10, aes(x=nc, y=NPP_M, col="M"), linetype="dotted") +            
#        geom_line(data=cdf2.r12, aes(x=nc, y=NPP_photo, col="C800")) +     
#        #geom_line(data=cdf2.r10, aes(x=nc, y=NPP_photo, col="C800"), linetype="dotted") +     
#        geom_point(data=edf1.r12, aes(x=nc_VL, y=NPP_I, fill="B"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf1.r10, aes(x=nc_act_VL, y=NPP_I), 
#                   shape=19, color=cbPalette[2], size=5) +
#        geom_point(data=edf1.r12, aes(x=nc_VL, y=NPP_VL, fill="A"), 
#                color="black", shape=21, size=5) + 
#        geom_point(data=edf1.r10, aes(x=nc_act_VL, y=NPP_act_VL), 
#                color=cbPalette[1], shape=19, size=5) + 
#        geom_point(data=edf2.r12, aes(x=nc_L, y=NPP_L, fill="D"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf2.r10, aes(x=nc_act_L, y=NPP_act_L), 
#                   shape=19, color=cbPalette[4], size=5) +
#        #geom_point(data=edf2.r12, aes(x=nc_M, y=NPP_M, fill="C"), 
#        #           shape=21, color="black", size=5) +
#        #geom_point(data=edf2.r10, aes(x=nc_act_M, y=NPP_act_M), 
#        #           shape=19, color=cbPalette[3], size=5) +
#        #geom_point(data=edf2.r12, aes(x=nc_VL, y=NPP_VL, fill="E"), 
#        #           shape=21, color="black", size=5) +
#        #geom_point(data=edf2.r10, aes(x=nc_act_VL, y=NPP_act_VL), 
#        #           shape=19, color=cbPalette[7], size=5) +
#        ylim(0.0, 2.0) + 
#        xlim(0.001, 0.01) +
#        labs(x="Leaf N:C Ratio", 
#             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
#        theme_linedraw() +
#        theme(panel.grid.minor=element_blank(),
#              axis.text=element_text(size=14),
#              axis.title=element_text(size=16),
#              legend.text=element_text(size=14),
#              legend.title=element_text(size=16),
#              panel.grid.major=element_line(color="grey")) +
#              scale_colour_manual(name="Constraint line", 
#                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = #cbbPalette[7],
#                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
#              scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = #cbPalette[2], "C" = cbPalette[3],
#                                                            "D" = cbPalette[4], "E" = #cbPalette[7])) 
#
#    
#    p3<-ggplot() + 
#        geom_line(data=cdf1.r12, aes(x=nc, y=NPP_photo, col="C400")) +   
#        #geom_line(data=cdf1.r10, aes(x=nc, y=NPP_photo, col="C400"), linetype="dotted") +   
#        #geom_line(data=cdf1.r12, aes(x=nc, y=NPP_VL, col="VL")) + 
#        #geom_line(data=cdf1.r10, aes(x=nc, y=NPP_VL, col="VL"), linetype="dotted") + 
#        #geom_line(data=cdf1.r12, aes(x=nc, y=NPP_L, col="L")) +  
#        #geom_line(data=cdf1.r10, aes(x=nc, y=NPP_L, col="L"), linetype="dotted") +  
#        geom_line(data=cdf1.r12, aes(x=nc, y=NPP_M, col="M")) +    
#        geom_line(data=cdf1.r10, aes(x=nc, y=NPP_M, col="M"), linetype="dotted") +            
#        geom_line(data=cdf2.r12, aes(x=nc, y=NPP_photo, col="C800")) +     
#        #geom_line(data=cdf2.r10, aes(x=nc, y=NPP_photo, col="C800"), linetype="dotted") +     
#        geom_point(data=edf1.r12, aes(x=nc_VL, y=NPP_I, fill="B"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf1.r10, aes(x=nc_act_VL, y=NPP_I), 
#                   shape=19, color=cbPalette[2], size=5) +
#        geom_point(data=edf1.r12, aes(x=nc_VL, y=NPP_VL, fill="A"), 
#                color="black", shape=21, size=5) + 
#        geom_point(data=edf1.r10, aes(x=nc_act_VL, y=NPP_act_VL), 
#                color=cbPalette[1], shape=19, size=5) + 
#        #geom_point(data=edf2.r12, aes(x=nc_L, y=NPP_L, fill="D"), 
#        #           shape=21, color="black", size=5) +
#        #geom_point(data=edf2.r10, aes(x=nc_act_L, y=NPP_act_L), 
#        #           shape=19, color=cbPalette[4], size=5) +
#        geom_point(data=edf2.r12, aes(x=nc_M, y=NPP_M, fill="C"), 
#                   shape=21, color="black", size=5) +
#        geom_point(data=edf2.r10, aes(x=nc_act_M, y=NPP_act_M), 
#                   shape=19, color=cbPalette[3], size=5) +
#        #geom_point(data=edf2.r12, aes(x=nc_VL, y=NPP_VL, fill="E"), 
#        #           shape=21, color="black", size=5) +
#        #geom_point(data=edf2.r10, aes(x=nc_act_VL, y=NPP_act_VL), 
#        #           shape=19, color=cbPalette[7], size=5) +
#        ylim(0.0, 2.0) + 
#        xlim(0.001, 0.01) +
#        labs(x="Leaf N:C Ratio", 
#             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
#        theme_linedraw() +
#        theme(panel.grid.minor=element_blank(),
#              axis.text=element_text(size=14),
#              axis.title=element_text(size=16),
#              legend.text=element_text(size=14),
#              legend.title=element_text(size=16),
#              panel.grid.major=element_line(color="grey")) +
#              scale_colour_manual(name="Constraint line", 
#                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], "VL" = #cbbPalette[7],
#                                       "L" = cbbPalette[4], "M" = cbbPalette[3])) +
#              scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = #cbPalette[2], "C" = cbPalette[3],
#                                                            "D" = cbPalette[4], "E" = #cbPalette[7]))
#    
#    plot(p1)
#    
#    plot(p2)
#    
#    plot(p3)
#        
#    if (save.plot) {
#        pdf("Plots/CLM_vs_CLM_FUN_variable_wood_1.pdf", width=8, height=6)
#        plot(p1)
#        dev.off()
#        
#        pdf("Plots/CLM_vs_CLM_FUN_variable_wood_2.pdf", width=8, height=6)
#        plot(p2)
#        dev.off()
#        
#        pdf("Plots/CLM_vs_CLM_FUN_variable_wood_3.pdf", width=8, height=6)
#        plot(p3)
#        dev.off()
#    } 
    

```

### 5.7 Medlyn and Dewar allocation coupling
``` {r, echo=FALSE, warning=F, message=F}

    sensDF <- matrix(nrow=4*6, ncol=6)
    colnames(sensDF) <- c("coupling", "af", "ar", "aw", "value", "Time")
    sensDF <- as.data.frame(sensDF)
    sensDF$coupling <- rep(c("no_coupling", "linear_coupling"), each = 3*4)
    sensDF$Time <- rep(c("I", "M", "L", "VL"), by = 6)
    
    sensDF$af <- c(rep(0.18, each=4), rep(0.2, each=4), rep(0.22, each=4),
                   rep(0.18, each=4), rep(0.2, each=4), rep(0.22, each=4))
    sensDF$aw <- c(rep(0.6, each=3*4), rep(0.54, each=4), rep(0.6, each=4), rep(0.66, each=4))
    sensDF$ar <- 1.0 - sensDF$aw - sensDF$af

    sensDF[sensDF$af==0.18&sensDF$aw==0.6&sensDF$Time=="I","value"] <- (r20$eDF$NPP_I[2]-r20$eDF$NPP_VL[1])/r20$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.18&sensDF$aw==0.6&sensDF$Time=="M","value"] <- (r20$eDF$NPP_M[2]-r20$eDF$NPP_VL[1])/r20$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.18&sensDF$aw==0.6&sensDF$Time=="L","value"] <- (r20$eDF$NPP_L[2]-r20$eDF$NPP_VL[1])/r20$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.18&sensDF$aw==0.6&sensDF$Time=="VL","value"] <- (r20$eDF$NPP_VL[2]-r20$eDF$NPP_VL[1])/r20$eDF$NPP_VL[1]

    sensDF[sensDF$af==0.2&sensDF$aw==0.6&sensDF$coupling=="no_coupling"&sensDF$Time=="I","value"] <- (r21$eDF$NPP_I[2]-r21$eDF$NPP_VL[1])/r21$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.2&sensDF$aw==0.6&sensDF$coupling=="no_coupling"&sensDF$Time=="M","value"] <- (r21$eDF$NPP_M[2]-r21$eDF$NPP_VL[1])/r21$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.2&sensDF$aw==0.6&sensDF$coupling=="no_coupling"&sensDF$Time=="L","value"] <- (r21$eDF$NPP_L[2]-r21$eDF$NPP_VL[1])/r21$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.2&sensDF$aw==0.6&sensDF$coupling=="no_coupling"&sensDF$Time=="VL","value"] <- (r21$eDF$NPP_VL[2]-r21$eDF$NPP_VL[1])/r21$eDF$NPP_VL[1]
    
    sensDF[sensDF$af==0.22&sensDF$aw==0.6&sensDF$Time=="I","value"] <- (r22$eDF$NPP_I[2]-r22$eDF$NPP_VL[1])/r22$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.22&sensDF$aw==0.6&sensDF$Time=="M","value"] <- (r22$eDF$NPP_M[2]-r22$eDF$NPP_VL[1])/r22$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.22&sensDF$aw==0.6&sensDF$Time=="L","value"] <- (r22$eDF$NPP_L[2]-r22$eDF$NPP_VL[1])/r22$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.22&sensDF$aw==0.6&sensDF$Time=="VL","value"] <- (r22$eDF$NPP_VL[2]-r22$eDF$NPP_VL[1])/r22$eDF$NPP_VL[1]
    
    sensDF[sensDF$af==0.18&sensDF$aw==0.54&sensDF$Time=="I","value"] <- (r23$eDF$NPP_I[2]-r23$eDF$NPP_VL[1])/r23$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.18&sensDF$aw==0.54&sensDF$Time=="M","value"] <- (r23$eDF$NPP_M[2]-r23$eDF$NPP_VL[1])/r23$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.18&sensDF$aw==0.54&sensDF$Time=="L","value"] <- (r23$eDF$NPP_L[2]-r23$eDF$NPP_VL[1])/r23$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.18&sensDF$aw==0.54&sensDF$Time=="VL","value"] <- (r23$eDF$NPP_VL[2]-r23$eDF$NPP_VL[1])/r23$eDF$NPP_VL[1]
    
    sensDF[sensDF$af==0.2&sensDF$coupling=="linear_coupling"&sensDF$Time=="I","value"] <- (r24$eDF$NPP_I[2]-r24$eDF$NPP_VL[1])/r24$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.2&sensDF$coupling=="linear_coupling"&sensDF$Time=="M","value"] <- (r24$eDF$NPP_M[2]-r24$eDF$NPP_VL[1])/r24$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.2&sensDF$coupling=="linear_coupling"&sensDF$Time=="L","value"] <- (r24$eDF$NPP_L[2]-r24$eDF$NPP_VL[1])/r24$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.2&sensDF$coupling=="linear_coupling"&sensDF$Time=="VL","value"] <- (r24$eDF$NPP_VL[2]-r24$eDF$NPP_VL[1])/r24$eDF$NPP_VL[1]
    
    sensDF[sensDF$af==0.22&sensDF$aw==0.66&sensDF$Time=="I","value"] <- (r25$eDF$NPP_I[2]-r25$eDF$NPP_VL[1])/r25$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.22&sensDF$aw==0.66&sensDF$Time=="M","value"] <- (r25$eDF$NPP_M[2]-r25$eDF$NPP_VL[1])/r25$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.22&sensDF$aw==0.66&sensDF$Time=="L","value"] <- (r25$eDF$NPP_L[2]-r25$eDF$NPP_VL[1])/r25$eDF$NPP_VL[1]
    sensDF[sensDF$af==0.22&sensDF$aw==0.66&sensDF$Time=="VL","value"] <- (r25$eDF$NPP_VL[2]-r25$eDF$NPP_VL[1])/r25$eDF$NPP_VL[1]
    
    sensDF$af <- as.character(sensDF$af)
    sensDF$aw <- as.character(sensDF$aw)
    sensDF$value <- sensDF$value * 100
    
    plotDF1 <- subset(sensDF, coupling=="no_coupling")
    plotDF2 <- subset(sensDF, coupling=="linear_coupling")
    plotDF3 <- subset(sensDF, Time=="VL")

    ### plot
    p1<-ggplot() + 
        geom_bar(data=plotDF1, aes(x=af, y=value, fill=Time), 
                 stat="identity", position=position_dodge()) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              axis.text.x = element_text(angle = 90, hjust=1),
              panel.grid.major=element_line(color="grey")) + 
        ylab("% change") + xlab("") + 
        scale_fill_manual(name="Timescale", values = c("I" = cbPalette[2], "M" = cbPalette[3],
                                                      "L" = cbPalette[4], "VL" = cbPalette[7])) 

    p2<-ggplot() + 
        geom_bar(data=plotDF2, aes(x=af, y=value, fill=Time), 
                 stat="identity", position=position_dodge()) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              axis.text.x = element_text(angle = 90, hjust=1),
              panel.grid.major=element_line(color="grey")) + 
        ylab("% change") + xlab("") + 
        scale_fill_manual(name="Timescale", values = c("I" = cbPalette[2], "M" = cbPalette[3],
                                                      "L" = cbPalette[4], "VL" = cbPalette[7])) 
    
    p3<-ggplot() + 
        geom_bar(data=plotDF3, aes(x=af, y=value, fill=coupling), 
                 stat="identity", position=position_dodge()) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              axis.text.x = element_text(angle = 90, hjust=1),
              panel.grid.major=element_line(color="grey")) + 
        ylab("% change") + xlab("") + 
        scale_fill_manual(name="Coupling", values = c("no_coupling" = cbPalette[2], "linear_coupling" = cbPalette[3])) 

    if (save.plot) {
        pdf("Plots/sensitivity_Medlyn_Dewar_coupling.pdf", width=10, height=6)
        plot(p3)
        dev.off()
    } 
    

```


### 5.8 Medlyn and Dewar allocation coupling simplified
``` {r, echo=FALSE, warning=F, message=F}

    ### split the data.list into two data frames
    cdf.r26 <- r26$cDF
    edf.r26 <- r26$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r26 <- cdf.r26[cdf.r26$CO2 == 400, ]
    cdf2.r26 <- cdf.r26[cdf.r26$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r26 <- edf.r26[edf.r26$CO2 == 400, ]
    edf2.r26 <- edf.r26[edf.r26$CO2 == 800, ]
    
    ### split the data.list into two data frames
    cdf.r27 <- r27$cDF
    edf.r27 <- r27$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r27 <- cdf.r27[cdf.r27$CO2 == 400, ]
    cdf2.r27 <- cdf.r27[cdf.r27$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r27 <- edf.r27[edf.r27$CO2 == 400, ]
    edf2.r27 <- edf.r27[edf.r27$CO2 == 800, ]
    
    
    ### split the data.list into two data frames
    cdf.r28 <- r28$cDF
    edf.r28 <- r28$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r28 <- cdf.r28[cdf.r28$CO2 == 400, ]
    cdf2.r28 <- cdf.r28[cdf.r28$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r28 <- edf.r28[edf.r28$CO2 == 400, ]
    edf2.r28 <- edf.r28[edf.r28$CO2 == 800, ]
    
    
    ### split the data.list into two data frames
    cdf.r29 <- r29$cDF
    edf.r29 <- r29$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r29 <- cdf.r29[cdf.r29$CO2 == 400, ]
    cdf2.r29 <- cdf.r29[cdf.r29$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r29 <- edf.r29[edf.r29$CO2 == 400, ]
    edf2.r29 <- edf.r29[edf.r29$CO2 == 800, ]
    
    ### plot
    p1<-ggplot() + 
        geom_line(data=cdf1.r26, aes(x=nc, y=NPP_photo, col="C400")) +   
        geom_line(data=cdf1.r26, aes(x=nc, y=NPP_L, col="L")) +  
        geom_line(data=cdf2.r26, aes(x=nc, y=NPP_photo, col="C800")) +   
        geom_line(data=cdf1.r28, aes(x=nc, y=NPP_L, col="L_10")) +  
        geom_point(data=edf2.r26, aes(x=nc_L, y=NPP_L, fill="B"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf1.r26, aes(x=nc_L, y=NPP_L, fill="A"), 
                   shape=21, color="black", size=5) +
                geom_point(data=edf2.r28, aes(x=nc_L, y=NPP_L, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf1.r28, aes(x=nc_L, y=NPP_L, fill="C"), 
                   shape=21, color="black", size=5) +
        ylim(1.0, 2.5) + 
        xlim(0.001, 0.1) +
        labs(x="Leaf N:C Ratio", 
             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], 
                                       "L" = cbbPalette[4], "L_10" = cbbPalette[3]),
                            labels = c("C400", "C800", "L", "L+10%")) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4])) 
    
    
    p2<-ggplot() + 
        geom_line(data=cdf1.r27, aes(x=nc, y=NPP_photo, col="C400")) +   
        geom_line(data=cdf1.r27, aes(x=nc, y=NPP_L, col="L")) +  
        geom_line(data=cdf2.r27, aes(x=nc, y=NPP_photo, col="C800")) +   
        geom_line(data=cdf1.r29, aes(x=nc, y=NPP_L, col="L_10")) +  
        geom_point(data=edf2.r27, aes(x=nc_L, y=NPP_L, fill="B"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf1.r27, aes(x=nc_L, y=NPP_L, fill="A"), 
                   shape=21, color="black", size=5) +
                geom_point(data=edf2.r29, aes(x=nc_L, y=NPP_L, fill="D"), 
                   shape=21, color="black", size=5) +
        geom_point(data=edf1.r29, aes(x=nc_L, y=NPP_L, fill="C"), 
                   shape=21, color="black", size=5) +
        ylim(1.0, 2.5) + 
        xlim(0.001, 0.1) +
        labs(x="Leaf N:C Ratio", 
             y=expression(paste("NPP [kg C ", m^-2, " ", yr^-1, "]"))) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              panel.grid.major=element_line(color="grey")) +
        scale_colour_manual(name="Constraint line", 
                            values = c("C400" = cbbPalette[1], "C800" = cbbPalette[2], 
                                       "L" = cbbPalette[4], "L_10" = cbbPalette[3]),
                            labels = c("C400", "C800", "L", "L+10%")) +
        scale_fill_manual(name="QE point", values = c("A" = cbPalette[1], "B" = cbPalette[2], "C" = cbPalette[3],
                                                      "D" = cbPalette[4])) 

    plot(p2)
  
    if (save.plot) {
        pdf("Plots/sensitivity_Medlyn_Dewar_no_coupling.pdf", width=10, height=6)
        plot(p1)
        dev.off()
        
        pdf("Plots/sensitivity_Medlyn_Dewar_linear_coupling.pdf", width=10, height=6)
        plot(p2)
        dev.off()
    } 
    

```


### 5.9 Kirschbaum et al. 1994
``` {r, echo=FALSE, warning=F, message=F}

    ### split the data.list into two data frames
    cdf.r1.1 <- r1.1$cDF
    edf.r1.1 <- r1.1$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r1.1 <- cdf.r1.1[cdf.r1.1$CO2 == 400, ]
    cdf2.r1.1 <- cdf.r1.1[cdf.r1.1$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r1.1 <- edf.r1.1[edf.r1.1$CO2 == 400, ]
    edf2.r1.1 <- edf.r1.1[edf.r1.1$CO2 == 800, ]
    
    ### split the data.list into two data frames
    cdf.r1.2 <- r1.2$cDF
    edf.r1.2 <- r1.2$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r1.2 <- cdf.r1.2[cdf.r1.2$CO2 == 400, ]
    cdf2.r1.2 <- cdf.r1.2[cdf.r1.2$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r1.2 <- edf.r1.2[edf.r1.2$CO2 == 400, ]
    edf2.r1.2 <- edf.r1.2[edf.r1.2$CO2 == 800, ]
    
    ### split the data.list into two data frames
    cdf.r1 <- r1$cDF
    edf.r1 <- r1$eDF

    ### extract CO2 conc. from constraint DF
    cdf1.r1 <- cdf.r1[cdf.r1$CO2 == 400, ]
    cdf2.r1 <- cdf.r1[cdf.r1$CO2 == 800, ]

    ### extract CO2 conc. from eq DF
    edf1.r1 <- edf.r1[edf.r1$CO2 == 400, ]
    edf2.r1 <- edf.r1[edf.r1$CO2 == 800, ]
  
    ### Create a storage df
    sensDF <- matrix(nrow=4*3, ncol=3)
    colnames(sensDF) <- c("n_trans", "value", "Time")
    sensDF <- as.data.frame(sensDF)
    sensDF$n_trans <- rep(c(0.1, 0.5, 0.9), each = 4)
    sensDF$Time <- rep(c("I", "M", "L", "VL"), by = 3)
    

    sensDF[sensDF$n_trans==0.1&sensDF$Time=="I","value"] <- (r1.1$eDF$NPP_I[2]-r1.1$eDF$NPP_VL[1])/r1.1$eDF$NPP_VL[1]
    sensDF[sensDF$n_trans==0.1&sensDF$Time=="M","value"] <- (r1.1$eDF$NPP_M[2]-r1.1$eDF$NPP_VL[1])/r1.1$eDF$NPP_VL[1]
    sensDF[sensDF$n_trans==0.1&sensDF$Time=="L","value"] <- (r1.1$eDF$NPP_L[2]-r1.1$eDF$NPP_VL[1])/r1.1$eDF$NPP_VL[1]
    sensDF[sensDF$n_trans==0.1&sensDF$Time=="VL","value"] <- (r1.1$eDF$NPP_VL[2]-r1.1$eDF$NPP_VL[1])/r1.1$eDF$NPP_VL[1]
    
    
    sensDF[sensDF$n_trans==0.5&sensDF$Time=="I","value"] <- (r1$eDF$NPP_I[2]-r1$eDF$NPP_VL[1])/r1$eDF$NPP_VL[1]
    sensDF[sensDF$n_trans==0.5&sensDF$Time=="M","value"] <- (r1$eDF$NPP_M[2]-r1$eDF$NPP_VL[1])/r1$eDF$NPP_VL[1]
    sensDF[sensDF$n_trans==0.5&sensDF$Time=="L","value"] <- (r1$eDF$NPP_L[2]-r1$eDF$NPP_VL[1])/r1$eDF$NPP_VL[1]
    sensDF[sensDF$n_trans==0.5&sensDF$Time=="VL","value"] <- (r1$eDF$NPP_VL[2]-r1$eDF$NPP_VL[1])/r1$eDF$NPP_VL[1]
    
    sensDF[sensDF$n_trans==0.9&sensDF$Time=="I","value"] <- (r1.2$eDF$NPP_I[2]-r1.2$eDF$NPP_VL[1])/r1.2$eDF$NPP_VL[1]
    sensDF[sensDF$n_trans==0.9&sensDF$Time=="M","value"] <- (r1.2$eDF$NPP_M[2]-r1.2$eDF$NPP_VL[1])/r1.2$eDF$NPP_VL[1]
    sensDF[sensDF$n_trans==0.9&sensDF$Time=="L","value"] <- (r1.2$eDF$NPP_L[2]-r1.2$eDF$NPP_VL[1])/r1.2$eDF$NPP_VL[1]
    sensDF[sensDF$n_trans==0.9&sensDF$Time=="VL","value"] <- (r1.2$eDF$NPP_VL[2]-r1.2$eDF$NPP_VL[1])/r1.2$eDF$NPP_VL[1]
    
    sensDF$n_trans <- as.character(sensDF$n_trans)
    
    ### Plotting
    p1<-ggplot() + 
        geom_bar(data=sensDF, aes(x=n_trans, y=value, fill=Time), 
                 stat="identity", position=position_dodge()) +
        theme_linedraw() +
        theme(panel.grid.minor=element_blank(),
              axis.text=element_text(size=14),
              axis.title=element_text(size=16),
              legend.text=element_text(size=14),
              legend.title=element_text(size=16),
              axis.text.x = element_text(angle = 90, hjust=1),
              panel.grid.major=element_line(color="grey")) + 
        ylab("% change") + xlab("") + 
        scale_fill_manual(name="Timescale", values = c("I" = cbPalette[2], "M" = cbPalette[3],
                                                      "L" = cbPalette[4], "VL" = cbPalette[7])) 
    

    if (save.plot) {
        pdf("Plots/sensitivity_Kirshbaum_retranslocation.pdf", width=10, height=6)
        plot(p1)
        dev.off()
    
    } 
    

```

# 6. To-do list:
Below is a continuous update of the to-do list of this project:

    - Check N fixation (McMurtrie 2000)
    - Implement C cost of N uptake properly (FUN model), right now problem is L and VL term NPP response ratio is higher under fixed wood assumption as compared to the model without FUN model. We introduced c cost to acquire N, and at L term, c cost is reduced hence more response? Does it make sense?
    - NUE plotting, understand difference compared to McMurtrie 2000 and Zaehle 2014

<br/>
<br/>