% COMPUTE CRPS TO SEE IF DATING IS CORRECT
p0=reshape(Xa_t2m_sig(:,ir),length(lat),length(lon),length(ir));
phsg=reshape(p0(pt,pn,:),length(pn)*length(pt),length(ir));

spt_crps0=zeros(length(ir)-2,1);
spt_crpsm1=zeros(length(ir)-2,1);
spt_crpsp1=zeros(length(ir)-2,1);
k=1;
for i=2:(length(ir)-1)
   spt_crps0(k)=crps_m(ph_nm1(:,i),phsg(:,i),be_nm1(:,i));
   spt_crpsm1(k)=crps_m(ph_nm1(:,i-1),phsg(:,i-1),be_nm1(:,i));
   spt_crpsp1(k)=crps_m(ph_nm1(:,i+1),phsg(:,i+1),be_nm1(:,i));
   k=k+1;
end

figure
hold on
plot(ir(2:end-1),spt_crps0,'-ok','linewidth',2)
plot(ir(2:end-1),spt_crpsm1,'-ob','linewidth',2)
plot(ir(2:end-1),spt_crpsp1,'-or','linewidth',2)
ax=gca;
line('ydata',ax.YLim,'xdata',[1883 1883],...
'linestyle','-','color','green','linewidth',2)
ax.YGrid = 'on';
box on
xlim([1871 1891])
ylim([0.5 3])
xlabel('Year','fontsize',16)
ylabel('Spatial CRPS','fontsize',16)
set(gca,'FontSize',14)
legend('0','-1','+1')
title('CRPS(PHYDA,BEarth) Eurasia')

print(['/home/nsteiger/Dropbox/enso-volcanoes/phyda_bearth_djf_crps_18711891_eurasia.png'],'-dpng','-r300');



