% DROUGHT ANALYSIS FIGURES

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


%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';
%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=800:1600;
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); % 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;

%=====================================
% LOAD ENSEMBLE MEGADROUGHT ESTIMATES
%=====================================
load('svmd_ens_90pagree.mat')


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


%===========================================================================
% LOAD FORCED CLIMATE RESPONSE ESTIMATES

% Full forcing
load('../megad-forcing/output_twoboxCM_full_06-Aug-2018_15:09:13.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_f=gmt_est_mn-mean(gmt_est_mn(ig2));

% solar 
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')
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')
gmt_est_mn=mean(gmt_fr,2);
gmt_est_g=gmt_est_mn-mean(gmt_est_mn(ig2));


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

load distinct_colors_20.mat

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

% 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],'state_tp','xmeta','Xa_m','Xa_ens','reconYrs','p_yrs')


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

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


% Load DJF for Nino analysis
%efl='cesm_lme010_r12000_p998_state_2ozn_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:16:10.mat';
%efl='cesm_lme010_r12000_p998_state_2ozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:24:35.mat';
%load([pth,efl],'state_tp','xmeta','Xa_ens','reconYrs','p_yrs')

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);


% Precentiles to check
%pct=[1 5 10 25 50 75 90 95 99];
pct=[5 10 25 50 75 90 95];


%==========
% FORCINGS
%==========

% Reference PDF
gmtf_ref_pdf=gmt_est_f(ig2);


%ds=ds1;
ds=svmd;
ird_all=[];

% Compute percentiles of reference
gmtf_ref_p=prctile(gmtf_ref_pdf(:),pct);

gmtf_d_p=zeros(size(ds,1),length(pct));

for i=1:size(ds,1)

  [~,~,irab] = intersect((ds(i,1)+analyrs1(1)-1):(ds(i,2)+analyrs1(1)-1),myrs); % indices of recon segment

  % Collect all years of drought
  ird_all=cat(1,ird_all,irab);
  
  % Pull out distributions during this drought
  gmtf_pdf=gmt_est_f(irab,:); 
  %amo_pdf=xa_amo_ens(irab,:)-mean(mean(xa_amo_ens(irab,:),2)); % adjust mean state to reference 

  gmtf_d_p(i,:)=prctile(gmtf_pdf(:),pct);

end

gmtf_d_pch=bsxfun(@minus,gmtf_d_p,gmtf_ref_p);


% Colormaps for lines
load('./colormaps/MPL_Reds.mat')
%load('./colormaps/MPL_viridis.mat')
cidx=flip(round(linspace(10,length(cmap),size(ds,1))));

figure 
hold on
for j=1:size(ds,1)
   plot(pct,gmtf_d_pch(j,:),'-o','linewidth',1.5,'color',cmap(cidx(j),:));
end
box on
plot(pct,zeros(length(pct),1),'--','linewidth',2,'color',clrs.grey)
set(gca,'fontsize',16)
title('Forcing')
ax=gca;
xlabel('Percentiles','fontsize',18)
ylabel(['Temperature change (' char(176) 'C)'],'fontsize',18)
ax.XTick=pct;
ax.XTickLabel=num2cell(pct);
xtickangle(90)

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


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

% Get an offset value for anomalies
amo_off=mean(mean(xa_amo_ens(ir2,:),2));
% Reference PDF
amo_ref_pdf=xa_amo_ens(ir2,:)-amo_off;


%ds=ds1;
ds=svmd;
ird_all=[];

% Compute percentiles of reference
amo_ref_p=prctile(amo_ref_pdf(:),pct);

amo_d_p=zeros(size(ds,1),length(pct));

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,:)-amo_off; 
  %amo_pdf=xa_amo_ens(irab,:)-mean(mean(xa_amo_ens(irab,:),2)); % adjust mean state to reference 

  amo_d_p(i,:)=prctile(amo_pdf(:),pct);

end

amo_d_pch=bsxfun(@minus,amo_d_p,amo_ref_p);


% Colormaps for lines
%load('./colormaps/MPL_Reds.mat')
%cidx=flip(round(linspace(10,128,size(ds,1))));

figure 
hold on
for j=1:size(ds,1)
   plot(pct,amo_d_pch(j,:),'-o','linewidth',1.5,'color',cmap(cidx(j),:));
end
box on
plot(pct,zeros(length(pct),1),'--','linewidth',2,'color',clrs.grey)
set(gca,'fontsize',16)
title('AMO')
ax=gca;
xlabel('Percentiles','fontsize',18)
ylabel(['Temperature change (' char(176) 'C)'],'fontsize',18)
ax.XTick=pct;
ax.XTickLabel=num2cell(pct);
xtickangle(90)

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


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


% Get an offset value for anomalies
nino_off=mean(mean(xa_nino_ens(irm2_a:irm2_b,:),2));

% Reference PDF
nino_ref_pdf=xa_nino_ens(irm2_a:irm2_b,:)-nino_off;


%ds=ds1;
ds=svmd;
ird_all=[];
ird_ann=[];

% Compute percentiles of reference
nino_ref_p=prctile(nino_ref_pdf(:),pct);

nino_d_p=zeros(size(ds,1),length(pct));

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);
  
  % Pull out distributions during this drought
  nino_pdf=xa_nino_ens(irm_a:irm_b,:)-nino_off; 
  %nino_pdf=xa_nino_ens(irm_a:irm_b,:)-mean(mean(xa_nino_ens(irm_a:irm_b,:),2)); % adjust mean state to reference 

  nino_d_p(i,:)=prctile(nino_pdf(:),pct);

end

nino_d_pch=bsxfun(@minus,nino_d_p,nino_ref_p);


% Colormaps for lines
%load('./colormaps/MPL_Reds.mat')
%cidx=flip(round(linspace(10,128,size(ds,1))));

figure 
hold on
for j=1:size(ds,1)
   plot(pct,nino_d_pch(j,:),'-o','linewidth',1.5,'color',cmap(cidx(j),:));
end
box on
plot(pct,zeros(length(pct),1),'--','linewidth',2,'color',clrs.grey)
set(gca,'fontsize',16)
title('Nino3.4')
ax=gca;
xlabel('Percentiles','fontsize',18)
ylabel(['Temperature change (' char(176) 'C)'],'fontsize',18)
ax.XTick=pct;
ax.XTickLabel=num2cell(pct);
xtickangle(90)

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

% 
% 
% %----------
% % ENSO PDF
% %----------
% 
% % Get modern pdf and mean offset value
% enso_mdrn_pdf=xa_nino_ens(irm2_a:irm2_b,:);
% enp=struct([]);
% 
% %ds=ds1;
% ird_all=[];
% ird_ann=[];
% figure; hold on;
% 
% % PDF of the contemporary period
% enso_mdrn_pdf=enso_mdrn_pdf(:)-nino_off;
% 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,'color',clrs.blue,'linewidth',2.5)
% area(x,Yk1,'facecolor',clrs.blue,'linestyle','none')
% 
% 
% 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);
%   
%   % Pull out distributions during this drought
%   enso_pdf=xa_nino_ens(irm_a:irm_b,:)-nino_off; 
%   enp{i}.pdf=enso_pdf(:);
% 
%   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',1.75)
%   plot(x,Yk1,'color',clrs.red,'linewidth',1.75)
%   title('Nino3.4')
%   %xlabel('Temp. anom. (deg C)','fontsize',18)
%   %ylabel('Probability density','fontsize',18)
%   %set(gca,'fontsize',18)
%   box on
% 
% 
% 
% end
% 
% 
% % PDF of the non-drought years
% irw = setdiff(irm1_a:irm1_b,ird_all);
% ensow_pdf=xa_nino_ens(irw,:)-nino_off;
% 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,'k','linewidth',3)
% 
% 
% % PDF of all megadroughts
% xa_nino_ann=mon2ann(xa_nino_ens',4,3)';
% %xa_nino_ann=mon2ann(xa_nino_ens',12,2)';% DJF nino
% xa_nino_ann(end+1,:)=xa_nino_ann(end,:); % repeat end year because of annualizing
% nino_mdrght_pdf=xa_nino_ann(ird_ann,:)-nino_off;
% 
% 
% 
% % Remove the labels, squash down the figure
% xlim([-4.75 4.75])
% set(gca,'xticklabel',[])
% set(gca,'yticklabel',[])
% set(gca,'XDir','reverse'); % Reverse direction for flipping on its side
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 4.5 2.5]; % [left bottom width height]
% 
% if prntfg=='y';print(['./figs/enso_droughtpdf_shifts.png'],'-dpng','-r300');end
% %print(['./figs/enso_droughtpdf_shifts.png'],'-dpng','-r300');
% 
% 
% % Get annualized versions of the other periods as well
% nino_mdrn_pdf=xa_nino_ann(ir2,:)-nino_off;
% irwa = setdiff(ir1,ird_ann);
% nino_wt_pdf=xa_nino_ann(irwa,:)-nino_off;
% 
% 
% 
% 
% 
% %----------------------
% % FILTERED TIME SERIES
% %----------------------
% 
% 
% 
% figure
% hold on
% %plot(analyrs3,amo_p(ir3,2)-amo_off,'linewidth',2,'color',[0 0.4470 0.7410])
% h2=plot(analyrs3,amo_p(ir3,1)-amo_off,'--','color',[0.5 0.5 0.5]);
% plot(analyrs3,amo_p(ir3,3)-amo_off,'--','color',[0.5 0.5 0.5])
% h1=plot(analyrs3,amo_p(ir3,2)-amo_off,'linewidth',2,'color',[0 0 0]);
% h1=plot(analyrs3,smooth(amo_p(ir3,2)-amo_off,20,'lowess'),'linewidth',2,'color',clrs.yellow);
% %legend([h1,h2],{'AMO (NASST)','5-95th prcntl'},'location','southwest')
% set(gca,'fontsize',16)
% ylim([-0.8 1])
% xlim([analyrs3(1) analyrs3(end)])
% xlabel('Year','fontsize',18)
% ylabel(['Temperature anomaly (' char(176) 'C)'],'fontsize',18)
% box on
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]
% 
% if  prntfg=='y'
%  %print('-dpdf',['./figs/timeseries_amo_recon_analyrs_' num2str(analyrs3(1)) '_' num2str(analyrs3(end)) '.pdf'])
%  print('-dpng',['./figs/timeseries_filter_amo_recon_analyrs_' num2str(analyrs3(1)) '_' num2str(analyrs3(end)) '.png'],'-r300')
% end
% 
% 
% 
% figure
% hold on
% h2=plot(mon_r(irm3_a:irm3_b),nino_p(irm3_a:irm3_b,1)-nino_off,'--','color',[0.5 0.5 0.5]);
% plot(mon_r(irm3_a:irm3_b),nino_p(irm3_a:irm3_b,3)-nino_off,'--','color',[0.5 0.5 0.5])
% h1=plot(mon_r(irm3_a:irm3_b),nino_p(irm3_a:irm3_b,2)-nino_off,'linewidth',2,'color',[0 0 0]);
% h1=plot(mon_r(irm3_a:irm3_b),smooth(nino_p(irm3_a:irm3_b,2)-nino_off,12*20,'lowess'),'linewidth',2,'color',clrs.yellow);
% %legend([h1,h2],{'Nino3.4 (mon)','5-95th prcntl'},'location','northeast')
% set(gca,'fontsize',16)
% ylim([-4.75 4.75])
% xlim([analyrs3(1) analyrs3(end)])
% xlabel('Year','fontsize',18)
% ylabel(['Temperature anomaly (' char(176) 'C)'],'fontsize',18)
% box on
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]
% 
% 
% 
% if  prntfg=='y'
%  %print('-dpdf',['./figs/timeseries_nino_recon_analyrs_' num2str(analyrs3(1)) '_' num2str(analyrs3(end)) '.pdf'])
%  print('-dpng',['./figs/timeseries_filter_nino_recon_analyrs_' num2str(analyrs3(1)) '_' num2str(analyrs3(end)) '.png'],'-r300')
% end
% 
% 
% 
% 
% 
% %================================
% %   JOINT PDF OF ENSO AND AMO
% %================================
% 
% 
% X1=cat(2,nino_mdrn_pdf(:),amo_mdrn_pdf(:));
% %X2=cat(2,nino_mdrght_pdf(:),amo_mdrght_pdf(:));
% 
% %bwn=[0.0422,0.0097];
% %[f1,x1] = ksdensity(X1,'bandwidth',bwn);
% [f1,x1] = ksdensity(X1);
% %[f2,x2] = ksdensity(X2);
% 
% xl=[-4,4];yl=[-0.9,0.9];
% %[xq,yq] = meshgrid(linspace(-4.75,4.75),linspace(-1,1));
% [xq,yq] = meshgrid(linspace(xl(1),xl(2)),linspace(yl(1),yl(2)));
% z1 = griddata(x1(:,1),x1(:,2),f1,xq,yq);
% %z2 = griddata(x2(:,1),x2(:,2),f2,xq,yq);
% 
% figure
% hold on
% %[c1,h1]=contour(xq,yq,z1,'fill','on');
% %[c1,h1]=contour(xq,yq,z1);
% %[c1,h1]=contourf(xq,yq,z1,0.05:0.15:1.4);
% [c1,h1]=contourf(xq,yq,z1,0.1:0.2:1.4,'showtext','on');
% colormap(flip(gray))
% clb=colorbar;
% clb.Label.String='Probability density';
% clb.Label.FontSize=16;
% %clb.Label.HorizontalAlignment='right';
% %[c2,h2]=contour(xq,yq,z2,'showtext','on','linecolor',clrs.red,'linestyle','-.');
% %[c2,h2]=contour(xq,yq,z2,'linecolor',clrs.red,'linestyle','-.');
% %clabel(c2,h2,'Color',clrs.red)
% 
% % Get a colormap for the errorbars
% load('./colormaps/MPL_Reds.mat')
% cidx=flip(round(linspace(6,128,size(ds,1))));
% 
% % Plot each drought pdf 
% 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; 
%   enso_pdf=xa_nino_ann(irab,:)-nino_off;
%   ep=enso_pdf(:);
%   amo_pdf=xa_amo_ens(irab,:)-amo_off; 
%   ap=amo_pdf(:);
% 
%   %errorbar(x,y,yneg,ypos,xneg,xpos)
%   %errorbar(mean(ep),mean(ap),prctile(ap,5),prctile(ap,95),prctile(ep,5),prctile(ep,95),...
%   % 'o','markeredgecolor',clrs.red,'color',clrs.red,'linewidth',2)
%   b=errorbar(median(ep),median(ap),prctile(ap,5),prctile(ap,95),prctile(ep,5),prctile(ep,95),...
%    'o','linewidth',2,'color',cmap(cidx(i),:));
% 
% end
% 
% 
% line([xl(1),xl(2)],[0,0],'color','k','linestyle','--')
% line([0,0],[yl(1),yl(2)],'color','k','linestyle','--')
% %set(gca,'tickdir','both')
% %xlim([-3.75 3.75])
% %ylim([-0.7 0.7])
% xlim([-3 3])
% ylim([-0.7 0.7])
% set(gca,'fontsize',16)
% xlabel(['Nino3.4 index (annual, ' char(176) 'C)'],'fontsize',18)
% ylabel(['AMO (NASST) index (' char(176) 'C)'],'fontsize',18)
% box on
% 
% %prntfg='y';
% if prntfg=='y';print(['./figs/jointpdf_megaD_shifts_errbar.png'],'-dpng','-r300'); end
% %print(['./figs/jointpdf_megaD_shifts_clrs.png'],'-dpng','-r300'); 
% %print(['./figs/jointpdf_megaD_shifts_pdf.png'],'-dpng','-r300'); 
% 
% %prntfg='n';
% 
% 
% 
% %=====================================================
% %   JOINT PDF OF ENSO AND AMO WITH MEAN STATE SHIFT
% %=====================================================
% 
% 
% figure
% hold on
% %[c1,h1]=contour(xq,yq,z1,'fill','on');
% %[c1,h1]=contour(xq,yq,z1);
% %[c1,h1]=contourf(xq,yq,z1,0.05:0.15:1.4);
% [c1,h1]=contourf(xq,yq,z1,0.1:0.2:1.4,'showtext','on');
% colormap(flip(gray))
% clb=colorbar;
% clb.Label.String='Probability density';
% clb.Label.FontSize=16;
% %clb.Label.HorizontalAlignment='right';
% %[c2,h2]=contour(xq,yq,z2,'showtext','on','linecolor',clrs.red,'linestyle','-.');
% %[c2,h2]=contour(xq,yq,z2,'linecolor',clrs.red,'linestyle','-.');
% %clabel(c2,h2,'Color',clrs.red)
% 
% % Get a colormap for the errorbars
% load('./colormaps/MPL_Blues.mat')
% cidx=flip(round(linspace(6,128,size(ds,1))));
% 
% % Plot each drought pdf 
% 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; 
%   enso_pdf=xa_nino_ann(irab,:)-nino_off;
%   edf=xa_nino_ann(ir1,:)-nino_off;
%   ep=enso_pdf(:)-mean(edf(:));
% 
%   amo_pdf=xa_amo_ens(irab,:)-amo_off; 
%   adf=xa_amo_ens(ir1,:)-amo_off; 
%   ap=amo_pdf(:)-mean(adf(:));
% 
%   %errorbar(x,y,yneg,ypos,xneg,xpos)
%   %errorbar(mean(ep),mean(ap),prctile(ap,5),prctile(ap,95),prctile(ep,5),prctile(ep,95),...
%   % 'o','markeredgecolor',clrs.red,'color',clrs.red,'linewidth',2)
%   b=errorbar(median(ep),median(ap),prctile(ap,5),prctile(ap,95),prctile(ep,5),prctile(ep,95),...
%    'o','linewidth',2,'color',cmap(cidx(i),:));
% 
% end
% 
% 
% line([xl(1),xl(2)],[0,0],'color','k','linestyle','--')
% line([0,0],[yl(1),yl(2)],'color','k','linestyle','--')
% %set(gca,'tickdir','both')
% %xlim([-3.75 3.75])
% %ylim([-0.7 0.7])
% xlim([-3 3])
% ylim([-0.7 0.7])
% set(gca,'fontsize',16)
% xlabel(['Nino3.4 index (annual, ' char(176) 'C)'],'fontsize',18)
% ylabel(['AMO (NASST) index (' char(176) 'C)'],'fontsize',18)
% box on
% 
% %prntfg='y';
% if prntfg=='y';print(['./figs/jointpdf_megaD_shifts_errbar_offset.png'],'-dpng','-r300'); end
% %print(['./figs/jointpdf_megaD_shifts_clrs.png'],'-dpng','-r300'); 
% %print(['./figs/jointpdf_megaD_shifts_pdf.png'],'-dpng','-r300'); 
% 
% %prntf='n';
% 
% 
% 
% 
% %================================================
% %   JOINT PDF OF ENSO AND AMO: FOR EACH DROUGHT
% %================================================
% 
% 
% X1=cat(2,nino_mdrn_pdf(:),amo_mdrn_pdf(:));
% 
% %bwn=[0.0422,0.0097];
% %[f1,x1] = ksdensity(X1,'bandwidth',bwn);
% [f1,x1] = ksdensity(X1);
% 
% xl=[-4,4];yl=[-0.9,0.9];
% %[xq,yq] = meshgrid(linspace(-4.75,4.75),linspace(-1,1));
% [xq,yq] = meshgrid(linspace(xl(1),xl(2),200),linspace(yl(1),yl(2),200));
% z1 = griddata(x1(:,1),x1(:,2),f1,xq,yq);
% %z2 = griddata(x2(:,1),x2(:,2),f2,xq,yq);
% 
% %clb.Label.HorizontalAlignment='right';
% %[c2,h2]=contour(xq,yq,z2,'showtext','on','linecolor',clrs.red,'linestyle','-.');
% %[c2,h2]=contour(xq,yq,z2,'linecolor',clrs.red,'linestyle','-.');
% %clabel(c2,h2,'Color',clrs.red)
% 
% 
% % Plot each drought pdf 
% for i=1:size(ds,1)
% 
%    figure
%    hold on
%    [c1,h1]=contourf(xq,yq,z1,0.1:0.2:1.5,'showtext','on');
%    colormap(flip(gray))
%    %colormap(gray)
%    clb=colorbar;
%    clb.Label.String='Probability Density';
%    clb.Label.FontSize=16;
% 
%    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; 
%    enso_pdf=xa_nino_ann(irab,:)-nino_off;
%    amo_pdf=xa_amo_ens(irab,:)-amo_off; 
% 
%    % PLOT EACH MEGADROUGHT PDF INDIVIDUALLY
%    X2=cat(2,enso_pdf(:),amo_pdf(:));
%    [f2,x2] = ksdensity(X2);
%    z2 = griddata(x2(:,1),x2(:,2),f2,xq,yq);
%    
%    [c2,h2]=contour(xq,yq,z2,0.1:0.2:1.5,'showtext','on','linecolor',clrs.red,'linestyle','-','linewidth',1.25);
%    clabel(c2,h2,'Color',clrs.red)
% 
%    % Find the max and plot it
%    [~,m1]=max(max(z2,[],1));[~,m2]=max(max(z2,[],2));
%    plot(xq(m2,m1),yq(m2,m1),'*','linewidth',2,'color',clrs.red,'markersize',10,'markerfacecolor',clrs.red)
% 
%    line([xl(1),xl(2)],[0,0],'color','k','linestyle','--','linewidth',1.5)
%    line([0,0],[yl(1),yl(2)],'color','k','linestyle','--','linewidth',1.5)
%    %set(gca,'tickdir','both')
%    %xlim([-3.75 3.75])
%    xlim([-2.5 2.5])
%    ylim([-0.7 0.7])
%    set(gca,'fontsize',16)
%    xlabel('Nino3.4 (annual temp. anomaly ^{\circ}C)','fontsize',18)
%    ylabel('AMO (temp. anomaly ^{\circ}C)','fontsize',18)
%    box on
% 
%    if prntfg=='y';print(['./figs/jointpdf_megaD' num2str(i) '_shift_labels.png'],'-dpng','-r300'); end
%    %errorbar(x,y,yneg,ypos,xneg,xpos)
%    %errorbar(mean(ep),mean(ap),prctile(ap,5),prctile(ap,95),prctile(ep,5),prctile(ep,95),...
%    % 'o','markeredgecolor',clrs.red,'color',clrs.red,'linewidth',2)
%    %b=errorbar(mean(ep),mean(ap),prctile(ap,5),prctile(ap,95),prctile(ep,5),prctile(ep,95),...
%    %'o','linewidth',2,'color',cmap(cidx(i),:));
% 
% end
% 
% 
% 
% % Statistical tests of the distributions for Fig 1
% 
% for i=1:length(ds)
%    disp(['ENSO stats for megadrought = ' num2str(i)])
%    disp('KS-test')
%    [h,p,stats]=kstest2(enp{i}.pdf,enso_mdrn_pdf)
%    disp('Welchs t-test')
%    [h,p,ci,stats]=ttest2(enp{i}.pdf,enso_mdrn_pdf,'vartype','unequal')
% end
% 

% g_ref=ones(length(amo_mdrn_pdf(:)),1);
% g_MD=2*ones(length(amo_mdrght_pdf(:)),1);
% g_WT=3*ones(length(amow_pdf(:)),1);
% G=cat(1,g_ref,g_MD,g_WT);
% 
% v_amo=cat(1,amo_mdrn_pdf(:),amo_mdrght_pdf(:),amow_pdf(:));
% v_nino=cat(1,nino_mdrn_pdf(:),nino_mdrght_pdf(:),nino_wt_pdf(:));
% 
% 
% % Plot of joint pdf scatter plot
% figure
% scatterhist(v_nino,v_amo,'Group',G,'Kernel','on','location','southeast','direction','out')
% 
% 
% % LOOK AT JUST ENSEMBLE MEAN
% g_ref=cell(size(amo_mdrn_pdf,1),1);
% g_ref(:)={'1601-1925'};
% g_MD=cell(size(amo_mdrght_pdf,1),1);
% g_MD(:)={'Megadrought'};
% g_WT=cell(size(amow_pdf,1),1);
% g_WT(:)={'Non-drought'};
% 
% G=cat(1,g_ref,g_WT,g_MD);
% 
% amo=cat(1,mean(amo_mdrn_pdf,2),mean(amow_pdf,2),mean(amo_mdrght_pdf,2));
% nino=cat(1,mean(nino_mdrn_pdf,2),mean(nino_wt_pdf,2),mean(nino_mdrght_pdf,2));
% 
% 
% % Plot of joint pdf scatter plot
% figure
% scatterhist(nino,amo,'Group',G,'Kernel','on','location','southeast','direction','out','marker','.')




% %======================================================
% % Correlation between Nino3.4 and spatial PDSI and T2m
% %======================================================
% 
% 
% % Get nino index from temperature field
% xatnino=nino(Xt2m,lat,lon,3,'n');
% 
% analyrs5=1870:2000;
% [~,~,ir5] = intersect(analyrs5,r_o:r_f);
% 
% pMap=zeros(length(lat)*length(lon),1);
% 
% tMap=zeros(length(lat)*length(lon),1);
% for i=1:length(lat)*length(lon)
% %   pMap(i)=nancorr(xa_nino_mn_ann(ir5),xp1(i,ir5));
% %   tMap(i)=nancorr(xa_nino_mn_ann(ir5),xt1(i,ir5));
%    pMap(i)=nancorr(xatnino(ir4),xp1(i,ir4));
%    tMap(i)=nancorr(xatnino(ir4),xt1(i,ir4));
% %   pMap(i)=nancorr(xatnino(ird_ann),xp1(i,ird_ann));
% %   tMap(i)=nancorr(xatnino(ird_ann),xt1(i,ird_ann));
% end
% 
% 
% s=load('coast');mlat=s.lat;mlon=s.long;
% 
% figure;
% h=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
% pcolorm(lat,lon,reshape(tMap,length(lat),length(lon)));
% %pcolor(reshape(tMap,length(lat),length(lon)))
% C=load('./colormaps/NCV_blue_red.mat');
% colormap(C.cmap)
% caxis([-1 1])
% colorbar
% plotm(mlat,mlon,'color',[0.5 0.5 0.5])
% title('Corr(Nino3.4, T2m A2M), 800-1925')
% if prntfg=='y';print(['./figs/corrmap_nino_t2m.png'],'-dpng','-r300'); end
% 
% 
% figure;
% h=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
% pcolorm(lat,lon,reshape(pMap,length(lat),length(lon)));
% %pcolor(reshape(pMap,length(lat),length(lon)))
% colormap(C.cmap)
% caxis([-1 1])
% colorbar
% plotm(mlat,mlon,'color',[0.5 0.5 0.5])
% title('Corr(Nino3.4, PDSI JJA), 800-1925')
% if prntfg=='y';print(['./figs/corrmap_nino_pdsi.png'],'-dpng','-r300'); end




% %==================================================================================================r
% %                        REGRESSION MODELS (LINEAR & LOGISTIC)
% %===================================================================================================
% 
% 
% ird_ann=[];
% for i=1:size(ds,1)
%   [~,~,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);
% end
% 
% % Bivariate logistic regression
% X(:,1)=mean(xa_nino_ann(ir4,:)-nino_off,2);
% X(:,2)=mean(xa_amo_ens(ir4,:)-amo_off,2);
% 
% y0=zeros(reconYrs,1);
% y0(ird_ann)=1;
% Y=categorical(y0(ir4));
% [B,dev,stats] = mnrfit(X,Y);
% 
% % Bivariate logistic regression with interaction
% Xi(:,1)=mean(xa_nino_ann(ir4,:)-nino_off,2);
% Xi(:,2)=mean(xa_amo_ens(ir4,:)-amo_off,2);
% Xi(:,3)=Xi(:,1).*Xi(:,2);
% [Bi,devi,statsi] = mnrfit(Xi,Y);
% 
% 
% 
% % Explore the predicted probabilities using the range of ENSO and AMO
% 
% % Define a grid of all possible Nino and AMO values
% %ns=xa_nino_ann(ir4,:)-nino_off;as=xa_amo_ens(ir4,:)-amo_off;
% % Just pick the mean values over which the model is built
% ns=mean(xa_nino_ann(ir4,:)-nino_off,2);as=mean(xa_amo_ens(ir4,:)-amo_off,2);
% [mx,my] = meshgrid(min(ns(:)):0.01:max(ns(:)),min(as(:)):0.01:max(as(:)));
% Xm(:,1)=mx(:);
% Xm(:,2)=my(:);
% 
% % Get probability for random observations
% pihat = mnrval(B,Xm);
% [pihat,dlow,dhi] = mnrval(B,Xm,stats);% get uncertainties
% % NOTE THAT UPPER AND LOWER BOUNDS ARE THE SAME!!
% 
% 
% % Reshape to match the grid
% pM=reshape(pihat(:,2),size(mx));
% pL=reshape(dlow(:,2),size(mx));
% %pU=reshape(dhi(:,2),size(mx));
% 
% % figure
% % [Cc,h]=contour(mx,my,pM,'ShowText','on','fill','on','linecolor','r');
% % colormap(flip(gray))
% % clabel(Cc,h,'FontSize',14,'Color','red')
% % title('Predicted megadrought probability','fontsize',20)
% % set(gca,'fontsize',16)
% % xlabel('Nino3.4 (annual anom. deg. C)','fontsize',18)
% % ylabel('AMO (anom. deg. C)','fontsize',18)
% % 
% % if prntfg=='y';print(['./figs/logreg_mdrought_prob_nino_amo.png'],'-dpng','-r300');end
% % 
% % 
% % figure
% % [Cc,h]=contour(mx,my,pL,'ShowText','on','linecolor','k');
% % clabel(Cc,h,'FontSize',14,'Color','black')
% % title('Probability 95% conf.','fontsize',20)
% % set(gca,'fontsize',16)
% % xlabel('Nino3.4 (annual anom. deg. C)','fontsize',18)
% % ylabel('AMO (anom. deg. C)','fontsize',18)
% % 
% % if prntfg=='y';print(['./figs/logreg_mdrought_prob_nino_amo_unc.png'],'-dpng','-r300');end
% 
% 
% % Alternative plot with both probability and uncertainty
% figure
% hold on
% [C1,h1]=contour(mx,my,pM,'fill','on');
% set(gca,'tickdir','both')
% colormap(flip(gray))
% c=colorbar;
% c.Label.String = 'Drought probability';
% [Cc,h]=contour(mx,my,pL,'showtext','on','linecolor',clrs.red,'linestyle','-.');
% clabel(Cc,h,'Color',clrs.red)
% set(gca,'fontsize',16)
% xlabel('Nino3.4 (annual anom. deg. C)','fontsize',18)
% ylabel('AMO (anom. deg. C)','fontsize',18)
% box on
% 
% if prntfg=='y';print(['./figs/logreg_mdrought_prob_nino_amo_data&unc.png'],'-dpng','-r300');end
% 
% 
% %-------------------------
% % TEST GMT INSTEAD OF AMO
% %-------------------------
% 
% clear X Xm Xi
% 
% % Bivariate logistic regression
% X(:,1)=mean(xa_nino_ann(ir4,:)-nino_off,2);
% X(:,2)=mean(xa_gmt_ens(ir4,:)-gmt_off,2);
% 
% y0=zeros(reconYrs,1);
% y0(ird_ann)=1;
% Y=categorical(y0(ir4));
% [B,dev,stats] = mnrfit(X,Y);
% 
% % Bivariate logistic regression with interaction
% Xi(:,1)=mean(xa_nino_ann(ir4,:)-nino_off,2);
% Xi(:,2)=mean(xa_gmt_ens(ir4,:)-gmt_off,2);
% Xi(:,3)=Xi(:,1).*Xi(:,2);
% [Bi,devi,statsi] = mnrfit(Xi,Y);
% 
% 
% 
% % Explore the predicted probabilities using the range of ENSO and AMO
% 
% % Define a grid of all possible Nino and AMO values
% %ns=xa_nino_ann(ir4,:)-nino_off;as=xa_gmt_ens(ir4,:)-gmt_off;
% % Just pick the mean values over which the model is built
% ns=mean(xa_nino_ann(ir4,:)-nino_off,2);as=mean(xa_gmt_ens(ir4,:)-gmt_off,2);
% [mx,my] = meshgrid(min(ns(:)):0.01:max(ns(:)),min(as(:)):0.01:max(as(:)));
% Xm(:,1)=mx(:);
% Xm(:,2)=my(:);
% 
% % Get probability for random observations
% pihat = mnrval(B,Xm);
% [pihat,dlow,dhi] = mnrval(B,Xm,stats);% get uncertainties
% % NOTE THAT UPPER AND LOWER BOUNDS ARE THE SAME!!
% 
% % Reshape to match the grid
% pM=reshape(pihat(:,2),size(mx));
% pL=reshape(dlow(:,2),size(mx));
% 
% 
% % Alternative plot with both probability and uncertainty
% figure
% hold on
% [C1,h1]=contour(mx,my,pM,'fill','on');
% set(gca,'tickdir','both')
% colormap(flip(gray))
% c=colorbar;
% c.Label.String = 'Drought probability';
% [Cc,h]=contour(mx,my,pL,'showtext','on','linecolor',clrs.red,'linestyle','-.');
% clabel(Cc,h,'Color',clrs.red)
% set(gca,'fontsize',16)
% xlabel('Nino3.4 (annual anom. deg. C)','fontsize',18)
% ylabel('GMT (anom. deg. C)','fontsize',18)
% box on
% 
% if prntfg=='y';print(['./figs/logreg_mdrought_prob_nino_gmt.png'],'-dpng','-r300');end
% 
% 
% 
% 
% 
% 








