%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simple integration of the advection-diffusion equation.%
% Includes correction to Kz for convective instabilities.
% Coastal and Shelf Oceanography
% February 2003
% Jonathan Sharples
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear			% Matlab command to clear all variables
close all		% Matlab command to make sure all graphics windows are closed

% set up the initial conditions
N=10;			% Number of elements in the vertical
ds_dx=1.0e-3;		% Horizontal salinity gradient
dt=10;			% time step (seconds)
dz=2.0;			% depth step (metres)
s_old(1:N)=33.0;	% Initial value for salinity through all of the depth
time=0;			% Start time = 0.0 seconds
Kz1=0.01;		% vertical eddy diffusivity (m2 s-1)
Kz2=0.15;		% vertical eddy diffusivity (m2 s-1) for convective instabilities
itotal=fix((12.42*3600)/dt);	% total number of timesteps in 1 full tidal cycle (fixed to nearest integer)

height=N*dz;		% total height of water column is N x depth cell thickness
for i=1:N,		% These lines calculate the depths of each of the depth cells
  depths(i)=-((N*dz)-(i*dz-dz/2));
end,

u0=1;			% Horizontal current amplitude (metres per sec) for use in tidal current equation
omega=1.40525e-4;	% tidal frequency (M2 tide, radians per second)

% The next lines provide a loop that works through the integration of the advection equation.
% At each step the current profile is calculated, the effect of advection on the salinity is calculated,
% the old salinity is updated with the new salinity, and the time is incremented by the time step.

for i=1:itotal,		
  
  % Calculate the depth profile of the horizontal current
  for k=1:N,
    b=depths(k)/height;
    u(k)=u0*sin(omega*time)*(1.15-0.425*b^2);	% A simple quadratic representation of the depth-variability of a current
  end,
  
  % Then advect the water horizontally at each position through the depth
  for k=1:N,
    delta_s(k)=-u(k)*ds_dx*dt;		% Calculate the change in salinity over the time step
    s_new(k)=s_old(k)+delta_s(k);	% Calculate the new salinity, knowing the old salinity and the salinity change
  end,
  s_old=s_new;				% update the salinity
  
  % A solution to the convective instability of the water:
  % Increase Kz if the density increases with distance above the bed
    Kz=Kz1;
    if s_old(1) < s_old(N),
      Kz=Kz2;
    end,

  % And mix the salinity vertically using the vertical eddy diffusivity
  
  % First do the interior points where d2s/dz2 is well-defined
  for k=2:N-1
    deltas(k)=dt*(Kz*((s_old(k+1)-s_old(k))-(s_old(k)-s_old(k-1)))/dz^2);
    s_new(k)=s_old(k)+deltas(k);
  end,

  % Do the bottom depth cell with the condition of no s flux through the bed
  deltas(1)=dt*(Kz*(s_old(2)-s_old(1))/dz^2);
  s_new(1)=s_old(1)+deltas(1);

  % Do the surface depth cell with the condition of no s flux through the surface
  deltas(N)=dt*(Kz*(s_old(N-1)-s_old(N))/dz^2);
  s_new(N)=s_old(N)+deltas(N);
  
  time=time+dt;			% Update the time by one time step
  s_old=s_new;			% Update the old salinity with the new salinity
  
  itest=MOD(i,100);
  if itest == 0,						% We only want to plot the data every 100 time steps
    plot(s_old,depths)		% Plot the data as a salinity profile
    xlabel('salinity (PSS78)')	% Add a label to x axis
    ylabel('depth (metres)')    % Add a label to y axis
    smin=mean(s_old)-0.75;	% set a minimum value for the salinity axis
    smax=mean(s_old)+0.75;	% set a maximum value for the salinity axis
    AXIS([smin smax -height 0])	% control the axes limits
    pause(0.25)			% wait for 0.25 sec before continuing
  end,
  
end,

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5