
%--------------------------------------------------------------
% COMPARE THE SPECTRA OF THE CESM PRIOR WITH THE RECONSTRUCTION
%--------------------------------------------------------------





%=============================================================================================
%           LOAD THE UNCHANGED VARIABLES FROM CESM
%=============================================================================================


% Load reconstructed indices
mdl_prior='cesm_lme010';
gcm_eyr=1849; gcm_syr=850;
mon_avg_o=4;mon_avg_f=3;
p_yrs=[851:1848]-gcm_syr+1;

% Enter state vector data types (no spaces)
%  surface (t)emperature, 500hPa (g)eopotential,
%  (p)recip, surface p(r)essure, (e) P-E
%  (a)moc ind., am(o), ocean (h)eat OR N.(h)eat trans, (u)wind200mb
%  (s)ea level pressure or (s)st, (w)ind, soi(l) moisture, (z)onal precip
%  OR itc(z), stream (f)unction, 
%  PDSI (r) [full field], PDSI (q) [land only]
%  SPEI12 (i) [full field], SPEI12 (j) [land only] 
%  PDSI NAS(w) index, PDSI (e)ast Africa index
%  (n)ino 3.4, sea surface salinit(y), gmt (2)m
state_tp = 'on';

% Load AMO and Nino variables at annual resolution
[Xb_o,xmeta]=load_S_nobc(mdl_prior,state_tp,mon_avg_o,mon_avg_f,p_yrs);

k=find(state_tp=='n');i=3; % = Nino 3.4
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=='o'); % AMO
xb_amo=Xb_o((xmeta{k}.id_X(1)),:);

% Load JJA PDSI
mon_avg_o=6;mon_avg_f=8;
state_tp = 'q';

[Xb_o,xmeta]=load_S_nobc(mdl_prior,state_tp,mon_avg_o,mon_avg_f,p_yrs);

% PDSI
k=find(state_tp=='q');
Xvar=Xb_o(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);
Xpdsi=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),length(p_yrs));
Xpdsi(xmeta{k}.lndidx,:)=Xvar;
Xpdsi=reshape(Xpdsi,length(xmeta{k}.lat),length(xmeta{k}.lon),length(p_yrs));

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)]));
Xb_di=wmean_a(Xv_msk,A); % RECONSTRUCTION DROUGHT INDEX


%=============================================================================================
%           LOAD THE RECONSTRUCTION VARIABLES
%=============================================================================================

% ALL PROXY TYPES, LMR DATABASE
efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:33:07.mat';
pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';load([pth,efl])


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

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


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

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

% Load colors for plotting
load distinct_colors_20.mat

% What years to look at?
analyrs1=850:1849; % period consistent with CESM
analyrs2=800:2000;
analyrs3=1870: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
% 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 OBSERVATIONAL DATA
load('/d1/nsteiger/climate-data/climate-indices/AMO_KaplanSST.mat')
a0=amo_mon';
amo=mon2ann(a0(:),4,3);
ayrs=ayrs(1:end-1);% last year removed because of averaging
[~,~,ia3] = intersect(analyrs3,ayrs); % indices of recon segment

load('/d1/nsteiger/climate-data/climate-indices/nino34_NOAA_HadISST_1870_2016.mat')
n0=nino34(:,2:end)';
nino0=n0(:);
%[~,~,in3] = intersect(analyrs3,yrs); % indices of recon segment
ovrlprng=rng_intr([yrs(1),yrs(end)],[analyrs3(1),analyrs3(end)]);
mon_r=yrs(1):(1/12):(yrs(end)+11/12);
inm3_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
inm3_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;


%======================
%======  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
[~,~,id5] = intersect(analyrs5,yrs_da); % indices of recon segment

[~,~,ir5] = intersect(analyrs5,r_o:r_f); % indices of recon segment


% PLOT PHYDA AND NADA DROUGHT TIME SERIES

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

% Compute percentiles
%amo_p=prctile(xa_amo_ens,[5 50 95],2);

figure
hold on
plot(analyrs2,DA_di(id2),'linewidth',2,'color',clrs.black)
plot(analyrs2,Xa_di(ir2),'linewidth',2,'color',clrs.red)
legend('NADA','PHYDA','location','southeast')
set(gca,'fontsize',16)
%ylim([-0.8 1])
xlim([analyrs2(1) analyrs2(end)])
xlabel('Year','fontsize',18)
ylabel(['NASW PDSI'],'fontsize',18)
box on
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]

%print('-dpng',['./figs/timeseries_phyda_nada_' num2str(analyrs3(1)) '_' num2str(analyrs3(end)) '.png'],'-r300')

%===============================
%         NINO 3.4
%===============================


% [Pxy,F] = cpsd(x,y,window,noverlap,nfft,fs)
% [pxx,f] = periodogram(x,window,f,fs)
% [pxx,f] = pwelch(x,window,noverlap,f,fs)
%[p_t,~]=periodogram(detrend(t_nino(ih)),hanning(length(t_nino(ih))),[],12);
%[p_r,f]=periodogram(detrend(xa_nino(ir)),hanning(length(t_nino(ih))),[],12);


[p_o,fo,p_oc]=pwelch(detrend(nino0(inm3_a:inm3_b)),round(length(nino0(inm3_a:inm3_b))/8),[],[],12,'ConfidenceLevel',0.9);
[p_0,f0,p_0c]=pwelch(detrend(xb_nino),round(length(xb_nino)/8),[],[],12,'ConfidenceLevel',0.9);
[p_r,fr,p_rc]=pwelch(detrend(xa_nino),round(length(xa_nino)/8),[],[],12,'ConfidenceLevel',0.9);
[p_2,f2,p_2c]=pwelch(detrend(xa_nino(irm2_a:irm2_b)),round(length(xa_nino(irm2_a:irm2_b))/8),[],[],12,'ConfidenceLevel',0.9);
[p_1,f1,p_1c]=pwelch(detrend(xa_nino(irm1_a:irm1_b)),round(length(xa_nino(irm1_a:irm1_b))/8),[],[],12,'ConfidenceLevel',0.9);
[p_3,f3,p_3c]=pwelch(detrend(xa_nino(irm3_a:irm3_b)),round(length(xa_nino(irm3_a:irm3_b))/8),[],[],12,'ConfidenceLevel',0.9);


% Validation figures, OBS and PHYDA
figure
hold on
ho=plot(fo,real(p_o),'linewidth',2,'color',clrs.black);
plot(fo,real(p_oc),'--','color',clrs.black);
h3=plot(f3,real(p_3),'linewidth',2,'color',clrs.red);
plot(fo,real(p_3c),'--','color',clrs.red);
legend([ho,h3],{'HadISST 1870-2000',['PHYDA ' num2str(analyrs3(1)) '-' num2str(analyrs3(end))]})
xlim([0 0.5]) % limit view to 100 years
ax = gca;
ax.XTick = [0 0.1 0.2 0.3 0.4 0.5];
ax.XTickLabel = {'0','0.1 [10]','0.2 [5]','0.3 [3.3]','0.4 [2.5]','0.5 [2]'};
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 7 4]; % [left bottom width height]
%title('Nino3.4 index')
box on
xlabel('Frequency (1/yrs), Period [yrs]','fontsize',14)
ylabel('Power spectral density','fontsize',14)
set(gca,'fontsize',14)
hold off

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


% How does PHYDA compare with CESM?
figure
hold on
h0=plot(f0,real(p_0),'linewidth',2,'color',clrs.grey);
plot(f0,real(p_0c),'--','color',clrs.grey);
h1=plot(f1,real(p_1),'linewidth',2,'color',clrs.blue);
plot(f1,real(p_1c),'--','color',clrs.blue);
h3=plot(f3,real(p_3),'linewidth',2,'color',clrs.red);
plot(f3,real(p_3c),'--','color',clrs.red);
legend([h0,h1,h3],{'CESM 850-1849',['PHYDA ' num2str(analyrs1(1)) '-' num2str(analyrs1(end))],...
['PHYDA ' num2str(analyrs3(1)) '-' num2str(analyrs3(end))]})
%,['PHYDA ' num2str(analyrs2(1)) '-' num2str(analyrs2(end))]...
%,['PHYDA ' num2str(analyrs3(1)) '-' num2str(analyrs3(end))])
xlim([0 0.5]) % limit view to 100 years
ax = gca;
ax.XTick = [0 0.1 0.2 0.3 0.4 0.5];
ax.XTickLabel = {'0','0.1 [10]','0.2 [5]','0.3 [3.3]','0.4 [2.5]','0.5 [2]'};
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 7 4]; % [left bottom width height]
%title('Nino3.4 index')
box on
xlabel('Frequency (1/yrs), Period [yrs]','fontsize',14)
ylabel('Power spectral density','fontsize',14)
set(gca,'fontsize',14)
hold off

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



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



[p_o,fo,p_oc]=pwelch(detrend(amo(ia3)),round(length(amo(ia3))/8),[],[],1,'ConfidenceLevel',0.9);
[p_0,f0,p_0c]=pwelch(detrend(xb_amo),round(length(xb_amo)/8),[],[],1,'ConfidenceLevel',0.9);
[p_r,fr,p_rc]=pwelch(detrend(xa_amo),round(length(xa_amo)/8),[],[],1,'ConfidenceLevel',0.9);
[p_2,f2,p_2c]=pwelch(detrend(xa_amo(ir2)),round(length(xa_amo(ir2))/8),[],[],1,'ConfidenceLevel',0.9);
[p_1,f1,p_1c]=pwelch(detrend(xa_amo(ir1)),round(length(xa_amo(ir1))/8),[],[],1,'ConfidenceLevel',0.9);
[p_3,f3,p_3c]=pwelch(detrend(xa_amo(ir3)),round(length(xa_amo(ir3))/8),[],[],1,'ConfidenceLevel',0.9);

figure
hold on
ho=plot(fo,real(p_o),'linewidth',2,'color',clrs.black);
plot(fo,real(p_oc),'--','color',clrs.black);
h3=plot(f3,real(p_3),'linewidth',2,'color',clrs.red);
plot(fo,real(p_3c),'--','color',clrs.red);
legend([ho,h3],{'HadISST 1870-2000',['PHYDA ' num2str(analyrs3(1)) '-' num2str(analyrs3(end))]})
xlim([0 0.5]) % limit view to 100 years
ax = gca;
ax.XTick = [0 0.1 0.2 0.3 0.4 0.5];
ax.XTickLabel = {'0','0.1 [10]','0.2 [5]','0.3 [3.3]','0.4 [2.5]','0.5 [2]'};
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 7 4]; % [left bottom width height]
%title('AMO index')
box on
xlabel('Frequency (1/yrs), Period [yrs]','fontsize',14)
ylabel('Power spectral density','fontsize',14)
set(gca,'fontsize',14)
hold off

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

figure
hold on
h0=plot(f0,real(p_0),'linewidth',2,'color',clrs.grey);
plot(f0,real(p_0c),'--','color',clrs.grey);
h1=plot(f1,real(p_1),'linewidth',2,'color',clrs.blue);
plot(f1,real(p_1c),'--','color',clrs.blue);
h3=plot(f3,real(p_3),'linewidth',2,'color',clrs.red);
plot(f3,real(p_3c),'--','color',clrs.red);
legend([h0,h1,h3],{'CESM 850-1849',['PHYDA ' num2str(analyrs1(1)) '-' num2str(analyrs1(end))],...
['PHYDA ' num2str(analyrs3(1)) '-' num2str(analyrs3(end))]})
xlim([0 0.5]) % limit view to 100 years
ax = gca;
ax.XTick = [0 0.1 0.2 0.3 0.4 0.5];
ax.XTickLabel = {'0','0.1 [10]','0.2 [5]','0.3 [3.3]','0.4 [2.5]','0.5 [2]'};
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 7 4]; % [left bottom width height]
%title('AMO index')
box on
xlabel('Frequency (1/yrs), Period [yrs]','fontsize',14)
ylabel('Power spectral density','fontsize',14)
set(gca,'fontsize',14)
hold off

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

%===============================
%         PDSI
%===============================



[p_o,fo,p_oc]=pwelch(detrend(DA_di(id3)),round(length(DA_di(id3))/8),[],[],1,'ConfidenceLevel',0.9);
[p_0,f0,p_0c]=pwelch(detrend(Xb_di),round(length(Xb_di)/8),[],[],1,'ConfidenceLevel',0.9);
[p_r,fr,p_rc]=pwelch(detrend(Xa_di),round(length(Xa_di)/8),[],[],1,'ConfidenceLevel',0.9);
[p_2,f2,p_2c]=pwelch(detrend(Xa_di(ir2)),round(length(Xa_di(ir2))/8),[],[],1,'ConfidenceLevel',0.9);
[p_1,f1,p_1c]=pwelch(detrend(Xa_di(ir1)),round(length(Xa_di(ir1))/8),[],[],1,'ConfidenceLevel',0.9);
[p_3,f3,p_3c]=pwelch(detrend(Xa_di(ir3)),round(length(Xa_di(ir3))/8),[],[],1,'ConfidenceLevel',0.9);

figure
hold on
ho=plot(fo,real(p_o),'linewidth',2,'color',clrs.black);
plot(fo,real(p_oc),'--','color',clrs.black);
h3=plot(f3,real(p_3),'linewidth',2,'color',clrs.red);
plot(fo,real(p_3c),'--','color',clrs.red);
legend([ho,h3],{['NADA ' num2str(analyrs3(1)) '-' num2str(analyrs3(end))]...
,['PHYDA ' num2str(analyrs3(1)) '-' num2str(analyrs3(end))]})
xlim([0 0.5]) % limit view to 100 years
ax = gca;
ax.XTick = [0 0.1 0.2 0.3 0.4 0.5];
ax.XTickLabel = {'0','0.1 [10]','0.2 [5]','0.3 [3.3]','0.4 [2.5]','0.5 [2]'};
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 7 4]; % [left bottom width height]
%title('NASW PDSI')
box on
xlabel('Frequency (1/yrs), Period [yrs]','fontsize',14)
ylabel('Power spectral density','fontsize',14)
set(gca,'fontsize',14)
hold off

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


figure
hold on
h0=plot(f0,real(p_0),'linewidth',2,'color',clrs.grey);
plot(f0,real(p_0c),'--','color',clrs.grey);
h1=plot(f1,real(p_1),'linewidth',2,'color',clrs.blue);
plot(f1,real(p_1c),'--','color',clrs.blue);
h3=plot(f3,real(p_3),'linewidth',2,'color',clrs.red);
plot(f3,real(p_3c),'--','color',clrs.red);
legend([h0,h1,h3],{'CESM 850-1849',['PHYDA ' num2str(analyrs1(1)) '-' num2str(analyrs1(end))],...
['PHYDA ' num2str(analyrs3(1)) '-' num2str(analyrs3(end))]})
xlim([0 0.5]) % limit view to 100 years
ax = gca;
ax.XTick = [0 0.1 0.2 0.3 0.4 0.5];
ax.XTickLabel = {'0','0.1 [10]','0.2 [5]','0.3 [3.3]','0.4 [2.5]','0.5 [2]'};
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 7 4]; % [left bottom width height]
%title('AMO index')
box on
xlabel('Frequency (1/yrs), Period [yrs]','fontsize',14)
ylabel('Power spectral density','fontsize',14)
set(gca,'fontsize',14)
hold off

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


%=============================================================================================
%           COMPUTE THE CROSS CORRELATIONS
%=============================================================================================


% AMO AND PDSI

[ro,lo]=xcorr(detrend(DA_di(id5)),detrend(xa_amo(ir5))',30,'coeff');
[r0,l0]=xcorr(detrend(Xb_di),detrend(xb_amo)',30,'coeff');
[rr,lr]=xcorr(detrend(Xa_di),detrend(xa_amo)',30,'coeff');
[r1,l1]=xcorr(detrend(Xa_di(ir1)),detrend(xa_amo(ir1))',30,'coeff');
[r2,l2]=xcorr(detrend(Xa_di(ir2)),detrend(xa_amo(ir2))',30,'coeff');
[r3,l3]=xcorr(detrend(Xa_di(ir3)),detrend(xa_amo(ir3))',30,'coeff');

figure
hold on
plot(lo,ro,'linewidth',2,'color',clrs.black)
%plot(l0,r0,'--','linewidth',2,'color',[0.4 0.4 0.4])
%plot(lr,rr,'linewidth',2)
%plot(l1,r1,'linewidth',2)
%plot(l2,r2,'linewidth',2)
plot(l3,r3,'linewidth',2,'color',clrs.red)
legend(['NADA & HadISST ' num2str(analyrs5(1)) '-' num2str(analyrs5(end))],...
['PHYDA ' num2str(analyrs3(1)) '-' num2str(analyrs3(end))],'location','northwest')
xlim([l0(1) l0(end)])
xlabel('Lag (yrs)','fontsize',18)
ylabel('Cross-correlation','fontsize',18)
%title('NASW PDSI & AMO')
box on
grid on
ylim([-0.3 0.7])
set(gca,'fontsize',16)
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 6 4]; % [left bottom width height]
%legend('CESM',['Xa ' num2str(r_o) '-' num2str(r_f)])

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

figure
hold on
%plot(lo,ro,'linewidth',2,'color',clrs.black)
plot(l0,r0,'linewidth',2,'color',clrs.grey)
%plot(lr,rr,'linewidth',2)
plot(l1,r1,'linewidth',2,'color',clrs.blue)
%plot(l2,r2,'linewidth',2)
plot(l3,r3,'linewidth',2,'color',clrs.red)
legend(['CESM 850-1849'],...
['PHYDA ' num2str(analyrs1(1)) '-' num2str(analyrs1(end))],...
['PHYDA ' num2str(analyrs3(1)) '-' num2str(analyrs3(end))],'location','northwest')
xlim([l0(1) l0(end)])
xlabel('Lag (yrs)','fontsize',18)
ylabel('Cross-correlation','fontsize',18)
%title('NASW PDSI & AMO')
box on
grid on
ylim([-0.3 0.7])
set(gca,'fontsize',16)
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 6 4]; % [left bottom width height]
%legend('CESM',['Xa ' num2str(r_o) '-' num2str(r_f)])

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



% NINO 3.4 AND PDSI


% Annualize Nino for the computations
xb_ninoA=mon2ann(xb_nino,4,3);
xa_ninoA=mon2ann(xa_nino,4,3);

id5n=id5(1:end-1);
ir5n=ir5(1:end-1);
ir3n=ir3(1:end-1);
ir2n=ir2(1:end-1);

[ro,lo]=xcorr(detrend(DA_di(id5n)),detrend(xa_ninoA(ir5n))',10,'coeff');
[r0,l0]=xcorr(detrend(Xb_di(1:end-1)),detrend(xb_ninoA),10,'coeff');
[rr,lr]=xcorr(detrend(Xa_di(1:end-1)),detrend(xa_ninoA),10,'coeff');
[r1,l1]=xcorr(detrend(Xa_di(ir1)),detrend(xa_ninoA(ir1)),10,'coeff');
[r2,l2]=xcorr(detrend(Xa_di(ir2n)),detrend(xa_ninoA(ir2n)),10,'coeff');
[r3,l3]=xcorr(detrend(Xa_di(ir3n)),detrend(xa_ninoA(ir3n)),10,'coeff');

figure
hold on
plot(lo,ro,'linewidth',2,'color',clrs.black)
%plot(l0,r0,'linewidth',2,'color',[0.4 0.4 0.4])
%plot(lr,rr,'linewidth',2)
%plot(l1,r1,'linewidth',2)
%plot(l2,r2,'linewidth',2)
plot(l3,r3,'linewidth',2,'color',clrs.red)
legend(['NADA & HadISST ' num2str(analyrs5(1)) '-' num2str(analyrs5(end))],...
['PHYDA ' num2str(analyrs3(1)) '-' num2str(analyrs3(end))],'location','northwest')
xlim([l0(1) l0(end)])
ylim([-0.3 0.7])
xlabel('Lag (yrs)','fontsize',18)
ylabel('Cross-correlation','fontsize',18)
title('NASW PDSI & Nino3.4')
box on
grid on
set(gca,'fontsize',16)
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 6 4]; % [left bottom width height]
%legend('CESM',['Xa ' num2str(r_o) '-' num2str(r_f)])

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



figure
hold on
%plot(lo,ro,'linewidth',2,'color',clrs.black)
plot(l0,r0,'linewidth',2,'color',clrs.grey)
%plot(lr,rr,'linewidth',2)
plot(l1,r1,'linewidth',2,'color',clrs.blue)
%plot(l2,r2,'linewidth',2)
plot(l3,r3,'linewidth',2,'color',clrs.red)
legend(['CESM 850-1849'],...
['PHYDA ' num2str(analyrs1(1)) '-' num2str(analyrs1(end))],...
['PHYDA ' num2str(analyrs3(1)) '-' num2str(analyrs3(end))],'location','northwest')
xlim([l0(1) l0(end)])
ylim([-0.3 0.7])
xlabel('Lag (yrs)','fontsize',18)
ylabel('Cross-correlation','fontsize',18)
title('NASW PDSI & Nino3.4')
box on
grid on
set(gca,'fontsize',16)
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 6 4]; % [left bottom width height]
%legend('CESM',['Xa ' num2str(r_o) '-' num2str(r_f)])

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








