Global Index (short | long) | Local contents | Local Index (short | long)
Get data Start with SLP
| This script calls | |
|---|---|
clear
cd /home/disk/tao/data/nmc.reanalysis/monthly
tim = (12*(1950-1948)+1):(12*(1990-1948)+12); % we'll seasonally average
lims = [117.5 272.5 20 89];
[slp, lat, lon] = getnc2('slp.mon.mean.nc', 'slp', lims, 1, tim);
slp = squeeze(slp);
[ntim, nlat1, nlon1] = size(slp);
slp2 = repmat(NaN, [ntim, nlat1/2, nlon1/3]);
for i = 1:(nlat1/2);
ind1 = 2*(i-1)+[1:2];
for j = 1:(nlon1/3);
ind2 = 3*(j-1)+[1:3];
slp2(:,i,j) = squeeze(mean(mean(shiftdim(slp(:,ind1,ind2), 1))));
end
end
lat1 = repmat(NaN, [nlat1/2], 1);
for i = 1:(nlat1/2);
ind1 = 2*(i-1)+[1:2];
lat1(i) = mean(lat(ind1));
end
lon1 = repmat(NaN, [nlon1/3], 1);
for j = 1:(nlon1/3);
ind2 = 3*(j-1)+[1:3];
lon1(j) = mean(lon(ind2));
end
slp = slp2; clear slp2 lat lon;
% Next, SST
%lims2 = [122.5 280 -25 25]
lims2 = [30 292 -30 30]
[sst, lat, lon] = getnc2('skt.mon.mean.nc', 'skt', lims2, 1, tim);
sst = squeeze(sst);
[ntim, nlat2, nlon2] = size(sst);
sst2 = repmat(NaN, [ntim, nlat2/2, nlon2/4]);
% eliminate land points
lm = getnc2('land.sfc.gauss.nc', 'land', lims2, 1, 1);
land = find(lm);
ocean = find(1-lm);
sst(:,land) = NaN;
for i = 1:(nlat2/2);
ind1 = 2*(i-1)+[1:2];
for j = 1:(nlon2/4);
ind2 = 4*(j-1)+[1:4];
sst2(:,i,j) = squeeze(mean2(mean2(shiftdim(sst(:,ind1,ind2), 1))));
end
end
lat2 = repmat(NaN, [nlat2/2], 1);
for i = 1:(nlat2/2);
ind1 = 2*(i-1)+[1:2];
lat2(i) = mean(lat(ind1));
end
lon2 = repmat(NaN, [nlon2/4], 1);
for j = 1:(nlon2/4);
ind2 = 4*(j-1)+[1:4];
lon2(j) = mean(lon(ind2));
end
sst = sst2; clear sst2 lat lon;
[ntim, nlat1, nlon1] = size(slp);
[ntim, nlat2, nlon2] = size(sst);
% Now, get to business
sst = annave(sst);
slp = annave(slp);
% 3-month running mean
sst = rave(sst, 3);
slp = rave(slp, 3);
% Start by just computing the pops:
slp = cosweight(slp, lat1);
sst = cosweight(sst, lat2);
% Transform into first 15 eof's
nmode = 15;
% SLP
slp = reshape(slp, ntim, nlat1*nlon1);
c = slp'*slp/(ntim-1);
[lds, lam] = eig(c);
l = diag(lam);
[lam, k] = sort(l'); lds = lds(:,k);
lam = fliplr(lam);
lds = fliplr(lds);
pcs = slp*lds;
slp2 = pcs(:,1:nmode);
slplam = lam;
slppcs = pcs;
slplds = lds;
clear lam lds pcs l;
% SST
tem = mean(sst);
kp = find(~isnan(tem));
sst = sst(:, kp);
c = sst'*sst/(ntim-1);
[lds, lam] = eig(c);
l = diag(lam);
[lam, k] = sort(l'); lds = lds(:,k);
lam = fliplr(lam);
lds = fliplr(lds);
pcs = sst*lds;
sst2 = pcs(:,1:nmode);
sstlam = lam;
sstpcs = pcs;
sstlds = lds;
clear lam lds pcs l;
% Now, compute B: Start with just SST
lag = 7;
tem = sst2; % detrend(sst2);
c0 = tem'*tem/(ntim-1);
ctau = sst2(1:(ntim-lag),:)'*sst2([1:(ntim-lag)]+lag,:)/(ntim-lag-1);
B = logm(ctau*inv(c0))/lag;
[ldsb, lamb] = eig(B);
l = diag(lamb);
[tem, k] = sort(real(l)');
lamb = diag(lamb); lamb = lamb(k); ldsb = ldsb(:,k);
lamb = fliplr(lamb'); ldsb = fliplr(ldsb);
% Now, reconstruct modes from sstlds
figure_tall(1); clf;
colormap default
global_axes(6, 2, 1/3, 0.5, 1);
global_latlon(lat2, lon2, lims2);
num = 1; wgt = 0.15;
num = 4; wgt = 0.15;
tem = ldsb(:,num)'*sstlds(:,1:15)';
tem2 = repmat(NaN, nlat2, nlon2);
tem2(kp) = tem;
% Case 1: lamb is real
subplot2(1,1); cla;
map_axis('giso');
map_surface_interp(real(exp(sqrt(-1)*0)*tem2), -0.5);
fill_landmap('over');
tightmap
caxis([-1 1]*wgt);
t1 = title(['Damp: ' num2str(round(-100/real(lamb(num)))/100) ...
', Period: ' num2str(round(-100/imag(lamb(num)))/100)]);
set(t1, 'fontsize', 10);
% Case 2: lamb is complex
num = 2; wgt = 0.15;
num = 3; wgt = 0.15;
num = 4; wgt = 0.15;
num = 7; wgt = 0.15;
num = 11; wgt = 0.15;
num = 13; wgt = 0.15;
eta = 0.5*atan2(-2*real(ldsb(:,num))'*imag(ldsb(:,num)), ...
real(ldsb(:,num))'*real(ldsb(:,num)) - ...
imag(ldsb(:,num))'*imag(ldsb(:,num)));
ldstem = exp(sqrt(-1)*eta)*ldsb(:,num);
tem = conj(ldstem)'*sstlds(:,1:15)';
tem2 = repmat(NaN, nlat2, nlon2);
tem2(kp) = tem;
subplot2(1,1); cla;
map_axis('giso');
map_surface_interp(real(exp(sqrt(-1)*0)*tem2), -0.5);
fill_landmap('over');
tightmap
caxis([-1 1]*wgt);
t1 = title(['Damp: ' num2str(round(-100/real(lamb(num)))/100) ...
', Period: ' num2str(round(-100/imag(lamb(num)))/100)]);
set(t1, 'fontsize', 10);
subplot2(1,2); cla;
map_axis('giso');
map_surface_interp(real(exp(sqrt(-1)*pi/4)*tem2), -0.5);
fill_landmap('over');
tightmap
caxis([-1 1]*wgt);
subplot2(1,3); cla;
map_axis('giso');
map_surface_interp(real(exp(sqrt(-1)*2*pi/4)*tem2), -0.5);
fill_landmap('over');
tightmap
caxis([-1 1]*wgt);
subplot2(1,4); cla;
map_axis('giso');
map_surface_interp(real(exp(sqrt(-1)*3*pi/4)*tem2), -0.5);
fill_landmap('over');
tightmap
caxis([-1 1]*wgt);
% Now, redo this with G
lag = 7;
c0 = sst2'*sst2/(ntim-1);
ctau = sst2(1:(ntim-lag),:)'*sst2([1:(ntim-lag)]+lag,:)/(ntim-lag-1);
G = ctau*inv(c0);
[ldsb, lamb] = eig(G'*G);
l = diag(lamb);
[tem, k] = sort(real(l)');
lamb = diag(lamb); lamb = lamb(k); ldsb = ldsb(:,k);
lamb = fliplr(lamb'); ldsb = fliplr(ldsb);
ldsb = G*ldsb(:,1);
[u, s, v] = svd(G);
ldsb = u;
ldsb = v;