Rewrite-Log-Rick.Rmd

---
title: "Rewrite Log Rick"
author: "Rick Venema"
date: "`r Sys.Date()`"
output:
  pdf_document
---

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


# Introduction
The model we have recieved from Remi Cardinael from Cardinael et al 2018 was based around the SOC stocks in an Agroforestry plot containing walnut trees and durum wheat. 

## References
1. Cardinael, R., Guenet, B., Chevallier, T., Dupraz, C., Cozzi, T., and Chenu, C.: High organic imputs explain shallow and deep SOC storage in a long-term agroforestry system - combining experimental and modeling approaches, Biogeosciences, 15, 297-317, https://di.org/10.5194/bg-15-297-2018, 2018.

# 7-6-2018
This was the first day we started rewriting the code given by Cardinael. This code was around 948 lines long. 
The model was a hard coded model of the reality of that specific tree. We decided that we wanted to rewrite the model. 
This was because of the first intention to add the nitrogen cycle to the model was not easy because of the bad coding we recieved. So in order to add the Nitrogen model, the model needs to be rewritten so it can be used to generate a model based on parameters instead of fixed data.
```{r 7-6-18a, eval=F}
modelp3difft <- function(t, initial_state, parms){
  with (as.list(parms),{
          A <- initial_state[1:dim(z)[1]]
          S <- initial_state[(dim(z)[1]+1):((2*dim(z)[1]))]
          P <- initial_state[((2*dim(z)[1])+1):((3*dim(z)[1]))]  
          
          #Fluxes in z direction
          FluxA <- Dt * (c(0,A))/dz -D * diff(c(0,A,0)) / dz - c(0,Dmix) * diff(c(0,A,0)) / dz 
          FluxS <- Dt * (c(0,S))/dz -D_slow*diff(c(0,S,0)) / dz - c(0,Dmix) * diff(c(0,S,0)) / dz
          FluxP <- Dt * (c(0,P))/dz -D_slow*diff(c(0,P,0)) / dz - c(0,Dmix) * diff(c(0,P,0)) / dz
          FluxA[1]=0.
          FluxS[1]=0.
          FluxP[1]=0.
          
          #Reaction       
          Import<-import_tree_be*mr_tree + import_grass_be*mr_grass + import_crop_be*mr_crop
          dA=-diff(FluxA) + (e*ks*S*frac_SA + e*kp*P - kf* A* clay_func)* mf* tf  + Import
          dS=-diff(FluxS) + (frac_AS*kf*e*A*clay_func - ks*S)* mf* tf
          dP=-diff(FluxP) + ((1-frac_SA)*e*ks*S + (1-frac_AS)*kf*e*A*clay_func - kp*P)* mf* tf
          
          return(list(c(dA=dA,dS=dS,dP=dP)))
        })
}
```
Above, the model is given. This is the model what generates the output.
  
In the original code `PARAM=read.table('run_options_p0_3pools.def', header=FALSE)` can be found. This file contains only numbers and thus can be added to the model as a vector to adjust parameter even more. The PARAM part can also be ommited to give the parameters their respective values directly instead of being loaded from a file and then passed as a dataframe. The original code called the PARAM variable by `parameter = PARAM[i,1]/10000. `. This is unreadable for himself when he tries to debug this code. It is easier to give the parameters directly with their correct values. An example follows below:

```{r 7-6-18b, eval=F}
### Code with PARAM as hard coded vector
PARAM <-  c(0.01, 0.83, 5.24, 21.60, 0.34, 0.99, 0.94)
Dt=PARAM[4]/10000.
D_slow = PARAM[3]/10000.

### Code with values directly added to variable
Dt = 21.60/10000.
D_slow = 5.24/10000.
```

# 8-6-2018
The next part I tackled was the crop_be part of the code. This is the carbon input of the tree, grass, and crop. 
My idea was that it can be different with each situation. But how we could fix it, that its different on each location, is unclear at the moment, I could create a simple model that can be specified to different inputs

To implement this, I first need to understand the code and know what the form of the profile is. If I can implement a type of different input to the model, the model can be used on different locations/situations.
The model can than be used to an agriculture plot with potatoes for example.
  
My first thought is to rewrite the existing input matrixes. The inputs can later be modified to be corrected on the situation.
If there are different types of input models, the model can be tweaked a lot more. The first thing I tackled was crop input, because that part of input I was interested by the most, because I wanted the model to be 
able to be used in agriculture research with wheat (future research by me I hope)s

```{r 8-6-18a, eval=F}

######### NEW CODE ##########
profil_CR_R<-matrix(ncol=dim(d)[1],nrow=dim(z)[1])

## Formula to calculate profile 26.443*exp((-2.6)*(-z[i-1,1]))
create_crop_profile <- function(z){
  return(26.433*exp((-2.6)*(-z))/100)
}

profil_CR_R <- create_crop_profile(z)
profil_CR_R <- as.data.frame(rep(profil_CR_R, dim(d)[1]))
colnames(profil_CR_R) <- as.character(d[,1])
rownames(profil_CR_R) <- as.character(z[,1])
profil_CR_R[z[,1]<= -1.5,] <- 0
rofil_CR_R_SPIN<-profil_CR_R
####### OLD CODE #####


#Roots profil of crop roots (% of the total root mass)
profil_CR_R<-matrix(ncol=dim(d)[1]+1,nrow=dim(z)[1]+1)
for (i in 1:dim(d)[1]) {profil_CR_R[1,i+1]<-d[i,]}
for (i in 1:dim(z)[1]) {profil_CR_R[i+1,1]<-z[i,]}

for (j in seq(z[1,1],z[dim(z)[1],1], by=-step_depth)){ 
  for (i in 1:dim(z)[1]+1) {
    if (profil_CR_R[i,1]==as.character(j)){
      for (k in 1:dim(d)[1]+1)     {
        profil_CR_R[i,k]<-26.443*exp((-2.6)*(-z[i-1,1])) 
        profil_CR_R[i,k]<-profil_CR_R[i,k]/100 #Conversion from % to proportion
        if (profil_CR_R[i,1]<limit_root_crop) {profil_CR_R[i,k]<-0} # no more crop roots below 1.5m
      }
    }
  }
}
profil_CR_R_SPIN<-profil_CR_R
for (i in 1:dim(z)[1]+1) {
  for (k in 1:dim(d)[1]+1)     {
    if (profil_CR_R[1,k]<=limit_grass) {profil_CR_R[i,k]<-0} # no crop on the tree line
  }
}		

```

I have rewritten the input of crop carbon. The old code had too many for loops, this is reduced to none. The build in functions of dataframes can be used to do the same things that Cardinael did with his 7 for loops.
  
In the model the yield of decomposed FOM that goes to SOM is defined by `e = PARAM[5,1]`. This can be rewritten to `e = 0.34`. This is better to read and is better to adjust.