clear
clc


load sstpy2_anns_oisst
data0=sstpy2*1e12;
% load v2_anns_oisst
% data0=v2*1024;
poslon=lon>=200&lon<=260;
poslat=lat>=-2&lat<=2;
data1=squeeze(mean(mean(data0(:,poslat,poslon),3),2));
tim1=tim;
postim=tim>=1993&tim<=2021;
data1m=nanmean(data1(postim)); 
data1=(data1-data1m)*2+30;
tim1=tim1(2:end);
data1=data1(2:end);

% load eke_anns_reof_sst.mat
% data0=eke*1e4;
% poslon=lon>=200&lon<=260;
% poslat=lat>=-2&lat<=2;
% data1=squeeze(mean(mean(data0(:,poslat,poslon),3),2));
% s1=data1(tim>=1982&tim<=2021);
% tim1=tim;
% postim=tim>=1993&tim<=2021;
% data1m=nanmean(data1(postim)); 
% data1=(data1-data1m)*1+2400;

% load eke0_anns_cmems
load eke_anns_model
data0=eke*1024;
poslon=lon>=200&lon<=260;
poslat=lat>=-2&lat<=2;
data2=squeeze(nanmean(nanmean(data0(:,poslat,poslon),3),2));
tim2=tim;
postim=tim>=1993&tim<=2021;
data2m=nanmean(data2(postim));
data2=(data2-data2m)+90;

load eke0_anns_ofes2.mat
% eke=squeeze(nanmean(eke(:,dep>=-3,:,:),2));
data0=eke*1024;
poslon=lon>=200&lon<=260;
poslat=lat>=-2&lat<=2;
data3=squeeze(nanmean(nanmean(data0(:,poslat,poslon),3),2));
s3=data3(tim>=1982&tim<=2021);
tim3=tim;
postim=tim>=1993&tim<=2019;
data3m_em=nanmean(data3(postim));
postim=tim>=1993&tim<=2021;
data3m=nanmean(data3(postim));
data3=(data3-data3m);

load eke_anns_tao
eke=squeeze(nanmean(eke(:,dep>=25&dep<=35,:,:),2));
data0=eke*1024;
poslon=lon>=220&lon<=220;
poslat=lat>=-.2&lat<=.2;
data4=squeeze(nanmean(nanmean(data0(:,poslat,poslon),3),2));
tim4=tim;
% plot(tim4,data4)
postim=tim>=1993&tim<=2020;
data4m=nanmean(data4(postim));
data4=(data4-data4m)/2+60;
s=data4(postim);
% s(isnan(s))=300;
trend(s,1,'omitnan')*10/data4m

% (nanmean(data4(tim4>=1982&tim4<=2006)))/data4m
% (nanmean(data4(tim4>=2007&tim4<=2021)))/data4m

%% plot
tname='Time series of TIW intensity (different datasets)';
ax1=1956;
ax2=2023;
ax3=-25;
ax4=125;

figure(1)
set(gcf,'unit','centimeters','position',[0 0 9 13])
% plot([ax1 ax2],[0 0],':m')
plot([ax1 ax2],[0 0],':','color',[0.00 0.45 0.74])
hold on
plot([ax1 ax2],[30 30],':','color',[0.47 0.67 0.19])
plot([ax1 ax2],[60 60],'k:')
plot([ax1 ax2],[90 90],':','color',[1 0 0])
plot(tim1,movmean(data1,1),'-','color',[0.47 0.67 0.19],'linewidth',.3);
plot(tim2,movmean(data2,1),'-','color',[1 0 0],'linewidth',.3);
plot(tim4,movmean(data4,1),'-k','linewidth',.3);
plot(tim3,movmean(data3,1),'-','color',[0.00 0.45 0.74],'linewidth',.3);

smt=3;
load sshpy2_anns_ofes2_10.mat
data0=sshpy2*1e14;
poslat=lat>=-1&lat<=0.5;
poslon=lon>=200&lon<=260;
data10=squeeze(nanmean(nanmean(data0(:,:,poslat,poslon),3),4));
data20=movmean(data10,smt-2,2);
data10m=squeeze(nanmean(data20,1));
data1mm=nanmean(data10m(tim>=1993&tim<=2019));
tim5=tim;
data1a=1*squeeze(nanstd(data20,0,1))/data1mm*data3m_em;
lat=tim5((smt-1)/2:end-(smt-1)/2+1);
s1a=data1a((smt-1)/2:end-(smt-1)/2+1);
s1am=data3(tim3>=lat(1)&tim3<=lat(end))';
s1=s1am-s1a;
s2=s1am+s1a;
h2=fill([lat,fliplr(lat)],[s1,fliplr(s2)],'r','linestyle','none');
set(h2,'facealpha',0.3,'facecolor',[0.1 0.1 0.1]*5)

aap=0.05;
postim1=tim1>=1982&tim1<=2021;
postim3=tim3>=1982&tim3<=2021;
[rr,p,dof,r05]=corr_free(data1(postim1),data3(postim3),aap);
postim2=tim2>=1993&tim2<=2021;
postim3=tim3>=1993&tim3<=2021;
[rr,p,dof,r05]=corr_free(data2(postim2),data3(postim3),aap);
postim4=tim4>=1983&tim4<=2020;
postim3=tim3>=1983&tim3<=2020;
[rr,p,dof,r05]=corr_free(data4(postim4),data3(postim3),aap);

smt=9;
data1=movmean(data1,smt);
data2=movmean(data2,smt);
data3=movmean(data3,smt);
data40=data4;
data4=movmean(data4,smt,'omitnan');
data4(isnan(data40))=nan;


trend(data1(tim1>=1997&tim1<=2017),1,'omitnan')*10/data1m/2000
trend(data2(tim2>=1997&tim2<=2017),1,'omitnan')*10/data2m
trend(data3(tim3>=1997&tim3<=2017),1,'omitnan')*10/data3m
% trend(data3(tim3>=1997&tim3<=2017),1,'omitnan')*10/data3m


f1=plot(tim1((smt-1)/2:end-(smt-1)/2+1),data1((smt-1)/2:end-(smt-1)/2+1),'-','color',[0.47 0.67 0.19],'linewidth',1.5);
f2=plot(tim2((smt-1)/2:end-(smt-1)/2+1),data2((smt-1)/2:end-(smt-1)/2+1),'-','color',[1 0 0],'linewidth',1.5);
f3=plot(tim3((smt-1)/2:end-(smt-1)/2+1),data3((smt-1)/2:end-(smt-1)/2+1),'-','color',[0.00 0.45 0.74],'linewidth',1.5);
% plot(tim1([1:(smt-1)/2]),data1([1:(smt-1)/2]),':','color',[0.47 0.67 0.19],'linewidth',1.5);
% plot(tim2([1:(smt-1)/2]),data2([1:(smt-1)/2]),':','color',[1 0 0],'linewidth',1.5);
% plot(tim4(tim4<=1986),data4(tim4<=1986),':k','linewidth',1.5);
% plot(tim3([1:(smt-1)/2]),data3([1:(smt-1)/2]),':','color',[0.00 0.45 0.74],'linewidth',1.5);
% plot(tim1([end-(smt-1)/2+1:end]),data1([end-(smt-1)/2+1:end]),':','color',[0.47 0.67 0.19],'linewidth',1.5);
% plot(tim2([end-(smt-1)/2+1:end]),data2([end-(smt-1)/2+1:end]),':','color',[1 0 0],'linewidth',1.5);
% plot(tim4(tim4>=2010),data4(tim4>=2010),':k','linewidth',1.5);
% plot(tim3([end-(smt-1)/2+1:end]),data3([end-(smt-1)/2+1:end]),':','color',[0.00 0.45 0.74],'linewidth',1.5);
f4=plot(tim4(tim4>=1986&tim4<=2010),data4(tim4>=1986&tim4<=2010),'k-','linewidth',1.5);

% [p,s1]=polyfit(tim1(tim1>=1997&tim1<=2017),data1(tim1>=1997&tim1<=2017),1);
% y=p(2)+p(1)*tim1(tim1>=1997&tim1<=2017);
% plot(tim1(tim1>=1997&tim1<=2017),y,'--','color',[0.47 0.67 0.19],'linewidth',1.5)
% [p,s2]=polyfit(tim2(tim2>=1997&tim2<=2017),data2(tim2>=1997&tim2<=2017),1);
% y=p(2)+p(1)*tim2(tim2>=1997&tim2<=2017);
% plot(tim2(tim2>=1997&tim2<=2017),y,'--','color',[1 0 0],'linewidth',1.5)
% data40(tim4==2014)=1600;
% [p,s4]=polyfit(tim4(tim4>=1993&tim4<=2020),data40(tim4>=1993&tim4<=2020),1);
% y=p(2)+p(1)*tim4(tim4>=1993&tim4<=2020);
% plot(tim4(tim4>=1993&tim4<=2020),y,'--k','linewidth',.2)
% [p,s3]=polyfit(tim3(tim3>=1997&tim3<=2017),data3(tim3>=1997&tim3<=2017),1);
% y=p(2)+p(1)*tim3(tim3>=1997&tim3<=2017);
% plot(tim3(tim3>=1997&tim3<=2017),y,'--','color',[0.00 0.45 0.74],'linewidth',1.5)

load sshpy2_anns_ofes2_10.mat
data0=sshpy2*1e14;
poslat=lat>=-1&lat<=0.5;
poslon=lon>=200&lon<=260;
data1=squeeze(nanmean(nanmean(data0(:,:,poslat,poslon),3),4));
data2=movmean(data1,smt,2);
data1m=squeeze(nanmean(data2,1));
data1mm=nanmean(data1m(tim>=1993&tim<=2019));
tim5=tim;
data1a=1*squeeze(nanstd(data2,0,1))/data1mm*data3m_em;
lat=tim5((smt-1)/2:end-(smt-1)/2+1);
s1a=data1a((smt-1)/2:end-(smt-1)/2+1);
s1am=data3(tim3>=lat(1)&tim3<=lat(end))';
s1=s1am-s1a;
s2=s1am+s1a;
h2=fill([lat,fliplr(lat)],[s1,fliplr(s2)],'r','linestyle','none');
set(h2,'facealpha',0.3,'facecolor',[0.1 0.1 0.1]*2)


hold off

% hl=legend([f2 f4 f1 f3],['Altimeter EKE'],'TAO EKE',['SST meander intensity'],['OFES2 EKE']);
% set(hl,'location','west','box','off','orientation','vertical','fontsize',7)
% hm=text(1960,1150,['Time series of TIW intensity',sprintf('\n'),'(1993–2021 mean removed)'],...
%     'rotation',0,'fontsize',7,'HorizontalAlignment','left');

ylabel('J m-3')
xlabel('Year')
set(gca,'tickdir','out')
set(gca,'ytick',[-50:10:ax4])
set(gca,'xtick',[1960:10:2020])
xticklabel=char(get(gca,'XTickLabel'));
set(gca,'xticklabel',{xticklabel})
set(gca,'xgrid','on')
set(gca,'xminortick','on')
set(gca,'ticklength',[0.005 0.005])

axis([ax1 ax2 ax3 ax4])
la=axis;
% title(tname)

set(gca,'fontsize',7)
box on