clear
clc

load rouw_anns_ofes2
datestr(tim);
poslat=lat>=-7&lat<=10;
poslon=lon>=200&lon<=260;
postim=tim>=1993&tim<=2021;
data0=-squeeze(nanmean(rouw(postim,:,poslat,poslon),4))*9.8*1e6;
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
data3=trend(data2,1)*10;
data4=mann_kendall(data2,0.05);

[x,y]=meshgrid(lon,lat);

% plot
figure(1)
set(gcf,'unit','centimeters','position',[5 5 7 5])
set(gca,'position',[0.1 0.35 0.86 0.6])
cax1=-15;
cax2=15;  
ax1=-2.5;
ax2=6;
ax3=-150;
ax4=0;
nl=40;
data2t=data3;
data2t=movmean(data2t,1,2,'omitnan');
data2t(data2t>=cax2)=cax2;
data2t(data2t<=cax1)=cax1;
% data2t(isnan(data2t))=0;

contourf(x,y,data2t,nl,'linestyle','none');
hold on 
plot([0 0],[ax3 ax4],'k:','linewidth',.5)
stipple(x,y,data4,'markersize',1,'density',150,'color',[0.1 0.1 0.1]*5)

level=5:10:150;
datam=squeeze(nanmean(data2,3));
[c,h]=contour(x,y,datam,level,'-','color',[0.1 0.1 0.1]*5,'linewidth',0.5);
clabel(c,h,'labelspacing',720,'color',[0.1 0.1 0.1]*3,'fontsize',6)
[c,h]=contour(x,y,datam,-level,'--','color',[0.1 0.1 0.1]*5,'linewidth',0.5);
clabel(c,h,'labelspacing',720,'color',[0.1 0.1 0.1]*3,'fontsize',6)
[c,h]=contour(x,y,datam,-[0 0],':','color',[0.1 0.1 0.1]*5,'linewidth',1);
clabel(c,h,'labelspacing',720,'color',[0.1 0.1 0.1]*3,'fontsize',6)

hold off

h1=colorbar("southoutside");
set(h1,'position',[0.305 0.2 0.3 0.03])
ticks=[cax1:cax2:cax2];
set(h1,'ticks',ticks,'fontsize',6,'TickLength',0.03,'tickdirection','out')
caxis([cax1 cax2])
% h2=get(h1,'title');
% set(h2,'string','(^oC/100km)^2 / decade','rotation',90)
hold on
hm=text(4.2,-188,['BCR trend',sprintf('\n'),'x 10^{-6} W m^{-3} / decade'],...
    'rotation',0,'fontsize',6,'HorizontalAlignment','center');
% hm=m_text(178,9,['Satellite'],...
%     'rotation',0,'fontsize',6,'HorizontalAlignment','center');
% hm=m_text(220.,-1,['TAO'],...
%     'rotation',0,'fontsize',6,'HorizontalAlignment','center');
hold off
% load br_sparse
% colormap(br_sparse)
colormap(othercolor('BuOr_12'))
hc=cbarrow;
uistack(hc,'top')

xtick=-2:2:6;
ytick=-150:50:0;

yticklabel={'150','100','50','0'};
xticklabel={'2°S','0°','2°N','4°N','6°N'};

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.006 0.025])

set(gca,'fontsize',6)

axis([ax1 ax2 ax3 ax4])