
ens='010';
varpth=['/d6/haibo/lme/cam5/BLMTRC5CN/ALLFORCEINGS/' ens];
lon=ncread([varpth '/TREFHT/b.e11.BLMTRC5CN.f19_g16.' ens '.cam.h0.TREFHT.085001-184912.nc'],'lon');
lat=ncread([varpth '/TREFHT/b.e11.BLMTRC5CN.f19_g16.' ens '.cam.h0.TREFHT.085001-184912.nc'],'lat');

%======================
%======  SADA  ========
%======================

fln='/home/nsteiger/Dropbox/sa-phyda/SADA_1400_2000.nc';
da=flip(permute(ncread(fln,'scpdsi'),[2 1 3]),1);
%nada(nada<-99)=NaN; % put nans in fill spots
lat_da=flip(ncread(fln,'latitude'),1);
lon_da=ncread(fln,'longitude');
yrs_da=ncread(fln,'Time');


% Coarsen SADA to match climate model grid
ndm=load('sada_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.sadamsk==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


% Land mask over Argentina
msk=load('argentina_mask.mat');
lndmsk=msk.argntn;
% Land mask over Chile
%msk=load('chile_mask.mat');
%lndmsk=msk.chl;


% Latitude weighting
A=cosd(repmat(lat,[1 length(lon)]));

% Compute PDSI time series based on masked region
v_msk=bsxfun(@times,da_i,lndmsk);
bx_sp=wmean_a(v_msk,A); % RECONSTRUCTION DROUGHT INDEX


%=====================================
% LOAD ENSEMBLE MEGADROUGHT ESTIMATES

%load('pampas_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('pampas_megaD_ens_noNAproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('pampas_megaD_ens_noNASWproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
load('argentina_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('chile_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('chlarg_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
xa_di_sa=Xa_di;
xa_di_ens_sa=Xa_di_ens;
ds_sa=svmd;
mdagr_sa=ens_drght;
sa_yrs=analyrs;




% Indices of the drought atlas
analyrs1=1400:2000;
[~,ir1,id1] = intersect(1:2000,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



load distinct_colors_20.mat

% time series to compare
sadatm=bx_sp(id1);
sadatm=sadatm-mean(sadatm);
phydatm=xa_di_sa(ir1);
phydatm=phydatm-mean(phydatm);

% PHYDA AND HIGHLIGHT MEGADROUGHTS
figure
hold on
h1=plot(analyrs1,sadatm,'linewidth',1,'color',[clrs.blue 0.5]); % extra number on the end adds transparency
h2=plot(analyrs1,phydatm,'linewidth',1,'color',[clrs.purple 0.5]);
h3=plot(analyrs1,smooth(sadatm,11),'linewidth',2,'color',clrs.blue);
h4=plot(analyrs1,smooth(phydatm,11),'linewidth',2,'color',clrs.purple);
plot(analyrs1,zeros(length(analyrs1),1),'--','linewidth',2,'color',clrs.yellow)
set(gca,'fontsize',16)
%ylim([-0.8 1])
xlim([analyrs1(1) analyrs1(end)])
xlabel('Year','fontsize',18)
ylabel(['PDSI'],'fontsize',18)
box on
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]

%y=Xa_di(ir3); % using 'l' from above
%yMin=min(y)-0.5;yMax=max(y)+0.5;
ax=gca;yMin=ax.YLim(1);yMax=ax.YLim(2);
ylim([yMin yMax])


legend([h3,h4],{'SADA','PHYDA'},'location','northwest')
hold off

%print(['/home/nsteiger/Dropbox/sa-phyda/phyda_sada_chile.png'],'-dpng','-r300');
print(['/home/nsteiger/Dropbox/sa-phyda/phyda_sada_argentina.png'],'-dpng','-r300');











