Global Index (short | long) | Local contents | Local Index (short | long)
Next, SST lims2 = [122.5 280 -25 25]
| This script calls | |
|---|---|
clear
cd /home/disk/tao/data/coads/coads1c
tim = (12*(1950-1800)+1):(12*(1990-1800)+12); % we'll seasonally average
lims2 = [30 290 -30 30]
[sst, lat, lon] = getnc2('sst.mean.nc', 'sst', lims2, 1, tim);
sst = squeeze(sst);
[ntim, nlat2, nlon2] = size(sst);
sst2 = repmat(NaN, [ntim, nlat2/2, nlon2/5]);
for tind = 1:ntim;
for i = 1:(nlat2/2);
ind1 = 2*(i-1)+[1:2];
for j = 1:(nlon2/5);
ind2 = 5*(j-1)+[1:5];
if length(find(~isnan(squeeze(sst(tind,ind1,ind2))))) > 2;
sst2(tind,i,j) = ...
squeeze(mean2(mean2(squeeze(sst(tind,ind1,ind2)))));
end
end
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/5], 1);
for j = 1:(nlon2/5);
ind2 = 5*(j-1)+[1:5];
lon2(j) = mean(lon(ind2));
end
sst = sst2; clear sst2 lat lon;
[ntim, nlat2, nlon2] = size(sst);
% 3-month running mean
sst = rave(sst, 3);
% Now, get to business
sst = annave(sst);
% Start by just computing the pops:
sst = cosweight(sst, lat2);
% Transform into first 15 eof's
nmode = 15;
% SST
tem = mean2(sst);
kp = find(~isnan(tem));
sst = sst(:, kp);
c = covar_nan3(sst, sst, 10, 'show');
c(isnan(c)) = 0;
[lds, lam] = eig(c);
l = diag(lam);
[lam, k] = sort(l'); lds = lds(:,k);
lam = fliplr(lam);
lds = fliplr(lds);
pcs = covar_nan3(sst', lds(:,1:nmode));
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;