% vertical_pred_prey is a Matlab script for calculating
% the interactions among particle flux, flux feeder, and prey
%
%
% Time-stamp: <2004-05-17 16:53:03 adrian>
%
%
% USAGE:
% Type vertical_pred_prey at the Matlab prompt.
%
% HISTORY:
% **-**-01: Original code by GAJ.
% 25-04-03: Code re-written & re-organized by ABB
%
% NOTES:
% For theory, see
% Jackson, G.A. & A.B. Burd (2002) Deep-Sea Res. II 49:193-218
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% George A. Jackson
% Department of Oceanography, Texas A&M University.
%
% Adrian B. Burd
% Department of Marine Sciences, University of Georgia
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Clear workspace
close all
clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% MODEL PARAMETERS
%
% Define dimensional (and some dimensionless) parameters and from
% them, calculate required dimensionless parameter values.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define dimensional parameters
a = 0.3; % efficiency of P eating particles
b = 123; % clearance rate for P [m^3/gC_p/d]
c = 0.3; % efficiency of Q eating P
d = 0.01; % Predator specific mortality [1/d]
F_surf = 0.1; % Surface particle flux [gC/m^2/d]
sigma_p = 3.95; % Interaction cross section [m^2/gC_p]
lambda = 6.25e-04; % Bacterial degradation rate [1/d]
v_f = 30; % Particle settling vely [m/d]
omega_f = 2*pi/365; % Flux forcing frequency [1/d]
% Define some non-dimensional parameters
gamma = 0.25; % relative amplitude of seasonal flux oscillations
% Calculate some dimensional constants
beta_0 = sqrt(a*sigma_p*d*F_surf); % Lotka Voleterra frequency
p_ss = d/(b*c); % Steady-state P [gC/m^3]
% Calculate non-dimensional parameters
psi = sigma_p*v_f/(b*c);
nu = b*c*lambda/(d*v_f*sigma_p);
beta_0 = beta_0/d;
omega_f = omega_f/d;
beta_0 = omega_f/0.5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% MODEL GRID
%
% Set up a grid from 100 -> 1000 m depth with 100 layers and 0 -> 4
% years in time with the time step determined by the CFL
% condition. Note convention that dimensional depths and time have
% "depth" and "t" in their names, dimensionless version have "z"
% and "tau".
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t_start = 0; % Starting day [d]
t_final = 12*365; % End day [d]
top_depth = 100; % Depth of top layer [m]
bottom_depth = 1000; % Depth of bottom layer [m]
n_layers = 100; % Number of layers
tau_start = t_start*d; % Dimensionless start time
tau_final = t_final*d; % Dimensionless final time
z_top = top_depth/(v_f/d); % Dimensionless top layer
z_bot = bottom_depth/(v_f/d); % Dimensionless bottom layer
z = linspace(z_top, z_bot, n_layers); % Dimensionless depths
z = z';
delta_z = (z_bot - z_top)/(n_layers-1); % Dimensionless delta z
delta_tau = delta_z;
n_tau = floor((tau_final - tau_start)/delta_tau);
tau = (0:(n_tau-1))*delta_tau;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% THE INTEGRATION
%
% Following G, use a method of lines by discretizing space (use
% backward difference for flux and central difference for
% diffusion) and an Euler scheme for the time derivative.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some functions we need - these are the time variations of the
% particle flux at the surface and the time variations of the
% vertical migration: NOTE, the definition of ekz is different from
% that in G's code, my version is 1/his.
f0 = 1 + gamma*sin(omega_f*tau); % Flux forcing at surface
ekz = exp(-psi*(1 + nu)*z); % Dimensionless exp(-k*z)
beta_sq = beta_0^2.*ekz;
% Define a matrix of coefficients for the central differences used
% in the second order derivatives.
ckz = 1e-4;
Akz = -2*diag(ones(n_layers,1)) + diag(ones((n_layers-1),1),1) + ...
diag(ones(n_layers-1,1),-1);
Akz = sparse((delta_z^(-2)*ckz)*Akz);
% The dimensionless initial values - i.e. the initial depth
% distributions. Here we choose the steady state sans migration
% case for f and q and a Gaussian for p. This translates into a
% constant for the dimensionless f
f_init = ones(n_layers,1);
q_init = f_init.*ekz;
p_init = f_init;
% Set up indexing vectors for the finite difference
z_range_1 = 1:n_layers-1;
z_range_2 = 2:n_layers;
% Set up storage for the solutions and put the initial conditions
% in the storage
ft = zeros(n_layers, n_tau);
pt = ft;
qt = ft;
ft(:,1) = f_init;
pt(:,1) = p_init;
qt(:,1) = q_init;
% Now loop through each time step
for it = 2 : n_tau
dfdt = zeros(n_layers,1);
dpdt = dfdt;
dqdt = dfdt;
% The equation for the derivative of the flux
dfdt = psi*ft(:,it-1).*(1 - pt(:,it-1))*delta_tau;
% The prey equation. This includes a diffusive term. Note the
% difference in the formulation of the (f-q) term when we don't
% have the diffusion in.
dpdt = beta_sq.*pt(:,it-1).*(ft(:,it-1) - qt(:,it-1)./ekz) ...
+ Akz*pt(:,it-1);
% The predator equations
dqdt = qt(:,it-1).*(pt(:,it-1) - 1) + Akz*qt(:,it-1);
% Update the flux, prey and predator depth profiles
ft(1,it) = f0(it);
ft(z_range_2,it) = ft(z_range_1,it-1) + dfdt(z_range_1);
pt(:,it) = pt(:,it-1) + delta_tau*dpdt;
qt(:,it) = qt(:,it-1) + delta_tau*dqdt;
% Impose no-flux conditions at the boundaries
pt([1, n_layers],it) = pt([2, n_layers-1],it);
qt([1, n_layers],it) = qt([2, n_layers-1],it);
end
% Plot some things
figure(1)
subplot(3,1,1)
plot(tau/d/365, ft(1,:), 'b', tau/d/365, ft(n_layers,:), 'r--')
ylabel('particle flux')
title(['\psi = ' num2str(psi) ' \nu = ', num2str(nu), ' \gamma =' ...
' ' num2str(gamma)])
subplot(3,1,2)
plot(tau/d/365, pt(1,:), 'b', tau/d/365, pt(n_layers,:), 'r--')
ylabel('particle feeder p')
subplot(3,1,3)
plot(tau/d/365, qt(1,:), 'b', tau/d/365, qt(n_layers,:), 'r--')
ylabel('predator q');
xlabel('time (yr)')
figure(2)
plot(pt(1,:), qt(1,:), 'b')
plot(pt(n_layers,:), qt(n_layers,:), 'r')
xlabel('particle feeder p')
ylabel('predator q')
title(['\psi = ' num2str(psi) ' \nu = ', num2str(nu), ' \gamma =' ...
' ' num2str(gamma)])
%------------ Last line of vertical_pred_prey.m