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

% Which drought information to use?
ds=svmd;

% 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;
% Compute annual mean ENSO
[nino_ann,nmeta] = nino(Xt2m,lat,lon,3,'n');
 
% 
% % 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')
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));


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

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


% 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);
% Nino offset with respect to reference period
nino_off=mean(mean(xa_nino_ens(irm2_a:irm2_b,:),2));



% VOLCANO YEARS
vlcyrs=[1258 1458 1230 1641 1171 1695 1345 1286 1182 1276];


% Pull out ENSO chunks
chkln=10*12;% chunk of time to pull out
enso_mon=zeros(chkln,length(p_yrs),length(vlcyrs));
enso_ann=zeros(chkln/12-1,length(p_yrs),length(vlcyrs));

for i=1:length(vlcyrs)

   ovrlprng=rng_intr([r_o,r_f],[vlcyrs(i)-3,vlcyrs(i)+6]);
   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;

   enso_mon(:,:,i)=xa_nino_ens(irm3_a:irm3_b,:);
   % annual mean, ensemble mean
   e1=mean(mon2ann(xa_nino_ens(irm3_a:irm3_b,:)',4,3),1);
   % remove reference
   enso_nma(:,i)=e1-mean(e1(1:3));
   enso_ann(:,:,i)=mon2ann(xa_nino_ens(irm3_a:irm3_b,:)',4,3)';

end

% random draws
%vlcyrs_r=randsample(1000:1900,1000,true);
vlcyrs_r=1000:1900;
enso_nma_r=zeros(9,length(vlcyrs_r));
for i=1:length(vlcyrs_r)

   ovrlprng=rng_intr([r_o,r_f],[vlcyrs_r(i)-3,vlcyrs_r(i)+6]);
   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;

   enso_mon(:,:,i)=xa_nino_ens(irm3_a:irm3_b,:);
   % annual mean, ensemble mean
   e1=mean(mon2ann(xa_nino_ens(irm3_a:irm3_b,:)',4,3),1);
   % remove reference
   enso_nma_r(:,i)=e1-mean(e1(1:3));

end

% randomly subselect sets of years
enso_mr=zeros(9,1000);
for i=1:1000
   ri=randsample(length(vlcyrs_r),length(vlcyrs));
   %enso_mr(:,i)=mean(enso_nma_r(:,ri),2);
   enso_mr(:,i)=median(enso_nma_r(:,ri),2);
end


% get percentiles of null
eprc=prctile(enso_mr,[5 95],2);

figure;plot(-3:5,enso_nma)
hold on
%plot(-3:5,mean(enso_nma,2),'linewidth',3)
plot(-3:5,median(enso_nma,2),'linewidth',3)
plot(-3:5,0*[1:9],'k--','linewidth',2)
plot(-3:5,eprc,'--','color',clrs.grey,'linewidth',2)%uncertainties
title('Nino3.4 response to volcanic eruptions')

fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 7 4]; % [left bottom width height]

%print(['/home/nsteiger/Dropbox/enso_volcanoe_sea_median.png'],'-dpng','-r300');












