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/nmc.reanalysis/monthly
%cd /home/disk/tao/data/reynolds/eof
filin = ['sst.mon.mean.nc'];
%filin = ['ssteof5092.nc'];
nc = netcdf(filin, 'nowrite');
ts = nc{'air'}(:,:,:);
datets = nc{'time'}(:);
lat = nc{'lat'}(:);
lon = nc{'lon'}(:);
add_offset = nc{'air'}.add_offset(:);
scale_factor = nc{'air'}.scale_factor(:);
nc = close(nc)
ts = ts * scale_factor;
ts = ts + add_offset;
[ntim, nlat, nlon] = size(ts);
% Load landmask
cd /home/disk/tao/data/nmc.reanalysis
filin = ['landmask.nc'];
nc = netcdf(filin, 'nowrite');
lm = nc{'lsmask'}(:);
lmlat = nc{'lat'}(:);
lmlon = nc{'lon'}(:);
nc = close(nc);
[m,n] = size(lm);
loncyc = [[lmlon((n-9):n) - 360]; lmlon; [lmlon(1:10) + 360]];
lmcyc = [lm(:,(n-9):n) lm lm(:,1:10)];
lm2 = sign(interp2(loncyc, lmlat, lmcyc, lon, lat'));
lndpt = find(reshape(lm2, 1, nlat*nlon) == -1);
ocpt = find(reshape(lm2, 1, nlat*nlon) == 1);
ts = reshape(ts, ntim, nlat*nlon);
ts = ts(:, ocpt);
% Eliminate ice points:
ticept = [];
for i = 1:length(ocpt);
icept = find(ts(:,i) <= -1.8);
seapt = find(~ismember([1:ntim], icept));
if length(seapt) <= 120;
ticept = [ticept i];
elseif ~isempty(icept);
ts(icept, i) = mean(ts(seapt, i)) * ones(size(icept));
end
end;
maskpt = union(lndpt, ocpt(ticept));
kpt = ocpt(~ismember(ocpt, ocpt(ticept)));
ts = ts(:, kpt);
% Get highpass filtered and lowpass filtered datasets if called for:
gethp = 0; % Set to 1 for hp/lp stuff. This is to save time.
if gethp;
[ts, clim] = annave(ts);
tslp = myrunning_ave(ts, 25);
tslp = myrunning_ave(tslp, 37);
tshp = ts - tslp;
cd /home/disk/hayes2/dvimont/data
save nmc_hp_sst.mat tshp
save nmc_lp_sst.mat tslp
save nmc_noice_sst.mat ts
clear tslp ts tshp
save lat_lon_kpt.mat lat lon kpt
end
% Calculate CT*:
cd /home/disk/hayes2/dvimont/data
load nmc_hp_sst.mat
load nmc_logistics.mat
tkp = 1:ntim;
tshp = tshp(tkp,:);
ntim = size(tshp, 1);
tsnew = NaN * ones(ntim, nlat*nlon);
tsnew(:, kpt) = tshp;
tsnew = reshape(tsnew, ntim, nlat, nlon);
ctlim = [180 270 -6 6];
[xk, yk] = keep_var(ctlim, lon, lat);
ctstar = squeeze(mean2(mean2(shiftdim(tsnew(:, yk, xk), 1))));
clear tshp
load nmc_noice_sst.mat
ts = ts(tkp, :);
tsnew = reshape(tsnew, ntim, nlat*nlon);
tsnew(:, kpt) = ts;
ctpat = ((ctstar - mean(ctstar)) / std(ctstar))' * tsnew / ntim;
ctpat = reshape(ctpat, nlat, nlon);
cd /home/disk/hayes2/dvimont/data
ctstar = (ctstar - mean(ctstar)) / std(ctstar);
%save ct_gr.mat ctstar ctpat
load ct_gr.mat
get_global
define_global
figure(1)
sd(1)
gcont(ctpat, [-5:.2:5])
dc
title(['NMC: SST regressed on CT*'])
xlabel(['Contour Interval: 0.2 K (std)^-^1'])
sd(2)
plot((ctstar - mean(ctstar)) / std(ctstar))
set(gca, 'XTick',[49:60:409],'XTickLabel',[65:5:100],'YTick',[-3:3])
axis([0 409 -3 3.5])
xlabel('Year: ticks correspond to January')
ylabel('STD')
title('CT* Time Series')
grid
cd /home/disk/hayes2/dvimont/data/Plots
% Now, calculate G time series:
cd /home/disk/hayes2/dvimont/data
load nmc_noice_sst.mat
ntim = size(ts, 1);
%[ts, clim] = annave(ts, clim, 1);
%[ts, clim] = annave(ts(37:(ntim-36),:));
%ntim = size(ts, 1);
tsnew = NaN * ones(ntim, nlat*nlon);
tsnew(:, kpt) = ts;
tsnew = reshape(tsnew, ntim, nlat, nlon);
tsnew = cosweight(tsnew, lat);
tsnew = reshape(tsnew, ntim, nlat*nlon);
c = tsnew(:, kpt) * tsnew(:, kpt)' / length(kpt);
[lam, pc, per] = eof(c);
pc10 = ones(ntim,1) * (1./(std(pc(:,1:10)))) .* pc(:,1:10);
lds = pc10'*ts / ntim;
num = 1;
cint = 0.2;
tem = NaN * ones(1, nlat*nlon);
tem(1, kpt) = lds(num, :);
tem = reshape(tem, nlat, nlon);
sd(1)
gcont(tem, [-5:cint:5])
dc
title(['NMC: EOF' num2str(num) ' of raw SST'])
xlabel(['Contour Interval: ' num2str(cint) ' K (std)^-^1'])
sd(2)
plot(pc10(:,num))
set(gca, 'XTick',[49:60:409],'XTickLabel',[65:5:100],'YTick',[-3:3])
axis([0 409 -3 3.5])
xlabel('Year: ticks correspond to January')
ylabel('STD')
title(['PC' num2str(num)])
grid
cd /home/disk/hayes2/dvimont/data/Plots
% Now, calculate residual time series:
a1 = ctstar' * pc10(:,1) / ntim;
gr = pc10(:,1) - a1*ctstar;
gr = (gr - mean(gr)) / std(gr);
grpat = gr' * ts / ntim;
cint = 0.1;
tem = NaN * ones(1, nlat*nlon);
tem(1, kpt) = grpat;
tem = reshape(tem, nlat, nlon);
sd(1)
gcont(tem, [-5:cint:5])
dc
title(['NMC: GR pattern'])
xlabel(['Contour Interval: ' num2str(cint) ' K (std)^-^1'])
sd(2)
plot(gr)
set(gca, 'XTick',[49:60:409],'XTickLabel',[65:5:100],'YTick',[-3:3])
axis([0 409 -3 3.5])
xlabel('Year: ticks correspond to January')
ylabel('STD')
title(['GR time series'])
grid
cd /home/disk/hayes2/dvimont/data/Plots
% Check out other possibilities:
tkp = 37:ntim-36;
tkp = 1:ntim;
grlp = myrunning_ave(gr, 25);
grlp = myrunning_ave(grlp, 37);
grlp = (grlp(tkp) - mean(grlp(tkp))) / std(grlp(tkp));
grlppat = grlp' * ts(tkp, :) / length(tkp);
cint = 0.1;
tem = NaN * ones(1, nlat*nlon);
tem(1, kpt) = grlppat;
tem = reshape(tem, nlat, nlon);
tem = myrunning_ave(myrunning_ave(tem', 3)', 3);
sd(1)
gcont(tem, [-5:cint:5])
dc
title(['NMC: Low Pass GR pattern'])
xlabel(['Contour Interval: ' num2str(cint) ' K (std)^-^1'])
sd(2)
plot(grlp)
set(gca, 'XTick',[49:60:409],'XTickLabel',[65:5:100],'YTick',[-3:3])
axis([0 409 -2 2])
xlabel('Year: ticks correspond to January')
ylabel('STD')
title(['Low Pass GR time series'])
grid
cd /home/disk/hayes2/dvimont/data/Plots
grhp = gr - myrunning_ave(myrunning_ave(gr, 25), 37);
grhp = (grhp - mean(grhp)) / std(grhp);
grhppat = grhp' * ts / ntim;
cint = 0.1;
tem = NaN * ones(1, nlat*nlon);
tem(1, kpt) = grhppat;
tem = reshape(tem, nlat, nlon);
sd(1)
gcont(tem, [-5:cint:5])
dc
title(['NMC: High Pass GR pattern'])
xlabel(['Contour Interval: ' num2str(cint) ' K (std)^-^1'])
sd(2)
plot(grhp)
set(gca, 'XTick',[49:60:409],'XTickLabel',[65:5:100],'YTick',[-3:3])
axis([0 409 -3 3.5])
xlabel('Year: ticks correspond to January')
ylabel('STD')
title(['High Pass GR time series'])
grid
cd /home/disk/hayes2/dvimont/data/Plots
% Check lag/lead relationship:
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])
cd /home/disk/hayes2/dvimont/data/Plots
%
% Now, add GRLPPAT to sst to force CCM:
cd /home/disk/hayes/dvimont/ccm3.6/data/GR
filin = 'T42M5079.nc'
nc = netcdf(filin, 'nowrite');
clat = nc{'lat'}(:);
clon = nc{'lon'}(:);
sst = nc{'SST'}(:);
nc = close(nc);
lonb = [lon((nlon-4):nlon)-360; lon; lon(1:5)+360];
temb = [tem(:,(nlon-4):nlon) tem tem(:,1:5)];
cgrlp = interp2(lonb, lat, temb, clon, clat');
figure(2)
XAX = clon;
YAX = clat;
gcont(cgrlp, [-5:cint:5]);
dc
cgrlp(find(isnan(cgrlp))) = zeros(size(find(isnan(cgrlp))));
icept = min(min(min(sst)));
sstw = icept * ones(size(sst));
for i = 1:12;
kpsst = find(squeeze(sst(i,:)) > icept);
sstw(i,kpsst) = cgrlp(kpsst) + squeeze(sst(i,kpsst));
end;
filin = 'gr_warm_sst.nc';
nc = netcdf(filin, 'write');
nc{'SST'}(:) = sstw;
nc = close(nc);
sstc = icept * ones(size(sst));
for i = 1:12;
kpsst = find(squeeze(sst(i,:)) > icept);
sstc(i,kpsst) = -1*cgrlp(kpsst) + squeeze(sst(i,kpsst));
end;
filin = 'gr_cold_sst.nc';
nc = netcdf(filin, 'write');
nc{'SST'}(:) = sstc;
nc = close(nc);