% DROUGHT ANALYSIS FIGURES

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

addpath('../')


% % ALL PROXY TYPES, LMR DATABASE
% %efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_AprMar_prxydtst_5_prxtp_tca_2605_swtchbld100_27-Sep-2017_15:41:21.mat';
% %efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2601_swtchbld100_27-Sep-2017_15:43:37.mat';
% %efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2590_swtchbld100_27-Sep-2017_15:45:13.mat';
% % ALL ENSEMBLES SAVED
% efl='cesm_lme010_r12000_p998_state_2ozn_avg_AprMar_prxydtst_5_prxtp_tca_2605_swtchbld100_13-Oct-2017_15:13:59.mat';
% 
% pth='/d2/nsteiger/output-da/hydroclimate/'; 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');
% 
% 
% xa_nino=reshape(Xa_m(id_nino(1):id_nino(2),:),12*reconYrs,1);
% xa_nino_ens=reshape(permute(Xa_ens(id_nino(1):id_nino(2),:,:),[1 3 2]),12*reconYrs,length(p_yrs));
% 
% xa_amo=Xa_m(id_amo(1):id_amo(2),:)';
% xa_amo_ens=permute(squeeze(Xa_ens(id_amo(1):id_amo(2),:,:)),[2 1]);
% 
% % Get rid of all previous variables from annual reconstruction
% clearvars -except xa_nino xa_nino_ens xa_amo xa_amo_ens
% 
% %% DROUGHT ATLASES: USE JJA INSTEAD OF ANNUAL MEAN RECONSTRUCTION FOR PDSI
% 
% 
% efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2590_swtchbld100_27-Sep-2017_15:45:13.mat';
% pth='/d2/nsteiger/output-da/hydroclimate/'; load([pth,efl])
% 
% 
% % What years to look at?
% analyrs1=800:1600;
% analyrs2=1601:1900;
% [~,~,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;
% 
% 
% 
% % Expand out PDSI land-only fields from reconstruction
% % NEED TO MAKE 'LNDIDX' AVAILABLE IN THE OUTPUT FILE!!! MAYBE LOAD IT AS A GLOBAL VARIABLE SOMEWHERE? OR JUST LOAD IT IN PARAMS SCRIPT?
% load(['/d1/nsteiger/climate-data/cesm-lme/cesm_lme_010_pdsi_output_biascorr_AWC_c_7_u2_gcm.mat'],'lndidx')
% 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]);
% Xa_pdsi_lndonly=Xa_m(id_pdsilnd(1):id_pdsilnd(2),:);
% Xa_pdsi_f=NaN(length(lat)*length(lon),reconYrs);
% Xa_pdsi_f(lndidx,:)=Xa_pdsi_lndonly;
% %Xa_pdsi_f=reshape(Xa_pdsi_f,length(lat),length(lon),reconYrs);
% 
% % Analysis options
% anal_da='y'; % drought atlas verification?
% plot_p='y'; % plot proxies?
% make_figs='n'; % Save figures?
% 
% 
% %==============================
% % ANALYZE THE DROUGHT ATLASES?
% %==============================
% if anal_da=='y'
% 
% %------------------------------------------
% % Load each drought atlas and interpolate to the resolution of the reconstruction
% %------------------------------------------
% 
% 
% if mon_avg_o==6 && mon_avg_f==8
% 
%    %======================
%    %======  NADA  ========
%    %======================
% 
%    da=permute(ncread('/d1/nsteiger/drought-atlas/nada_hd2_cl.nc','pdsi'),[2 3 1]);
%    %nada(nada<-99)=NaN; % put nans in fill spots
%    lat_da=ncread('/d1/nsteiger/drought-atlas/nada_hd2_cl.nc','lat');
%    lon_da=ncread('/d1/nsteiger/drought-atlas/nada_hd2_cl.nc','lon');
%    yrs_da=ncread('/d1/nsteiger/drought-atlas/nada_hd2_cl.nc','time');
% 
%    if yrs_da(1)<=analyrs1(1)
% 
%    % Coarsen NADA to match climate model grid
%    ndm=load('nada_mask.mat');
%    % create grids for interpolation
%    nlns=find(lon_da<0);lon_da(nlns)=lon_da(nlns)+360;
%    [xn,yn]=meshgrid(lon_da,lat_da);
%    [xc,yc]=meshgrid(lon,lat);
%    % note masks of both data types
%    nm=find(~isnan(da(:,:,end)));
%    cm=find(ndm.sftlf==1);
%    da(da<-99)=NaN; % put nans in fill spots but after determining mask
% 
%    da_i=nan(length(lat),length(lon),size(da,3));
%    for i=1:size(da,3)
%       % linear interpolation and extrapolation for irregular data
%       da1=da(:,:,i);
%       F = scatteredInterpolant(xn(nm),yn(nm),da1(nm),'linear','linear');
%       vq=F(xc(cm),yc(cm)); % do the interpolation
%       % expand out the interpolated data points onto the climate model mask
%       vd=nan(length(lat),length(lon));
%       vd(cm)=vq; % full interpolated field
%       da_i(:,:,i)=vd;
%    end
% 
% 
%    % get years of overlap
%    %[~,ir,ih] = intersect(r_o:r_f,time);
%    % Only look at a certain segment of the reconstructions
%    %analyrs=1800:1899;
%    [~,~,ir] = intersect(analyrs1,r_o:r_f); % indices of recon segment
%    [~,~,id] = intersect(analyrs1,yrs_da); % indices of DA segment
% 
%    % PULL OUT SPATIAL MEAN OVER AMERICAN SOUTHWEST
%    load('/d1/nsteiger/climate-data/cesm-lme/Nrth_A_Swthwst_mask.mat')
%    nasw(nasw==0)=NaN;
%    Xv=reshape(Xa_pdsi_f,length(lat),length(lon),reconYrs);
%    Xv_msk=bsxfun(@times,Xv,nasw);
%    A=cosd(repmat(lat,[1 length(lon)]));
%    Xa_di=wmean_a(Xv_msk,A); % RECONSTRUCTION DROUGHT INDEX
% 
%    DA_msk=bsxfun(@times,da_i,nasw);
%    DA_di=wmean_a(DA_msk,A); % DROUGHT ATLAS DROUGHT INDEX
% 
%    % Simple plot of spatial means (corr = 0.875)
%    %figure
%    %plot(analyrs,DA_di(id),'linewidth',2)
%    %hold on
%    %plot(analyrs,Xa_di(ir),'linewidth',2)
%    %legend('NADA','Recon')
%    %xlim([analyrs(1) analyrs(end)])
%    %
%    %if  make_figs=='y'
%    % print('-dpng',['./figs/timeseries_nada_recon_analyrs_' num2str(analyrs(1)) '_' num2str(analyrs(end)) '_' fln(1:end-13) '.png',])
%    %end
% 
% 
% %   xa_v1=Xa_di(ir1);
% %   xa_v2=Xa_di(ir2);
% %
% %   p1=gevfit(xa_v1-mean(xa_v2)); % subtract mean of post-1600
% %   x=-3:0.01:3;
% %   Yc1 = gevpdf(x,p1(1),p1(2),p1(3));
% %
% %   p2=gevfit(xa_v2-mean(xa_v2)); % subtract mean of post-1600
% %   Yc2 = gevpdf(x,p2(1),p2(2),p2(3));
% %
% %
% %   figure
% %   histogram(xa_v1,30)
% %   hold on
% %   histogram(xa_v2,30)
% %   title('PDSI NASW')
% %   legend('Pre-1600','Post-1600')
% %
% %   figure
% %   plot(x,Yc1,'linewidth',2)
% %   hold on
% %   plot(x,Yc2,'linewidth',2)
% %   title('PDSI NASW')
% %   legend('Pre-1600','Post-1600')
% 
% 
%    else
%    disp('NADA skipped...')
%    end % if within analysis range
% 
% end % if JJA
% 
% end % if: analyze drought atlases
% 
% 
% 
% 
% 
% %===========================================================================
% 
% 
% %% COMPUTE THE DROUGHTS
% 
% % Drought = negative anomalies around 11 year running mean
% 
% 
% % Find the droughts
% drghts=struct([]);
% 
% disp(['Computing megadroughts over the interval ' num2str(ir1(1)) '-' num2str(ir1(end)) ])
% dr_v=smooth(Xa_di(ir1),11); % extract drought variable & smooth it at 11 years running mean
% dr_v=dr_v-mean(dr_v); % Remove time mean
% neg_ids=find(dr_v<0);
% negpos=ones(length(analyrs1),1); % positive years will be ones
% negpos(neg_ids)=0;% set negative years as zeros
% 
% i=1;j=1;flag=0;
% while i+1 <= length(analyrs1)
%   % 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(analyrs1)
% 	      d_end=i-1;
% 	      drghts(j).strtend=[d_start d_end];
% 	      drghts(j).dsum=sum(dr_v(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];
%       drghts(j).dsum=sum(dr_v(d_start:d_end));
%       
%       j=j+1;
%   end
%   
%   i=i+1;
% end
% 
% 
% 
% % Rank the droughts, find the most severe
% d=zeros(length(drghts),3);
% for kk=1:length(drghts)
%   d(kk,1:2)=drghts(kk).strtend;
%   d(kk,3)=drghts(kk).dsum;
% 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,:)=[];  
% 
% 
% % Sort the droughts according to the most severe drought density (sum of
% % length and intensity)
% ds=sortrows(d,3);

%save('drought_analysis_out1.mat','-v7.3')
% Uncomment above to compute all the prep work...
load('drought_analysis_out1.mat')
make_figs='y';

%--------------------------------
% PROBABILISTIC DROUGHT ANALYSIS
%--------------------------------


% What years to look at?
analyrs3=1601:1925;
[~,~,ir3] = intersect(analyrs3,r_o:r_f); % indices of recon segment
% monthly indices
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;

% Pull out an extreme percentile
enso_mdrn_pdf=xa_nino_ens(irm3_a:irm3_b,:);
%em5=prctile(enso_mdrn_pdf(:),1);
em50=prctile(enso_mdrn_pdf(:),50);

amo_mdrn_pdf=xa_amo_ens(ir3,:);
%em5=prctile(enso_mdrn_pdf(:),1);
am50=prctile(amo_mdrn_pdf(:),50);

cldfrac_d=zeros(size(ds,1),1);
cldfrac_d_a=zeros(size(ds,1),1);
smz=200; % sample size
cldfrac_d_r=zeros(size(ds,1),smz);
cldfrac_d_a_r=zeros(size(ds,1),smz);
for i=1:size(ds,1)

  ovrlprng=rng_intr([r_o,r_f],[ds(i,1)+analyrs1(1)-1,ds(i,2)+analyrs1(1)-1]);
  irm_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
  irm_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;

  [~,~,irab] = intersect((ds(i,1)+analyrs1(1)-1):(ds(i,2)+analyrs1(1)-1),r_o:r_f); 
  
  % Pull out distributions during this drought
  enso_pdf=xa_nino_ens(irm_a:irm_b,:); 
  amo_pdf=xa_amo_ens(irab,:); 


  %x=(min(enso_pdf(:))-1):0.01:(max(enso_pdf(:))+1);
  %pk1=fitdist(enso_pdf(:),'kernel');
  %Yk1=pdf(pk1,x); 
  %figure;plot(x,Yk1)
  %hold on

  % Fraction of ENSO and AMO states that are below threshold
  cldfrac_d(i)=length(find(enso_pdf(:)<em50))/length(enso_pdf(:));
  cldfrac_d_a(i)=length(find(amo_pdf(:)>am50))/length(amo_pdf(:));


  % Randomly sample droughts of this same length from a reference period
  rsmn1=randi([irm3_a,irm3_b-length(irm_a:irm_b)+1],smz,1);
  rsmn=bsxfun(@plus,rsmn1,repmat(0:(length(irm_a:irm_b)-1),smz,1));
  rsmn11=randi([ir3(1),ir3(end)-length(irab)+1],smz,1);
  rsmn2=bsxfun(@plus,rsmn11,repmat(0:(length(irab)-1),smz,1));
  for j=1:smz
     renso_pdf=xa_nino_ens(rsmn(j,:),:);
     ramo_pdf=xa_amo_ens(rsmn2(j,:),:);
     
     cldfrac_d_r(i,j)=length(find(renso_pdf(:)<em50))/length(renso_pdf(:));
     cldfrac_d_a_r(i,j)=length(find(ramo_pdf(:)>am50))/length(ramo_pdf(:));
    

     %pk2=fitdist(renso_pdf(:),'kernel');
     %Yk2=pdf(pk2,x);
     %plot(x,Yk2)
  end

end

% p05=prctile(cldfrac_d_r,5,2);
% p95=prctile(cldfrac_d_r,95,2);
% 
% figure
% plot(cldfrac_d,'*')
% hold on
% plot(p05,'sk')
% plot(p95,'sk')
% title(['Fraction of ENSO states < median ' num2str(analyrs3(1)) '-' num2str(analyrs3(end))])
% xlabel('Drought rank')
% ylabel('Fraction')
% legend('Drought','rand 5th','rand 95th')
% 
% 
% p05a=prctile(cldfrac_d_a_r,5,2);
% p95a=prctile(cldfrac_d_a_r,95,2);
% 
% figure
% plot(cldfrac_d_a,'*')
% hold on
% plot(p05a,'sk')
% plot(p95a,'sk')
% title(['Fraction of AMO states > median ' num2str(analyrs3(1)) '-' num2str(analyrs3(end))])
% xlabel('Drought rank')
% ylabel('Fraction')
% legend('Drought','rand 5th','rand 95th')




figure
plot(cldfrac_d,'ro','markerfacecolor','r')
hold on
boxplot(cldfrac_d_r','PlotStyle','compact','symbol','','labelorientation','horizontal','color',[0.7 0.7 0.7])
plot(cldfrac_d,'ro','markerfacecolor','r')
title(['Nino3.4 colder than ' num2str(analyrs3(1)) '-' num2str(analyrs3(end)) ' median'])
xlabel('Drought severity rank','fontsize',18)
ylabel('Fraction of states in PDF','fontsize',18)
ylim([0.2 0.8])
hold off
set(gca,'fontsize',18)
if make_figs=='y';print(['./figs/enso_coldfrac_drought.png'],'-dpng','-painters');end



figure
plot(cldfrac_d_a,'ro','markerfacecolor','r')
hold on
boxplot(cldfrac_d_a_r','PlotStyle','compact','symbol','','labelorientation','horizontal','color',[0.7 0.7 0.7])
plot(cldfrac_d_a,'ro','markerfacecolor','r')
title(['AMO warmer than ' num2str(analyrs3(1)) '-' num2str(analyrs3(end)) ' median'])
xlabel('Drought severity rank','fontsize',18)
ylabel('Fraction of states in PDF','fontsize',18)
ylim([0.075 0.925])
hold off
set(gca,'fontsize',18)

if make_figs=='y';print(['./figs/amo_warmfrac_drought.png'],'-dpng','-painters');end





% What years to look at?
analyrs3=1601:1925;
[~,~,ir3] = intersect(analyrs3,r_o:r_f); % indices of recon segment
% monthly indices
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;

% Pull out an extreme percentile
enso_mdrn_pdf=xa_nino_ens(irm3_a:irm3_b,:);
enso_anc_pdf=xa_nino_ens(irm1_a:irm1_b,:);
enso_offset=mean(enso_mdrn_pdf(:))-mean(enso_anc_pdf(:));
%em5=prctile(enso_mdrn_pdf(:),1);
em50=prctile(enso_mdrn_pdf(:),50);

amo_mdrn_pdf=xa_amo_ens(ir3,:);
amo_anc_pdf=xa_amo_ens(ir1,:);
amo_offset=mean(amo_mdrn_pdf(:))-mean(amo_anc_pdf(:));
%em5=prctile(enso_mdrn_pdf(:),1);
am50=prctile(amo_mdrn_pdf(:),50);




cldfrac_d=zeros(size(ds,1),1);
cldfrac_d_a=zeros(size(ds,1),1);
smz=200; % sample size
cldfrac_d_r=zeros(size(ds,1),smz);
cldfrac_d_a_r=zeros(size(ds,1),smz);
for i=1:size(ds,1)

  ovrlprng=rng_intr([r_o,r_f],[ds(i,1)+analyrs1(1)-1,ds(i,2)+analyrs1(1)-1]);
  irm_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
  irm_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;

  [~,~,irab] = intersect((ds(i,1)+analyrs1(1)-1):(ds(i,2)+analyrs1(1)-1),r_o:r_f); 
  
  % Pull out distributions during this drought
  enso_pdf=xa_nino_ens(irm_a:irm_b,:);enso_pdf=enso_pdf(:); 
  amo_pdf=xa_amo_ens(irab,:);amo_pdf=amo_pdf(:); 

  % Adjust mean of distribution to 800-1600 period
  enso_pdf=enso_pdf+enso_offset;
  amo_pdf=amo_pdf+amo_offset;



  %x=(min(enso_pdf(:))-1):0.01:(max(enso_pdf(:))+1);
  %pk1=fitdist(enso_pdf(:),'kernel');
  %Yk1=pdf(pk1,x); 
  %figure;plot(x,Yk1)
  %hold on

  % Fraction of ENSO and AMO states that are below threshold
  cldfrac_d(i)=length(find(enso_pdf(:)<em50))/length(enso_pdf(:));
  cldfrac_d_a(i)=length(find(amo_pdf(:)>am50))/length(amo_pdf(:));


  % Randomly sample droughts of this same length from a reference period
  rsmn1=randi([irm3_a,irm3_b-length(irm_a:irm_b)+1],smz,1);
  rsmn=bsxfun(@plus,rsmn1,repmat(0:(length(irm_a:irm_b)-1),smz,1));
  rsmn11=randi([ir3(1),ir3(end)-length(irab)+1],smz,1);
  rsmn2=bsxfun(@plus,rsmn11,repmat(0:(length(irab)-1),smz,1));
  for j=1:smz
     renso_pdf=xa_nino_ens(rsmn(j,:),:);
     ramo_pdf=xa_amo_ens(rsmn2(j,:),:);
     
     cldfrac_d_r(i,j)=length(find(renso_pdf(:)<em50))/length(renso_pdf(:));
     cldfrac_d_a_r(i,j)=length(find(ramo_pdf(:)>am50))/length(ramo_pdf(:));
  end

end



figure
plot(cldfrac_d,'ro','markerfacecolor','r')
hold on
boxplot(cldfrac_d_r','PlotStyle','compact','symbol','','labelorientation','horizontal','color',[0.7 0.7 0.7])
plot(cldfrac_d,'ro','markerfacecolor','r')
title(['Nino3.4 colder than ' num2str(analyrs3(1)) '-' num2str(analyrs3(end)) ' median'])
xlabel('Drought severity rank','fontsize',18)
ylabel('Fraction of states in PDF','fontsize',18)
ylim([0.2 0.8])
hold off
set(gca,'fontsize',18)

if make_figs=='y';print(['./figs/enso_coldfrac_drought_offset.png'],'-dpng','-painters');end



figure
plot(cldfrac_d_a,'ro','markerfacecolor','r')
hold on
boxplot(cldfrac_d_a_r','PlotStyle','compact','symbol','','labelorientation','horizontal','color',[0.7 0.7 0.7])
plot(cldfrac_d_a,'ro','markerfacecolor','r')
title(['AMO warmer than ' num2str(analyrs3(1)) '-' num2str(analyrs3(end)) ' median'])
xlabel('Drought severity rank','fontsize',18)
ylabel('Fraction of states in PDF','fontsize',18)
ylim([0.075 0.925])
hold off
set(gca,'fontsize',18)

if make_figs=='y';print(['./figs/amo_warmfrac_drought_offset.png'],'-dpng','-painters');end
















% %-------------------------------------------------------------------
% % ENSO mode identification based on deviations from the X year mean
% %-------------------------------------------------------------------
% 
% % Mean determined by smoothed time series
% %xa_nino_mn=smooth(xa_nino,30*12,'rlowess');
% %xa_nino_mn=smooth(xa_nino,100*12,'lowess'); % centennial time scale mean state changes
% xa_nino_mn=smooth(xa_nino,300*12,'lowess'); % centennial time scale mean state changes
% 
% % Mean determined by time average over some period
% %ovrlprng=rng_intr([r_o,r_f],[1601,2000]);
% %ovrlprng=rng_intr([r_o,r_f],[analyrs(1),analyrs(end)]);
% %irm_a=find(abs(ovrlprng(1)-mon_r)<1e-6);irm_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;
% %xa_nino_mn=mean(xa_nino(irm_a:irm_b));
% 
% % Remove long-time scale mean state
% xa_nino_NM=xa_nino-xa_nino_mn;
% 
% mon_r=r_o:(1/12):(r_f+11/12);
% 
% 
% % PLOT NINO TIME SERIES AND TREND REMOVAL
% figure
% hold on
% plot(mon_r,xa_nino,'linewidth',2)
% plot(mon_r,xa_nino_mn,'linewidth',2)
% %xlim([r_o,r_f])
% xlim([analyrs(1),analyrs(end)])
% %legend([xt,h.mainLine],{'NADA','Recon'})
% %text(0.05,0.05,['r = ' sprintf('%0.2f',rp) ', CE = ' sprintf('%0.2f',cet)...
% %    ', CRPSS = ' sprintf('%0.2f',m_crpss)],'units','normalized','fontsize',14)
% %title('North American Southwest Drought')
% xlabel('Years','fontsize',18)
% ylabel('Nino 3.4 temperature anomalies','fontsize',18)
% set(gca,'fontsize',18)
% hold off
% box on
% %fig = gcf;fig.PaperUnits = 'inches';
% %fig.PaperPosition = [0.1 2 8 4]; % [left bottom width height]
% 
% 
% % above and below thresholds
% %ln_ids=find(xa_nino<-0.5);
% %en_ids=find(xa_nino>0.5);
% %nt_ids=find(xa_nino<=0.5 & xa_nino>=-0.5);
% 
% % just positive or negative (EL or LA)
% ln_ids=find(xa_nino_NM<0);
% en_ids=find(xa_nino_NM>0);
% 
% 
% lnino=zeros(length(xa_nino_NM),1);
% lnino(ln_ids)=1;
% 
% enino=zeros(length(xa_nino_NM),1);
% enino(en_ids)=1;
% 
% prct_ln_d=zeros(size(ds,1),1);
% prct_en_d=zeros(size(ds,1),1);
% for i=1:size(ds,1)
% 
%   ovrlprng=rng_intr([r_o,r_f],[ds(i,1)+analyrs(1)-1,ds(i,2)+analyrs(1)-1]);
%   irm_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
%   irm_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;
% 
%   prct_ln_d(i)=sum(lnino(irm_a:irm_b))/length(irm_a:irm_b);
%   prct_en_d(i)=sum(enino(irm_a:irm_b))/length(irm_a:irm_b);
% 
%   %prct_ln_d(i)=sum(lnino(ds(i,1)*12:((ds(i,2)+1)*12-1)))/length(ds(i,1)*12:((ds(i,2)+1)*12-1));
%   %prct_en_d(i)=sum(enino(ds(i,1)*12:((ds(i,2)+1)*12-1)))/length(ds(i,1)*12:((ds(i,2)+1)*12-1));
% end
% 
% %%
% 
% %dl=15; % top 'd' droughts
% dl=size(ds,1);
% 
% figure
% %plot(prct_ln(1:d)*100,'o','color',[0 0.4471 0.7412],'linewidth',2)
% hold on
% plot(prct_ln_d(1:dl)*100,'s','color',[0 0.4471 0.7412],'markerfacecolor',[0 0.4471 0.7412],'linewidth',2)
% plot([1 dl],[length(ln_ids)/length(xa_nino) length(ln_ids)/length(xa_nino)]*100,'--','color',[0 0.4471 0.7412])
% %plot(prct_en(1:d)*100,'o','color',[0.8510 0.3255 0.0980],'linewidth',2)
% plot(prct_en_d(1:dl)*100,'s','color',[0.8510 0.3255 0.0980],'markerfacecolor',[0.8510 0.3255 0.0980],'linewidth',2)
% plot([1 dl],[length(en_ids)/length(xa_nino) length(en_ids)/length(xa_nino)]*100,'--','color',[0.8510 0.3255 0.0980])
% xlim([1 dl])
% legend('La Nina','Recon avg','El Nino','Recon avg','location','eastoutside')
% xlabel('Drought severity rank','fontsize',14)
% ylabel('Monthly percent occurrence','fontsize',14)
% set(gca,'fontsize',14)
% box on
% %fig = gcf;fig.PaperUnits = 'inches';
% %fig.PaperPosition = [0.5 2 7 6]; % [left bottom width height]
% 
% 
% 













% %=============================================
% % Add a patches for top five times from recon
% %=============================================
% 
% x=r_o:r_f;
% y=smooth(nasw_d(l).indx,11); % using 'l' from above
% %h=shadedErrorBar(r_o:r_f,y,2*sigma_t,{'color',[0.494 0.184 0.556]},0);
% %set(h.mainLine,'linewidth',2)
% hx=plot(x,y,'color',[0.494 0.184 0.556],'linewidth',2);
% hold off
% set(gca,'fontsize',14)
% 
% legend([ht,hx],...
%     {'Truth','Recon.'},'location','southeast')
% 
% xlabel('Years','fontsize',14)
% ylabel('N. American Southwest PDSI','fontsize',14)
% 
% for i=1:5
%     p=patch(r_o+[ds(i,1) ds(i,2) ds(i,2) ds(i,1)],[yMin yMin yMax yMax],[0.494 0.184 0.556]);
%     p.FaceAlpha=0.25;p.EdgeColor='none';
% end
% 
% %axis([r_o,r_f,-0.8,0.8])
% xlim([r_o,r_f])
% 
% rp=corr(y,yt);
% cet=nan_ce(yt,y);
% text(0.05,0.05,['r = ' sprintf('%0.2f',rp) ', CE = ' sprintf('%0.2f',cet)],'units','normalized','fontsize',14)
% 
% box on
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0.1 2 8 4]; % [left bottom width height]
% 
% yr1=(ds(1,1)+r_o);
% yr2=(ds(1,2)+r_o);
% 
% 
% print(['/home/nsteiger/MATLAB/annual-da/hydro/pdsi_5droughts_11yravg_' num2str(yr1) '_' num2str(yr2) fln(1:end-4) '.pdf'],'-dpdf','-painters')
% 


