function [] = precompute_ghosts_uvwp_ucf(hucf,dout,nghost,varargin)
  % [] = precompute_ghosts_uvwp_ucf(hucf,dout,nghost,varargin)
  % Constructs a given number of ghost cells per processor and saves them to
  % mat-files.
  % The ghost cell generation is performed with a single read operation (input sweep),
  % which extracts the points needed by the neighbors, and a subsequent processing
  % sweep, which exchanges the ghost cell data. The data is hold in memory using
  % a sparse data structure (c.f. https://de.mathworks.com/matlabcentral/fileexchange/29832-n-dimensional-sparse-arrays)
  % and is also written to file in that way, so the 'ndSparse' class is a prerequisite.
  % This way, the ghost cells can be added to a chunk of data with a simple add
  % operation while avoiding large consumption of disk space.
  % Input
  %   hucf              UCF handle (ustar,ucfmulti)
  %   dout              output directory
  %   nghost            number of ghost cells to compute
  %                       periodic:       take data from neighboring processor
  %                       non-periodic:   duplicate last datapoint
  % ? verbosity         verbose output? (includes estimation of the memory footprint) [default: 0]
  % ? filebase          basename of the output files [default: ghost<nghost>]
  % ? iseq              sequence number of the output files [default: 0]

  % Parse input
  par = inputParser;
  addParamValue(par,'verbosity',0,@isnumeric);
  addParamValue(par,'filebase',sprintf('ghost%d',nghost),@ischar);
  addParamValue(par,'iseq',0,@isnumeric);
  parse(par,varargin{:});
  flag_verb = par.Results.verbosity;
  filebase  = par.Results.filebase;
  iseq      = par.Results.iseq;
  
  % Read info from tar-archive
  [ibegu,iendu,jbegu,jendu,kbegu,kendu,...
   ibegv,iendv,jbegv,jendv,kbegv,kendv,...
   ibegw,iendw,jbegw,jendw,kbegw,kendw,...
   ibegp,iendp,jbegp,jendp,kbegp,kendp] = read_procgrid_ucf(hucf);
  [params] = read_parameters_ucf(hucf);
  nxp = params.mesh.nxp;
  nyp = params.mesh.nyp;
  nzp = params.mesh.nzp;
  nxprocs = params.parallel.nxprocs;
  nyprocs = params.parallel.nyprocs;
  nzprocs = params.parallel.nzprocs;
  xperiodic = params.geometry.xperiodic;
  yperiodic = params.geometry.yperiodic;
  zperiodic = params.geometry.zperiodic;
  
  % Buffer for boundary data of each chunk
  ubd = cell(nxprocs,nyprocs,nzprocs);
  vbd = cell(nxprocs,nyprocs,nzprocs);
  wbd = cell(nxprocs,nyprocs,nzprocs);
  pbd = cell(nxprocs,nyprocs,nzprocs);
  
  % Estimate memory requirement (conservative simplification)
  if flag_verb
    nxl_estm = floor(nxp/nxprocs);
    nyl_estm = floor(nyp/nyprocs);
    nzl_estm = floor(nzp/nzprocs);
    imem = 0;
    for ii=0:nghost-1
      imem = imem + 4*8*(...
              2*(nxl_estm-ii)*(nyl_estm-ii) + ...
              2*(nyl_estm-ii)*(nzl_estm-ii) + ...
              2*(nxl_estm-ii)*(nzl_estm-ii)) * ...
              (nxprocs*nyprocs*nzprocs);
    end
    fprintf('Estimated memory requirement: %7.0f MiB\n',imem/(1024*1024));
  end
  
  % Read boundary data of each chunk into memory
  for ixproc=0:nxprocs-1
    for iyproc=0:nyprocs-1
      for izproc=0:nzprocs-1
        ichunk = ixproc*nyprocs*nzprocs+iyproc*nzprocs+izproc;
        if flag_verb
          fprintf('Reading original chunk %5d\n',ichunk);
        end
        [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_ucf(hucf,'rank',ichunk,'ghost',0);
        % Set interior to zero
        u(1+nghost:end-nghost,1+nghost:end-nghost,1+nghost:end-nghost) = zeros(nxul-2*nghost,nyul-2*nghost,nzul-2*nghost);
        v(1+nghost:end-nghost,1+nghost:end-nghost,1+nghost:end-nghost) = zeros(nxvl-2*nghost,nyvl-2*nghost,nzvl-2*nghost);
        w(1+nghost:end-nghost,1+nghost:end-nghost,1+nghost:end-nghost) = zeros(nxwl-2*nghost,nywl-2*nghost,nzwl-2*nghost);
        p(1+nghost:end-nghost,1+nghost:end-nghost,1+nghost:end-nghost) = zeros(nxpl-2*nghost,nypl-2*nghost,nzpl-2*nghost);
        % Save ghost points as sparse matrix
        ubd{ixproc+1,iyproc+1,izproc+1} = ndSparse(u);
        vbd{ixproc+1,iyproc+1,izproc+1} = ndSparse(v);
        wbd{ixproc+1,iyproc+1,izproc+1} = ndSparse(w);
        pbd{ixproc+1,iyproc+1,izproc+1} = ndSparse(p);
      end
    end
  end
  
  % Construct ghost cells for each chunk
  for ixproc=0:nxprocs-1
    nxul = iendu(ixproc+1)-ibegu(ixproc+1)+1;
    nxvl = iendv(ixproc+1)-ibegv(ixproc+1)+1;
    nxwl = iendw(ixproc+1)-ibegw(ixproc+1)+1;
    nxpl = iendp(ixproc+1)-ibegp(ixproc+1)+1;
    for iyproc=0:nyprocs-1
      nyul = jendu(iyproc+1)-jbegu(iyproc+1)+1;
      nyvl = jendv(iyproc+1)-jbegv(iyproc+1)+1;
      nywl = jendw(iyproc+1)-jbegw(iyproc+1)+1;
      nypl = jendp(iyproc+1)-jbegp(iyproc+1)+1;
      for izproc=0:nzprocs-1
        nzul = kendu(izproc+1)-kbegu(izproc+1)+1;
        nzvl = kendv(izproc+1)-kbegv(izproc+1)+1;
        nzwl = kendw(izproc+1)-kbegw(izproc+1)+1;
        nzpl = kendp(izproc+1)-kbegp(izproc+1)+1;
        
        ichunk = ixproc*nyprocs*nzprocs+iyproc*nzprocs+izproc;
        if flag_verb
          fprintf('Constructing ghost cells for chunk %5d\n',ichunk);
        end
        
        % Create sparse arrays which hold the results
        ugh = ones(nxul+2*nghost,nyul+2*nghost,nzul+2*nghost);
        ugh(1+nghost:end-nghost,1+nghost:end-nghost,1+nghost:end-nghost) = zeros(nxul,nyul,nzul);
        ugh = ndSparse(ugh);
        vgh = ones(nxvl+2*nghost,nyvl+2*nghost,nzvl+2*nghost);
        vgh(1+nghost:end-nghost,1+nghost:end-nghost,1+nghost:end-nghost) = zeros(nxvl,nyvl,nzvl);
        vgh = ndSparse(vgh);
        wgh = ones(nxwl+2*nghost,nywl+2*nghost,nzwl+2*nghost);
        wgh(1+nghost:end-nghost,1+nghost:end-nghost,1+nghost:end-nghost) = zeros(nxwl,nywl,nzwl);
        wgh = ndSparse(wgh);
        pgh = ones(nxpl+2*nghost,nypl+2*nghost,nzpl+2*nghost);
        pgh(1+nghost:end-nghost,1+nghost:end-nghost,1+nghost:end-nghost) = zeros(nxpl,nypl,nzpl);
        pgh = ndSparse(pgh);
        
        % Fill them with ghost cell data
        % u
        % 1st faces (6)
        ugh(1:nghost,           1+nghost:end-nghost,1+nghost:end-nghost) = ubd{mod(ixproc-1,nxprocs)+1,iyproc+1,izproc+1}(end-nghost+1:end,:,               :               );
        ugh(end-nghost+1:end,   1+nghost:end-nghost,1+nghost:end-nghost) = ubd{mod(ixproc+1,nxprocs)+1,iyproc+1,izproc+1}(1:nghost,        :,               :               );
        ugh(1+nghost:end-nghost,1:nghost,           1+nghost:end-nghost) = ubd{ixproc+1,mod(iyproc-1,nyprocs)+1,izproc+1}(:,               end-nghost+1:end,:               );
        ugh(1+nghost:end-nghost,end-nghost+1:end,   1+nghost:end-nghost) = ubd{ixproc+1,mod(iyproc+1,nyprocs)+1,izproc+1}(:,               1:nghost,        :               );
        ugh(1+nghost:end-nghost,1+nghost:end-nghost,1:nghost           ) = ubd{ixproc+1,iyproc+1,mod(izproc-1,nzprocs)+1}(:,               :,               end-nghost+1:end);
        ugh(1+nghost:end-nghost,1+nghost:end-nghost,end-nghost+1:end   ) = ubd{ixproc+1,iyproc+1,mod(izproc+1,nzprocs)+1}(:,               :,               1:nghost        );
        % 2nd edges (12)
        ugh(1:nghost,           1:nghost,           1+nghost:end-nghost) = ubd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,izproc+1               }(end-nghost+1:end,end-nghost+1:end,:);
        ugh(1:nghost,           1+nghost:end-nghost,1:nghost           ) = ubd{mod(ixproc-1,nxprocs)+1,iyproc+1,               mod(izproc-1,nzprocs)+1}(end-nghost+1:end,:,               end-nghost+1:end);
        ugh(1+nghost:end-nghost,1:nghost,           1:nghost           ) = ubd{ixproc+1,               mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(:,               end-nghost+1:end,end-nghost+1:end);
        ugh(end-nghost+1:end,   end-nghost+1:end,   1+nghost:end-nghost) = ubd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,izproc+1               }(1:nghost,        1:nghost,        :               );
        ugh(end-nghost+1:end,   1+nghost:end-nghost,end-nghost+1:end   ) = ubd{mod(ixproc+1,nxprocs)+1,iyproc+1,               mod(izproc+1,nzprocs)+1}(1:nghost,        :,               1:nghost        );
        ugh(1+nghost:end-nghost,end-nghost+1:end,   end-nghost+1:end   ) = ubd{ixproc+1,               mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(:,               1:nghost,        1:nghost        );
        ugh(1:nghost,           end-nghost+1:end,   1+nghost:end-nghost) = ubd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,izproc+1               }(end-nghost+1:end,1:nghost,        :               );
        ugh(1:nghost,           1+nghost:end-nghost,end-nghost+1:end   ) = ubd{mod(ixproc-1,nxprocs)+1,iyproc+1,               mod(izproc+1,nzprocs)+1}(end-nghost+1:end,:,               1:nghost        );
        ugh(1+nghost:end-nghost,1:nghost,           end-nghost+1:end   ) = ubd{ixproc+1,               mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(:,               end-nghost+1:end,1:nghost        );
        ugh(end-nghost+1:end,   1:nghost,           1+nghost:end-nghost) = ubd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,izproc+1               }(1:nghost,        end-nghost+1:end,:               );
        ugh(end-nghost+1:end,   1+nghost:end-nghost,1:nghost           ) = ubd{mod(ixproc+1,nxprocs)+1,iyproc+1,               mod(izproc-1,nzprocs)+1}(1:nghost,        :,               end-nghost+1:end);
        ugh(1+nghost:end-nghost,end-nghost+1:end,   1:nghost           ) = ubd{ixproc+1,               mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(:,               1:nghost,        end-nghost+1:end);
        % 3rd corners (8)
        ugh(1:nghost,        1:nghost,        1:nghost        ) = ubd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(end-nghost+1:end,end-nghost+1:end,end-nghost+1:end);
        ugh(end-nghost+1:end,1:nghost,        1:nghost        ) = ubd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(1:nghost,        end-nghost+1:end,end-nghost+1:end);
        ugh(1:nghost,        end-nghost+1:end,1:nghost        ) = ubd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(end-nghost+1:end,1:nghost,        end-nghost+1:end);
        ugh(1:nghost,        1:nghost,        end-nghost+1:end) = ubd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(end-nghost+1:end,end-nghost+1:end,1:nghost        );
        ugh(end-nghost+1:end,end-nghost+1:end,1:nghost        ) = ubd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(1:nghost,        1:nghost,        end-nghost+1:end);
        ugh(end-nghost+1:end,1:nghost,        end-nghost+1:end) = ubd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(1:nghost,        end-nghost+1:end,1:nghost        );
        ugh(1:nghost,        end-nghost+1:end,end-nghost+1:end) = ubd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(end-nghost+1:end,1:nghost,        1:nghost        );
        ugh(end-nghost+1:end,end-nghost+1:end,end-nghost+1:end) = ubd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(1:nghost,        1:nghost,        1:nghost        );
        % Correct if not periodic
        if ~xperiodic && ixproc==0
          ugh(1:nghost,        :,                  :                  ) = repmat(ugh(nghost+1,:,:),[nghost,1,1]);
          ugh(1:nghost,        1+nghost:end-nghost,1+nghost:end-nghost) = repmat(ubd{ixproc+1,iyproc+1,izproc+1}(1,:,:),[nghost,1,1]);
        end
        if ~xperiodic && ixproc==nxprocs-1
          ugh(end-nghost+1:end,:,                  :                  ) = repmat(ugh(end-nghost,:,:),[nghost,1,1]);
          ugh(end-nghost+1:end,1+nghost:end-nghost,1+nghost:end-nghost) = repmat(ubd{ixproc+1,iyproc+1,izproc+1}(end,:,:),[nghost,1,1]);
        end
        if ~yperiodic && iyproc==0
          ugh(:,                  1:nghost,        :                  ) = repmat(ugh(:,nghost+1,:),[1,nghost,1]);
          ugh(1+nghost:end-nghost,1:nghost,        1+nghost:end-nghost) = repmat(ubd{ixproc+1,iyproc+1,izproc+1}(:,1,:),[1,nghost,1]);
        end
        if ~yperiodic && iyproc==nyprocs-1
          ugh(:,                  end-nghost+1:end,:                  ) = repmat(ugh(:,end-nghost,:),[1,nghost,1]);
          ugh(1+nghost:end-nghost,end-nghost+1:end,1+nghost:end-nghost) = repmat(ubd{ixproc+1,iyproc+1,izproc+1}(:,end,:),[1,nghost,1]);
        end
        if ~zperiodic && izproc==0
          ugh(:,                  :,                  1:nghost        ) = repmat(ugh(:,:,nghost+1),[1,1,nghost]);
          ugh(1+nghost:end-nghost,1+nghost:end-nghost,1:nghost        ) = repmat(ubd{ixproc+1,iyproc+1,izproc+1}(:,:,1),[1,1,nghost]);
        end
        if ~zperiodic && izproc==nzprocs-1
          ugh(:,                  :,                  end-nghost+1:end) = repmat(ugh(:,:,end-nghost),[1,1,nghost]);
          ugh(1+nghost:end-nghost,1+nghost:end-nghost,end-nghost+1:end) = repmat(ubd{ixproc+1,iyproc+1,izproc+1}(:,:,end),[1,1,nghost]);
        end
        
        % Fill them with ghost cell data
        % v
        % 1st faces (6)
        vgh(1:nghost,           1+nghost:end-nghost,1+nghost:end-nghost) = vbd{mod(ixproc-1,nxprocs)+1,iyproc+1,izproc+1}(end-nghost+1:end,:,               :               );
        vgh(end-nghost+1:end,   1+nghost:end-nghost,1+nghost:end-nghost) = vbd{mod(ixproc+1,nxprocs)+1,iyproc+1,izproc+1}(1:nghost,        :,               :               );
        vgh(1+nghost:end-nghost,1:nghost,           1+nghost:end-nghost) = vbd{ixproc+1,mod(iyproc-1,nyprocs)+1,izproc+1}(:,               end-nghost+1:end,:               );
        vgh(1+nghost:end-nghost,end-nghost+1:end,   1+nghost:end-nghost) = vbd{ixproc+1,mod(iyproc+1,nyprocs)+1,izproc+1}(:,               1:nghost,        :               );
        vgh(1+nghost:end-nghost,1+nghost:end-nghost,1:nghost           ) = vbd{ixproc+1,iyproc+1,mod(izproc-1,nzprocs)+1}(:,               :,               end-nghost+1:end);
        vgh(1+nghost:end-nghost,1+nghost:end-nghost,end-nghost+1:end   ) = vbd{ixproc+1,iyproc+1,mod(izproc+1,nzprocs)+1}(:,               :,               1:nghost        );
        % 2nd edges (12)
        vgh(1:nghost,           1:nghost,           1+nghost:end-nghost) = vbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,izproc+1               }(end-nghost+1:end,end-nghost+1:end,:);
        vgh(1:nghost,           1+nghost:end-nghost,1:nghost           ) = vbd{mod(ixproc-1,nxprocs)+1,iyproc+1,               mod(izproc-1,nzprocs)+1}(end-nghost+1:end,:,               end-nghost+1:end);
        vgh(1+nghost:end-nghost,1:nghost,           1:nghost           ) = vbd{ixproc+1,               mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(:,               end-nghost+1:end,end-nghost+1:end);
        vgh(end-nghost+1:end,   end-nghost+1:end,   1+nghost:end-nghost) = vbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,izproc+1               }(1:nghost,        1:nghost,        :               );
        vgh(end-nghost+1:end,   1+nghost:end-nghost,end-nghost+1:end   ) = vbd{mod(ixproc+1,nxprocs)+1,iyproc+1,               mod(izproc+1,nzprocs)+1}(1:nghost,        :,               1:nghost        );
        vgh(1+nghost:end-nghost,end-nghost+1:end,   end-nghost+1:end   ) = vbd{ixproc+1,               mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(:,               1:nghost,        1:nghost        );
        vgh(1:nghost,           end-nghost+1:end,   1+nghost:end-nghost) = vbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,izproc+1               }(end-nghost+1:end,1:nghost,        :               );
        vgh(1:nghost,           1+nghost:end-nghost,end-nghost+1:end   ) = vbd{mod(ixproc-1,nxprocs)+1,iyproc+1,               mod(izproc+1,nzprocs)+1}(end-nghost+1:end,:,               1:nghost        );
        vgh(1+nghost:end-nghost,1:nghost,           end-nghost+1:end   ) = vbd{ixproc+1,               mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(:,               end-nghost+1:end,1:nghost        );
        vgh(end-nghost+1:end,   1:nghost,           1+nghost:end-nghost) = vbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,izproc+1               }(1:nghost,        end-nghost+1:end,:               );
        vgh(end-nghost+1:end,   1+nghost:end-nghost,1:nghost           ) = vbd{mod(ixproc+1,nxprocs)+1,iyproc+1,               mod(izproc-1,nzprocs)+1}(1:nghost,        :,               end-nghost+1:end);
        vgh(1+nghost:end-nghost,end-nghost+1:end,   1:nghost           ) = vbd{ixproc+1,               mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(:,               1:nghost,        end-nghost+1:end);
        % 3rd corners (8)
        vgh(1:nghost,        1:nghost,        1:nghost        ) = vbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(end-nghost+1:end,end-nghost+1:end,end-nghost+1:end);
        vgh(end-nghost+1:end,1:nghost,        1:nghost        ) = vbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(1:nghost,        end-nghost+1:end,end-nghost+1:end);
        vgh(1:nghost,        end-nghost+1:end,1:nghost        ) = vbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(end-nghost+1:end,1:nghost,        end-nghost+1:end);
        vgh(1:nghost,        1:nghost,        end-nghost+1:end) = vbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(end-nghost+1:end,end-nghost+1:end,1:nghost        );
        vgh(end-nghost+1:end,end-nghost+1:end,1:nghost        ) = vbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(1:nghost,        1:nghost,        end-nghost+1:end);
        vgh(end-nghost+1:end,1:nghost,        end-nghost+1:end) = vbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(1:nghost,        end-nghost+1:end,1:nghost        );
        vgh(1:nghost,        end-nghost+1:end,end-nghost+1:end) = vbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(end-nghost+1:end,1:nghost,        1:nghost        );
        vgh(end-nghost+1:end,end-nghost+1:end,end-nghost+1:end) = vbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(1:nghost,        1:nghost,        1:nghost        );
        % Correct if not periodic
        if ~xperiodic && ixproc==0
          vgh(1:nghost,        :,                  :                  ) = repmat(vgh(nghost+1,:,:),[nghost,1,1]);
          vgh(1:nghost,        1+nghost:end-nghost,1+nghost:end-nghost) = repmat(vbd{ixproc+1,iyproc+1,izproc+1}(1,:,:),[nghost,1,1]);
        end
        if ~xperiodic && ixproc==nxprocs-1
          vgh(end-nghost+1:end,:,                  :                  ) = repmat(vgh(end-nghost,:,:),[nghost,1,1]);
          vgh(end-nghost+1:end,1+nghost:end-nghost,1+nghost:end-nghost) = repmat(vbd{ixproc+1,iyproc+1,izproc+1}(end,:,:),[nghost,1,1]);
        end
        if ~yperiodic && iyproc==0
          vgh(:,                  1:nghost,        :                  ) = repmat(vgh(:,nghost+1,:),[1,nghost,1]);
          vgh(1+nghost:end-nghost,1:nghost,        1+nghost:end-nghost) = repmat(vbd{ixproc+1,iyproc+1,izproc+1}(:,1,:),[1,nghost,1]);
        end
        if ~yperiodic && iyproc==nyprocs-1
          vgh(:,                  end-nghost+1:end,:                  ) = repmat(vgh(:,end-nghost,:),[1,nghost,1]);
          vgh(1+nghost:end-nghost,end-nghost+1:end,1+nghost:end-nghost) = repmat(vbd{ixproc+1,iyproc+1,izproc+1}(:,end,:),[1,nghost,1]);
        end
        if ~zperiodic && izproc==0
          vgh(:,                  :,                  1:nghost        ) = repmat(vgh(:,:,nghost+1),[1,1,nghost]);
          vgh(1+nghost:end-nghost,1+nghost:end-nghost,1:nghost        ) = repmat(vbd{ixproc+1,iyproc+1,izproc+1}(:,:,1),[1,1,nghost]);
        end
        if ~zperiodic && izproc==nzprocs-1
          vgh(:,                  :,                  end-nghost+1:end) = repmat(vgh(:,:,end-nghost),[1,1,nghost]);
          vgh(1+nghost:end-nghost,1+nghost:end-nghost,end-nghost+1:end) = repmat(vbd{ixproc+1,iyproc+1,izproc+1}(:,:,end),[1,1,nghost]);
        end
        
        % Fill them with ghost cell data
        % w
        % 1st faces (6)
        wgh(1:nghost,           1+nghost:end-nghost,1+nghost:end-nghost) = wbd{mod(ixproc-1,nxprocs)+1,iyproc+1,izproc+1}(end-nghost+1:end,:,               :               );
        wgh(end-nghost+1:end,   1+nghost:end-nghost,1+nghost:end-nghost) = wbd{mod(ixproc+1,nxprocs)+1,iyproc+1,izproc+1}(1:nghost,        :,               :               );
        wgh(1+nghost:end-nghost,1:nghost,           1+nghost:end-nghost) = wbd{ixproc+1,mod(iyproc-1,nyprocs)+1,izproc+1}(:,               end-nghost+1:end,:               );
        wgh(1+nghost:end-nghost,end-nghost+1:end,   1+nghost:end-nghost) = wbd{ixproc+1,mod(iyproc+1,nyprocs)+1,izproc+1}(:,               1:nghost,        :               );
        wgh(1+nghost:end-nghost,1+nghost:end-nghost,1:nghost           ) = wbd{ixproc+1,iyproc+1,mod(izproc-1,nzprocs)+1}(:,               :,               end-nghost+1:end);
        wgh(1+nghost:end-nghost,1+nghost:end-nghost,end-nghost+1:end   ) = wbd{ixproc+1,iyproc+1,mod(izproc+1,nzprocs)+1}(:,               :,               1:nghost        );
        % 2nd edges (12)
        wgh(1:nghost,           1:nghost,           1+nghost:end-nghost) = wbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,izproc+1               }(end-nghost+1:end,end-nghost+1:end,:);
        wgh(1:nghost,           1+nghost:end-nghost,1:nghost           ) = wbd{mod(ixproc-1,nxprocs)+1,iyproc+1,               mod(izproc-1,nzprocs)+1}(end-nghost+1:end,:,               end-nghost+1:end);
        wgh(1+nghost:end-nghost,1:nghost,           1:nghost           ) = wbd{ixproc+1,               mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(:,               end-nghost+1:end,end-nghost+1:end);
        wgh(end-nghost+1:end,   end-nghost+1:end,   1+nghost:end-nghost) = wbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,izproc+1               }(1:nghost,        1:nghost,        :               );
        wgh(end-nghost+1:end,   1+nghost:end-nghost,end-nghost+1:end   ) = wbd{mod(ixproc+1,nxprocs)+1,iyproc+1,               mod(izproc+1,nzprocs)+1}(1:nghost,        :,               1:nghost        );
        wgh(1+nghost:end-nghost,end-nghost+1:end,   end-nghost+1:end   ) = wbd{ixproc+1,               mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(:,               1:nghost,        1:nghost        );
        wgh(1:nghost,           end-nghost+1:end,   1+nghost:end-nghost) = wbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,izproc+1               }(end-nghost+1:end,1:nghost,        :               );
        wgh(1:nghost,           1+nghost:end-nghost,end-nghost+1:end   ) = wbd{mod(ixproc-1,nxprocs)+1,iyproc+1,               mod(izproc+1,nzprocs)+1}(end-nghost+1:end,:,               1:nghost        );
        wgh(1+nghost:end-nghost,1:nghost,           end-nghost+1:end   ) = wbd{ixproc+1,               mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(:,               end-nghost+1:end,1:nghost        );
        wgh(end-nghost+1:end,   1:nghost,           1+nghost:end-nghost) = wbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,izproc+1               }(1:nghost,        end-nghost+1:end,:               );
        wgh(end-nghost+1:end,   1+nghost:end-nghost,1:nghost           ) = wbd{mod(ixproc+1,nxprocs)+1,iyproc+1,               mod(izproc-1,nzprocs)+1}(1:nghost,        :,               end-nghost+1:end);
        wgh(1+nghost:end-nghost,end-nghost+1:end,   1:nghost           ) = wbd{ixproc+1,               mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(:,               1:nghost,        end-nghost+1:end);
        % 3rd corners (8)
        wgh(1:nghost,        1:nghost,        1:nghost        ) = wbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(end-nghost+1:end,end-nghost+1:end,end-nghost+1:end);
        wgh(end-nghost+1:end,1:nghost,        1:nghost        ) = wbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(1:nghost,        end-nghost+1:end,end-nghost+1:end);
        wgh(1:nghost,        end-nghost+1:end,1:nghost        ) = wbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(end-nghost+1:end,1:nghost,        end-nghost+1:end);
        wgh(1:nghost,        1:nghost,        end-nghost+1:end) = wbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(end-nghost+1:end,end-nghost+1:end,1:nghost        );
        wgh(end-nghost+1:end,end-nghost+1:end,1:nghost        ) = wbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(1:nghost,        1:nghost,        end-nghost+1:end);
        wgh(end-nghost+1:end,1:nghost,        end-nghost+1:end) = wbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(1:nghost,        end-nghost+1:end,1:nghost        );
        wgh(1:nghost,        end-nghost+1:end,end-nghost+1:end) = wbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(end-nghost+1:end,1:nghost,        1:nghost        );
        wgh(end-nghost+1:end,end-nghost+1:end,end-nghost+1:end) = wbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(1:nghost,        1:nghost,        1:nghost        );
        % Correct if not periodic
        if ~xperiodic && ixproc==0
          wgh(1:nghost,        :,                  :                  ) = repmat(wgh(nghost+1,:,:),[nghost,1,1]);
          wgh(1:nghost,        1+nghost:end-nghost,1+nghost:end-nghost) = repmat(wbd{ixproc+1,iyproc+1,izproc+1}(1,:,:),[nghost,1,1]);
        end
        if ~xperiodic && ixproc==nxprocs-1
          wgh(end-nghost+1:end,:,                  :                  ) = repmat(wgh(end-nghost,:,:),[nghost,1,1]);
          wgh(end-nghost+1:end,1+nghost:end-nghost,1+nghost:end-nghost) = repmat(wbd{ixproc+1,iyproc+1,izproc+1}(end,:,:),[nghost,1,1]);
        end
        if ~yperiodic && iyproc==0
          wgh(:,                  1:nghost,        :                  ) = repmat(wgh(:,nghost+1,:),[1,nghost,1]);
          wgh(1+nghost:end-nghost,1:nghost,        1+nghost:end-nghost) = repmat(wbd{ixproc+1,iyproc+1,izproc+1}(:,1,:),[1,nghost,1]);
        end
        if ~yperiodic && iyproc==nyprocs-1
          wgh(:,                  end-nghost+1:end,:                  ) = repmat(wgh(:,end-nghost,:),[1,nghost,1]);
          wgh(1+nghost:end-nghost,end-nghost+1:end,1+nghost:end-nghost) = repmat(wbd{ixproc+1,iyproc+1,izproc+1}(:,end,:),[1,nghost,1]);
        end
        if ~zperiodic && izproc==0
          wgh(:,                  :,                  1:nghost        ) = repmat(wgh(:,:,nghost+1),[1,1,nghost]);
          wgh(1+nghost:end-nghost,1+nghost:end-nghost,1:nghost        ) = repmat(wbd{ixproc+1,iyproc+1,izproc+1}(:,:,1),[1,1,nghost]);
        end
        if ~zperiodic && izproc==nzprocs-1
          wgh(:,                  :,                  end-nghost+1:end) = repmat(wgh(:,:,end-nghost),[1,1,nghost]);
          wgh(1+nghost:end-nghost,1+nghost:end-nghost,end-nghost+1:end) = repmat(wbd{ixproc+1,iyproc+1,izproc+1}(:,:,end),[1,1,nghost]);
        end
        
        % Fill them with ghost cell data
        % p
        % 1st faces (6)
        pgh(1:nghost,           1+nghost:end-nghost,1+nghost:end-nghost) = pbd{mod(ixproc-1,nxprocs)+1,iyproc+1,izproc+1}(end-nghost+1:end,:,               :               );
        pgh(end-nghost+1:end,   1+nghost:end-nghost,1+nghost:end-nghost) = pbd{mod(ixproc+1,nxprocs)+1,iyproc+1,izproc+1}(1:nghost,        :,               :               );
        pgh(1+nghost:end-nghost,1:nghost,           1+nghost:end-nghost) = pbd{ixproc+1,mod(iyproc-1,nyprocs)+1,izproc+1}(:,               end-nghost+1:end,:               );
        pgh(1+nghost:end-nghost,end-nghost+1:end,   1+nghost:end-nghost) = pbd{ixproc+1,mod(iyproc+1,nyprocs)+1,izproc+1}(:,               1:nghost,        :               );
        pgh(1+nghost:end-nghost,1+nghost:end-nghost,1:nghost           ) = pbd{ixproc+1,iyproc+1,mod(izproc-1,nzprocs)+1}(:,               :,               end-nghost+1:end);
        pgh(1+nghost:end-nghost,1+nghost:end-nghost,end-nghost+1:end   ) = pbd{ixproc+1,iyproc+1,mod(izproc+1,nzprocs)+1}(:,               :,               1:nghost        );
        % 2nd edges (12)
        pgh(1:nghost,           1:nghost,           1+nghost:end-nghost) = pbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,izproc+1               }(end-nghost+1:end,end-nghost+1:end,:);
        pgh(1:nghost,           1+nghost:end-nghost,1:nghost           ) = pbd{mod(ixproc-1,nxprocs)+1,iyproc+1,               mod(izproc-1,nzprocs)+1}(end-nghost+1:end,:,               end-nghost+1:end);
        pgh(1+nghost:end-nghost,1:nghost,           1:nghost           ) = pbd{ixproc+1,               mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(:,               end-nghost+1:end,end-nghost+1:end);
        pgh(end-nghost+1:end,   end-nghost+1:end,   1+nghost:end-nghost) = pbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,izproc+1               }(1:nghost,        1:nghost,        :               );
        pgh(end-nghost+1:end,   1+nghost:end-nghost,end-nghost+1:end   ) = pbd{mod(ixproc+1,nxprocs)+1,iyproc+1,               mod(izproc+1,nzprocs)+1}(1:nghost,        :,               1:nghost        );
        pgh(1+nghost:end-nghost,end-nghost+1:end,   end-nghost+1:end   ) = pbd{ixproc+1,               mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(:,               1:nghost,        1:nghost        );
        pgh(1:nghost,           end-nghost+1:end,   1+nghost:end-nghost) = pbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,izproc+1               }(end-nghost+1:end,1:nghost,        :               );
        pgh(1:nghost,           1+nghost:end-nghost,end-nghost+1:end   ) = pbd{mod(ixproc-1,nxprocs)+1,iyproc+1,               mod(izproc+1,nzprocs)+1}(end-nghost+1:end,:,               1:nghost        );
        pgh(1+nghost:end-nghost,1:nghost,           end-nghost+1:end   ) = pbd{ixproc+1,               mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(:,               end-nghost+1:end,1:nghost        );
        pgh(end-nghost+1:end,   1:nghost,           1+nghost:end-nghost) = pbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,izproc+1               }(1:nghost,        end-nghost+1:end,:               );
        pgh(end-nghost+1:end,   1+nghost:end-nghost,1:nghost           ) = pbd{mod(ixproc+1,nxprocs)+1,iyproc+1,               mod(izproc-1,nzprocs)+1}(1:nghost,        :,               end-nghost+1:end);
        pgh(1+nghost:end-nghost,end-nghost+1:end,   1:nghost           ) = pbd{ixproc+1,               mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(:,               1:nghost,        end-nghost+1:end);
        % 3rd corners (8)
        pgh(1:nghost,        1:nghost,        1:nghost        ) = pbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(end-nghost+1:end,end-nghost+1:end,end-nghost+1:end);
        pgh(end-nghost+1:end,1:nghost,        1:nghost        ) = pbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(1:nghost,        end-nghost+1:end,end-nghost+1:end);
        pgh(1:nghost,        end-nghost+1:end,1:nghost        ) = pbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(end-nghost+1:end,1:nghost,        end-nghost+1:end);
        pgh(1:nghost,        1:nghost,        end-nghost+1:end) = pbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(end-nghost+1:end,end-nghost+1:end,1:nghost        );
        pgh(end-nghost+1:end,end-nghost+1:end,1:nghost        ) = pbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(1:nghost,        1:nghost,        end-nghost+1:end);
        pgh(end-nghost+1:end,1:nghost,        end-nghost+1:end) = pbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(1:nghost,        end-nghost+1:end,1:nghost        );
        pgh(1:nghost,        end-nghost+1:end,end-nghost+1:end) = pbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(end-nghost+1:end,1:nghost,        1:nghost        );
        pgh(end-nghost+1:end,end-nghost+1:end,end-nghost+1:end) = pbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(1:nghost,        1:nghost,        1:nghost        );
        % Correct if not periodic
        if ~xperiodic && ixproc==0
          pgh(1:nghost,        :,                  :                  ) = repmat(pgh(nghost+1,:,:),[nghost,1,1]);
          pgh(1:nghost,        1+nghost:end-nghost,1+nghost:end-nghost) = repmat(pbd{ixproc+1,iyproc+1,izproc+1}(1,:,:),[nghost,1,1]);
        end
        if ~xperiodic && ixproc==nxprocs-1
          pgh(end-nghost+1:end,:,                  :                  ) = repmat(pgh(end-nghost,:,:),[nghost,1,1]);
          pgh(end-nghost+1:end,1+nghost:end-nghost,1+nghost:end-nghost) = repmat(pbd{ixproc+1,iyproc+1,izproc+1}(end,:,:),[nghost,1,1]);
        end
        if ~yperiodic && iyproc==0
          pgh(:,                  1:nghost,        :                  ) = repmat(pgh(:,nghost+1,:),[1,nghost,1]);
          pgh(1+nghost:end-nghost,1:nghost,        1+nghost:end-nghost) = repmat(pbd{ixproc+1,iyproc+1,izproc+1}(:,1,:),[1,nghost,1]);
        end
        if ~yperiodic && iyproc==nyprocs-1
          pgh(:,                  end-nghost+1:end,:                  ) = repmat(pgh(:,end-nghost,:),[1,nghost,1]);
          pgh(1+nghost:end-nghost,end-nghost+1:end,1+nghost:end-nghost) = repmat(pbd{ixproc+1,iyproc+1,izproc+1}(:,end,:),[1,nghost,1]);
        end
        if ~zperiodic && izproc==0
          pgh(:,                  :,                  1:nghost        ) = repmat(pgh(:,:,nghost+1),[1,1,nghost]);
          pgh(1+nghost:end-nghost,1+nghost:end-nghost,1:nghost        ) = repmat(pbd{ixproc+1,iyproc+1,izproc+1}(:,:,1),[1,1,nghost]);
        end
        if ~zperiodic && izproc==nzprocs-1
          pgh(:,                  :,                  end-nghost+1:end) = repmat(pgh(:,:,end-nghost),[1,1,nghost]);
          pgh(1+nghost:end-nghost,1+nghost:end-nghost,end-nghost+1:end) = repmat(pbd{ixproc+1,iyproc+1,izproc+1}(:,:,end),[1,1,nghost]);
        end
        
        fout = sprintf('%s/%s_%04d.%05d.mat',dout,filebase,iseq,ichunk);
        save(fout,'ugh','vgh','wgh','pgh');
      end
    end
  end
end