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 = ['air.mon.mean.srfc.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(:);
mv = nc{'air'}.missing_value(:);
nc = close(nc)
ts(find(ts == mv)) = NaN*ones(size(find(ts == mv)));
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 = interp2(loncyc, lmlat, lmcyc, lon, lat');
% Interpolate to a more coarse grid.
tem = NaN*ones(ntim,((nlat-1)/2+1),(nlon/4));
temlnd = NaN*ones(((nlat-1)/2+1),(nlon/4));
temlon = NaN*ones((nlon/4),1);
temlat = NaN*ones((nlat-1)/2+1, 1);
for lati = 1:((nlat-1)/4);
latind = 2*(lati-1)+[1:2];
for loni = 1:nlon/4;
lonind = 4*(loni-1)+[1:4];
tem(:,lati,loni) = mean2(mean2(shiftdim(ts(:,latind,lonind),1)));
temlnd(lati,loni) = mean2(mean2(lm2(latind,lonind)));
end
temlat(lati) = mean2(lat(latind));
end
for loni = 1:nlon/4;
lati = ((nlat-1)/4+1);
latind = 37;
lonind = 4*(loni-1)+[1:4];
tem(:,lati,loni) = mean2(shiftdim(ts(:,latind,lonind),2));
temlnd(lati,loni) = mean2(mean2(lm2(latind,lonind)));
temlon(loni) = squeeze(mean2(lon(lonind)));
end
temlat((nlat-1)/4+1) = mean2(lat(latind));
for lati = ((nlat-1)/4 + 2):((nlat-1)/2 + 1);
latind = 2*(lati-1)+[0:1];
for loni = 1:nlon/4;
lonind = 4*(loni-1)+[1:4];
tem(:,lati,loni) = mean2(mean2(shiftdim(ts(:,latind,lonind),1)));
temlnd(lati,loni) = mean2(mean2(lm2(latind,lonind)));
end
temlat(lati) = mean2(lat(latind));
end
ts = tem;
lat = temlat;
lon = temlon;
lm2 = temlnd;
[ntim, nlat, nlon] = size(ts);
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) <= (ntim-1);
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))));
cd /home/disk/hayes2/dvimont/data
%save sst_1958-1999_coarse.mat ts kpt lat lon
%
%
% Calculate ctstar:
%
%
cd /home/disk/hayes2/dvimont/data
load sst_1958-1999_coarse.mat
[ntim, nkp] = size(ts);
nlat = length(lat); nlon = length(lon);
[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(mean2(mean2(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);
% 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);
%tsnct = ts;
% Method 2: calculate eofs, then remove ctstar
tsnew = zeros(ntim, nlat*nlon);
tsnew(:, kpt) = tsnct;
tsnew = cosweight(tsnew, lat);
tsnct = tsnew(:, kpt);
clear tsnew
[lam, per, lds, pcs] = eof_dan(tsnct);
if 0
cd /home/disk/tao/data/nmc.reanalysis/monthly
filin = ['air.mon.mean.srfc.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 = lat(yk);
[ntim, nlat, nlon] = size(ts);
ts = reshape(ts, ntim, nlat*nlon);
[ts, clim] = annave(ts);
end
% Try some plots:
default_global
FRAME = [0 360 -60 60];
num = 1;
cint = 0.1;
tem = NaN * ones(nlat*nlon,1);
tem(kpt,1) = lds(:, num);
%tem = lds(:, num);
tem = reshape(tem, nlat, nlon);
sd(1)
gcont(1*tem, [-5:cint:5])
dc
title(['NMC: EOF' num2str(num) ' of raw SST'])
xlabel(['Contour Interval: ' num2str(cint) ' K (std)^-^1'])
sd(2)
plot(1*pcs(:,num))
set(gca, 'XTick',[25:60:500],'XTickLabel',[60:5:100],'YTick',[-3:3])
axis([0 ntim+1 -2.5 2])
xlabel('Year: ticks correspond to January')
ylabel('STD')
title(['PC' num2str(num)])
grid
cd /home/disk/tao/dvimont/matlab/NMC/Plots2
gr1 = pcs(:,1) - corr(pcs(:,1),ctstar)*ctstar;
gr1 = (gr1 - mean(gr1)) / std(gr1);
cd /home/disk/hayes2/dvimont/data
%save gr_ct_pacific_only.mat gr ctstar
cd /home/disk/tao/data/nmc.reanalysis/monthly
filin = ['sst.mon.mean.nc'];
nc = netcdf(filin, 'nowrite');
time = nc{'time'}(:);
nc = close(nc);
dpy = (365 * 400 + 97) / 400;
yr = floor( time / (dpy * 24) + 1);