function [ibegu,iendu,jbegu,jendu,kbegu,kendu,...
          ibegv,iendv,jbegv,jendv,kbegv,kendv,...
          ibegw,iendw,jbegw,jendw,kbegw,kendw,...
          ibegp,iendp,jbegp,jendp,kbegp,kendp]=generate_procgrid(...
              nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,...
              nxp,nyp,nzp,nxprocs,nyprocs,nzprocs)
  %[ibegu,iendu,jbegu,jendu,kbegu,kendu,...
  % ibegv,iendv,jbegv,jendv,kbegv,kendv,...
  % ibegw,iendw,jbegw,jendw,kbegw,kendw,...
  % ibegp,iendp,jbegp,jendp,kbegp,kendp]=generate_procgrid(...
  %     nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,...
  %     nxp,nyp,nzp,nxprocs,nyprocs,nzprocs)
  % Generates a processore grid for a given grid.
  % Input
  %   nxu,nyu,nzu,...     staggered grids
  %   nxprocs,...         number of processors
  
  chi='_';
  % U:
  [ibegu,iendu]=mpe_decomp1d(nxu,nxprocs,chi);
  [jbegu,jendu]=mpe_decomp1d(nyu,nyprocs,chi);
  [kbegu,kendu]=mpe_decomp1d(nzu,nzprocs,chi);
  % V:
  [ibegv,iendv]=mpe_decomp1d(nxv,nxprocs,chi);
  [jbegv,jendv]=mpe_decomp1d(nyv,nyprocs,chi);
  [kbegv,kendv]=mpe_decomp1d(nzv,nzprocs,chi);
  % W:
  [ibegw,iendw]=mpe_decomp1d(nxw,nxprocs,chi);
  [jbegw,jendw]=mpe_decomp1d(nyw,nyprocs,chi);
  [kbegw,kendw]=mpe_decomp1d(nzw,nzprocs,chi);
  % P:
  [ibegp,iendp]=mpe_decomp1d(nxp,nxprocs,chi);
  [jbegp,jendp]=mpe_decomp1d(nyp,nyprocs,chi);
  [kbegp,kendp]=mpe_decomp1d(nzp,nzprocs,chi);
end

function [s,e]=mpe_decomp1d(n,numprocs,chi)
  %
  % determines 1d decomposition (as in fortran dns code)
  %
  for myid=0:numprocs-1
    nlocal  = floor(n/ numprocs);
    s(myid+1) = myid * nlocal + 1;
    deficit = mod(n,numprocs);
    s(myid+1) = s(myid+1) + min(myid,deficit);
    if myid<deficit
      nlocal = nlocal + 1;
    end
    e(myid+1) = s(myid+1) + nlocal - 1;
    if (e(myid+1)> n | myid==numprocs-1)
      e(myid+1) = n;
    end
    if chi=='p'
      %    /* finally: augment all counts by one */
      s(myid+1) = s(myid+1) + 1;
      e(myid+1) = e(myid+1) + 1;
    end
  end
end