
% Make spatial composites during megadroughts 

addpath('../')


%===========================================================================
% LOAD CLIMATE INDICES

% Load full ensembles for the nino and amo variables
efl='cesm_lme010_r12000_p998_state_2ozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:24:35.mat';
% LOAD JJA FOR AMO ANALYSIS
%efl='cesm_lme010_r12000_p998_state_2ozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:11:15.mat';

pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';
%load([pth,efl],'r_o','r_f','state_tp','xmeta','Xa_m','Xa_ens','reconYrs','p_yrs')
load([pth,efl])


% What years to look at?
analyrs1=1000:1925;
analyrs2=800: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;
% 
% 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;



% k=find(state_tp=='o');
% xa_amo_ens=permute(squeeze(Xa_ens(xmeta{k}.id_X(1),:,:)),[2 1]);
% % Ensemble mean of AMO
% xa_amo_mn=mean(xa_amo_ens,2);
% amo_off=mean(mean(xa_amo_ens(ir2,:),2));
% 
% 
% k=find(state_tp=='2');
% xa_gmt_ens=permute(squeeze(Xa_ens(xmeta{k}.id_X(1),:,:)),[2 1]);
% 
% k=find(state_tp=='n'); % specific nino variable
% 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_nino_mn=reshape(Xa_m((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:),12*reconYrs,1);
% % Nino offset with respect to reference period
% nino_off=mean(mean(xa_nino_ens(irm2_a:irm2_b,:),2));
% 
% % Annual nino from monthly data (leaves year ends not quite correct)
% xa_nino_ann=mon2ann(xa_nino_mn,4,3);
% % Annual nino from spatial field
% %xa_nino_ann=nino(Xt2m,lat,lon,3,'n');



% Get annual Nino 3.4 value from annual spatial field
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:33:07.mat';
efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:40:23.mat';disp('Loading DJF data...')
load([pth,efl],'state_tp','xmeta','Xa_m')
k=find(state_tp=='t');
lat=xmeta{k}.lat;
lon=xmeta{k}.lon;
Xt2m=reshape(Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:),length(lat),length(lon),reconYrs);

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

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

% Get the ensemble mean; identical to the original reconstruction
Xv_msk=bsxfun(@times,Xt2m,lndmsk);
t2m_casa=wmean_a(Xv_msk,A); % RECONSTRUCTION DROUGHT INDEX




% 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';
%pth='./output-recon/';load([pth,efl],'xmeta')


% Load DJF PDSI
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:40:23.mat';
%load([pth,efl],'state_tp','xmeta','Xa_m')

% Load PDSI
k=find(state_tp=='q');
Xvar=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);


%--------------------------------------
% 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
[~,~,ip2] = intersect(analyrs2,result.times); % indices of recon segment

paicoE=result.resample.signals';
paicoM=mean(result.resample.signals);
xa_di_sa=paicoM;

%paicotm=paicotm-mean(paicotm)+mean(phydatm);
%paicote=paicote0-mean(paicotm)+mean(phydatm); % same time mean ensemble mean for comparison



%=====================================
% LOAD ENSEMBLE MEGADROUGHT ESTIMATES

% Load both south american and north america drought indices
%load('svmd_ens_90pagree_sa.mat')
%load('svmd_ens_90pagree_pampas.mat')
%load('pampas_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('pampas_megaD_ens_noNAproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('pampas_megaD_ens_noNASWproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')

% top here is the original, full proxy network
load('chile_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('chile_megaD_ens_s800_noNAproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('chile_megaD_ens_s800_noNASWproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('chile_megaD_ens_s800_chileonly.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')

% xa_di_sa=Xa_di;
% xa_di_ens_sa=Xa_di_ens;
% ds_sa=svmd;
% mdagr_sa=ens_drght;
% sa_yrs=analyrs;




%load('svmd_ens_50pagree.mat')
%load('svmd_ens_90pagree.mat')
load('nasw_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
xa_di_na=Xa_di;
xa_di_ens_na=Xa_di_ens;
ds_na=svmd;
mdagr_na=ens_drght;
na_yrs=analyrs;

%ds=svmd;


% %===========================================================================
% % LOAD FORCED CLIMATE RESPONSE ESTIMATES
% 
% % Full forcing
% %load('../megad-forcing/output_twoboxCM_full_06-Aug-2018_15:09:13.mat','gmt_fr','myrs')
% load('../ebm-held/output_twoboxCM_full_21-Sep-2018_13:29:49.mat','gmt_fr','myrs')
% 
% 
% % Plot mean forced response
% gmt_est_mn=mean(gmt_fr,2);
% [~,~,ig2] = intersect(analyrs2,myrs);% 1601-1925 reference 
% [~,~,ig3] = intersect(analyrs3,myrs); 
% [~,~,ig4] = intersect(analyrs4,myrs); 
% 
% % remove mean from ref
% gmt_est_ens=gmt_fr-mean(gmt_est_mn(ig2));
% gmt_est_f=gmt_est_mn-mean(gmt_est_mn(ig2));
% 
% % solar 
% load('../ebm-held/output_twoboxCM_sol_21-Sep-2018_13:29:03.mat','gmt_fr')
% %load('../megad-forcing/output_twoboxCM_solar_08-Aug-2018_13:47:28.mat','gmt_fr')
% gmt_est_mn=mean(gmt_fr,2);
% gmt_est_s=gmt_est_mn-mean(gmt_est_mn(ig2));
% 
% % volcanic 
% %load('../megad-forcing/output_twoboxCM_volc_08-Aug-2018_13:38:26.mat','gmt_fr')
% load('../ebm-held/output_twoboxCM_volc_21-Sep-2018_13:35:31.mat','gmt_fr')
% gmt_est_mn=mean(gmt_fr,2);
% gmt_est_v=gmt_est_mn-mean(gmt_est_mn(ig2));
% 
% % ghg
% %load('../megad-forcing/output_twoboxCM_ghg_08-Aug-2018_14:05:19.mat','gmt_fr')
% load('../ebm-held/output_twoboxCM_ghg_21-Sep-2018_13:32:35.mat','gmt_fr')
% gmt_est_mn=mean(gmt_fr,2);
% gmt_est_g=gmt_est_mn-mean(gmt_est_mn(ig2));
% 
% % ghg + solar
% load('../ebm-held/output_twoboxCM_solghg_26-Sep-2018_14:06:47.mat','gmt_fr')
% gmt_est_mn=mean(gmt_fr,2);
% gmt_est_sg=gmt_est_mn-mean(gmt_est_mn(ig2));
% 
% % volc + solar
% load('../ebm-held/output_twoboxCM_volcsol_26-Sep-2018_15:37:39.mat','gmt_fr')
% gmt_est_mn=mean(gmt_fr,2);
% gmt_est_vs=gmt_est_mn-mean(gmt_est_mn(ig2));
% 
% %===========================================================================



% Find continguous years
%[~,~,irn] = intersect(analyrs2,r_o:r_f); % indices of recon segment
%[ds_na]=contig(dyrs_na,irn);

%[~,~,irs] = intersect(sa_yrs,r_o:r_f); % indices of recon segment
%[~,~,irs] = intersect(sa_yrs,pr_o:pr_f); % indices of recon segment
%[ds_sa]=contig(dyrs_sa,irs);



% Droughts in the real data mean
ds0_n=mdrought(xa_di_na(ir1),1,'y');
ds0_s=mdrought(xa_di_sa(ip1),1,'y');
%ds0_n=mdrought(xa_di_na(ir2),1,'y');
%ds0_s=mdrought(xa_di_sa(ip2),1,'y');

% Droughts too frequent, throw out based on relative threshold
dcut=50;
np67=prctile(ds0_n(:,3),dcut);
rmv=find(ds0_n(:,3)>np67);
ds0_n(rmv,:)=[];

sp67=prctile(ds0_s(:,3),dcut);
rmv=find(ds0_s(:,3)>sp67);
ds0_s(rmv,:)=[];






% % AR null for mean
% addpath('../ar-fit/')
% ntr=1000;
% d0=xa_di_na(analyrs);
% [w,A,C,~,~,~]=arfit(d0,1,5); % explore up to 5, but it looks like all are 2nd order
% na_ar_ens=squeeze(arsim(w,A,C,[length(analyrs),ntr]));% generate base ensemble AR noise time series 
% 
% d0=xa_di_sa(analyrs);
% [w,A,C,~,~,~]=arfit(d0,1,5); % explore up to 5, but it looks like all are 2nd order
% sa_ar_ens=squeeze(arsim(w,A,C,[length(analyrs),ntr]));% generate base ensemble AR noise time series 
% 
% % Ensemble estimate of drought periods
% dsens_n=struct();
% dsens_s=struct();
% for i=1:ntr
%    % Compute the megadroughts
%    dsens_n.ds{i}=mdrought(na_ar_ens(:,i),1,'y');
%    dsens_s.ds{i}=mdrought(sa_ar_ens(:,i),1,'y');
% 
%    % Throw out droughts below Xth percentile from real data
%    rmv=find(dsens_n.ds{i}(:,3)>np67);
%    dsens_n.ds{i}(rmv,:)=[];
%    rmv=find(dsens_s.ds{i}(:,3)>sp67);
%    dsens_s.ds{i}(rmv,:)=[];
% end
% 
% 
% % Find which years there is agreement
% pd_nthr=zeros(ntr,1);
% pd_na=zeros(ntr,1);
% pd_sa=zeros(ntr,1);
% pd_bh=zeros(ntr,1);
% for i=1:ntr
% 
%    % North America
%    dn0=dsens_n.ds{i}(:,1:2); % pull out dates of droughts
%    dn0=dn0';dn0=dn0(:);
% 
%    % Fill out all years of droughts in vector
%    k=1;dn1=[];
%    while k<length(dn0)
%       dn1=cat(2,dn1,dn0(k):dn0(k+1));
%       k=k+2;
%    end
% 
%    % South America
%    ds0=dsens_s.ds{i}(:,1:2); % pull out dates of droughts
%    ds0=ds0';ds0=ds0(:);
% 
%    % Fill out all years of droughts in vector
%    k=1;ds1=[];
%    while k<length(ds0)
%       ds1=cat(2,ds1,ds0(k):ds0(k+1));
%       k=k+2;
%    end
% 
%    % 1/0 for kinds of drought
%    dtm=zeros(length(analyrs),1); 
%    dtm(dn1)=1;dtm(ds1)=1;ibth=intersect(dn1,ds1);dtm(ibth)=1;
%    df0=find(dtm==0); % no drought will still be zero
%    
%    % Probabilities for differention possibilities
%    pd_nthr(i)=length(df0)/length(analyrs);
%    pd_na(i)=length(dn1)/length(analyrs);
%    pd_sa(i)=length(ds1)/length(analyrs);
%    pd_bh(i)=length(ibth)/length(analyrs);
% 
% end


%----------------------------
% Probabilites for real data
%----------------------------

% North America
dn0=ds0_n(:,1:2)';dn0=dn0(:);

% Fill out all years of droughts in vector
k=1;dn1=[];
while k<length(dn0)
   dn1=cat(2,dn1,dn0(k):dn0(k+1));
   k=k+2;
end

% South America
ds0=ds0_s(:,1:2)';ds0=ds0(:);

% Fill out all years of droughts in vector
k=1;ds1=[];
while k<length(ds0)
   ds1=cat(2,ds1,ds0(k):ds0(k+1));
   k=k+2;
end

% Figure out indices of neither, NA, SA, both
inthr=setdiff(1:length(analyrs),union(dn1,ds1)); % neither
[~,ina,isa]=setxor(dn1,ds1); % indices of NA only and SA only
ibth=intersect(dn1,ds1); % both


% Probabilities for differention possibilities
pd_nthr0=length(inthr)/length(analyrs);
pd_na0=length(ina)/length(analyrs);
pd_sa0=length(isa)/length(analyrs);
pd_bh0=length(ibth)/length(analyrs);


% % Combine noise null
% pnl(:,1)=pd_bh;
% pnl(:,2)=pd_nthr;
% pnl(:,3)=pd_na;
% pnl(:,4)=pd_sa;
% 
% c0={'Both','None','NASW','Pampas'};
% c=categorical(c0,c0);
% 
% % Colors for boxplots
% cm2=cat(1,clrs.red,clrs.grey,clrs.blue,clrs.purple);
% 
% 
% figure
% hold on
% hl=boxplot(pnl,c,'symbol','','colors',cm2,'whisker',1); % whisker=1 => 2 sigma coverage
% set(hl,'linew',2)
% set(findobj(gcf,'-regexp','Tag','\w*Whisker'),'LineStyle','-')% whisker lines solid
% ax=gca;
% ylim([0 0.7])
% %line('xdata',ax.XLim,'ydata',[0 0],'linestyle',':','color',clrs.grey,'linewidth',2)
% ylabel('Probability','fontsize',16)
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% 
% h1=plot(1,pd_bh0,'*','color',clrs.black,'linewidth',2,'MarkerSize',10);
% h2=plot(2,pd_nthr0,'*','color',clrs.black,'linewidth',2,'MarkerSize',10);
% h3=plot(3,pd_na0,'*','color',clrs.black,'linewidth',2,'MarkerSize',10);
% h4=plot(4,pd_sa0,'*','color',clrs.black,'linewidth',2,'MarkerSize',10);
% 
% %legend([h1,h2,h3,h4],c0,'location','northeast')
% %legend([h1],'PHYDA','location','northeast')
% 
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 4 4]; % [left bottom width height]
% 
% 
% %print(['/home/nsteiger/Dropbox/sa-phyda/drght_prob_null_na_pampas_phyda_p' num2str(dcut) '.png'],'-dpng','-r300');



%-------------------------------------------------
% Composites over the different kinds of  periods
%-------------------------------------------------

% define input and output variable
inputf='cesm_latlon.nc';
%outputf=['phyda_megaD_tdjf_pdjfNoRef_cluster_comp.nc'];
%outputf=['phyda_megaD_tdjf_pdjf_cluster_comp.nc'];
%outputf=['phyda_noNA_megaD_tdjf_pdjf_cluster_comp.nc'];
%outputf=['phyda_noNASW_megaD_tdjf_pdjf_cluster_comp.nc'];
%outputf=['phyda_chileonly_megaD_tdjf_pdjf_cluster_comp.nc'];
outputf=['phyda_paico_megaD_tdjf_pdjf_cluster_comp_.nc'];
% grab lat/lon from model, put it into the output file
eval(['!ncks -A -v lat,lon ',inputf,' ',outputf]);

% Take analysis years
Xt2ma=Xt2m(:,:,analyrs);
Xpdsia=Xpdsi(:,:,analyrs);

cnm={'Neither','NA-only','SA-only','Both','NA','SA'};

% Composites over the possibilities
Xt2maM=bsxfun(@minus,Xt2ma,mean(Xt2ma,3));
X1=zeros(length(lat),length(lon),length(cnm));
X1(:,:,1)=mean(Xt2maM(:,:,inthr),3);
X1(:,:,2)=mean(Xt2maM(:,:,ina),3);
X1(:,:,3)=mean(Xt2maM(:,:,isa),3);
X1(:,:,4)=mean(Xt2maM(:,:,ibth),3);
X1(:,:,5)=mean(Xt2maM(:,:,dn1),3);
X1(:,:,6)=mean(Xt2maM(:,:,ds1),3);

XpaM=bsxfun(@minus,Xpdsia,mean(Xpdsia,3));
%XpaM=Xpdsia; disp('Dont remove PDSI mean...')
X2=zeros(length(lat),length(lon),length(cnm));
X2(:,:,1)=mean(XpaM(:,:,inthr),3);
X2(:,:,2)=mean(XpaM(:,:,ina),3);
X2(:,:,3)=mean(XpaM(:,:,isa),3);
X2(:,:,4)=mean(XpaM(:,:,ibth),3);
X2(:,:,5)=mean(XpaM(:,:,dn1),3);
X2(:,:,6)=mean(XpaM(:,:,ds1),3);


X2(isnan(X2))=1e20;


%figure
%s=load('coast');mlat=s.lat;mlon=s.long;
for i=1:length(cnm)

 %   figure
 %   %subplot(length(rnc),1,i)
 %   hold on
 %   h=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
 %   %h=contour(X,Y,cmp);
 %   h=pcolorm(lat,lon,X1(:,:,i));
 %   %set(h,'EdgeColor','none');
 %   plotm(mlat,mlon,'color',[0.5 0.5 0.5])
 %   %caxis([-3e2 3e2]) % Make the axis a log-axis?
 %   %colormap(cmap_br)
 %   %set(gca,'xtick',[],'ytick',[])
 %   title(['T2m ' num2str(i)])

   % SAVE COMPOSITE
   nccreate(outputf,['t_' cnm{i}],'Dimensions',{'lon',length(lon),'lat',length(lat)});
   ncwrite(outputf,['t_' cnm{i}],X1(:,:,i)');

   nccreate(outputf,['p_' cnm{i}],'Dimensions',{'lon',length(lon),'lat',length(lat)});
   ncwrite(outputf,['p_' cnm{i}],X2(:,:,i)');
   ncwriteatt(outputf,['p_' cnm{i}],'missing_value',1e20);

%   print(['/home/nsteiger/Dropbox/easterisle/ei_cluster_cmp' num2str(i) '_mslp.png'],'-dpng','-r300');
end




return





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

addpath('../')

% ALL PROXY TYPES, LMR DATABASE
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:49:45.mat';
efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:33:07.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])
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');



%-------------------------
% LOAD SPATIAL FIELDS
%-------------------------

% Temperature
k=find(state_tp=='t');
Xvar=Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);
Xt2m=reshape(Xvar,length(xmeta{k}.lat),length(xmeta{k}.lon),reconYrs);

k=find(state_tp=='n');i=3; % = Nino 3.4
xa_nino=reshape(Xa_m((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:),12*reconYrs,1);
xb_nino=reshape(Xb_o((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:),12*length(p_yrs),1);

k=find(state_tp=='z');i=3; % = Pacific (Schneider)
xa_itcz_pac=Xa_m((xmeta{k}.id_X(1)+i-1),:);
xb_itcz_pac=Xb_o((xmeta{k}.id_X(1)+i-1),:);
i=8; % = South America
xa_itcz_sa=Xa_m((xmeta{k}.id_X(1)+i-1),:);
xb_itcz_sa=Xb_o((xmeta{k}.id_X(1)+i-1),:);
i=1; % = Atlantic
xa_itcz_atl=Xa_m((xmeta{k}.id_X(1)+i-1),:);
xb_itcz_atl=Xb_o((xmeta{k}.id_X(1)+i-1),:);

k=find(state_tp=='o'); % AMO
xa_amo=Xa_m((xmeta{k}.id_X(1)),:);

clearvars -except xa_amo xb_itcz_sa xb_itcz_atl xb_itcz_pac xa_itcz_sa xa_itcz_atl xa_itcz_pac Xt2m xa_nino xb_nino


% 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';
efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:40:23.mat';
pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';load([pth,efl])
%pth='/d2/nsteiger/output-da/hydroclimate/'; load([pth,efl])


% What years to look at?
analyrs1=1000:1925;
%analyrs2=1601:1925;
%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);
% 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]);

% PDSI
k=find(state_tp=='q');
Xvar=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);

lat=xmeta{k}.lat;
lon=xmeta{k}.lon;


% Load full ensembles for the nino and amo variables
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/subrecons_5Jan18/';
load([pth,efl],'state_tp','xmeta','Xa_ens')


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

k=find(state_tp=='n'); % specific nino variable
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));

% Analysis options
prntfg='n'; % Save figures?



% %======================
% %======  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');
% 
% 
% % 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
% 
% 
% 
% % 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');
% % NASW = box of 235:255 E (125:105 W); 31:42 N (just the US, not Mexico)
% % Find nearest lat/lon equal to or outside bounds
% lonext=[235 255];latext=[31 42];
% ln_a=find((abs(lon-lonext(1))-min(abs(lon-lonext(1))))<1e-5);
% ln_b=find((abs(lon-lonext(2))-min(abs(lon-lonext(2))))<1e-5);
% lt_a=find((abs(lat-latext(1))-min(abs(lat-latext(1))))<1e-5);
% lt_b=find((abs(lat-latext(2))-min(abs(lat-latext(2))))<1e-5);
% nasw=zeros(size(sftlf));
% nasw(lt_a:lt_b,ln_a:ln_b)=sftlf(lt_a:lt_b,ln_a:ln_b);
% 
% nasw(nasw==0)=NaN;
% Xv_msk=bsxfun(@times,Xpdsi,nasw);
% A=cosd(repmat(xmeta{k}.lat,[1 length(xmeta{k}.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
% % Indices of the drought atlas
% [~,~,id1] = intersect(analyrs1,yrs_da); % indices of recon segment
% [~,~,id2] = intersect(analyrs2,yrs_da); % indices of recon segment
% [~,~,id3] = intersect(analyrs3,yrs_da); % indices of recon segment

% Load colors for plotting
load distinct_colors_20.mat


%=====================================
% LOAD ENSEMBLE MEGADROUGHT ESTIMATES

load('svmd_ens_90pagree_pampas.mat')

%load('svmd_ens_90pagree.mat')
%load('svmd_ens_50pagree.mat')
ds=svmd;



%=============================
% COMPOSITE OVER THE DROUGHTS
%=============================

% Pull out the SST values
omsk=sftlf(:);
oids=find(omsk==0);
lids=find(omsk==1);
Xsst=reshape(Xt2m,length(lat)*length(lon),reconYrs);
Xsst(lids,:)=NaN;
Xsst=reshape(Xsst,length(lat),length(lon),reconYrs);


% Remove referene period for composites
irc=ir1;
%irc=ir2;
T0=bsxfun(@minus,Xsst,mean(Xsst(:,:,irc),3));
P0=bsxfun(@minus,Xpdsi,mean(Xpdsi(:,:,irc),3));


% Get all droughts and composite over them
ds=svmd;
ird_all=[];
ird_ann=[];
Csst=zeros(length(lat),length(lon),length(ds));
Cpdsi=zeros(length(lat),length(lon),length(ds));
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;

  % Collect all months of drought
  ird_all=cat(2,ird_all,irm_a:irm_b);

  [~,~,irab] = intersect((ds(i,1)+analyrs1(1)-1):(ds(i,2)+analyrs1(1)-1),r_o:r_f); 
  % Collect all years of drought
  ird_ann=cat(1,ird_ann,irab);

  % Save composites over each drought
  Csst(:,:,i)=mean(T0(:,:,irab),3);
  Cpdsi(:,:,i)=mean(P0(:,:,irab),3);


end

% Composites over all the droughts
Csst_a=mean(T0(:,:,ird_ann),3);
Cpdsi_a=mean(P0(:,:,ird_ann),3);

savenc='n';

if savenc=='y'

   inputf='cesm_latlon.nc';

   % define the output variable
   outputf=['composite_megaD_mnref_' num2str(irc(1)) '_' num2str(irc(end)) '.nc'];
   %outputf='test_output.nc';inputf='ccsm4_latlon.nc';
   % grab lat/lon from model, put it into the output file
   eval(['!ncks -A -v lat,lon ',inputf,' ',outputf]);

   % Get composites of PDSI over the SOMs
   cmp_pdsi=permute(Cpdsi_a,[2 1]);
   cmp_pdsi(cmp_pdsi==NaN)=1e20;
   cnm=['cmp_pdsi'];

   tmp_ptrn=permute(Csst_a,[2 1]);
   tmp_ptrn(tmp_ptrn==NaN)=1e20;
   pnm=['cmp_sst'];

   nccreate(outputf,cnm,'Dimensions',{'lon',length(lon),'lat',length(lat)},'FillValue',1e20);
   ncwrite(outputf,cnm,cmp_pdsi);
   ncwriteatt(outputf,cnm,'missing_value',1e20);

   nccreate(outputf,pnm,'Dimensions',{'lon',length(lon),'lat',length(lat)},'FillValue',1e20);
   ncwrite(outputf,pnm,tmp_ptrn);
   ncwriteatt(outputf,pnm,'missing_value',1e20);

end

plotfigs='y';

if plotfigs=='y'
   % Plot the composites
   for j=1:length(ds)
      varX=Csst(:,:,j);
      figure
      hold on
      s=load('coast');mlat=s.lat;mlon=s.long;
      h=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
      %h=axesm('MapProjection','wagner4','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[120 33]);
      varZ=cat(2,varX,varX(:,1));
      pcolorm(lat,[lon; 360],varZ);
      %C=importdata('redblue_256_rgb.txt');
      %colormap(flip(C,1)/255)
      C=load('./colormaps/NCV_blue_red.mat');
      colormap(C.cmap)
      caxis([-1 1])
      set(gca,'fontsize',16)
      set(get(colorbar('location','southoutside'),'xlabel'),'string','Temperature anom. (deg C)','fontsize',16)
      cbh = findobj( 0, 'tag', 'Colorbar' );
      %set(cbh,'YTick',[-1:0.2:1])
      plotm(mlat,mlon,'color',[0.5 0.5 0.5])
      tightmap
      if prntfg=='y';print(['figs/temp_megaD_comp.png'],'-dpng');end
   end


   for j=1:length(ds)
      varX=Cpdsi(:,:,j);
      figure
      hold on
      s=load('coast');mlat=s.lat;mlon=s.long;
      h=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
      %h=axesm('MapProjection','wagner4','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[120 33]);
      varZ=cat(2,varX,varX(:,1));
      pcolorm(lat,[lon; 360],varZ);
      %cptcmap('GMT_red2green')
      %cptcmap('GMT_drywet')
      %C=importdata('MPL_BrBG.txt');
      C=load('./colormaps/MPL_BrBG.mat');
      colormap(C.cmap)
      caxis([-2.3 2.3])
      set(gca,'fontsize',16)
      set(get(colorbar('location','southoutside'),'xlabel'),'string','PDSI','fontsize',16)
      cbh = findobj( 0, 'tag', 'Colorbar' );
      %set(cbh,'YTick',[-1:0.2:1])
      plotm(mlat,mlon,'color',[0.5 0.5 0.5])
      tightmap
      if prntfg=='y';print(['figs/pdsi_megaD_comp.png'],'-dpng');end
   end

end


% %------------------
% % BOXPLOT ANALYSIS
% %------------------
% 
% % Pull out an extreme percentile
% enso_mdrn_pdf=xa_nino_ens(irm2_a:irm2_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(ir2,:);
% 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);
% 
% 
% %ds=ds1;
% ds=svmd;
% cldfrac_d=zeros(size(ds,1),1);
% cldfrac_d_a=zeros(size(ds,1),1);
% cldfrac_do=zeros(size(ds,1),1);
% cldfrac_d_ao=zeros(size(ds,1),1);
% smz=500; % 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_pdfo=enso_pdf+enso_offset;
%   amo_pdfo=amo_pdf+amo_offset;
% 
% 
%   % 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(:));
% 
%   % OFFSET: Fraction of ENSO and AMO states that are below threshold
%   cldfrac_do(i)=length(find(enso_pdfo(:)<em50))/length(enso_pdfo(:));
%   cldfrac_d_ao(i)=length(find(amo_pdfo(:)>am50))/length(amo_pdfo(:));
% 
%   % Randomly sample droughts of this same length from a reference period
%   rsmn1=randi([irm2_a,irm2_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([ir2(1),ir2(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
% hold on
% boxplot(cldfrac_d_r','PlotStyle','compact','symbol','','labelorientation','horizontal','color',[0.1 0.1 0.1])
% plot(cldfrac_d,'o','color',clrs.magenta,'markerfacecolor',clrs.magenta,'markersize',8)
% plot(cldfrac_do,'^','color',clrs.blue,'markerfacecolor',clrs.blue,'markersize',8)
% %title(['Nino3.4 colder than ' num2str(analyrs2(1)) '-' num2str(analyrs2(end)) ' median'])
% xlabel('Drought severity rank','fontsize',18)
% ylabel({'Fraction of Nino3.4 states ',...
% ['colder than ' num2str(analyrs2(1)) '-' num2str(analyrs2(end)) ' median']},'fontsize',18)
% ylim([0.225 0.8])
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 6 5]; % [left bottom width height]
% hold off
% set(gca,'fontsize',18)
% %prntfg='y';
% if prntfg=='y';print(['./figs/enso_coldfrac_drought.png'],'-dpng','-r300');end
% 
% 
% 
% figure
% hold on
% boxplot(cldfrac_d_a_r','PlotStyle','compact','symbol','','labelorientation','horizontal','color',[0.1 0.1 0.1])
% plot(cldfrac_d_a,'o','color',clrs.magenta,'markerfacecolor',clrs.magenta,'markersize',8)
% plot(cldfrac_d_ao,'^','color',clrs.blue,'markerfacecolor',clrs.blue,'markersize',8)
% %title(['AMO warmer than ' num2str(analyrs2(1)) '-' num2str(analyrs2(end)) ' median'])
% xlabel('Drought severity rank','fontsize',18)
% %ylabel('Fraction of states in PDF','fontsize',18)
% ylabel({'Fraction of AMO states',...
% ['warmer than ' num2str(analyrs2(1)) '-' num2str(analyrs2(end)) ' median']},'fontsize',18)
% ylim([0.2 0.925])
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 6 5]; % [left bottom width height]
% hold off
% set(gca,'fontsize',18)
% 
% if prntfg=='y';print(['./figs/amo_warmfrac_drought.png'],'-dpng','-r300');end
% 
% prntfg='n';




% %============================================
% % VIOLIN PLOTS OF ENSO AND AMO DISTRIBUTIONS
% %============================================
% 
% 
% addpath('/home/nsteiger/Dropbox/code/Violinplot')
% 
% ep=struct([]);
% ap=struct([]);
% epr=struct([]);
% apr=struct([]);
% smz=10; % sample size
% 
% 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,:)-nino_off; 
%   ep{i}=enso_pdf(:);
%   amo_pdf=xa_amo_ens(irab,:)-amo_off; 
%   ap{i}=amo_pdf(:);
% 
%   % Randomly sample droughts of this same length from a reference period
%   rsmn1=randi([irm2_a,irm2_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([ir2(1),ir2(end)-length(irab)+1],smz,1);
%   rsmn2=bsxfun(@plus,rsmn11,repmat(0:(length(irab)-1),smz,1));
%   re=[];ra=[];
%   for j=1:smz
%      renso_pdf=xa_nino_ens(rsmn(j,:),:)-nino_off;
%      ramo_pdf=xa_amo_ens(rsmn2(j,:),:)-amo_off;
%      re=cat(1,re,renso_pdf(:));
%      ra=cat(1,ra,ramo_pdf(:));
% 
%   end
%  
%   epr{i}=randsample(re,length(enso_pdf(:)));
%   apr{i}=randsample(ra,length(amo_pdf(:)));
% 
% end
% 
% 
% figure
% violinplot(ep,'histOri','left','widthDiv',[2 1],'showMM',0,'color',clrs.red);
% violinplot(epr,'histOri','right','widthDiv',[2 2],'showMM',0,'color',clrs.grey);
% xlabel('Drought severity rank','fontsize',18)
% ylabel('Nino3.4 (anom. deg. C)','fontsize',18)
% set(gca,'fontsize',16)
% box on
% 
% %print(['./figs/violin_nino34_megaD.png'],'-dpng','-r300');
% 
% figure
% violinplot(ap,'histOri','left','widthDiv',[2 1],'showMM',0,'color',clrs.blue);
% violinplot(apr,'histOri','right','widthDiv',[2 2],'showMM',0,'color',clrs.grey);
% xlabel('Drought severity rank','fontsize',18)
% ylabel('AMO (anom. deg. C)','fontsize',18)
% set(gca,'fontsize',16)
% box on
% 
% %print(['./figs/violin_amo_megaD.png'],'-dpng','-r300');



% %===============================
% % PLOT PDFS OF SELECTED PERIODS
% %===============================
% 
% %----------
% % ENSO PDF
% %----------
% 
% % Get modern pdf and mean offset value
% enso_mdrn_pdf=xa_nino_ens(irm2_a:irm2_b,:);
% mnoff=mean(enso_mdrn_pdf(:));
% 
% ds=ds1;
% %ds=ds2;
% ird_all=[];
% figure; hold on;
% 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;
% 
%   % Collect all months of drought
%   ird_all=cat(2,ird_all,irm_a:irm_b);
% 
%   [~,~,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,:)-mnoff; 
% 
%   x=(min(enso_pdf(:))-1):0.01:(max(enso_pdf(:))+1);
%   pk1=fitdist(enso_pdf(:),'kernel');
%   Yk1=pdf(pk1,x); 
%   plot(x,Yk1,'color',[0.5 0.5 0.5],'linewidth',2)
%   title('Nino3.4')
%   xlabel('Temp. anom. (deg C)','fontsize',18)
%   ylabel('Probability density','fontsize',18)
%   set(gca,'fontsize',18)
%   box on
% 
%   %x2=(min(amo_pdf(:))-1):0.01:(max(amo_pdf(:))+1);
%   %pk2=fitdist(amo_pdf(:),'kernel');
%   %Yk2=pdf(pk2,x2); 
%   %plot(x2,Yk2)
%   %title('AMO PDF')
%   %xlabel('Temperature (deg C)','fontsize',18)
%   %ylabel('Probability density','fontsize',18)
%   %set(gca,'fontsize',18)
% 
% end
% 
% % PDF of the contemporary period
% enso_mdrn_pdf=enso_mdrn_pdf(:)-mnoff;
% x=(min(enso_mdrn_pdf(:))-1):0.01:(max(enso_mdrn_pdf(:))+1);
% pk1=fitdist(enso_mdrn_pdf(:),'kernel');
% Yk1=pdf(pk1,x); 
% plot(x,Yk1,'k','linewidth',3)
% 
% % PDF of the 800-1601 period
% %enso1_pdf=xa_nino_ens(irm1_a:irm1_b,:);
% %x=(min(enso1_pdf(:))-1):0.01:(max(enso1_pdf(:))+1);
% %pk1=fitdist(enso1_pdf(:),'kernel');
% %Yk1=pdf(pk1,x); 
% %plot(x,Yk1,'k','linewidth',3)
% 
% 
% % PDF of the non-drought years
% irw = setdiff(irm1_a:irm1_b,ird_all);
% ensow_pdf=xa_nino_ens(irw,:)-mnoff;
% x_w=(min(ensow_pdf(:))-1):0.01:(max(ensow_pdf(:))+1);
% pk_w=fitdist(ensow_pdf(:),'kernel');
% Yk_w=pdf(pk_w,x_w); 
% plot(x_w,Yk_w,'b','linewidth',3)
% 
% xlim([-5 5])
% 
% if prntfg=='y';print(['./figs/enso_droughtpdf_shifts.png'],'-dpng','-r300');end
% 
% 
% % Plot both wet and drought PDFs
% %ensod_pdf=xa_nino_ens(ird_all,:);
% %x_d=(min(ensod_pdf(:))-1):0.01:(max(ensod_pdf(:))+1);
% %pk_d=fitdist(ensod_pdf(:),'kernel');
% %Yk_d=pdf(pk_d,x_d); 
% 
% %figure
% %hold on
% %plot(x_w,Yk_w,'b','linewidth',3)
% %plot(x_d,Yk_d,'k','linewidth',3)
% 
% 
% %----------
% % AMO PDF
% %----------
% 
% amo_mdrn_pdf=xa_amo_ens(ir2,:);
% mnoff=mean(amo_mdrn_pdf(:));
% 
% ds=ds1;
% %ds=ds2;
% ird_all=[];
% figure; hold on;
% 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); % indices of recon segment
% 
%   % Collect all years of drought
%   ird_all=cat(1,ird_all,irab);
%   
%   % Pull out distributions during this drought
%   amo_pdf=xa_amo_ens(irab,:)-mnoff; 
% 
%   x2=(min(amo_pdf(:))-1):0.01:(max(amo_pdf(:))+1);
%   pk2=fitdist(amo_pdf(:),'kernel');
%   Yk2=pdf(pk2,x2); 
%   plot(x2,Yk2,'color',[0.5 0.5 0.5],'linewidth',2)
%   title('AMO')
%   xlabel('Temp. anom. (deg C)','fontsize',18)
%   ylabel('Probability density','fontsize',18)
%   set(gca,'fontsize',18)
%   box on
% 
% end
% 
% % PDF of the contemporary period
% amo_mdrn_pdf=amo_mdrn_pdf-mnoff;
% x=(min(amo_mdrn_pdf(:))-1):0.01:(max(amo_mdrn_pdf(:))+1);
% pk1=fitdist(amo_mdrn_pdf(:),'kernel');
% Yk1=pdf(pk1,x); 
% plot(x,Yk1,'k','linewidth',3)
% 
% % PDF of the 800-1601 period
% %amo1_pdf=xa_amo_ens(ir1,:);
% %x=(min(amo1_pdf(:))-1):0.01:(max(amo1_pdf(:))+1);
% %pk1=fitdist(amo1_pdf(:),'kernel');
% %Yk1=pdf(pk1,x); 
% %plot(x,Yk1,'k','linewidth',3)
% 
% 
% % PDF of the non-drought years
% irw = setdiff(ir1,ird_all);
% amow_pdf=xa_amo_ens(irw,:)-mnoff;
% x_w=(min(amow_pdf(:))-1):0.01:(max(amow_pdf(:))+1);
% pk_w=fitdist(amow_pdf(:),'kernel');
% Yk_w=pdf(pk_w,x_w); 
% plot(x_w,Yk_w,'b','linewidth',3)
% 
% xlim([-1 1])
% 
% if prntfg=='y';print(['./figs/amo_droughtpdf_shifts.png'],'-dpng','-r300');end
% 




% %======================================
% %     PLOT HYDROCLIMATE INDICES
% %======================================
% 
% 
% rp=nancorr(Xa_di(ir3),DA_di(id3));
% 
% figure
% hold on
% p1=plot(id3,DA_di(id3),'linewidth',2);
% p2=plot(ir3,Xa_di(ir3),'linewidth',2);
% xlim([ir3(1) ir3(end)])
% %plot(ir3,zeros(length(ir3),1),'--k')
% text(0.05,0.05,['r = ' sprintf('%0.2f',rp)],'units','normalized','fontsize',18)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 8.5 4.5]; % [left bottom width height]
% legend('NADA','PHYDA')
% box on
% xlabel('Years','fontsize',18)
% ylabel('NASW PDSI','fontsize',18)
% set(gca,'fontsize',16)
% hold off
% 
% if prntfg=='y';print(['figs/naswpdsi_phyda_nada.png'],'-dpng','-r300'); end












% 
% % Do next plot
% figure
% hold on
% plot(mon_r(irm2_a:irm2_b),zscore(xa_nino(irm2_a:irm2_b)))
% plot(ir2,zscore(Xa_di(ir2)))
% legend('Nino','PDSI')
% y=zscore(Xa_di(ir2));
% yMin=min(y)-0.5;yMax=max(y)+0.5;
% 
% % ADD NEXT IN THE PLOT VERTICAL BARS TO SEE PERIODS OF DROUGHT
% % Add a patches for top five times from truth
% for i=1:length(ds)
%     p=patch(ir2(1)+[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
% xlim([ir2(1) ir2(end)])
% 
% if prntfg=='y';print(['figs/nino_pdsi_megaD_highlight_prd2.png'],'-dpng');end






% %========================
% % CROSS CORRELATIONS
% %========================
% 
% %[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)')
% 
% 
% %----------
% % NINO 3.4
% %----------
% 
% % Compare the xcorr for the two different time periods
% a_nino=mon2ann(xa_nino,4,3);
% [r3,l3]=xcorr(zscore(Xa_di(ir1)),zscore(a_nino(ir1)),10,'unbiased');
% [r4,l4]=xcorr(zscore(Xa_di(ir2)),zscore(a_nino(ir2)),10,'unbiased');
% 
% 
% % Cross correlation using prior data
% 
% % PDSI
% k=find(state_tp=='q');
% Xvar=Xb_o(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);
% Xbpdsi=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),length(p_yrs));
% Xbpdsi(xmeta{k}.lndidx,:)=Xvar;
% Xbpdsi=reshape(Xbpdsi,length(xmeta{k}.lat),length(xmeta{k}.lon),length(p_yrs));
% Xv_msk=bsxfun(@times,Xbpdsi,nasw);
% Xb_di=wmean_a(Xv_msk,A); % RECONSTRUCTION DROUGHT INDEX
% 
% ab_nino=mon2ann(xb_nino,4,3);
% [r1,l1]=xcorr(zscore(Xb_di(1:end-1)),zscore(ab_nino),10,'unbiased');
% 
% figure;
% hold on;
% plot(l1,r1,'-o')
% plot(l3,r3,'-o')
% plot(l4,r4,'-o')
% xlabel('Lag (yrs)')
% title('xcorr(PDSI,Nino 3.4)')
% legend('Prior',[num2str(analyrs1(1)) '-' num2str(analyrs1(end))],[num2str(analyrs2(1)) '-' num2str(analyrs2(end))])
% 
% if prntfg=='y';print(['figs/xcorr_nino_pdsi.png'],'-dpng');end
% 
% 
% %----------
% % ITCZ
% %----------
% 
% [r5,l5]=xcorr(zscore(Xa_di(ir1)),zscore(xa_itcz_pac(ir1)),10,'unbiased');
% [r6,l6]=xcorr(zscore(Xa_di(ir2)),zscore(xa_itcz_pac(ir2)),10,'unbiased');
% [r7,l7]=xcorr(zscore(Xb_di),zscore(xb_itcz_pac),10,'unbiased');
% 
% figure;
% hold on;
% plot(l7,r7,'-o')
% plot(l5,r5,'-o')
% plot(l6,r6,'-o')
% xlabel('Lag (yrs)')
% title('xcorr(PDSI, ITCZ Pac)')
% legend('Prior',[num2str(analyrs1(1)) '-' num2str(analyrs1(end))],[num2str(analyrs2(1)) '-' num2str(analyrs2(end))])
% 
% if prntfg=='y';print(['figs/xcorr_itczPac_pdsi.png'],'-dpng');end
% 
% [r8,l8]=xcorr(zscore(Xa_di(ir1)),zscore(xa_itcz_sa(ir1)),10,'unbiased');
% [r9,l9]=xcorr(zscore(Xa_di(ir2)),zscore(xa_itcz_sa(ir2)),10,'unbiased');
% [r10,l10]=xcorr(zscore(Xb_di),zscore(xb_itcz_sa),10,'unbiased');
% 
% figure;
% hold on;
% plot(l10,r10,'-o')
% plot(l8,r8,'-o')
% plot(l9,r9,'-o')
% xlabel('Lag (yrs)')
% title('xcorr(PDSI, ITCZ SA)')
% legend('Prior',[num2str(analyrs1(1)) '-' num2str(analyrs1(end))],[num2str(analyrs2(1)) '-' num2str(analyrs2(end))])
% 
% if prntfg=='y';print(['figs/xcorr_itczSA_pdsi.png'],'-dpng');end
% 
% 
% [r11,l11]=xcorr(zscore(Xb_di),zscore(xb_itcz_atl),10,'unbiased');
% [r12,l12]=xcorr(zscore(Xa_di(ir1)),zscore(xa_itcz_atl(ir1)),10,'unbiased');
% [r13,l13]=xcorr(zscore(Xa_di(ir2)),zscore(xa_itcz_atl(ir2)),10,'unbiased');
% 
% figure;
% hold on;
% plot(l11,r11,'-o')
% plot(l12,r12,'-o')
% plot(l13,r13,'-o')
% xlabel('Lag (yrs)')
% title('xcorr(PDSI, ITCZ Atl)')
% legend('Prior',[num2str(analyrs1(1)) '-' num2str(analyrs1(end))],[num2str(analyrs2(1)) '-' num2str(analyrs2(end))])
% 
% if prntfg=='y';print(['figs/xcorr_itczAtl_pdsi.png'],'-dpng');end





% %=================================
% % Composite over the top droughts
% %=================================
% plot_ind='n';
% cpyrs=[];
% for i=1:length(ds);
%    c0=(ds(i,1):ds(i,2))+analyrs1(1)-1;
%    if plot_ind=='y'
%    % Plot composites over individual droughts
%    
%    % Remove referene period, then composite
%    T0=bsxfun(@minus,Xt2m,mean(Xt2m(:,:,ir2),3));
%    varX=mean(T0(:,:,c0),3);
% 
%    % Plot the composites
%    figure
%    hold on
%    s=load('coast');mlat=s.lat;mlon=s.long;
%    %h=axesm('MapProjection','wagner4','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
%    h=axesm('MapProjection','wagner4','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[120 33]);
%    varZ=cat(2,varX,varX(:,1));
%    pcolorm(lat,[lon; 360],varZ);
%    %C=importdata('redblue_256_rgb.txt');
%    %colormap(flip(C,1)/255)
%    C=load('./colormaps/NCV_blue_red.mat');
%    colormap(C.cmap)
%    caxis([-1 1])
%    set(gca,'fontsize',16)
%    set(get(colorbar('location','southoutside'),'xlabel'),'string','Temperature anom. (deg C)','fontsize',16)
%    cbh = findobj( 0, 'tag', 'Colorbar' );
%    %set(cbh,'YTick',[-1:0.2:1])
%    plotm(mlat,mlon,'color',[0.5 0.5 0.5])
%    tightmap
%    if prntfg=='y';print(['figs/temp_megaD_comp_d' num2str(i) '.png'],'-dpng');end
% 
%    end
%    cpyrs=cat(1,cpyrs,c0(:));
% end
% 
% % Get composite relative to mean reference period
% %Ct2m=mean(Xt2m(:,:,cpyrs),3)-mean(Xt2m(:,:,ir2),3);
% %Cpdsi=mean(Xpdsi(:,:,cpyrs),3)-mean(Xpdsi(:,:,ir2),3);
% 
% % Remove referene period, then composite
% T0=bsxfun(@minus,Xt2m,mean(Xt2m(:,:,ir2),3));
% Ct2m=mean(T0(:,:,cpyrs),3);
% P0=bsxfun(@minus,Xpdsi,mean(Xpdsi(:,:,ir2),3));
% Cpdsi=mean(P0(:,:,cpyrs),3);
% 
% 
% % Plot the composites
% varX=Ct2m;
% figure
% hold on
% s=load('coast');mlat=s.lat;mlon=s.long;
% %h=axesm('MapProjection','wagner4','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
% h=axesm('MapProjection','wagner4','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[120 33]);
% varZ=cat(2,varX,varX(:,1));
% pcolorm(lat,[lon; 360],varZ);
% %C=importdata('redblue_256_rgb.txt');
% %colormap(flip(C,1)/255)
% C=load('./colormaps/NCV_blue_red.mat');
% colormap(C.cmap)
% caxis([-1 1])
% set(gca,'fontsize',16)
% set(get(colorbar('location','southoutside'),'xlabel'),'string','Temperature anom. (deg C)','fontsize',16)
% cbh = findobj( 0, 'tag', 'Colorbar' );
% %set(cbh,'YTick',[-1:0.2:1])
% plotm(mlat,mlon,'color',[0.5 0.5 0.5])
% tightmap
% if prntfg=='y';print(['figs/temp_megaD_comp.png'],'-dpng');end
% 
% 
% 
% addpath('/home/nsteiger/MATLAB/matlab-goody-bag/cptcmap-pkg-master/cptcmap')
% varX=Cpdsi;
% figure
% hold on
% s=load('coast');mlat=s.lat;mlon=s.long;
% %h=axesm('MapProjection','wagner4','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
% h=axesm('MapProjection','wagner4','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[120 33]);
% varZ=cat(2,varX,varX(:,1));
% pcolorm(lat,[lon; 360],varZ);
% %cptcmap('GMT_red2green')
% %cptcmap('GMT_drywet')
% %C=importdata('MPL_BrBG.txt');
% C=load('./colormaps/MPL_BrBG.mat');
% colormap(C.cmap)
% caxis([-2.3 2.3])
% set(gca,'fontsize',16)
% set(get(colorbar('location','southoutside'),'xlabel'),'string','PDSI','fontsize',16)
% cbh = findobj( 0, 'tag', 'Colorbar' );
% %set(cbh,'YTick',[-1:0.2:1])
% plotm(mlat,mlon,'color',[0.5 0.5 0.5])
% tightmap
% if prntfg=='y';print(['figs/pdsi_megaD_comp.png'],'-dpng');end
% 








































% 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')
% 


