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

%varpth='/d2/nsteiger/output-da/hydroclimate/';
varpth='/home/nsteiger/MATLAB/steiger-scidata-2018/output-recon/';

%efl='cesm_lme010_r10002000_p998_state_tj_avg_AprMar_prxydtst_5_prxtp_tca_2605_swtchbld100_30-Nov-2017_09:27:38.mat';
%efl='cesm_lme010_r10002000_p998_state_tj_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_16-Feb-2018_14:21:17.mat';

% ENSEMBLES FOR HUN'S PANCON DROUGHTS
%efl='cesm_lme010_r8002000_p998_state_tq_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_12-Jul-2018_14:50:03.mat';
%efl='cesm_lme010_r8002000_p998_state_tq_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_12-Jul-2018_14:59:32.mat';
efl='cesm_lme010_r8002000_p998_state_tq_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_06-Aug-2018_15:52:58.mat';

disp('Loading data...')

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

disp('Converting to netcdf...')


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



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

   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=['pdsi_mn_da_hydro_r.' num2str(r_o) '-' num2str(r_f) '_d.' timestp(1:end-9) '.nc'];
   outputf2=['pdsi_ens_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]);

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


   % 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_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,'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','degC');
   ncwriteatt(outputf1,'pdsi_mn','missing_value',NaN);


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




end











