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
cint = [-10 -15 1 1]
fnum = 4;
global RADUS
% Parameters
% Continuous
y = [-15:.2:15]'; % = y/re
x = [-180:2:180]; % = y/re
wn = 1; % 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
x = RADUS*pi/180*x/re;
eta = zeros(fnum, length(y), length(x));
u = zeros(fnum, length(y), length(x));
v = zeros(fnum, length(y), length(x));
% Start with Kelvin wave
u(1,:,:) = (9.8/c)*exp(-0.25 * y.^2)*cos(re*k*x);
eta(1,:,:) = exp(-0.25 * y.^2)*cos(re*k*x);
v(1,:,:) = zeros(1, length(y), length(x));
% 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;
temv = 2.^(-0.5*n).*hn.*exp(-0.25*y2.^2)*2^-0.5;
v((n+1)/2+1,:,:) = temv*cos(k*re*x);
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;
q = q*sin(re*k*x);
r = r*sin(re*k*x);
u((n+1)/2+1,:,:) = 0.5*(q-r);
eta((n+1)/2+1,:,:) = (c./(2*9.8))*(q+r);
end;
for i = 1:fnum;
scale = max(max(abs(eta(i,:,:))));
eta(i,:,:) = eta(i,:,:)./scale;
u(i,:,:) = u(i,:,:)./scale;
v(i,:,:) = v(i,:,:)./scale;
end
get_global;
XAX = re*180/(pi*RADUS)*x'; YAX = y*re*180/(RADUS*pi);
%FRAME = [min(XAX) max(XAX) -12 12];
FRAME = [-180 180 -20 20];
[xk, yk] = keep_var(FRAME, XAX, YAX);
XAX = XAX(xk); YAX = YAX(yk);
for i = 1:3;
figure(1); fo(1);
sptalk(fnum,1,i);
% [tem2, lat2, lon2] = sph_div1(squeeze(u(i,:,:)), ...
% squeeze(v(i,:,:)), YAX, XAX, 1);
tem = eta(i,yk,xk);
gcont(tem, .2);
% hold on;
% pncont(lon2, lat2, -1e6*tem2, (3/i)*[-.5:.05:-.05 .05:.05:.5], 0, 'r');
temu = u(i,yk,xk); temv = v(i,yk,xk);
hold on;
gquiv(temu, temv, cint(i), 5, ' ');
hold off;
axis([-180 180 -10 10]);
set(gca, 'XTick', [-180:45:180], 'YTick', -10:5:10, 'box', 'on');
grid off
if i == 1;
ylabel('Kelvin');
elseif i == 2;
ylabel('Rossby N=1')
elseif i == 3;
ylabel('Rossby N=3')
end
xlabel('');
end
cd /home/disk/tao/dvimont/Thesis/Wave
print -dps2 cont_k_r1_r2.ps
% Look at model solution
global RADUS
% Parameters
for biff = 1:2;
% Model 'y'
if biff == 1;
[lat, lon] = getll('temp', [-0.1 360 -30 30]);
else
[lat, lon] = getll('u', [-0.1 360 -30 30]);
end
x = [-180:10:180]; % = y/re
wn = 1; % 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
%x = RADUS*pi/180*(lon-180)'/re;
x = RADUS*pi/180*x/re;
% Model 'y'
y = lat*pi/180*RADUS/re;
if biff == 1;
eta = zeros(fnum, length(y), length(x));
else
u = zeros(fnum, length(y), length(x));
v = zeros(fnum, length(y), length(x));
end
% Start with Kelvin wave
if biff == 1;
eta(1,:,:) = exp(-0.25 * y.^2)*cos(re*k*x);
else
u(1,:,:) = (9.8/c)*exp(-0.25 * y.^2)*cos(re*k*x);
v(1,:,:) = zeros(1, length(y), length(x));
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;
temv = 2.^(-0.5*n).*hn.*exp(-0.25*y2.^2)*2^-0.5;
v((n+1)/2+1,:,:) = temv*cos(k*re*x);
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;
q = q*sin(re*k*x);
r = r*sin(re*k*x);
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;
for i = 1:fnum;
scale(i) = max(max(abs(eta(i,:,:))));
eta(i,:,:) = eta(i,:,:) ./ scale(i);
end
else
for i = 1:fnum
u(i,:,:) = u(i,:,:) ./ scale(i);
v(i,:,:) = v(i,:,:) ./ scale(i);
end
end
get_global;
XAX = re*180/(pi*RADUS)*x'; YAX = y*re*180/(RADUS*pi);
%FRAME = [min(XAX) max(XAX) -12 12];
FRAME = [-180 180 -30 30];
for i = 1:3;
figure(2); fo(1);
sptalk(fnum,1,i);
tem = eta(i,:,:);
if biff == 1;
gcont(tem, .2);
else
temu = u(i,:,:); temv = v(i,:,:);
hold on;
gquiv(temu, temv, cint(i), 1, ' ');
hold off;
end
axis([-180 180 -10 10]);
set(gca, 'XTick', [-180:45:180], 'YTick', -10:5:10, 'box', 'on');
grid off
if i == 1;
ylabel('Kelvin');
elseif i == 2;
ylabel('Rossby N=1')
elseif i == 3;
ylabel('Rossby N=3')
end
xlabel('');
end
end; % biff loop
figure(2);
cd /home/disk/tao/dvimont/Thesis/Wave
print -dps2 csiro_k_r1_r2.ps
figure(1);
sptalk(2,1,1);
title('Kelvin Wave');
figure(2);
sptalk(2,1,2);
title('First order Rossby Wave');
cd ~/Thesis/Talk
%figure(1); print -dps2 Kelvin_wave.ps
%figure(2); print -dps2 Rossby1_wave.ps