clear
clc

load eke_anns_ofes2
datestr(tim);
poslat=lat>=-2&lat<=2;
poslon=lon>=200&lon<=260;
timup=1993;
timend=2021;
postim=tim>=timup&tim<=timend;
data0=squeeze(nanmean(eke(postim,:,poslat,poslon),4))*1024;
lon=lat(poslat);
lat=dep;
tim=tim(postim);
nt=length(tim);
ny=length(lat);
nx=length(lon);
data1=nan(ny,nx,nt);
run=4;
data2=nan(ny,nx,nt-2*run);
for t=1:nt
    data1(:,:,t)=data0(t,:,:);
end
for j=1:ny 
for i=1:nx
    s=squeeze(data1(j,i,:));
    s1=movmean(s,2*run+1);
    s1(isnan(s))=nan;
    data2(j,i,:)=s1(1+run:end-run);
end
end
data300=trend(data2,1)*10;
tim=tim(1+run:end-run);
for t=1:timend-run-(timup+run)-9+1
postim=tim>=timup+run-1+t&tim<=timup+run-1+t+9;
data30(t,:,:)=trend(data2(:,:,postim),1)*10;
t
end
eke=squeeze(nanmean(data300(lat>=-50,abs(lon)<=2),'all'));
data=squeeze(nanstd(data30,1,1));
ekea=squeeze(nanmean(data(lat>=-50,abs(lon)<=2),'all'));

load uvuy_anns_ofes2
datestr(tim);
poslat=lat>=-2&lat<=2;
poslon=lon>=200&lon<=260;
timup=1993;
timend=2021;
postim=tim>=timup&tim<=timend;
data0=-squeeze(nanmean(uvuy(postim,:,poslat,poslon),4))*1024;
lon=lat(poslat);
lat=dep;
tim=tim(postim);
nt=length(tim);
ny=length(lat);
nx=length(lon);
data1=nan(ny,nx,nt);
run=4;
data2=nan(ny,nx,nt-2*run);
for t=1:nt
    data1(:,:,t)=data0(t,:,:);
end
for j=1:ny 
for i=1:nx
    s=squeeze(data1(j,i,:));
    s1=movmean(s,2*run+1);
    s1(isnan(s))=nan;
    data2(j,i,:)=s1(1+run:end-run);
end
end
data300=trend(data2,1)*10;
tim=tim(1+run:end-run);
for t=1:timend-run-(timup+run)-9+1
postim=tim>=timup+run-1+t&tim<=timup+run-1+t+9;
data30(t,:,:)=trend(data2(:,:,postim),1)*10;
t
end
btr=squeeze(nanmean(data300(lat>=-50,abs(lon)<=2),'all'));
data=squeeze(nanstd(data30,1,1));
btra=squeeze(nanmean(data(lat>=-50,abs(lon)<=2),'all'));

load rouw_anns_ofes2
datestr(tim);
poslat=lat>=-2&lat<=2;
poslon=lon>=200&lon<=260;
timup=1993;
timend=2021;
postim=tim>=timup&tim<=timend;
data0=-squeeze(nanmean(rouw(postim,:,poslat,poslon),4))*9.8;
lon=lat(poslat);
lat=dep;
tim=tim(postim);
nt=length(tim);
ny=length(lat);
nx=length(lon);
data1=nan(ny,nx,nt);
run=4;
data2=nan(ny,nx,nt-2*run);
for t=1:nt
    data1(:,:,t)=data0(t,:,:);
end
for j=1:ny 
for i=1:nx
    s=squeeze(data1(j,i,:));
    s1=movmean(s,2*run+1);
    s1(isnan(s))=nan;
    data2(j,i,:)=s1(1+run:end-run);
end
end
data300=trend(data2,1)*10;
tim=tim(1+run:end-run);
clear data30
for t=1:timend-run-(timup+run)-9+1
postim=tim>=timup+run-1+t&tim<=timup+run-1+t+9;
data30(t,:,:)=trend(data2(:,:,postim),1)*10;
t
end
bcr=squeeze(nanmean(data300(lat>=-50,abs(lon)<=2),'all'));
data=squeeze(nanstd(data30,1,1));
bcra=squeeze(nanmean(data(lat>=-50,abs(lon)<=2),'all'));

load uvuy_anns_ofes2
datestr(tim);
poslat=lat>=-2&lat<=2;
poslon=lon>=200&lon<=260;
timup=1993;
timend=2021;
postim=tim>=timup&tim<=timend;
data0=-squeeze(nanmean(uvuy(postim,:,poslat,poslon),4))*1024;
% data0=data0(:,1:end-1,:)+0.5*diff(data0,1,2);
load rouw_anns_ofes2
datestr(tim);
poslat=lat>=-2&lat<=2;
poslon=lon>=200&lon<=260;
timup=1993;
timend=2021;
postim=tim>=timup&tim<=timend;
data00=-squeeze(nanmean(rouw(postim,:,poslat,poslon),4))*9.8;
data0=data0+data00;
lon=lat(poslat);
lat=dep;
tim=tim(postim);
nt=length(tim);
ny=length(lat);
nx=length(lon);
data1=nan(ny,nx,nt);
run=4;
data2=nan(ny,nx,nt-2*run);
for t=1:nt
    data1(:,:,t)=data0(t,:,:);
end
for j=1:ny 
for i=1:nx
    s=squeeze(data1(j,i,:));
    s1=movmean(s,2*run+1);
    s1(isnan(s))=nan;
    data2(j,i,:)=s1(1+run:end-run);
end
end
data300=trend(data2,1)*10;
tim=tim(1+run:end-run);
clear data30
for t=1:timend-run-(timup+run)-9+1
postim=tim>=timup+run-1+t&tim<=timup+run-1+t+9;
data30(t,:,:)=trend(data2(:,:,postim),1)*10;
t
end
tcr=squeeze(nanmean(data300(lat>=-50,abs(lon)<=2),'all'));
data=squeeze(nanstd(data30,1,1));
tcra=squeeze(nanmean(data(lat>=-50,abs(lon)<=2),'all'));

load eket_anns_ofes2
datestr(tim);
poslat=lat>=-2&lat<=2;
poslon=lon>=200&lon<=260;
timup=1993;
timend=2021;
postim=tim>=timup&tim<=timend;
data0t=squeeze(nanmean(eket(postim,:,poslat,poslon),4))*1024;
% data0t=data0t(:,1:end-1,:)+0.5*diff(data0t,1,2);
load uvuy_anns_ofes2
datestr(tim);
poslat=lat>=-2&lat<=2;
poslon=lon>=200&lon<=260;
timup=1993;
timend=2021;
postim=tim>=timup&tim<=timend;
data0=-squeeze(nanmean(uvuy(postim,:,poslat,poslon),4))*1024;
% data0=data0(:,1:end-1,:)+0.5*diff(data0,1,2);
load rouw_anns_ofes2
datestr(tim);
poslat=lat>=-2&lat<=2;
poslon=lon>=200&lon<=260;
timup=1993;
timend=2021;
postim=tim>=timup&tim<=timend;
data00=-squeeze(nanmean(rouw(postim,:,poslat,poslon),4))*9.8;
data0=data0t-(data0+data00);
lon=lat(poslat);
lat=dep;
tim=tim(postim);
nt=length(tim);
ny=length(lat);
nx=length(lon);
data1=nan(ny,nx,nt);
run=4;
data2=nan(ny,nx,nt-2*run);
for t=1:nt
    data1(:,:,t)=data0(t,:,:);
end
for j=1:ny 
for i=1:nx
    s=squeeze(data1(j,i,:));
    s1=movmean(s,2*run+1);
    s1(isnan(s))=nan;
    data2(j,i,:)=s1(1+run:end-run);
end
end
data300=trend(data2,1)*10;
tim=tim(1+run:end-run);
clear data30
for t=1:timend-run-(timup+run)-9+1
postim=tim>=timup+run-1+t&tim<=timup+run-1+t+9;
data30(t,:,:)=trend(data2(:,:,postim),1)*10;
t
end
sink=squeeze(nanmean(data300(lat>=-50,abs(lon)<=2),'all'));
data=squeeze(nanstd(data30,1,1));
sinka=squeeze(nanmean(data(lat>=-50,abs(lon)<=2),'all'));

load eket_anns_ofes2
datestr(tim);
poslat=lat>=-2&lat<=2;
poslon=lon>=200&lon<=260;
timup=1993;
timend=2021;
postim=tim>=timup&tim<=timend;
data0=squeeze(nanmean(eket(postim,:,poslat,poslon),4))*1024;
lon=lat(poslat);
lat=dep;
tim=tim(postim);
nt=length(tim);
ny=length(lat);
nx=length(lon);
data1=nan(ny,nx,nt);
run=4;
data2=nan(ny,nx,nt-2*run);
for t=1:nt
    data1(:,:,t)=data0(t,:,:);
end
for j=1:ny 
for i=1:nx
    s=squeeze(data1(j,i,:));
    s1=movmean(s,2*run+1);
    s1(isnan(s))=nan;
    data2(j,i,:)=s1(1+run:end-run);
end
end
data300=trend(data2,1)*10;
tim=tim(1+run:end-run);
clear data30
for t=1:timend-run-(timup+run)-9+1
postim=tim>=timup+run-1+t&tim<=timup+run-1+t+9;
data30(t,:,:)=trend(data2(:,:,postim),1)*10;
t
end
eket=squeeze(nanmean(data300(lat>=-50,abs(lon)<=2),'all'));
data=squeeze(nanstd(data30,1,1));
eketa=squeeze(nanmean(data(lat>=-50,abs(lon)<=2),'all'));

%% plot
figure(1)
set(gcf,'unit','centimeters','position',[5 5 9 5])
set(gca,'position',[0.25 0.35 0.73 0.6])
ax1=0.;
ax2=7;
ax3=-32;
ax4=30;


xtick=[1 3:6];
xticklabel={'EKE','EKE_t','BTR','BCR','Sink'};

% tt=7.9*24*3600;
load ekem_tcrm

data0=[eke/ekem 0 eket/tcrm btr/tcrm bcr/tcrm sink/tcrm]*100;
err=[ekea/ekem 0 eketa/tcrm btra/tcrm bcra/tcrm sinka/tcrm]*100;
h=bar(data0(1),'r','barwidth',0.6);
hold on
plot([2 2],[ax3 ax4],':k')

data01=data0;
data01(:)=0;
data01(3)=data0(3);
bar(data01,'facecolor',[0.93,0.69,0.13],'barwidth',0.6)

data01=data0;
data01(:)=0;
data01(4)=data0(4);
bar(data01,'facecolor',[0.93,0.69,0.13],'barwidth',0.6)

data01=data0;
data01(:)=0;
data01(5)=data0(5);
bar(data01,'facecolor',[0.93,0.69,0.13],'barwidth',0.6)

data01=data0;
data01(:)=0;
data01(6)=data0(6);
bar(data01,'facecolor',[0 1 1],'barwidth',0.6)

for n=3:6
    data1=data0;
    err1=err;
    data1(:)=nan;
    data1(n)=data0(n);
    err1(:)=nan;
    err1(n)=err(n);
errorbar(data1,err1,'k');
end
for n=1
    data1=data0;
    err1=err;
    data1(:)=nan;
    data1(n)=data0(n);
    err1(:)=nan;
    err1(n)=err(n);
errorbar(data1,err1,'k');
end

for n=1:1
text(n-0.4,data0(n)+err(n)+4,[num2str(round(data0(n))),' ± ',num2str(round(err(n))),'%'],'fontsize',7)
end
for n=3:5
text(n-0.4,data0(n)+err(n)+4,[num2str(round(data0(n))),' ± ',num2str(round(err(n))),'%'],'fontsize',7)
end
n=6;
text(n-0.4,data0(n)-err(n)-2.5,[num2str(round(data0(n))),' ± ',num2str(round(err(n))),'%'],'fontsize',7)


hold off
axis([ax1 ax2 ax3 ax4])
set(gca,'tickdir','out')
set(gca,'xtick',xtick,'xticklabel',xticklabel)

% ylabel('x 10-4 cm2 s-3')
% la=axis;
% text(la(1),la(4)+0.08*(la(4)-la(3)),'b. Mean value and seasonality of BTR and BCR')

% set(gca,'tickdir','out')
% set(gca,'xtick',xtick)
% set(gca,'xticklabel',xticklabel)
% set(gca,'ytick',ytick)
% set(gca,'yticklabel',yticklabel)
% % set(gca,'yminortick','on')
set(gca,'ticklength',[0.01 0.025])

set(gca,'fontsize',7)

axis([ax1 ax2 ax3 ax4])