Nothing
### ============================================================================
### Demonstration of a 2box simulation
### Here we use the data set of Bautzen reservoir with
### flexible mixing depth (ZMIX)
### ============================================================================
library(rSALMO)
### ================ work in progress ======================
### Warning: transport, sedimentation and sediment not yet implemented
### ToDo:
### - add functions or data set to handle hypsographic function
### - create 2 cbind() forcing matrices for epi- and hypolimnion
### - implement details of function SALMO.2box
### Open Tasks
### - sediment area
### - depth
### - oxygen
### ========================================================
## Data set from workgroup limnology of TU Dresden
data(bautzen1997)
## Hypsographic table of Bautzen reservoir
data(bautzen_hypso)
hyps <- hypso_functions(bautzen_hypso)
#bautzen1997$s
## sediment area of hypo- and epilimnion
## upper boundaries of hypo- and epilimnion
levels <- with(bautzen1997, cbind(hypo = s - zmixreal, epi = s))
## calculate sediment area row-wise for pairs of hypo and epilimnion depths
## "1" means row wise, t() means that result needs to be transposed
## result: 1st colum: hypolimnic sediment area, 2nd: epilimnic sediment area
# ToDo: make this a function
sediment_area <- t(apply(levels, 1, hyps$sediment_area))
## same for pelagic ratio (i.e. ratio of water contact area to total area, rest is sediment)
pelagic_ratio <- t(apply(levels, 1, hyps$pelagic_ratio))
## check pelagic_ratio calculation
# areaE <- hyps$area(levels[,2])
# areaH <- hyps$area(levels[,1])
#
# sedH <- areaH
# sedE <- areaE - sedH
#
# ratioH <- 1 - sedH/areaH
# ratioE <- 1 - sedE/areaE
#
# pr_test <- cbind(ratioH, ratioE)
## Reformat old-style input data structure into new structure
forcE <- with(bautzen1997,
data.frame(
time = t, # simulation time (in days)
vol = ve, # volume (m^3)
depth = zmixreal, # absolute depth of the layer (m), required for resuspension depth
dz = zmix, # zmix, or layer thickness (m)
qin = qin, # water inflow (m^3 d^-1)
ased = sediment_area[,2], # sediment contact area of the layer (m^2); important
srf = srf, # strong rain factor, an empirical index of turbidity
iin = iin, # photosynthetic active radiation (J cm^2 d^-1); approx 50% of global irradiation
temp = temp, # water temperature (deg. C)
nin = nin, # DIN concentration in inflow (mg L^-1)
pin = pin, # DIP concentratioin in inflow (mug L^-1)
pomin = pomin, # particulate organic matter in inflow, wet weight (mg L^-1)
zin = zin, # zooplankton in inflow (w.w. mg L^-1)
oin = oin, # oxygen concentration in inflow (mg L^-1)
aver = pelagic_ratio[,2], # 1 - (ratio of sediment contact area to total area); (redundant if SF=0)
ad = ae, # downwards flux between layers (m^3 d^-1)
au = ah, # upwards flux between layers (m^3 d^-1)
diff = 0, # eddy diffusion coefficient
x1in = xin1, # phytoplankton import of group 1 (w.w. mg L^-1)
x2in = xin2, # phytoplankton import of group 2 (w.w. mg L^-1)
x3in = 0 # phytoplankton import of group 3 (w.w. mg L^-1)
)
)
forcH <- with(bautzen1997,
data.frame(
time = t,
vol = vh,
depth = s - 154,
dz = zhm,
qin = qhin,
ased = sediment_area[,1],
srf = srf,
iin = 0,
temp = temph,
nin = nhin,
pin = phin,
pomin = pomhin,
zin = zhin,
oin = ohin,
aver = 0, # is zero in hypolimnion
ad = ae, # check this !!
au = ah, # check this !!
diff = 0,
x1in = xhin1,
x2in = xhin2,
x3in = 0
)
)
## Matrices are faster than data frames
forc <- as.matrix(cbind(forcE, forcH))
## Read default parameter set. The order of the vector must not be changed!
data(pp) # Phytoplankton parameters of SALMO, matrix with 1 column per phytopl. species
data(cc) # other SALMO parameters, vector
## A few parameters that are specific for Bautzen Reservoir
cc[c("MOMIN", "MOT", "KANSF", "NDSMAX", "NDSSTART", "NDSEND", "KNDS", "KNDST")] <-
c(0.005, 0.002, 0, 0.095, 0, 365, 0.005, 1.03)
#cc["Zres"] <- 0
cc["SF"] <- 0 # switch internal sedimentation off
## N fixer on
pp["NFIX",1] <- 0.24
## Background light extinction is lake specific
cc["EPSMIN"] <- 0.7
## Check that internal sedimentation is switched off
cc["SF"] == 0
## Set of technical parameters, a.o. for dimensioning dynamic variables
nOfVar <- c(
numberOfInputs = 21,
numberOfOutputs = 14,
numberOfStates = 11,
numberOfParameters = NA,
numberOfAlgae = 3,
numberOfLayers = 2
)
## Total number of parameters passed to the model (very important)
nOfVar["numberOfParameters"] <- length(pp) / nOfVar["numberOfAlgae"]
## Put all 3 kinds of model parameters in a list
parms <- list(pp=pp, cc=cc, nOfVar=nOfVar)
## -----------------------------------------------------------------------------
## ice in Dec, Jan, Feb
ice <- ifelse(((forcE$time + 30) %% 365.24) < 90, 1, 0)
forcE$iin <- with(forcE, ifelse(ice, iin * cc["RLW"], iin * cc["RL"]))
## Initial values
## Xi = Phytoplankton biomass (mg/L w.w.)
## Z = Zooplankton Biomass (mg/L w.w.)
## D = allochthonous detritus (mg/L w.w.)
## O = Oxygen concetration (mg/L)
## Gi = Carbon:Chlorophyll ratio for Baumert's photosynthesis model
## (currently not used)
x0 <- c(N=4, P=80, X1=.3, X2=.3, X3=.3, Z=1, D=4, O=14, G1=0, G2=0, G3=0)
## 2 boxes -----------------
x0 <- rep(x0, 2)
## Call one time for testing
ret <- call_salmodll("SalmoCore", nOfVar, cc, pp, forc, x0)
## Simulation time steps
times <- seq(0, 365, 0.1)
#times <- seq(0, 200, 1)
## Model simulation with "lsoda" from package deSolve
#system.time(
out <- ode(x0, times, salmo_2box, parms = parms, method="rk4", inputs=forc)
#)
#plot(out, which=c("X1", "X2", "X3", "N", "P", "O", "iin"), mfrow=c(4,3))
plot(out, which=c("X1", "X2", "X3", "N", "P", "O"), mfrow=c(4,3))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.