Code covered by the BSD License  

Highlights from
a Coriolis tutorial

image thumbnail

a Coriolis tutorial

by

 

16 Feb 2003 (Updated )

Scripts that demonstrate aspects of rotating reference frames, Coriolis force and geostrophy.

gyre_plot.m
% gyre_plot.m
%
% A matlab plotting program to read and plot disk files
% generated by gyre.f 
%
% The first time an experiment is to be plotted, the data that are
% read are ascii files written by the Fortran.  This is really, really
% slow. After this has been done once, you can save the workspace into 
% a mat file, and thereafter you can read the same data in very quickly. 
%
% by Jim Price, 2003, 2010, 2011.
%

clear all
clear Um
close all
format compact

%  Set some default graphics things.
set(0,'DefaultTextFontSize',13)
set(0,'DefaultAxesFontSize',13)
set(0,'DefaultAxesLineWidth',1.5)
set(0,'DefaultLineLineWidth',1.2)

format compact


xscale = 1.0;

readnew = 0;
disp('  Shall we (1) read in the asacii files written by gyre.f, or,')
disp('           (2) read in a previously-saved workspace/mat file? ')
readnew = input('  Enter 1 or 2 here:  ')

if readnew == 1
%  if yes, then come here and start plodding thjough the ascii files

load ga_what.dat    %  some scalar data that define the experiment
what = ga_what;
nx = what(1);
ny = what(2);
nz = what(3);
dx = what(5);
dy = what(6);
x1 = what(7);
L = -x1/1000.;
y1 = what(9); 
f0 = what(11);
beta = what(12);
tau0 = what(13);
H = what(15);
delh = what(16);
efold = what(17);


delh = 100.    %   this is purely a guess

g = 9.8
delr = 2.0
C = sqrt(g*H*delr/1030.)
Rd = C/f0
Clong = -beta*Rd^2
xyear = 365*8.64e4*Clong/1000.



kc = 0;

x = zeros(nx,1); y = zeros(ny,1);

for i=1:nx;
x(i) = x1 + (i-1)*dx;
x(i) = x(i)/1000.;
end;
for i=1:ny;
y(i) = y1 + (i-1)*dy;
y(i) = y(i)/1000.;
f(i) = f0 + y(i)*1000.*beta;
end;

for i=1:nx
    for j=1:ny
        fm(j,i) = f(j);
    end
end


% x = x - mean(x);
% y = y - mean(y);
       
load ga_time.dat    % the times of the matrix data read below
[ntim, kkk] = size(ga_time);

etax = zeros(ntim, nx);
vorx = etax; 
betav = etax;

h1 = zeros(ny, nx); h2 = h1; h3 = h1;  
us = h1, vs = h1;


disp('  reading the eta file...')
load ga_eta.dat;              %  this is the interface displacement
disp('  reading the u file... this is slow')
load ga_u.dat;                %  u component of velocity
disp('  reading the v file...')
load ga_v.dat;                %  v component of velocity
% disp('  reading the tracer file... oh my, this is very slow ...')
% load ga_t.dat;                %  a passive tracer
% disp('  reading the q file...')
% load ga_q.dat;                %  potential vorticity
% disp('  reading the float track file...')
% load ga_float.dat;            %  float tracks
disp('reading the y, f and tau file ...')
load ga_yftau.dat;            %  f and wind stress

%  now ask whether this workspace should be saved for future analysis 
savethis = input('Shall we save this workspace?   (1 = yes)')
if savethis == 1
    uisave 
end

end    %  if on whether to read in a fresh data set

%  come directly here if you want to read in a saved workspace
if readnew == 2
 disp(' Select a previously saved mat file ')
     thisfile = uigetfile 
     load(thisfile)
end

nnh = 0;
nm = 0;

ii = 0;

imov = 0;
% imov = input('  Shall we save avi movie files?  (1 = yes)')

ntime = length(ga_time)


for jj=1:ny
    y6(jj) = ga_yftau(jj,1)/1000.;
    f6(jj) = ga_yftau(jj,2);
    taux(jj) = ga_yftau(jj,3)*1030.;
end
% 
% figure(88)
% plot(f6,y6, 'linewidth', 1.8)

figure(89)
clf reset
hf1 = figure(89);
set(hf1,'paperposition', [0 0 680 500])
hold on
plot([-0.14 0.14], [L/2 L/2], 'r')
plot([-0.14 0.14], -[L/2 L/2], 'r')
plot(taux,y6, 'linewidth', 1.8)
for s=1:10:ny
    plot([0 taux(s)], [y(s) y(s)], 'b')
end
hold on
plot([0 0], [-4000 4000], 'k', 'linewidth', 1.8)
grid on
xlabel('zonal wind stress, Pa')
ylabel('y, north, km')
axis([-0.15 0.15, -3600 3600])
text(0.02, 1800., '{\it westerlies}', 'background', 'w')
text(-0.08, -1800., '{\sl easterlies}', 'background', 'w')
title('zonal wind stress')
ht = text(-0.1, 2600., {'subpolar gyre'; '\nabla \times \tau > 0'})
set(ht, 'color', 'r')
ht = text(-0.1, 0., {'subtropical gyre'; '\nabla \times \tau < 0'})
set(ht, 'color', 'r')
ht = text(0.031, -2600., {'tropical gyre'; '\nabla \times \tau > 0'})
set(ht, 'color', 'r')


   nf1 = 680;
   nf2 = 500;
  
nm = 0;
     for mm=1:10
     nm = nm + 1;
  Um(nm) = getframe(hf1, [0 0 nf1 nf2]); 
     end


%%%%%%%%%%%%% begin stepping through the times here  %%%%%%%%%%%%%%%%%%%%
%   there are so many possibilities here that it is impossible to explain
%   them all........

nt9 = ntime;    %  the total number of saved times
dt = 10         % the interval between plotting 
                  
iplot = [2 6 101]  %   use this to plot specific times 

%  loop on j, the time
for j = nt9              %  j = nt9 does the last time only
    
    %  is is used as the actual time index
 i = nt9       %  this just does the last time 
    
 %  for i = 1:10:nt9     %  for the seasonal case
 % i = iplot(j);
 
 i = nt9       %  this just does the last time 
                      
 n1 = (i-1)*ny + 1;
 n2 = n1 + ny - 1;

 h1(:,:) = ga_eta(n1:n2,:);      %  eta
 h2(:,:) = ga_u(n1:n2,:);        %  east-west velocity component
 h3(:,:) = ga_v(n1:n2,:);        %  the north-south velocity component
% h5(:,:) = ga_t(n1:n2,:);        %  the passive tracer
%  h6(:,:) = ga_q(n1:n2,:);        %  potential vorticity
% vor = h6.*(H + h1) - fm;        %  relative vorticity
 
 us = h2;
 vs = h3;
 
 h19 = h1;
 



% check that mass is conserved

time9(i) = ga_time(i);
mass(i) = mean(mean(h1));



figure(20)
subplot(2,1,1)
hold on
eta3 = h1(100,:) - mean(mean(h1));
plot(x, eta3/(1030/delr))

subplot(2,1,2)
hold on
ht = h1(100,:) + H;
plot(x, ht)

imod = mod(i,5)
if imod == 1 
    if i == 1
        time9(1) = 0.;
    end   
nyc = ceil(ny/2);


curltau = (tau0/1030)*pi/(L*1000);

Hactual = 580.;    %  this the actual H for the three gyre case.......
Svslope = f0*curltau/(Hactual*(g*beta))

 Svslope1 = f0*curltau/((Hactual+80)*(g*beta));
 Svslope2 = f0*curltau/((Hactual-80)*(g*beta));
% 

if i <= 50
figure(27)
hold on
eta3 = h1(nyc,:) - mean(mean(h1));
plot(x, eta3/(1030/delr), 'linewidth', 1.4)
xlabel('x distance, km')
ylabel('SSHA, m')
if time9(i) <= 2000
text(-1900., eta3(1,25)/(1030/delr), num2str(time9(i),4))
end

xs1 = -500.
xs2 = 1.5e3
ys1 = .27
ys2 = ys1 + (Svslope*(xs1 - xs2)*1000.)
plot([xs1 xs2], [ys1 ys2], 'r', 'linewidth', 1.8)
 ys2 = ys1 + (Svslope1*(xs1 - xs2)*1000.)
 plot([xs1 xs2], [ys1 ys2], 'g', 'linewidth', 1.5)
 ys2 = ys1 + (Svslope2*(xs1 - xs2)*1000.)
 plot([xs1 xs2], [ys1 ys2], 'g', 'linewidth', 1.5)
 
 
% 
% ys2 = ys1 + (-Svslope2*xs2*1000.)
% plot([xs1 xs2], [ys1 ys2], 'g', 'linewidth', 1.5)

text(500., 0.24, 'Sverdrup geo slope')

end

end

iwb = 1;   %  set = 1 to plot just the western bc region
   nf1 = 680;
   nf2 = 500;
  
%    Figure 1 starts here 
%  now plot velocity and eta contours
%  now subsample velocity to put the data on a domain better for plotting
xp1 = -L;
dx1 = 360.;
xp2 = L;
yp1 = -L;
dy1 = 360.;
if iwb == 1; dx1 = 20; dy1 = 20. ; end
yp2 = L;
xp = xp1+dx1/2:dx1:xp2;
yp = yp1:dy1:yp2;
xp = xp*xscale;
yp = yp*xscale;
[xp3, yp3] = meshgrid(xp, yp);

hp1 = interp2(x,y,h1,xp3,yp3);
hp2 = interp2(x,y,h2,xp3,yp3);
hp3 = interp2(x,y,h3,xp3,yp3);

speed3 = (hp2.^2 + hp3.^2).^0.5;
vmax = max(max(speed3))

%  now for eta (to add contours to velocity)
dx1 = 60.;
dy1 = 60.;
if iwb == 1; dx1 = 20; dy1 = 20. ; end
xph = xp1+dx1/2:dx1:xp2;
yph = yp1:dy1:yp2;
xph = xph*xscale;
yph = yph*xscale;
[xp3, yp3] = meshgrid(xph, yph);

hp9 = interp2(x,y,h19,xp3,yp3);
 
figure(1)
clf reset
hf1 = figure(1);
set(hf1,'paperposition', [0 0 680 500])

delh = 50.;
etamax = delh;

%  contours of eta
hold on
minmax9 = 500.;
hconts = [-500 -400 -300 -200:20:200 300 400 500];
if iwb == 1; hconts = -160:20:160; end
caxis = ([-minmax9 minmax9]);
% caxis manual

% colm = hp9/delh;
% colm(colm<=-1) = -1;
% colm(colm>=1.) = 1.;

if imov == 1; hp9(1,1) = 2.0; end

hold on
plot([-3400 3400], [L/2 L/2], 'r')
plot([-3400 3400], [-L/2 -L/2], 'r')
[c5, hc5] = contour(xp3, yp3, hp9, hconts); % [-minmax9:20:minmax9]);
set(hc5, 'linewidth', 1.6);
if imov == 0 && i >= 2; clabel(c5, hc5, 'manual'); end
if imov == 1; 
    clear c5 hc5
[c5, hc5] = contour(xp3, yp3, hp9, [-minmax9:100:minmax9]);
clabel(c5, hc5, 'LabelSpacing', 50000);
end
box on
caxis = ([-minmax9 minmax9]);
%colorbar

C = 3.1;
% try scaling velocity to give the little ones a chance
% v = (h2.^2 + h3.^2).^0.5;
% vmax = max(max(v)) + 0.00000001;
% if vmax <= 0.0001; vmax = C*etamax/H; end
vmax = 0.05; 
% % hp2(5,3) = vmax;   % 0.5*C*abs(etamax)/H;
% % hp3(5,3) = 0.;

vv = (h2.^2 + h3.^2).^0.5;
vmax = max(max(vv)); 
disp('  the maximum speed is '); vmax

usc = 10000.;  
if iwb == 1; usc = 15.; end
hq = quiver(xp, yp, usc*hp2, usc*hp3, 0.);
set(hq, 'color', 'b','linewidth', 1.0);
speed9 = vmax;
% % htt = text(xp(3), yp(3)+225*xscale, [num2str(speed9,2), ' m sec^{-1}']);
% % set(htt, 'color', 'w', 'fontsize', 12)

%  the long wave speed at y = 0
 xrbeta = ga_time(i)*8.64e4*Clong/1000. + L;
%% plot([xrbeta; xrbeta], [-L L], 'r', 'linewidth', 2.4)
 
Clongy = Clong*(f/f0).^-2;
Clongy(Clongy<-C) = -C;
 xrbetay = ga_time(i)*8.64e4*Clongy/1000. + L;
 plot(xrbetay', y, 'linewidth', 2.8, 'color', 0.6*[1 1 1])
 
 xx = [xrbetay L L];
 yy = [y' L -L];
 
 fillon = 1
 if imov == 1
 fillon = 0     % to make a movie, you will have to turn fill off
 end
 
 if fillon == 1
hf = fill(xx, yy, 0.9*[1 1 1])
set(hf, 'facealpha', 0.8)
 end

 


%  find the max and min of eta to track the peak
clear a1 a2;
if etamax >= 0.01
maxh = max(max(h1));
end
if etamax <= 0.01
 maxh = min(min(h1));
end
   
[a1] = find(h1 == maxh);
[xm ym] = meshgrid(x, y);
xx = xm(a1(1)); yy = ym(a1(1));

if etamax >= 0.5; 
ht =  text(xx, yy,'H');
end
if etamax <= 0.5; 
ht =  text(xx, yy,'L');
end
set(ht, 'fontsize', 14, 'fontweight', 'bold', 'color', 'k'); 

xpeak(i) = xx;
ypeak(i) = yy;


xlabel('x, east, km')
ylabel('y, north, km') 

  axis([-L L -L L])
  
if iwb == 1;  axis([-3600 -3400 -100 100.]);
    clabel(c5, hc5, 'manual'); end

  
 tim = ga_time(i);
% %  tim = round(tim/100)*100;
% %  tim = (i-1)*100;                 %  this only works for deltat = 100 days!!!!!
%  if tim >= 98 & tim <= 105; tim = 100; end
%   if tim >= 480 & tim <= 520; tim = 500; end
%   if tim >= 99401  & tim <= 10050; tim = 10000; end 
%  if tim >= 4.55; tim = floor(tim) + 1; end
 title(['shallow water beta-plane, wind-driven gyre, time = ', ...
     num2str(tim,5), ' days ']);
  
  
   nf1 = 680;
   nf2 = 500;
   
% % if i == 1
% %   hfig = openfig('tauofy.fig');  
% % set(hfig,'paperposition', [0 0 680 500])
% %     for mm=1:10
% %     nm = nm + 1;
% %  Um(nm) = getframe(hfig, [0 0 nf1 nf2]); 
% %     end
% % end

if imov == 1
nm = nm+1; 
Um(nm) = getframe(hf1, [0 0 nf1 nf2]); 
    pause(1.0)
else
    pause
end

%% compute and plot the streamfunction

if i == 101
psi = zeros(ny, nx);
for jj = nx-1:-1:1
    for kk = 1:ny
        psi(jj,kk) = psi(jj+1,kk) + dx*(h3(jj,kk) + H)*h2(jj,kk);
    end
end

mx = (hp1 + H).*hp2;
my = (hp1 + H).*hp3;
figure(90)
clf reset
[ch, hh] = contour(x, y, -psi/1.e6, [-50:5:50])
set(hh, 'color', [0 1 0])
set(hh, 'linewidth', 1.4)
clabel(ch, hh, 'manual')
xlabel('east distance, km')
ylabel('north distance, km')
hold on
plot([-4000 4000], [-1800 -1800], 'r', 'linewidth', 1.4)
plot([-4000 4000], [1800 1800], 'r', 'linewidth', 1.4)
hq = quiver(xp, yp, 0.7*usc*mx/H, 0.7*usc*my/H, 0.);
set(hq, 'color', 'b','linewidth', 1.4);
title('streamfunction, Sv,  time = 10000 days')
end

nyc = ceil(ny/2);
ysample = 0.
ny8 = nyc + round((ysample*1000.)/dy);
figure(82); 
subplot(2,1,1)
hold on
plot(x, h3(ny8,:))
axis([-3600 3200 -0.05 0.05])
subplot(2,1,2)
hold on
plot(x, h1(ny8,:))
axis([-3600 3200 -5 5])



end    %  time step loop ends here


if imov == 1 
    movie2avi(Um, 'moviefile.avi', 'FPS', 3)
end

%  the next threee lines allow to selecet a single y
%  jylast = 10    %  doe this to select specific latitudes
%  for jy =1:jylast
%  ysample = input(' enter y to sample for PV balance  ')
%  
%  use these next two lines to sweep through all y for the dot plot
 for jy = 1:71     % to 71 does all
 ysample = -3600. + jy*100;   %  to sweep through all the y  
% 

ny8 = nyc + round((ysample*1000.)/dy);

r = 1/(8.64e4*efold);
Sv = pi*tau0/(beta*L*1030*1000);   %  scale for transport per unit width
Cd = 0.01;
H6 = 500.;

trans = 0*x;
curltau = trans; gradh = trans; drag = trans; curlv = trans; 
v6 = trans; u6 = trans; h6 = trans; vgeo = trans; both = trans;
gradh = trans; all = trans; speed6 = trans;
for k = 1: nx 
    v6(k) = h3(ny8,k); u6(k) = h2(ny8,k); h6(k) = h1(ny8,k);
    speed6(k) = sqrt(u6(k)^2 + v6(k)^2);
   %  r6(k) = Cd*speed6(k)/H6;
    r6(k) = r;
    trans(k) = h3(ny8,k)*(h1(ny8,k)+H)/Sv;
    hfac = H/(h1(ny8,k) + H);
    curltau(k) = -(taux(ny8+1) - taux(ny8-1))/(beta*2*dy*1030*Sv);
end
for k=2:nx-1
    gradh(k) = -( h3(ny8,k)*(h1(ny8+1,k) - h1(ny8-1,k))/(2*dy) + ...
         h2(ny8,k)*(h1(ny8,k+1) - h1(ny8,k-1))/(2*dx) ); 
     gradh(k) = gradh(k)*(f(ny8)/(h1(ny8,k) + H))*(H/(beta*Sv));   %%%
     curlv(k) = (h3(ny8,k+1) - h3(ny8,k-1))/(2.*dx) - ...
              (h2(ny8+1,k) - h2(ny8-1,k))/(2.*dy);
%  the next is, I think, a better C grid version.......
          curlv(k) = (h3(ny8,k+1) - h3(ny8,k))/(dx) - ...
              (h2(ny8+1,k) - h2(ny8-1,k))/(2.*dy);
 hmid = h1(ny8,k)+ H;    
     drag(k) = -r6(k)*curlv(k)*H/(beta*Sv);   
    all(k) = curltau(k) + drag(k) + gradh(k);
    vgeo(k) = g*delr*((h1(ny8,k+1) - h1(ny8,k-1))/...
           (f(ny8)*2*dy*1030.));
end
trans(1) = trans(2); trans(nx) = trans(nx-1);
curltau(1) = curltau(2); drag(1) = drag(2);
curltau(nx) = curltau(nx-1); drag(nx) = drag(nx-1);
all(1) = all(2); all(nx) = all(nx-1);
epsy = trans - (curltau + drag);
figure(80), clf reset, hold on
hp = plot(x, trans, 'g.', 'MarkerSize', 6)

hp4 = plot(x, curltau, 'm', x, drag, 'b', x, epsy, 'r');
% plot(x, gradh, 'm')
set(hp, 'linewidth', 1.7)
set(hp4, 'linewidth', 1.7)
xlabel('x distance, km'), ylabel('vorticity tendency/\nabla \times \tau_o')
legend('\beta V h', '\nabla  \times \tau', 'drag', 'epsy')
%   '\nabla  \times  \tau + drag', 'grad h')
legend boxoff
axis([-4000 4000 -1.5 2.5 ])
plot([-4000 4000], [0 0], 'k', 'linewidth', 1.4)
% plot(x, gradh, 'm--')
title(['vorticity balance along y = ',num2str(ysample,4), ' km'])
grid on

figure(60)
clf reset

subplot(2,1,1)
plot(x(x<=-3300), v6(x<=-3300), x(x<=-3300), u6(x<=-3300))
ylabel('speed, m/sec')

subplot(2,1,2)
hold on
hp = plot(x(x<=-3300), trans(x<=-3300), 'g.', 'MarkerSize', 6)
hp4 = plot(x(x<=-3300), curltau(x<=-3300), 'm', ...
    x(x<=-3300), drag(x<=-3300), 'b', x(x<=-3300), epsy(x<=-3300), 'r');
% plot(x, gradh, 'm')
set(hp4, 'linewidth', 1.7)
xlabel('x distance, km'), ylabel('vorticity tendency/\nabla \times \tau_o')
legend('\beta V h', '\nabla  \times \tau', 'drag', 'epsy')
%   '\nabla  \times  \tau + drag', 'grad h')
legend boxoff
% axis([-4000 4000 -1.5 2.5 ])
plot([-3600 -3400], [0 0], 'k', 'linewidth', 1.4)


for m=1:nx;
    isv(m) = 0.;
end

varq = trans.^2 + curltau.^2 + drag.^2 + epsy.^2;

% nonlin balance?
for m=1:nx-1
    epsyn(m) = 0.;
    if abs(epsy(m)) >= abs(curltau(m)); epsyn(m) = epsyn(m) + 1; end
    if abs(epsy(m)) >= abs(drag(m)); epsyn(m) = epsyn(m) + 1; end
    if abs(epsy(m)) >= abs(trans(m)); epsyn(m) = epsyn(m) + 1; end
    if epsyn(m) >= 2; isv(m) = 4; end
end

% % varqnl = varq - epsy.^2;
% % for m=1:nx-1
% % if varqnl/varq <= 0.5; isv(m) = 4; end;
% % end


%  set the criteria for saying a mode is present;  
%  modecrit = 0.2 means 80% var accounted for (kind of low?)
modecrit = 0.2    

% Sverdrup balance?
varqn = varq - (curltau.^2 + trans.^2);
varqn = varqn./varq;
for m=1:nx
    if varqn(m) <= modecrit; isv(m) = 1; end
end

% wbc balance?
varqn = varq - (drag.^2 + trans.^2);
varqn = varqn./varq;
for m=1:nx
    if varqn(m) <= modecrit; isv(m) = 2; end
end

% zonal bc balance?
varqn = varq - (drag.^2 + curltau.^2);
varqn = varqn./varq;
for m=1:nx
    if varqn(m) <= modecrit; isv(m) = 3; end
end



% varqn = varq - (drag.^2 + curltau.^2 + trans.^2);
% varqn = varqn./varq;
% for m=1:nx
%     if varqn(m) >= 0.4; isv(m) = 4; end
% end
% 

figure(81) 
hold on

radb = 40.
axis([-3600 3600 -3600 3600])

xlabel('east distance, km')
ylabel('north distance, km')

if ysample == 1800. | ysample == -1800.
 plot( [-3600. 3600.], [-1800. -1800], 'r')
 plot( [-3600. 3600.], [1800. 1800], 'r')  
else
    
for m=3:5:nx-2
    
if isv(m) == 1; filledCircle( [x(m), ysample], radb, 500, 'g'); end
if isv(m) == 2; filledCircle( [x(m), ysample], radb, 500, 'k'); end
if isv(m) == 3; filledCircle( [x(m), ysample], radb, 500, 'b'); end
if isv(m) == 4; filledCircle( [x(m), ysample], radb, 500, 'r'); end

end

end   %   the if on ysample = +-1800


end   %  loop on j (pick the y to sample
box on
title('modes of potential vorticity balance')

%  add the zonal bl scale wdith
%  note that these are linear estimates, and ignore the very large 
%  thickness changes in the subpolar gyre
zbllin = (1/1.e3)*(r*(L - x)*1.e3/beta).^0.5;
zblnorth = (1/1.e3)*(r*(H/250.)*(L - x)*1.e3/beta).^0.5;
plot(x, (-3600 + zbllin), 'linewidth', 1.6)
plot(x, (3600 - zbllin), 'linewidth', 1.6)
help plot



%  abnalyze the y dependence of Sverdrup balance
trans9 = 0*y; svertrans = 0*y;
for j=2:ny-1
ny8 = j;
nwb = 20;
trans8 = 0.;
for k=nx:-1:nwb 
    v6(k) = h3(ny8,k); u6(k) = h2(ny8,k); h6(k) = h1(ny8,k);
    trans8 = trans8 + h3(ny8,k)*(h1(ny8,k)+H)*dx;
end
trans9(j) = trans8/1.e6;
curltau8 = -(taux(ny8+1) - taux(ny8-1))/(beta*2*dy*1030);
svertrans(j) = (1/1.e6)*curltau8*2*L*1000*(nx-nwb)/nx;
end

svertrans(1) = svertrans(2); svertrans(ny) = svertrans(ny-1);
figure(39)
clf reset
hp = plot(trans9, y, svertrans, y, 'r')
set(hp, 'linewidth', 1.7)
hold on
plot([0 0], [-4000 4000], 'k', 'linewidth', 1.6)
xlabel('volume transport, 10^6 m^3 sec^{-1}')
ylabel(' north,y,  km')
grid on
legend('meridional transport', 'Sverdrup transport', 'location', 'northwest')
axis([-40 40 -3600 3600])


%  plot the pv

% pv0 = f(nyc)/H
% for j=1:nx
%     for k=1:ny        
%         pv(k,j) = (f(k)/(h1(k,j) + H))/pv0;
%         pvlin(k,j) = (f(k)/H)/pv0;
%     end
% % end
% % figure(70)
% clf reset
% [c, hc1] = contour(x,y,pvlin,[-2:0.2:2])
% clabel(c, hc1, 'labelspacing', 500)
% title('linear pv')

% figure(71)
% clf reset
% clear c hc1
% [c, hc1] = contour(x,y,pv,[-4:0.2:8])
% clabel(c, hc1, 'labelspacing', 500)
% title('pv')
% 
% u = h2; v = h3;
% gradq = u*0.; gradqlin = u*0;
% 
% for j=2:ny-1
% for k=2:nx-1    
%     gradq(j,k) = pv0*v(j,k)*(pv(j+1,k) - pv(j-1,k))/(2*dy) ...
%                  +  u(j,k)*(pv(j,k+1) - pv(j,k+1))/(2*dx); 
%     gradqlin(j,k) = pv0*v(j,k)*(pvlin(j+1,k) - pvlin(j-1,k))/(2*dy) ...
%                  +  u(j,k)*(pvlin(j,k+1) - pvlin(j,k+1))/(2*dx);            
%            
% end
% end
% 
% Sv = pi*tau0/(L*1030*1000);   %  scale for transport per unit width
% qscale = Sv/(H^2)
% 
% 
% figure(72)
% clf reset
% [c, hc1] = contour(x,y,gradqlin/qscale, [-10 -5 -2:0.2:2.0 5 10])
% clabel(c, hc1, 'labelspacing', 500)
% caxis([-1 1])
% title('V dot grad linear  pv ')
% 
% figure(73)
% clf reset
% clear c hcl
% [c, hc1] = contour(x,y,gradq/qscale,  [-10 -5 -2:0.2:2.0 5 10])
% clabel(c, hc1, 'labelspacing', 500)
% title('V dot grad pv ')
% 
% 
% 
% figure(28)
% hold on
% subplot(2,1,1)
% eta9 = h3(100,:);
% meaneta = mean(mean(h3));
% eta9 = eta9 - meaneta;
% plot(x, h3(100,:)/(1030/delr))
% subplot(2,1,2)
% hold on
% ht = h1(100,:) + H;
% plot(x, ht)

ct = cumsum(h3(100,:).*ht)*dx;
figure(21)
plot(x, ct, '.')
sver = 2000000*(-0.1/1.e6)/(1000.*beta);



% figure(99)
% clf reset
% zbl = (1/1.e3)*(r*(L - x)*1.e3/beta).^0.;
% axis([-3600 3600 -3600 3600])
% plot(x, (-3600 + zbl))
% hold on
% plot(x, (3600 - zbl))
% axis([-3600 3600 -3600 3600]) 

Contact us