

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=800:1600;
analyrs2=1601:1925;
analyrs3=800:2000;
analyrs4=800:1925;
analyrs5=1900:2000;
[~,~,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
[~,~,ir5] = intersect(analyrs5,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';
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);


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

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

Xv_msk=bsxfun(@times,Xt2m,nasw);
t2m_nasw=wmean_a(Xv_msk,A); % RECONSTRUCTION DROUGHT INDEX



%========
% T 2m  
%========

[Xobs,xlat,xlon,x_yrs] = load_obs_S('bearth',mon_avg_o,mon_avg_f);
da=reshape(Xobs,length(xlat),length(xlon),length(x_yrs));

% Interpolate to CESM resolution
%ndm=load('dai_pdsi_mask.mat');
% create grids for interpolation
[xn,yn]=meshgrid(xlon,xlat);
[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

disp('Interpolating obs to climate model res...')
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(:),yn(:),da1(:),'linear','linear');
   vq=F(xc,yc); % do the interpolation
   da_i(:,:,i)=vq;
end

%bet2m=reshape(da_i,length(lat)*length(lon),size(da,3));


bv_msk=bsxfun(@times,da_i,nasw);
bet2m_nasw=wmean_a(bv_msk,A);

% Indices of the obs
%[~,~,ib3] = intersect(analyrs3,x_yrs); 
[~,~,ib5] = intersect(analyrs5,x_yrs); 

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




%======================
%======  PHYDA  =======
%======================

% k=find(state_tp=='q');
% 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;
% 
% % Latitude weighting
% A=cosd(repmat(xmeta{k}.lat,[1 length(xmeta{k}.lon)]));
% 
% % PDSI ensemble
% Xa_di_ens=zeros(reconYrs,sub_ens);
% for i=1:sub_ens
%    Xvar=squeeze(Xa_ens(xmeta{k}.id_X(1):xmeta{k}.id_X(2),i,:));
%    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);
%    Xv_msk=bsxfun(@times,Xpdsi,nasw);
%    Xa_di_ens(:,i)=wmean_a(Xv_msk,A); % RECONSTRUCTION DROUGHT INDEX
% end
% 
% % Get the ensemble mean; identical to the original reconstruction
% Xvar=squeeze(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);
% Xv_msk=bsxfun(@times,Xpdsi,nasw);
% 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 ENSEMBLE MEGADROUGHT ESTIMATES

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

%======================
%======  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;
% A=cosd(repmat(xmeta{k}.lat,[1 length(xmeta{k}.lon)]));
% 
% DA_msk=bsxfun(@times,da_i,nasw);
% DA_di=wmean_a(DA_msk,A); % DROUGHT ATLAS DROUGHT INDEX
% % Indices of the drought atlas
% analyrs5=1870:2000;
% [~,~,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
% [~,~,id4] = intersect(analyrs4,yrs_da); % indices of recon segment
% [~,~,id5] = intersect(analyrs5,yrs_da); % indices of recon segment
% 
% [~,~,ir5] = intersect(analyrs5,r_o:r_f); % indices of recon segment
% 
% 
% % MEGADROUGHTS IN NADA	
% nada_ds4=mdrought(DA_di(ir4),1,'y'); % 'y' indicates that the mean is removed


%===========================================================================
% 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); 
[~,~,ig5] = intersect(analyrs5,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));

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


load distinct_colors_20.mat

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


% How useful is the EBM for temperature in the NASW? Correlations between different quantities
gf5=gmt_est_f(ig5);
gmt5=mean(xa_gmt_ens(ir5,:),2);
be5=bet2m_nasw(ib5);

%disp(['corr(forc,xa_gmt) ' num2str(analyrs5(1)) '-' num2str(analyrs5(end)) ' = ' num2str(corr(gf5,gmt5))])

disp(['corr(forc,be_nasw) ' num2str(analyrs5(1)) '-' num2str(analyrs5(end)) ' = ' num2str(corr(gf5,be5))])

gf3=gmt_est_f(ig3);
gmt3=mean(xa_gmt_ens(ir3,:),2);

%disp(['corr(forc,xa_gmt) ' num2str(analyrs3(1)) '-' num2str(analyrs3(end)) ' = ' num2str(corr(gf3,gmt3))])

gf310=ann2xt(gmt_est_f(ig3),10);
gmt310=ann2xt(mean(xa_gmt_ens(ir3,:),2),10);

%disp(['corr(forc_10,xa_gmt_10) ' num2str(analyrs3(1)) '-' num2str(analyrs3(end)) ' = ' num2str(corr(gf310,gmt310))])


t3=t2m_nasw(ir3);

disp(['corr(forc,t2m_nasw) ' num2str(analyrs3(1)) '-' num2str(analyrs3(end)) ' = ' num2str(corr(gf3,t3))])

disp(['corr(forc_10,t2m10) ' num2str(analyrs3(1)) '-' num2str(analyrs3(end)) ' = ' num2str(corr(gf310,ann2xt(t3,10)))])


% Are drought years responding stronger to forcing than non-drought years?
xD=Xa_di(ir4);
gf=gmt_est_f(ig4);

figure
plot(gf,xD,'.')
xlabel('Forc')
ylabel('NASW PDSI')


[xr,lgs]=xcorr(gf,xD,20,'coeff');
figure;plot(lgs,xr)
title('Cross correlation with all data')



volyr=find(isoutlier(gf)==1);

xD(volyr)=[];
gf(volyr)=[];

figure
plot(gf,xD,'.r')
xlabel('Forc')
ylabel('NASW PDSI')


[xr,lgs]=xcorr(gf,xD,20,'coeff');
figure;plot(lgs,xr)
title('Cross correlation removing volcanic eruptions outliers')


% The negative correlations indicate a positive feedback between forcing, temperature, and PDSI?



