clear all

% read data 
IDIR='/cluster/shared/noresm/inputdata/atm/cam/ozone';
%
FPATH=[IDIR '/tracer_cnst_CAM6chem_piControl_cycle_3D_monthly_c170901.nc'];
ozone1=flip(mean(ncread(FPATH,'O3'),4),3);
ps1=mean(ncread(FPATH,'PS'),3);
lon1=ncread(FPATH,'lon');
lat1=ncread(FPATH,'lat');
lev1=flipud(ncread(FPATH,'lev')); 
hyam1=flipud(ncread(FPATH,'hyam')); 
hybm1=flipud(ncread(FPATH,'hybm')); 
%
FPATH=[IDIR '/vmro3_input4MIPs_ozone_CMIP_UReading-CCMI-1-0_gn_185001-185012.nc']; 
ozone2=mean(ncread(FPATH,'vmro3'),4); 
lon2=ncread(FPATH,'lon'); 
lat2=ncread(FPATH,'lat'); 
plev2=ncread(FPATH,'plev'); 

% vertical interpolation 
plev1=zeros(size(ozone1)); 
for k=1:numel(lev1) 
 plev1(:,:,k)=1000*hyam1(k)+ps1/100*hybm1(k);  
end
plev1=reshape(plev1,numel(lon1)*numel(lat1),[])';
ozone1=reshape(ozone1,numel(lon1)*numel(lat1),[])';
for n=1:size(ozone1,2) 
  n
  ozone1on2(:,n)=interp1(plev1(:,n),ozone1(:,n),plev2); 
end
%ozone1on2=interp1(plev2,ozone1on2,plev2,'nearest','extrap'); 
ozone1on2=reshape(ozone1on2',numel(lon1),numel(lat1),numel(plev2));

% compute zonal averages 
ozone1zon=squeeze(nanmean(ozone1on2));
ozone2zon=squeeze(mean(ozone2));

% interpolate in lat direction
ozone1zon=interp1(lat1,ozone1zon,lat2);

% plot 
if true
figure(1)
clf
subplot(1,3,1) 
h=pcolor(lat2,plev2,ozone1zon');
set(h,'edgecolor','none')
colormap(jet)
caxis([0 1e-5])
set(gca,'ydir','reverse','yscale','log','ylim',[5 600],...
'ytick',[5 10 20 50 100 200 500 1000],'yticklabel',num2str([5 10 20 50 100 200 500 1000]'),'xtick',[-90:30:90])
colorbar('southoutside')
title('NCAR-CAM6')

subplot(1,3,2) 
h=pcolor(lat2,plev2,ozone2zon');
set(h,'edgecolor','none')
caxis([0 1e-5])
colormap(jet)
set(gca,'ydir','reverse','yscale','log','ylim',[5 600],...
'ytick',[5 10 20 50 100 200 500 1000],'yticklabel',num2str([5 10 20 50 100 200 500 1000]'),'xtick',[-90:30:90])
colorbar('southoutside')
title('CMIP6')

subplot(1,3,3)
h=pcolor(lat2,plev2,ozone1zon'-ozone2zon');
set(h,'edgecolor','none')
%cmap=[fliplr(pink);flipud(pink)];
%colormap(cmap)
caxis([0 1e-5])
set(gca,'ydir','reverse','yscale','log','ylim',[5 600],...
'ytick',[5 10 20 50 100 200 500 1000],'yticklabel',num2str([5 10 20 50 100 200 500 1000]'),'xtick',[-90:30:90])
colorbar('southoutside')
title('NCAR-CAM6 - CMIP6')
print(gcf,'-dpng','-r300','o3_ncar-cam6_vs_cmip6_a.png') 

end 

if false
figure(2)
cmap=[fliplr(pink);flipud(pink)];
colormap(cmap)
clf
clev=[0:0.1:2]*1e-5; 
clevdiff=[-2:0.1:2]*1e-6;
contourf(lat2,plev2,ozone1zon'-ozone2zon',clevdiff,'edgecolor','none'); 
hold on 
[c,h]=contour(lat2,plev2,ozone1zon',clev,'color',[0 1 0]*1,'linewidth',1.5);
%clabel(c,h);
[c,h]=contour(lat2,plev2,ozone2zon',clev,'color',[1 1 1]*0.3,'linewidth',1.5);
clabel(c,h);
hold off 
caxis([-2 2]*1e-6)
set(gca,'ydir','reverse','yscale','log','ylim',[5 600],...
'ytick',[5 10 20 50 100 200 500 1000],'yticklabel',num2str([5 10 20 50 100 200 500 1000]'),'xtick',[-90:30:90])
colorbar('southoutside')
ylabel('Pressure (hPa)')
xlabel('Latitude (degN)') 
title('Mole fraction of ozone in air - NCAR-CAM6 (green) less CMIP6 (black)')
print(gcf,'-dpng','-r300','o3_ncar-cam6_vs_cmip6_b.png') 
end
