clear
clc


f1='/Volumes/Backup Plus/PAROCE/siconc_all_icemod_1980_2023.nc';% 1980.1-2023.2
sic=ncread(f1,'siconc')*100; 
sic22=sic(:,:,505:518);
sic22(sic22<5)=nan;

f2='/Volumes/Backup Plus/PAROCE/snvolu_all_icemod_1980_2023.nc';% 1980.1-2023.2
sn=ncread(f2,'snvolu'); 
sn22=sn(:,:,505:518);
sn22(isnan(sic22)==1)=nan;
for i=13:26
    snclim(:,:,i-12)=nanmean(sn(:,:,i:12:374),3);
end
snclim(isnan(sic22)==1)=nan;
lon=ncread(f1,'nav_lon');
lat=ncread(f1,'nav_lat');

% set(gcf,'Position',[350 200 1000 500]);
% ha = tight_subplot(2,4,[0.05 0.02],[.04 .04],[.02 .02]);
% for i=1:7
%     axes(ha(i))
%     m_proj('azimuthal equal-area','latitude',-90,'radius',55,'rectbox','on');
%     hold on;
%     m_pcolor(lon,lat,snano(:,:,i));
%     m_gshhs_i('color','k');
%     m_grid('xaxislocation','top','xticklabels',[],'yticklabels',[],'fontsize',8);
%     colortable = textread('/Users/jinfei/Documents/MATLAB/colorbar/ncl_default.txt');
%     colormap(colortable);
%     caxis([-0.25 0.25]);
% end
% delete(ha(8))
% hcb1 = colorbar('FontSize',12,'YTick',[-0.25 0 0.25]);
% set(hcb1,'position',[0.765 0.055 0.02 0.4])
% title(hcb1,'m');

ind=find(lon>-180&lon<-120&lat>-80&lat<-70);
for i=1:14
    sn220=sn22(:,:,i);
    snclim0=snclim(:,:,i);
    sn22t(i)=nanmean(sn220(ind));
    snct(i)=nanmean(snclim0(ind));
end
snano=sn22t-snct;

plot(snano,'k','LineWidth',1.5)
set(gca,'xtick',1:2:14,'xticklabel',{'Jan' 'Mar' 'May' 'July' 'Sep' 'Nov' 'Jan'},'fontname','Arial','fontsize',14);
axis([1 14 -0.15 0])
ylabel('Snow depth anomaly (m)','fontname','Arial','fontsize',16)

