Names of the state variables, parameters and their standard values are shown in Table
4.1.
Box 4.3 Matlab script for an NPZ model similar to
that of Franks et al. (1986)
(^)
%Program Franks
clear all
ndays=60; nsteps=120;
%Set up storage vectors for results
Pct=ones(nstepsndays,1); Tct=Pct; Nct=Pct; Zct=Pct;
%Set values of all parameters
%for standard run. Daily rates
%are reduced for small time steps
Vm=2./nsteps; m=0.1/nsteps;
Rm=1.5/nsteps; d=0.2/nsteps;
Mix=0.02/nsteps;
%Try Vm=0.69; try other changes
Gamma=0.3; Ks=1.; Lmda=1.; NatZ=10.6;
%Set values of starting conditions
NIT=1.6; P=0.3; Z=0.1; P0=0.;
%P0 initially set to zero; try 1.0
ct=0;
%Main loop starts here,
%one cycle per model day
for i=1:ndays
%To add autumn mixing use
%the following statements;
%make ndays=200.
%if i>120
%Mix=Mix+0.02/nsteps; end
%Subloop to allow nstep
%time steps per day
for j=1:nsteps
ct=ct+1; Tct(ct)=ct/nsteps;
Nct(ct)=NIT; Pct(ct)=P; Zct(ct)=Z;
UPTAKE=VmNIT/(Ks+NIT);
if P>P0
Ivlev=Rm(1-exp(-Lmda(P-P0)));
else Ivlev=0.;
end
delP=UPTAKEP-mP-ZIvlev;
delZ=GammaZIvlev-dZ;
delN=-UPTAKEP+mP+...
(1-Gamma)ZIvlev+dZ;
%To mix with deeper water and sink
%some organic matter use the
%following delN= statement instead:
%delN=-UPTAKEP+0.6mP+0.4
% (1-Gamma)ZIvlev+0.4dZ
% +Mix(NatZ-NIT);
%Calculate new values of P, Z, NIT:
P=P+delP; Z=Z+delZ;
NIT=NIT+delN;
end
end
figure
plot(Tct(1:end), Nct(1:end),′k′);
hold on;
plot(Tct(1:end), Pct(1:end),′g′);
plot(Tct(1:end), Zct(1:end),′r′);
(^)