Documentation of svd_intraseasonal

Global Index (short | long) | Local contents | Local Index (short | long)

Help text

  Get landmask

Cross-Reference Information

Listing of script svd_intraseasonal

clean
cd /home/disk/tao/data/nmc.reanalysis
load landmask.mat
lm(lm == -1) = NaN;

cd /home/disk/tao/data/nmc.reanalysis/monthly
uwnlim = [120 285 -20 20];
slplim = [110 270 20 90];

filin1 = 'uwnd.mon.mean.nc'; var1 = 'uwnd';
filin2 = 'slp.mon.mean_48_Current.nc'; var2 = 'slp';
tim = 11:622;

[uwn, lat1, lon1] = getnc2(filin1, var1, uwnlim, 925, tim);
[slp, lat2, lon2] = getnc2(filin2, var2, slplim, 500, tim);

uwn = squeeze(uwn); slp = squeeze(slp);

[uwn, uwnc] = annave(uwn);
[slp, slpc] = annave(slp);

[ntim, nlat1, nlon1] = size(uwn);
[ntim, nlat2, nlon2] = size(slp);

uwn2 = reshape(uwn, ntim, nlat1*nlon1);
slp2 = reshape(slp, ntim, nlat2*nlon2);

%  Get time indices - 0NDJFMA for SLP, NDJFMAJ for UWN
ind1 = sort([1:12:ntim 2:12:ntim 3:12:ntim 4:12:ntim ...
	    5:12:ntim 6:12:ntim]);
ind2 = ind1+6;
nmo = length(11:16);

uwn2 = uwn2(ind2,:);
slp2 = slp2(ind2,:);

ntim = length(ind1);

%  Remove winter time mean from data
nyr = floor(ntim/nmo);

for i = 1:nyr;
  ind = nmo*(i-1)+[1:nmo];
  uwn2(ind,:) = uwn2(ind,:)-ones(nmo,1)*mean(uwn2(ind,:));
  slp2(ind,:) = slp2(ind,:)-ones(nmo,1)*mean(slp2(ind,:));
end

%  Normalize, as in Deser and Timlin
if 0;
  uwnstd = std(uwn2);
  slpstd = std(slp2);

  uwn2 = uwn2 ./ (ones(ntim,1)*uwnstd);
  slp2 = slp2 ./ (ones(ntim,1)*slpstd);
end

%  Weight by cosine latitude
uwn2 = reshape(uwn2, ntim, nlat1, nlon1);
uwn2 = cosweight(uwn2, lat1);
uwn2 = reshape(uwn2, ntim, nlat1*nlon1);
slp2 = reshape(slp2, ntim, nlat2, nlon2);
slp2 = cosweight(slp2, lat2);
slp2 = reshape(slp2, ntim, nlat2*nlon2);

%  Get rid of land points for UWN
if 0;
  [xk, yk] = keep_var(uwnlim, lon, lat);
  kp = find(~isnan(lm(yk, xk)));
  uwn2 = uwn2(:,kp);
end

%  Form covariance matrix
c = slp2'*uwn2./(ntim-1);

%  Do SVD
[slpu, s, uwnv] = svd(c);
slpx = slp2*slpu;
uwny = uwn2*uwnv;

rmsc = (sum(sum(c.^2))./(sum(var(slp2))*sum(var(uwn2))))^0.5;
scv = diag(s).^2./sum(diag(s).^2);