clear all 

% global parameters
FAC=100; % convert m to cm
LW=3;
FS=16;
AR=0.65;
XLIM=[1950 2010];XTICK=[XLIM(1):5:XLIM(2)];
YLIM=[-0.3 0.3];YTICK=[YLIM(1):0.1:YLIM(2)];
YLABEL='0-2000m T anomaly (K)';
%TAG='48N65N65W10E';
TAG='48N65N60W15W';
ZMAX=2000;
STARTYEARS=[1960:1:2000];
%STARTYEARS=[];

% read data
load(sprintf('Data/TSBOX_cmip6_%s_0-%4.4i_test.mat',TAG,ZMAX));
load(sprintf('Data/TSBOX_cmip6_%s_0-%4.4i.mat',TAG,ZMAX));
%load(sprintf('Data/TSBOX_obs_%s_0-%4.4i.mat',TAG,ZMAX));

% compute climatologies 
f5010=mean(boxtf((1950-YEAR1F+1):(2010-YEAR1F)));
a5010=mean(boxta((1950-YEAR1A+1):(2010-YEAR1A)));
o5010=mean(boxto((1950-YEAR1O+1):(2010-YEAR1O)));

% plot 
figure(1)
clf
figure_height_scale=0.65;
set(gcf,'paperunits','inches','papersize',[8 6*figure_height_scale],'inverthardcopy','off','color',[1 1 1],...
        'paperposition',[0 0 8 6*figure_height_scale],'renderer','painters')
hold on 
% free run
X=[YEAR1F:YEARNF fliplr(YEAR1F:YEARNF)];
Y=[bmxtf' fliplr(bmntf')]-f5010
patch(X,Y,[0.1 0.9 1],'edgecolor','none')
hf=plot(YEAR1F:YEARNF,boxtf-f5010,'color',[0 0 1],'linewidth',LW);
% analysis
X=[YEAR1A:YEARNA fliplr(YEAR1A:YEARNA)];
Y=[bmxta' fliplr(bmnta')]-f5010;
patch(X,Y,[1 0.7 0.5],'edgecolor','none')
ha1=plot(YEAR1A:YEARNA,boxta-f5010,'color',[1 0 0],'linewidth',LW);
% obs
ho=plot(YEAR1O:YEARNO,boxto-f5010-o5010+a5010,'color','k','linewidth',LW);
% hindcasts 
if ~isempty(STARTYEARS)
  for n=numel(YEARS):-1:1
    if YEARS(n)>=STARTYEARS(1) & YEARS(n)<=STARTYEARS(end)
      xvec=[YEARS(n)-1 YEARS(n):(YEARS(n)+size(boxth,1)-1)];
      yvec=[boxta(YEARS(n)-YEAR1A) boxth(:,n)']-f5010; 
      hh1=plot(xvec,yvec,'color',[1 0 1],'linewidth',LW); 
    end 
  end
end
      xvec=[YEARSTEST(1)-1 YEARSTEST(1):(YEARSTEST(1)+size(boxtt,1)-1)];
      yvec=[boxta(YEARSTEST(1)-YEAR1A) boxtt(:,n)']-f5010;
      ht1=plot(xvec,yvec,'color',['g'],'linewidth',LW+0.5);
% misc 
plot([1995 1995],YLIM,'--','color',[1 1 1]*0.5)
hold off 
set(gca,'xlim',XLIM,'ylim',YLIM,'xtick',XTICK,'ytick',YTICK,'fontsize',FS,'fontweight','bold','xcolor','k');
ylabel(YLABEL,'fontsize',FS,'fontweight','bold');
legend([ho,hf,ha1,hh1],{'EN4','FREE','ANALYSIS','HINDCAST'},'location','northwest')
if isempty(STARTYEARS)
  print(gcf,'-depsc2',sprintf('SPNAtemp0-%im_cmip6_analysis.eps',ZMAX));
else
  print(gcf,'-depsc2',sprintf('SPNAtemp0-%im_cmip6_hindcast%4.4i-%4.4i_test.eps',ZMAX,STARTYEARS(1),STARTYEARS(end)));
end
