%  bore;  1-D gravitational adjustment of a 
%  single shallow layer of constant density fluid.
%  Useful for showing the formation of a bore.
%  Uses an A-grid (velocity and height at the
%  same grid points) and upwind differencing on advection
%  terms to preserve a shock. Time-stepping is by the  
%  leap-frog/trapezoidal time method.  This may be about the 
%  simplest PE model that works reasonably well for non-linear
%  cases, i.e., you can recognize the major features of the 
%  solution, including a fairly reasonable shock (or bore). 
%  This is related to the geostrophic adjustment problem 
%  but is plotted with different scales, etc.  
%
%  Jim Price,  January, 2001.

clear
close all
format compact

set(0,'DefaultLineLineWidth',1.5)
set(0,'DefaultTextFontSize',12)
set(0,'DefaultAxesLineWidth',1.4)
set(0,'DefaultAxesFontSize',12)

%  display the text string str2
str2(1) =  {'bore:' };
str2(2) =  {'  '};
str2(3) =  {'Adjustment under gravity in one dimension.  '};
str2(4) =  {'Solved by a primitive equation numerical '};
str2(5) =  {'model representing one (shallow) layer.'};
str2(6) =  {'The initial thickness represents a double '}; 
str2(7) =  {'front having zero initial current.  '}; 
str2(8) =  {' '}; 
str2(9) =  {' '};
str2(10) = {'This code is in the public domain.'};
str2(11) = {'Hit any key to continue.'};
str2(12) = {''};
str2(13) = {''};
str2(14) = {'Jim Price, January, 2001.'};

hf3 = figure(1);
clf
% set(hf3,'Position',[50 200 500 500])
set(gca,'Visible','off')
text(0, 0.50, str2,'FontSize',12,'HorizontalAlignment','left')

pause
close(1)

g = 9.8;
%  set important inital conditions here
H = 100.;                  %  the nominal layer thickness, m
delH = 10.;                %  amplitude of the thickness change across
                           %    an unbalanced front (initially at rest)
delr = 1.e-3;              %  fractional density change          
c = sqrt(delr*g*H);        %  shallow water wave phase speed


% this code started off as geo-adjustment, and so we have
%   a few lines of relic code
lat = 0.;
% lat = input(' Enter the latitude (degrees, 0 is OK)'); 
f = 2*7.29e-5*sin(lat*pi/180);
IP = 2*pi/f;

ndays = 15;                %  number of days to integrate

ndays = 12.

dt = 0.005*8.64e4;         %  time step, seconds
nstep = round(ndays*8.64e4/dt);  %  number of time steps to take

Rd = c/(f + 1.e-10);       %  Rossby radius of deformation

dx = 10.;                %  grid interval, km
L = 2000.;               %  domain half-width, km  (the larger, the better)
xu = -L:dx:L;            %  grid points for u,v,h (input as km)
xu = 1000.*xu;           %  convert to m
dx = 1000*dx;         
 
CFL = c*dt/dx            %  CFL number (needs to be kept quite small, < 0.2)
  
K = 1.e2;                %  horizontal diffusion (keep this small)
Kd = sqrt(K*ndays*8.64e4)/1000.;  
                         
udecay = (1./(5.*8.64e4));    %  Rayleigh drag  (not used at present)
udecay = 0.;

nonlinu = 1.;            %  set these both to 0 to get a linear model
nonlinh = 1.;            %   (for advection of momentum and thickness)

[m nu] = size(xu);
xh = ones(1,nu-1);
xh = cumsum(xh)*dx;
nh = nu - 1;

mov = 1;
%mov = input(' Would you like to make a movie  ? (1 = yes) ')
if mov == 1 
M = moviein(round(nstep/20));
end

%  pre-allocate some arrays 
u = zeros(1,nu); unew = u; uold = u; uold2 = u;  ut = u;
v = u; vnew = u; vold = u; vold2 = u; vt = u;
h = u; ht = u; hnew = u; hold = u; hold2 = u; 
na = 0; ha = u; ua = u; va = u;

%  pre-allocate some arrays for storage of diagnostics
uc = zeros(1,nstep);
hc = uc; vc = uc; he = uc; ue = uc; ve = vc;
ns9 = floor(nstep/20);    % the constant here is for decimation of time
tss = zeros(1, ns9); hs = zeros(nu, ns9); vs = zeros(nu, ns9);
ns = 0; hs = zeros(nu,round(nstep/20)); vs = hs; 
ts = zeros(1,round(nstep/20));

%  Specify the IC on h, the interface displacement.

Hscale = 30.e3;   %  scale over which h changes, 10-30 km is OK
xl = 500e3;
for m = 1:nh
   xp = xh(m) - L*1000;
   if xp <= -10; h(m) = delH*(0.5 + 0.5*tanh((xp + xl)/Hscale)); end;
   if xp > -10; h(m) = delH*(0.5  + 0.5*tanh(-(xp - xl)/Hscale)); end;
   %   h(m) = 3.*delH*exp(-xu(m)^2/Hscale^2);
   h0(m) = h(m);
end
hold = h;  hold2 = h;

icvmode = menu('Want a velocity IC?', 'No thanks', ...
   'Yes, vortical mode', ' Yes, divergent mode')

%  add a vortical or divergent mode to the initial field

if icvmode == 2
   for m=1:nu
      if h(m)> 0.01
         v(m) = (c/3)*(xu(m)/xl);
      end
   end
   vold = v;
   vold2 = v; 
end
if icvmode == 3
   for m=1:nu
      if h(m)> 0.01
         u(m) = (c/3)*(xu(m)/xl);
      end
   end
   uold = u;
   uold2 = u; 
end


%  store the initial h and PV

hi = h;
PVi = f./(H + h);


for m = 1:nstep   

%  time step this baby; use a leap-frog method with an occasional 
%   trapezoidal correction to supress time splitting

if m == 1

%  actually, the first time step is Euler forward
   du = 0.5*(diff([0 u]) + diff([u 0]));
   dv = 0.5*(diff([0 v]) + diff([v 0]));
   dh = 0.5*(diff([0 h]) + diff([h 0]));

   du(1) = 0.; du(nu) = 0.; 
   dv(1) = 0.; dv(nu) = 0.; 
   dh(1) = 0.; dh(nu) = 0.;

   ut = f*v - (delr*g*dh/dx + udecay*u + nonlinu*u.*du/dx);
   vt = -(f*u + udecay*v + nonlinu*u.*dv/dx);
   ht = -(H*du/dx + nonlinh*h.*du/dx + nonlinh*u.*dh/dx);
  
   uold = u;
   vold = v;
   hold = h;
  
   u = uold + dt*ut;
   v = vold + dt*vt;
   h = hold + dt*ht;   

else  
   
   mtype = mod(m,100);  %  decide whether to do a trapezoidal step
   
   if mtype >= 1       %  come here for a simple leap frog 
          
   du = 0.5*(diff([0 u]) + diff([u 0]));
   dv = 0.5*(diff([0 v]) + diff([v 0]));
   dh = 0.5*(diff([0 h]) + diff([h 0]));
   du(1) = 0.; du(nu) = 0.; 
   dv(1) = 0.; dv(nu) = 0.; 
   dh(1) = 0.; dh(nu) = 0.;

%  velocity and derivatives needed for an upwind difference

   ul = 0.5*(u + abs(u));
   ur = 0.5*(u - abs(u));

   dhl = diff([0 h]); dhl(1) = 0.; dhl(nu) = 0.;
   dhr = diff([h 0]); dhr(1) = 0.; dhr(nu) = 0.;
   dul = diff([0 u]); dul(1) = 0.; dul(nu) = 0.;
   dur = diff([u 0]); dur(1) = 0.; dur(nu) = 0.;
   dvl = diff([0 v]); dvl(1) = 0.; dvl(nu) = 0.;
   dvr = diff([v 0]); dvr(1) = 0.; dvr(nu) = 0.;


   du2 = [0 diff(u,2) 0]/dx^2;
   dv2 = [0 diff(v,2) 0]/dx^2;
   dh2 = [0 diff(h,2) 0]/dx^2;
      
   ut = f*v + K*du2 - (delr*g*dh/dx + udecay*u + ...
           nonlinu*(ur.*dur/dx + ul.*dul/dx));

   vt = K*dv2 - (f*u + udecay*v + ...
           nonlinu*(ur.*dvr/dx + ul.*dvl/dx));

   ht = K*dh2 - (H*du/dx + nonlinh*h.*du/dx + ...
           nonlinh*(ur.*dhr/dx + ul.*dhl/dx));
   
   uold2 = uold;
   vold2 = vold;
   hold2 = hold;
   
   uold = u;
   vold = v;
   hold = h;
  
   u = uold2 + 2*dt*ut;
   v = vold2 + 2*dt*vt;
   h = hold2 + 2*dt*ht;   
   
   else    %  do a leap-frog/trapezoidal time step
    
   du = 0.5*(diff([0 u]) + diff([u 0]));
   dv = 0.5*(diff([0 v]) + diff([v 0]));
   dh = 0.5*(diff([0 h]) + diff([h 0])); 
   du(1) = 0.; du(nu) = 0.; 
   dv(1) = 0.; dv(nu) = 0.; 
   dh(1) = 0.; dh(nu) = 0.;

%  velocity and derivatives needed for an upwind difference

   ul = 0.5*(u + abs(u));
   ur = 0.5*(u - abs(u));

   dhl = diff([0 h]); dhl(1) = 0.; dhl(nu) = 0.;
   dhr = diff([h 0]); dhr(1) = 0.; dhr(nu) = 0.;
   dul = diff([0 u]); dul(1) = 0.; dul(nu) = 0.;
   dur = diff([u 0]); dur(1) = 0.; dur(nu) = 0.;
   dvl = diff([0 v]); dvl(1) = 0.; dvl(nu) = 0.;
   dvr = diff([v 0]); dvr(1) = 0.; dvr(nu) = 0.;

   du2 = [0 diff(u,2) 0]/dx^2;
   dv2 = [0 diff(v,2) 0]/dx^2;
   dh2 = [0 diff(h,2) 0]/dx^2;
      
   ut = f*v + K*du2 - (delr*g*dh/dx + udecay*u + ...
           nonlinu*(ur.*dur/dx + ul.*dul/dx));

   vt = K*dv2 - (f*u + udecay*v + ...
           nonlinu*(ur.*dvr/dx + ul.*dvl/dx));

   ht = K*dh2 - (H*du/dx + nonlinh*h.*du/dx + ...
           nonlinh*(ur.*dhr/dx + ul.*dhl/dx));
 
   u1 = u + 0.5*dt*ut;
   v1 = v + 0.5*dt*vt;
   h1 = h + 0.5*dt*ht; 
  
   
   du = 0.5*(diff([0 u1]) + diff([u1 0]));
   dv = 0.5*(diff([0 v1]) + diff([v1 0]));
   dh = 0.5*(diff([0 h1]) + diff([h1 0]));
   du(1) = 0.; du(nu) = 0.; 
   dv(1) = 0.; dv(nu) = 0.; 
   dh(1) = 0.; dh(nu) = 0.;

%  velocity and derivatives needed for an upwind difference

   ul = 0.5*(u1 + abs(u1));
   ur = 0.5*(u1 - abs(u1));

   dhl = diff([0 h1]); dhl(1) = 0.; dhl(nu) = 0.;
   dhr = diff([h1 0]); dhr(1) = 0.; dhr(nu) = 0.;
   dul = diff([0 u1]); dul(1) = 0.; dul(nu) = 0.;
   dur = diff([u1 0]); dur(1) = 0.; dur(nu) = 0.;
   dvl = diff([0 v1]); dvl(1) = 0.; dvl(nu) = 0.;
   dvr = diff([v1 0]); dvr(1) = 0.; dvr(nu) = 0.;
 
   du2 = [0 diff(u1,2) 0]/dx^2;
   dv2 = [0 diff(v1,2) 0]/dx^2;
   dh2 = [0 diff(h1,2) 0]/dx^2;
      
   ut = f*v + K*du2 - (delr*g*dh/dx + udecay*u + ...
           nonlinu*(ur.*dur/dx + ul.*dul/dx));

   vt = K*dv2 - (f*u + udecay*v + ...
           nonlinu*(ur.*dvr/dx + ul.*dvl/dx));

   ht = K*dh2 - (H*du/dx + nonlinh*h.*du/dx + ...
           nonlinh*(ur.*dhr/dx + ul.*dhl/dx));
     
  
   u = uold + dt*ut;
   v = vold + dt*vt;
   h = hold + dt*ht;
    
   uold = u;
   vold = v;
   hold = h;
 
   ut = 0.*ut;
   vt = 0.*vt;
   ht = 0.*ht;
   
end  %  if on mtype


end    %  if on m ==1

 
%  Apply a radiation BC to the grid ends.  This seems to be  
%   acceptable but not great since the phase speed changes in time.
   
   rbc = c*dt/dx;
   
   u(1) = u(1) + rbc*(u(2) - u(1));
   v(1) = v(1) + rbc*(v(2) - v(1));
   h(1) = h(1) + rbc*(h(2) - h(1));
   u(nu) = u(nu) + rbc*(u(nu-1) - u(nu));
   v(nu) = v(nu) + rbc*(v(nu-1) - v(nu));
   h(nu) = h(nu) + rbc*(h(nu-1) - h(nu));
   
% store some data for plotting later
 
   ts(m) = m*dt/8.64e4;       %  the time (days) 
 
   hc(m) = h(round(nu/2));    %  central h,u,v
   uc(m) = u(round(nu/2));
   vc(m) = v(round(nu/2));
   
   he(m) = h(round(2*nu/3));  %  h,u,v at some other point
   ue(m) = u(round(2*nu/3));
   ve(m) = v(round(2*nu/3));

   na = na + 1;               %  sum for averages
   ha = ha + h;
   ua = ua + u;
   va = va + v;

%  evaluate the energy balance

ke(m) = sum(0.5*(H + h).*(u.^2 + v.^2));
pe(m) = sum(0.5*g*delr*(h.^2));
   
if mod(m,20) == 1    %  save the entire grid occasionally
   ns = ns+1;
   hs(:,ns) = h(:);
   vs(:,ns) = sqrt(v(:).^2 + u(:).^2);
   tss(ns) = ts(m);

%  make a movie frame, occasionally

if mov == 1
hf1 = figure(1);
clf reset

subplot(2,1,1)
plot(xu/1000, u/c, xu/1000, v/c)
axis([-2.e6/1000. 2.e6/1000. -0.3 0.3])  %  may need to change these 
legend('across-front', 'along-front',2)
ylabel('speed/C')
title('1-D adjustment under gravity')

subplot(2,1,2)
plot(xu/1000., h/delH); 
fill([xu(1)/1000 xu/1000 max(xu)/1000], [-9 h/delH -9], [0.6 0.6 0.8]);
axis([-2.e6/1000. 2.e6/1000. -0.5 1.5])  %  may need to change these for each case
ylabel('displacement/\deltah')
xlabel('distance, km ')
timed = m*dt/8.64e4;
text(0., 1.2, ['days = ', num2str(timed)])

M(:,ns) = getframe;


end    %  if on movie or not

end    %  if on plot this time step ot not
   
end    %  end of the time stepping loop

%  compute time averages

ha = ha/na;
ua = ua/na;
va = va/na;

%  plot some things

figure(2)
clf reset

subplot(2,1,1)
plot(ts, uc/c, 'b', ts, vc/c, '--r')
ylabel('U(0), m/sec')
legend('U','V')
title('Currents at x = 0 and 2/3')

subplot(2,1,2)
plot(ts, ue/c, 'b', ts, ve/c, '--r')
xlabel('time, days'); ylabel('U(2/3), m/sec')
legend('U','V')

figure(3)
clf reset
mesh(tss, xu/1000., hs/delH)
xlabel('time, days')
ylabel('distance, km')
zlabel('height/\deltah')
axis([0 20 -2000 2000 -0.5 1.5])
title('Adjustment under gravity; height') 


figure(5)
clf reset
Esc = g*delr*delH^2*xl/dx;
plot(ts, ke/Esc, ts, pe/Esc, ts, (ke+pe)/Esc)
legend('KE', 'PE', 'KE+PE')
xlabel('time, days')
ylabel('energy/PE(0)')
title('Energy budget')


% reshuffle the figures
figure(5)
figure(3)
figure(2)

% replay the movie 
if mov == 1
figure(1)
movie(M, 1)
end

