% Seasonal cycle analysis of CMIP3 models, reanalysis data and IMD rainfall
% April 2009

clear
%% IMD data
load /Users/andy/data/India/IMD/IMD-Daily1DegRain-1901-2004/IMD-Daily1DegRain-1901-2004.mat
nx=35; ny=33; ndays=dayN-day1+1;
A=reshape(A,ndays,ny*nx);
x_axis=X(1,:)'; y_axis=Y(:,1);
X=reshape(X,1,ny*nx); Y=reshape(Y,1,ny*nx);
dvec=datevec(day1:dayN);

% select peninsular India, land points
% xysel=find(Y<=30 & X<=88 & A(end,:)~=-999);
xysel=find(A(end,:)~=-999);
A1=A(:,xysel);
X1=Y(xysel); Y1=Y(xysel);
nxy=length(xysel);

% Monthly averages,averaged over years 1971-2000
nmon=12;
rain=NaN*ones(nmon,nxy);
imon=0; 
for mon=1:12 
  imon=imon+1;
  dsel=find(dvec(:,2)==mon & dvec(:,1)>=1971 & dvec(:,1)<=2000);
  rain(imon,:)=mean(A1(dsel,:));
end
clear A A1

%% plot 12 monthly maps
smon={'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'};
figure(1); clf
for mon=1:12
  subplot(4,3,mon)
  map=NaN*ones(ny,nx);
  map(xysel)=rain(mon,:);
  plot_map2(x_axis,y_axis,map,[0 15])
  title(smon{mon})
end
  
%% EOF analysis via SVD: 
% [U,S,V]=svd(rain) where U are the PCs, V' are the EOFs, S**2/(nmon-1) 
% are the variances

rains=rain-repmat(mean(rain),nmon,1); % remove time mean
[U,S,V]=svd(rains,'econ');
PCvar=diag(S).^2/sum(diag(S).^2);

% standardize the PCs to unit variance
U=U(:,1:11); % 12th PC is a constant 
U=U./repmat(std(U),nmon,1);

% weight PCs by sqrt(PCvar)
U=U.*repmat(sqrt(PCvar(1:11))',nmon,1);

%% plot PCs
figure(2); clf
plot(U(:,1:4))
legend('PC1','PC2','PC3','PC4')
grid on
title('IMD Seasonal cycle PCs')

%% plot EOFs
figure(3); clf
for k=1:4
  subplot(2,2,k)
  map=NaN*ones(ny,nx);
  map(xysel)=V(:,k);
  plot_map2(x_axis,y_axis,map,[-0.2 0.2])
  title(['EOF ',num2str(k)])
end
