%
first_year=1900;last_year=2005;
first_year1=701;last_year1=1200;
first_year3=1501;last_year3=2000;
nx=360;ny=180;

sfirst_year=sprintf('%4.4d',first_year);
slast_year=sprintf('%4.4d',last_year);

fff=15;

% 1) HadISST
load('Data/ninoindexes_HadISST.mat');

first_ind=(first_year-1870)*12+1;
last_ind=(last_year-1870+1)*12;

xi=nino34(first_ind:last_ind);
x=detrend(xi);
x=x-reshape(mean(reshape(x,12,[]),2)*ones(1,last_year-first_year+1),1,[]);

% 2) N20TRAERCN
load('Data/ninoindexes_NAER1850CNOC_f19_g16_06_0701-1200.mat');
xi1=nino34;
x1=detrend(xi1);
x1=x1-reshape(mean(reshape(x1,12,[]),2)*ones(1,last_year1-first_year1+1),1,[]);

% 3) NBF1850OC
load('Data/ninoindexes_NBF1850OC_f19_tn11_02_1501-2000.mat');
xi3=nino34;
x3=detrend(xi3);
x3=x3-reshape(mean(reshape(x3,12,[]),2)*ones(1,last_year3-first_year3+1),1,[]);




%
figure;
orient tall;
nn=3,
lw=0.2;
% 1) NINO3.4 SST (monthly evolution)
subplot(nn,1,1);
xi_mon=mean(reshape(xi,12,[]),2);
xi1_mon=mean(reshape(xi1,12,[]),2);
xi3_mon=mean(reshape(xi3,12,[]),2);

pp=plot(1:12,[xi_mon xi1_mon xi3_mon]);
set(pp(1),'color',0.7*ones(3,1),'linewidth',2);
set(pp(2),'color','g');set(pp(3),'color','r');
%legend(pp,{'HadISST' 'historical (NorESM1-M)' 'historical (BCCR fast)' 'PI'});

set(gca,'xtick',1:12);xlim([1 12]);
set(gca,'fontsize',fff);
ylabel('SST (^{\circ}C)','fontsize',fff);
text(2,26,'a)','fontsize',fff)

% 2) standard deviation of interannual NINO3.4 anomalies
subplot(nn,1,2);
xj=std(reshape(x,12,[]),0,2);
xj1=std(reshape(x1,12,[]),0,2);
xj3=std(reshape(x3,12,[]),0,2);

pp=plot(1:12,[xj xj1 xj3]);
set(pp(1),'color',0.6*ones(3,1),'linewidth',2);
set(pp(2),'color','g');set(pp(3),'color','r');

set(gca,'xtick',1:12);xlim([1 12]);
set(gca,'fontsize',fff);
ylabel('\sigma (^{\circ}C)','fontsize',fff);
text(2,0.35,'b)','fontsize',fff)


% 3) skewness of interannual NINO3.4 anomalies
subplot(nn,1,3);
xk=skewness(reshape(x,12,[]),1,2);
xk1=skewness(reshape(x1,12,[]),1,2);
xk3=skewness(reshape(x3,12,[]),1,2);

pp=plot(1:12,[xk xk1 xk3]);
set(pp(1),'color',0.6*ones(3,1),'linewidth',2);
set(pp(2),'color','g');set(pp(3),'color','r');

set(gca,'xtick',1:12);xlim([1 12]);
set(gca,'fontsize',fff);
xlabel('months','fontsize',fff);ylabel('skewness','fontsize',fff);
text(2,-0.6,'c)','fontsize',fff)

legend(pp,{'HadISST' 'PI (NorESM1-M)' 'PI (NorESM1-F)'},'location','southeast');

if 0
% 4) NINO3.4 power spectrum
subplot(nn,1,4);
%[Pxx,F]=pmtm(x/std(x),[],[],1);
[Pxx1,F]=pmtm(x1/std(x1),[],[],1);
[Pxx3,F]=pmtm(x3/std(x3),[],[],1);
%PXX=[Pxx Pxx1 Pxx2 Pxx3];
PXX=[Pxx1 Pxx3];

pp=loglog((F*12),PXX,'-','LineWidth',1);
%set(pp(1),'color',0.6*ones(3,1),'linewidth',2);
set(pp(1),'color','g');set(pp(2),'color','r');

grid on;
xtick=[60 30 15 8 4 2 1 0.5 0.25];
set(gca,'xtick',1./xtick,'xticklabel',{'60','30','15','8','4','2','1','0.5','0.25'},'xminortick','off','xminorgrid','off')
%xtick=[60 30 15 8 6 4 3 2 1 0.5 0.25];
%set(gca,'xtick',1./xtick,'xticklabel',{'60','30','15','8','6','4','3','2','1','0.5','0.25'},'xminortick','off','xminorgrid','off')


set(gca,'ytick',[0.01 0.1 1 10 100],'yminortick','off','yminorgrid','off')
xlim([1/60 1/0.25]);
ylim([0.01 100]);


%h=legend(ph,{'HadISST' 'historical' 'PI control'},'location','Best');set(h,'FontSize',12)
%legend boxoff
%box on
set(gca,'LineWidth',2,'FontSize',12)
xlabel('Period (years)','FontSize',12)
ylabel('Power spectral density','FontSize',12)
text(1/0.33,50,'d)','fontsize',12)
end

%return;
%print('-depsc','Figs/plot_nino34index_various_2exp.eps');
%print('-dpng','-r300','Figs/plot_nino34index_various_2exp.png');
print('-depsc','Figs/plot_nino34index_various_2exp_3panels.eps');
print('-dpng','-r300','Figs/plot_nino34index_various_2exp_3panels.png');
