Models of Soil Organic Matter Decomposition

231

0

This function creates a numerical model for n independent (parallel) pools that can be queried afterwards.
It is used by the convinient wrapper functions `TwopParallelModel`

and `ThreepParallelModel`

but can also be used independently.

1 2 3 | ```
ParallelModel(times,
coeffs_tm, startvalues, inputrates, solverfunc = deSolve.lsoda.wrapper,
pass = FALSE)
``` |

`times` |
A vector containing the points in time where the solution is sought. |

`coeffs_tm` |
A TimeMap object consisting of a vector valued function containing the decay rates for the n pools as function of time and the time range where this function is valid. The length of the vector is equal to the number of pools. |

`startvalues` |
A vector containing the initial amount of carbon for the n pools. <<The length of this vector is equal to the number of pools and thus equal to the length of k. This is checked by the function. |

`inputrates` |
An object consisting of a vector valued function describing the inputs to the pools as funtions of time |

`solverfunc` |
The function used to actually solve the ODE system. This can be |

`pass` |
if TRUE forces the constructor to create the model even if it is invalid |

a model object

Carlos A. Sierra <csierra@bgc-jena.mpg.de>, Markus Mueller <mamueller@bgc-jena.mpg.de>

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 | ```
t_start=0
t_end=10
tn=50
timestep=(t_end-t_start)/tn
t=seq(t_start,t_end,timestep)
k=TimeMap.new(t_start,t_end,function(times){c(-0.5,-0.2,-0.3)})
c0=c(1, 2, 3)
#constant inputrates
inputrates=BoundInFlux(
function(t){matrix(nrow=3,ncol=1,c(1,1,1))},
t_start,
t_end
)
mod=ParallelModel(t,k,c0,inputrates)
Y=getC(mod)
lt1=1 ;lt2=2 ;lt3=3
col1=1; col2=2; col3=3
plot(t,Y[,1],type="l",lty=lt1,col=col1,
ylab="C stocks",xlab="Time")
lines(t,Y[,2],type="l",lty=lt2,col=col2)
lines(t,Y[,3],type="l",lty=lt3,col=col3)
legend(
"topleft",
c("C in pool 1",
"C in 2",
"C in pool 3"
),
lty=c(lt1,lt2,lt3),
col=c(col1,col2,col3)
)
Y=getAccumulatedRelease(mod)
plot(t,Y[,1],type="l",lty=lt1,col=col1,ylab="C release",xlab="Time")
lines(t,Y[,2],lt2,type="l",lty=lt2,col=col2)
lines(t,Y[,3],type="l",lty=lt3,col=col3)
legend("topright",c("R1","R2","R3"),lty=c(lt1,lt2,lt3),col=c(col1,col2,col3))
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

Please suggest features or report bugs with the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.