
%
% historical run
CNAM='NBF20TROC_f19_tn11_01';
YEAR1=1976;
YEARN=2005;
SEAS='ann';

% PI run
CNAM='NBF1850OC_f19_tn11_02';
YEAR1=1801;
YEARN=2000;
SEAS='ann';


% mis3 run
CNAM='NBF1850_f19_tn11_test_mis3b';
YEAR1=1801;
YEARN=2000;
SEAS='ann';


switch SEAS
    case 'ann'
        mon1=1; mon2 =12;
    case 'djf'
        mon1=1; mon2 =3;
end


%
fname=['Data/' CNAM '.cam2.h0.' num2str(YEAR1) '_' num2str(YEARN) 'y_' SEAS '.nc'];
fname_eof=['Data/' CNAM '.cam2.h0.' num2str(YEAR1) '_' num2str(YEARN) 'y_' SEAS '_eof_AO.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,'PSL')/100; % Pa-->hPa

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

% get the monthly anomalies
data = zeros(ny,nx,nr);
nn=0;
for iy=YEAR1:YEARN
  for im=mon1:mon2
   nn=nn+1;
   irec = (iy-YEAR1)*nmon + 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,:);

% Preparing data
[d1,d2,d3]=size(data1);
data1=reshape(data1,d1,d2*d3);
%[pcs,eofs,v,l,w]=eof(data,areaatm);

% EOF
[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;

%caxis([-6 6])
% [xi2,yi2]=m_ll2xy(xi,yi);
% [C,H] = contourf(xi2',yi2',eof1',vDr);
% [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)

 [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(['NorESM ('  num2str(round(var(imod)*100), '%2.2d') '%) ' upper(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




