% 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.
%
% Revision history:
% June 2004 to refine the plotting.
% March 2010 to add linear drag and gui interface.
% September 2010 to fix an error in input of eta0.
% February 2011 to allow more complete specification of the
% initial conditions, i.e., negative eta0, etc.
% October 2011 to add floats.
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. Six parameters '};
str2(6) = {'may be modified or defaulted:'};
str2(7) = {' '};
str2(8) = {'latitude, -90 -- 90 degrees'};
str2(9) = {'layer thickness, 50 -- 5000 meters'};
str2(10)= {'density difference, 0.1 - 10 kg/m^3'};
str2(11)= {'half-width of the ridge, 20 -- 500 kilometers'};
str2(12)= {'height of the ridge, eta_0, -200 -- 200 meters'};
str2(13)= {'e-folding by Rayleigh drag, 1 -- 99999 days'};
str2(14)= {'ndays to integrate, 3 - 20 days. '};
hf3 = figure(9);
clf
set(hf3,'Position',[50 200 600 800])
set(gca,'Visible','off')
text(-.05, 0.80, str2,'FontSize',14,'HorizontalAlignment','left');
latitude = 30.;
hlat = uicontrol('Style', 'Edit', 'Units', 'Normalized', ...
'Position', [0.05 0.44 0.40 0.06], ...
'FontSize', 14,'fontweight', 'demi', ...
'Callback', 'eval(get(hlat,''string''))', ...
'String', 'latitude = 30 ');
thickness = 500.;
hH = uicontrol('Style', 'Edit', 'Units', 'Normalized', ...
'Position', [0.05 0.37 0.40 0.06], ...
'FontSize', 14, 'fontweight', 'demi',...
'Callback', 'eval(get(hH,''string''))', ...
'String', 'thickness = 500 ');
delrho = 2.0;
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', 'delrho = 2.0 ');
width = 100.;
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 = 100 ');
eta0 = 50.;
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', 'eta0 = 50 ');
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 = 4;
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 = 4 ');
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 = delrho/1000.; % normalized density difference
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 + eps); % an inertial period, seconds
% set a flag to indicate whether rotation is 'on'
rota = 1;
if abs(latitude) <= 0.00001
rota = 0;
end
% calculate the Radius of deformation
Rd = Inf;
if rota == 1
Rd = c/abs((f + eps)); % Rossby radius of deformation
disp(' the Rossby radius of deformation, km ')
disp(Rd/1000)
end
W = width/Rd
disp(' the non-d with, L/Rd')
disp(W)
abseta0 = abs(eta0)
% the velocity scale
Ugeo = g*delr*abseta0/(f*Rd);
Usc = (abseta0/H)*c;
disp(' the velocity scale, m/s ')
disp(Usc)
Mach = Usc/c
disp(' the Mach number, Usc/c ')
disp(Mach)
% 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 since we will allow f = 0
% if rota == 0
dx = 3.e3; % grid interval
L = 1200.e3; % calculate on this model domain (bigger is better)
Ld = 500.e3; % plot and diagnose on this smaller domain
% 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
% some initial float postions
dxf = 100.e3;
nf = round(2*Ld/dxf)+1;
for jf=1:nf
xfo(jf) = (jf-1)*dxf - Ld;
zf(jf) = -1.1; % float 'depth' for plotting purposes
end
xf = xfo; % time-dependent float position
% define the Rayleigh drag from the given e-folding time
udecay = (1./(efold*8.64e4)); % Rayleigh drag
E = udecay/abs(f);
disp(' The Ekman number ')
disp(E)
% use CFL to set the time step
CFL = 0.05; % specify 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
% store and plot the data at this interval
nplot = 75;
% 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 = 30; % horizontal diffusion, mks, 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;
% 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
nsp = 0; % another counter for the movie storage
% 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 and any initial v a little
for jj=1:6
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;
end % loop on smoothing the IC
% 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.1; % 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(400*1000/dx);
hff(m) = h(nf); % far field, x = 400 km, h,u,v
uff(m) = u(nf);
vff(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 ridge
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.
% the h here is eta of the paper........
rad = dt*g*delr*( u(n1e)*h(n1e)*(H + h(n1e)) ...
- u(n2e)*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
% integrate the float postions
uf = interp1(xu, u, xf, 'linear'); % find x-velocity at the float postions
xf = xf + uf*dt; % step ahaead (this is awfully simple.....)
% 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, every 2*nplot steps
figure(1)
clf
set(gcf,'DoubleBuffer','on')
set(gcf, 'position', [250 250 800 400])
axis([-Ld/1000. Ld/1000. -1.25 3.25])
vv = axis;
xyscale = (vv(2) - vv(1))/(vv(4) - vv(3));
for q=1:2:nu
hold on
% xa = xu(q)/1000.;
if rota == 1
xa = xu(q)/1000.;
else
xa = xu(q)/1.e3;
end
Lplot = Ld; % 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/abseta0, 'g', 'linewidth', 1.0) % the IC
if rota == 0
text(-425., 3.0, 'these currents are rotated 90^o ccw')
end
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([-300 (-300 + Rd/1000)], [1.5 1.5], 'k')
text(-370., 1.5, 'Rd')
end
hold on
if rota == 0
E = NaN;
end
if rota == 1
plot(xu/1000, holda/abseta0);
fill([xu(1)/1000 xu/1000 max(xu)/1000], [-9 holda/abseta0 -9], [0.6 0.6 0.8]);
else
plot(xu/1.e3, holda/abseta0, 'EraseMode', 'none');
fill([xu(1)/1.e3 xu/1.e3 max(xu)/1.e3], [-9 holda/abseta0 -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(200., 1.5, ['time, days = ', num2str(timed,2)])
tIP = timed*8.64e4/IP;
timeIP = (floor(100*tIP))/100.;
text(200., 1.1, ['time, IP = ', num2str(timeIP,2)])
else
text(200., 1.5, ['time, days = ', num2str(timed,2)])
end
plot(xu/1000., hi/abseta0, 'g', 'linewidth', 1.0) % the IC
% plot the float present and initial positions
plot(xf/1000, zf, 'b*', 'markersize', 5)
plot(xfo/1000, zf, 'g+', 'markersize', 5)
% title1 = ['geoadj\_1d: lat = ', num2str(latitude,2), '^o, C = ', ...
% num2str(c,2), ...
% ' m s^{-1}, R_d = ', num2str(Rd/1000, 2), ...
% ' km, H = ',num2str(thickness,3), ' m, \eta_0 = ', ...
% num2str(eta0, 2), ' m'];
%
W = width*1000/Rd;
title2 = ['C \eta_0/H = ', num2str(Usc, 2), ' m s^{-1},' ...
' L/R_d = ', num2str(W,2), ', Mach no. = ', ...
num2str(Mach,2), ', Ekman no. = ', num2str(E,2)];
% title([{title1, title2}])
drawnow
mov = 1;
if mov == 1 & mod(ns,2) == 1
nsp = nsp + 1;
M(:,nsp) = 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, uff/Usc, 'b', ts, vff/Usc, '--r')
ylabel('U far/(C \eta_0/H)')
legend('U','V')
title('Currents in the far field (x=400 km) and at 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])
% interpolate eta onto a reasonably-sized grid before plotting
nt9 = length(ts9);
xd9 = -Ld:2000.:Ld;
nxd9 = length(xd9);
clear etas9 eta9
etas9 = zeros(nxd9, nt9);
for j=1:nt9
bc = hs(:,j);
eta9(:,j) = interp1(xu, bc', xd9);
end
mesh(ts9, xd9/1000., eta9/abseta0);
xlabel('time, days')
ylabel('X distance, km')
zlabel('\eta/\eta_0')
axis([0 ndays -Ld/1000 Ld/1000 -1.1 1.1])
% title('geostrophic/gravitational adjustment in 1-d')
view([-70 60])
% add the gravity wave speed
t4 = (Ld/c)/8.64e4;
tx = 0:0.1:t4;
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')
% 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])
if rota == 1; % if latitude approx 0, skip the things below
% start and end PV....
% 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
% plot some average currents
% 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 -abseta0/H abseta0/H])
% plot the PV budget here
L = 500000;
figure(10)
clf reset
qscale = f/H;
set(gcf,'papersize', [7 14])
subplot(3,1,3)
%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
xlabel('distance, km')
for k = 1:2:length(xu)
plot([xu(k)/1000 xu(k)/1000], [0 va(k)/Usc], 'b')
end
% compute and plot the geostrophic velocity
nxu = length(xu);
ugeo = zeros(1,nxu);
for k=2:nxu-1
ugeo(k) = g*delr*(h(k+1) - h(k-1))/(2*dx*f);
end
ugeo(1) = 0.;
ugeo(nxu) = 0.;
plot(xu/1000., ugeo/Usc, 'r--', 'linewidth', 2.0)
subplot(3,1,1)
plot(xu/1000., PVi/qscale,'g', xu/1000., PVf/qscale, 'b') %, ...
% xu/1000., PVflin/qscale, 'r');
axis([-L/2000 L/2000 0.88 1.12]);
ylabel('q/(f/H)')
hl = legend('q_0', 'q_{final}') %,'q_{Sverdrup}');
set(hl,'fontsize', 14)
legend boxoff
grid on
% radius of deformation
g = 9.8;
dr = 2.;
r0 = 1000.;
H = thickness;
lat = latitude;
f = 2.*7.292e-5*sin(lat*pi/180);
C = sqrt(g*dr*H/r0);
Rd = C/abs(f);
Ly = width;
y = -2.5*Ly:2:2.5*Ly;
y = 1000*y;
Ly = Ly*1000.;
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,2)
plot(xu/1000, hi/abseta0, 'g', xu/1000, h/abseta0, 'b', y/1000, eta/abseta0, 'r')
hl = legend('\eta_0', '\eta_{final}', '\eta_{geo}');
set(hl,'fontsize', 14)
legend boxoff
axis([-L/2000 L/2000 -1.2 1.2])
grid on
ylabel('\eta/\eta_0')
end % if on rotation or not
% 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
%
% % use this to save an avi movie
% aviobj = avifile('mymovie.avi')
% aviobj = addframe(aviobj,M);
% aviobj = close(aviobj)
%
