Global Index (short | long) | Local contents | Local Index (short | long)
clear cd /home/disk/hayes2/dvimont/csm/300yr/srf load yrs151-200.mat sst1 = sst; oro1 = oro; load yrs201-250.mat sst2 = sst; oro2 = oro; load yrs251-300.mat sst = [sst1; sst2; sst]; clear sst1 sst2 oro = [oro1; oro2; oro]; clear oro1 oro2 save sst151-300 sst save oro151-300 oro
| This script calls | |
|---|---|
% The folowing is for the EOF analysis. The regression part can be
% blocked out to get the same analysis as figures 3 and 4
clear
cd /home/disk/hayes2/dvimont/csm/300yr/srf
load sst151-300.mat
load oro151-300.mat
load lims1-150.mat
whos
lims = [0 360 -90 90];
%sst = sst(360:1800,:,:);
%oro = oro(360:1800,:,:);
sst = subset(sst, nxa, nya, lims);
oro = subset(oro, nxa, nya, lims);
nxkeep = nxa(find(nxa >= lims(1) & nxa <= lims(2)));
nykeep = nya(find(nya >= lims(3) & nya <= lims(4)));
[ntmstp, nypts, nxpts] = size(sst);
whos
sst = reshape(sst, ntmstp, nypts*nxpts);
oro = reshape(oro, ntmstp, nypts*nxpts);
oro = mean(oro);
dokeep = find(oro == 0);
sst = sst(:,dokeep);
%
% Alter sst for ice points. Notice, the annual cycle feature was
% not included for figures 1 and 2 in ~dvimont/pacific
%
noelim = [];
for i = 1:length(dokeep);
noelimtmp = [];
for j = 1:12
icept = (find(sst(j:12:ntmstp,i) <= 273.14)-1).*12 + j;
noice = (find(sst(j:12:ntmstp,i) >= 273.14)-1).*12 + j;
if (~isempty(icept) & ~isempty(noice));
sst(icept,i) = mean(sst(noice,i)) * ones(length(icept),1);
end
if ~isempty(noice);
noelimtmp = [noelimtmp j];
else
sst(icept,i) = mean(sst(icept,i)) * ones(length(icept),1);
end
end
if ~isempty(noelimtmp);
noelim = [noelim i];
end
end
sst = sst(:,noelim);
dokeep = dokeep(noelim);
% Remove any point that ever had ice
%noelim = [];
%for i = 1:length(dokeep);
% icept = find(sst(1:ntmstp,i) <= 273.14);
% if isempty(icept);
% noelim = [noelim i];
% end
%end
%sst = sst(:,noelim);
%dokeep = dokeep(noelim);
whos
%
% Remove linear trend from each point's time series
%
%x = 1:ntmstp;
%xprime = x - mean(x);
%a1 = xprime * (sst - ones(ntmstp,1)*mean(sst)) / (xprime * xprime');
%a0 = mean(sst) - a1*mean(x);
%yhat = ones(ntmstp,1)*a0 + x'*a1;
%sst = sst - yhat;
%trendpat = zeros(1, length(dokeep));
%for i = 1:length(dokeep);
% yhat(:,i) = yhat(:,i) - mean(yhat(:,i)) * ones(size(yhat(:,i)));
% yhat(:,i) = yhat(:,i) / sqrt(yhat(:,i)' * yhat(:,i));
% trendpat(1,i) = yhat(:,i)' * (sst(:,i) - mean(sst(:,i)) * ones(size(sst(:,i))));
%end
%temptrend = zeros(1, nypts*nxpts);
%temptrend(1,dokeep) = trendpat;
%gcont(squeeze(reshape(temptrend(1,:),1,nypts,nxpts)),lims)
%clear yhat;
% Remove climatology
[sst, clim] = annave(sst);
sstemp = zeros(ntmstp, nxpts*nypts);
sstemp(:,dokeep) = sst;
sst = sstemp;
clear sstemp
sst = cosweight(sst, nya);
sst = sst(:, dokeep);
% Running Mean
sst3yr = myrunning_ave(sst, 36);
% This might make things faster
sst3yr = sst3yr(1:3:ntmstp,:);
% Take EOF's. (Might alter the dimensions if need be...)
biff = ['starting correllation']
c = sst3yr * sst3yr' / (length(dokeep));
biff = ['starting EOF']
[lam, eof, per] = eof(c);
pc1 = eof(:,1:10);
eof1 = pc1'*sst3yr;
eof1 = eof1' * diag(1 ./ sqrt(diag(eof1*eof1'/length(dokeep))));
eof1 = eof1';
rescale = std(pc1);
pc1 = pc1 * diag(1./rescale);
eof1 = diag(rescale) * eof1;
temp = zeros(10, nypts*nxpts);
temp(1:10,dokeep) = eof1(1:10,:);
for i = 1:length(dokeep);
a(i) = 0.5*(corr(sst(:,i),sst(:,i),1) + sqrt(abs(corr(sst(:,i),sst(:,i),2))));
end
a = mean(a);
dof = -1 * 600 * log(a) / 2;
dlam = lam * sqrt(2/dof);