% Jul 31, 2001
% MK Tippett
% matlab script to read CCM data from DL and compute S2N eofs.
%
m=4; %number of ensemble members
nt=204; %number of time points
nx=128; %number of x-points
ny=64; %number of y-points
ens=zeros(ny,nx,nt,m);

loaddods('http://ingrid.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.IRI/.MP/.RESEARCH/.ATMOSPHERE/.GLOBAL/.CCM/.mem01/.ORO/0/maskgt%5Btime%5Daverage/dods')

mask=~isnan(ORO(:));


% read data
ldstr1='http://ingrid.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.IRI/.MP/.RESEARCH/.ATMOSPHERE/.GLOBAL/.CCM/.mem0';
ldstr2='/time/%28Jan%201979%29%28Dec%201995%29RANGE/.TAUX/dods'
for ii=1:m
  loaddods([ldstr1,int2str(ii),ldstr2])
  ens(:,:,:,ii)=TAUX+ORO(:,:,ones(1,nt));
end

% compute ensmeble mean and deviations
ens_mean = mean(ens,4);
ens = ens - ens_mean(:,:,:,ones(1,m));

% compute yearly-anomalies
W = reshape(ens_mean,ny*nx,nt);
W_mean=repmat(NaN,ny*nx,12);
Wa=repmat(NaN,ny*nx,nt);
for ii=1:12
  W_mean(mask,ii)=mean(W(mask,ii:12:nt),2);
end
for ii=1:12
  ind = ii:12:nt;
  Wa(mask,ind) = W(mask,ind) - W_mean(mask,ii)*ones(1,length(ind));
end
ens_meana = reshape(Wa,ny,nx,nt);

[lx,ly,lt,m] = size(ens);


nmodes = 100; % number of modes to keep.
pmodes = 10;  % number of modes to plot

% compute ensemble mean covariance

W = reshape(ens_meana,nx*ny,nt)./sqrt(nt-1);
k=1:nmodes;
ts_covar;
D.ens=diag(d);
ts.ens=time_series;
Z.ens = z(:,1:pmodes);
clear W z d time_series

pcolor(lon,lat,reshape(Z.ens(:,1),ny,nx));shading('flat')
drawnow

% compute inter-ensemble covariance

W = reshape(ens,nx*ny,nt*m);
W = W./sqrt(nt*m-1);
k=1:nmodes;
ts_covar;
D.inter_ens=diag(d);
ts.inter_ens=time_series;
Z.inter_ens = z(:,1:pmodes);

pcolor(lon,lat,reshape(Z.inter_ens(:,1),ny,nx));shading('flat')
drawnow


%F = z(mask(:),:)*diag(sqrt(m)./sqrt(diag(d)));
F = z(mask(:),:)*diag(1./sqrt(diag(d)));
tmp = reshape(ens_meana,nx*ny,nt);
W = F'*tmp(mask(:),:);
k=1:nmodes;
ts_covar;
D.sn=diag(d);
ts.sn=time_series;
Z.sn = z(:,1:pmodes);

flt=zeros(nx*ny,nmodes);
flt(mask(:),:) = F*z; flt(~mask(:),:)=NaN;

regres = inv(time_series*time_series')*time_series;
sn = reshape(ens_meana,nx*ny,nt)*regres';

pcolor(lon,lat,reshape(sn(:,1),ny,nx));shading('flat')
drawnow

signal = reshape(sn(:,1:pmodes),nx*ny,pmodes)*time_series(1:pmodes,:);
pcolor(lon,lat,reshape(signal(:,1),ny,nx));shading('flat')
drawnow

clear W








