
%============================================================
% EXTRACT RECONSTRUCTION VARIABLES AND SAVE TO NETCDF FORMAT
%============================================================

varpth='./output-recon/ens_100_21nov18/';

%efl='cesm_lme010_r12000_p998_state_2ozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_21-Nov-2018_10:46:49.mat';
%efl='cesm_lme010_r12000_p998_state_2ozn_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_21-Nov-2018_10:43:26.mat';
%efl='cesm_lme010_r12000_p998_state_2ozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_21-Nov-2018_10:42:13.mat';
%efl='cesm_lme010_r12000_p998_state_j_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_21-Nov-2018_20:56:57.mat';
%efl='cesm_lme010_r12000_p998_state_j_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_21-Nov-2018_21:12:22.mat';
%efl='cesm_lme010_r12000_p998_state_j_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_21-Nov-2018_21:20:05.mat';
%efl='cesm_lme010_r12000_p998_state_q_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_21-Nov-2018_20:55:55.mat';
%efl='cesm_lme010_r12000_p998_state_q_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_21-Nov-2018_20:48:05.mat';
%efl='cesm_lme010_r12000_p998_state_q_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_21-Nov-2018_20:47:28.mat';
%efl='cesm_lme010_r12000_p998_state_t_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_21-Nov-2018_10:59:46.mat';
%efl='cesm_lme010_r12000_p998_state_t_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_21-Nov-2018_11:05:22.mat';
efl='cesm_lme010_r12000_p998_state_t_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_21-Nov-2018_11:06:59.mat';


load([varpth,efl])

% Use lat/lon from CESM for it's variable attributes
% varpth=['/d6/haibo/lme/cam5/BLMTRC5CN/ALLFORCEINGS/' ens];
% lat=ncread([varpth '/TREFHT/b.e11.BLMTRC5CN.f19_g16.' ens '.cam.h0.TREFHT.085001-184912.nc'],'lat');
% ncks -A -v lat,lon b.e11.BLMTRC5CN.f19_g16.010.cam.h0.TREFHT.085001-184912.nc cesm_latlon.nc
inputf='./input-recon/cesm_latlon.nc';



% % define the output variable
% outputf1=['da_hydro_' mon_o mon_f '_r.' num2str(r_o) '-' num2str(r_f) '_d.' timestp(1:end-9) '.nc'];
% 
% % grab lat/lon from model, put it into the output file
% eval(['!ncks -A -v lat,lon ',inputf,' ',outputf1]);
% 
% % Years
% nccreate(outputf1,'time','Dimensions',{'time',reconYrs});
% ncwrite(outputf1,'time',r_o:r_f);
% ncwriteatt(outputf1,'time','long_name','Years CE');
% 
% if strcmp([mon_o mon_f],'AprMar')
%    ncwriteatt(outputf1,'time','units','Annual mean starting from April to the next calendar year March');
% elseif strcmp([mon_o mon_f],'DecFeb')
%    ncwriteatt(outputf1,'time','units','December-January-February average');
% elseif strcmp([mon_o mon_f],'JunAug')
%    ncwriteatt(outputf1,'time','units','June-July-August average');
% end

% Months
mon_r=r_o:(1/12):(r_f+11/12);

% nccreate(outputf1,'tmon','Dimensions',{'tmon',length(mon_r)});
% ncwrite(outputf1,'tmon',mon_r);
% ncwriteatt(outputf1,'tmon','long_name','Fractional years CE');
% ncwriteatt(outputf1,'tmon','units','Monthly mean');


%----------- SPATIAL VARIABLES -----------------


if ismember('t',state_tp)
   k=find(state_tp=='t');

   % define the output variable
   outputf1=['phyda_ens_tas_' mon_o mon_f '_r.' num2str(r_o) '-' num2str(r_f) '_d.' timestp(1:end-9) '.nc'];

   % grab lat/lon from model, put it into the output file
   eval(['!ncks -A -v lat,lon ',inputf,' ',outputf1]);

   % Years
   nccreate(outputf1,'time','Dimensions',{'time',reconYrs});
   ncwrite(outputf1,'time',r_o:r_f);
   ncwriteatt(outputf1,'time','long_name','Years CE');

   if strcmp([mon_o mon_f],'AprMar')
      ncwriteatt(outputf1,'time','units','Annual mean starting from April to the next calendar year March');
   elseif strcmp([mon_o mon_f],'DecFeb')
      ncwriteatt(outputf1,'time','units','December-January-February average');
   elseif strcmp([mon_o mon_f],'JunAug')
      ncwriteatt(outputf1,'time','units','June-July-August average');
   end

   % Extract ensemble	  
   X_ens0=Xa_ens(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:,:);
   X_ens=permute(reshape(X_ens0,length(xmeta{k}.lat),length(xmeta{k}.lon),sub_ens,reconYrs),[2 1 4 3]);
   clear X_ens0

   nccreate(outputf1,'tas_ens','Dimensions',{'lon',length(xmeta{k}.lon),'lat',length(xmeta{k}.lat),'time',reconYrs,'sub_ens',sub_ens});
   ncwrite(outputf1,'tas_ens',X_ens);
   ncwriteatt(outputf1,'tas_ens','long_name',[xmeta{k}.varnmlng ' (' num2str(sub_ens) ' member sub-ensemble)']);
   ncwriteatt(outputf1,'tas_ens','units','degC');


end


if ismember('q',state_tp)
   k=find(state_tp=='q');

   % define the output variable
   outputf1=['phyda_ens_pdsi_' mon_o mon_f '_r.' num2str(r_o) '-' num2str(r_f) '_d.' timestp(1:end-9) '.nc'];

   % grab lat/lon from model, put it into the output file
   eval(['!ncks -A -v lat,lon ',inputf,' ',outputf1]);

   % Years
   nccreate(outputf1,'time','Dimensions',{'time',reconYrs});
   ncwrite(outputf1,'time',r_o:r_f);
   ncwriteatt(outputf1,'time','long_name','Years CE');

   if strcmp([mon_o mon_f],'AprMar')
      ncwriteatt(outputf1,'time','units','Annual mean starting from April to the next calendar year March');
   elseif strcmp([mon_o mon_f],'DecFeb')
      ncwriteatt(outputf1,'time','units','December-January-February average');
   elseif strcmp([mon_o mon_f],'JunAug')
      ncwriteatt(outputf1,'time','units','June-July-August average');
   end
   
   % Extract ensemble
   X_ens0=Xa_ens(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:,:);
   X_ens0f=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),sub_ens,reconYrs);
   X_ens0f(xmeta{k}.lndidx,:,:)=X_ens0;
   clear X_ens0
   X_ens=permute(reshape(X_ens0f,length(xmeta{k}.lat),length(xmeta{k}.lon),sub_ens,reconYrs),[2 1 4 3]);
   clear X_ens0f

   nccreate(outputf1,'pdsi_ens','Dimensions',{'lon',length(xmeta{k}.lon),'lat',length(xmeta{k}.lat),'time',reconYrs,'sub_ens',sub_ens},'FillValue',NaN);
   ncwrite(outputf1,'pdsi_ens',X_ens);
   ncwriteatt(outputf1,'pdsi_ens','long_name',[xmeta{k}.varnmlng ' (' num2str(sub_ens) ' member sub-ensemble)']);
   ncwriteatt(outputf1,'pdsi_ens','units','Standardized units of relative dry and wet');
   ncwriteatt(outputf1,'pdsi_ens','missing_value',NaN);

end


if ismember('j',state_tp)
   k=find(state_tp=='j');

   % define the output variable
   outputf1=['phyda_ens_spei_' mon_o mon_f '_r.' num2str(r_o) '-' num2str(r_f) '_d.' timestp(1:end-9) '.nc'];

   % grab lat/lon from model, put it into the output file
   eval(['!ncks -A -v lat,lon ',inputf,' ',outputf1]);

   % Years
   nccreate(outputf1,'time','Dimensions',{'time',reconYrs});
   ncwrite(outputf1,'time',r_o:r_f);
   ncwriteatt(outputf1,'time','long_name','Years CE');

   if strcmp([mon_o mon_f],'AprMar')
      ncwriteatt(outputf1,'time','units','Annual mean starting from April to the next calendar year March');
   elseif strcmp([mon_o mon_f],'DecFeb')
      ncwriteatt(outputf1,'time','units','December-January-February average');
   elseif strcmp([mon_o mon_f],'JunAug')
      ncwriteatt(outputf1,'time','units','June-July-August average');
   end
   
   % Extract ensemble
   X_ens0=Xa_ens(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:,:);
   X_ens0f=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),sub_ens,reconYrs);
   X_ens0f(xmeta{k}.lndidx,:,:)=X_ens0;
   clear X_ens0
   X_ens=permute(reshape(X_ens0f,length(xmeta{k}.lat),length(xmeta{k}.lon),sub_ens,reconYrs),[2 1 4 3]);
   clear X_ens0f

   nccreate(outputf1,'spei_ens','Dimensions',{'lon',length(xmeta{k}.lon),'lat',length(xmeta{k}.lat),'time',reconYrs,'sub_ens',sub_ens},'FillValue',NaN);
   ncwrite(outputf1,'spei_ens',X_ens);
   ncwriteatt(outputf1,'spei_ens','long_name',[xmeta{k}.varnmlng ' (' num2str(sub_ens) ' member sub-ensemble)']);
   ncwriteatt(outputf1,'spei_ens','units','Standardized units of relative dry and wet');
   ncwriteatt(outputf1,'spei_ens','missing_value',NaN);


end



%-------- INDEX VARIABLES ------------

% index variables are all together in a recon, so assume they're all together

if strcmp('2ozn',state_tp)

   % define the output variable
   outputf1=['phyda_ens_indexvars_' mon_o mon_f '_r.' num2str(r_o) '-' num2str(r_f) '_d.' timestp(1:end-9) '.nc'];

   % grab lat/lon from model, put it into the output file
   eval(['!ncks -A -v lat,lon ',inputf,' ',outputf1]);

   % Years variable
   nccreate(outputf1,'time','Dimensions',{'time',reconYrs});
   ncwrite(outputf1,'time',r_o:r_f);
   ncwriteatt(outputf1,'time','long_name','Years CE');

   % Months variable
   nccreate(outputf1,'tmon','Dimensions',{'tmon',length(mon_r)});
   ncwrite(outputf1,'tmon',mon_r);
   ncwriteatt(outputf1,'tmon','long_name','Fractional years CE');
   ncwriteatt(outputf1,'tmon','units','Monthly mean');

   if strcmp([mon_o mon_f],'AprMar')
      ncwriteatt(outputf1,'time','units','Annual mean starting from April to the next calendar year March');
   elseif strcmp([mon_o mon_f],'DecFeb')
      ncwriteatt(outputf1,'time','units','December-January-February average');
   elseif strcmp([mon_o mon_f],'JunAug')
      ncwriteatt(outputf1,'time','units','June-July-August average');
   end


% NINO INDICES	
if ismember('n',state_tp)
   k=find(state_tp=='n');

   for i=1:length(xmeta{k}.indinfo)

      % Extract ensemble	  
      X_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,sub_ens);

      vrnm0=strrep(strrep(xmeta{k}.indinfo{i}.indxnm,' ','_'),'+',''); % remove plus and space
      vlngnm=xmeta{k}.indinfo{i}.indxnm;
      if i==5; vrnm0='PacDelSST'; vlngnm=strrep(vlngnm,':',''); end % simple kludge

      nccreate(outputf1,[vrnm0 '_ens'],'Dimensions',{'tmon',length(mon_r),'sub_ens',sub_ens});
      ncwrite(outputf1,[vrnm0 '_ens'],X_ens);
      ncwriteatt(outputf1,[vrnm0 '_ens'],'long_name',[vlngnm ' (' num2str(sub_ens) ' member sub-ensemble)']);
      ncwriteatt(outputf1,[vrnm0 '_ens'],'units','degC');

   end

end



% ITCZ INDICES	
if ismember('z',state_tp)
   k=find(state_tp=='z');
   vrnm0={'Atl','Pac130','Pac160','EPac','WPac','SAsia','IndOcn','SAmer','Ind','Afr','EAfr'};

   for i=1:length(xmeta{k}.indinfo)

      % Extract ensemble	  
      X_ens=permute(squeeze(Xa_ens(xmeta{k}.id_X(1)+i-1:xmeta{k}.id_X(1)+i-1,:,:)),[2,1]);

      nccreate(outputf1,[vrnm0{i} '_ens'],'Dimensions',{'time',reconYrs,'sub_ens',sub_ens});
      ncwrite(outputf1,[vrnm0{i} '_ens'],X_ens);
      ncwriteatt(outputf1,[vrnm0{i} '_ens'],'long_name',['ITCZ ' xmeta{k}.indinfo{i}.varnmlng ...
      ' (' num2str(sub_ens) ' member sub-ensemble)']);
      ncwriteatt(outputf1,[vrnm0{i} '_ens'],'units','lat (deg)');
    
   end

end


% GLOBLA MEAN INDEX		
if ismember('2',state_tp)
   k=find(state_tp=='2');

   % Extract ensemble	  
   X_ens=permute(squeeze(Xa_ens(xmeta{k}.id_X(1),:,:)),[2,1]);

   nccreate(outputf1,['gmt_ens'],'Dimensions',{'time',reconYrs,'sub_ens',sub_ens});
   ncwrite(outputf1,['gmt_ens'],X_ens);
   ncwriteatt(outputf1,['gmt_ens'],'long_name',[xmeta{k}.varnmlng ...
      ' (' num2str(sub_ens) ' member sub-ensemble)']);
   ncwriteatt(outputf1,['gmt_ens'],'units','degC');

end


% AMO INDEX		
if ismember('o',state_tp)
   k=find(state_tp=='o');

   % Extract ensemble	  
   X_ens=permute(squeeze(Xa_ens(xmeta{k}.id_X(1),:,:)),[2,1]);

   nccreate(outputf1,['amo_ens'],'Dimensions',{'time',reconYrs,'sub_ens',sub_ens});
   ncwrite(outputf1,['amo_ens'],X_ens);
   ncwriteatt(outputf1,['amo_ens'],'long_name',[xmeta{k}.varnmlng ...
      ' (' num2str(sub_ens) ' member sub-ensemble)']);
   ncwriteatt(outputf1,['amo_ens'],'units','degC');

end

end % END of index variables




%This example creates a new netCDF file, specifies a global attribute, and assigns a value to the attribute.
%ncid = netcdf.create('myfile.nc','CLOBBER');
%varid = netcdf.getConstant('GLOBAL');
%netcdf.putAtt(ncid,varid,'creation_date',datestr(now));
%netcdf.close(ncid);


% ADD GLOBAL ATTRIBUTES TO THE FILE
eval(['!ncatted -O --glb_att_add author="Nathan Steiger, LDEO Columbia University, Nov 2018" ' outputf1]);
eval(['!ncatted -O --glb_att_add citation="http://doi.org/10.1038/sdata.2018.86" ' outputf1]);

% DELETE THE HISTORY OF THE FILES
eval(['!ncatted -hO -a history_of_appended_files,global,d,, ' outputf1]);
eval(['!ncatted -hO -a history,global,d,, ' outputf1]);






























% BELOW IS CODE FOR EXTRACTING THE ENSEMBLES FOR BO LI TO USE

% if ismember('t',state_tp)
%    k=find(state_tp=='t');
% 
%    % define the output variable
%    outputf1=['tas_mn_da_hydro_r.' num2str(r_o) '-' num2str(r_f) '_d.' timestp(1:end-9) '.nc'];
%    outputf2=['tas_ens_da_hydro_r.' num2str(r_o) '-' num2str(r_f) '_d.' timestp(1:end-9) '.nc'];
%    outputf3=['tas_prior_da_hydro_r.' num2str(r_o) '-' num2str(r_f) '_d.' timestp(1:end-9) '.nc'];
% 
%    % grab lat/lon from model, put it into the output file
%    eval(['!ncks -A -v lat,lon ',inputf,' ',outputf1]);
%    eval(['!ncks -A -v lat,lon ',inputf,' ',outputf2]);
%    eval(['!ncks -A -v lat,lon ',inputf,' ',outputf3]);
% 
%    nccreate(outputf1,'time','Dimensions',{'time',reconYrs});
%    ncwrite(outputf1,'time',r_o:r_f);
%    ncwriteatt(outputf1,'time','units','Annual mean starting from April to the next calendar year March');
%    ncwriteatt(outputf1,'time','long_name','Years CE');
% 
%    nccreate(outputf2,'time','Dimensions',{'time',reconYrs});
%    ncwrite(outputf2,'time',r_o:r_f);
%    ncwriteatt(outputf2,'time','units','Annual mean starting from April to the next calendar year March');
%    ncwriteatt(outputf2,'time','long_name','Years CE');
% 
%    nccreate(outputf3,'time2','Dimensions',{'time2',length(p_yrs)});
%    ncwrite(outputf3,'time2',p_yrs);
%    ncwriteatt(outputf3,'time2','units','Annual mean starting from April to the next calendar year March');
%    ncwriteatt(outputf3,'time2','long_name','Model years');
% 
% 
%    % Extract the mean, prior, and ensembles
% 
%    X_mn0=Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);
%    X_mn=permute(reshape(X_mn0,length(xmeta{k}.lat),length(xmeta{k}.lon),reconYrs),[2 1 3]);
% 
%    X_xb0=Xb_o(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);
%    X_xb=permute(reshape(X_xb0,length(xmeta{k}.lat),length(xmeta{k}.lon),length(p_yrs)),[2 1 3]);
% 	  
%    X_ens0=Xa_ens(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:,:);
%    X_ens=permute(reshape(X_ens0,length(xmeta{k}.lat),length(xmeta{k}.lon),sub_ens,reconYrs),[2 1 4 3]);
%    clear X_ens0
% 
%    nccreate(outputf1,'tas_mn','Dimensions',{'lon',length(xmeta{k}.lon),'lat',length(xmeta{k}.lat),'time',reconYrs});
%    ncwrite(outputf1,'tas_mn',X_mn);
%    ncwriteatt(outputf1,'tas_mn','long_name',[xmeta{k}.varnmlng ' (reconstruction mean)']);
%    ncwriteatt(outputf1,'tas_mn','units','degC');
% 
% 
%    nccreate(outputf2,'tas_ens','Dimensions',{'lon',length(xmeta{k}.lon),'lat',length(xmeta{k}.lat),'time',reconYrs,'sub_ens',sub_ens});
%    ncwrite(outputf2,'tas_ens',X_ens);
%    ncwriteatt(outputf2,'tas_ens','long_name',[xmeta{k}.varnmlng ' (' num2str(sub_ens) ' member sub-ensemble)']);
%    ncwriteatt(outputf2,'tas_ens','units','degC');
% 
%    nccreate(outputf3,'tas_xb','Dimensions',{'lon',length(xmeta{k}.lon),'lat',length(xmeta{k}.lat),'time2',length(p_yrs)});
%    ncwrite(outputf3,'tas_xb',X_xb);
%    ncwriteatt(outputf3,'tas_xb','long_name',[xmeta{k}.varnmlng ' (reconstruction prior ensemble members)']);
%    ncwriteatt(outputf3,'tas_xb','units','degC');
% 
% 
% 
% end
% 
% 
% % SPEI
% if ismember('j',state_tp)
%    k=find(state_tp=='j');
% 
%    load(['/d1/nsteiger/climate-data/cesm-lme/cesm_lme_010_pdsi_output_biascorr_AWC_c_7_u2_gcm.mat'],'lndidx')
% 
%    % define the output variable
%    outputf1=['spei_mn_da_hydro_r.' num2str(r_o) '-' num2str(r_f) '_d.' timestp(1:end-9) '.nc'];
%    outputf2=['spei_ens_da_hydro_r.' num2str(r_o) '-' num2str(r_f) '_d.' timestp(1:end-9) '.nc'];
%    outputf3=['spei_prior_da_hydro_r.' num2str(r_o) '-' num2str(r_f) '_d.' timestp(1:end-9) '.nc'];
% 
%    % grab lat/lon from model, put it into the output file
%    eval(['!ncks -A -v lat,lon ',inputf,' ',outputf1]);
%    eval(['!ncks -A -v lat,lon ',inputf,' ',outputf2]);
%    eval(['!ncks -A -v lat,lon ',inputf,' ',outputf3]);
% 
%    nccreate(outputf1,'time','Dimensions',{'time',reconYrs});
%    ncwrite(outputf1,'time',r_o:r_f);
%    ncwriteatt(outputf1,'time','units','Annual mean starting from April to the next calendar year March');
%    ncwriteatt(outputf1,'time','long_name','Years CE');
% 
%    nccreate(outputf2,'time','Dimensions',{'time',reconYrs});
%    ncwrite(outputf2,'time',r_o:r_f);
%    ncwriteatt(outputf2,'time','units','Annual mean starting from April to the next calendar year March');
%    ncwriteatt(outputf2,'time','long_name','Years CE');
% 
%    nccreate(outputf3,'time2','Dimensions',{'time2',length(p_yrs)});
%    ncwrite(outputf3,'time2',p_yrs);
%    ncwriteatt(outputf3,'time2','units','Annual mean starting from April to the next calendar year March');
%    ncwriteatt(outputf3,'time2','long_name','Model years');
% 
% 
%    % Extract the mean, prior, and ensembles
% 
%    X_mn0=Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);
%    X_mn0f=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),reconYrs);
%    X_mn0f(lndidx,:)=X_mn0;
%    X_mn=permute(reshape(X_mn0f,length(xmeta{k}.lat),length(xmeta{k}.lon),reconYrs),[2 1 3]);
% 
%    X_xb0=Xb_o(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);
%    X_xb0f=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),length(p_yrs));
%    X_xb0f(lndidx,:)=X_xb0;
%    X_xb=permute(reshape(X_xb0f,length(xmeta{k}.lat),length(xmeta{k}.lon),length(p_yrs)),[2 1 3]);
% 	  
%    X_ens0=Xa_ens(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:,:);
%    X_ens0f=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),sub_ens,reconYrs);
%    X_ens0f(lndidx,:,:)=X_ens0;
%    clear X_ens0
%    X_ens=permute(reshape(X_ens0f,length(xmeta{k}.lat),length(xmeta{k}.lon),sub_ens,reconYrs),[2 1 4 3]);
%    clear X_ens0f
% 
%    nccreate(outputf1,'spei_mn','Dimensions',{'lon',length(xmeta{k}.lon),'lat',length(xmeta{k}.lat),'time',reconYrs},'FillValue',NaN);
%    ncwrite(outputf1,'spei_mn',X_mn);
%    ncwriteatt(outputf1,'spei_mn','long_name',[xmeta{k}.varnmlng ' (reconstruction mean)']);
%    ncwriteatt(outputf1,'spei_mn','units','degC');
%    ncwriteatt(outputf1,'spei_mn','missing_value',NaN);
% 
% 
%    nccreate(outputf2,'spei_ens','Dimensions',{'lon',length(xmeta{k}.lon),'lat',length(xmeta{k}.lat),'time',reconYrs,'sub_ens',sub_ens},'FillValue',NaN);
%    ncwrite(outputf2,'spei_ens',X_ens);
%    ncwriteatt(outputf2,'spei_ens','long_name',[xmeta{k}.varnmlng ' (' num2str(sub_ens) ' member sub-ensemble)']);
%    ncwriteatt(outputf2,'spei_ens','units','degC');
%    ncwriteatt(outputf2,'spei_ens','missing_value',NaN);
% 
%    nccreate(outputf3,'spei_xb','Dimensions',{'lon',length(xmeta{k}.lon),'lat',length(xmeta{k}.lat),'time2',length(p_yrs)},'FillValue',NaN);
%    ncwrite(outputf3,'spei_xb',X_xb);
%    ncwriteatt(outputf3,'spei_xb','long_name',[xmeta{k}.varnmlng ' (reconstruction prior ensemble members)']);
%    ncwriteatt(outputf3,'spei_xb','units','degC');
%    ncwriteatt(outputf3,'spei_xb','missing_value',NaN);
% 
% 
% 
% end














