
%  geoadj_1d.m; adjustment under gravity, rotation and friction in 1-D
%
%  A 1-D numerical ocean model of a single
%  shallow layer of rotating, incompressible fluid.
%  Useful for showing the numerical version of geostrophic
%  adjustment. The numerical methods used here are adequate
%  but not great:  an A-grid (velocity and height at the
%  same grid points) and upwind differencing on advection
%  terms help preserve a shock, but at the expense of high 
%  numerical viscosity. High frequency ripples on 
%  ridge height should not be taken literally. Time-stepping is by the  
%  leap-frog/trapezoidal time method.  In all, 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) in
%  the case of low latitude adjustment to a large initial ridge.  
%
%  by Jim Price,  January, 2001.  
%  Revised in June 2004 to refine the plotting.
%  In March, 2010 to add linear drag and gui interface.

clear all
close all
format compact
set(0,'DefaultLineLineWidth',1.5)
set(0,'DefaultTextFontSize',14)
set(0,'DefaultAxesLineWidth',1.4)
set(0,'DefaultAxesFontSize',14)

%  text to display: 
str2(1) =  {'geoadj\_1d.m' };
str2(2) =  {'Geostrophic adjustment in one dimension  '};
str2(3) =  {'solved by a primitive equation numerical '};
str2(4) =  {'model representing one (shallow) layer of'};
str2(5) =  {'fluid having reduced gravity.  Five parameters '}; 
str2(6) =  {'may be modified below:'};
str2(7) =  {'  '};
str2(8) =  {'latitude, 0 -- 90 degrees'}; 
str2(9) =  {'layer thickness, 50 -- 5000 meters'}; 
str2(10)=  {'half-width of the ridge, 20 -- 500 kilometers'};
str2(11)=  {'height of the ridge, eta_0, 1 -- 500 meters'};
str2(12)=  {'e-folding by Rayleigh drag, 1 -- 99999 days'};
str2(13)=  {'ndays to integrate, 3 - 20 days. '};

hf3 = figure(9);
clf
set(hf3,'Position',[50 200 600 700])
set(gca,'Visible','off')
text(-.05, 0.75, str2,'FontSize',14,'HorizontalAlignment','left');

latitude = 20.;
hlat = uicontrol('Style', 'Edit', 'Units', 'Normalized', ...
	'Position',  [0.05 0.37 0.40 0.06], ...
	'FontSize', 14,'fontweight', 'demi', ...
	'Callback', 'eval(get(hlat,''string''))', ...
	'String', 'latitude  = 20 ');

thickness = 500.;
hH = uicontrol('Style', 'Edit', 'Units', 'Normalized', ...
	'Position',  [0.05 0.30 0.40 0.06], ...
	'FontSize', 14, 'fontweight', 'demi',...
	'Callback', 'eval(get(hH,''string''))', ...
	'String', 'thickness  = 500  ');

width = 200.;
hwid = uicontrol('Style', 'Edit', 'Units', 'Normalized', ...
	'Position',  [0.05 0.23 0.40 0.06], ...
	'FontSize', 14, 'fontweight', 'demi',...
	'Callback', 'eval(get(hwid,''string''))', ...
	'String', 'width  = 200  ');

eta0 = 20.;
heta = uicontrol('Style', 'Edit', 'Units', 'Normalized', ...
	'Position',  [0.05 0.16 0.40 0.06], ...
    'FontSize', 14,'fontweight', 'demi', ...
	'Callback', 'eval(get(heta,''string''))', ...
	'String', 'etao = 20 ');

efold = 999.;
hefold = uicontrol('Style', 'Edit', 'Units', 'Normalized', ...
	'Position',  [0.05 0.09 0.40 0.06], ...
	'Callback', 'eval(get(hefold,''string''))', ...
	'FontSize', 14, 'fontweight', 'demi',...
	'String', 'efold = 999  ');

ndays = 6.;
hndays = uicontrol('Style', 'Edit', 'Units', 'Normalized', ...
	'Position',  [0.05 0.02 0.40 0.06], ...
	'Callback', 'eval(get(hndays,''string''))', ...
	'FontSize', 14,'fontweight', 'demi', ...
	'String', 'ndays = 6 ');
hgo = uicontrol('Style', 'Pushbutton', 'Units', 'Normalized', ...
	'Position',  [0.80 0.10 0.1 0.06], ...
	'Callback', 'uiresume(9)', ...
	'FontSize', 14,'fontweight', 'demi', ...
	'String', 'start');
hquit = uicontrol('Style', 'Pushbutton', 'Units', 'Normalized', ...
	'Position',  [0.80 0.02 0.1 0.06], ...
	'Callback', 'lquit = 1; uiresume(9)', ...
	'FontSize', 14,'fontweight', 'demi', ...
	'String', 'end');

lquit = 0;
while (lquit == 0)
disp  ('  ... waiting for control panel action.')
uiwait(9)
if lquit == 1; break; end;


clear h u 

%  define the basic statification
H = thickness;
g = 9.8;
delr = 2/1030;            %  fractional density change          
c = sqrt(delr*g*H);       %  shallow water phase speed
disp(' the non-disersive phase speed, m/s is ')
disp(c)

% ICtype = menu('What sort of IC?', ...
%   'initial interface displacement', 'initial current')
ICtype = 1;   % initial interface displacement 

% lat = input(' Enter the latitude (degrees, 0 < lat < 90) '); 
f = 2*7.29e-5*sin(latitude*pi/180);
IP = 2*pi/(f + 1.e-10);   % an inertial period, seconds

%  set a flag to indicate whether rotation is 'on'
rota = 1;
if abs(latitude) <= 0.1
   rota = 0;
end

%  calculate the Radius of deformation
Rd = 1.e10;
if rota == 1
Rd = c/(f + 1.e-10);        %  Rossby radius of deformation
disp(' the Rossby radius of deformation, km  ')
disp(Rd/1000)
end

%  the velocity scale
Ugeo = g*delr*eta0/(f*Rd);
Usc = (eta0/H)*c;
disp(' the velocity scale, m/s ')
disp(Usc)

disp(' the Mach number, Usc/c ')
disp(Usc/c)

%  define the grid in x
% dx = Rd/10;               %  grid interval, in units of Rd
% L = 20.*Rd;               %  domain half-width, Rd  (the larger, the better)

%  in general, probably better to do this.....
% if rota == 0
    dx = 5.e3;
    L = 1100.e3;
    Ld = 1000.e3;
% end
xu = -L:dx:L;            %  grid points for u,v,h
nx = length(xu);
n1e =  round((L - Ld)/dx);
n2e = nx - round((L - Ld)/dx);

for k=1:nx
    Ldindex(k) = 1;
    if k <= n1e;  Ldindex(k) = 0; end
    if k >= n2e;  Ldindex(k) = 0; end
end

%  define the Rayleigh drag from the given efolding time
udecay = (1./(efold*8.64e4));  %  Rayleigh drag
Ekman = udecay/abs(f);
disp(' The Ekman number    ')
disp(Ekman)

% use CFL to set the time step
CFL = 0.05;                 %  set the CFL number, used to set dt
dt = CFL*dx/c;              %  time step, seconds
nstep = round(ndays*8.64e4/dt);  %  number of time steps to take

%  set a flag to make a movie (1) or not (0)
mov = 1;
if mov == 1 
M = moviein(round(nstep/30));
end

% define horizontal diffusion;  helps reduce small scale numerical noise
K = 1.e2;   %  horizontal diffusion; keep as small as possible 
disp('  the diffusion length scale (ndays), km   ') 
Kd = sqrt(K*ndays*8.64e4)/1000.;
disp(Kd)

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

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

%  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; holda = u; hold2 = u; 
na = 0; ha = u; ua = u; va = u;

%  store and plot the data at this interval 
nplot = 50;

%  pre-allocate some arrays for storage of diagncostics
uc = zeros(1,nstep);
hc = uc; vc = uc; he = uc; ue = uc; ve = vc; ts = uc;
hf = he; uf = ue; vf = ue; ke = ue; pe = ue; fw = pe; rw = pe;
ns9 = floor(nstep/nplot) + 1;    % for data that are decimated in time
ts9 = zeros(1, ns9);
hs = zeros(nu, ns9); 
% vs = hs; 
ns = 0;    %  a counter used to keep track of how many stored  

%  Specify the IC, first on h, the interface displacement above the
%       background.
xl = width*1000.;      %  convert width to m 
for m = 1:nh
   xp = xh(m) - L;
   if abs(xp) <= xl;  h(m) = eta0;   %  this makes a box-shaped ridge
   else
       h(m) = 0;
   end
%  if xp < 0; h(m) = delH*(0.5 + 0.5*tanh((xp + xl)/Hscale)); end;
%  if xp > 0; h(m) = delH*(0.5  + 0.5*tanh(-(xp - xl)/Hscale)); end;
end


%  add an IC other than zero on the current 
ICtype = 0
if ICtype == 2
    for j=1:nh
  v(j) = (eta0*f/H)*xu(j);
  if xu(j) <= -width*1000; v(j) = -eta0*f*width*1000/H; end;
  if xu(j) >= width*1000; v(j) = eta0*f*width*1000/H; end;
  vic(j) = v(j);
    end
end

%  now smooth the ridge a little
a4 = 0*h;
a5 = a4;
 for m=2:nh        
    a4(m) = 0.33*(h(m-1) + h(m+1) + h(m));
    a5(m) = 0.33*(v(m-1) + v(m) + v(m+1));
 end
 a4(1) = h(1);  a4(nh+1) = h(nh+1);
 h = a4;
 v = a5;

%  initialize time levels
holda = h;
hold2 = h;
h0 = h;
vold = v;
vold2 = v;
v0 = v;

%  store the initial h and PV
hi = h;
PVi = f./(H + h);

%
%  now we are ready to start time stepping
%  the first step has to be 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;
   holda = h;
  
   u = uold + dt*ut;
   v = vold + dt*vt;
   h = holda + dt*ht;  
 
hf1 = figure(1);
clf reset  
     
%  now for the real integration.....
for m = 1:nstep   

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

   trap = 20;
   mtype = mod(m,trap);  %  decide whether to do a trapezoidal step (20 - 100)
   
   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 = holda;
   
   uold = u;
   vold = v;
   holda = 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 = holda + dt*ht;
    
   uold = u;
   vold = v;
   holda = h;
 
   ut = 0.*ut;
   vt = 0.*vt;
   ht = 0.*ht;
   
end  %  if on mtype (whether leap or trap)

%  Apply a radiation BC to the grid ends.  This is only moderately 
%   succesful because the phase speed changes in time.
   
   bcfactor = 1.2;     % an ad hoc 'correction' for the radiation BC
   rbc = bcfactor*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)
   np = length(u);
   nf = np - round(xl/dx);
   hf(m) = h(nf);    %  far field h,u,v
   uf(m) = u(nf);
   vf(m) = v(nf);
   
   ne = round(np/2 + round(xl/dx));
   he(m) = h(ne);  %  h,u,v at a point near the edge of the front
   ue(m) = u(ne);
   ve(m) = v(ne);

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

%  evaluate the energy balance at each time step over the
%      interval -Ld to Ld

keall = 0.5*(H + h).*(u.^2 + v.^2);
peall = 0.5*g*delr*(h.^2); 
ke(m) = dx*sum(Ldindex.*keall);
pe(m) = dx*sum(Ldindex.*peall);
fw3 =  udecay*(u.^2 + v.^2).*(H + h);
fs1 = -dt*dx*sum(Ldindex.*fw3);

%  the radiation (pressure work) evaluated at +- 1000 km.
rad = dt*( u(n1e)*g*delr*h(n1e)*(H + h(n1e)) ...
    - u(n2e)*g*delr*h(n2e)*(H + h(n2e)) );

% sum the friction and radiation
if m==1  
    fw(m) = fs1;
    rw(m) = rad;   
else
    fw(m) = fw(m-1) + fs1;
    rw(m) = rw(m-1) + rad;  
end

% other diagnostic data are decimated in time

if mod(m,nplot) == 1    %  save the entire grid occasionally
   ns = ns+1;
   hs(:,ns) = h(:);
%  vs(:,ns) = sqrt(v(:).^2 + u(:).^2);
   ts9(ns) = ts(m);

% plot the solution, occasionally
figure(1)
clf 

set(gcf,'DoubleBuffer','on')
set(gcf, 'position', [250 250 800 400])
axis([-1000 1000 -0.5 3.5])

vv = axis;
xyscale = (vv(2) - vv(1))/(vv(4) - vv(3));
for q=1:4:nu
   hold on
   % xa = xu(q)/1000.; 
   if rota == 1
   xa = xu(q)/1000.;
else 
    xa = xu(q)/1.e3;
   end
   Lplot = 1000.e3;    % this is the scale of the plot
   ya = 2.5;
   xb = xa + (2*Lplot/(1000*4))*uold(q)/Usc;
%   xb = xa + 7.5*uold(q)/Usc;    %  this 5 depends upon the axis set below
   yb = ya + vold(q)/Usc;
if rota == 0            %  switch the components if there is no rotation
   xb = xa + 7.5*vold(q)/Usc;
   yb = ya + uold(q)/Usc;
end
   plot([xa xb], [ya yb], 'LineWidth', 1.0)
end


plot(xu/1000., hi/eta0,  'g', 'linewidth', 1.0)   % the IC


if rota == 0
 %   title(['gravitational adjustment, latitude = ', num2str(latitude,2)])
else
% title(['geostrophic adjustment; latitude = ', num2str(latitude,2), ...
 %     ', Ekman = ', num2str(Ekman,2)]) 
end

%  add the Rd scale
if abs(latitude) >= 0.1
plot([-500 (-500 + Rd/1000)], [1.5 1.5], 'k')
text(-650, 1.5, 'Rd')
end

hold on

if rota == 1
 plot(xu/1000, holda/eta0); 
 fill([xu(1)/1000 xu/1000 max(xu)/1000], [-9 holda/eta0 -9], [0.6 0.6 0.8]);
else
   plot(xu/1.e3, holda/eta0, 'EraseMode', 'none'); 
 fill([xu(1)/1.e3 xu/1.e3 max(xu)/1.e3], [-9 holda/eta0 -9], [0.6 0.6 0.8]);
end
   ylabel('\eta/\eta_0  and  U/(C \eta_0/H)')
   if rota == 1
      ylabel('\eta/\eta_0  and  V/(C \eta_0/H)') 
  end

xlabel('X distance, km', 'EraseMode', 'none')

timed = (floor(100*m*dt/8.64e4))/100.;
if rota == 1
	text(250., 1.5, ['time = ', num2str(timed), ' days'])
   tIP = timed*8.64e4/IP;
   timeIP = (floor(100*tIP))/100.;
   text(250., 1.2, ['        = ', num2str(timeIP), ' IP'])
else 
	text(250., 1.5, ['time = ', num2str(timed), ' days'])
end
  % text(15., 1.3, ['latitude = ', num2str(lat)])
  
drawnow


mov = 0;

if mov == 1
M(:,ns) = getframe(gcf, [0 0 800 400]);
end


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


%  plot some current components
figure(2)
clf reset

tip = ts;
subplot(2,1,1)
plot(ts, uf/Usc, 'b', ts, vf/Usc, '--r')
ylabel('U far/(C \eta_0/H)')
legend('U','V')
title('Currents in the far field (x=900 km) and on the ridge edge')

subplot(2,1,2)
plot(ts, ue/Usc, 'b', ts, ve/Usc, '--r')
xlabel('time, days'); ylabel('Uedge/(C \eta_0/H)')
legend('U','V')


%  3-d plot of eta
figure(3)
clf reset
set(gcf, 'paperposition', [1 1 6 6])
mesh(ts9, xu/1000., hs/eta0);
xlabel('time, days')
ylabel('X distance, km')
zlabel('\eta/\eta0')
axis([0 ndays -Ld/1000 Ld/1000 -0.5 1.5])
% title('geostrophic/gravitational adjustment in 1-d') 
view([-70 60])
% add the gravity wave speed 
tx = 0:0.1:5;
cx = tx*8.64e4*c;
zx = 0.55*ones(1,length(tx));
hold on
plot3(tx, cx/1000., zx, 'r', tx, -cx/1000., zx, 'r')


if rota == 1
%  compute the final PV

dv = 0.5*(diff([0 va]) + diff([va 0]));
dv(1) = 0.; dv(nu) = 0.; 
rvort = dv/dx;
PVf = (f + rvort)./(H + h);
PVflin = (f + 0*rvort)./(H + h);
%
% figure(4)
% clf reset
% subplot(2,1,1)
% plot(xu/1000., hi/eta0, 'g', xu/1000., h/eta0,'r'); 
% axis([-Ld/1000. Ld/1000 -1.2 1.2]);
% ylabel('h/\eta_0')
% legend('initial', 'final')
% legend boxoff
% title('Initial and final thickness and PV')
% 
% subplot(2,1,2)
% plot(xu/1000., PVi,'g', xu/1000., PVf, 'r', xu/1000., PVflin, 'b'); 
% axis([-L/1000 L/1000  0.95e-7 1.02e-7]); 
% xlabel('distance, km')
% ylabel('Potential Vorticity')
% legend('initial', 'final','f/ H + \eta')
% legend boxoff
% end

% average current and rel vorticity

figure(5)
clf reset
subplot(2,1,1)
plot(xu/1000., ua/Usc, xu/1000., va/Usc, 'g')
ylabel('U, V avg speeds/(C \eta_0/H)')
legend('Uavg', 'Vavg')
axis([-L/1000 L/1000 -0.5 0.5]) 
title('Time average currents, etc.')

subplot(2,1,2)
plot(xu/1000., rvort/f)
xlabel('distance, km')
ylabel('rel vorticity/f')
axis([-L/1000 L/1000 -eta0/H eta0/H])

end    % if on rotation or not

%  energy budget
figure(8)
clf reset
set(gcf, 'paperposition', [0 0 5 5])
Esc = pe(1); 
plot(ts, ke/Esc, 'r', ts, pe/Esc, 'b', ts, (ke + pe)/Esc, 'g', ...
    ts, fw/Esc, 'r--', ts, rw/Esc, 'm')
hl = legend('KE', 'PE','KE + PE', 'FW', 'PW', ...
       'Location', 'SouthWest')
legend boxoff
set(hl, 'fontsize', 14)
xlabel('time, days')
ylabel('energy/PE_0')
hold on
plot([0 max(ts)], [0 0], 'k--')
axis([0 max(ts) -1.2 1.2])


L = 1000000;
figure(10)
clf reset
qscale = f/H;
set(gcf,'paperposition', [0 0 5 8.5])

subplot(3,1,2)
%plot(xu/1000., va/Usc) % , xu/1000., ua/Usc, 'r')
ylabel('V/(C \eta_0/H)')
%legend('Vavg', 'Uavg')
legend boxoff
axis([-L/2000 L/2000 -0.5 0.5])
grid on; box on
hold on
for k = 1:3:length(xu)
    plot([xu(k)/1000 xu(k)/1000], [0 va(k)/Usc], 'b')
end


subplot(3,1,3)
plot(xu/1000., PVi/qscale,'g', xu/1000., PVf/qscale, 'b', ...
         xu/1000., PVflin/qscale, 'r'); 
axis([-L/2000 L/2000  0.95 1.02]); 
ylabel('q/(f/H)')
hl = legend('q_0', 'q_{final}','q_{Sverdrup}')
set(hl,'fontsize', 14)
legend boxoff
xlabel('distance, km')
grid on

%  radius of deformation
g = 9.8;
dr = 2.;
r0 = 1000.;
H = 500.;
lat = 20.;
f = 2.*7.292e-5*sin(lat*pi/180);

C = sqrt(g*dr*H/r0);
Rd = C/f;

Ly = 200.;
y = -2.5*Ly:2:2.5*Ly;
y = 1000*y;
Ly = Ly*1000.;

eta0 = 20.
for j=1:length(y)

if y(j) >= 0    
eta(j) = (eta0/2)*(1 - exp((y(j) - Ly)/Rd));
if (y(j) - Ly) >= 0
eta(j) = eta0/2*(exp(-(y(j) - Ly)/Rd) - 1);
end
end

if y(j) <= 0
eta(j) = (eta0/2)*(exp((y(j) + Ly)/Rd) - 1);
if (y(j) + Ly) >= 0
eta(j) = eta0/2*(1 - exp(-(y(j) + Ly)/Rd));
end
end

end
eta = eta + eta0/2

subplot(3,1,1) 
plot(xu/1000, hi/eta0, 'g',  xu/1000, h/eta0, 'b', y/1000, eta/eta0, 'r')
hl = legend('\eta_0', '\eta_{final}', '\eta_{q conserv}'); 
set(hl,'fontsize', 14)
legend boxoff
axis([-L/2000 L/2000 -0.2 1.2])
grid on
ylabel('\eta/\eta_0')


figure(9)

% % if you want to see the movie again
% replayn = 0;
% replayn = input('  To replay the movie, enter 1, 2 or 3 (0 = no)  ');
% if replayn >= 1
%     figure(1); 
%     set(gcf, 'paperposition', [0 0 850 400])
%     set(gca, 'position', [0 0 1 1])
%     movie(M,replayn,6)
% end

end    %  end of the function


% 
% 
% imov = input('  Shall we save an mpeg movie file?  (1 = yes)');
% if imov == 1
% if  exist('geoadjmovie.m') == 1; 
%     delete 'geoadjmovie.mpg';  
% end
%     mpgwrite(M, colormap, 'geoadjmovie.mpg')
% end

