diff --git a/matlab/read_field_chunk_legacy.m b/matlab/read_field_chunk_legacy.m
new file mode 100644
index 0000000..563d26a
--- /dev/null
+++ b/matlab/read_field_chunk_legacy.m
@@ -0,0 +1,157 @@
+function [data,ib,jb,kb,nxl,nyl,nzl,ighost] = read_field_chunk_legacy(file,fproc,ighost,field,varargin)
+    % [data,ib,jb,kb,nxl,nyl,nzl,ighost] = read_field_chunk_legacy(file,fproc,ighost,field,varargin)
+    % Reads arbitrary field from one processor chunk.
+    % Input
+    %   file                path to chunk
+    %   fproc               path to 'current_proc_grid.tmp'
+    %   ighost              data written with ghost cells?
+    %   field               field to be read 
+    %                         {'u','v','w','p','s1','s2',...}
+    % ? 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
+    %   data                partial field with dim(nxl+2*ighost,nyl+2*ighost,nzl+2*ighost)
+    %   ib,jb,kb            global index of first grid point of the partial field w/o ghost cells
+    %   nxl,nyl,nzl         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;
+    
+    % Parse field
+    if ischar(field)
+      switch field(1)
+      case 'u'; ifield=1; fbase='uvwp';
+      case 'v'; ifield=2; fbase='uvwp';
+      case 'w'; ifield=3; fbase='uvwp';
+      case 'p'; ifield=4; fbase='uvwp';
+      case 's'; ifield=str2double(field(2:end)); fbase='scal';
+      end
+    else
+      error('field must be of type char');
+    end
+    
+    % 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(file,'[\d]+$');
+    if isempty(idxbeg)
+      error('Invalid file: does not contain rank. %s',fuvwp);
+    end
+    iproc = str2double(file(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
+    
+    % Determine header size
+    nheader = 4*reclen+(4+nparam)*nprecision;
+    
+    % Skip datasets until correct one is reached
+    for iset=1:ifield
+      switch fbase
+      case 'uvwp'
+        if iset==1; nxl=nxul; nyl=nyul; nzl=nzul; end
+        if iset==2; nxl=nxvl; nyl=nyvl; nzl=nzvl; end
+        if iset==3; nxl=nxwl; nyl=nywl; nzl=nzwl; end
+        if iset==4; nxl=nxpl; nyl=nypl; nzl=nzpl; end
+      case 'scal'
+        nxl=nxpl; nyl=nypl; nzl=nzpl;
+      end
+      ndata = (nxl+2*ighost)*(nyl+2*ighost)*(nzl+2*ighost);
+      if iset~=ifield
+        nskip = nheader+2*reclen+ndata*nprecision;
+        fseek(fid,nskip,'cof');
+      else
+        fseek(fid,nheader+reclen,'cof');
+        data = fread(fid,ndata,precision);
+        data = reshape(data,nxl,nyl,nzl);
+      end
+    end
+      
+    % Get processor bounds
+    switch fbase
+    case 'uvwp'
+      if ifield==1; ib=ibu; jb=jbu; kb=kbu; nxl=nxul; nyl=nyul; nzl=nzul; end
+      if ifield==2; ib=ibv; jb=jbv; kb=kbv; nxl=nxvl; nyl=nyvl; nzl=nzvl; end
+      if ifield==3; ib=ibw; jb=jbw; kb=kbw; nxl=nxwl; nyl=nywl; nzl=nzwl; end
+      if ifield==4; ib=ibp; jb=jbp; kb=kbp; nxl=nxpl; nyl=nypl; nzl=nzpl; end
+    case 'scal'
+      ib=ibp; jb=jbp; kb=kbp; nxl=nxpl; nyl=nypl; nzl=nzpl;
+    end
+    
+    % Close file
+    fclose(fid);
+        
+    % Remove ghosts if necessary
+    if ighost && ~keepghost
+      data = data(2:end-1,2:end-1,2:end-1);
+      ighost = 0;
+    end
+end