% ------------------------------------------------------------
% Eulerian phytoplankton water column model 
% (based on Huisman et al.(2006))
% ------------------------------------------------------------
%Programmed by C Lindemann & DL Aksnes (02102018)

clear all
close all

%% Parameter

% Temporal and spatial discretization
nZ=100;               % Number of depth cells 
DZ=1;                 % Length of a depth cell (m)
nT=10;                % Number of time steps  
DT=60*10;             % Length of a time step(s)'
t= (0:nT)/(60*60*24/DT);    % for plotting

%Water column properties 
kappa=10/(60*60*24) ; % Turbulent diffusivity (m2 s-1)
I0=100;                 % Irradiance at surface (mmol photons m-2 s-1)
Ksw=0.15;               % Background light attenuation (m-1)
NB=5;               % Fixed nutrient concentration in last depth cell (mmol N m-3)     

%Phytoplankton traits 
um=1/(60*60*24);     % Max growth rate (s-1)
m=0.25/(60*60*24);   % loss (death)rate (s-1)
Hi=20;               % Half saturation irradiance (micromol photons m-2 s-1)
Hn=0.025;            % Half saturation nutrient (mmol m-3)
v=1/(60*60*24);      % Sinking velocity of phytoplankton (m s-1) 
alpha=0.5;           % Nutrient recycling coefficient

%% testing variable senstivity analysis

varKsw=[0,0,0,0,0];

Nall=zeros(length(varKsw),nZ,nT+1);
Pall=zeros(length(varKsw),nZ,nT+1);

for Ktest=1:length(varKsw)
  
    Ksw=varKsw(Ktest);
    
%% Initial conditions

% define initial conditions
phy1=linspace(0.3,0.05,nZ);
nit1=linspace(0.1,NB,nZ);

% create storage structures
P=zeros(nZ,nT+1);
N=zeros(nZ,nT+1);

% save initial conditions to storage
P(:,1)=phy1;
N(:,1)=nit1;

%% Discretisation & Constant rates

umax = um*DT;          % constant maximum growth rate
Mort=m*DT;             % constant Mortality
Diff=kappa*DT/(DZ*DZ); % the fraction of a depth cell that is exchanged 
                       % with the neighbour cells due to turbulent diffusivity
Sink=v*DT/DZ;          % The fraction of phytoplankton in a depth cell that sinks 
                       % to the cell below during DT

% time-independent light in the water column
z = 0.5*DZ:DZ:DZ*nZ-0.5*DZ; % depth array
Iz=I0*exp(-Ksw.*z);         % water column irradiance
Ilim=Iz./(Iz+Hi);           % Light limitation 
        
%% Time loop
       
for tstep=1:nT 
   
    %Phytoplankton growth 
    Nlim=nit1./(nit1+Hn); %nutrient limitation of depth cell
    u=umax*min(Ilim,Nlim);     %phytoplankton growth rate 

    %Update first (surface) depth cell - no exchange through surface
    nit2(1)=phy1(1)*(-u(1)+Mort*alpha)               + nit1(1)*(1-Diff) +nit1(2)*Diff; 
    phy2(1)=phy1(1)*(u(1)-Mort)        -phy1(1)*Sink + phy1(1)*(1-Diff) +phy1(2)*Diff;
  
    %Update last depth cell. Nutrient concentration is fixed in last cell! 
    %Phytoplankton is not transported (by diffusion/sinking) below the last depth cell.
    nit2(nZ)=NB; 
    phy2(nZ)=phy1(end)*(u(end)-Mort) +phy1(end-1)*Sink + phy1(end)*(1-Diff) +phy1(end-1)*(Diff);
   
   % Update all the intermediate depth cells
    for n=2:nZ-1
        nit2(n)= phy1(n)*(-u(n)+Mort*alpha)                                   + nit1(n)*(1-Diff-Diff) +nit1(n-1)*Diff  +nit1(n+1)*Diff;
        phy2(n)= phy1(n)*(u(n)-Mort)         + phy1(n-1)*Sink - phy1(n)*Sink  + phy1(n)*(1-Diff-Diff) +phy1(n-1)*Diff +phy1(n+1)*Diff;
    end
    
    %Store results
    N(:,tstep+1)=nit2(:); 
    P(:,tstep+1)=phy2(:);
                                   
    %Proceed to next time step with updated state variables
    nit1=nit2;                
    phy1=phy2;
end 

Nall(Ktest,:,:)=N;
Pall(Ktest,:,:)=P;

%% plot output

figure;
subplot(2,1,1)
contourf(t,z,N);
title(['Nutrient concentration (mmol m^-^3) for K=' num2str(Ksw)]);
xlabel('Time (days)');
ylabel('Depth (m)');
colorbar
set(gca,'YDir','reverse');

subplot(2,1,2)
contourf(t,z,P);
title(['Phytoplankton concentration (mmol m^-^3) for K=' num2str(Ksw)]);
xlabel('Time (days)');
ylabel('Depth (m)');
colorbar
set(gca,'YDir','reverse');


%clear N P
end

Nfinal=squeeze(Nall(:,:,end));
Pfinal=squeeze(Pall(:,:,end));