function [nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,...
          xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp,dx,dy,dz]=...
            generate_grid(a,b,c,d,e,f,nxp,nyp,nzp,xperiodic,yperiodic,zperiodic)
  % [nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,...
  %  xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp,dx,dy,dz]=...
  %          generate_grid(a,b,c,d,e,f,nxp,nyp,nzp,xperiodic,yperiodic,zperiodic)
  % Generates a staggered grid.
  % Input
  %   a,b,c,d,e,f       domain bounds
  %   nxp,nyp,nzp       number of points for pressure grid
  %   x/y/zperiodic     domain periodicity
  
  if xperiodic==1
    nxu=nxp;
    nxv=nxp;
    nxw=nxp;
  else
    nxu=nxp+1;
    nxv=nxp;
    nxw=nxp;
  end
  if yperiodic==1
    nyu=nyp;
    nyv=nyp;
    nyw=nyp;
  else
    nyu=nyp;
    nyv=nyp+1;
    nyw=nyp;
  end
  if zperiodic==1
    nzu=nzp;
    nzv=nzp;
    nzw=nzp;
  else
    nzu=nzp;
    nzv=nzp;
    nzw=nzp+1;
  end
  %grid step:
  if xperiodic==1
    dx=(b-a)/(nxp);
  else
    dx=(b-a)/(nxp-1);
  end
  if yperiodic==1
    dy=(d-c)/(nyp);
  else
    dy=(d-c)/(nyp-1);
  end
  if zperiodic==1
    dz=(f-e)/(nzp);
  else
    dz=(f-e)/(nzp-1);
  end
  %
  if xperiodic==1
    xu=a+((1:nxu)-1)*dx;
  else
    xu=a+((1:nxu)-3/2)*dx;
  end
  if yperiodic==1
    yu=c+((1:nyu)-1/2)*dy;
  else
    yu=c+((1:nyu)-1)*dy;
  end
  if zperiodic==1
    zu=e+((1:nzu)-1/2)*dz;
  else
    zu=e+((1:nzu)-1)*dz;
  end
  %
  if xperiodic==1
    xv=a+((1:nxv)-1/2)*dx;
  else
    xv=a+((1:nxv)-1)*dx;
  end
  if yperiodic==1
    yv=c+((1:nyv)-1)*dy;
  else
    yv=c+((1:nyv)-3/2)*dy;
  end
  if zperiodic==1
    zv=e+((1:nzv)-1/2)*dz;
  else
    zv=e+((1:nzv)-1)*dz;
  end
  %
  if xperiodic==1
    xw=a+((1:nxw)-1/2)*dx;
  else
    xw=a+((1:nxw)-1)*dx;
  end
  if yperiodic==1
    yw=c+((1:nyw)-1/2)*dy;
  else
    yw=c+((1:nyw)-1)*dy;
  end
  if zperiodic==1
    zw=e+((1:nzw)-1)*dz;
  else
    zw=e+((1:nzw)-3/2)*dz;
  end
  %
  if xperiodic==1
    xp=a+((1:nxp)-1/2)*dx;
  else
    xp=a+((1:nxp)-1)*dx;
  end
  if yperiodic==1
    yp=c+((1:nyp)-1/2)*dy;
  else
    yp=c+((1:nyp)-1)*dy;
  end
  if zperiodic==1
    zp=e+((1:nzp)-1/2)*dz;
  else
    zp=e+((1:nzp)-1)*dz;
  end
end