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.

twolayer_eig.m

%  some symbolic math for modal calculations
%  one, two and three layers.....
%
%  Jim Price, Nov. 2010. 


clear all
close all

set(0,'DefaultTextFontSize',18)
set(0,'DefaultAxesFontSize',18)
set(0,'DefaultAxesLineWidth',1.3)
set(0,'DefaultLineLineWidth',1.7)

format compact

%  start with a one layer model.....
syms A w k H1 H2 rrat g gp w1 Ar f kx ky

A = [-w g*k; H1*k -w];
pretty(simplify(det(A)))

%  with rotation, f = constant;   must assume waves ~ exp(i(kx - wt))

Ar =     [-i*w     -f    i*g*kx; 
           f      -i*w   i*g*ky; 
        i*H1*kx  i*H1*ky  -i*w];
pretty(simplify(det(Ar)))

clear A w k H1 H2 rrat g gp w1 Ar f

g = 9.8*2/1030;
H1 = 500.;
f = 2*7.292e-5*sin(30*pi/180.);
c = sqrt(g*H1)
Rd = c/f

for m=1:100
    lambda = 50.e3 + (m-1)*20000.;
    k(m) = 2*pi/lambda;
    kn = k(m);

w1(m) = sqrt(g*H1)*kn;
wr(m) = sqrt(g*H1*kn^2 + f^2);
wnd(m) = sqrt(kn^2*Rd^2 + 1);
end

figure(20)
clf reset

plot(k*Rd, w1/f, k*Rd, wr/f, 'r', k*Rd, wnd,'g' )
axis([0 5 0 5])
xlabel('k*Rd')
ylabel('freq/f')


%  for j=1:100
%      lambda = j*5.e3;
%      k = 2*pi/lambda;
%      w1val(j) = subs(w1);
%      wcbtr(j) =  k*sqrt(g*H1);
%      kval(j) = k;
%  end
%  
%  figure(1)
%  clf reset
%  
%  plot(kval, w1val, 'b', kval, wcbtr, 'b*')
%  xlabel('k, 1/m')
%  ylabel('freq')
%  grid on
%  
%  now for the two layer model.....
%  w is the frequency,  rrat is rho1/rho2 < 1

clear A g k w H1 H2 rrat
syms g k H1 H2 rrat w H3  rrat1 rrat2 

A = [ -w   g*k       0       g*k; 
      H1*k -w        0        0; 
       0  g*rrat*k  -w       g*k;
       0    0       H2*k     -w ]
   
 clear V D 
 [V,D]  = eig(A);
 disp('  the two layer eigenvectors...')
 pretty(simplify(V))
 disp('  the two layer eigenvalues.....')
 pretty(simplify(D))
 
  %  oceanic startification
 H1 = 500.;
 H2 = 3500.;
 rrat = 1025/1027;
 
 H1bcl = (1/(2*H2*rrat))*(H1 - (H1^2 + H2^2 - 2*H1*H2 + 4*H1*H2*rrat)^(1/2) -H2)
 
 g = 9.8;
 %  oceanic startification
 H1 = 500.;
 H2 = 3500.;
 rrat = 1025/1027;
 %  have to specify some k, w
 Pdays = 1.
 w = 1./(Pdays*8.64e4); 
 lambda = 100.e3
 k = 2*pi/lambda
 
 disp('  the numerical two layer eigen vectors')
 Vval = subs(V)
 disp('  the numerical two layer eigen values....')
 Dval = subs(D)

 
 
%  compute the  characteristic equation 
    
Adet = det(A);
disp('  the determinant: ')
pretty(simplify(Adet))



syms  a b c w1 w2
 
g = 9.8
 %  oceanic stratification
 H1 = 500.;
 H2 = 3500.;
 rrat = 1025/1027;

%  just type these in.......


 w1 = (-b + sqrt(b^2 - 4*a*c))/(2*a)   %  the barotropic mode
 % w1lim = limit(w1,rrat,1,'left')
 
 w2 =  (-b - sqrt(b^2 - 4*a*c))/(2*a)  %  the baroclinic mode
 % w2lim = limit(w2, H2, inf)
 

 
 for m=1:1   
 
 
 %  a rough guess at atmospheric startification 
 %  the baroclinic wave seems reasonable;  the barotropic mode,  doubtful.
 if m == 2
 H1 = 8000;
 H2 = 8000;
 rrat = 0.1/0.6;
 end

 w1val = zeros(1, 200);
 w2val = w1val; 
 wcbcl = w1val;
 wcbcl2 = w1val;
 wcbtr = w1val; 
 kval = w1val;
 
 for j=1:200
     lambda = (j-1)*200.e3 + 10.;
     k = 2*pi/lambda;
     
 a = 1;
 b =  -k^2*g*(H1 + H2); 
 c = (1 - rrat)*H1*H2*g^2*k^4;
     
     w1val(j) = sqrt(subs(w1));
     w2val(j) = sqrt(subs(w2));
     wcbtr(j) = k*sqrt(g*(H1 + H2));
     wcbcl1(j) = k*sqrt(g*(1 - rrat)*H1); 
     wcbcl2(j) = k*sqrt(g*(1 - rrat)*H1*H2/(H1 + H2));
     kval(j) = k;
 end
 
 kval = kval*1000;     %  convert to km
 w1val = w1val*3600;   %  hours
 w2val = w2val*3600;
 wcbtr = wcbtr*3600;
 wcbcl1 = wcbcl1*3600;
 wcbcl2 = wcbcl2*3600;
 
 figure(m)
 clf reset
 hold on
 
 set(gcf, 'papersize', [7 3.5])
 plot(kval, w1val, 'b') % ,  kval, wcbtr, 'b*')
 plot(-kval, w1val, 'b')
 hold on
 plot(kval, w2val, 'r', kval, wcbcl1, 'r--') % , kval, wcbcl2, 'ro')
  plot(-kval, w2val, 'r', kval, wcbcl1, 'r--')
 xlabel('wavenumber, k, km^{-1}')
 ylabel('frequency, \omega,  hour^{-1}')
 grid on

 cbtr = w1val(100)/(kval(100)*3.6)    %  C in m/sec
 cbcl = w2val(100)/(kval(100)*3.6)    %  

 if m == 1
 axis([-0.1 0.1 0 10])
 if cbtr >= 100; cbtr = 200;   end;
 if cbcl <= 4; cbcl = 3; end;
 ss = {'barotropic'; ' C = 200 m sec^{-1}'}
 st = {'baroclinic'; ' C = 3 m sec^{-1}'}
 text(-0.08, 7., ss )    
 text(0.02, 1.7, st)
 % title('two-layer oceanic gravity waves')
 plot([0 0 ], [0 100], 'k')
 end
 if m == 2;
 axis([0 0.1 0 100])
 text(0.003, 90., ['barotropic, C = ', num2str(cbtr, 3), ' m sec^{-1}'])    
 text(0.03, 10., ['baroclinic, C = ', num2str(cbcl, 3), ' m sec^{-1}'])
 title('atmospheric gravity waves')
 end
 end

 g = 9.8;
 H1 = 500.;
 H2 = 3500.;
 rrat = 1025/1027;
 Pdays = 1.
 w = 1./(Pdays*8.64e4); 
 lambda = 100.e3
 k = 2*pi/lambda
 Vval = subs(V)
 
 Vval(1,1)*H1
 Vval(3,1)*H2
 
 disp('   the two layer eigen vectors and values....')

 clear V D 
 [V,D]  = eig(A);
 pretty(simplify(V))
 pretty(simplify(D));
 
%  check out Gamm1(rrat)
  H1n = 500.
  H2n = 3500.
  
  
  for j=1:50
      rhor(j) = 1 - 0.0001*(j-1);
      Gamma1(j) = (1/(2*rhor(j)))*...
          ((H1n - ((H1n^2 + H2n^2 - 2*H1n*H2n + 4*H1n*H2n*rhor(j))^0.5))/H2n - 1);
      G7(j) = -1/rhor(j);  
  end
  figure(10)
  clf reset
  plot(rhor, Gamma1, rhor, G7, 'r*')
  xlabel('rhor')
  ylabel('Gamma1')
  grid on
  
 
%  try three layers.....  
%  I can't seem to get the purely symbolic math to work,
%  but this does compute the numerical values that are just as useful.

clear A3 g k w H1 H2 rrat V3val D3val 
syms A3 g k H1 H2 rrat w H3  rrat1 rrat2 

 A3 = [ -w    g*k      0     g*k     0    g*k; 
      H1*k   -w       0      0      0     0; 
       0  g*rrat1*k  -w     g*k     0    g*k;
       0     0       H2*k   -w      0     0;
       0  g*rrat2*k   0  g*rrat1*k  -w   g*k;   
       0     0        0      0     H3*k   -w ]
  
 g = 9.8;
 H1 = 1000.;
 H2 = 1000.;
 H3 = 1000.;
 rrat1 = 0.99;   % rho1/rho2; also rho2/rho3
 rrat2 = 0.98;   % rho1/rho3
 Pdays = 1.
 w = 1./(Pdays*8.64e4); 
 lambda = 100.e3
 k = 2*pi/lambda
   
 A3val = subs(A3);     % the numerical A    
 [V3val,D3val]  = eig(A3val)  % the numerical eigenvectors and values.
 


Contact us