% DROUGHT ANALYSIS FIGURES

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LOAD RECONSTRUCTION DATA FIRST!!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

addpath('../')


% ALL PROXY TYPES, LMR DATABASE

% ALL ENSEMBLES SAVED
%efl='cesm_lme010_r12000_p998_state_2ozn_avg_AprMar_prxydtst_5_prxtp_tca_2605_swtchbld100_16-Nov-2017_15:57:59.mat';
efl='cesm_lme010_r12000_p998_state_2ozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:24:35.mat';

%pth='/d2/nsteiger/output-da/hydroclimate/'; load([pth,efl])
pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';load([pth,efl])

%varpth=['/d6/haibo/lme/cam5/BLMTRC5CN/ALLFORCEINGS/010'];
%sftlf=permute(ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'landmask'),[2 1]);
%lat=ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'lat');
%lon=ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'lon');



% EXTRACT EACH VARIABLE AND EXPLORE WHETHER OR NOT THERE WAS A MEAN STATE SHIFT IN THAT VARIABLE
% look for shift in ENSO, AMO, local T, 



make_figs='n';



% EXTRACT EACH VARIABLE AND EXPLORE WHETHER OR NOT THERE WAS A MEAN STATE SHIFT IN THAT VARIABLE
% look for shift in ENSO, AMO, local T, 

% What years to look at?
analyrs1=800:1600;
analyrs2=1601:1925;
[~,~,ir1] = intersect(analyrs1,r_o:r_f); % indices of recon segment
[~,~,ir2] = intersect(analyrs2,r_o:r_f); % indices of recon segment
% monthly indices
ovrlprng=rng_intr([r_o,r_f],[analyrs1(1),analyrs1(end)]);
mon_r=r_o:(1/12):(r_f+11/12);
irm1_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
irm1_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;

ovrlprng=rng_intr([r_o,r_f],[analyrs2(1),analyrs2(end)]);
mon_r=r_o:(1/12):(r_f+11/12);
irm2_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
irm2_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;


%========
% AMO
%========

k=find(state_tp=='o');
 
xa_amo_ens=permute(squeeze(Xa_ens(xmeta{k}.id_X(1),:,:)),[2 1]);


xa_v1=xa_amo_ens(ir1,:);
xa_v2=xa_amo_ens(ir2,:);

% Mean offset for AMO is the 1601-1925 period
amo_off=mean(xa_v2(:));


p1=fitdist(xa_v1(:),'kernel');
x1=min(xa_v1(:)):0.01:max(xa_v1(:));
Yc1=pdf(p1,x1);

p2=fitdist(xa_v2(:),'kernel');
x2=min(xa_v2(:)):0.01:max(xa_v2(:));
Yc2=pdf(p2,x2);


figure
plot(x2-amo_off,Yc2,'k','linewidth',3)
hold on
plot(x1-amo_off,Yc1,'linewidth',3,'color',[0 0.4470 0.7410])
%legend(['AMO ' num2str(analyrs2(1)) '-' num2str(analyrs2(end))],['AMO ' num2str(analyrs1(1)) '-' num2str(analyrs1(end))])
set(gca,'fontsize',16)
ylabel('Prob. density','fontsize',18)
xlabel('Temp. anom. (deg C)','fontsize',18)
box on

if  make_figs=='y'
 %print('-dpdf',['./figs/pdf_amo_shift_noleg.pdf'])
 print('-dpng',['./figs/pdf_amo_shift_noleg.png'],'-r300')
end


%===========
% NINO
%===========

k=find(state_tp=='n');
if strcmp(xmeta{k}.indinfo{3}.indxnm,'Nino 3.4'); i=3; else; error('Wrong Nino specification'); end

xa_nino_ens=reshape(permute(Xa_ens((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:,:),[1 3 2]),12*reconYrs,length(p_yrs));
    
xa_v1=xa_nino_ens(irm1_a:irm1_b,:);
xa_v2=xa_nino_ens(irm2_a:irm2_b,:);

% Mean offset for AMO is the 1601-1925 period
nino_off=mean(xa_v2(:));

p1=fitdist(xa_v1(:),'kernel');
x1=min(xa_v1(:)):0.01:max(xa_v1(:));
Yc1=pdf(p1,x1);

p2=fitdist(xa_v2(:),'kernel');
x2=min(xa_v2(:)):0.01:max(xa_v2(:));
Yc2=pdf(p2,x2);


figure
plot(x2-nino_off,Yc2,'k','linewidth',3)
hold on
plot(x1-nino_off,Yc1,'linewidth',3,'color',[0.8500 0.3250 0.0980])
%legend(['Nino3.4 ' num2str(analyrs2(1)) '-' num2str(analyrs2(end))],['Nino3.4 ' num2str(analyrs1(1)) '-' num2str(analyrs1(end))])
set(gca,'fontsize',16)
ylabel('Prob. density','fontsize',18)
xlabel('Temp. anom. (deg C)','fontsize',18)
box on
xlim([-6 5])

if  make_figs=='y'
 %print('-dpdf',['./figs/pdf_nino_shift_noleg.pdf'])
 print('-dpng',['./figs/pdf_nino_shift_noleg.png'],'-r300')
end


%------------------------------------------------------------------
% PLOT TIME SERIES


%========
% AMO
%========

k=find(state_tp=='o');

pltyrs=800:2000;
[~,~,ip] = intersect(pltyrs,r_o:r_f); % indices of recon segment

figure
plot(pltyrs,squeeze(Xa_prct(xmeta{k}.id_X(1),2,ip))-amo_off,'linewidth',2,'color',[0 0.4470 0.7410])
hold on
plot(pltyrs,squeeze(Xa_prct(xmeta{k}.id_X(1),1,ip))-amo_off,'--','color',[0.5 0.5 0.5])
plot(pltyrs,squeeze(Xa_prct(xmeta{k}.id_X(1),3,ip))-amo_off,'--','color',[0.5 0.5 0.5])
legend('AMO (NASST)','5-95th prcntl')
set(gca,'fontsize',16)
xlim([pltyrs(1) pltyrs(end)])
xlabel('Years','fontsize',18)
ylabel('Temp. anom. (deg C)','fontsize',18)
box on


if  make_figs=='y'
 %print('-dpdf',['./figs/timeseries_amo_recon_analyrs_' num2str(pltyrs(1)) '_' num2str(pltyrs(end)) '.pdf'])
 print('-dpng',['./figs/timeseries_amo_recon_analyrs_' num2str(pltyrs(1)) '_' num2str(pltyrs(end)) '.png'],'-r300')
end


%========
% TEMP
%========

k=find(state_tp=='2');

pltyrs=800:2000;
[~,~,ip] = intersect(pltyrs,r_o:r_f); % indices of recon segment

gmtmd=squeeze(Xa_prct(xmeta{k}.id_X(1),2,ip));
gmto=mean(gmtmd); % mean offset

figure
plot(pltyrs,squeeze(Xa_prct(xmeta{k}.id_X(1),2,ip))-gmto,'linewidth',2)
hold on
plot(pltyrs,squeeze(Xa_prct(xmeta{k}.id_X(1),1,ip))-gmto,'--','color',[0.5 0.5 0.5])
plot(pltyrs,squeeze(Xa_prct(xmeta{k}.id_X(1),3,ip))-gmto,'--','color',[0.5 0.5 0.5])
set(gca,'fontsize',16)
xlim([pltyrs(1) pltyrs(end)])
xlabel('Years','fontsize',18)
ylabel('Temp. anom. (deg C)','fontsize',18)
box on
% Plot line showing end of analysis period
yl=ylim;line([1925 1925],yl,'color','k','linewidth',2)
legend('GMT','5-95th prcntl')

make_figs='y';
if  make_figs=='y'
 %print('-dpdf',['./figs/timeseries_gmt_recon_analyrs_' num2str(pltyrs(1)) '_' num2str(pltyrs(end)) '.pdf'])
 print('-dpng',['./figs/timeseries_gmt_recon_analyrs_' num2str(pltyrs(1)) '_' num2str(pltyrs(end)) '.png'],'-r300')
end
make_figs='n';

%==========
% Nino 3.4
%==========

k=find(state_tp=='n');
if strcmp(xmeta{k}.indinfo{3}.indxnm,'Nino 3.4'); i=3; else; error('Wrong Nino specification'); end

ovrlprng=rng_intr([r_o,r_f],[pltyrs(1),pltyrs(end)]);
mon_r=r_o:(1/12):(r_f+11/12);
irm_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
irm_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;


x1=reshape(permute(Xa_prct((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:,:),[2 1 3]),3,12*reconYrs);

figure
%plot(mon_r(irm_a:irm_b),squeeze(x1(2,irm_a:irm_b)),'linewidth',2)
hold on
hp=plot(mon_r(irm_a:irm_b),squeeze(x1(1,irm_a:irm_b))-nino_off,'--','color',[0.5 0.5 0.5]);
plot(mon_r(irm_a:irm_b),squeeze(x1(3,irm_a:irm_b))-nino_off,'--','color',[0.5 0.5 0.5])
hn=plot(mon_r(irm_a:irm_b),squeeze(x1(2,irm_a:irm_b))-nino_off,'linewidth',2,'color',[0.8500 0.3250 0.0980]);
legend([hn,hp],'Nino3.4 (mon)','5-95th prcntl')
set(gca,'fontsize',16)
xlim([pltyrs(1) pltyrs(end)])
xlabel('Years','fontsize',18)
ylabel('Temp. anom. (deg C)','fontsize',18)
box on


if  make_figs=='y'
 %print('-dpdf',['./figs/timeseries_nino_recon_analyrs_' num2str(pltyrs(1)) '_' num2str(pltyrs(end)) '.pdf'])
 print('-dpng',['./figs/timeseries_nino_recon_analyrs_' num2str(pltyrs(1)) '_' num2str(pltyrs(end)) '.png'],'-r300')
end





































