%% DJF 81-23
clear
clc
f1='/Volumes/Backup Plus/PAROCE/siconc_all_icemod_1980_2023.nc';% 1980.1-2023.2
sic=ncread(f1,'siconc')*100;
X=ncread(f1,'nav_lon');
X=double(X);
index=find(X >= -180 & X <= -120);
sic12=sic(:,:,12:12:516);
sic01=sic(:,:,13:12:517);
sic02=sic(:,:,14:12:518);
f2='/Volumes/Backup Plus/PAROCE/so025_eORCA025_cellarea.nc';
ca=ncread(f2,'cell_area');
for i=1:43
    E01=sum(ca(intersect(index,find(sic01(:,:,i)>15))))/1e12;
    E02=sum(ca(intersect(index,find(sic02(:,:,i)>15))))/1e12;
    E03=sum(ca(intersect(index,find(sic12(:,:,i)>15))))/1e12;
    E0(i)=(E01+E02+E03)/3;
%     E0(i)=E03;
end
modEano=E0-mean(E0(1:30));
mod=normalize(modEano);

load ntSIC
filepath='/Volumes/Backup Plus/backup/polar/dailySICv2_nrt/';
filelist=dir([filepath,'N*.nc']);
k=length(filelist);
for i=1:59
  filename=[filepath,filelist(i).name];
  SIC=ncread(filename,'F17_ICECON')*100;
  SIC(SIC>100)=nan;
  C0(:,:,i)=SIC;
end
C(:,:,529)=nanmean(C0(:,:,1:31),3);
C(:,:,530)=nanmean(C0(:,:,32:59),3);
id=fopen('pss25area_v3.dat','rb');
area=fread(id,[316,332],'int32');
load lon
index=find(lon >= -180 & lon <= -120);
for i=1:530
    E(i)=sum(area(intersect(index,find(C(:,:,i)>15))))/1e9;
end
for i=1:43
    E1(i)=nanmean(E(12+i*12:14+i*12));
%     E1(i)=E(12+i*12);
end
E1(8)=nan;
obsEano=E1-nanmean(E1(1:30));
obs=normalize(obsEano);


plot(1:43,ones(43),'--','color','k','LineWidth',1)
hold on
plot(1:43,-1*ones(43),'--','color','k','LineWidth',1)
p1=plot(1:43,obs,'color',sixteen2ten('#E85472')/255,'LineWidth',2);
p2=plot(1:43,mod,'color',sixteen2ten('#2D59D2')/255,'LineWidth',2);
ylabel('Standardized SIE anomaly','fontname','Arial','fontsize',16)
set(gca,'xtick',1:10:41,'xticklabel',{'1981' '1991' '2001' '2011' '2021'},'fontname','Arial','fontsize',14);
axis([1 43 -4 4])
legend([p1 p2],{'Observation','Model'},'fontsize',14)
text(25,3.5,'R=0.51','fontname','Arial','fontsize',16)

mod(8)=[];
obs(8)=[];
% allstats(mod,obs)
[r,p] = corr(mod',obs');
%% SON 81-22
clear
clc
f1='/Volumes/Backup Plus/PAROCE/siconc_all_icemod_1980_2023.nc';% 1980.1-2023.2
sic=ncread(f1,'siconc')*100;
X=ncread(f1,'nav_lon');
X=double(X);
index=find(X >= -180 & X <= -120);
sic12=sic(:,:,18:12:518);
sic01=sic(:,:,19:12:518);
sic02=sic(:,:,20:12:518);
f2='/Volumes/Backup Plus/PAROCE/so025_eORCA025_cellarea.nc';
ca=ncread(f2,'cell_area');
for i=1:42
    E01=sum(ca(intersect(index,find(sic01(:,:,i)>15))))/1e12;
    E02=sum(ca(intersect(index,find(sic02(:,:,i)>15))))/1e12;
    E03=sum(ca(intersect(index,find(sic12(:,:,i)>15))))/1e12;
    E0(i)=(E01+E02+E03)/3;
end
modEano=E0-mean(E0(1:30));
mod=normalize(modEano);

load ntSIC
id=fopen('pss25area_v3.dat','rb');
area=fread(id,[316,332],'int32');
load lon
index=find(lon >= -180 & lon <= -120);
for i=1:528
    E(i)=sum(area(intersect(index,find(C(:,:,i)>15))))/1e9;
end
for i=1:42
    E1(i)=nanmean(E(18+i*12:20+i*12));
end
obsEano=E1-nanmean(E1(1:30));
obs=normalize(obsEano);

plot(1:42,ones(42),'--','color','k','LineWidth',1)
hold on
plot(1:42,-1*ones(42),'--','color','k','LineWidth',1)
p1=plot(1:42,obs,'color',sixteen2ten('#E85472')/255,'LineWidth',2);
p2=plot(1:42,mod,'color',sixteen2ten('#2D59D2')/255,'LineWidth',2);
ylabel('Standardized SIE anomaly','fontname','Arial','fontsize',16)
set(gca,'xtick',1:10:41,'xticklabel',{'1981' '1991' '2001' '2011' '2021'},'fontname','Arial','fontsize',14);
axis([1 42 -4 4])
legend([p1 p2],{'Observation','Model'},'fontsize',14)
text(25,3.5,'R=0.95','fontname','Arial','fontsize',16)

allstats(mod,obs)

%% corrmap
clear
clc

% D/JF
f1='/Volumes/Backup Plus/PAROCE/siconc_all_icemod_1980_2023.nc';% 1980.1-2023.2
sic=ncread(f1,'siconc')*100;
X=ncread(f1,'nav_lon');
Y=ncread(f1,'nav_lat');
X=double(X);
Y=double(Y);
% sic10=sic(:,:,25:12:518);
% sic11=sic(:,:,26:12:518);
% C0=(sic10+sic11)/2;
C0=sic(:,:,26:12:518);
modCano=C0-repmat(nanmean(C0(:,:,1:30),3),1,1,42);

load ntSIC
load lon
load lat
filepath='/Volumes/Backup Plus/backup/polar/dailySICv2_nrt/';
filelist=dir([filepath,'N*.nc']);
k=length(filelist);
for i=1:59
  filename=[filepath,filelist(i).name];
  SIC=ncread(filename,'F17_ICECON')*100;
  SIC(SIC>100)=nan;
  C2(:,:,i)=SIC;
end
C(:,:,529)=nanmean(C2(:,:,1:31),3);
C(:,:,530)=nanmean(C2(:,:,32:59),3);
% C10=C(:,:,37:12:530);
% C11=C(:,:,38:12:530);
% C1=(C10+C11)/2;
C1=C(:,:,38:12:530);
obsCano=C1-repmat(nanmean(C1(:,:,1:30),3),1,1,42);

for m=1:42
    modCanogrid(:,:,m)=griddata(X,Y,modCano(:,:,m),lon,lat);
end

for i=1:316
    for j=1:332
    for m=1:42
        mod(m)=squeeze(modCanogrid(i,j,m));
        obs(m)=squeeze(obsCano(i,j,m));
    end
    valid_indices = ~isnan(mod) & ~isnan(obs);
    mod_valid = mod(valid_indices);
    obs_valid = obs(valid_indices);
    if isempty(mod_valid)
        r(i,j)=nan;
        p(i,j)=nan;
    else
    [r(i,j),p(i,j)] = corr(mod_valid',obs_valid');
    end
    end
end

% h_int_sum(314:316,332)=0;
sig=find(p<0.05);
% sig=find(h_int_sum==1);

r0=r;
r0(p>0.05)=nan;

m_proj('azimuthal equal-area','latitude',-90,'radius',55,'rectbox','on');
hold on;
m_pcolor(lon,lat,r0);
% m_plot(lon(sig(1:20:length(sig))),lat(sig(1:20:length(sig))),'.','color',[0.3 0.3 0.3],'MarkerSize',5);
colortable = textread('/Volumes/Backup Plus/backup/MATLAB/colorbar/cmp_flux.txt');
colormap(colortable);
m_coast('patch',[.7 .7 .7],'edgecolor',[.7 .7 .7],'MarkerSize',0.1);
m_grid('xaxislocation','top','xticklabels',[],'yticklabels',[],'fontsize',8);
caxis([-1 1]);
colorbar



