

addpath('../')



% Drought atlases use JJA instead of annual
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:49:45.mat';
%pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';load([pth,efl])
%pth='/d2/nsteiger/output-da/hydroclimate/'; load([pth,efl])

% PDSI ONLY, ENSEMBLE RECONSTRUCTION
%efl='cesm_lme010_r12000_p998_state_q_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_06-Aug-2018_17:32:52.mat';
%efl='cesm_lme010_r12000_p998_state_q_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_06-Aug-2018_17:38:44.mat';
%pth='./output-recon/';load([pth,efl])
%efl='cesm_lme010_r12000_p998_state_q_avg_DecFeb_prxydtst_5_prxtp_tca_1228_swtchbld100_noNAproxies_04-Jun-2019_12:02:57.mat';
%efl='cesm_lme010_r12000_p998_state_q_avg_DecFeb_prxydtst_5_prxtp_tca_2119_swtchbld100_noNASWproxies_04-Jun-2019_12:43:27.mat';
efl='cesm_lme010_r12000_p998_state_q_avg_DecFeb_prxydtst_5_prxtp_tca_75_swtchbld100_chileonly_31-Jan-2020_12:12:57.mat';
pth='/d2/nsteiger/output-da/hydroclimate/nhsh_recontests/';load([pth,efl])


% What years to look at?
%analyrs1=1000:1925;
analyrs1=800:1925;
%analyrs2=1601:1925;
%analyrs3=1000:2000;
analyrs3=800:2000;
%analyrs4=800:1925;
[~,~,ir1] = intersect(analyrs1,r_o:r_f); % indices of recon segment
%[~,~,ir2] = intersect(analyrs2,r_o:r_f); % indices of recon segment
[~,~,ir3] = intersect(analyrs3,r_o:r_f); % indices of recon segment
%[~,~,ir4] = intersect(analyrs4,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;

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


%======================
%======  PHYDA  =======
%======================

k=find(state_tp=='q');
lat=xmeta{k}.lat;
lon=xmeta{k}.lon;

% PULL OUT SPATIAL MEAN OVER AMERICAN SOUTHWEST
%load('/d1/nsteiger/climate-data/cesm-lme/Nrth_A_Swthwst_mask.mat')
%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');

% Load mask for NASW
%load('NASW_mask.mat','nasw');
%lndmsk=nasw;

% Load mask for Chile/Argentina
%msk=load('ChlArg_mask.mat');
%lndmsk=msk.sftlf;
%lndmsk(lndmsk==0)=NaN;

% Land mask over Pampas (+Andes)
%msk=load('Pampas_mask.mat');
%lndmsk=msk.sftlf;
%lndmsk(lndmsk==0)=NaN;

% Land mask over Argentina
%msk=load('argentina_mask.mat');
%lndmsk=msk.argntn;

% Land mask over Chile 
msk=load('chile_mask.mat');
lndmsk=msk.chl;

% Land mask over Chile+Argentina
%msk=load('chile+argentina_mask.mat');
%lndmsk=msk.chlarg;

% Latitude weighting
A=cosd(repmat(xmeta{k}.lat,[1 length(xmeta{k}.lon)]));

% PDSI ensemble
Xa_di_ens=zeros(reconYrs,sub_ens);
for i=1:sub_ens
   Xvar=squeeze(Xa_ens(xmeta{k}.id_X(1):xmeta{k}.id_X(2),i,:));
   Xpdsi=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),reconYrs);
   Xpdsi(xmeta{k}.lndidx,:)=Xvar;
   Xpdsi=reshape(Xpdsi,length(xmeta{k}.lat),length(xmeta{k}.lon),reconYrs);
   Xv_msk=bsxfun(@times,Xpdsi,lndmsk);
   Xa_di_ens(:,i)=wmean_a(Xv_msk,A); % RECONSTRUCTION DROUGHT INDEX
end

% Get the ensemble mean; identical to the original reconstruction
Xvar=squeeze(Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:));
Xpdsi=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),reconYrs);
Xpdsi(xmeta{k}.lndidx,:)=Xvar;
Xpdsi=reshape(Xpdsi,length(xmeta{k}.lat),length(xmeta{k}.lon),reconYrs);
Xv_msk=bsxfun(@times,Xpdsi,lndmsk);
Xa_di=wmean_a(Xv_msk,A); % RECONSTRUCTION DROUGHT INDEX

% How to compute the megadroughts? 1 => 11 year running mean, 2 => 2S2E
%ds1=mdrought(Xa_di(ir1),1,'y'); % 'y' indicates that the mean is removed
%ds2=mdrought(Xa_di(ir1),2,'y');

% How to compute the megadroughts? 1 => 11 year running mean, 2 => 2S2E
%ds4=mdrought(Xa_di(ir4),1,'y'); % 'y' indicates that the mean is removed
% Most severe droughts have severity < -10 and all occur before 1600 (last one ends in 1590)
%si=find(ds4(:,3)<-10);
%si=find(ds4(:,2)+analyrs1(1)<analyrs1(end));
%svmd=ds4(si,:);


%--------------------------------------
% Load Paico-based PDSI reconstruction
%--------------------------------------

fln='chile_phyda_r.800_1995_prxyN.47_03-Oct-2019_16:36:59.mat';% all PDSI-like proxies in Chile box
pth='/home/nsteiger/MATLAB/PaiCo/phyda/results/';
load([pth,fln])

[~,~,ip1] = intersect(analyrs1,result.times); % indices of recon segment
phydatm=mean(Xa_di_ens(ir1,:),2);
phydae=Xa_di_ens(ir1,:);

%paicote0=result.resample.signals(:,ip1)';
paicote0=result.resample.signals(:,:)';
%paicote=bsxfun(@minus,paicote0,phydatm); % same mean for comparison

%paicotm=mean(result.resample.signals(:,ip1));
paicotm=mean(result.resample.signals(:,:));
paicotm=paicotm-mean(paicotm)+mean(phydatm);
paicote=paicote0-mean(paicotm)+mean(phydatm); % same time mean ensemble mean for comparison


% Ensemble estimate of drought periods
sub_ens=1000; % CHANGE FOR PAICO
dsens=struct();
for i=1:sub_ens
%   dsens.ds{i}=mdrought(Xa_di_ens(ir1,i),1,'y');
   dsens.ds{i}=mdrought(paicote(ip1,i),1,'y');
end


% Find which years there is agreement
mds=zeros(length(ir1),1);
for i=1:sub_ens
   for j=1:length(dsens.ds{i})
      m0=zeros(length(ir1),1);
      m0(dsens.ds{i}(j,1):dsens.ds{i}(j,2))=1;
      mds=mds+m0; % add up droughts
   end
end





% figure
% plot(analyrs1,mds,'o')
% hold on
% plot(analyrs1,90*ones(length(analyrs1),1),'r','linewidth',2)
% %title('Ensemble agreement about drought year')
% ylabel('Percent agreement','fontsize',18)
% xlabel('Year','fontsize',18)
% xlim([analyrs1(1) analyrs1(end)])
% box on
% set(gca,'fontsize',16)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]

%print(['~/Dropbox/sa-na-coupled-mega/ens_drght_agreement_chile_paico.png'],'-dpng','-r300');


% Find years of 90% agreement
dyrs=find(mds>=sub_ens*0.9);
%dyrs=find(mds>=sub_ens*0.5);


% NOW FIGURE OUT WHICH YEARS ARE CONTINUOUS CHUNKS
drghts=struct([]); % initialize structure
negpos=ones(length(ir1),1);
negpos(dyrs)=0;
% Find the timing of droughts
i=1;j=1;flag=0;
while i+1 <= length(ir1)
% Find the start of a drought
if negpos(i)==0
   d_start=i;
   i=i+1;
   % Keep going to the end of a drought
   while negpos(i)~=1
       i=i+1;
       if i+1 > length(ir1)
	   d_end=i-1;
	   drghts(j).strtend=[d_start d_end];
	   flag=1;
	   break
       end
   end
   if flag==1;break;end % get out of nested loops
   
   d_end=i-1;
   drghts(j).strtend=[d_start d_end];
   
   j=j+1;
end

i=i+1;
end
   
% Pull out years from structure
d=zeros(length(drghts),2);
for kk=1:length(drghts)
  d(kk,1:2)=drghts(kk).strtend;
end

% Throw out droughts that aren't at least 10 years long
rmv=[];
ii=1;
for jj=1:length(drghts)
   if (d(jj,2)-d(jj,1)) < 10
      rmv(ii)=jj;
      ii=ii+1;
   end
end
d(rmv,:)=[];  

% FOR RIGHT NOW, JUST USE ALL DROUGHTS; NEED TO SORT OUT LATER...
svmd=d;
ens_drght=mds;
analyrs=analyrs1;

%save('svmd_ens_50pagree.mat','svmd','Xa_di','Xa_di_ens')
%save('svmd_ens_90pagree_sa.mat','svmd','Xa_di','Xa_di_ens')

%save('nasw_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%save('pampas_megaD_ens_noNASWproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')

%save('argentina_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%save('chile_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%save('chlarg_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')

%save('chile_megaD_ens_s800.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%save('chile_megaD_ens_s800_noNAproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%save('chile_megaD_ens_s800_noNASWproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%save('chile_megaD_ens_s800_chileonly.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')

% for PaiCo
Xa_di=paicotm';
Xa_di_ens=paicote;
pr_o=result.times(1);
pr_f=result.times(end);
save('chile_megaD_ens_s800_PaiCo_chile.mat','pr_o','pr_f','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')

return % end here without the figures...


%=============
%   FIGURES
%=============



load distinct_colors_20.mat

% Analysis options
plot_p='y'; % plot proxies?
prntfg='n'; % Save figures?




% PHYDA AND HIGHLIGHT MEGADROUGHTS
figure
hold on
%plot(analyrs3,DA_di(id3),'linewidth',1,'color',clrs.grey)
plot(analyrs3,Xa_di(ir3),'linewidth',1,'color',clrs.grey)
%plot(analyrs3,smooth(DA_di(id3),11),'linewidth',2,'color',clrs.red)
plot(analyrs3,smooth(Xa_di(ir3),11),'linewidth',2,'color',clrs.purple)
plot(analyrs3,zeros(length(analyrs3),1)+mean(Xa_di(ir1)),'--','linewidth',2,'color',clrs.yellow)
set(gca,'fontsize',16)
%ylim([-0.8 1])
xlim([analyrs3(1) analyrs3(end)])
xlabel('Year','fontsize',18)
ylabel(['PDSI'],'fontsize',18)
box on
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]

%y=Xa_di(ir3); % using 'l' from above
%yMin=min(y)-0.5;yMax=max(y)+0.5;
ax=gca;yMin=ax.YLim(1);yMax=ax.YLim(2);
ylim([yMin yMax])
% Add a patches for top five times from truth
for i=1:length(svmd)
    p=patch(analyrs3(1)+[d(i,1) d(i,2) d(i,2) d(i,1)],[yMin yMin yMax yMax],clrs.blue);
    p.FaceAlpha=0.35;p.EdgeColor='none';
end

legend({'PHYDA','PHYDA-11'},'location','southeast')
hold off

%print(['figs_sa/drght_pampas_phyda_11mean.png'],'-dpng','-r300');


% PHYDA + NADA AND HIGHLIGHT MEGADROUGHTS
figure
hold on
plot(analyrs3,prctile(Xa_di_ens(ir3,:)',5),'linewidth',1,'color',clrs.grey)
h3=plot(analyrs3,prctile(Xa_di_ens(ir3,:)',95),'linewidth',1,'color',clrs.grey);
h2=plot(analyrs3,Xa_di(ir3),'linewidth',1.5,'color',clrs.purple);
%plot(analyrs3,smooth(DA_di(id3),11),'linewidth',2,'color',clrs.red)
%plot(analyrs3,smooth(Xa_di(ir3),11),'linewidth',2,'color',clrs.purple)
plot(analyrs3,zeros(length(analyrs3),1)+mean(Xa_di(ir1)),'--','linewidth',2,'color',clrs.yellow)
%legend({'NADA','PHYDA'},'location','southeast')
%legend({'PHYDA','NADA'},'location','southeast')
legend([h2,h3],{'PHYDA','5^{th}-95^{th} pctl'},'location','southeast')
set(gca,'fontsize',16)
%ylim([-0.8 1])
xlim([analyrs3(1) analyrs3(end)])
xlabel('Year','fontsize',18)
ylabel(['PDSI'],'fontsize',18)
box on
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]

% y=DA_di(id2); % using 'l' from above
% yMin=min(y)-0.5;yMax=max(y)+0.5;
%ylim([yMin yMax])
% % Add a patches for top five times from truth
% for i=1:9
%     p=patch(analyrs3(1)+[ds4(i,1) ds4(i,2) ds4(i,2) ds4(i,1)],[yMin yMin yMax yMax],clrs.blue);
%     p.FaceAlpha=0.40;p.EdgeColor='none';
% end
% 
% hold off

%print(['figs/drght_cmp_annual.png'],'-dpng','-r300');
%print(['figs_sa/drght_pampas_phyda_ann_unc.png'],'-dpng','-r300');




