
% PARAMETERS FOR DATA ASSIMILATION BASED RECONSTRUCTIONS

%rng('shuffle') % Shuffle random number seed

%=================================
% Choose the model for the prior
%=================================

% CESM option: 002 through 010
mdl_prior='cesm_lme002';

disp(['Prior model = ' mdl_prior])

%----------------------------------------
% Base time scale for the reconstruction
%----------------------------------------
% Jan=1, Feb=2, Mar=3, Apr=4, May=5, Jun=6
% Jul=7, Aug=8, Sep=9, Oct=10, Nov=11, Dec=12

% Annual average (respecting the climate system, not the calendar)
%mon_avg_o=4;mon_avg_f=3;
%mon_o='Apr';mon_f='Mar';

% DJF seasonal average
%mon_avg_o=12;mon_avg_f=2; 
%mon_o='Dec';mon_f='Feb'; 

% JJA seasonal average
mon_avg_o=6;mon_avg_f=8; 
mon_o='Jun';mon_f='Aug'; 

disp(['Seasonality = ' mon_o '--' mon_f])

% 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 = 't2qjozn';
%state_tp = 'q';
%state_tp = '2oznw';
%state_tp = '2ozn';
disp(['State type = ' state_tp ])

% Get model information
[gcm_syr,gcm_eyr]=modelyrs(mdl_prior);

% Number of reconstruction iterations
%itr=5;
%itr=25;
itr=1;
%disp(['MC iterations = ' num2str(itr)])

% Switchblade proxies? y/n
%swtchbld=30;
%swtchbld=50;
%swtchbld=75;
swtchbld=100; % use all proxies but random ordering
disp(['Switchblade ' num2str(swtchbld) '% of proxies.'])

% Percentiles from posterior
xa_prct=[5 50 95];
disp(['Will save ' num2str(xa_prct) ' percentiles' ])

% Save analysis ensemble?
%ens_save= 'a'; % Save all posterior ensembles: ONLY USE THIS FOR INDEX RECONS!!!
%ens_save= 'y';
ens_save= 'n';
%ens_save='s'; sub_ens=100; % Save subsample of posterior
disp(['Save analysis ensemble? ' ens_save])

% Save output?  y/n
outsave = 'y';
%outsave = 'n';
disp(['Save output? ' outsave])
% Set location of output
outpath='/d2/nsteiger/output-da/hydroclimate/';
%outpath='/d2/nsteiger/output-da/pages2kv2/';

if outsave=='y'; disp(['Output will be saved to ' outpath]); end

% USE THIS OPTION FOR ALL CESM_LME REAL RECONSTRUCTIONS
if strcmp(mdl_prior(1:8),'cesm_lme') 
    
    %rng('default');disp('Using consistent random prior...')
    %rng('shuffle');disp('Using different random prior...')
    %p_yrs=randsample((gcm_syr:gcm_eyr-1)-gcm_syr+1,900); % rand sample
    
    % almost all 1000 years (can't use the beg/end years b/c of year definition and var quirks)
    p_yrs=[851:1848]-gcm_syr+1; 
    
    % Years of the reconstruction and calibration period
    r_o = p_yrs(1)+gcm_syr-1;r_f = p_yrs(end)+gcm_syr-1; clb_prd=[1920 1999]; % do a verif of 1870:1919 (50 yrs)
    %r_o = 1;r_f = 2000; clb_prd=[1920 1999]; % do a verif of 1870:1919 (50 yrs)
    %r_o = 1;r_f = 2000; clb_prd=[1911 1995];disp('Calibration period for PAGES2k recons...')
    %r_o = 1000;r_f = 2000; clb_prd=[1920 1999]; % do a verif of 1870:1919 (50 yrs)
    disp(['PSM calib period = ' num2str(clb_prd)]) 
    
    reconYrs=length(r_o:r_f);
else
    error('Must add more options for prior type')
end

disp(['Prior size = ' num2str(length(p_yrs))])

disp(['Reconstruction years = ' num2str(r_o) '-' num2str(r_f) ])

disp('Loading data...')

% LOAD THE MODEL DATA
%[Xb_o,xmeta]=load_S(mdl_prior,state_tp,mon_avg_o,mon_avg_f,p_yrs);
[Xa_m,xmeta]=load_S_nobc(mdl_prior,state_tp,mon_avg_o,mon_avg_f,p_yrs);

% State vector length
svl=size(Xa_m,1);

if ens_save=='a' && svl > 100
   error('State vector is too large to save all posterior ensembles!')
elseif ens_save=='s' && svl*sub_ens*reconYrs*8 > 2e10
   error('State vector will be over 20 GB! Cannot save posterior ensembles.')
end 

% Inflation/deflation factor (if = 1, prior is not inflated)
infl = 1; if infl~=1; disp(['Inflation = ' num2str(infl)]); end

% file name corresponding to current time stamp
timestp=datestr(now);
timestp=strrep(timestp,' ','_'); % remove space

%fln=['da_recon_' fdata(1:end-4) '.mat'];

fln=[mdl_prior '_r' num2str(r_o) num2str(r_f) '_p' num2str(length(p_yrs))...
  '_state_' state_tp '_avg_' mon_o mon_f '_' timestp '.mat'];


disp(['outfile = ' fln])







% 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='cesm_latlon.nc';



% define the output variable
outputf1=[mdl_prior '_' mon_o mon_f '_d.' num2str(r_o) '-' num2str(r_f) '.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]);









