clear
clc

vfx_lon1=-179.75:0.25:-120;
vfx_lat1=-70;

vfx_lon2=-179.75;
vfx_lat2=-78:0.07:-70;

vfx_lon3=-120;
vfx_lat3=-73.6:0.07:-70;

f1='/Volumes/Backup Plus/PAROCE/simotion_all_icemod_1980_2023.nc';
f2='/Volumes/Backup Plus/PAROCE/sivolu_all_icemod_1980_2023.nc';
f3='/Volumes/Backup Plus/PAROCE/siconc_all_icemod_1980_2023.nc';
lon=ncread(f1,'nav_lon');
lat=ncread(f1,'nav_lat');

for i=1:518
    
%     v=ncread(f1,'sivelv',[1 1 i],[1442 559 1],[1 1 1]);
    u=ncread(f1,'sivelu',[1 1 i],[1442 559 1],[1 1 1]);
    t=ncread(f2,'sivolu',[1 1 i],[1442 559 1],[1 1 1]);
    c=ncread(f3,'siconc',[1 1 i],[1442 559 1],[1 1 1]);
    SIT_gate = griddata(double(lon),double(lat),double(t),double(vfx_lon2),double(vfx_lat2));
    SIC_gate = griddata(double(lon),double(lat),double(c),double(vfx_lon2),double(vfx_lat2));
%     V_gate = griddata(double(lon),double(lat),double(v),double(vfx_lon1),double(vfx_lat1));
    U_gate = griddata(double(lon),double(lat),double(u),double(vfx_lon2),double(vfx_lat2));
%     L_gate = distance([-70,-120.25],[-70,-120])*6371000*3.14/180;
    L_gate = distance([-70.07,-120],[-70,-120])*6371000*3.14/180;
    Qx = SIT_gate.*U_gate.*SIC_gate*L_gate*(60*60*24*30);
    Qx_all2(i)= sum(Qx)/1e9;

end

%% transport
clear
clc

load Qx_all1
load Qx_all2
load Qx_all3

for i=1:12
    Qx_clim1(i)=mean(Qx_all1(12+i:12:372));
    Qx_clim2(i)=mean(Qx_all2(12+i:12:372));
    Qx_clim3(i)=mean(Qx_all3(12+i:12:372));
end
Qx_c1=[Qx_clim1(3:12) Qx_clim1(1:2)];
Qx_c2=[Qx_clim2(3:12) Qx_clim2(1:2)];
Qx_c3=[Qx_clim3(3:12) Qx_clim3(1:2)];

% 2022.3-2023.2
Qx_ano1=Qx_all1(507:518)-Qx_c1;
Qx_ano2=Qx_all2(507:518)-Qx_c2;
Qx_ano3=Qx_all3(507:518)-Qx_c3;

Qx_ano=Qx_ano1-Qx_ano2+Qx_ano3;

% area
f2='/Volumes/Backup Plus/PAROCE/so025_eORCA025_cellarea.nc';
ca=ncread(f2,'cell_area');
f='/Volumes/Backup Plus/PAROCE/siconc_all_icemod_1980_2023.nc';
lon=ncread(f,'nav_lon');
lat=ncread(f,'nav_lat');
ind=find(lon>-180&lon<-120&lat>-80&lat<-70);
E=sum(ca(ind))/1e6;





%% wind
f='/Volumes/Backup Plus/backup/polar/ERA5/surface7922.nc';
% ncdisp(f)
for m=27:38
    vclim(:,:,m-26)=mean(squeeze(ncread(f,'v10',[1 1 1 m],[1440 361 1 30],[1 1 1 12])),3);
end
x0=ncread(f,'longitude');
y0=ncread(f,'latitude');
[X,Y]=meshgrid(x0,y0);
X=double(X);Y=double(Y);

f0='/Volumes/Backup Plus/backup/polar/ERA5/v23.nc';
v0(:,:,1)=permute(squeeze(ncread(f,'v10',[1 1 1 519],[1440 361 1 1],[1 1 1 1]))-vclim(:,:,1),[2 1]);
v0(:,:,2)=permute(squeeze(ncread(f,'v10',[1 1 1 520],[1440 361 1 1],[1 1 1 1]))-vclim(:,:,2),[2 1]);
v0(:,:,3)=permute(squeeze(ncread(f,'v10',[1 1 1 521],[1440 361 1 1],[1 1 1 1]))-vclim(:,:,3),[2 1]);
v0(:,:,4)=permute(squeeze(ncread(f,'v10',[1 1 1 522],[1440 361 1 1],[1 1 1 1]))-vclim(:,:,4),[2 1]);
v0(:,:,5)=permute(squeeze(ncread(f,'v10',[1 1 1 523],[1440 361 1 1],[1 1 1 1]))-vclim(:,:,5),[2 1]);
v0(:,:,6)=permute(squeeze(ncread(f,'v10',[1 1 1 524],[1440 361 1 1],[1 1 1 1]))-vclim(:,:,6),[2 1]);
v0(:,:,7)=permute(squeeze(ncread(f,'v10',[1 1 1 525],[1440 361 1 1],[1 1 1 1]))-vclim(:,:,7),[2 1]);
v0(:,:,8)=permute(squeeze(ncread(f,'v10',[1 1 2 526],[1440 361 1 1],[1 1 1 1]))-vclim(:,:,8),[2 1]);
v0(:,:,9)=permute(squeeze(ncread(f,'v10',[1 1 2 527],[1440 361 1 1],[1 1 1 1]))-vclim(:,:,9),[2 1]);
v0(:,:,10)=permute(squeeze(ncread(f,'v10',[1 1 2 528],[1440 361 1 1],[1 1 1 1]))-vclim(:,:,10),[2 1]);
v0(:,:,11)=permute(ncread(f0,'v10',[1 1 1],[1440 361 1],[1 1 1])-vclim(:,:,11),[2 1]);
v0(:,:,12)=permute(ncread(f0,'v10',[1 1 2],[1440 361 1],[1 1 1])-vclim(:,:,12),[2 1]);

ind=find(X>-180&X<-120&Y>-80&Y<-70);
for i=1:12
    va=v0(:,:,i);
    vave(i)=mean(va(ind));
end


%% ASL

f='/Users/jinfei/Documents/polar/2023extreme/amundsen-sea-low-index-master/asli_era5_v3-latest.xlsx';
asli=readtable(f);
lon23=str2double(table2array(asli(784:795,2))); %2022.3-2023.2
act23=str2double(table2array(asli(784:795,4)));

lonc0=str2double(table2array(asli(292:651,2)));
actc0=str2double(table2array(asli(292:651,4)));
for i=1:12
    lonc(i)=mean(lonc0(i:12:360));
    actc(i)=mean(actc0(i:12:360));
end

lona=lon23-lonc';
acta=act23-actc';

%% albedo
f='/Volumes/Backup Plus/PAROCE/albedo_all_icemod_1980_2023.nc';
% ncdisp(f)
alb=ncread(f,'alb_ice');
lon=ncread(f,'nav_lon');
lat=ncread(f,'nav_lat');
for i=20:24
    albclim(:,:,i-19)=nanmean(alb(:,:,i:12:372),3);
end
alb22=alb(:,:,512:516);
albano=alb22-albclim;
ind=find(lon>-180&lon<-120&lat>-80&lat<-70);
for i=1:5
    aba=albano(:,:,i);
    abave(i)=nanmean(aba(ind));
end


%%

set(gcf,'Position',[350 200 650 900]);
ha = tight_subplot(5,1,[0.05 0],[.05 .05],[.13 .03]);

axes(ha(1))
plot(1:5,acta(6:10),'-','color',sixteen2ten('#359E5A')/255,'LineWidth',2)
ylabel({'ASL depth';'(hPa)'},'fontsize',14)
set(gca,'xtick',1:5,'xticklabel',{'Aug' 'Sep' 'Oct' 'Nov' 'Dec'},'fontname','Arial','fontsize',14);
xlim([1 5])
ylim([-10 0])

axes(ha(2))
plot(1:5,lona(6:10),'-','color',sixteen2ten('#FFA726')/255,'LineWidth',2)
ylabel({'ASL position';'(degree)'},'fontsize',14)
set(gca,'xtick',1:5,'xticklabel',{'Aug' 'Sep' 'Oct' 'Nov' 'Dec'},'fontname','Arial','fontsize',14);
xlim([1 5])
hold on
plot(1:5,zeros(1,5),'--k')

axes(ha(3))
plot(1:5,vave(6:10),'-','color',sixteen2ten('#363472')/255,'LineWidth',2)
ylabel({'Meridional wind';'(m/s)'},'fontsize',14)
set(gca,'xtick',1:5,'xticklabel',{'Aug' 'Sep' 'Oct' 'Nov' 'Dec'},'fontname','Arial','fontsize',14);
xlim([1 5])
hold on
plot(1:5,zeros(1,5),'--k')

axes(ha(4))
plot(1:5,Qx_ano(6:10),'-','color',sixteen2ten('#A14170')/255,'LineWidth',2)
ylabel({'SIV export';'(km^3/mon)'},'fontsize',14)
set(gca,'xtick',1:5,'xticklabel',{'Aug' 'Sep' 'Oct' 'Nov' 'Dec'},'fontname','Arial','fontsize',14);
xlim([1 5])
hold on
plot(1:5,zeros(1,5),'--k')

axes(ha(5))
plot(1:5,abave(1:5),'-','color',sixteen2ten('#469DE9')/255,'LineWidth',2)
ylabel('Albedo','fontsize',14)
set(gca,'xtick',1:5,'xticklabel',{'Aug' 'Sep' 'Oct' 'Nov' 'Dec'},'fontname','Arial','fontsize',14);
xlim([1 5])

