Global Index (short | long) | Local contents | Local Index (short | long)
Load landmask
| This script calls | |
|---|---|
clear
get_global
cd /home/disk/tao/data/nmc.reanalysis/monthly
filin = ['sst.mon.mean.nc'];
nc = netcdf(filin, 'nowrite');
lat = nc{'lat'}(:);
lon = nc{'lon'}(:);
[xk, yk] = keep_var([0 360 -60 60], lon, lat);
ts = nc{'air'}(1:480,yk,xk);
datets = nc{'time'}(:);
add_offset = nc{'air'}.add_offset(:);
scale_factor = nc{'air'}.scale_factor(:);
nc = close(nc)
ts = ts * scale_factor;
ts = ts + add_offset;
lon = lon(xk);
lat = lat(yk);
[ntim, nlat, nlon] = size(ts);
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);
kpt = ocpt;
ts = reshape(ts, ntim, nlat*nlon);
ts = ts(:, ocpt);
% Eliminate ice points:
ticept = [];
for i = 1:length(ocpt);
icept = find(ts(:,i) <= 0.);
seapt = find(~ismember([1:ntim], icept));
if length(seapt) <= 479;
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(:, find(~ismember(ocpt, ocpt(ticept))));
% Calculate ctstar:
[ctx, cty] = keep_var([180 270 -6 6], lon, lat);
tsnew = zeros(ntim, nlat*nlon);
tsnew(:,kpt) = ts;
tsnew = reshape(tsnew, ntim, nlat, nlon);
ct = squeeze(mean(mean(shiftdim(tsnew(:,cty,ctx),1))));
[ct, tem] = annave(ct);
ctstar = ct - myrunning_ave(myrunning_ave(ct, 25), 37);
ctstar = (ctstar - mean(ctstar)) / std(ctstar);
% Now, calculate GR time series.
[ts, clim] = annave(ts);
if 0
tslp = rave(ts, 25);
tslp = rave(tslp, 37);
[tslp, clim] = annave(tslp);
cd /home/disk/hayes2/dvimont/data
save nmc_lp_sst.mat tslp lat lon ntim nlat nlon kpt
end
% Method 1: remove ctstar time series from all points:
clear tsnew
a1 = ctstar' * ts / ntim;
tsnct = ts - ctstar*a1;
tsnct = tsnct - ones(ntim, 1) * mean(tsnct);
if 0
tslpnct = rave(ts, 25);
tslpnct = rave(tslpnct, 17);
[tslpnct, clim] = annave(tslpnct);
cd /home/disk/hayes2/dvimont/data
save nmc_lpnoct_sst.mat tslpnct lat lon ntim nlat nlon kpt
end
%tsnct = ts;
tsnct = tslpnct;
clear tslpnct
% Method 2: calculate eofs, then remove ctstar
tsnew = zeros(ntim, nlat*nlon);
tsnew(:, kpt) = tsnct;
tsnew = cosweight(tsnew, lat);
tsnct = tsnew(:, kpt);
clear tsnew
tim = [25:456]; ntim = length(tim);
tsnct = tsnct(tim,:);
c = tsnct * tsnct' / (length(kpt));
[lam, pcs, per] = eof(c);
if 0
cd /home/disk/tao/data/nmc.reanalysis/monthly
filin = ['sst.mon.mean.nc'];
nc = netcdf(filin, 'nowrite');
lat = nc{'lat'}(:);
lon = nc{'lon'}(:);
[xk, yk] = keep_var([0 360 -90 90], lon, lat);
ts = nc{'air'}(:,yk,xk);
datets = nc{'time'}(:);
add_offset = nc{'air'}.add_offset(:);
scale_factor = nc{'air'}.scale_factor(:);
nc = close(nc)
ts = ts * scale_factor;
ts = ts + add_offset;
lon = lon(xk);
lat = lat(yk);
[ntim, nlat, nlon] = size(ts);
ts = reshape(ts, ntim, nlat*nlon);
[ts, clim] = annave(ts);
elseif 0
cd /home/disk/hayes2/dvimont/data
load nmc_lpnoct_sst.mat
ts = tslpnct; clear tslpnct;
end
pc10 = pcs(:,1:10);
%pc10 = myrunning_ave(pcs(:,1:10), 9);
%pc10 = myrunning_ave(pc10, 5);
pc10 = ones(ntim,1) * (1./(std(pc10(:,1:10)))) .* pc10(:,1:10);
%ts = ts(:, kpt);
tsnct = ts;
a1 = ctstar' * ts / ntim;
tsnct = ts - ctstar*a1;
tsnct = tsnct - ones(ntim, 1) * mean(tsnct);
lds = pc10'*tsnct / ntim;
default_global
FRAME = [0 360 -60 60];
orient landscape
for num = 1:4;
cint = 0.2;
tem = NaN * ones(1, nlat*nlon);
%tem(1, kpt) = lds(num, :);
%tem(1, kpt) = tem2(num, :);
tem(1, kpt) = grpat(num, :);
%tem = lds(num, :);
tem = reshape(tem, nlat, nlon);
%sd(2)
% plot(1*pc10(:,num))
% plot(gr2)
% set(gca, 'XTick',[25:60:500],'XTickLabel',[60:5:100],'YTick',[-3:3])
% set(gca, 'XTick',[49:60:500],'XTickLabel',[65:5:100],'YTick',[-3:3])
% axis([0 (length(yrs)+1) -2.5 2])
% xlabel('Year: ticks correspond to January')
% ylabel('STD')
% title(['PC' num2str(num)])
% grid
%subplot(2,2,num-0)
% if ismap(gca); clma; end
sd(1)
gcont(1*tem, [-5:cint:5])
dc
title(['NMC: LP EOF1, MON = ' num2str(num)])
xlabel(['Contours: ' num2str(cint) ' K (std)^-^1'])
end
cd /home/disk/tao/dvimont/matlab/CCM/GR/GR_Plots
gr1 = pc10(:,1) - corr(pc10(:,1),ctstar)*ctstar;
gr1 = (gr1 - mean(gr1)) / std(gr1);
cd /home/disk/hayes2/dvimont/data
tem1 = gr1' * ts / ntim;
yrs = (37:444);
gr2 = rave(rave(gr1(yrs), 25), 37);
gr2 = (gr2 - mean(gr2)) / std(gr2);
tem2 = gr2' * tsnct(yrs,:) / length(yrs);
gr1 = pc10(:,1);
for i = 1:12;
yrs1 = (36 + [i:12:408]);
% tem = rave(rave(gr1, 25), 37);
tem = gr1(yrs1);
tem = (tem - mean(tem)) / std(tem);
grpat(i,:) = tem' * tsnct(yrs1,:) / length(yrs1);
end
cd /home/disk/hayes2/dvimont/ccm/ccm3.6/data
save seas_grpat gr2 grpat
% Check out highpass pc's:
hppc1 = pc10_withct - rave(rave(pc10_withct,25),37);
hppc2 = pc10 - rave(rave(pc10,25),37);
hppc1 = (hppc1 - ones(ntim, 1)*mean(hppc1)) ./ (ones(ntim, 1)*std(hppc1));
hppc2 = (hppc2 - ones(ntim, 1)*mean(hppc2)) ./ (ones(ntim, 1)*std(hppc2));
%
% I think the best method to do this is the following:
%
% 1. Load region [0 360 -60 60].
% 2. Eliminate all points that ever had ice.
% 3. Calculate HP CT (CT*).
% 4. Regress TS onto CT*, resulting in n_spatial_points regression
% coefficients.
% 5. Subtract CT* times the regression coefficient from each spatial
% point's time series.
% 6. Low Pass filter the data to remove all high-frequency phenomenon
% associated with high frequency component of ENSO.
% 7. Compute EOF's on this new, low-pass filtered, CT-free data.
%