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

varpth='./output-recon/';

efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:33:07.mat';
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:40:23.mat';
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:49:45.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');


   % Extract the mean and uncertainties

   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_s0=Xa_sigma(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);
   X_s=permute(reshape(X_s0,length(xmeta{k}.lat),length(xmeta{k}.lon),reconYrs),[2 1 3]);
	  
   X_p0=Xa_prct(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:,:);
   X_p=permute(reshape(X_p0,length(xmeta{k}.lat),length(xmeta{k}.lon),3,reconYrs),[2 1 4 3]);


%    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_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,'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(outputf1,'tas_sg','Dimensions',{'lon',length(xmeta{k}.lon),'lat',length(xmeta{k}.lat),'time',reconYrs});
   ncwrite(outputf1,'tas_sg',X_s);
   ncwriteatt(outputf1,'tas_sg','long_name',[xmeta{k}.varnmlng ' (1 standard deviation of ensemble)']);
   ncwriteatt(outputf1,'tas_sg','units','degC');

   nccreate(outputf1,'tas_pc','Dimensions',{'lon',length(xmeta{k}.lon),'lat',length(xmeta{k}.lat),'time',reconYrs,'prctl',3});
   ncwrite(outputf1,'tas_pc',X_p);
   ncwriteatt(outputf1,'tas_pc','long_name',[xmeta{k}.varnmlng ' (ensemble percentiles: 5th, 50th, 95th)']);
   ncwriteatt(outputf1,'tas_pc','units','degC');

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


end


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


   % Extract the mean and uncertainties
   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(xmeta{k}.lndidx,:)=X_mn0;
   X_mn=permute(reshape(X_mn0f,length(xmeta{k}.lat),length(xmeta{k}.lon),reconYrs),[2 1 3]);

   X_s0=Xa_sigma(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);
   X_s0f=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),reconYrs);
   X_s0f(xmeta{k}.lndidx,:)=X_s0;
   X_s=permute(reshape(X_s0f,length(xmeta{k}.lat),length(xmeta{k}.lon),reconYrs),[2 1 3]);

   X_p0=Xa_prct(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:,:);
   X_p0f=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),3,reconYrs);
   X_p0f(xmeta{k}.lndidx,:,:)=X_p0;
   X_p=permute(reshape(X_p0f,length(xmeta{k}.lat),length(xmeta{k}.lon),3,reconYrs),[2 1 4 3]);


   nccreate(outputf1,'pdsi_mn','Dimensions',{'lon',length(xmeta{k}.lon),'lat',length(xmeta{k}.lat),'time',reconYrs},'FillValue',NaN);
   ncwrite(outputf1,'pdsi_mn',X_mn);
   ncwriteatt(outputf1,'pdsi_mn','long_name',[xmeta{k}.varnmlng ' (reconstruction mean)']);
   ncwriteatt(outputf1,'pdsi_mn','units','Standardized units of relative dry and wet');
   ncwriteatt(outputf1,'pdsi_mn','missing_value',NaN);

   nccreate(outputf1,'pdsi_sg','Dimensions',{'lon',length(xmeta{k}.lon),'lat',length(xmeta{k}.lat),'time',reconYrs},'FillValue',NaN);
   ncwrite(outputf1,'pdsi_sg',X_s);
   ncwriteatt(outputf1,'pdsi_sg','long_name',[xmeta{k}.varnmlng ' (1 standard deviation of ensemble)']);
   ncwriteatt(outputf1,'pdsi_sg','units','Standardized units of relative dry and wet');
   ncwriteatt(outputf1,'pdsi_sg','missing_value',NaN);

   nccreate(outputf1,'pdsi_pc','Dimensions',{'lon',length(xmeta{k}.lon),'lat',length(xmeta{k}.lat),'time',reconYrs,'prctl',3},'FillValue',NaN);
   ncwrite(outputf1,'pdsi_pc',X_p);
   ncwriteatt(outputf1,'pdsi_pc','long_name',[xmeta{k}.varnmlng ' (ensemble percentiles: 5th, 50th, 95th)']);
   ncwriteatt(outputf1,'pdsi_pc','units','Standardized units of relative dry and wet');
   ncwriteatt(outputf1,'pdsi_pc','missing_value',NaN);




end


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


   % Extract the mean and uncertainties
   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(xmeta{k}.lndidx,:)=X_mn0;
   X_mn=permute(reshape(X_mn0f,length(xmeta{k}.lat),length(xmeta{k}.lon),reconYrs),[2 1 3]);

   X_s0=Xa_sigma(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);
   X_s0f=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),reconYrs);
   X_s0f(xmeta{k}.lndidx,:)=X_s0;
   X_s=permute(reshape(X_s0f,length(xmeta{k}.lat),length(xmeta{k}.lon),reconYrs),[2 1 3]);

   X_p0=Xa_prct(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:,:);
   X_p0f=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),3,reconYrs);
   X_p0f(xmeta{k}.lndidx,:,:)=X_p0;
   X_p=permute(reshape(X_p0f,length(xmeta{k}.lat),length(xmeta{k}.lon),3,reconYrs),[2 1 4 3]);


   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','Standardized units of relative dry and wet');
   ncwriteatt(outputf1,'spei_mn','missing_value',NaN);

   nccreate(outputf1,'spei_sg','Dimensions',{'lon',length(xmeta{k}.lon),'lat',length(xmeta{k}.lat),'time',reconYrs},'FillValue',NaN);
   ncwrite(outputf1,'spei_sg',X_s);
   ncwriteatt(outputf1,'spei_sg','long_name',[xmeta{k}.varnmlng ' (1 standard deviation of ensemble)']);
   ncwriteatt(outputf1,'spei_sg','units','Standardized units of relative dry and wet');
   ncwriteatt(outputf1,'spei_sg','missing_value',NaN);

   nccreate(outputf1,'spei_pc','Dimensions',{'lon',length(xmeta{k}.lon),'lat',length(xmeta{k}.lat),'time',reconYrs,'prctl',3},'FillValue',NaN);
   ncwrite(outputf1,'spei_pc',X_p);
   ncwriteatt(outputf1,'spei_pc','long_name',[xmeta{k}.varnmlng ' (ensemble percentiles: 5th, 50th, 95th)']);
   ncwriteatt(outputf1,'spei_pc','units','Standardized units of relative dry and wet');
   ncwriteatt(outputf1,'spei_pc','missing_value',NaN);

end




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


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


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

      X_mn=reshape(Xa_m((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:),12*reconYrs,1);
      X_s=reshape(Xa_sigma((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:),12*reconYrs,1);
      X_p=reshape(permute(Xa_prct((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:,:),[1 3 2]),12*reconYrs,3);

      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 '_mn'],'Dimensions',{'tmon',length(mon_r)});
      ncwrite(outputf1,[vrnm0 '_mn'],X_mn);
      ncwriteatt(outputf1,[vrnm0 '_mn'],'long_name',[vlngnm ' (reconstruction mean)']);
      ncwriteatt(outputf1,[vrnm0 '_mn'],'units','degC');

      nccreate(outputf1,[vrnm0 '_sg'],'Dimensions',{'tmon',length(mon_r)});
      ncwrite(outputf1,[vrnm0 '_sg'],X_s);
      ncwriteatt(outputf1,[vrnm0 '_sg'],'long_name',[vlngnm ' (1 standard deviation of ensemble)']);
      ncwriteatt(outputf1,[vrnm0 '_sg'],'units','degC');

      nccreate(outputf1,[vrnm0 '_pc'],'Dimensions',{'tmon',length(mon_r),'prctl',3});
      ncwrite(outputf1,[vrnm0 '_pc'],X_p);
      ncwriteatt(outputf1,[vrnm0 '_pc'],'long_name',[vlngnm ' (ensemble percentiles: 5th, 50th, 95th)']);
      ncwriteatt(outputf1,[vrnm0 '_pc'],'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)

      X_mn=reshape(Xa_m(xmeta{k}.id_X(1)+i-1:xmeta{k}.id_X(1)+i-1,:),reconYrs,1);
      X_s=reshape(Xa_sigma(xmeta{k}.id_X(1)+i-1:xmeta{k}.id_X(1)+i-1,:),reconYrs,1);
      X_p=permute(squeeze(Xa_prct(xmeta{k}.id_X(1)+i-1:xmeta{k}.id_X(1)+i-1,:,:)),[2 1]);


      nccreate(outputf1,[vrnm0{i} '_mn'],'Dimensions',{'time',reconYrs});
      ncwrite(outputf1,[vrnm0{i} '_mn'],X_mn);
      ncwriteatt(outputf1,[vrnm0{i} '_mn'],'long_name',['ITCZ ' xmeta{k}.indinfo{i}.varnmlng ' (reconstruction mean)']);
      ncwriteatt(outputf1,[vrnm0{i} '_mn'],'units','lat (deg)');

      nccreate(outputf1,[vrnm0{i} '_sg'],'Dimensions',{'time',reconYrs});
      ncwrite(outputf1,[vrnm0{i} '_sg'],X_s);
      ncwriteatt(outputf1,[vrnm0{i} '_sg'],'long_name',['ITCZ ' xmeta{k}.indinfo{i}.varnmlng ' (1 standard deviation of ensemble)']);
      ncwriteatt(outputf1,[vrnm0{i} '_sg'],'units','lat (deg)');

      nccreate(outputf1,[vrnm0{i} '_pc'],'Dimensions',{'time',reconYrs,'prctl',3});
      ncwrite(outputf1,[vrnm0{i} '_pc'],X_p);
      ncwriteatt(outputf1,[vrnm0{i} '_pc'],'long_name',['ITCZ ' xmeta{k}.indinfo{i}.varnmlng ' (ensemble percentiles: 5th, 50th, 95th)']);
      ncwriteatt(outputf1,[vrnm0{i} '_pc'],'units','lat (deg)');

    
   end


end


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

   X_mn=reshape(Xa_m(xmeta{k}.id_X(1),:),reconYrs,1);
   X_s=reshape(Xa_sigma(xmeta{k}.id_X(1),:),reconYrs,1);
   X_p=permute(squeeze(Xa_prct(xmeta{k}.id_X(1),:,:)),[2 1]);


   nccreate(outputf1,['gmt_mn'],'Dimensions',{'time',reconYrs});
   ncwrite(outputf1,['gmt_mn'],X_mn);
   ncwriteatt(outputf1,['gmt_mn'],'long_name',[xmeta{k}.varnmlng ' (reconstruction mean)']);
   ncwriteatt(outputf1,['gmt_mn'],'units','degC');

   nccreate(outputf1,['gmt_sg'],'Dimensions',{'time',reconYrs});
   ncwrite(outputf1,['gmt_sg'],X_s);
   ncwriteatt(outputf1,['gmt_sg'],'long_name',[xmeta{k}.varnmlng ' (1 standard deviation of ensemble)']);
   ncwriteatt(outputf1,['gmt_sg'],'units','degC');

   nccreate(outputf1,['gmt_pc'],'Dimensions',{'time',reconYrs,'prctl',3});
   ncwrite(outputf1,['gmt_pc'],X_p);
   ncwriteatt(outputf1,['gmt_pc'],'long_name',[xmeta{k}.varnmlng ' (ensemble percentiles: 5th, 50th, 95th)']);
   ncwriteatt(outputf1,['gmt_pc'],'units','degC');

end


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

   X_mn=reshape(Xa_m(xmeta{k}.id_X(1),:),reconYrs,1);
   X_s=reshape(Xa_sigma(xmeta{k}.id_X(1),:),reconYrs,1);
   X_p=permute(squeeze(Xa_prct(xmeta{k}.id_X(1),:,:)),[2 1]);


   nccreate(outputf1,['amo_mn'],'Dimensions',{'time',reconYrs});
   ncwrite(outputf1,['amo_mn'],X_mn);
   ncwriteatt(outputf1,['amo_mn'],'long_name',[xmeta{k}.varnmlng ' (reconstruction mean)']);
   ncwriteatt(outputf1,['amo_mn'],'units','degC');

   nccreate(outputf1,['amo_sg'],'Dimensions',{'time',reconYrs});
   ncwrite(outputf1,['amo_sg'],X_s);
   ncwriteatt(outputf1,['amo_sg'],'long_name',[xmeta{k}.varnmlng ' (1 standard deviation of ensemble)']);
   ncwriteatt(outputf1,['amo_sg'],'units','degC');

   nccreate(outputf1,['amo_pc'],'Dimensions',{'time',reconYrs,'prctl',3});
   ncwrite(outputf1,['amo_pc'],X_p);
   ncwriteatt(outputf1,['amo_pc'],'long_name',[xmeta{k}.varnmlng ' (ensemble percentiles: 5th, 50th, 95th)']);
   ncwriteatt(outputf1,['amo_pc'],'units','degC');

end



%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, Mar 2018" ' outputf1]);
eval(['!ncatted -O --glb_att_add citation="https://doi.org/10.5281/zenodo.1154913" ' 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














