## How do we create a soil pool model in R?

soc <- function(carbon.input, soc.labile, soc.slow, soc.passive, k1 = 0.06, k2 = 0.006){

  ## First we consider the carbon input
  ## Possibly there is a partition to the different pools
  q1 <- 0.7
  q2 <- 0.2
  q3 <- 0.1
  
  add.carbon.input.soc.labile <- carbon.input * q1
  add.carbon.input.soc.slow <- carbon.input * q2
  add.carbon.input.soc.passive <- carbon.input * q3

  soc.labile <- soc.labile + add.carbon.input.soc.labile
  soc.slow <- soc.slow + add.carbon.input.soc.slow
  soc.passive <- soc.passive + add.carbon.input.soc.passive

  ## Let's assume that there is some flow of C from labile to slow

  out.labile <- soc.labile * k1
  soc.labile <- soc.labile - out.labile

  if(soc.labile < 0) stop('k1 might be too high since soc.labile became negative')

  ## Add it to the slow

  soc.slow <- soc.slow + out.labile

  ## Some soc.slow moves to the passive pool

  out.slow <- soc.slow * k2
  soc.slow <- soc.slow - out.slow

  if(soc.slow < 0) stop('k2 might be too high since soc.labile became negative')

  ## Add it to the passive

  soc.passive <- soc.passive + out.slow

  ans <- c(soc.labile, soc.slow, soc.passive)
  names(ans) <- c('soc.labile', 'soc.slow', 'soc.passive')

  ans

}

##

## Testing the function

cinput <- 100 # g/m2
clabile <- 300
cslow <- 200
cpassive <- 100

nyears <- 2000

col.res <- matrix(nrow = nyears, ncol = 3)

for(i in 1:nyears){

  if(i == 1){
    res <- soc(cinput, clabile, cslow, cpassive)
  }else{
    res <- soc(0, clabile, cslow, cpassive)
  }

  clabile <- res[1]
  cslow <- res[2]
  cpassive <- res[3]

  col.res[i,] <- res

}
  
## plotting results

xyplot(col.res[,1] + col.res[,2] + col.res[,3] ~ 1:nyears,
       type = 'l',
       lwd = 3,
       col = c('red', 'purple', 'black'),
       ylab = 'SOC',
       xlab = 'years',
       key = list(text = list(c('labile','slow','passive')),
                    col = c('red', 'purple', 'black'), cex=1.5))
  

## Trying a model with repeated inputs of C  
cinput <- 100
clabile <- 300
cslow <- 200
cpassive <- 100

nyears <- 2000

col.res <- matrix(nrow = nyears, ncol = 3)

for(i in 1:nyears){

    res <- soc(cinput, clabile, cslow, cpassive)

  clabile <- res[1]
  cslow <- res[2]
  cpassive <- res[3]

  col.res[i,] <- res

}


xyplot(col.res[,1] + col.res[,2] + col.res[,3] ~ 1:nyears,
       type = 'l',
       lwd = 3,
#       ylim = c(0,10000),
#       xlim = c(0,2000),
       col = c('red', 'purple', 'black'),
       ylab = 'SOC',
       xlab = 'years',
       key = list(text = list(c('labile','slow','passive')),
                    col = c('red', 'purple', 'black'), cex=1.5))