
%
CNAM='NCEP';
YEAR1=1979;
YEARN=2008;
SEAS='ANN';
switch SEAS
    case 'ANN'
        mon1=1; mon2 =12;
    case 'DJF'
        mon1=12; mon2 =14;
    case 'DJFM'
        mon1=11; mon2 =14;
    case 'JFM'
        mon1=11; mon2 =14;
    case 'JJA'
        mon1=6; mon2=8;
end


%
fname=['Data/mslp.mon.mean.nc'];
fname_eof=['Data/mslp.mon.mean_eof_AAO.mat'];
if exist(fname_eof);
%   load(fname_eof);
%   disp('Reading from existing data...')
%else

%
lon=ncread(fname,'lon');
lat=ncread(fname,'lat');
[xi yi] = meshgrid(lon,lat);
[nx ny] = size(xi);
dataa=ncread(fname,'mslp')/100; % Pa-->hPa

% calculate climatology
dataa_clim=zeros(ny,nx,12);
nr=YEARN-YEAR1+1;
for mon=1:12
 data_clim(:,:,mon)=mean(dataa(:,:,mon:12:nr*12),3);
end

% get the monthly anomalies
data = zeros(ny,nx,nr);
nn=0;
for iy=YEAR1:YEARN
  for im=1:12
   nn=nn+1;
   irec = (iy-YEAR1)*12 + im;
   data(:,:,nn)=dataa(:,:,irec)-data_clim(:,:,im);
  end
end
data=permute(data,[3 2 1]);

% for data north of 20N
indx=find(lat<-20);
lat1=lat(indx);
lon1=lon;
data1=data(:,indx,:);

% get EOF
[d1,d2,d3]=size(data1);
data1=reshape(data1,d1,d2*d3);
[pc,eof,var]=is_svdcov(lon1,lat1,data1,2,3,'linear');

% NAM/AO
pc1=-pc(1,:);
var1=var(1);

%
save(fname_eof,  'pc', 'eof','var','xi','yi','data');

end



%
% ---------------- PLOT ----------------
%

% Preparing data
PSL1 = data;
[d1,d2,d3]=size(PSL1);
PSL1=reshape(PSL1,d1,d2*d3);

% Remove mean from time series
PSL1=detrend(PSL1,'constant');
PSL1=reshape(PSL1,d1,d2,d3);

vDr = [-25:5:25]/5;

% calculate for the first few EOFs 
for imod =1:1
    eof1=ones(d2, d3);
    for i=1:d2
      for j=1:d3
        eof1(i,j)=squeeze(PSL1(:,i,j))'*(pc(imod,:))';
      end
    end
    eof1 =eof1/(d1-1);
    pc1 = pc(1,:);
    pc2 = pc(2,:);

    disp(['Range: ' num2str(min(min(eof1)),'%5.2f') ' to ' num2str(max(max(eof1)),'%5.2f') ] );

    [nx ny] = size(xi);
    xi(:,ny+1)=xi(:,1);
    yi(:,ny+1)=yi(:,1);
    eof1(:,ny+1)=eof1(:,1);

figure;
m_proj('Stereographic','lon',0,'lat',-90,'rad',70);
%m_pcolor(xi,yi,eof1);shading flat;

 %[C,H] = m_contourf(xi',yi',eof1',vDr);
 %set (H, 'LineWidth' ,8, 'EdgeColor', 'none')
 %clabel(C,H,'FontSize',12,'FontWeight','Bold','Color','b','Rotation',0)
 %caxis([-8 8])
 %colorbar;

 [C,H] = m_contour(xi,yi,eof1,vDr(find(vDr<0)), 'b--');
 set (H, 'LineWidth' ,2, 'EdgeColor', 'b')
 clabel(C,H,'FontSize',12,'FontWeight','Bold','Color','b','Rotation',0)
% 
 hold on
 [C,H] = m_contour(xi,yi,eof1,vDr(find(vDr>0)), 'r-');
 set (H, 'LineWidth' ,2, 'EdgeColor', 'r')
 clabel(C,H,'FontSize',12,'FontWeight','Bold','Color','b','Rotation',0)


m_grid('xtick',[],'ytick',[],'xticklabel',[],'yticklabel',[]);
m_coast('color',[.4 .4 .4],'linewidth',1);


title(['NCEP ('  num2str(round(var(imod)*100), '%2.2d') '%) ' SEAS], 'FontSize',16,'FontWeight','Bold')

print('-depsc',['Figs/' CNAM '_eof_PSL_AAO_' num2str(YEAR1) '-' num2str(YEARN) '.eps']);
print('-dpng','-r300',['Figs/' CNAM '_eof_PSL_AAO_' num2str(YEAR1) '-' num2str(YEARN) '.png']);
end




