function Var3
% 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. 4a, top left 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] %note: this value is actually given for N not C in Moore's paper resulting in a slight bias
c2=0.44; % grazing half saturation [mmole^2 m^-6] %note: this value is actually given for N not C in Moore's paper resulting in a slight bias
c3=0.5;  % ingestion efficiency [unitles]
c4=0.06; % herbivore non-predation loss [day^-1]
c5=1.0; %herbivore quadratic mortality [mmole^-1 m^2 day^-1] %note: this value is actually given for N not C in Moore's paper resulting in a slight bias
N0=10;  %deep nutrients [mmole m^-3]
J=0.5;  %uptake half saturation [mmole m^-3]
b=0.07; %phytoplankton respiration [day^-1]
lat=47; %latitude [degrees]
m=3; %diffusion velocity [m day^-1]

tspan=[0 365*9:1:365*10-1];%integrate for 10 years get data, for year 10.
y0=[10 0.26 0.26]; %initial condition

%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(:,2)=y(:,2)*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(:,3)=y(:,3)*14*6;

figure(1)

II=find(365*9+1<tspan & tspan < 365*10); %plot year 10
[AX,H1,H2] = plotyy(tspan(II)-min(tspan(II)),y(II,2),tspan(II)-min(tspan(II)),y(II,3),'plot');
set(get(AX(1),'Ylabel'),'String','P [mg m^{-3}] of Chl')
set(AX(1),'YLim',[0 2.5])
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.5:2.3])
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 VariationNo3.txt out -ascii
%-------------------------------------------------------

 function dydt = f(t,y)
  % Derivative function
  year_day=mod(t,365);
  MLD=75+50*cos(2*pi*(year_day-60)/365);  %mixed layer depth. Affects only the light level.
   w_e=-50*2*pi/365*sin(2*pi*(year_day-60)/365); %entrainment velocity. zero when MLD shallows
   w_e_pos=w_e;
   if w_e_pos<0 
       w_e_pos=0 
   end
 
   
   dydt = [ 
      -(alpha(lat,year_day,MLD)*y(1)/(J+y(1))-b)*y(2)+(m+w_e_pos)*(N0-y(1))/MLD;
      (alpha(lat,year_day,MLD)*y(1)/(J+y(1))-b)*y(2)-c1*y(3)*((y(2)*y(2))/((y(2)*y(2))+c2));%-(m+w_e_pos)*y(2)/MLD;
      c3*c1*y(3)*(y(2)*y(2))/((y(2)*y(2))+c2)-c4*y(3)-c5*(y(3)*y(3))];%-w_e*y(3)/MLD];
 end

%-------------------------------------------------------
  

function growth_rate=alpha(lat, year_day, MLD)
% function to compute MLD averaged growth rate
% inputs: lat [degrees], year_day [Julian day], MLD [m]
% based on Evans and Parslow, 1985.

cloud_cover=0.9;
trans_loss=3/8;
solar_const=1.94; %[ly]
alpha_J=0.04*1343/1.94; %Low light photosynthetic slope [ly^-1]
Q=2; %maximum photosynthetic rate [d^-1]
kk=0.1; %light attenuation by water [m^-1]

delta=-0.406*cos(2*pi*(year_day+10)/365); % solar declination
tau=acos(-tan(delta)*tan(lat*pi/180)); %half day length in radians
tau_days=tau/pi/2; %half day length in days
Daily_radiation_TOA=(sin(delta)*sin(lat*pi/180).*tau+cos(delta)*cos(lat*pi/180).*sin(tau))*solar_const/pi*trans_loss*.7*cloud_cover; %Jxtau
beta=Q*tau_days/alpha_J/(Daily_radiation_TOA/tau_days); %has units of days
yy=beta.*exp(kk*MLD); 
F1=sqrt(yy.^2+tau_days.^2)-tau_days.*log((tau_days+sqrt(yy.^2+tau_days.^2))./yy);
F2=sqrt(beta.^2+tau_days.^2)-tau_days.*log((tau_days+sqrt(beta.^2+tau_days.^2))./beta);
F3=yy;
F4=beta;
growth_rate=2*Q./MLD/kk.*[F1-F2-F3+F4];
end
%-------------------------------------------------------

end
