clear
clc

%% slp
f='/Volumes/Backup Plus/backup/polar/ERA5/surface7922.nc';
for m=27:38
    mclim(:,:,m-26)=mean(squeeze(ncread(f,'msl',[1 1 1 m],[1440 361 1 30],[1 1 1 12])/100),3);
end
x0=ncread(f,'longitude');
y0=ncread(f,'latitude');
[X,Y]=meshgrid(x0,y0);
X=double(X);Y=double(Y);
ma17spr(:,:,1)=permute(squeeze(ncread(f,'msl',[1 1 1 453],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,7),[2 1]);
ma17spr(:,:,2)=permute(squeeze(ncread(f,'msl',[1 1 1 454],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,8),[2 1]);
ma17spr(:,:,3)=permute(squeeze(ncread(f,'msl',[1 1 1 455],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,9),[2 1]);
ma22spr(:,:,1)=permute(squeeze(ncread(f,'msl',[1 1 1 513],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,7),[2 1]);
ma22spr(:,:,2)=permute(squeeze(ncread(f,'msl',[1 1 1 514],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,8),[2 1]);
ma22spr(:,:,3)=permute(squeeze(ncread(f,'msl',[1 1 1 515],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,9),[2 1]);

ma17sum(:,:,1)=permute(squeeze(ncread(f,'msl',[1 1 1 456],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,10),[2 1]);
ma17sum(:,:,2)=permute(squeeze(ncread(f,'msl',[1 1 1 457],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,11),[2 1]);
ma17sum(:,:,3)=permute(squeeze(ncread(f,'msl',[1 1 1 458],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,12),[2 1]);
ma22sum(:,:,1)=permute(squeeze(ncread(f,'msl',[1 1 1 516],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,10),[2 1]);
ma22sum(:,:,2)=permute(squeeze(ncread(f,'msl',[1 1 1 517],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,11),[2 1]);
ma22sum(:,:,3)=permute(squeeze(ncread(f,'msl',[1 1 1 518],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,12),[2 1]);

ma23spr(:,:,1)=permute(squeeze(ncread(f,'msl',[1 1 1 525],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,7),[2 1]);
ma23spr(:,:,2)=permute(squeeze(ncread(f,'msl',[1 1 2 526],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,8),[2 1]);
ma23spr(:,:,3)=permute(squeeze(ncread(f,'msl',[1 1 2 527],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,9),[2 1]);
f0='/Volumes/Backup Plus/backup/polar/ERA5/surface23.nc';
ma23sum(:,:,1)=permute(squeeze(ncread(f,'msl',[1 1 2 528],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,10),[2 1]);
ma23sum(:,:,2)=permute(squeeze(ncread(f0,'msl',[1 1 1 1],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,11),[2 1]);
ma23sum(:,:,3)=permute(squeeze(ncread(f0,'msl',[1 1 1 2],[1440 361 1 1],[1 1 1 1])/100)-mclim(:,:,12),[2 1]);

mslspr17=mean(ma17spr,3,"omitmissing");
mslspr22=mean(ma22spr,3,"omitmissing");
mslspr23=mean(ma23spr,3,"omitmissing");
mslsum17=mean(ma17sum,3,"omitmissing");
mslsum22=mean(ma22sum,3,"omitmissing");
mslsum23=mean(ma23sum,3,"omitmissing");

%% SIC
load lon
load lat
load ntSIC
for i=27:38
    Cclim(:,:,i-26)=mean(C(:,:,i:12:386),3,"omitmissing");
end
filepath='/Volumes/Backup Plus/backup/polar/dailySICv2_nrt/';
filelist=dir([filepath,'N*.nc']);
k=length(filelist);
for i=1:167
  filename=[filepath,filelist(i).name];
  SIC=ncread(filename,'F18_ICECON')*100;
  SIC(SIC>100)=nan;
  C0(:,:,i)=SIC;
end

Ca17spr(:,:,1)=C(:,:,453)-Cclim(:,:,7);
Ca17spr(:,:,2)=C(:,:,454)-Cclim(:,:,8);
Ca17spr(:,:,3)=C(:,:,455)-Cclim(:,:,9);
Ca17sum(:,:,1)=C(:,:,456)-Cclim(:,:,10);
Ca17sum(:,:,2)=C(:,:,457)-Cclim(:,:,11);
Ca17sum(:,:,3)=C(:,:,458)-Cclim(:,:,12);

Ca22spr(:,:,1)=C(:,:,513)-Cclim(:,:,7);
Ca22spr(:,:,2)=C(:,:,514)-Cclim(:,:,8);
Ca22spr(:,:,3)=C(:,:,515)-Cclim(:,:,9);
Ca22sum(:,:,1)=C(:,:,516)-Cclim(:,:,10);
Ca22sum(:,:,2)=C(:,:,517)-Cclim(:,:,11);
Ca22sum(:,:,3)=C(:,:,518)-Cclim(:,:,12);

Ca23spr(:,:,1)=C(:,:,525)-Cclim(:,:,7);
Ca23spr(:,:,2)=C(:,:,526)-Cclim(:,:,8);
Ca23spr(:,:,3)=C(:,:,527)-Cclim(:,:,9);
Ca23sum(:,:,1)=C(:,:,528)-Cclim(:,:,10);
Ca23sum(:,:,2)=mean(C0(:,:,1:31),3,"omitmissing")-Cclim(:,:,11);
Ca23sum(:,:,3)=mean(C0(:,:,32:59),3,"omitmissing")-Cclim(:,:,12);

Cspr17=mean(Ca17spr,3,"omitmissing");
Cspr22=mean(Ca22spr,3,"omitmissing");
Cspr23=mean(Ca23spr,3,"omitmissing");
Csum17=mean(Ca17sum,3,"omitmissing");
Csum22=mean(Ca22sum,3,"omitmissing");
Csum23=mean(Ca23sum,3,"omitmissing");

%% plot
set(gcf,'Position',[300 400 800 400]);
ha = tight_subplot(1,2,[.03 0.05],[.12 .04],[.05 0.06]);
axes(ha(1))
m_proj('azimuthal equal-area','latitude',-90,'radius',55,'rectbox','on');
hold on;
m_pcolor(X,Y,mslsum22);
colortable = textread('/Users/jinfei/Documents/MATLAB/colorbar/cmp_flux.txt');
colormap(ha(1),colortable);
m_gshhs_i('color','k');
m_grid('xaxislocation','top','xtick',[-120 -60 0 60 120 180],'ytick',[-90 -80 -70 -60],'fontsize',11);
clim([-15 15]);
% hcb1 = colorbar('southoutside','FontSize',14');
% set(hcb1,'position',[0.11 0.05 0.3 0.04],'YTick',[-15 0 15],'FontSize',12)
% t1=get(hcb1,'Yticklabel');
% t1=strcat(t1,'hPa');
% set(hcb1,'Yticklabel',t1);

axes(ha(2))
m_proj('azimuthal equal-area','latitude',-90,'radius',55,'rectbox','on');
hold on;
m_gshhs_i('patch',[.8 .8 .8],'edgecolor',[.8 .8 .8],'MarkerSize',0.1);
m_pcolor(lon,lat,Csum22);
colortable = textread('/Users/jinfei/Documents/MATLAB/colorbar/sunshine_diff_12lev.txt');
colormap(ha(2),colortable);
m_grid('xaxislocation','top','xtick',[-120 -60 0 60 120 180],'ytick',[-90 -80 -70 -60],'fontsize',11);
clim([-50 50]);
% hcb3 = colorbar('southoutside','FontSize',14);
% set(hcb3,'position',[0.58 0.05 0.3 0.04],'YTick',[-50 0 50],'FontSize',12)
% t1=get(hcb3,'Yticklabel');
% t1=strcat(t1,'%');
% set(hcb3,'Yticklabel',t1);