function [q,x,y,z] = read_field_complete_ucfmf(fdir,iseq,field,varargin)
    % [q,x,y,z] = read_field_complete_ucfmf(file,field,varargin)
    % Reads a specific field from all processor chunks.(Multi-file tar version)
    % Input
    %   fdir                directory with files
    %   iseq                sequence index
    %   field               field to be read 
    %                         {'u','v','w','p','s1','s2',...}
    % ? step                index of step to be read (default: 1)
    % ? verbosity           verbose output? (default: 0)
    % ? debug               debug output? (default: 0)
    % Output
    %   q                   complete field
    %   x,y,z               corresponding grid

    % Parse optional input arguments
    par = inputParser;
    addParamValue(par,'step',1,@isnumeric);
    addParamValue(par,'verbosity',0,@isnumeric);
    addParamValue(par,'debug',0,@isnumeric);
    parse(par,varargin{:});
    istep = par.Results.timestep;

    % Parse field
    if ischar(field)
      switch field(1)
      case 'u'; ifield=1; fbase='uvwp'; cmesh='u';
      case 'v'; ifield=2; fbase='uvwp'; cmesh='v';
      case 'w'; ifield=3; fbase='uvwp'; cmesh='w';
      case 'p'; ifield=4; fbase='uvwp'; cmesh='p';
      case 's'; ifield=str2double(field(2:end)); fbase='scal'; cmesh='p';
      otherwise; error('Invalid field: %s',field);
    end
    
    % Read parameters
    fparam = sprintf('%s/parameters_%04d.asc',fdir,iseq);
    params = read_parameters_ucf(fparam,'tarmode',0);
    
    % Construct an array with final mesh size
    nx = params.mesh.(sprintf('nx%s',cmesh));
    ny = params.mesh.(sprintf('ny%s',cmesh));
    nz = params.mesh.(sprintf('nz%s',cmesh));
    q = zeros(nx,ny,nz);
    
    % Read grid
    fgrid = sprintf('%s/grid_%04d.bin',fdir,iseq);
    [x.u,y.u,z.u,x.v,y.v,z.v,x.w,y.w,z.w,x.p,y.p,z.p] = read_grid_ucf(fgrid,'tarmode',0);
    x = x.(cmesh);
    y = y.(cmesh);
    z = z.(cmesh);
    
    % Get number of processors
    nprocs = params.parallel.nprocs;
    
    % Now read chunk by chunk
    for iproc=0:nprocs-1
      fq = sprintf('%s/%s_%04d.%05d',fpath,iseq,iproc);
      [data,ib,jb,kb,nxl,nyl,nzl] = read_field_chunk_ucf(fq,ifield,'tarmode',0,'ghost',0,'step',istep);
      q(ib:ib+nxl-1,jb:jb+nyl-1,kb:kb+nzl-1) = data;
    end
end