Hi,
I am trying to run a 1D nutrient-phytoplankton-zooplankton model in R
using the package 'deSolve'. The code is shown below:
DEPTH = seq(2.5, 147.5, 5)
NPZ = function(t, state, params){
with(as.list(params), {
P <- state[1:NB]
Z <- state[(NB + 1): (2*NB)]
N <- state[(2*NB + 1): (3*NB)]
F.I = function(z, hr){
I0 = function(hr){
if (hr <=16){
450*sin(hr*pi/16)
} else {
0
}
}
i = which(z == DEPTH)
I = I0(hr=hr)*0.43*exp(-z*Kw - Kc*(cumsum(delz*P)[i]
- delz*P[i]/2))
FI = alpha*I/(um^2 + (alpha*I)^2)^.5
return(FI)
}
J1 = function(x){
n = 10^3
FI2 <- function(y)F.I(z = x, y)
J1 = sum(sapply(runif(n, 0, 16), FI2))*16/24/n
return(J1)
}
J = sapply(DEPTH, J1) #the factor of light limitation
deltaz <- c(rep(5, NB - 1), 2.5)
P.Flux <- -D*diff(c(P,0))/deltaz
Z.Flux <- -D*diff(c(Z,0))/deltaz
N.Flux <- -D*diff(c(N,N0))/deltaz
dP.dt = P*um*N/(N+Kn)*J - Z*Im*P^2/(P^2+Kp^2) - diff(c(0,P.Flux))/delz
dZ.dt = Z*(eps*Im*P^2/(P^2+Kp^2) - g*Z) - diff(c(0,Z.Flux))/delz
dN.dt = -P*um*N/(N+Kn)*J + Z*(1-eps)*Im*P^2/(P^2+Kp^2) + g*Z^2
- diff(c(0,N.Flux))/delz
return(list(c(dP.dt, dZ.dt, dN.dt)))
})
}
params1 = c(um = 1, #maximal phytoplankton growth rate, d-1
Kn = 0.2, #nitrogen half-saturation constant uM/L
Kw = 0.04, #light attenuation due to water, unit: m^-1
Kc = 0.03, #light attenuation by phytoplankton,
#unit: m^2 (mmol N)^-1
alpha = 0.025,
Im = 0.6, #zooplankton maximal ingestion rate, unit: /d
Kp = 1, #zooplankton half-saturation for
ingestion, unit: mmol N/m^3
eps = .3, #Growth efficiency of zooplankton
delz = 5,
D = 2.6, #diffusion coefficiency, unit: m^2/d
g = 0.025, #zooplankton mortality rate, unit: d^-1
(mmol N)^-1
N0 = 14, #Nitrate concentration at 150 m
NB = length(DEPTH)) #Number of depth intervals
P.int = c(rep(0.07, 8),.12, .18,rep(0.21,3),rep(0.35,2),
rep(0.42,2),rep(0.07, 8), rep(0.01, 5))
Z.int = P.int/7
N.int = c(rep(0.1, 13), rep(4.55, 5), rep(10,12))
state = c(P.int, Z.int, N.int)
Time = seq(0,100, by = 1)
NPZ.out <- ode.1D(state, Time, NPZ, params1, nspec = 3, names =
c("P","Z","N"))
After running for about 7 hours, the result returns:
DLSODA- At current T (=R1), MXSTEP (=I1) steps
taken on this call before reaching TOUT
In above message, I [1] 5000
In above message, R [1] 0.0498949
Warning?
1: In lsoda(y[ii], times, func = bmod, parms = parms, bandup = nspec * :
an excessive amount of work (> maxsteps ) was done, but integration
was not successful - increase maxsteps
2: In lsoda(y[ii], times, func = bmod, parms = parms, bandup = nspec * :
Returning early. Results are accurate, as far as they go
And the diagnostics(NPZ.out) shows:
lsoda return code
--------------------
return code (idid) = -1
Excess work done on this call. (Perhaps wrong Jacobian type MF.)
--------------------
INTEGER values
--------------------
1 The return code : -1
2 The number of steps taken for the problem so far: 5000
3 The number of function evaluations for the problem so far: 44438
5 The method order last used (successfully): 1
6 The order of the method to be attempted on the next step: 1
7 If return flag =-4,-5: the largest component in error vector 0
8 The length of the real work array actually required: 1732
9 The length of the integer work array actually required: 110
14 The number of Jacobian evaluations and LU decompositions so far: 4834
15 The method indicator for the last succesful step,
1=adams (nonstiff), 2= bdf (stiff): 2
16 The current method indicator to be attempted on the next step,
1=adams (nonstiff), 2= bdf (stiff): 2
--------------------
RSTATE values
--------------------
1 The step size in t last used (successfully): 2.189586e-06
2 The step size to be attempted on the next step: 1.039085e-05
3 The current value of the independent variable which the solver has
reached: 0.0498949
4 Tolerance scale factor > 1.0 computed when requesting too much accuracy:
0
5 The value of t at the time of the last method switch, if any: 0.005314453
So I wonder what are the exact problems with this model. Is it just
because of inappropriate parameters? I also wonder how I can speed up
the calculation time in R.
Thanks,
--
Bingzhang Chen
Ph. D.,
State Key Lab of Marine Environmental Science,
College of Oceanography and Environmental Science,
Xiamen University,
Xiamen, Fujian 361005
P. R. China
