


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


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


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

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




%==============================
% CALCULATE THE HURST EXPONENT
%==============================

% Using R/S method
H(1)=hurst(Xa_di(ir1),19);
H(2)=hurst(Xa_di(ir2),19);
H(3)=hurst(xa_nino_mn(irm1_a:irm1_b),19*12);
H(4)=hurst(xa_nino_mn(irm2_a:irm2_b),19*12);
H(5)=hurst(mean(xa_amo_ens(ir1,:),2),19);
H(6)=hurst(mean(xa_amo_ens(ir2,:),2),19);
H(7)=hurst(gmt_est_f(ir1),19);
H(8)=hurst(gmt_est_f(ir2),19);

disp('Hurst exponent: ')
H


% % ADJUST THE YEAR ENDPOINTS FOR SENSITIVITIES
% % **** The coe below shows that there are no sensitivities to the particular years chosen
% for i=0:5:50;
% 
%    analyrs1=800:(1575+i);
%    analyrs2=(1576+i):1925;
%    [~,~,ir1] = intersect(analyrs1,r_o:r_f);[~,~,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;
% 
% 
%    H2(1)=hurst(Xa_di(ir1),19);H2(2)=hurst(Xa_di(ir2),19);
%    H2(3)=hurst(xa_nino(irm1_a:irm1_b),19*12);H2(4)=hurst(xa_nino(irm2_a:irm2_b),19*12);
% 
%    H2
% 
% end
% 
% 
% % Same differences if the aggregated variance method is used
% hurst_aggvar(Xa_di(ir1));
% hurst_aggvar(Xa_di(ir2));
% hurst_aggvar(xa_nino(irm1_a:irm1_b));
% hurst_aggvar(xa_nino(irm2_a:irm2_b));

