% Laboratory in Oceanography - Data and Methods
% SMAST, UMass Dartmouth
%
% Example of simple numerical model for simulating Allens pond tide 
% range behavior
% 
% Written by Miles A. Sundermeyer, 4/9/08
% Last modified 2/23/09

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assume the following 1-d momentum balance:
% du/dt + u*du/dx = -1/rho * dP/dx + Cd*|u|^2

g = 9.81;			% (m/s^2)
rho_0 = 1024;			% (kg/m^3)
vel_scalefac = 1000;		% to scale down u for basin size

Cd_all = 0.25*[0.001 0.01 0.1 1.0];			% (1/m) quadratic drag coeff

mm=2;							% set which of these Cd's to run

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Cd = Cd_all(mm);

color = ['bgkr'];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define the basin geometry
pond_length = 0.5;		% (km)
channel_length = 0.3;		% (km)
channel_depth = 2*0.3048;	% (m)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% start with a tidal signal outside the estuary
M2 = 0.0805114 * 24;		% M2 tide in cycles per day
S2 = 0.0833333 * 24;

time = [0:60:30*86400];		% (s)
time = [0:60*24:30*86400];		% (s)

A_outer = 4*0.3048;		% (m)
stage_outer = A_outer * ( sin(2*pi*time*M2/86400) ...
			  + sin(2*pi*time*S2/86400) );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(mm)
subplot(2,1,1)

plot(time/86400,stage_outer,'c:')
hold on
xlabel('time (days)')
ylabel('stage_{outer} (m)')

set(gca,'xlim',[0 max(time/86400)])
set(gca,'ylim',A_outer*2.1*[-1 1]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
stage_inner_0 = 2*0.3048;	% (m)

clear dPdx u stage_inner

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% step through and compute flow and plot
for n=1:length(time)
  if n==1
    stage_inner(n) = stage_inner_0;
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % compute the velocity and step the stage
  if (stage_inner(n) > stage_outer(n))
    if (stage_inner(n) > channel_depth)
      if (stage_outer(n) >= channel_depth)
        dPdx(n) = -rho_0 * g * (stage_inner(n) - stage_outer(n));
      elseif (stage_outer(n) < channel_depth)
        dPdx(n) = -rho_0 * g * (stage_inner(n) - channel_depth);
      end
    else
      dPdx(n) = 0;
    end
  elseif (stage_inner(n) <= stage_outer(n))
    if (stage_outer(n) > channel_depth)
      if (stage_inner(n) >= channel_depth)
        dPdx(n) = rho_0 * g * (stage_outer(n) - stage_inner(n));
      elseif (stage_inner(n) < channel_depth )
        dPdx(n) = rho_0 * g * (stage_outer(n) - channel_depth);
      end
    else
      dPdx(n) = 0;
    end
  end

  if dPdx(n) < 0	% low tide
    u(n) = -sqrt(-1/rho_0 * dPdx(n) / Cd);
  else			% high tide
    u(n) = sqrt(1/rho_0 * dPdx(n) / Cd);
  end

  stage_inner(n+1) = stage_inner(n) + u(n)/vel_scalefac;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  figure(mm)

  xlim = [0 max(time/86400)];
  xlim = [time(n)/86400 + 3.5*[-2 2]];

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  subplot(4,1,1)
  cla
  plot(time/86400,stage_outer,'c-')
  hold on
  plot(time(1:n)/86400,stage_outer(1:n),'b-')

  ylabel('stage_{outer} (m)')

  set(gca,'xlim',xlim)
  set(gca,'ylim',A_outer*2.1*[-1 1]);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  subplot(4,1,2)
  plot(time(1:n)/86400,stage_inner(1:n),[color(mm),'-'])

  ylabel('stage_{inner} (m)')

  set(gca,'xlim',xlim)
  set(gca,'ylim',[0.9*channel_depth A_outer*2.05]);
  %set(gca,'ylim',A_outer*2.1*[-1 1]);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  subplot(4,1,3)
  plot(time(1:n)/86400,u(1:n)/100,[color(mm),'-'])

  xlabel('time (days)')
  ylabel('channel vel (m/s)')

  set(gca,'xlim',xlim)
  set(gca,'ylim',1.7*A_outer*[-1 1]);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % draw bottom bottom
  subplot(4,1,4)
  cla

  fill([0 0 pond_length*[1 1] (pond_length + channel_length)*[1 1] 1 1], ...
		[-A_outer*4.0 0 0 channel_depth*[1 1] -A_outer*2.1*[1 1 4/2.1]],'k');
  hold on

  % draw the channel bottom 
  fill([pond_length*[1 1] (pond_length + channel_length)*[1 1]], ...
		channel_depth*[0 1 1 0],'k');
  
  % draw the estuary height
  h1 = fill([0 0 pond_length*[1 1]], stage_inner(n)*[0 1 1 0],'b');

  % draw the ocean height
  h2 = fill([(pond_length + channel_length)*[1 1] 1.0 1.0],  ...
		[-4 stage_outer(n)*[1 1] -4],'b');

  % draw the channel water height
  if (stage_outer(n) > channel_depth)
    h3 = fill([pond_length*[1 1] (pond_length + channel_length)*[1 1]], ...
	[channel_depth stage_inner(n) stage_outer(n) channel_depth],'b');
  elseif (stage_outer(n) < channel_depth)
    h3 = fill([pond_length*[1 1] (pond_length + channel_length)*[1 1]], ...
	[channel_depth stage_inner(n) channel_depth channel_depth],'b');
  end

  set(h1,'edgecolor','none');
  set(h2,'edgecolor','none');
  set(h3,'edgecolor','none');

  xlabel('Along-estuary dist (km)')
  set(gca,'xlim',[0 1]);
  set(gca,'ylim',A_outer*2.1*[-1 1]);

  drawnow

end