function [] = interp_uvwp_ucf(hucf,dout,nxpn,nypn,nzpn,nxprocsn,nyprocsn,nzprocsn,varargin)
  % [] = interp_uvwp_ucf(hucf,dout,nxpn,nypn,nzpn,nxprocsn,nyprocsn,nzprocsn,varargin)
  % Interpolates a given flow field onto a new grid. 
  % This function reads every chunk of the original field only once, while
  % trying to minimize the memory footprint by writing the new chunks as
  % soon as possible to disk. Memory for the new chunks is only allocated if
  % they contain partial data. 
  % Input
  %   hucf              UCF handle (ustar,ucfmulti)
  %   dout              output directory
  %   nxpn,...          number of grid points of the new pressure grid
  %   nxprocs,...       number of processors of the new processor grid
  % ? verbosity         verbose output? (includes estimation of the current memory footprint) [default: 0]
  % ? warnings          do some sanity checks and throw warnings? [default: 0]
  % ? filebase          basename of the output files [default: uvwp]
  % ? iseq              sequence number of the output files [default: 0]
  % ? time              simulation time of the output files [default: 0.0]
  
  % Parse input
  par = inputParser;
  addParamValue(par,'verbosity',0,@isnumeric);
  addParamValue(par,'warnings',0,@isnumeric);
  addParamValue(par,'filebase','uvwp',@ischar);
  addParamValue(par,'iseq',0,@isnumeric);
  addParamValue(par,'time',0.0,@isnumeric);
  parse(par,varargin{:});
  flag_verb = par.Results.verbosity;
  flag_warn = par.Results.warning;
  filebase  = par.Results.filebase;
  iseq      = par.Results.iseq;
  simtime   = par.Results.time;

  % Read info from tar-archive
  [xuo,yuo,zuo,xvo,yvo,zvo,xwo,ywo,zwo,xpo,ypo,zpo] = read_grid_ucf(hucf);
  [params] = read_parameters_ucf(hucf);
  nprocso = params.parallel.nprocs;
  a = params.domain.a;
  b = params.domain.b;
  c = params.domain.c;
  d = params.domain.d;
  e = params.domain.e;
  f = params.domain.f;
  xperiodic = params.domain.xperiodic;
  yperiodic = params.domain.yperiodic;
  zperiodic = params.domain.zperiodic;

  [nxun,nyun,nzun,nxvn,nyvn,nzvn,nxwn,nywn,nzwn,...
   xun,yun,zun,xvn,yvn,zvn,xwn,ywn,zwn,xpn,ypn,zpn,dxn,dyn,dzn]=...
              generate_grid(a,b,c,d,e,f,nxpn,nypn,nzpn,xperiodic,yperiodic,zperiodic);
  [ibegun,iendun,jbegun,jendun,kbegun,kendun,...
   ibegvn,iendvn,jbegvn,jendvn,kbegvn,kendvn,...
   ibegwn,iendwn,jbegwn,jendwn,kbegwn,kendwn,...
   ibegpn,iendpn,jbegpn,jendpn,kbegpn,kendpn]=generate_procgrid(...
                nxun,nyun,nzun,nxvn,nyvn,nzvn,nxwn,nywn,nzwn,...
                nxpn,nypn,nzpn,nxprocsn,nyprocsn,nzprocsn);
                
  if flag_warn 
    if (~xperiodic && mod(nxpn,2)==0) || (~yperiodic && mod(nypn,2)==0) || (~zperiodic && mod(nzpn,2)==0)
      warning('Pressure grid should contain extra point in non-periodic directions.')
    end
    if abs((dxn-dyn)/dxn)>1e-8 || abs((dyn-dzn)/dyn)>1e-8
      warning('New mesh is not equidistant: %f, %f, %f.',dxn,dyn,dzn);
    end
    if mod(nxpn-~xperiodic,64)~=0
      warning('nxpn is not a multiple of 64: MGD performance might be bad.');
    end
    if mod(nypn-~yperiodic,64)~=0
      warning('nypn is not a multiple of 64: MGD performance might be bad.');
    end
    if mod(nzpn-~xperiodic,64)~=0
      warning('nzpn is not a multiple of 64: MGD performance might be bad.');
    end
  end

  % Cell arrays which hold interpolation results
  un = cell(nxprocsn,nyprocsn,nzprocsn);
  vn = cell(nxprocsn,nyprocsn,nzprocsn);
  wn = cell(nxprocsn,nyprocsn,nzprocsn);
  pn = cell(nxprocsn,nyprocsn,nzprocsn);

  % Cell arrays which indicate if points have been interpolated yet
  un_ind = cell(nxprocsn,nyprocsn,nzprocsn);
  vn_ind = cell(nxprocsn,nyprocsn,nzprocsn);
  wn_ind = cell(nxprocsn,nyprocsn,nzprocsn);
  pn_ind = cell(nxprocsn,nyprocsn,nzprocsn);

  % Logical arrays which indicate whether memory has been allocated
  un_init = false(nxprocsn,nyprocsn,nzprocsn);
  vn_init = false(nxprocsn,nyprocsn,nzprocsn);
  wn_init = false(nxprocsn,nyprocsn,nzprocsn);
  pn_init = false(nxprocsn,nyprocsn,nzprocsn);

  % Logical array which indicate whether a chunk has been finalized
  chunk_final = false(nxprocsn,nyprocsn,nzprocsn);

  imem = 0;
  maxmem = 0;
  for iproco=0:nprocso-1
    if flag_verb
      fprintf('Processing original chunk %5d ',iproco);
    end
    [uo,ibuo,jbuo,kbuo,nxuol,nyuol,nzuol,...
     vo,ibvo,jbvo,kbvo,nxvol,nyvol,nzvol,...
     wo,ibwo,jbwo,kbwo,nxwol,nywol,nzwol,...
     po,ibpo,jbpo,kbpo,nxpol,nypol,nzpol,ighost] = ...
              read_uvwp_chunk_ucf(hucf,'rank',iproco);
    imem = imem+(numel(uo)+numel(vo)+numel(wo)+numel(po))*8;
    if flag_verb
      fprintf('(current memory footprint: %7.0f MiB)\n',imem/(1024*1024));
    end
    maxmem = max(imem,maxmem);
    
    [xuol,yuol,zuol] = grid_chunk(xuo,yuo,zuo,ibuo,jbuo,kbuo,nxuol,nyuol,nzuol,ighost);
    [xvol,yvol,zvol] = grid_chunk(xvo,yvo,zvo,ibvo,jbvo,kbvo,nxvol,nyvol,nzvol,ighost);
    [xwol,ywol,zwol] = grid_chunk(xwo,ywo,zwo,ibwo,jbwo,kbwo,nxwol,nywol,nzwol,ighost);
    [xpol,ypol,zpol] = grid_chunk(xpo,ypo,zpo,ibpo,jbpo,kbpo,nxpol,nypol,nzpol,ighost);
    
    Fu = griddedInterpolant({xuol,yuol,zuol},uo);
    Fv = griddedInterpolant({xvol,yvol,zvol},vo);
    Fw = griddedInterpolant({xwol,ywol,zwol},wo);
    Fp = griddedInterpolant({xpol,ypol,zpol},po);
    
    for izprocn=0:nzprocsn-1
      for iyprocn=0:nyprocsn-1
        for ixprocn=0:nxprocsn-1
          if chunk_final(ixprocn+1,iyprocn+1,izprocn+1)
            continue;
          end
          ibun = ibegun(ixprocn+1); nxunl = iendun(ixprocn+1)-ibegun(ixprocn+1)+1;
          jbun = jbegun(iyprocn+1); nyunl = jendun(iyprocn+1)-jbegun(iyprocn+1)+1;
          kbun = kbegun(izprocn+1); nzunl = kendun(izprocn+1)-kbegun(izprocn+1)+1;
          ibvn = ibegvn(ixprocn+1); nxvnl = iendvn(ixprocn+1)-ibegvn(ixprocn+1)+1;
          jbvn = jbegvn(iyprocn+1); nyvnl = jendvn(iyprocn+1)-jbegvn(iyprocn+1)+1;
          kbvn = kbegvn(izprocn+1); nzvnl = kendvn(izprocn+1)-kbegvn(izprocn+1)+1;
          ibwn = ibegwn(ixprocn+1); nxwnl = iendwn(ixprocn+1)-ibegwn(ixprocn+1)+1;
          jbwn = jbegwn(iyprocn+1); nywnl = jendwn(iyprocn+1)-jbegwn(iyprocn+1)+1;
          kbwn = kbegwn(izprocn+1); nzwnl = kendwn(izprocn+1)-kbegwn(izprocn+1)+1;
          ibpn = ibegpn(ixprocn+1); nxpnl = iendpn(ixprocn+1)-ibegpn(ixprocn+1)+1;
          jbpn = jbegpn(iyprocn+1); nypnl = jendpn(iyprocn+1)-jbegpn(iyprocn+1)+1;
          kbpn = kbegpn(izprocn+1); nzpnl = kendpn(izprocn+1)-kbegpn(izprocn+1)+1;
          [xunl,yunl,zunl] = grid_chunk(xun,yun,zun,ibun,jbun,kbun,nxunl,nyunl,nzunl,0);
          [xvnl,yvnl,zvnl] = grid_chunk(xvn,yvn,zvn,ibvn,jbvn,kbvn,nxvnl,nyvnl,nzvnl,0);
          [xwnl,ywnl,zwnl] = grid_chunk(xwn,ywn,zwn,ibwn,jbwn,kbwn,nxwnl,nywnl,nzwnl,0);
          [xpnl,ypnl,zpnl] = grid_chunk(xpn,ypn,zpn,ibpn,jbpn,kbpn,nxpnl,nypnl,nzpnl,0);
          
          % u
          ixu = find(xunl>=xuol(1) & xunl<=xuol(end));
          iyu = find(yunl>=yuol(1) & yunl<=yuol(end));
          izu = find(zunl>=zuol(1) & zunl<=zuol(end));
          if ~isempty(ixu) && ~isempty(iyu) && ~isempty(izu)
            if ~un_init(ixprocn+1,iyprocn+1,izprocn+1)
              un{ixprocn+1,iyprocn+1,izprocn+1} = zeros(nxunl,nyunl,nzunl);
              un_ind{ixprocn+1,iyprocn+1,izprocn+1} = false(nxunl,nyunl,nzunl);
              un_init(ixprocn+1,iyprocn+1,izprocn+1) = true;
              imem = imem+(nxunl*nyunl*nzunl)*(8+1);
              maxmem = max(imem,maxmem);
              if ~xperiodic 
                % staggered grid goes beyond boundary in non-periodic direction:
                % if field is coarsened, no interpolation point will be found!
                % Thus, just leave it and mark as processed
                un_ind{ixprocn+1,iyprocn+1,izprocn+1}(xunl<xuo(1),:,:) = true;
                un_ind{ixprocn+1,iyprocn+1,izprocn+1}(xunl>xuo(end),:,:) = true;
              end
            end
            un{ixprocn+1,iyprocn+1,izprocn+1}(ixu,iyu,izu) = Fu({xunl(ixu),yunl(iyu),zunl(izu)});
            un_ind{ixprocn+1,iyprocn+1,izprocn+1}(ixu,iyu,izu) = true;
          end
          
          % v
          ixv = find(xvnl>=xvol(1) & xvnl<=xvol(end));
          iyv = find(yvnl>=yvol(1) & yvnl<=yvol(end));
          izv = find(zvnl>=zvol(1) & zvnl<=zvol(end));
          if ~isempty(ixv) && ~isempty(iyv) && ~isempty(izv)
            if ~vn_init(ixprocn+1,iyprocn+1,izprocn+1)
              vn{ixprocn+1,iyprocn+1,izprocn+1} = zeros(nxvnl,nyvnl,nzvnl);
              vn_ind{ixprocn+1,iyprocn+1,izprocn+1} = false(nxvnl,nyvnl,nzvnl);
              vn_init(ixprocn+1,iyprocn+1,izprocn+1) = true;
              imem = imem+(nxunl*nyunl*nzunl)*(8+1);
              maxmem = max(imem,maxmem);
              if ~yperiodic 
                % staggered grid goes beyond boundary in non-periodic direction:
                % if field is coarsened, no interpolation point will be found!
                % Thus, just leave it and mark as processed
                vn_ind{ixprocn+1,iyprocn+1,izprocn+1}(:,yvnl<yvo(1),:) = true;
                vn_ind{ixprocn+1,iyprocn+1,izprocn+1}(:,yvnl>yvo(end),:) = true;
              end
            end
            vn{ixprocn+1,iyprocn+1,izprocn+1}(ixv,iyv,izv) = Fv({xvnl(ixv),yvnl(iyv),zvnl(izv)});
            vn_ind{ixprocn+1,iyprocn+1,izprocn+1}(ixv,iyv,izv) = true;
          end
          
          % w
          ixw = find(xwnl>=xwol(1) & xwnl<=xwol(end));
          iyw = find(ywnl>=ywol(1) & ywnl<=ywol(end));
          izw = find(zwnl>=zwol(1) & zwnl<=zwol(end));
          if ~isempty(ixw) && ~isempty(iyw) && ~isempty(izw)
            if ~wn_init(ixprocn+1,iyprocn+1,izprocn+1)
              wn{ixprocn+1,iyprocn+1,izprocn+1} = zeros(nxwnl,nywnl,nzwnl);
              wn_ind{ixprocn+1,iyprocn+1,izprocn+1} = false(nxwnl,nywnl,nzwnl);
              wn_init(ixprocn+1,iyprocn+1,izprocn+1) = true;
              imem = imem+(nxunl*nyunl*nzunl)*(8+1);
              maxmem = max(imem,maxmem);
              if ~zperiodic 
                % staggered grid goes beyond boundary in non-periodic direction:
                % if field is coarsened, no interpolation point will be found!
                % Thus, just leave it and mark as processed
                wn_ind{ixprocn+1,iyprocn+1,izprocn+1}(:,:,zwnl<zwo(1)) = true;
                wn_ind{ixprocn+1,iyprocn+1,izprocn+1}(:,:,zwnl>zwo(end)) = true;
              end
            end
            wn{ixprocn+1,iyprocn+1,izprocn+1}(ixw,iyw,izw) = Fw({xwnl(ixw),ywnl(iyw),zwnl(izw)});
            wn_ind{ixprocn+1,iyprocn+1,izprocn+1}(ixw,iyw,izw) = true;
          end
          
          % p
          ixp = find(xpnl>=xpol(1) & xpnl<=xpol(end));
          iyp = find(ypnl>=ypol(1) & ypnl<=ypol(end));
          izp = find(zpnl>=zpol(1) & zpnl<=zpol(end));
          if ~isempty(ixp) && ~isempty(iyp) && ~isempty(izp)
            if ~pn_init(ixprocn+1,iyprocn+1,izprocn+1)
              pn{ixprocn+1,iyprocn+1,izprocn+1} = zeros(nxpnl,nypnl,nzpnl);
              pn_ind{ixprocn+1,iyprocn+1,izprocn+1} = false(nxpnl,nypnl,nzpnl);
              pn_init(ixprocn+1,iyprocn+1,izprocn+1) = true;
              imem = imem+(nxunl*nyunl*nzunl)*(8+1);
              maxmem = max(imem,maxmem);
            end
            pn{ixprocn+1,iyprocn+1,izprocn+1}(ixp,iyp,izp) = Fp({xpnl(ixp),ypnl(iyp),zpnl(izp)});
            pn_ind{ixprocn+1,iyprocn+1,izprocn+1}(ixp,iyp,izp) = true;
          end
          
          if un_init(ixprocn+1,iyprocn+1,izprocn+1) && ...
             vn_init(ixprocn+1,iyprocn+1,izprocn+1) && ...
             wn_init(ixprocn+1,iyprocn+1,izprocn+1) && ...
             pn_init(ixprocn+1,iyprocn+1,izprocn+1) && ...
             all(un_ind{ixprocn+1,iyprocn+1,izprocn+1}(:)) && ...
             all(vn_ind{ixprocn+1,iyprocn+1,izprocn+1}(:)) && ...
             all(wn_ind{ixprocn+1,iyprocn+1,izprocn+1}(:)) && ...
             all(pn_ind{ixprocn+1,iyprocn+1,izprocn+1}(:))
            ichunkn = ixprocn*nyprocsn*nzprocsn+iyprocn*nzprocsn+izprocn;
            fout = sprintf('%s/%s_%04d.%05d',dout,filebase,iseq,ichunkn);
            write_uvwp_chunk_ucf(fout,...
               un{ixprocn+1,iyprocn+1,izprocn+1},ibun,jbun,kbun,nxunl,nyunl,nzunl,...
               vn{ixprocn+1,iyprocn+1,izprocn+1},ibvn,jbvn,kbvn,nxvnl,nyvnl,nzvnl,...
               wn{ixprocn+1,iyprocn+1,izprocn+1},ibwn,jbwn,kbwn,nxwnl,nywnl,nzwnl,...
               pn{ixprocn+1,iyprocn+1,izprocn+1},ibpn,jbpn,kbpn,nxpnl,nypnl,nzpnl,...
                   simtime,0);
            un{ixprocn+1,iyprocn+1,izprocn+1} = [];
            vn{ixprocn+1,iyprocn+1,izprocn+1} = [];
            wn{ixprocn+1,iyprocn+1,izprocn+1} = [];
            pn{ixprocn+1,iyprocn+1,izprocn+1} = [];
            un_ind{ixprocn+1,iyprocn+1,izprocn+1} = [];
            vn_ind{ixprocn+1,iyprocn+1,izprocn+1} = [];
            wn_ind{ixprocn+1,iyprocn+1,izprocn+1} = [];
            pn_ind{ixprocn+1,iyprocn+1,izprocn+1} = [];
            chunk_final(ixprocn+1,iyprocn+1,izprocn+1) = true;
            imem = imem-(...
              (nxunl*nyunl*nzunl)+...
              (nxvnl*nyvnl*nzvnl)+...
              (nxwnl*nywnl*nzwnl)+...
              (nxpnl*nypnl*nzpnl) ...
              )*(8+1);
          end
        end
      end
    end 
    imem = imem-(numel(uo)+numel(vo)+numel(wo)+numel(po))*8;
  end

  % Check if everything has been finalized
  for izprocn=0:nzprocsn-1
    for iyprocn=0:nyprocsn-1
      for ixprocn=0:nxprocsn-1
        if ~chunk_final(ixprocn+1,iyprocn+1,izprocn+1)
          ichunkn = ixprocn*nyprocsn*nzprocsn+iyprocn*nzprocsn+izprocn;
          fout = sprintf('%s/%s_%04d.%05d',dout,filebase,iseq,ichunkn);
          if flag_warn
            warning('Chunk %2d,%2d,%2d (rank: %5d) has not been finalized! Writing anyway...',ixprocn,iyprocn,izprocn,ichunkn);
          end
          write_uvwp_chunk_ucf(fout,...
             un{ixprocn+1,iyprocn+1,izprocn+1},ibun,jbun,kbun,nxunl,nyunl,nzunl,...
             vn{ixprocn+1,iyprocn+1,izprocn+1},ibvn,jbvn,kbvn,nxvnl,nyvnl,nzvnl,...
             wn{ixprocn+1,iyprocn+1,izprocn+1},ibwn,jbwn,kbwn,nxwnl,nywnl,nzwnl,...
             pn{ixprocn+1,iyprocn+1,izprocn+1},ibpn,jbpn,kbpn,nxpnl,nypnl,nzpnl,...
                 simetime,0);
        end
      end
    end
  end
  if flag_verb
    fprintf('Done. (maximum memory footprint: %7.0f MiB)\n',maxmem/(1024*1024));
  end
end