% 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_16-Nov-2017_16:03:43.mat';
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2601_swtchbld100_16-Nov-2017_16:11:22.mat';
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2590_swtchbld100_16-Nov-2017_16:20:57.mat';
% INDEX ONLY RECONSTRUCTION WITH FULL ENSEMBLES
%efl='cesm_lme010_r12000_p998_state_2ozn_avg_AprMar_prxydtst_5_prxtp_tca_2605_swtchbld100_16-Nov-2017_15:57: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');



% 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: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;


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


xa_v1=Xa_m(id_amo(1),ir1);
xa_v2=Xa_m(id_amo(1),ir2);

p1=gevfit(xa_v1-mean(xa_v2)); % subtract mean of post-1600
x=-1:0.01:1;
Yc1 = gevpdf(x,p1(1),p1(2),p1(3));

p2=gevfit(xa_v2-mean(xa_v2)); % subtract mean of pre-1600
Yc2 = gevpdf(x,p2(1),p2(2),p2(3));


figure
histogram(xa_v1,30)
hold on
histogram(xa_v2,30)
title('AMO')
legend('Pre-1600','Post-1600')

figure
plot(x,Yc1,'linewidth',2)
hold on
plot(x,Yc2,'linewidth',2)
title('AMO')
legend('Pre-1600','Post-1600')

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


xa_V=reshape(Xa_m(id_nino(1):id_nino(2),:),12*reconYrs,1);
    
xa_v1=xa_V(irm1_a:irm1_b);
xa_v2=xa_V(irm2_a:irm2_b);


p1=gevfit(xa_v1-mean(xa_v2)); % subtract mean of post-1600
x=-2.5:0.01:2.5;
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('Nino 3.4')
legend('Pre-1600','Post-1600')

figure
plot(x,Yc1,'linewidth',2)
hold on
plot(x,Yc2,'linewidth',2)
title('Nino 3.4')
legend('Pre-1600','Post-1600')



%======================
% LOCAL TEMPERATURE
%======================

load('/d1/nsteiger/climate-data/cesm-lme/Nrth_A_Swthwst_mask.mat')
nasw(nasw==0)=NaN;
Xv=reshape(Xa_m(id_ts(1):id_ts(2),:),length(lat),length(lon),reconYrs);
Xv_msk=bsxfun(@times,Xv,nasw);
A=cosd(repmat(lat,[1 length(lon)]));
Xv_ti=wmean_a(Xv_msk,A); % RECONSTRUCTION DROUGHT INDEX


xa_v1=Xv_ti(ir1);
xa_v2=Xv_ti(ir2);

p1=gevfit(xa_v1-mean(xa_v2)); % subtract mean of pre-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 pre-1600
Yc2 = gevpdf(x,p2(1),p2(2),p2(3));


figure
histogram(xa_v1,30)
hold on
histogram(xa_v2,30)
title('Temperature NASW')
legend('Pre-1600','Post-1600')

figure
plot(x,Yc1,'linewidth',2)
hold on
plot(x,Yc2,'linewidth',2)
title('Temperature NASW')
legend('Pre-1600','Post-1600')


%figure
%plot(analyrs,Xv_ti(ir),'linewidth',2)
%xlim([analyrs(1) analyrs(end)])
%title('NASW temperature')
%if  make_figs=='y'
% print('-dpng',['./figs/timeseries_nada_recon_analyrs_' num2str(analyrs(1)) '_' num2str(analyrs(end)) '_' fln(1:end-13) '.png',])
%end


% GATHER INDICES FOR ALL THE DYNAMIC VARIABLES
%Xa_ti=Xv_ti;
Xa_amo=Xa_m(id_amo(1),:);
xa_nino=reshape(Xa_m(id_nino(1):id_nino(2),:),12*reconYrs,1);

% Get rid of all previous variables from annual reconstruction
clearvars -except Xv_ti Xa_amo xa_nino

%% 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';
efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2590_swtchbld100_16-Nov-2017_16:20:57.mat';
pth='/d2/nsteiger/output-da/hydroclimate/'; load([pth,efl])


analyrs=800:2000;
[~,~,ir] = intersect(analyrs,r_o:r_f); % indices of recon segment
disp(['Drought analysis years ' num2str(analyrs(1)) '-' num2str(analyrs(end))])


% 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)<=analyrs(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(analyrs,r_o:r_f); % indices of recon segment
   [~,~,id] = intersect(analyrs,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



%===========================
% BOTH TEMPERATURE AND PDSI
%===========================

figure
plot(analyrs,zscore(Xv_ti(ir)),'linewidth',2)
hold on
plot(analyrs,zscore(Xa_di(ir)),'linewidth',2)
legend('T','PDSI')
xlim([analyrs(1) analyrs(end)])
title('NASW T & PDSI')



%========================
% T, PDSI, AMO, NINO
%========================

sm=24;mthd='lowess'; % 2 year smooth
Xa_n_sm=smooth(xa_nino,sm,mthd);
ovrlprng=rng_intr([r_o,r_f],[analyrs(1),analyrs(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;

figure
subplot(2,1,1)
plot(analyrs,zscore(Xv_ti(ir)),'linewidth',2)
hold on
plot(analyrs,zscore(Xa_di(ir)),'linewidth',2)
plot(analyrs,zscore(Xa_amo(ir)),'linewidth',2)

subplot(2,1,2)
plot(zscore(Xa_n_sm(irm1_a:irm1_b)),'linewidth',2)




%========================
% CROSS CORRELATION
%========================

[r1,l1]=xcorr(zscore(Xa_di(ir)),zscore(Xa_amo(ir))',30);
figure;plot(l1,r1)
xlabel('Lag (yrs)')
title('xcorr(PDSI,AMO)')



[r2,l2]=xcorr(zscore(Xa_di(ir)),zscore(Xv_ti(ir)),30);
figure;plot(l2,r2)
xlabel('Lag (yrs)')
title('xcorr(PDSI,T)')

a_nino=mon2ann(xa_nino,4,3);
[r3,l3]=xcorr(zscore(Xa_di(ir(1:end-1))),zscore(a_nino(ir(1:end-1))),30);
figure;plot(l3,r3)
xlabel('Lag (yrs)')
title('xcorr(PDSI,Nino 3.4)')



%===========================================================================


%% COMPUTE THE DROUGHTS

% Drought = negative anomalies around 11 year running mean

analyrs=800:1600; % Just look at megadrought period
[~,~,ir] = intersect(analyrs,r_o:r_f); % indices of recon segment

% Find the droughts
drghts=struct([]);

dr_v=smooth(Xa_di(ir),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(analyrs),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(analyrs)
  % 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(analyrs)
	      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('southwest_droughts_output.mat','-v7.3')

error('stop here')

%-------------------------------------------------------------------
% 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]






































































































% Compare the DA reconstruction with the drought atlases over a common interval



% % Pull out drought index
% Xa_gmt_m=Xa_m(id_nasw(1),:)';
% sigma_t=Xa_sigma(id_nasw(1),:);
% 
% %------------------------------
% % Scaled reconstruction state
% %------------------------------
% 
% % Global mean from update step
% xa_gmt=Xa_gmt_m-mean(Xa_gmt_m);
% 
% 
% % North American Southwest region: lat = 25-->42.5, lon = 125-->105 W
% 
% % Load data and manipulate
% % flip lat and time
% pdsi=flip(flip(permute(ncread('/d1/nsteiger/drought-atlas/NADAv2-2008.nc','PDSI'),[2 1 3]),1),3);
% lat=flip(ncread('/d1/nsteiger/drought-atlas/NADAv2-2008.nc','lat'));
% lon=ncread('/d1/nsteiger/drought-atlas/NADAv2-2008.nc','lon');
% time=flip(ncread('/d1/nsteiger/drought-atlas/NADAv2-2008.nc','time'));
% 
% % Just take over the region
% naswpdsi=pdsi(4:11,6:14,:);
% A=cosd(repmat(lat(4:11),[1 length(6:14)]));
% % Compute the index
% nasw_pdsi_idx=wmean_a(naswpdsi,A);
% 
% % make recon have same reference mean as observational data
% xa_gmt=xa_gmt+mean(nasw_pdsi_idx);
% %Xa_gmt_rm=Xa_gmt_rm+mean(tas_had);
% 
% % Bias-correct mean and variance of PDSI?
% %disp('Variance matching the PDSI variables...')
% %xa_gmt=bias_correct(nasw_pdsi_idx',xa_gmt','n');
% %xa_gmt=detrend(xa_gmt);
% 
% % get years of overlap
% %[~,ir,ih] = intersect(r_o:r_f,time);
% % Only look at a certain segment of the reconstructions
% % get years of overlap
% %reconYrs=1;
% analyrs=1000:2000;
% [~,~,ia] = intersect(analyrs,r_o:r_f); % indices of recon analysis segment
% [~,~,in] = intersect(analyrs,time); % indices of NADA
% 
% % TIME SERIES CHECK
% rp=nancorr(detrend(xa_gmt(ia)),nasw_pdsi_idx(in));
% cet=nan_ce(nasw_pdsi_idx(in),detrend(xa_gmt(ia)));
% 
% % compute mean crpss
% frcst_c=crps_m(xa_gmt(ia),sigma_t(ia),nasw_pdsi_idx(in));
% 
% xb_gmt=Xb_o(id_nasw(1),:);
% xb_gmt=xb_gmt-mean(xb_gmt);
% ref_c=crps_m(mean(xb_gmt),std(xb_gmt),nasw_pdsi_idx(in));
% % compute skill score
% m_crpss=1-frcst_c/ref_c;
% 
% %%
% 
% figure
% hold on
% h=shadedErrorBar(r_o:r_f,xa_gmt,2*sigma_t,'b',0);
% set(h.mainLine,'linewidth',2)
% xt=plot(time,nasw_pdsi_idx,'k','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('NASW JJA-PDSI','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]
% 
% if make_figs=='y';print(['/home/nsteiger/MATLAB/annual-da/hydro/nasw_pdsi_' fln(1:end-4) '.pdf'],'-dpdf','-painters');end
% 
% %%
% 
% figure
% hold on
% h=shadedErrorBar(r_o:r_f,xa_gmt,2*sigma_t,'b',0);
% set(h.mainLine,'linewidth',2)
% %xt=plot(time,nasw_pdsi_idx,'k','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('NASW JJA-PDSI','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]
% 
% if make_figs=='y';print(['/home/nsteiger/MATLAB/annual-da/hydro/nasw_pdsi_nonada_' fln(1:end-4) '.pdf'],'-dpdf','-painters');end
% 
% 
% %%
% % skill of smoothed versions
% sm=10;
% rp=nancorr(smooth(detrend(xa_gmt(ia)),sm),smooth(nasw_pdsi_idx(in),sm));
% cet=nan_ce(smooth(nasw_pdsi_idx(in),sm),smooth(detrend(xa_gmt(ia)),sm));
% 
% figure
% hold on
% h=shadedErrorBar(r_o:r_f,smooth(xa_gmt,sm),2*smooth(sigma_t,sm),'b',0);
% set(h.mainLine,'linewidth',2)
% xt=plot(time,smooth(nasw_pdsi_idx,sm),'k','linewidth',2);
% %xlim([r_o,r_f])
% xlim([analyrs(1),analyrs(end)])
% title('North American Southwest Drought')
% 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)
% xlabel('Years','fontsize',18)
% ylabel('NASW JJA-PDSI','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]
% 
% if make_figs=='y';print(['/home/nsteiger/MATLAB/annual-da/hydro/nasw_pdsi_sm_' fln(1:end-4) '.pdf'],'-dpdf','-painters');end
% 
% % COMPUTE THE DROUGHTS
% %%
% % skill of smoothed versions
% 
% figure
% hold on
% h=shadedErrorBar(r_o:r_f,smooth(xa_gmt,sm),2*smooth(sigma_t,sm),'b',0);
% set(h.mainLine,'linewidth',2)
% %xt=plot(time,smooth(nasw_pdsi_idx,sm),'k','linewidth',2);
% %xlim([r_o,r_f])
% xlim([analyrs(1),analyrs(end)])
% title('North American Southwest Drought')
% %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)
% xlabel('Years','fontsize',18)
% ylabel('NASW JJA-PDSI','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]
% 
% if make_figs=='y';print(['/home/nsteiger/MATLAB/annual-da/hydro/nasw_pdsi_sm_nonada_' fln(1:end-4) '.pdf'],'-dpdf','-painters');end
% 
% %%
% % COHERENCE
% 
% % [s, c, ph, ci, phi] = cmtm(x,y,dt,NW,qbias,confn,qplot);
% [f,coh,~,ci,~] = cmtm(xa_gmt(ia),nasw_pdsi_idx(in),1,8,0,0,0);
% 
% figure
% hold on
% plot(f,coh,'linewidth',2,'color',[0.4940 0.1840 0.5560])
% plot(f,ci,'--','linewidth',2,'color',[0.5 0.5 0.5])
% %plot(f,coh,'linewidth',2,'color',[0.4940 0.1840 0.5560])
% %plot(f,ci,'--','linewidth',2,'color',[0.5 0.5 0.5])
% legend('C(NADA,Recon.)','95% conf.')
% xlim([0 0.5])
% ylim([0 1])
% set(gca,'XTick',[0 0.1 0.2 0.3 0.4 0.5])
% set(gca,'XTickLabel',{'0','0.1 [10]','0.2 [5]','0.3 [3.3]','0.4 [2.5]','0.5 [2]'})
% box on
% title('NASW drought index','fontsize',18,'fontweight','bold')
% xlabel('Frequency (1/yrs), Period [yrs]','fontsize',14)
% ylabel('Coherence','fontsize',14)
% set(gca,'fontsize',14)
% hold off
% 
% if make_figs=='y';print(['/home/nsteiger/MATLAB/annual-da/hydro/nasw_pdsi_coh_' fln(1:end-4) '.pdf'],'-dpdf','-painters');end
% 
% %%
% 
% % COMPUTE THE DROUGHTS
% 
% % Drought = negative anomalies around 11 year running mean
% 
% % Find the droughts
% drghts=struct([]);
% 
% dr_v=smooth(xa_gmt(ia),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(analyrs),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(analyrs)
%     % 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(analyrs)
%                 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
% 
% % Sort the droughts according to the most severe drought density (sum of
% % length and intensity)
% ds=sortrows(d,3);
% 
% 
% %% Droughts from observations
% 
% % Find the droughts
% drghts=struct([]);
% 
% dr_v=smooth(nasw_pdsi_idx(in),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(analyrs),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(analyrs)
%     % 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(analyrs)
%                 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
% 
% % Sort the droughts according to the most severe drought density (sum of
% % length and intensity)
% ds_obs=sortrows(d,3);
% 
% 
% 
% 
% %% ENSO FREQUENCY OF OCCURENCE ALONG WITH DROUGHTS
% 
% % % Use different PSM-based Nino reconstruction (since Nino with VS-lite is
% % % way better than T and PDSI)
% % efl='cesm_lme009_r12000_p998_state_2wnez_avg_AprMar_prxydtst_4_prxtp_tca_stnd_y_itr15_swtchbld75_01-Jan-2017_21:32:40.mat';
% % pth='/d2/nsteiger/output-da/';
% % load([pth,efl])
% 
% 
% %------------------------------
% % Scaled reconstruction state
% %------------------------------
% 
% Xa_nino_m=reshape(Xa_m(id_nino(1):id_nino(2),:),12*reconYrs,1);
% xa_nino=Xa_nino_m-mean(Xa_nino_m);
% 
% sigma_t=reshape(Xa_sigma(id_nino(1):id_nino(2),:),12*reconYrs,1);
% 
% % monthly indices
% ovrlprng=rng_intr([r_o,r_f],[analyrs(1),analyrs(end)]);
% mon_r=r_o:(1/12):(r_f+11/12);
% mon_a=analyrs(1):(1/12):(analyrs(end)+11/12);
% 
% % get years of overlap
% ir1=find(abs(ovrlprng(1)-mon_r)<1e-6);
% ir2=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;
% 
% xa_nino=xa_nino(ir1:ir2);
% 
% %%
% % % ENSO mode identification
% % ln_ids=find(t_nino<-0.5);
% % en_ids=find(t_nino>0.5);
% % nt_ids=find(t_nino<=0.5 & t_nino>=-0.5);
% % 
% % lnino=zeros(length(t_nino),1);
% % lnino(ln_ids)=1;
% % 
% % enino=zeros(length(t_nino),1);
% % enino(en_ids)=1;
% % 
% % prct_ln=zeros(size(ds,1),1);
% % prct_en=zeros(size(ds,1),1);
% % for i=1:size(ds,1)
% %     prct_ln(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(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
% 
% % ENSO mode identification based on deviations from the X year mean
% %xa_nino_mn=smooth(xa_nino,30*12,'rlowess');
% xa_nino_mn=smooth(xa_nino,100*12,'lowess'); % centennial time scale mean state changes
% figure;plot(xa_nino);hold on; plot(xa_nino_mn)
% 
% figure
% hold on
% plot((ir1:ir2)/12,xa_nino,'linewidth',2)
% plot((ir1:ir2)/12,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]
% 
% if make_figs=='y';print(['/home/nsteiger/MATLAB/annual-da/hydro/nino_trend_removed_' fln(1:end-4) '.pdf'],'-dpdf','-painters');end
% 
% %%
% xa_nino=xa_nino-xa_nino_mn;
% 
% 
% 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);
% 
% lnino=zeros(length(xa_nino),1);
% lnino(ln_ids)=1;
% 
% enino=zeros(length(xa_nino),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)
%     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
% 
% %%
% 
% d=15; % top 'd' droughts
% 
% figure
% %plot(prct_ln(1:d)*100,'o','color',[0 0.4471 0.7412],'linewidth',2)
% hold on
% plot(prct_ln_d(1:d)*100,'s','color',[0 0.4471 0.7412],'markerfacecolor',[0 0.4471 0.7412],'linewidth',2)
% plot([1 d],[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:d)*100,'s','color',[0.8510 0.3255 0.0980],'markerfacecolor',[0.8510 0.3255 0.0980],'linewidth',2)
% plot([1 d],[length(en_ids)/length(xa_nino) length(en_ids)/length(xa_nino)]*100,'--','color',[0.8510 0.3255 0.0980])
% xlim([1 d])
% 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]
% %%
% if make_figs=='y';print(['/home/nsteiger/MATLAB/annual-da/hydro/nino_prctg_occurence_' fln(1:end-4) '.pdf'],'-dpdf','-painters');end
% 



%% PLOT PROXY AVAILABILITY

% clear
% 
% load(['/d1/nsteiger/proxy-data/PAGES2k_v2/v1_13_1/'...
%     'pages2kTSv1_13_1_unpack_AMA_allproxies_njs_31-Dec-2016_14:29:02.mat'])
% disp('Using PAGES2k temperature proxy data v 1.13.1')
% 
% % Convert proxy lons into 0 to 360
% nlns=find(p_lon<0);
% p_lon(nlns)=p_lon(nlns)+360;
% prx_i=proxy_ama;
% 
% prx_im=bsxfun(@minus,prx_i,nanmean(prx_i,1));
% prx_d=bsxfun(@rdivide,prx_im,nanstd(prx_im,0,1));
% 
% p_idx = find(ismember(archive,{'tree'}));
% prxs=prx_d(:,p_idx);
% avl_t=sum(~isnan(prxs),2);
% 
% p_idx = find(ismember(archive,{'coral','sclerosponge'}));
% prxs=prx_d(:,p_idx);
% avl_c=sum(~isnan(prxs),2);
% 
% p_idx = find(ismember(archive,{'glacier ice'}));
% prxs=prx_d(:,p_idx);
% avl_i=sum(~isnan(prxs),2);
% 
% p_idx = find(~ismember(archive,{'coral','tree','sclerosponge','marine sediment'}));
% prxs=prx_d(:,p_idx);
% avl_o=sum(~isnan(prxs),2);
% %%
% X=repmat(year,[4 1])';
% %Y=cat(2,avl_t,avl_c,avl_i,avl_o);
% Y=cat(2,avl_o,avl_i,avl_c,avl_t);
% 
% figure
% area(X,Y)
% xlim([0 2005])
% %xlim([1600 2005])
% %legend('Trees','Corals','Ice','Other','location','northwest')
% legend('Other','Ice','Corals','Trees','location','north')
% xlabel('Year','fontsize',14)
% ylabel('Proxy availability','fontsize',14)
% set(gca,'fontsize',14)
% box on
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0.1 2 8 4]; % [left bottom width height]
% 
% print(['/home/nsteiger/MATLAB/annual-da/proxy_availability_PAGES2kv1_13_1_v2.png'],'-dpng','-painters')

% 
% %%
% 
% figure
% 
% x=r_o:r_f;
% yt=smooth(nasw_d(l).indx,11); % using 'l' from above
% yMin=min(yt)-0.5;yMax=max(yt)+0.5;
% ht=plot(x,yt,'k','linewidth',2);
% hold on
% %title('PDSI droughts: red = truth, blue = recon, purple = overlap')
% ylim([yMin yMax])
% % Add a patches for top five times from truth
% for i=1:5
%     p=patch(r_o+[ds(i,1) ds(i,2) ds(i,2) ds(i,1)],[yMin yMin yMax yMax],'k');
%     p.FaceAlpha=0.25;p.EdgeColor='none';
% end
% 
% 
% 
% nasw_d=struct([]);
% jj=1;
% for j=1:length(state_tp)
%     
%     switch state_tp(j)
%         case 't' % temperature
%             continue
%         case 'u'
%             continue
%         case 'p'
%             continue
%         case 'g'
%             continue
%         case 'l' % normalize soil moisture (otherwise no need to normalize PDSI or SPEI)
%             id_y=id_slms;
%             var='slms';
%             od=Xma_m(id_y(1):id_y(2),:);
%             odt=Xt(id_y(1):id_y(2),:);
%             % standardize data
%             Xma_m(id_y(1):id_y(2),:)=zscore(od,0,2);
%             Xt(id_y(1):id_y(2),:)=zscore(odt,0,2);
%         case 'i'
%             id_y=id_spei12;
%             var='spei12';
%             % remove time mean
%             od=Xma_m(id_y(1):id_y(2),:);
%             odt=Xt(id_y(1):id_y(2),:);
%             Xma_m(id_y(1):id_y(2),:)=bsxfun(@minus,od,mean(od,2));
%             Xt(id_y(1):id_y(2),:)=bsxfun(@minus,odt,mean(odt,2));
%         case 'r'
%             id_y=id_pdsi;
%             var='pdsi';
%             % remove time mean
%             od=Xma_m(id_y(1):id_y(2),:);
%             odt=Xt(id_y(1):id_y(2),:);
%             Xma_m(id_y(1):id_y(2),:)=bsxfun(@minus,od,mean(od,2));
%             Xt(id_y(1):id_y(2),:)=bsxfun(@minus,odt,mean(odt,2));
%         case 'o'
%             continue
%         case 'n'
%             continue
%     end
%     
%     
%     %     % NaN fill non-land areas
%     %     Xf=NaN(length(lat)*length(lon),reconYrs);
%     %     Xf(lndidx,:)=Xma_m(id_y(1):id_y(2),:);
%     %     %Xf_r=reshape(Xf,length(lat),length(lon));
%     %
%     %     % APPLY SPATIAL MASK
%     %     nasw(nasw==0)=NaN;
%     %     Xf_msk=bsxfun(@times,Xf,nasw(:));
%     %     Xf_msk=reshape(Xf_msk,length(lat),length(lon),reconYrs);
%     %
%     %     % COMPUTE SPATIAL MEAN
%     %     A=cosd(repmat(lat,[1 cols]));
%     %     NASW_x=wmean_a(Xf_msk,A);
%     
%     
%     % NaN fill non-land areas
%     Xft=NaN(length(lat)*length(lon),reconYrs);
%     %Xft(lndidx,:)=Xt(id_y(1):id_y(2),:); % true state
%     Xft(lndidx,:)=Xma_m(id_y(1):id_y(2),:); % recon state
%     %Xf_r=reshape(Xf,length(lat),length(lon));
%     
%     % APPLY SPATIAL MASK
%     nasw(nasw==0)=NaN;
%     Xft_msk=bsxfun(@times,Xft,nasw(:));
%     Xft_msk=reshape(Xft_msk,length(lat),length(lon),reconYrs);
%     
%     % COMPUTE SPATIAL MEAN
%     A=cosd(repmat(lat,[1 cols]));
%     nasw_d(jj).var=var;
%     nasw_d(jj).indx=wmean_a(Xft_msk,A);
%     
%     jj=jj+1;
% end
% 
% % COMPUTE THE DROUGHTS FOR THE MC ITERATIONS
% 
% % WESTERN MEAN DROUGHT INDEX FROM PDSI
% 
% id_y=id_pdsi;
% var='pdsi';
% % NaN fill non-land areas
% Xf=NaN(length(lat)*length(lon),reconYrs,size(Xma_ens_f,3));
% Xf(lndidx,:,:)=Xma_ens_f(id_y(1):id_y(2),:,:);
% 
% %Xf_r=reshape(Xf,length(lat),length(lon),length(p_yrs),reconYrs);
% 
% % APPLY SPATIAL MASK
% nasw(nasw==0)=NaN;
% Xft_msk=bsxfun(@times,Xf,nasw(:));
% Xft_msk=reshape(Xft_msk,length(lat),length(lon),reconYrs,size(Xma_ens_f,3));
% 
% % COMPUTE SPATIAL MEAN
% A=cosd(repmat(lat,[1 cols]));
% nasw_d_ens=zeros(size(Xma_ens_f,3),reconYrs);
% for i=1:size(Xma_ens_f,3)
%     nasw_d_ens(i,:)=wmean_a(squeeze(Xft_msk(:,:,:,i)),A);
% end
% 
% sigma_mc_drght=std(nasw_d_ens,0,1)';
% 
% 
% % add uncertainty in quadrature
% sigma_t=sqrt(sigma_xa_drght.^2+sigma_mc_drght.^2);
% 
% 
% % We employ a drought definition similar to
% % that described in Herweijer et al. (2007), with a drought
% % commencing after two consecutive years of negative soil
% % moisture anomalies and continuing until two consecutive
% % years of positive anomalies (2S2E).
% % Droughts identified using the 2S2E definition were
% % ordered by creating a drought density rank. For each
% % drought period, the NASW index was summed from the
% % first to the last year of the drought. These values were
% % subsequently ranked by the negative value of the sum.
% % This drought density ranking was chosen over a purely
% % length-based ranking in order to incorporate both the
% % persistence and severity of each drought.
% 
% % Find droughts based on soil moisture variable
% %slms_id=find(strcmp('slms',{nasw_d.var})==1);
% % CORRELATIONS ARE USUALLY ABOVE 0.9, BUT DROUGHT DENSITIES DISAGREE BASED
% % ON THE INDEX
% 
% % Drought = negative anomalies around 11 year running mean
% 
% % Find the droughts
% drghts=struct([]);
% for k=1:length(nasw_d)
%     
%     drghts(k).varn=nasw_d(k).var;% note name of variable
%     dr_v=smooth(nasw_d(k).indx,11); % extract drought variable & smooth it at 11 years running mean
%     neg_ids=find(dr_v<0);
%     negpos=ones(reconYrs,1); % positive years will be ones
%     negpos(neg_ids)=0;% set negative years as zeros
%     
%     i=1;j=1;flag=0;
%     while i+1 <= reconYrs
%         % 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 > reconYrs
%                     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(k).var(j).strtend=[d_start d_end];
%             drghts(k).var(j).dsum=sum(dr_v(d_start:d_end));
%             
%             j=j+1;
%         end
%         
%         i=i+1;
%     end
%     
% end
% 
% 
% % Rank the droughts, find the most severe
% l=1; % NUMBERS DEPEND ON WHICH VARIABLE IS FIRST IN STATE VECTOR 1 = spei, 2 = pdsi
% d=zeros(length(drghts(l).var),3);
% for kk=1:length(drghts(l).var)
%     d(kk,1:2)=drghts(l).var(kk).strtend;
%     d(kk,3)=drghts(l).var(kk).dsum;
% end
% 
% % Sort the droughts according to the most severe drought density (sum of
% % length and intensity)
% ds=sortrows(d,3);
% 
% 
% 
% %=============================================
% % 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')
% 


