
function Var2
% Code used in Bherenfeld and Boss, 2014, 
% Resurrecting the Ecological Underpinnings of Ocean Plankton Blooms? 
% Annual Review of Marine Science. 
% This code provided the results shown in Fig. 4, top right panel
% 
% In this run the code is set with a prescribed phytoplankton growth-rate
% (alpha) which is varried every month starting in day 1000.
%
% The model is based on that of Evans and Parslow, 1985, with some
% modifications following Moore et al., 2002.


% clear the memory and close all figures
clear all
close all

%parameter used:
c1=3.24; % maximum herbivore growth rate  [day^-1]
c2=0.44; % grazing half saturation [mmole^2 m^-6]
c3=0.5;  % ingestion efficiency [unitles]
c4=0.06; % herbivore non-predation loss [day^-1]
c5=1.0;  % herbivore quadratic mortality[mmole-1 m3 day-1]

tspan=[0:1:1460];%integrate for 4 years. Use only 12mo where alpha was varying
y0=[0.26 0.26]; %initial condition for Phytoplantkon and Zooplankton biomass

%solve equations numerically by integrating ODE
[t,y]=ode23(@f,tspan,y0);

%convert phytoplankton from [mM m^{-3}] of Nitrogen to Chl
%Chl (mg/m3) = Phytoplankton N (mM/m3) * 14 mg/mM * 6 C/N * 1/60 Chl/C
y(:,1)=y(:,1)*14*6/60;

%convert herbivores from [mM m^{-3}] of Nitrogen to mg/m^3 of carbon
%grazer biomass (mg/m3) = grazer N (mM/m3) * 14 mg/mM * 6 C/N
y(:,2)=y(:,2)*14*6;

figure(1)
II=find(tspan>970 & tspan < 1330);
[AX,H1,H2] = plotyy(tspan(II)-min(tspan(II)),y(II,1),tspan(II)-min(tspan(II)),y(II,2),'plot');%red=phyto, blue=grazers
set(get(AX(1),'Ylabel'),'String','P [mg m^{-3}] of Chl')
set(AX(1),'YLim',[0 0.6])
set(AX(1),'XLim',[0 360])
set(get(AX(2),'Ylabel'),'String','H [mg m^{-3}] of C') 
set(AX(2),'Ylim',[0 60])
set(AX(2),'XLim',[0 360])
set(AX(1),'XTick',[0:30:360],...
        'YTick',[0:0.1:0.7])
set(AX(2),'XTick',[0:30:360],...
        'YTick',[0:10:60])
    
xlabel('year day')
legend('Phytoplankton','Herbivores')

%%% Prepare the data for out the first concatenates the 
%%% "t" and "y" variabiles so they are a single variable (3 columns wide)
%%% the second line save the data out, the syntax is: save filename variable_name -ascii
out=[t,y];
save VariationNo2.txt out -ascii
%-------------------------------------------------------

 function dydt = f(t,y)
  % Derivative function
  % alpha is the phytoplankton growth rate [day^-1]
  if t<1000 || t>1270
      alpha=0.05; 
  elseif (t>=1000 && t<1030) || (t>=1240 && t<1270)
      alpha=0.15;
  elseif (t>=1030 && t<1060) || (t>=1210 && t<1240)
      alpha=0.25;
  elseif (t>=1060 && t<1090) || (t>=1180 && t<1210)
      alpha=0.35;
  elseif (t>=1090 && t<1120) || (t>=1150 && t<1180)
      alpha=0.45;
  else
      alpha=0.55;
  end

  dydt = [ alpha*y(1)-c1*y(2)*((y(1)*y(1))/((y(1)*y(1))+c2));
           c1*c3*y(2)*((y(1)*y(1))/((y(1)*y(1))+c2))-(c4*y(2))-(c5*y(2)^2)];
 end
 
%-------------------------------------------------------
end

