clear
clc

%% a
f='/Users/jinfei/Documents/polar/2023extreme/S_seaice_extent_daily_v3.0.csv';
SIE=readtable(f);
SIE23=str2double(table2array(SIE(14183:14547,4))); %2022.3.1-2023.2.28
f1='/Users/jinfei/Documents/polar/2023extreme/S_seaice_extent_climatology_1981-2010_v3.0.csv';
SIE0=readtable(f1);
SIEc1=table2array(SIE0(1:59,2));%1.1-2.28
SIEc2=table2array(SIE0(61:366,2));%3.1-12.31
SIEc=[SIEc2;SIEc1];

SIEd=str2double(table2array(SIE(65:14182,1:4)));%3.1-2.28
year=SIEd(:,1);mon=SIEd(:,2);day=SIEd(:,3);
SIEd(mon==2&day==29,:)=[];
SIEd=[SIEd(1:1651,:);nan(41,4);SIEd(1652:end,:)];
for i=1:1546
    SI(2*i-1,:)=SIEd(i,:); 
    SI(2*i,:)=SIEd(i,:); 
end
for i=3093:15695
    SI(i,:)=SIEd(i-1546,:);
end
for i=1:43
   SIE00{i}=SI(365*i-364:365*i,4); 
end
% min&max
for i=1:365
    ma(i)=max(SI(i:365:15695,4));
    mi(i)=min(SI(i:365:15695,4));
end
% for i=1:365
%     ddp=SI(i:365:15695,4);
%     sd(i)=nanstd(ddp);
% end
% sd=sd';
%% a

set(gcf,'Position',[500 400 750 450]);
day=1:365;
p=fill([day fliplr(day)],[mi';flipud(ma')],[0.8 0.8 0.8],'edgealpha', '0', 'facealpha', '.5');
hold on
% grid on
% axis square
ylabel('SIE (million km^2)','fontname','Arial','fontsize',14)
% xlabel('Date','fontname','Arial','fontsize',12)
axis([1 365 0 20.5])
p2=plot(day,SIEc,'-','color','k','LineWidth',1.5);
p3=plot(day,cell2mat(SIE00(43)),'color',sixteen2ten('#061BEE')/255,'LineWidth',1.5);
p4=plot(day,cell2mat(SIE00(39)),'color',sixteen2ten('#075FE8')/255,'LineWidth',1.5);
p5=plot(day,cell2mat(SIE00(38)),'color',sixteen2ten('#089EE2')/255,'LineWidth',1.5);
p1=plot(day,SIE23,'-','color','r','LineWidth',2.5);
set(gca,'xtick',[1 32 62 93 123 154 185 215 246 276 307 338],'xticklabel',{'Mar' 'Apr' 'May' 'Jun' 'July' 'Aug' 'Sep' 'Oct' 'Nov' 'Dec' 'Jan' 'Feb'},'fontname','Arial','fontsize',14);
legend([p5 p4 p3 p1 p2],{'2016-2017','2017-2018','2021-2022','2022-2023','1981-2010 mean'},'position',[0.75 0.77 0.1 0.1],'fontname','Arial','fontsize',12);
text(12,19.5,"a",'FontSize',16,'FontWeight','bold','fontname','Arial')


%% b
f='/Volumes/Backup Plus/backup/polar/dailySICv2_nrt/NSIDC0081_SEAICE_PS_S25km_20230221_v2.0.nc';
SIC=ncread(f,'F18_ICECON')*100;
SIC(SIC>100)=nan;
load SIC_daily % start from 1979.2
n=1;
for j=1:30
for i=751
    C(:,:,n)=SIC_daily(:,:,i+365*(j-1));
    n=n+1;
end
end
CC=nanmean(C,3);

load lon
load lat
% rano=(SIC-CC)./CC;
% rano(CC<15)=nan;
rano=SIC-CC;

%% b
m_proj('azimuthal equal-area','latitude',-90,'radius',50,'rectbox','on');
hold on;
m_gshhs_i('patch',[.8 .8 .8],'edgecolor',[.8 .8 .8],'MarkerSize',0.1);
m_pcolor(lon,lat,rano);
m_contour(lon,lat,SIC,[15 15],'color',sixteen2ten('#4E23D1')/255,'LineWidth',1.5);
colortable = textread('/Users/jinfei/Documents/MATLAB/colorbar/MPL_bwr.txt');
colormap(flipud(colortable));
m_grid('xaxislocation','top','xtick',[-120 -60 0 60 120 180],'ytick',[-90 -80 -70 -60],'fontsize',11);
caxis([-50 50]);
c=colorbar;
set(c, 'Position', [0.88 0.15 0.03 0.7]);
title(c,'%','fontsize',14)
text(-0.58,0.56,"b",'FontSize',16,'FontWeight','bold','fontname','Arial')

%% c
s1=cell2mat(SIE00(38)); %2017
s2=cell2mat(SIE00(39)); %2018
s3=cell2mat(SIE00(43)); %2022
s4=SIE23; %2023
sc=SIEc; %clim

d1=smoothdata(diff(s1(1:365)),'gaussian',30);
d2=smoothdata(diff(s2(1:365)),'gaussian',30);
d3=smoothdata(diff(s3(1:365)),'gaussian',30);
d4=smoothdata(diff(s4(1:365)),'gaussian',30);
dc=smoothdata(diff(sc(1:365)),'gaussian',30);

d1(365)=nan;
d2(365)=nan;
d3(365)=nan;
d4(365)=nan;
dc(365)=nan;

%% c
set(gcf,'Position',[500 400 750 450]);
plot(1:365,zeros(1,365),'--','color','k','LineWidth',1)
hold on
ylabel('SIE tendency (million km^2/day)','fontname','Arial','fontsize',14)
% xlabel('Date','fontname','Arial','fontsize',12)
qc=plot(dc,'color','k','LineWidth',1.5);
q1=plot(d1,'color',sixteen2ten('#089EE2')/255,'LineWidth',1.5);
q2=plot(d2,'color',sixteen2ten('#075FE8')/255,'LineWidth',1.5);
q3=plot(d3,'color',sixteen2ten('#061BEE')/255,'LineWidth',1.5);
q4=plot(d4,'color','r','LineWidth',2.5);
set(gca,'xtick',[1 32 62 93 123 154 185 215 246 276 307 338],'xticklabel',{'Mar' 'Apr' 'May' 'Jun' 'July' 'Aug' 'Sep' 'Oct' 'Nov' 'Dec' 'Jan' 'Feb'},'fontname','Arial','fontsize',14);
axis([1 365 -0.32 0.15])
legend([q1 q2 q3 q4 qc],{'2016-2017','2017-2018','2021-2022','2022-2023','1981-2010 mean'},'position',[0.75 0.77 0.1 0.1],'fontname','Arial','fontsize',12);
text(10,0.13,"c",'FontSize',16,'FontWeight','bold','fontname','Arial')

%% d
filepath='/Volumes/Backup Plus/backup/polar/dailySICv2/';
filelist=dir([filepath,'NSIDC0051_SEAICE_PS_S25km_202212*.nc']);
k=length(filelist);
for i=1:k
  filename=[filepath,filelist(i).name];
  SIC=ncread(filename,'F17_ICECON')*100;
  SIC(SIC>100)=nan;
  SIC12(:,:,i)=SIC;
end
C2022=nanmean(SIC12,3);
for i=2:30
    R1(:,:,i)=SIC12(:,:,i+1)-SIC12(:,:,i-1);
end
R12_2022=nanmean(R1,3);

n=1;
for j=1:30
for i=1034:1064
    C(:,:,n)=SIC_daily(:,:,i+1+365*(j-1))-SIC_daily(:,:,i-1+365*(j-1));
    n=n+1;
end
end
R12c=nanmean(C,3);

R1ano=R12_2022-R12c;
R1ano(C2022==0)=0;

%% d
% str=["a","b","c","d"];
% set(gcf,'Position',[350 200 750 700]);
% ha = tight_subplot(2,2,[0.08 0.02],[.06 .04],[.08 .02]);

m_proj('azimuthal equal-area','latitude',-90,'radius',50,'rectbox','on');
hold on;
m_gshhs_i('patch',[.8 .8 .8],'edgecolor',[.8 .8 .8],'MarkerSize',0.1);
m_pcolor(lon,lat,R1ano);
colortable = textread('/Users/jinfei/Documents/MATLAB/colorbar/MPL_bwr.txt');
colormap(flipud(colortable));
m_grid('xaxislocation','top','xtick',[-120 -60 0 60 120 180],'ytick',[-90 -80 -70 -60],'fontsize',11);
caxis([-6 6]);
c=colorbar('FontSize',12);
set(c, 'Position', [0.88 0.15 0.03 0.7]);
title(c,'%/day','FontSize',14)
% text(-0.58,0.56,"d",'FontSize',16,'FontWeight','bold','fontname','Arial')

% print(gcf, '-depsc', '/Users/jinfei/Downloads/1d.eps');


%% polynya contribution feb
clear
clc
load ntSIC
load lon
load lat
Cclim=nanmean(C(:,:,26:12:384),3);

filepath='/Volumes/Backup Plus/backup/polar/dailySICv2_nrt/';
filelist=dir([filepath,'NSIDC0081_SEAICE_PS_S25km_202302*.nc']);
k=length(filelist);
for i=1:k
  filename=[filepath,filelist(i).name];
  SIC=ncread(filename,'F17_ICECON')*100;
  SIC(SIC>100)=nan;
  SIC02(:,:,i)=SIC;
end
C2023=nanmean(SIC02,3);

f1=fopen('/Users/jinfei/Documents/MATLAB/2023/pss25area_v3.dat','rb');
area=fread(f1,[316,332],'int32');


Eclim=sum(area(Cclim>15))/1e9;
E2023=sum(area(C2023>15))/1e9;


ind=find(lon>-180&lon<-120&lat>-80&lat<-70);
Eclimers=sum(area(intersect(find(Cclim>15),ind)))/1e9;
E2023ers=sum(area(intersect(find(C2023>15),ind)))/1e9;


p=(Eclimers-E2023ers)/(Eclim-E2023);






