%function plot_atmos2d_diff(VNAM,CNAM1,CNAM2,YEAR1,YEARN) 

VNAME='TREFHT';
%VNAME='PSL';
%VNAME='PRECT';
%VNAME='TOA';
%VNAME='CLDTOT';
%VNAME='FLNT';

YEAR1=1780;
YEARN=1799;
CNAM1='NBF1850_f19_tn11_test_mis3b_fwf3b';
CNAM2='NBF1850_f19_tn11_test_mis3b';

ODIR='Figs_atmos2d_diff';
OPATH=[ODIR '/atmos_' VNAME '_' CNAM1 '-' CNAM2 sprintf('_%4.4i-%4.4i_sterographic.png',YEAR1,YEARN)];

flagtoa=false
switch VNAME
case 'TREFHT'
  CVEC=[-6:0.5:6];
  CVEC=[-24:4:-4 -2 -1 0 1 2 4:4:24];
  CVEC=[-20 -15 -10 -5 -2 0 2 5 10 15 20];
  LNAME='2m air temperature';
  UNITS='K';
  FAC=1;
case 'PSL'
  CVEC=[-1:0.1:1];
  LNAME='Sea level pressure';
  UNITS='hPa';
  FAC=1e-2;
case 'PRECT'
  CVEC=[-1:0.1:1];
  LNAME='Total precipitation';
  UNITS='mm/d';
  FAC=3600*24*1000;
case 'TOA'
  CVEC=[-1:0.1:1]*2;
  LNAME='TAO net radiation';
  UNITS='W/m2';
  FAC=1;
  flagtoa=true;
case 'CLDTOT'
case 'RHREFHT'
case 'U200'
  CVEC=[-20:1:20];
  LNAME='U_200 speed';
  UNITS='m/s';
  FAC=1;

end 


addpath ~/matlab/mats
addpath ~/matlab/m_map
addpath ~/matlab/ingo
addpath ~/matlab/cmap

% read  data
area=ncread('../../Data/areaatm.nc','areacella');
weights=area/sum(area(:));
if flagtoa
  fld=(ncread(IPATH,'FSNT')-ncread(IPATH,'FLNT'))*FAC;
else
  load(['Data/atmos2d_ANN_' VNAME '_' CNAM1 '_' sprintf('%4.4i',YEAR1) '-' sprintf('%4.4i',YEARN) '.mat']);ff1=data;
  load(['Data/atmos2d_ANN_' VNAME '_' CNAM2 '_' sprintf('%4.4i',YEAR1) '-' sprintf('%4.4i',YEARN) '.mat']);ff2=data;
  fld=(ff1-ff2)*FAC;
end
fldglb=sum(fld(:).*weights(:));
fldglb=floor(sign(fldglb)*fldglb*100)/100*sign(fldglb)
%fld=cat(1,fld(find(lon>=30),:),fld(find(lon<30),:));
%lon=cat(1,lon(find(lon>=30)),360+lon(find(lon<30)));

mkdir(ODIR);

% plot 
[nx ny]=size(fld);
lon=lon*ones(1,ny);
lat=ones(nx,1)*lat';
maskmv=ones(size(fld));
maskmv(find(~isnan(fld)))=0;
fld=extrap(fld);
figure_height_scale=8/6;
cbar_width_scale=1;
fontsize=12;
cv=CVEC;
%
figure
clf
%m_proj('Equidistant Cylindrical','lon',[30 390],'lat',[-90 90]);
m_proj('stereographic','lat',90,'long',0,'radius',60);
set(gcf,'paperunits','inches','papersize',[8 6*figure_height_scale], ...
        'paperposition',[0 0 8 6*figure_height_scale],'renderer','painters')
set(gca,'outerposition',[0 0 1 1],'position',[0.1 0.1 0.8 0.75])
colormap(cbsafemap(511,'rdbu'))
[c,h]=m_contourf(lon,lat,-fld,[-cv(1) inf],'linecolor','none');
hold on
set(findobj(h,'type','patch'),'facecolor','flat','cdata',-inf);
[c,hc]=m_contourf(lon,lat,fld,cv,'linecolor','none');
%m_contour(lon,lat,fld,[0 0],'linecolor',[1 1 1]*0.6,'linewidth',3);
[c,h]=m_contourf(lon,lat,maskmv,[.5 inf],'linecolor','none');
%[c10,h10]=m_contour(lon,lat,fld,[10 10]);
%set(h10,'color','g','linewidth',3)
patchhole(findobj(h,'type','patch'),[1 1 1]*0.4);
m_grid
m_coast('color',[1 1 1]*0.01);
set(gca,'DataAspectRatioMode','auto')
hb=cbarfmb([-inf inf],cv,'vertical','nonlinear');
ytick=get(hb,'ytick');
yticklabel=get(hb,'yticklabel');
pos=get(hb,'position');
pos(3)=pos(3)*cbar_width_scale;
set(hb,'position',pos,'fontsize',fontsize)
cdatamidlev([hc;hb],cv,'nonlinear');
caxis([cv(1) cv(end)])
caxis(hb,[cv(1) cv(end)])
title([LNAME ' difference in ' UNITS sprintf(' (global=%5.2f)',fldglb) ],'fontsize',14)
print(gcf,'-dpng','-r500',OPATH);

