function [field time_step] = nctiles2aste_v2(varName,rootDir,tile_list,levels_list,times_list,tiles,flag_aste);   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Read specified netCDF tiles from disk, given selected levels and time steps, and return a compact format aste field
%This script requires the output of aste2tiles, a structure named tiles
%
%Input:
%	varName		Name of the diagnostic, e.g. in the files THETA.xxxx.data in binary tile format varName='THETA'
%	rootDir		Directory where diag tiles live (full path)
%	tile_list	Vector with a list of the tiles (tile number) being extracted, e.g. [8:16] selects the tiles that make up face 3 of ASTE
%	levels_list	Vector of vertical levels, e.g. [1:nz]
%	times_list	Time step ordinal being stracted, e.g. 1 <=times<= 12 for climatology or 1 <=times<= 192 for ASTE release 1 
%	tiles		Structure with fields tiles.tot,tiles.wet,...,tiles.index 
%			This is the output of aste2tiles.m
%       flag_aste       0 = leave outputs as individual tiles, save memory if want to analyze only a few tiles
%                       1 = put into aste compact format [270 1350 [1 or nz]] , for use if analyze the entire ASTE domain
%
%Output:
%	field		individual tile or ASTE compact format field with dimensions [270 1350 length(levels_list) length(times_list)]
%	time_step	Time step number corresponding to the times selected from the tiles 
%			For climatologies, first and last timesteps used top calculate thwe climatology
%
% V.Ocana 15.Nov.2019
% atn 01.Nov.2020: modify to output as double arrays instead of cell structures (flag_aste=0)
%                         OR output is for a individual tiles (flag_aste=1);
% atn 03.Feb.2021: modify to accomodate reading 2D or 3D fields (add if loop based on nc_info.dims)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%keyboard
%Get info from meta file
dataDir =[rootDir varName '/'];
fList  = dir([dataDir '*.nc']);
nc_info=read_nc_info(fList(1).name,dataDir);

varindex = find(ismember(nc_info.varName, varName));
varid = nc_info.varID(varindex);
timindex = find(ismember(nc_info.varName, 'timstep'));
timid = nc_info.varID(timindex);

%Set dimensions
nx = tiles.fsize(1);
ny = tiles.fsize(2);
nz=1; if (nc_info.ndims==5); nz= nc_info.dimSize(3); end;
dtilex = tiles.tsize(1);
dtiley = tiles.tsize(2);
fieldSize=[nx ny nz];

ntimes_max   = nc_info.dimSize(2); 
ntiles_max = length(fList);

ntiles  = length(tile_list); 
nlevels = length(levels_list); 
ntimes  = length(times_list);

%Check dimensions
if(dtilex~=nc_info.dimSize(end) || dtiley~=nc_info.dimSize(end-1));
   fprintf('nctiles2aste ERROR: Dimensions from tiles.tsize %dx%d are different from meta.dimList %dx%d',dtilex,dtiley,meta.dimList(1),meta.dimList(4));
   return;
end;
%Check that tiles, levels and times requested are available
if(ntiles_max~=tiles.wet)
   disp(['nctiles2aste WARNING: number of tiles in directory ' num2str(ntiles_max) ' is different from number of tiles from ASTE2tiles = ' num2str(tiles.wet)]);
end;
if(length(tile_list)>tiles.wet);
   disp(['nctiles2aste ERROR: Number of wet tiles in dataDir is ' num2str(tiles.wet) ' but you are requesting ' num2str(length(tile_list)) ' tiles']); 
   return;
end;
if(max(tile_list)>tiles.wet);
   disp(['nctiles2aste ERROR: Number of wet tiles in dataDir is ' num2str(tiles.wet) ' but you are requesting a tile number' num2str(max(tile_list))]); 
   return;
end;
if(nlevels>nz);
   disp(['nctiles2aste ERROR: Number of vertical levels in dataDir is ' num2str(nz) ' but you are requesting ' num2str(nlevels) ' levels']); 
   return;
end;
if(max(levels_list)>nz);
   disp(['nctiles2aste ERROR: Number of vertical levels in dataDir is ' num2str(nz) ' but you are requesting a level number ' num2str(max(levels_list))]); 
   return;
end;
if(max(times_list)>ntimes_max);
   disp(['nctiles2aste ERROR: Number of time steps in dataDir is ' num2str(ntimes_max) ' but you are requesting time step ' num2str(max(times_list))]); 
   return;
end;

%Extract the data from tiles (and re-assemble in compact form if flag_aste=1)
field_tiles=[];
for itile = 1:ntiles;
   tmp=nan(dtilex,dtiley,nlevels,ntimes);	%put into an array instead of cell structure
   tile_num = tile_list(itile);
   indices{itile} = tiles.index{tile_num};
   fName = [varName '.' sprintf('%04d',tile_num) '.nc'];
   ncid = netcdf.open([dataDir fName],'NC_NOWRITE');
   for itime=1:ntimes;
      time=times_list(itime);
      if(itile==1); time_step(itime,1) = netcdf.getVar(ncid,timid,[times_list(itime)-1],[1]);end;
      for ilev=1:nlevels;
         level=levels_list(ilev);
%for 2D fields we need to reduce dim of start and count compared to for 3D fields
         if(nc_info.ndims==4);
           start = [0 0 (time-1)];
           count = [dtilex dtiley 1];
         elseif(nc_info.ndims==5);
           start = [0 0 (level-1) (time-1)];
           count = [dtilex dtiley 1 1];
         end;
         tmp(:,:,ilev,itime) =netcdf.getVar(ncid,varid,start,count);
      end;
   end;
   field_tiles{itile} = tmp;
   netcdf.close(ncid);
   clear tmp;
end;


%note that this is memory intensive
%fprintf('flag_aste: %i\n',flag_aste);
if(flag_aste);
   fprintf('WARNING: putting fields into full ASTE compact format, memory intensive\n');
   field=nan(tiles.fsize(1),tiles.fsize(2),nlevels,ntimes);

   for itime=1:ntimes;
     for itile=1:ntiles
       tmp{itile}=field_tiles{itile}(:,:,:,itime);
     end;
     field(:,:,:,itime) = tiles2field(tmp,indices,fieldSize);
   end;
else;
   field=field_tiles;
end;
