clear
clc

load sst_anns_ersst
poslat1=lat>=-15&lat<=0;
poslon1=lon>=200&lon<=360-80;
poslat2=lat>=0&lat<=15;
poslon2=lon>=200&lon<=360-80;
postim=tim>=1993&tim<=2021;
sst1=squeeze(nanmean(nanmean(sst(:,poslat1,poslon1),2),3));
sst2=squeeze(nanmean(nanmean(sst(:,poslat2,poslon2),2),3));
% sst1=sst1-nanmean(sst1(postim));
% sst2=sst2-nanmean(sst2(postim));
% sst3=sst3-nanmean(sst3(postim));
data=sst2-sst1;
s4m=nanmean(data(postim));
data=movmean(movmean(data,9),1);
data=data-nanmean(data(postim));
s4=data;
tim4=tim;

load sst_anns_oisst_globe
poslat1=lat>=-15&lat<=0;
poslon1=lon>=200&lon<=360-80;
poslat2=lat>=0&lat<=15;
poslon2=lon>=200&lon<=360-80;
postim=tim>=1993&tim<=2021;
sst1=squeeze(nanmean(nanmean(sst(:,poslat1,poslon1),2),3));
sst2=squeeze(nanmean(nanmean(sst(:,poslat2,poslon2),2),3));
% sst1=sst1-nanmean(sst1(postim));
% sst2=sst2-nanmean(sst2(postim));
% sst3=sst3-nanmean(sst3(postim));
data=sst2-sst1;
s2m=nanmean(data(postim));
data=movmean(movmean(data,9),1);
data=data-nanmean(data(postim));
s2=data;
tim2=tim;

load sst_anns_ofes2_globe.mat
poslat1=lat>=-15&lat<=0;
poslon1=lon>=200&lon<=360-80;
poslat2=lat>=0&lat<=15;
poslon2=lon>=200&lon<=360-80;
postim=tim>=1993&tim<=2021;
sst1=squeeze(nanmean(nanmean(sst(:,poslat1,poslon1),2),3));
sst2=squeeze(nanmean(nanmean(sst(:,poslat2,poslon2),2),3));
% sst1=sst1-nanmean(sst1(postim));
% sst2=sst2-nanmean(sst2(postim));
% sst3=sst3-nanmean(sst3(postim));
data=sst2-sst1;
s3m=nanmean(data(postim));
data=movmean(movmean(data,9),1);
data=data-nanmean(data(postim));
s3=data;
tim3=tim;

load sst_anns_hadisst
poslat1=lat>=-15&lat<=0;
poslon1=lon>=200&lon<=360-80;
poslat2=lat>=0&lat<=15;
poslon2=lon>=200&lon<=360-80;
postim=tim>=1993&tim<=2021;
sst1=squeeze(nanmean(nanmean(sst(:,poslat1,poslon1),2),3));
sst2=squeeze(nanmean(nanmean(sst(:,poslat2,poslon2),2),3));
% sst1=sst1-nanmean(sst1(postim));
% sst2=sst2-nanmean(sst2(postim));
% sst3=sst3-nanmean(sst3(postim));
data=sst2-sst1;
s1m=nanmean(data(postim));
data=movmean(movmean(data,9),1);
data=data-nanmean(data(postim));
s1=data;
tim1=tim;

% save ssty_time_series s1 s2 s3 s4 tim1 tim2 tim3 tim4 s1m s2m s3m s4m


%%
% plot
tname='SSTy';
ax1=1957;
ax2=2022;
ax3=-.35;
ax4=.55;

figure(1)
set(gcf,'unit','centimeters','position',[5 5 6 7])
set(gca,'position',[0.15 0.11 0.775 0.815])

plot([ax1 ax2],[0 0],'k:')
hold on
% plot([1960 1960],[ax3 ax4],'k:','linewidth',1)
smt=9;

% s1(tim1<1960)=nan;
% s2(tim2<1960)=nan;
% s3(tim3<1960)=nan;
% s4(tim4<1960)=nan;

% f1=plot(tim1((smt-1)/2:end-(smt-1)/2+1),s1((smt-1)/2:end-(smt-1)/2+1),'-','color',[0.93 0.69 0.13],'linewidth',1);
f4=plot(tim4((smt-1)/2:end-(smt-1)/2+1),s4((smt-1)/2:end-(smt-1)/2+1),'-g','linewidth',1);
f3=plot(tim3((smt-1)/2:end-(smt-1)/2+1),s3((smt-1)/2:end-(smt-1)/2+1),'-','color',[0.00 0. 1],'linewidth',1);
f2=plot(tim2((smt-1)/2:end-(smt-1)/2+1),s2((smt-1)/2:end-(smt-1)/2+1),'-','color',[1 0 0],'linewidth',1);

% % plot(tim1([1:(smt-1)/2]),s1([1:(smt-1)/2]),':','color',[0.93 0.69 0.13],'linewidth',1);
% plot(tim4([1:(smt-1)/2]),s4([1:(smt-1)/2]),':g','linewidth',1);
% plot(tim3([1:(smt-1)/2]),s3([1:(smt-1)/2]),':','color',[0.00 0. 1],'linewidth',1);
% plot(tim2([1:(smt-1)/2]),s2([1:(smt-1)/2]),':','color',[1 0 0],'linewidth',1);
% % plot(tim1([end-(smt-1)/2+1:end]),s1([end-(smt-1)/2+1:end]),':','color',[0.93 0.69 0.13],'linewidth',1);
% plot(tim4([end-(smt-1)/2+1:end]),s4([end-(smt-1)/2+1:end]),'g:','linewidth',1);
% plot(tim3([end-(smt-1)/2+1:end]),s3([end-(smt-1)/2+1:end]),':','color',[0.00 0. 1],'linewidth',1);
% plot(tim2([end-(smt-1)/2+1:end]),s2([end-(smt-1)/2+1:end]),':','color',[1 0 0],'linewidth',1);

% load ssty_time_series
% s1(tim1>1960)=nan;
% s2(tim2>1960)=nan;
% s3(tim3>1960)=nan;
% s4(tim4>1960)=nan;
% tim1=(tim1-1960)/3+1960;
% tim2=(tim2-1960)/3+1960;
% tim3=(tim3-1960)/3+1960;
% tim4=(tim4-1960)/3+1960;
% f4=plot(tim4((smt-1)/2:end-(smt-1)/2+1),s4((smt-1)/2:end-(smt-1)/2+1),'-','color',[0.93 0.69 0.13],'linewidth',1);
% f1=plot(tim1((smt-1)/2:end-(smt-1)/2+1),s1((smt-1)/2:end-(smt-1)/2+1),'-','color',[0.47 0.67 0.19],'linewidth',1);
% f2=plot(tim3((smt-1)/2:end-(smt-1)/2+1),s3((smt-1)/2:end-(smt-1)/2+1),'-','color',[0.00 0.45 0.74],'linewidth',1);
% f3=plot(tim2((smt-1)/2:end-(smt-1)/2+1),s2((smt-1)/2:end-(smt-1)/2+1),'-','color',[1 0 0],'linewidth',1);
% plot(tim4([1:(smt-1)/2]),s4([1:(smt-1)/2]),':','color',[0.93 0.69 0.13],'linewidth',1);
% plot(tim1([1:(smt-1)/2]),s1([1:(smt-1)/2]),':','color',[0.47 0.67 0.19],'linewidth',1);
% plot(tim3([1:(smt-1)/2]),s3([1:(smt-1)/2]),':','color',[0.00 0.45 0.74],'linewidth',1);
% plot(tim2([1:(smt-1)/2]),s2([1:(smt-1)/2]),':','color',[1 0 0],'linewidth',1);
% plot(tim4([end-(smt-1)/2+1:end]),s4([end-(smt-1)/2+1:end]),':','color',[0.93 0.69 0.13],'linewidth',1);
% plot(tim1([end-(smt-1)/2+1:end]),s1([end-(smt-1)/2+1:end]),':','color',[0.47 0.67 0.19],'linewidth',1);
% plot(tim3([end-(smt-1)/2+1:end]),s3([end-(smt-1)/2+1:end]),':','color',[0.00 0.45 0.74],'linewidth',1);
% plot(tim2([end-(smt-1)/2+1:end]),s2([end-(smt-1)/2+1:end]),':','color',[1 0 0],'linewidth',1);


hold off

% hl=legend([f4 f1 f2 f3],'ERSSTv5','HadISST','OFES2','Satellite');
hl=legend([f2 f4 f3],'Satellite','ERSSTv5','OFES2');
set(hl,'location','northwest','box','off','orientation','vertical','fontsize',6)


% xlabel('Year') 
ylabel('°C')

set(gca,'tickdir','out')
set(gca,'xtick',[1840:10:2020])
% set(gca,'xtick',[((1870:30:1930)-1960)/3+1960,1960:10:2020])
% xticklabel=char(get(gca,'XTickLabel'));
% set(gca,'xticklabel',{'1870','1900','1930','1960','1970','1980','1990','2000','2010','2020'})
set(gca,'xgrid','on')
% set(gca,'xminortick','on')
set(gca,'ticklength',[0.008 0.025])

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

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