clear
clc

% load v2_anns_model
load v2_anns_ofes2
datestr(tim);
poslat=lat>=-7&lat<=10;
poslon=lon>=170&lon<=285;
timup=1993;
timend=2021;
postim=tim>=timup&tim<=timend;
data0=v2(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));

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 4.5 5])
ax1=-8;
ax2=38;
ax3=-4;
ax4=8;
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)
% 0.93,0.69,0.13

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.5);
s1am=s1am0;
s1am(s1am>=0)=nan;
f2=plot(s1am,lat,'color',[0 1 1],'linewidth',1.5);

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=-0:10:100;
ytick=-4:4:10;

xticklabel={'0','10','20','30','40','100'};
yticklabel={'4°S','0°','4°N','8°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))