One Phytoplankton One Resource Model

% basic_model.m
% one resource, one phytoplankton    

clear all; close all;

% initialize figure plotting...
figure
set(gcf,'Renderer','zbuffer') 


% timestep (days)
deltat = 0.1;
% number of days to timestep
max_days = 100.0 ;
nstepmax = max_days/deltat;

% nout = how often to plot data point (integer days)
nout = 0.01/deltat ;

% set phyto characteristics
% growth (day-1)
  mu1 = 0.50;
% diluation rate (day-1)
  dil1 =  0.2;
% half-saturation, K (micromol l-1)
  k1 =  0.1;

% incoming resource concentration (micromol l-1)
  nutin1 = 1.0 ;

% initialization 
% resources (micromol l-1)
  nut1 = nutin1;
% phytoplankton (micromol l-1)
  phy1 = 0.002;
% initial time = 0.0
tdays = 0.0;

% begin timestepping.....................
for nstep = 1:nstepmax,
  tdays = tdays + deltat;

% determine tendencies

% net nutrient supply rate (micromol l-1 day-1) 
  sn1 = dil1*(nutin1 - nut1) ;

% phytoplankton 
  grow1 = mu1*nut1/(nut1+k1)*phy1;
  loss1 = dil1*phy1;
  dphy1dt = grow1 - loss1;

% resources
%  dnut1dt = - (grow1+grow2) + sn1;
   dnut1dt = - grow1 + sn1;

% forward step
  phy1 = phy1 + dphy1dt*deltat;
  nut1 = nut1 + dnut1dt*deltat;

% if phytoplankton concentration gets very small, hold constant
  if (phy1 < 1.0e-30) ,  phy1 = 1.0e-30; end;


% update variables for plotting
% only plot at interval defined by nstep
if mod(nstep, nout) == 0
% phytoplankton
  phystore1(nstep) = phy1;
% resources 
  nutstore1(nstep) = nut1;
% time
  tstore(nstep) = tdays;
end

end
% end timestepping ......................


% plot results 
    figure(gcf);
    plot(tstore,phystore1,'.b-');
    set(gca,'FontSize',14) ;
    axis([0 tdays -0.5 1.5]);
    xlabel('time (days) ','FontSize',16);
    ylabel('concentration  (umol l^-^1) ','FontSize',16);
    grid;
    hold on;
    plot(tstore,nutstore1,'.r-');
    hleg1 = legend('phyto','nutrient');