clear
clc

%% a
load lon
load lat
index{1}=find(lon >= -60 & lon <= 20); %weddell
index{2}=find(lon >= 20 & lon <= 130); %indian ocean
index{3}=find(not(lon >= -160 & lon <= 130)); %wpo
index{4}=find(lon >= -160 & lon <= -130); %eastern ross
index{5}=find(lon >= -130 & lon <= -60); %abs
index{6}=1:length(lon(:));

load SIC_daily
for x=1:6
for i=759:1123
    Cclim=nanmean(SIC_daily(:,:,i:365:11708),3);
    pick=intersect(cell2mat(index(x)),find(Cclim>14&Cclim<16));
    edgeclim(x,i-758)=mean(lat(pick));
end
end

filepath='/Volumes/Backup Plus/backup/polar/dailySICv2/';
filelist=dir([filepath,'NSIDC0051_SEAICE_PS_S25km_2022*.nc']);
k=length(filelist);
for i=60:k
  filename=[filepath,filelist(i).name];
  SIC=ncread(filename,'F17_ICECON')*100;
  SIC(SIC>100)=nan;
  C0(:,:,i-59)=SIC;
end
filepath='/Volumes/Backup Plus/backup/polar/dailySICv2_nrt/';
filelist=dir([filepath,'N*.nc']);
for i=1:59
  filename=[filepath,filelist(i).name];
  SIC=ncread(filename,'F17_ICECON')*100;
  SIC(SIC>100)=nan;
  C0(:,:,i+306)=SIC;
end

for x=1:6
for i=1:365
    Csingle=C0(:,:,i);
    pick=intersect(cell2mat(index(x)),find(Csingle>14&Csingle<16));
    edge(x,i)=mean(lat(pick));
end
end
ea=edge-edgeclim;

for x=1:6
   eamov(x,:)=movmean(ea(x,:),9); 
end



%%
% str=["a","b"];
% set(gcf,'Position',[350 200 700 800]);
% ha = tight_subplot(2,1,[0.08 0],[.1 .08],[.08 .02]);
% axes(ha(1))

set(gcf,'Position',[350 200 700 400]);
color{1}=[232/256 137/256 189/256];
color{2}=[103/256 194/256 163/256];
color{3}=[252/256 140/256 99/256];
color{4}=[142/256 160/256 201/256];
color{5}=[230/256 175/256 12/256];
color{6}=[0.3 0.3 0.3];
plot(zeros(1,92),'--k','LineWidth',1)
hold on
for x=[1 2 3 5 6]
    q(x)=plot(eamov(x,215:306),'LineWidth',1.5,'color',cell2mat(color(x)));
end
q(4)=plot(eamov(4,215:306),'LineWidth',3,'color',cell2mat(color(4)));
ylabel('Sea ice edge anomaly (degree)','fontname','Arial','fontsize',14)
axis([1 92 -8 6])
legend([q(1) q(2) q(3) q(4) q(5) q(6)],{'WS','IO','WPO','ERS','ABS','SO'},'position',[0.8 0.71 0.08 0.15],'fontname','Arial','fontsize',12);
% 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);
set(gca,'xtick',[1 32 62],'xticklabel',{'Oct' 'Nov' 'Dec'},'fontname','Arial','fontsize',14);

% text(10,5.5,str(1),'FontSize',16,'FontWeight','bold','fontname','Arial')

% axes(ha(2))
% m_proj('azimuthal equal-area','latitude',-90,'radius',50,'rectbox','on');
% hold on;
% m_grid('xaxislocation','top','ytick',[],'xtick',[-60,20,130,-160,-130],'fontsize',15,'linest','-');
% m_gshhs_i('patch',[.8 .8 .8],'edgecolor',[.8 .8 .8],'MarkerSize',0.1);
% text(-0.3,0.43,'WS','FontSize',15,'FontWeight','Bold');
% text(0.38,0.3,'IO','FontSize',15,'FontWeight','Bold');
% text(0.02,-0.4,'WPO','FontSize',15,'FontWeight','Bold');
% text(-0.4,-0.4,'ERS','FontSize',15,'FontWeight','Bold');
% text(-0.5,-0.15,'ABS','FontSize',15,'FontWeight','Bold');
% text(-0.55,0.55,str(2),'FontSize',16,'FontWeight','bold','fontname','Arial')
% m_line([-120,-120],[-80,-70],'color','r','linewi',2);
% m_line([-180,-180],[-80,-70],'color','r','linewi',2);
% m_line(-180:-120,-70:-70,'color','r','linewi',2);
% m_line(-180:-120,-80:-80,'color','r','linewi',2);

%% model
clear
clc

f1='/Volumes/Backup Plus/PAROCE/OND_2022_sic.nc';
f2='/Volumes/Backup Plus/PAROCE/OND_clim_sic.nc';
sic22=ncread(f1,'siconc')*100; 
sic22(sic22<5)=nan;
sicclim=ncread(f2,'siconc')*100; 
sicclim(sicclim<5)=nan;

X=ncread(f1,'nav_lon');
Y=ncread(f1,'nav_lat');
X=double(X);
Y=double(Y);

index{1}=find(X >= -60 & X <= 20); %weddell
index{2}=find(X >= 20 & X <= 130); %indian ocean
index{3}=find(not(X >= -160 & X <= 130)); %wpo
index{4}=find(X >= -160 & X <= -130); %eastern ross
index{5}=find(X >= -130 & X <= -60); %abs
index{6}=1:length(X(:));

for x=1:6
    for i=1:92
        cc=sicclim(:,:,i);
        pick=intersect(cell2mat(index(x)),find(cc>14&cc<16));
        edgeclim(x,i)=mean(Y(pick));
    end
end


for x=1:6
    for i=1:92
        c22=sic22(:,:,i);
        pick=intersect(cell2mat(index(x)),find(c22>14&c22<16));
        edge22(x,i)=mean(Y(pick));
    end
end


ea=edge22-edgeclim;

for x=1:6
   eamov(x,:)=movmean(ea(x,:),9); 
end
%%


set(gcf,'Position',[350 200 700 400]);
color{1}=[232/256 137/256 189/256];
color{2}=[103/256 194/256 163/256];
color{3}=[252/256 140/256 99/256];
color{4}=[142/256 160/256 201/256];
color{5}=[230/256 175/256 12/256];
color{6}=[0.3 0.3 0.3];
plot(zeros(1,92),'--k','LineWidth',1)
hold on
for x=[1 2 3 5 6]
    q(x)=plot(eamov(x,:),'LineWidth',1.5,'color',cell2mat(color(x)));
end
q(4)=plot(eamov(4,:),'LineWidth',3,'color',cell2mat(color(4)));
ylabel('Sea ice edge anomaly (degree)','fontname','Arial','fontsize',14)
axis([1 92 -10 6])
legend([q(1) q(2) q(3) q(4) q(5) q(6)],{'WS','IO','WPO','ERS','ABS','SO'},'position',[0.8 0.71 0.08 0.15],'fontname','Arial','fontsize',12);
set(gca,'xtick',[1 32 62],'xticklabel',{'Oct' 'Nov' 'Dec'},'fontname','Arial','fontsize',14);


%% fig1 model
clear
clc

f1='/Volumes/Backup Plus/PAROCE/OND_2022_sic.nc';
f2='/Volumes/Backup Plus/PAROCE/OND_clim_sic.nc';
sic22=ncread(f1,'siconc')*100; 
sic22(sic22<5)=nan;
sicclim=ncread(f2,'siconc')*100; 
sicclim(sicclim<5)=nan;

X=ncread(f1,'nav_lon');
Y=ncread(f1,'nav_lat');
X=double(X);
Y=double(Y);
f3='/Volumes/Backup Plus/PAROCE/so025_eORCA025_cellarea.nc';
ca=ncread(f3,'cell_area');
for i=1:92
    ind1=find(sic22(:,:,i)>15);
    SIE22(i)=sum(ca(ind1))/1e12;
    ind2=find(sicclim(:,:,i)>15);
    SIEc(i)=sum(ca(ind2))/1e12;
end

% d0=smoothdata(diff(SIE22),'gaussian',30);
% dc=smoothdata(diff(SIEc),'gaussian',30);
d0=diff(SIE22);
dc=diff(SIEc);
d0(92)=nan;
dc(92)=nan;

set(gcf,'Position',[500 400 750 450]);
plot(1:92,zeros(1,92),'--','color','k','LineWidth',1)
hold on
ylabel('SIE tendency (million km^2/day)','fontname','Arial','fontsize',14)
qc=plot(dc,'color','k','LineWidth',1.5);
q0=plot(d0,'color','r','LineWidth',2.5);
set(gca,'xtick',[1 32 62],'xticklabel',{'Oct' 'Nov' 'Dec'},'fontname','Arial','fontsize',14);
axis([1 92 -1 0.15])
legend('2022','Clim','fontname','Arial','fontsize',14);

%% d
clear
clc

f1='/Volumes/Backup Plus/PAROCE/OND_2022_sic.nc';
f2='/Volumes/Backup Plus/PAROCE/OND_clim_sic.nc';
sic22=ncread(f1,'siconc')*100; 
sic22(sic22<5)=nan;
sicclim=ncread(f2,'siconc')*100; 
sicclim(sicclim<5)=nan;

X=ncread(f1,'nav_lon');
Y=ncread(f1,'nav_lat');
X=double(X);
Y=double(Y);

for i=62:91
    R1(:,:,i-61)=sic22(:,:,i+1)-sic22(:,:,i-1);
end
R12_2022=nanmean(R1,3);

for i=62:91
    R1c(:,:,i-61)=sicclim(:,:,i+1)-sicclim(:,:,i-1);
end
R12c=nanmean(R1c,3);

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

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(X,Y,R1ano);
colortable = textread('/Users/jinfei/Documents/MATLAB/colorbar/MPL_bwr.txt');
colormap(flipud(colortable));
m_grid('xaxislocation','top','xticklabels',[],'yticklabels',[],'fontsize',8);
caxis([-6 6]);
c=colorbar('FontSize',12);
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');




























