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.

geoadj_2d_plot.m
% geoadj_2d_plot.m
%
% A matlab plotting program to read and plot disk files
% generated by geoadj_2d.for 
%
% 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 ..... 2013.
%

clear all
close all
format compact

%  Set some default graphics things.

set(0,'DefaultTextFontSize',12)    %  these had been 14
set(0,'DefaultAxesFontSize',11)
set(0,'DefaultAxesLineWidth',1.5)
set(0,'DefaultLineLineWidth',1.0)

westend = 1500.;
xscale = 1.0;

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

if readnew == 1
%  come here to 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);
y1 = what(9); 
f0 = what(11);
beta = what(12);
H = what(15);
delh = what(16)

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

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; h5 = h1; h6 = 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... suggest that you save this workspace')
load ga_float.dat;            %  float tracks

%  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 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;

ntime = length(ga_time)

%  unpack the float data
[nft, nf] = size(ga_float)
nf = (nf-1)/2

for m = 1:nft
    tf(m) = ga_float(m,1);
    for n=1:nf
    xf(m,n) = ga_float(m, (n+1))/1000.;
    yf(m,n) = ga_float(m, (nf+1+n))/1000.;
    end
end


% the default is to run the Fortran model to 371 days

 nt9 = ntime;    %  set nt9 = ntime-1 you want a movie to stop at 365 days

 for i=1:nt9  %  this loop is on 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
 
 if i == 1; 
     etamax = max(max(h1));
     etamin = min(min(h1));
     if abs(etamin) >= etamax;
         etamax = etamin;
     end  
 end

%   using the full data set, compute the center of mass
sumhx = sum(h1);
grandsumhx = sum(sumhx);
xcm(i) = sum(sumhx'.*x)/grandsumhx;
sumhy = sum(h1,2);
grandsumhy = sum(sumhy);
ycm(i) = sum(sumhy.*y)/grandsumhy;


%   now plot the interface displacement, eta 
hf22 = figure(35);
clf reset
set(hf22,'Renderer','painters')
set(hf22,'paperposition', [0 0 680 500])

dx1 = 20.;
dy1 = 20.;
xp = -1500:dx1:500;
yp = -1000:dy1:1000;
xp = xp*xscale;
yp = yp*xscale;
[xp3, yp3] = meshgrid(xp, yp);
hp9 = interp2(x,y,h1,xp3,yp3);

caxis = ([-0.2 1.0])

colm = hp9/delh;
colm(colm<=-0.1) = -0.1;
colm(colm>=0.7) = 0.7;

mesh(xp3, yp3, hp9/etamax, colm);
view(-30., 50.);

axis([-1500*xscale xscale*500 -1000*xscale 1000*xscale -0.2 1.0 ]);
  
 xlabel('east, km ');  
 ylabel('north, km ');
 zlabel(['\eta/\eta_o']) % , \eta_{max} =' num2str(etamax,3), ' m']); 
 text(-1200.*xscale, 1000*xscale, 0.8,['\eta_o = ', num2str(delh,3), ' m'])

 hold on
 
 if i == 1
 for jr =1:1000
     thet = 2*pi*jr/1000;
     xr0(jr) = cos(thet);
     yr0(jr) = sin(thet);
     zr4(jr) = 0.2;
 end
 end
 
 r0 = 100.;       %  initial eddy radius, km
 rr = r0 + ga_time(i)*8.64e4*C/1000.;
 xr4 = xr0*rr;
 yr4 = yr0*rr;
 
 xr4(xr4>=500) = NaN;
 xr4(xr4<=-1500) = NaN;
 yr4(yr4>=1000) = NaN;
 yr4(yr4<=-1000) = NaN;
 
 plot3(xr4, yr4, zr4, 'r')
 
 xrbeta = ga_time(i)*8.64e4*Clong/1000.;
 plot3([xrbeta; xrbeta], [-1000; 1000], [0.1; 0.1], 'r')
 

 % add the time
 tim = ga_time(i);
 tim = round(100*tim)/100;
 if tim >= 4.55; tim = floor(tim) + 1; end
 text(-400*xscale, 600*xscale, 0.9, ['time = ',num2str(tim,3), ' days '])                            
 drawnow
 
 
  
  nf1 = 680
  nf2 = 500
  nnh = nnh + 1;
  Hm(nnh) = getframe(hf22, [0 0 nf1 nf2]);
   

%  now plot velocity and eta contours
%  now subsample velocity to put the data on a domain better for plotting
xp1 = -1500.;
dx1 = 50.;
xp2 = 500.;
yp1 = -1000.;
dy1 = 50.;
yp2 = 1000.;
xp = xp1:dx1:xp2;
yp = yp1:dy1:yp2;
xp = xp*xscale;
yp = yp*xscale;
[xp3, yp3] = meshgrid(xp, yp);

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

%  now for eta (to add contours to velocity)
dx1 = 20.;
dy1 = 20.;
xph = xp1:dx1:xp2;
yph = yp1:dy1:yp2;
xph = xph*xscale;
yph = yph*xscale;
[xp3, yp3] = meshgrid(xph, yph);

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


%  contours of eta
hold on
caxis = ([-0.1 0.7])
colm = hp9/delh;
colm(colm<=-0.1) = -0.1;
colm(colm>=0.7) = 0.7;

[c, hc] = contourf(xp3, yp3, colm, [-0.15:0.1:0.75]);
C = 3.0;
% 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
hp2(5,3) = 0.5*C*abs(etamax)/H;
hp3(5,3) = 0.;
usc = 100./vmax;   
hq = quiver(xp, yp, usc*hp2, usc*hp3, 0.);
set(hq, 'color', 'w','linewidth', 1.5);
speed9 = 0.5*C*abs(etamax)/H
htt = text(xp(3), yp(3)+225*xscale, [num2str(speed9,2), ' m sec^{-1}']),
set(htt, 'color', 'w', 'fontsize', 12)

%  plot the floats

if nf >=1   
    
    dft = abs(tf - ga_time(i));
    itf = find(dft <= 5);
    for jf = 1:nf
        plot(xf(itf,jf), yf(itf,jf), '.y', 'markersize', 2);
        tend = max(itf);
        plot(xf(tend,jf), yf(tend, jf), 'oy', ...
                              'markeredgecolor', 'y', ...
                               'markerfacecolor', 'r', ...
                           'markersize', 3) 
                       if jf >= 10       
                           plot(xf(tend,jf), yf(tend, jf), 'oy', ...
                              'markeredgecolor', 'y', ...
                               'markerfacecolor', 'g', ...
                           'markersize', 3) 
                       end                        
    end
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', 'w'); 

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


% 
% clear a1 a2;
% minh = min(min(h1));
% [a1] = find(h1 == minh);
% xx = xm(a1); yy = ym(a1);   
% a2 = a1(1);
% 
% if del <= 0
% if minh <= -0.5; 
%     ht = text(xx, yy,'L');
%     set(ht, 'fontsize', 12, 'fontweight', 'bold', 'color', 'w'); 
%     xpeak(i) = xx;
%     ypeak(i) = yy;
% end
% 
% 

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

  axis([-1500*xscale 500*xscale -1000*xscale 1000*xscale])
  axis('square')
 tim = ga_time(i);
  tim = round(100*tim)/100;
 if tim >= 4.55; tim = floor(tim) + 1; end
 text(-1200., xscale*750., ['time = ',num2str(tim,3), ' days '], 'color', 'w');

 nm = nm + 1;
 Um(nm) = getframe(hf1, [0 0 nf1 nf2]); 
 pause(0.2)
 

% % add the equator
% 
%  xx = [0. x(nx)];
%  yd2 = y(ny)/2.;
%  yy = [yd2 yd2];
%  hold on
% % line(xx,yy)

% if i == 6
% print -depsc eta_cover.eps     %  save the last frame
% end

if i == 2
    hq = h1(2,2);
end
h1 = h1 - h2;

%   now plot the passive tracer ...............
plottr = 1;    %  or not, if this flag set to 0
if plottr == 1
hf20 = figure(20);
clf reset
set(gcf,'paperposition', [0 0 500 500])
 
dx1 = 5.;
dy1 = 5.;
xp = -westend:dx1:xp2;
yp = yp1:dy1:yp2;
xp = 1.5*xp*xscale;
yp = yp*xscale;
[xp3, yp3] = meshgrid(xp, yp);

hp = interp2(x,y,h5,xp3,yp3);

[cc ch] = contourf(xp3, yp3, hp, [-.05:0.05:1.0]);
caxis = ([-0.05 1.05])
set(ch,'color', 'none')  
hcb = colorbar;
set(hcb, 'position', [0.91 0.1800 0.03 0.60])
title('a passive tracer')

 xlabel('east, km ');  
 ylabel('north, km ');
 zlabel('tracer');  
 
 % add the time
 tim = ga_time(i);
 tim = round(100*tim)/100;
 if tim >= 4.55;  tim = floor(tim) + 1; end
ht = text(-1200*xscale, xscale*700., 0.8, ['time = ',num2str(tim,3), ' days ']) 
set(ht,'color', 'w')
% axis([-1500*xscale 500*xscale -1000*xscale 1000*xscale -0.5 1.5])
axis([-1500*xscale 500*xscale -1000*xscale 1000*xscale])
axis('square')

%  add a circle that expands at the gravity wave speed
r0 = 1000.;
for p=1:1000
    th = 2*pi*0.001*p;
    xr(p) = r0*cos(th);
    yr(p) = r0*sin(th);   
end
hold on
% plot(xr, yr, 'w--', 'linewidth', 1.8)

nf1 = 418;
nf2 = 418;  
Tm(nnh) = getframe(hf20, [0 0 nf1 nf2]);
 
pause(0.1)
end

 
%   now plot the potential vorticity, q 

qplot = 1
if qplot == 1
hf40 = figure(40);
clf reset
set(gcf,'paperposition', [0 0 500 500])
 
dx1 = 20.;
dy1 = 20.;
xp = -westend:dx1:xp2;
yp = yp1:dy1:yp2;
[xp3, yp3] = meshgrid(xp, yp);
hpq = interp2(x,y,h6,xp3,yp3);

%  shall it be the actual q or the q difference from f/H0?
qprime = 0    % =0 for the full q;   =1 for the deviation from f/H
if qprime == 1   
if i == 1;
    q1 = 0*hpq;
end    
if i == 2
    [nny nnx] = size(hpq);
    for mm=1:nny
    for nn=1:nnx
    q1(mm,nn) = hpq(mm,1);
    end
    end
end
hpq = hpq - q1;
end   %  if on qprime

% mesh(xp3, yp3, hp/(f0/H));

hpq = hpq/(f0/H);
clear cc ch qlv
qlv = [0.7 0.85:0.001:1.15];  
caxis = ([0.85 1.15]);
if qprime == 1
    qlv = ([-0.1:0.01:0.1]);
    caxis = ([-0.1 0.1]);
end
[cc ch] = contourf(xp3, yp3, hpq, qlv );
set(ch,'color', 'none')  
hcb = colorbar;
set(hcb, 'position', [0.91 0.1800 0.03 0.60])

title('potential vorticity/(f_o/H)')
 xlabel('east, km ');  
 ylabel('north, km ');
 zlabel('q, mks');  
 

 % add the time
 tim = ga_time(i);
 tim = round(100*tim)/100;
 if tim >= 4.55;  tim = floor(tim) + 1; end
ht = text(-1200*xscale, xscale*700., 0.8, ['time = ',num2str(tim,3), ' days ']) 
set(ht,'color', 'w') 
 
 axis('square')

  nf1 = 418;
  nf2 = 418;  
  Qm(nnh) = getframe(hf40, [0 0 nf1 nf2]);
 
 pause(0.1)
end    %  if on plot tracer and q 


end    %  end of the first time step loop
%
%

%  compute the speed of the peak and the Clong for the average latitude


figure(80)
clf reset
plot(ga_time, xpeak, 'b', ga_time, ypeak, 'r')
xlabel('time, days')
ylabel('xpeak, ypeak, km')
grid on

meany = 0.5*min(ypeak)
Cfactor = (1 - 2*1000*beta*meany/mean(f))

meanxpeak = 1000.*min(xpeak)/(max(ga_time)*8.64e4)
meanspeednorm = meanxpeak/(Clong*Cfactor)




%  now cycle throuogh the data time-wsie again for q-conservation analysis 
% (this is somewhat wasteful, but....)
%

for jy = 1:10   %  loop through the data set up to 10 times to take  
                %   different y slices
5 
%   yline = 0 is the center line of the initial eddy 
%   yline = -90 is the eddy peak for eta0 = 50 m, and
%   yline = -200 is the eddy q-center in the etao = 50 m experiment

kc = 0;
yline = input...
 ('  enter the y (km) for q-conservation;  very large y to exit out  ') 

 ya = abs(y - yline);
 [yamin, nyc] = min(ya);
 nyc
    
 for i=1:ntime  %  this is a loop on 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
     
     for la9=1:nx
     dqdx(la9) = 0.;
     dqdy(la9) = 0.;
     end
     for lx=3:nx-2
         dqdx(lx) = (h6(nyc,lx+1) - h6(nyc,lx-1))/(2*dx);
         dqdy(lx) = (h6(nyc+1,lx) - h6(nyc-1,lx))/(2*dy);
     end
     
     kc = kc + 1;
     
     etax(kc,:) = h1(nyc,:);
     vorx(kc,:) = vor(nyc,:);
     betav(kc,:) = h3(nyc,:); 
     advq(kc,:) =(h1(nyc,:) + H).*(h2(nyc,:).*dqdx + h3(nyc,:).*dqdy);
     advqx(kc,:) = (h1(nyc,:) + H).*(h2(nyc,:).*dqdx);
     advqy(kc,:) = (h1(nyc,:) + H).*(h3(nyc,:).*dqdy);
end    %  time step loop ends here

%  evaluate the terms in linear q-conservation written as
%     relvort + stretch = beta
%       dvort +   deta  = betav5
   
nct = ntime - 1;
dt7 = 8.64e4*(ga_time(nct+1) - ga_time(nct-1));
eta0 = delh;
bu = beta*(C*eta0/H);    %  the q scale

dvor = (vorx(nct+1,:) - vorx(nct-1,:))/dt7;
betav5 = -beta*betav(nct,:);
deta = -(f0/H)*(etax(nct+1,:) - etax(nct-1,:))/dt7;
advq9 = -advq(nct,:);
advq9x = -advqx(nct,:);     %  the east-west part
advq9y = -advqy(nct,:);     %  north-south part
detanl =  -(f0 + vorx(nct,:)).*(etax(nct+1,:) - etax(nct-1,:))./...
    (dt7*(H + etax(nct,:)));

 eps1 = dvor + deta - betav5;  %  should be very small if the balance is linear

set(0,'DefaultTextFontSize',16)
set(0,'DefaultAxesFontSize',16)
figure(26+jy)
clf reset

epsbv =   dvor + deta;
hp = plot(x, dvor/bu, x, deta/bu, x, betav5/bu, x, epsbv/bu, 'k--' );
hold on
set(hp,'linewidth', 1.8)
legend('relative +', 'stretching', ' = beta', 'rel + str')
legend boxoff
hold on
plot(x, 0*x, 'k')
xlabel('east, km')
ylabel('\partial vorticity / \partial t/(\beta C \eta_0/H)')
 axis([-1500*xscale 500*xscale -0.1 0.15])
text(-1300, 0.1,  ['365 days'])
grid on
title(['q conservation along y = ', num2str(yline), ' km'])

%  plot the nonlinear terms too
figure(90+jy)
clf reset
hold on
plot(x, (advq9 - betav5)/bu, 'm', x, eps1/bu, 'k')
hold on
plot(x, advq9x/bu, 'r', x, (advq9y - betav5)/bu, 'g')
legend('adv - beta', 'linear misfit', 'east term', 'north term')
xlabel(' east, km')
ylabel('\partial vorticity / \partial t/(\beta C \eta_0/H)')
axis([-1500*xscale 500*xscale -0.2 0.4])
title(['q conservation along y = ', num2str(yline), ' km'])

end   %  loop on jy 
% all done with q


%  save some movies?

imov = input('  Shall we save avi movie files?  (1 = yes)')
if imov == 1 
    movie2avi(Hm, 'ga2d-h.avi')
    movie2avi(Um, 'ga2d-u.avi')
    movie2avi(Tm, 'ga2d-q.avi') 
    movie2avi(Qm, 'ga2d-q.avi')
end

Contact us