clear
clc

load eke0_anns_ofes2
datestr(tim);
poslat=abs(lat)<=10;
poslon=lon>=150&lon<=285;
timup=1993;
timend=2021;
postim=tim>=timup&tim<=timend;
data0=eke(postim,poslat,poslon)*1e4;
lon=lon(poslon);
lat=lat(poslat); 
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

data3=squeeze(nanstd(data30,1));

% %%
% [x,y]=meshgrid(lon,lat);
% 
% % plot
% figure(1)
% set(gcf,'unit','centimeters','position',[5 5 14 7])
% cax1=-150;
% cax2=150;  
% ax1=150;
% ax2=285;
% ax3=-10;
% ax4=10;
% nl=40;
% data2t=data3;
% data2t=movmean(data2t,1,2,'omitnan');
% data2t(data2t>=cax2)=cax2;
% data2t(data2t<=cax1)=cax1;
% 
% m_proj('Miller','lon',[ax1 ax2],'lat',[ax3 ax4]);
% m_contourf(x,y,data2t,nl,'linestyle','none');
% % m_contourf(x,y,data2t,nl,'linewidth',0.2,'color',[0.1 0.1 0.1]*6);
% m_coast('patch',[.1 .1 .1]*9,'linestyle','none');        
% hold on 
% % m_contour(x,y,data2t,[0. 0.02],'-','color',[0.1 0.1 0.1]*6,'linewidth',0.3);
% m_coast('line','color',[.1 .1 .1]*6);   
% m_plot([ax1 ax2],[0 0],'k:','linewidth',1)
% 
% % m_plot([200 200],[-2 2],'k--','linewidth',.5)
% % m_plot([260 260],[-2 2],'k--','linewidth',.5)
% % m_plot([200 260],[-2 -2],'k--','linewidth',.5)
% % m_plot([200 260],[2 2],'k--','linewidth',.5)
% hold off
% 
% h1=colorbar;
% % set(h1,'ticks',cax1:(cax2-cax1)/10:cax2)
% caxis([cax1 cax2])
% h2=get(h1,'title');
% set(h2,'string','cm^2s^{-2} / decade')
% % load br_sparse
% % colormap(br_sparse)
% colormap(othercolor('BuOr_12'))
% 
% title(['EKE'])
% 
% xlim=get(gca,'xlim');
% xtick=xlim(1):5/(ax2-ax1)*(xlim(end)-xlim(1)):xlim(end);
% ylim=get(gca,'ylim');
% ytick=ylim(1):1/(ax4-ax3)*(ylim(end)-ylim(1)):ylim(end);
% 
% xticklabel={'160E','180W','160W','140W','120W','100W','80W'};
% yticklabel={'10S','5S','EQ','5N','10N'};
% 
% set(gca,'tickdir','out')
% set(gca,'xtick',xtick(3:4:end))
% set(gca,'xticklabel',xticklabel)
% set(gca,'ytick',ytick(1:5:end))
% set(gca,'yticklabel',yticklabel)
% % set(gca,'yminortick','on')
% set(gca,'ticklength',[0.012 0.025])
% 
% set(gca,'fontsize',10)

lat0=lat;

%
s1m=squeeze(nanmean(nanmean(data2(:,lon>=200&lon<=260,:),2),3));
s1ma=squeeze(nanstd(nanmean(data2(:,lon>=200&lon<=260,:),2),1,3));
s1a=squeeze(nanstd(nanmean(data30(:,:,lon>=200&lon<=260),3),1));
s1am=squeeze(nanmean(data300(:,lon>=200&lon<=260),2));
lat1=[-6.7:0.1:9.7];
s1m=interp1(lat0,s1m,lat1,'linear');
s1ma=interp1(lat0,s1ma,lat1,'linear');
s1a=interp1(lat0,s1a,lat1,'linear');
s1am=interp1(lat0,s1am,lat1,'linear');
lat=lat1;
zzz=zeros(size(s1am));
zzz0=zeros(size(s1am));

figure(1)
set(gcf,'unit','centimeters','position',[5 5 3.5 5.01])
ax1=-10;
ax2=30;
ax3=-7;
ax4=10;
max1=nanmax(s1m);
% max1=1.556e-11*1e10;  % average between 1993-2021
s1m=s1m/max1*100;
s1ma=s1ma/max1*100;
s1a=s1a/max1*100;
s1am=s1am/max1*100;
s1a0=s1a;
s1am0=s1am;

% f1=plot(s1m,lat,'color',[0.1 0.1 0.1]*6,'linewidth',.5);
hold on
% s1am(s1am<0)=nan;
% f2=plot(s1am,lat,'color',[1 0 0],'linewidth',1);
% s1am=s1am0;
% s1am(s1am>=0)=nan;
% f2=plot(s1am,lat,'color',[0 1 1],'linewidth',1);
% h1=fill([s1m-s1ma,fliplr(s1m+s1ma)],[lat,fliplr(lat)],'k','linestyle','none');
% set(h1,'facealpha',0.3,'facecolor',[0.1 0.1 0.1]*6)

s1=s1am0-s1a;
s2=s1am0+s1a;
s2(s2<0)=0;
zzz(s1>0)=s1(s1>0);
h2=fill([zzz,fliplr(s2)],[lat,fliplr(lat)],'r','linestyle','none');
set(h2,'facealpha',0.3,'facecolor',[0.1 0.1 0.1]*6)

s1=s1am0-s1a;
s2=s1am0+s1a;
s1(s1>0)=0;
zzz=zzz0;
zzz(s2<0)=s2(s2<0);
h2=fill([s1,fliplr(zzz)],[lat,fliplr(lat)],'b','linestyle','none');
set(h2,'facealpha',0.3,'facecolor',[0.1 0.1 0.1]*6)

% s1am=s1am0;
% s1a=s1a0;
% s1a(s1am>=0)=0;
% s1am(s1am>=0)=0;
% h2=fill([s1am-s1a,fliplr(s1am+s1a)],[lat,fliplr(lat)],'b','linestyle','none');
% set(h2,'facealpha',0.5,'facecolor',[0 0. 1])

plot([ax1 ax2],[0 0],'k:','linewidth',.5)
plot([0 0],[ax3 ax4],'k:','linewidth',.5)
plot([100 100],[ax3 ax4],'k:','linewidth',.5)

s1am=s1am0;
s1am(s1am<0)=nan;
f2=plot(s1am,lat,'color',[1 0 0],'linewidth',1);
s1am=s1am0;
s1am(s1am>=0)=nan;
f2=plot(s1am,lat,'color',[0 1 1],'linewidth',1);
hold off

% lh=legend([f1 f2],'Mean','Trend');
% set(lh,'location','northeast','box','on','fontsize',7,'Color','none','TextColor',[.1 .1 .1]*0)

% title(['<SST''y2>'])

% xtick=-20:0.2:1;
ytick=-5:5:10;

% xticklabel={'-20','0','20','40','60','80','100'};
yticklabel={'5°S','0°','5°N','10°N'};

set(gca,'tickdir','out')
% set(gca,'xtick',xtick)
set(gca,'ytick',ytick)
% set(gca,'xticklabel',xticklabel)
set(gca,'yticklabel',yticklabel)
% set(gca,'yminortick','on')
set(gca,'ticklength',[0.015 0.025])
axis([ax1 ax2 ax3 ax4])
set(gca,'fontsize',7)
box on

nanmax(s1am0)
s1a0(s1am0==nanmax(s1am0))