%% budget(10-12)AS
clear
clc
f='/Volumes/Backup Plus/PAROCE/vfxice_1d_10_12_2022.nc';

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

ind=find(X>-180&X<-120&Y>-80&Y<-70);
% ind=1:(1442*559);
for i=1:92
    vfx22=ncread(f,'vfxice',[1 1 i],[1442 559 1],[1 1 1]);
    flux22(i)=nanmean(vfx22(ind));
end

f='/Volumes/Backup Plus/PAROCE/budget_OND.nc';
% ncdisp(f)
for i=1:92
    vfxbog=ncread(f,'vfxbog',[1 1 i],[1442 559 1],[1 1 1]);
    bog22(i)=nanmean(vfxbog(ind));
    vfxdyn=ncread(f,'vfxdyn',[1 1 i],[1442 559 1],[1 1 1]);
    dyn22(i)=nanmean(vfxdyn(ind));
    vfxopw=ncread(f,'vfxopw',[1 1 i],[1442 559 1],[1 1 1]);
    opw22(i)=nanmean(vfxopw(ind));
    vfxsni=ncread(f,'vfxsni',[1 1 i],[1442 559 1],[1 1 1]);
    sni22(i)=nanmean(vfxsni(ind));
    vfxsur=ncread(f,'vfxsum',[1 1 i],[1442 559 1],[1 1 1]);
    sur22(i)=nanmean(vfxsur(ind));
    vfxbom=ncread(f,'vfxbom',[1 1 i],[1442 559 1],[1 1 1]);
    bom22(i)=nanmean(vfxbom(ind));
%     vfxthin=ncread(f,'vfxthin',[1 1 i],[1442 559 1],[1 1 1]);
%     thin(i)=nanmean(vfxthin(ind));
    vfxres=ncread(f,'vfxres',[1 1 i],[1442 559 1],[1 1 1]);
    res22(i)=nanmean(vfxres(ind));
end


% f='/Volumes/Backup Plus/PAROCE/budget_OND_add.nc';
% ncdisp(f)
% for i=1:92
%     vfxsnw=ncread(f,'vfxsnw',[1 1 i],[1442 559 1],[1 1 1]);
%     snw(i)=nanmean(vfxsnw(ind));
%     vfxsub=ncread(f,'vfxsub',[1 1 i],[1442 559 1],[1 1 1]);
%     sub(i)=nanmean(vfxsub(ind));
%     vfxspr=ncread(f,'vfxspr',[1 1 i],[1442 559 1],[1 1 1]);
%     spr(i)=nanmean(vfxspr(ind));
%     vfxpnd=ncread(f,'vfxpnd',[1 1 i],[1442 559 1],[1 1 1]);
%     pnd(i)=nanmean(vfxpnd(ind));
% end

total22=bog22+dyn22+opw22+sni22+sur22+bom22+res22;
lam22=flux22-total22;

% dyn
f3='/Volumes/Backup Plus/PAROCE/OND_2022_sit.nc';
sit22=ncread(f3,'sivolu');
for i=1:91
    dsit=sit22(:,:,i+1)-sit22(:,:,i);
    vol22(i)=nanmean(dsit(ind));
end
vol22(92)=nan;

dyn22=-vol22-flux22;

dyng22=dyn22;
dyng22(dyn22>0)=0;
dynm22=dyn22;
dynm22(dyn22<0)=0;

%quantification
surq=-sum(sur22(48:71))/-sum(flux22(48:71));
lamq=-sum(lam22(48:71))/-sum(flux22(48:71));
bomq=-sum(bom22(48:71))/-sum(flux22(48:71));

surq=-sum(sur22(72:92))/-sum(flux22(72:92));
lamq=-sum(lam22(72:92))/-sum(flux22(72:92));
bomq=-sum(bom22(72:92))/-sum(flux22(72:92));

%% patch

p1=[sur22',lam22',bom22',dynm22'];
p2=[bog22',opw22',sni22',dyng22'];
a=area(-p1);
hold on
b=area(-p2);

a(1).FaceColor = sixteen2ten('#78CFF6')/255;
a(1).EdgeColor = sixteen2ten('#78CFF6')/255;
a(1).LineWidth = 2;
a(1).FaceAlpha = 0.4;
a(2).FaceColor = sixteen2ten('#7CACF7')/255;
a(2).EdgeColor = sixteen2ten('#7CACF7')/255;
a(2).LineWidth = 2;
a(2).FaceAlpha = 0.4;
a(3).FaceColor = sixteen2ten('#B79FFD')/255;
a(3).EdgeColor = sixteen2ten('#B79FFD')/255;
a(3).LineWidth = 2;
a(3).FaceAlpha = 0.4;
a(4).FaceColor = sixteen2ten('#808AF8')/255;
a(4).EdgeColor = sixteen2ten('#808AF8')/255;
a(4).LineWidth = 2;
a(4).FaceAlpha = 0.4;

b(1).FaceColor = sixteen2ten('#FF6D00')/255;
b(1).EdgeColor = sixteen2ten('#FF6D00')/255;
b(1).LineWidth = 2;
b(1).FaceAlpha = 0.4;
b(2).FaceColor = sixteen2ten('#FFAB40')/255;
b(2).EdgeColor = sixteen2ten('#FFAB40')/255;
b(2).LineWidth = 2;
b(2).FaceAlpha = 0.4;
b(3).FaceColor = sixteen2ten('#C2185B')/255;
b(3).EdgeColor = sixteen2ten('#C2185B')/255;
b(3).LineWidth = 2;
b(3).FaceAlpha = 0.4;
b(4).FaceColor = sixteen2ten('#F5ACC5')/255;
b(4).EdgeColor = sixteen2ten('#F5ACC5')/255;
b(4).LineWidth = 2;
b(4).FaceAlpha = 0.4;

plot(vol22,'-','color','k','LineWidth',2)

set(gca,'xtick',[1 32 62],'xticklabel',{'Oct' 'Nov' 'Dec'},'fontname','Arial','fontsize',12);
axis([1 92 -0.03 0.02])
ylabel('SIT tendency (m/day)','fontname','Arial','fontsize',14)
legend('Surface melt','Lateral melt','Basal melt','Dynamic melt','Basal growth','Lateral growth','Snowice','Dynamic growth','Total','Location','southwest')

% title('2022.10-12 (120-180W,70-80S)')


%% budget(10-12)AS clim
clear
clc
f='/Volumes/Backup Plus/PAROCE/vfxiceclim.nc';
X=ncread(f,'nav_lon');
Y=ncread(f,'nav_lat');
X=double(X);
Y=double(Y);

ind=find(X>-180&X<-120&Y>-80&Y<-70);
for i=1:92
    vfxc=ncread(f,'vfxice',[1 1 i],[1442 559 1],[1 1 1]);
    fluxc(i)=nanmean(vfxc(ind));
end

f1='/Volumes/Backup Plus/PAROCE/vfxbogclim.nc';
f2='/Volumes/Backup Plus/PAROCE/vfxopwclim.nc';
f3='/Volumes/Backup Plus/PAROCE/vfxsniclim.nc';
f4='/Volumes/Backup Plus/PAROCE/vfxsumclim.nc';
f5='/Volumes/Backup Plus/PAROCE/vfxbomclim.nc';
f6='/Volumes/Backup Plus/PAROCE/vfxdynclim.nc';
for i=1:92
    vfxbog=ncread(f1,'vfxbog',[1 1 i],[1442 559 1],[1 1 1]);
    bogc(i)=nanmean(vfxbog(ind));
    vfxopw=ncread(f2,'vfxopw',[1 1 i],[1442 559 1],[1 1 1]);
    opwc(i)=nanmean(vfxopw(ind));
    vfxsni=ncread(f3,'vfxsni',[1 1 i],[1442 559 1],[1 1 1]);
    snic(i)=nanmean(vfxsni(ind));
    vfxsur=ncread(f4,'vfxsum',[1 1 i],[1442 559 1],[1 1 1]);
    surc(i)=nanmean(vfxsur(ind));
    vfxbom=ncread(f5,'vfxbom',[1 1 i],[1442 559 1],[1 1 1]);
    bomc(i)=nanmean(vfxbom(ind));
    vfxdyn=ncread(f6,'vfxdyn',[1 1 i],[1442 559 1],[1 1 1]);
    dync(i)=nanmean(vfxdyn(ind));
end

totalc=bogc+opwc+snic+surc+bomc+dync;
lamc=fluxc-totalc;

% dyn
f4='/Volumes/Backup Plus/PAROCE/OND_clim_sit.nc';
sitc=ncread(f4,'sivolu');
for i=1:91
    dsit=sitc(:,:,i+1)-sitc(:,:,i);
    volc(i)=nanmean(dsit(ind));
end
volc(92)=nan;

dync=-volc-fluxc;

dyngc=dync;
dyngc(dync>0)=0;
dynmc=dync;
dynmc(dync<0)=0;


%quantification
surq=-sum(surc(48:71))/-sum(fluxc(48:71));
lamq=-sum(lamc(48:71))/-sum(fluxc(48:71));
bomq=-sum(bomc(48:71))/-sum(fluxc(48:71));

surq=-sum(surc(72:92))/-sum(fluxc(72:92));
lamq=-sum(lamc(72:92))/-sum(fluxc(72:92));
bomq=-sum(bomc(72:92))/-sum(fluxc(72:92));


sum(volc(1:15))
sum(volc(16:31))
sum(volc(32:46))
sum(volc(47:61))
sum(volc(62:76))
nansum(volc(77:92))
