Global Index (short | long) | Local contents | Local Index (short | long)
Load and configure NMC.REANAL data
| This script calls | |
|---|---|
clear
cd /home/disk/tao/data/sst_ldeo
filin = ['sstanomldeo18561991.nc'];
nc = netcdf(filin, 'nowrite');
ts = nc{'ssta'}(:,:,:);
datets = nc{'T'}(:);
lat = nc{'Y'}(:);
lon = nc{'X'}(:);
add_offset = nc{'ssta'}.add_offset(:);
scale_factor = nc{'ssta'}.scale_factor(:);
mv = nc{'ssta'}.missing_value(:);
nc = close(nc)
ts = ts * scale_factor;
ts = ts + add_offset;
mv = mv * scale_factor;
mv = mv + add_offset;
yr = floor(datets/12) + 1960;
mo = floor(datets + 1248) - (yr - 1856)*12 + 1;
[ntim, nlat, nlon] = size(ts);
ts = reshape(ts, ntim, nlat*nlon);
tsave = mean(ts);
kp = find(tsave ~= mv);
ts = ts(:, kp);
% Get climatology:
[ts, clim] = annave(ts);
% Get CT*:
getlp = 0;
if getlp;
tslp = myrunning_ave(ts, 25);
tslp = myrunning_ave(tslp, 37);
tshp = ts - tslp;
cd /home/disk/hayes2/dvimont/data
save ldeo_sst.mat tshp tslp ts;
else
cd /home/disk/hayes2/dvimont/data
load ldeo_sst.mat
end
ntim = size(ts, 1);
%kptime = 37:(ntim - 36);
kptime = 1:ntim;
tsnew = NaN * ones(ntim, nlat*nlon);
tsnew(:, kp) = tshp(kptime, :);
tsnew = reshape(tsnew, ntim, nlat, nlon);
ctlim = [180 270 -6 6];
[xk, yk] = keep_var(ctlim, (360 + lon), lat);
ctstar = squeeze(mean2(mean2(shiftdim(tsnew(:, yk, xk), 1))));
ctstar = (ctstar - mean(ctstar)) / std(ctstar);
ctpat = ctstar' * ts(kptime, :) / (length(ctstar));
tem = NaN * ones(1, nlat*nlon);
tem(kp) = ctpat;
tem = reshape(tem, nlat, nlon);
get_global
default_global
FRAME = [min(lon) (min(lon) + 360) -90 90];
figure(1)
sd(1)
gcont(tem, [-1:.1:1]);
dc
sd(2)
plot(ctstar)
set(gca, 'XTick', [12:120:1600], 'XTickLabel', [1860:10:2000])
axis([0 1561 -3 4])
grid
% Now, get EOFs:
clear tshp tslp tsnew
%[ts, clim] = annave(ts, clim, 1);
%[ts, clim] = annave(ts(kptime,:));
ntim = size(ts, 1);
% Check first years (1859-1949)
tkp = 1:(12*(1949-1856+1));
ntim = length(tkp);
[ts, dum] = annave(ts(tkp,:));
tsnew = NaN * ones(ntim, nlat*nlon);
tsnew(:, kp) = ts;
tsnew = reshape(tsnew, ntim, nlat, nlon);
tsnew = cosweight(tsnew, lat);
tsnew = reshape(tsnew, ntim, nlat*nlon);
tsnew = tsnew(:, kp);
c = tsnew' * tsnew;
[lam, lds, per] = eof(c);
lds10 = (ones(size(ts, 2), 1) * (1./(std(lds(:,1:10))))) .* lds(:, 1:10);
pc10 = ts * lds10 / size(ts, 2);
pc10 = ones(ntim,1) * (1./(std(pc10))) .* pc10;
lds10 = pc10'*ts / ntim;
num = 1;
tem = NaN * ones(1, nlat*nlon);
tem(:, kp) = lds10(num, :);
tem = reshape(tem, nlat, nlon);
figure(1)
sd(1)
gcont(-1 * tem, [-2:.1:2]);
dc
sd(2)
plot(-1 * pc10(:, num));
% Now, regress CT* against G and yyy...
ctstar = ctstar(1:ntim);
a1 = corr(ctstar, pc10(1:ntim, 1));
gr = pc10(1:ntim,1) - a1 * ctstar;
gr = (gr - mean(gr)) / std(gr);
grpat = gr' * ts / ntim;
tem = NaN * ones(1, nlat*nlon);
tem(:, kp) = grpat;
tem = reshape(tem, nlat, nlon);
figure(1)
sd(1)
gcont(-1 * tem, [-2:.1:2]);
dc
sd(2)
plot(-1 * gr);
% Check out LP and HP filtered GR:
grlp = myrunning_ave(myrunning_ave(gr, 25), 37);
grhp = gr - grlp;
grlp = (grlp - mean(grlp)) / std(grlp);
grhp = (grhp - mean(grhp)) / std(grhp);
grlppat = grlp' * ts / ntim;
grhppat = grhp' * ts / ntim;
tem = NaN * ones(1, nlat*nlon);
tem(:, kp) = grhppat;
tem = reshape(tem, nlat, nlon);
figure(1)
sd(1)
gcont(-1 * tem, [-2:.1:2]);
dc
sd(2)
plot(-1 * grhp);
% Check lag/lead relationship:
ctstar = -1 * ctstar;
corr_grct = [];
for i = -24:12;
corr_i = corr(ctstar, gr, i);
corr_grct = [corr_grct corr_i];
end
subplot(2,1,1)
bar(-24:12, corr_grct)
grid
title(['Lagged Correlations between GR and CT*'])
xlabel(['Months: Negative => CT* leads GR; Positive => CT* lags GR'])
ylabel(['Correlation Coefficient'])
axis([-25 13 -0.2 0.5])
corr_grct = [];
corr_hpgrct = [];
for i = -24:12;
corr_i = corr(ctstar, grhp, i);
corr_grct = [corr_grct corr_i];
end
subplot(2,1,2)
bar(-24:12, corr_grct)
grid
title(['Lagged Correlations between High Pass GR and CT*'])
xlabel(['Months: Negative => CT* leads GR; Positive => CT* lags GR'])
ylabel(['Correlation Coefficient'])
axis([-25 13 -0.5 0.8])