Global Index (short | long) | Local contents | Local Index (short | long)
Define stuff
| This script calls | |
|---|---|
clear
c = 289;
T = 750*24*3600;
dd = 100e5;
Lx = 180*dd/2;
Ly = 102*dd;
global RADUS OMEGA
rearth = RADUS*100;
x = [0:dd/10:(2*Lx)]-Lx;
y = [0:dd/10:Ly]'-Ly/2;
xlim1a = [0 Lx];
xlim1b = [-50 50]*(pi/180)*rearth;
ylim = [-6*pi*rearth/180 6*pi*rearth/180];
% Non-dimensionalize
beta = 2*OMEGA/rearth;
Ldim = sqrt(c/beta);
Tdim = 1./sqrt(c*beta);
lam = Tdim/T;
x2 = x./Ldim;
y2 = y./Ldim;
xlim2a = xlim1a./Ldim;
xlim2b = xlim1b./Ldim;
ylim2 = ylim./Ldim;
% Define windstress
wsx = zeros(1, length(x2));
wsy = zeros(length(y2), 1);
xkp = find(x >= xlim1b(1) & x <= xlim1b(2));
wsx(xkp) = cos(pi*x(xkp)/(100*(pi/180)*rearth));
wsy = exp(-1*beta*(y.^2)/4000.);
taux = wsy*wsx;
%subplot(3,1,1); plot(x, wsx);
%subplot(3,1,2); plot(y, wsy);
%subplot(3,2,5); contour(x(1:20:1801), y(1:20:1021), taux(1:20:1021, 1:20:1801));
% Start with locally forced Kelvin wave
expfun = exp(-0.5*y2.^2);
n = 0;
c0 = c;
pcfun = expfun.*hermite(y2, n);
pcfun = pcfun./sqrt(pcfun'*pcfun);
yk = find(y2 >= ylim2(1) & y2 <= ylim2(2));
kscale = sum(pcfun(yk))./sum(pcfun);
xk2 = find(x > xlim1a(1) & x < xlim1a(2));
xk1 = 2:(xk2(1)-1);
D1 = T*c0/Lx * exp(x(xk1)/(c0*T))*(1-exp(-1*Lx/(c0*T)));
D2 = T*c0./(Lx-x(xk2)) .* (1-exp((x(xk2)-Lx)/(c0*T)));
projk = pcfun'*taux;
q0l = kscale*(sum(D1.*projk(xk1))+sum(D2.*projk(xk2)));
% Reflected Rossby waves
n = 1;
c0 = c;
c2 = c./(2*n+1);
pcfun = expfun.*hermite(y2, n+1);
pcfun = pcfun./sqrt(pcfun'*pcfun);
r2scale = sum(pcfun(yk))./sum(pcfun);
Re1 = sqrt(n/(n+1));
D3 = c2*T/Lx * exp((x-Lx)/(c0*T))*(1-exp(-1*Lx/(c2*T)));
q2e = r2scale*Re1*sum(D3([xk1 xk2]).*projk([xk1 xk2]))
n = 3;
c0 = c;
c4 = c./(2*n+1);
pcfun = expfun.*hermite(y2, n+1);
pcfun = pcfun./sqrt(pcfun'*pcfun);
r4scale = sum(pcfun(yk))./sum(pcfun);
Re3 = Re1*sqrt(n/(n+1));
D3 = c4*T/Lx * exp((x-Lx)/(c0*T))*(1-exp(-1*Lx/(c4*T)));
q4e = r4scale*Re3*sum(D3([xk1 xk2]).*projk([xk1 xk2]))
n = 5;
c0 = c;
c6 = c./(2*n+1);
pcfun = expfun.*hermite(y2, n+1);
pcfun = pcfun./sqrt(pcfun'*pcfun);
r6scale = sum(pcfun(yk))./sum(pcfun);
Re5 = Re3*sqrt(n/(n+1));
D3 = c6*T/Lx * exp((x-Lx)/(c0*T))*(1-exp(-1*Lx/(c6*T)));
q6e = r6scale*Re3*sum(D3([xk1 xk2]).*projk([xk1 xk2]))
% Turns out that q4l and q6l add an insignificant amount to the solution.
% Now, look at reflections from western boundaries.
n = 1;
c0=c;
c2=c/(2*n+1);
pcfun = expfun.*hermite(y, n+1);
pcfun = pcfun./sqrt(pcfun'*pcfun);
pcfun_p1 = hermite(y, n+1);
pcfun_p1 = 2^(-(n+1)/2)*pcfun_p1.*expfun;
pcfun_p1 = pcfun_p1./sqrt(pcfun_p1'*pcfun_p1);
pcfun_m1 = hermite(y, n-1);
pcfun_m1 = 2^(-(n-1)/2)*pcfun_m1.*expfun;
pcfun_m1 = pcfun_m1./sqrt(pcfun_m1'*pcfun_m1);
projr0 = (1/(2*n+1))*(n*pcfun_p1-sqrt(n*(n+1))*pcfun_m1)'*taux;
D4 = (c0*T/Lx)*exp(-1*((Lx+x)/(c2*T) + Lx/(c0*T)))*...
(1-exp(-1*Lx/(c0*T)));
rw2 = 1/sqrt(2);
q2r = rw2*sum(D4([xk1 xk2]).*projr0([xk1 xk2]));
n = 3;
c0=c;
c4=c/(2*n+1);
pcfun = expfun.*hermite(y, n+1);
pcfun = pcfun./sqrt(pcfun'*pcfun);
pcfun_p1 = hermite(y, n+1);
pcfun_p1 = 2^(-(n+1)/2)*pcfun_p1.*expfun;
pcfun_p1 = pcfun_p1./sqrt(pcfun_p1'*pcfun_p1);
pcfun_m1 = hermite(y, n-1);
pcfun_m1 = 2^(-(n-1)/2)*pcfun_m1.*expfun;
pcfun_m1 = pcfun_m1./sqrt(pcfun_m1'*pcfun_m1);
projr2 = (1/(2*n+1))*(n*pcfun_p1-sqrt(n*(n+1))*pcfun_m1)'*taux;
D4 = (c0*T/Lx)*exp(-1*((Lx+x)/(c4*T) + Lx/(c0*T)))*...
(1-exp(-1*Lx/(c0*T)));
rw4 = (sqrt(fact(2*(n-1)))/(2^(n-1)*fact(n-1)))/sqrt((n+1)*n);
q4r = rw4*sum(D4([xk1 xk2]).*projr2([xk1 xk2]));
n = 5;
c0=c;
c6=c/(2*n+1);
pcfun = expfun.*hermite(y, n+1);
pcfun = pcfun./sqrt(pcfun'*pcfun);
pcfun_p1 = hermite(y, n+1);
pcfun_p1 = 2^(-(n+1)/2)*pcfun_p1.*expfun;
pcfun_p1 = pcfun_p1./sqrt(pcfun_p1'*pcfun_p1);
pcfun_m1 = hermite(y, n-1);
pcfun_m1 = 2^(-(n-1)/2)*pcfun_m1.*expfun;
pcfun_m1 = pcfun_m1./sqrt(pcfun_m1'*pcfun_m1);
projr4 = (1/(2*n+1))*(n*pcfun_p1-sqrt(n*(n+1))*pcfun_m1)'*taux;
D4 = (c0*T/Lx)*exp(-1*((Lx+x)/(c6*T) + Lx/(c0*T)))*...
(1-exp(-1*Lx/(c0*T)));
rw6 = (sqrt(fact(2*(n-1)))/(2^(n-1)*fact(n-1)))/sqrt((n+1)*n);
q6r = rw6*sum(D4([xk1 xk2]).*projr2([xk1 xk2]));
% Locally forced Rossby waves
n=1;
pcfun = expfun.*hermite(y, n+1);
pcfun = pcfun./sqrt(pcfun'*pcfun);
D5 = (c2*T./x).*(1-exp(-1*x/(c2*T)));
q2l = sum(D5([xk2]).*projr0([xk2]))*r2scale;
n=3;
pcfun = expfun.*hermite(y, n+1);
pcfun = pcfun./sqrt(pcfun'*pcfun);
D5 = (c4*T./x).*(1-exp(-1*x/(c4*T)));
q4l = sum(D5([xk2]).*projr2([xk2]))*r4scale;
n=5;
pcfun = expfun.*hermite(y, n+1);
pcfun = pcfun./sqrt(pcfun'*pcfun);
D5 = (c6*T./x).*(1-exp(-1*x/(c6*T)));
q6l = sum(D5([xk2]).*projr4([xk2]))*r6scale;
[c T Lx]
[[sum([q0l q2l q4l q6l q2e q4e q6e q2r q4r q6r]) ...
sum([q0l q2l q4l q6l q2e q4e q6e]) ...
sum([q2r q4r q6r])]/1000 ...
sum([q2r q4r q6r])/sum([q0l q2l q4l q6l q2e q4e q6e])]
300 25920000 900000000
7.1024 7.5128 -0.4104 -0.0183
289 64800000 900000000 --> 40m/s ATM speed
7.1002 7.7188 -0.6186 -0.08something
6.7392 7.4590 -0.7199 -0.0965 --> 20m/s ATM speed