function [u,ibu,jbu,kbu,nxul,nyul,nzul,...
    v,ibv,jbv,kbv,nxvl,nyvl,nzvl,...
    w,ibw,jbw,kbw,nxwl,nywl,nzwl,...
    p,ibp,jbp,kbp,nxpl,nypl,nzpl,ighost] = read_uvwp_chunk_legacy(fuvwp,fproc,ighost,varargin)
    % [u,ibu,jbu,kbu,nxul,nyul,nzul,...
    %  v,ibv,jbv,kbv,nxvl,nyvl,nzvl,...
    %  w,ibw,jbw,kbw,nxwl,nywl,nzwl,...
    %  p,ibp,jbp,kbp,nxpl,nypl,nzpl,ighost] = read_uvwp_chunk_legacy(fuvwp,fproc,ighost,varargin)
    % Reads u,v,w,p from one processor chunk.
    % Input
    %   fuvwp               path to uvwp chunk
    %   fproc               path to 'current_proc_grid.tmp'
    %   ighost              data written with ghost cells?
    % ? ghost               keep ghost cells? (default: 1)
    % ? verbosity           verbose output? (default: 0)
    % ? nparam              number of parameters in header (default: 10)
    % ? precision           precision of data (default: 'float64')
    %                        'float32'
    %                        'float64'
    % ? endian              endianess of file (default: 'a')
    % ? reclen              FORTRAN record length in bytes (default: 4)
    % Output
    %   u,v,w,p             partial fields with dim(nxl,nyl,nzl)+2*ighost
    %   ibu,jbu,kbu,...     global index of first grid point of the partial field w/o ghost cells
    %   nxul,nyul,nzul,...  local field size w/o ghost cells
    %   ighost              ghost cell flag

    % Parse optional input arguments
    par = inputParser;
    addParamValue(par,'ghost',1,@isnumeric);
    addParamValue(par,'verbosity',0,@isnumeric);
    addParamValue(par,'nparam',10,@isnumeric);
    addParamValue(par,'precision','float64',@ischar);
    addParamValue(par,'endian','a',@ischar);
    addParamValue(par,'reclen',4,@isnumeric);
    parse(par,varargin{:});
    keepghost = par.Results.ghost;
    verbosity = par.Results.verbosity;
    nparam    = par.Results.nparam;
    precision = par.Results.precision;
    endian    = par.Results.endian;
    reclen    = par.Results.reclen;
    
    % First read processor grid from tmp file
    [ibegu,iendu,jbegu,jendu,kbegu,kendu,...
     ibegv,iendv,jbegv,jendv,kbegv,kendv,...
     ibegw,iendw,jbegw,jendw,kbegw,kendw,...
     ibegp,iendp,jbegp,jendp,kbegp,kendp,...
     nxprocs,nyprocs,nzprocs] = read_procgrid_legacy(fproc);

    % Determine processor rank from filename
    [idxbeg,idxend] = regexp(fuvwp,'[\d+]$');
    if isempty(idxbeg)
      error('Invalid file: does not contain rank. %s',fuvwp);
    end
    iproc = str2double(fuvwp(idxbeg:idxend));
    
    % Determine local array sizes
    ixproc = floor(iproc/(nyprocs*nzprocs));
    iyproc = mod(floor(iproc/nzprocs),nyprocs);
    izproc = mod(iproc,nzprocs);
    
    ibu  = ibegu(ixproc+1);
    jbu  = jbegu(iyproc+1);
    kbu  = kbegu(izproc+1);
    nxul = iendu(ixproc+1)-ibu+1;
    nyul = jendu(iyproc+1)-jbu+1;
    nzul = kendu(izproc+1)-kbu+1;
    
    ibv  = ibegv(ixproc+1);
    jbv  = jbegv(iyproc+1);
    kbv  = kbegv(izproc+1);
    nxvl = iendv(ixproc+1)-ibv+1;
    nyvl = jendv(iyproc+1)-jbv+1;
    nzvl = kendv(izproc+1)-kbv+1;
    
    ibw  = ibegw(ixproc+1);
    jbw  = jbegw(iyproc+1);
    kbw  = kbegw(izproc+1);
    nxwl = iendw(ixproc+1)-ibw+1;
    nywl = jendw(iyproc+1)-jbw+1;
    nzwl = kendw(izproc+1)-kbw+1;
    
    ibp  = ibegp(ixproc+1);
    jbp  = jbegp(iyproc+1);
    kbp  = kbegp(izproc+1);
    nxpl = iendp(ixproc+1)-ibp+1;
    nypl = jendp(iyproc+1)-jbp+1;
    nzpl = kendp(izproc+1)-kbp+1;
    
    % Convert 'precision' into bytes
    switch precision
      case 'float32'; nprecision = 4;
      case 'float64'; nprecision = 8;
      otherwise; error('Invalid precision: %s',precision)
    end
    
    % Open uvwp file
    fid = fopen(fuvwp,'r',endian);
    if fid<0
      error('File not found: %s',fuvwp);
    end
    
    % Skip headers and read fields
    nheader = 4*reclen+(4+nparam)*nprecision;
    fseek(fid,nheader+reclen,'bof');
    nx = (nxul+2*ighost);
    ny = (nyul+2*ighost);
    nz = (nzul+2*ighost);
    u = fread(fid,nx*ny*nz,precision);
    u = reshape(u,nx,ny,nz);
    fseek(fid,nheader+2*reclen,'cof');
    nx = (nxvl+2*ighost);
    ny = (nyvl+2*ighost);
    nz = (nzvl+2*ighost);
    v = fread(fid,nx*ny*nz,precision);
    v = reshape(v,nx,ny,nz);
    fseek(fid,nheader+2*reclen,'cof');
    nx = (nxwl+2*ighost);
    ny = (nywl+2*ighost);
    nz = (nzwl+2*ighost);
    w = fread(fid,nx*ny*nz,precision);
    w = reshape(w,nx,ny,nz);
    fseek(fid,nheader+2*reclen,'cof');
    nx = (nxpl+2*ighost);
    ny = (nypl+2*ighost);
    nz = (nzpl+2*ighost);
    p = fread(fid,nx*ny*nz,precision);
    p = reshape(p,nx,ny,nz);
    fseek(fid,reclen,'cof');
    
    % Test for eof
    fread(fid,1,'char');
    if ~feof(fid)
      warning('End-of-file not reached: data might be corrupt.');
    end
    
    % Close file
    fclose(fid);
        
    % Remove ghosts if necessary
    if ighost && ~keepghost
      u = u(2:end-1,2:end-1,2:end-1);
      v = v(2:end-1,2:end-1,2:end-1);
      w = w(2:end-1,2:end-1,2:end-1);
      p = p(2:end-1,2:end-1,2:end-1);
      ighost = 0;
    end
end