Code covered by the BSD License

a Coriolis tutorial

James Price (view profile)

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

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

disp('  the determinant: ')

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.