Global Index (short | long) | Local contents | Local Index (short | long)
Note: equatorial rossby radius of deformation is given by:
R_e = sqrt(c/(2*beta)); Note again, however, that the wave
speed 'c' should be given by the equivelent depth, rather than
some depth of the ocean for a barotropic wave. (Gill, p437).
| This script calls | |
|---|---|
clear
global RADUS
% Parameters
% Continuous
y = [-20:.1:20]'; % = y/re
x = [-20:.1:20]; % = y/re
wn = 2; % Wavenumber
c = 2; % ==> He ~= 1m, or a 150m 2-layer system with g' = 0.06
% Also, this is Kelvin wave speed.
%c = sqrt(9.8 * 150);
beta = 2*7.292e-5./RADUS;
k = -2*pi*wn/(2*pi*RADUS); % Wavenumber 4
re = sqrt(c/(2*beta)) % re = 2.5606e+05
eta = zeros(length(y), 4);
u = zeros(length(y), 4);
v = zeros(length(y), 4);
% Start with Kelvin wave
u(:,1) = (9.8/c)*exp(-0.25 * y.^2);
eta(:,1) = exp(-0.25 * y.^2);
v(:,1) = zeros(length(y), 1);
% Planetary waves
y2 = sqrt(1)*y;
c = 1*c;
for n = 1:2:5; % Start with first mode
omega = -1*beta*k ./ (k^2 + (2*n+1)*beta/c);
hn = hermite((1/sqrt(2))*y2, n); % 2*y; % Hermite Polynomial (from table)
hnp1 = hermite((1/sqrt(2))*y2, n+1); % 4*y.^2-2;
hnm1 = hermite((1/sqrt(2))*y2, n-1); % 1;
v(:,(n+1)/2+1) = 2.^(-0.5*n).*hn.*exp(-0.25*y2.^2)*2^-0.5;
q = (c * k - omega)^-1 * (2*beta*c)^0.5 * 2^(-n/2 - 0.5) * ...
hnp1 .* exp(-0.25*y2.^2) * 2^-0.5;
r = (c * k + omega)^-1 * (2*beta*c)^0.5 * 2^(-n/2 + 0.5) * n * ...
hnm1 .* exp(-0.25*y2.^2) * 2^-0.5;
u(:,(n+1)/2+1) = 0.5*(q-r);
eta(:,(n+1)/2+1) = (c./(2*9.8))*(q+r);
end;
scale = max(abs(eta));
eta2 = eta ./ (ones(length(y), 1) * scale);
u2 = u ./ (ones(length(y), 1) * scale);
v2 = v ./ (ones(length(y), 1) * scale);
figure(1); fo(1);
subplot(3,1,1);
h1 = plot(...
y, eta2(:,1), '-k', ...
y, eta2(:,2), '--k', ...
y, eta2(:,3), '-.k', ...
y, eta2(:,4), '.:k')
set(h1(1), 'linewidth', 2);
axis([-8 8 -1.25 1.25]);
grid on;
title('Normalized Sea Level Anomaly');
legend(h1, 'K', 'R1', 'R3', 'R5');
subplot(3,1,2);
h1 = plot(...
y, u2(:,1), '-k', ...
y, u2(:,2), '--k', ...
y, u2(:,3), '-.k', ...
y, u2(:,4), '.:k')
set(h1(1), 'linewidth', 2);
axis([-8 8 -6 8]);
grid on;
title('Scaled Zonal Velocity');
legend(h1, 'K', 'R1', 'R3', 'R5');
subplot(3,1,3);
h1 = plot(...
y, v2(:,1), '-k', ...
y, v2(:,2), '--k', ...
y, v2(:,3), '-.k', ...
y, v2(:,4), '.:k')
set(h1(1), 'linewidth', 2);
axis([-8 8 -.5 .5]);
set(gca, 'YTick', [-.5:.25:.5]);
grid on;
title('Scaled Meridional Velocity');
% xlabel('Meridional Distance (Scaled)')
legend(h1, 'K', 'R1', 'R3', 'R5');
xlabel('Y / R_e; R_e = 209km = 1.9 deg');
cd ~/Thesis/Talk
%print -dps2 Rossby_Kelvin_discrete.ps
% Look at model solution
clear
global RADUS
% Parameters
for biff = 1:2;
% Model 'y'
if biff == 1;
lat = getll('temp', [0 20 -90 90]);
else
lat = getll('u', [0 20 -90 90]);
end
wn = 4; % Wavenumber
c = 2; % ==> He ~= 1m, or a 150m 2-layer system with g' = 0.06
% Also, this is Kelvin wave speed.
%c = sqrt(9.8 * 150);
beta = 2*7.292e-5./RADUS;
k = -2*pi*wn/(2*pi*RADUS); % Wavenumber
re = sqrt(c/(2*beta)) % re = 2.5606e+05
% Model 'y'
y = lat*pi/180*RADUS/re;
if biff == 1;
eta = zeros(length(y), 4);
else
u = zeros(length(y), 4);
v = zeros(length(y), 4);
end
% Start with Kelvin wave
if biff == 1;
eta(:,1) = exp(-0.25 * y.^2);
else
u(:,1) = (9.8/c)*exp(-0.25 * y.^2);
v(:,1) = zeros(length(y), 1);
end
% Planetary waves
y2 = sqrt(1)*y;
c = 1*c;
for n = 1:2:5; % Start with first mode
omega = -1*beta*k ./ (k^2 + (2*n+1)*beta/c);
hn = hermite((1/sqrt(2))*y2, n); % 2*y; % Hermite Polynomial (from table)
hnp1 = hermite((1/sqrt(2))*y2, n+1); % 4*y.^2-2;
hnm1 = hermite((1/sqrt(2))*y2, n-1); % 1;
if biff ~= 1;
v(:,(n+1)/2+1) = 2.^(-0.5*n).*hn.*exp(-0.25*y2.^2)*2^-0.5;
end
q = (c * k - omega)^-1 * (2*beta*c)^0.5 * 2^(-n/2 - 0.5) * ...
hnp1 .* exp(-0.25*y2.^2) * 2^-0.5;
r = (c * k + omega)^-1 * (2*beta*c)^0.5 * 2^(-n/2 + 0.5) * n * ...
hnm1 .* exp(-0.25*y2.^2) * 2^-0.5;
if biff ~= 1;
u(:,(n+1)/2+1) = 0.5*(q-r);
else
eta(:,(n+1)/2+1) = (c./(2*9.8))*(q+r);
end
end;
if biff == 1;
scale = max(abs(eta));
eta2 = eta ./ (ones(length(y), 1) * scale);
else
u2 = u ./ (ones(length(y), 1) * scale);
v2 = v ./ (ones(length(y), 1) * scale);
end
figure(2); fo(1);
if biff == 1;
subplot(3,1,1);
h1 = plot(...
y, eta2(:,1), '-k', ...
y, eta2(:,2), '--k', ...
y, eta2(:,3), '-.k', ...
y, eta2(:,4), '.:k')
set(h1(1), 'linewidth', 2);
axis([-8 8 -1.25 1.25]);
grid on;
title('Normalized Sea Level Anomaly');
legend(h1, 'K', 'R1', 'R3', 'R5');
else
subplot(3,1,2);
h1 = plot(...
y, u2(:,1), '-k', ...
y, u2(:,2), '--k', ...
y, u2(:,3), '-.k', ...
y, u2(:,4), '.:k')
set(h1(1), 'linewidth', 2);
axis([-8 8 -4 6]);
grid on;
title('Scaled Zonal Velocity');
legend(h1, 'K', 'R1', 'R3', 'R5');
subplot(3,1,3);
h1 = plot(...
y, v2(:,1), '-k', ...
y, v2(:,2), '--k', ...
y, v2(:,3), '-.k', ...
y, v2(:,4), '.:k')
set(h1(1), 'linewidth', 2);
axis([-8 8 -.75 .75]);
set(gca, 'YTick', [-.5:.25:.5]);
grid on;
title('Scaled Meridional Velocity');
% xlabel('Meridional Distance (Scaled)')
legend(h1, 'K', 'R1', 'R3', 'R5');
xlabel('Y / R_e; R_e = 256km = 2.3 deg');
end
end; % biff loop
cd ~/Thesis/Talk
%print -dps2 Rossby_Kelvin_CSIRO.ps