c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine make_set

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     This function allocates space for a set of cells and fills the set     
c     such that communication between cells on different nodes is only
c     nearest neighbor                                                   
c---------------------------------------------------------------------

      include 'header.h'
      include 'mpinpb.h'


      integer p, i, j, c, dir, size, excess, ierr,ierrcode

c---------------------------------------------------------------------
c     compute square root; add small number to allow for roundoff
c     (note: this is computed in setup_mpi.f also, but prefer to do
c     it twice because of some include file problems).
c---------------------------------------------------------------------
      ncells = dint(dsqrt(dble(no_nodes) + 0.00001d0))

c---------------------------------------------------------------------
c     this makes coding easier
c---------------------------------------------------------------------
      p = ncells
      
c---------------------------------------------------------------------
c     determine the location of the cell at the bottom of the 3D 
c     array of cells
c---------------------------------------------------------------------
      cell_coord(1,1) = mod(node,p) 
      cell_coord(2,1) = node/p 
      cell_coord(3,1) = 0

c---------------------------------------------------------------------
c     set the cell_coords for cells in the rest of the z-layers; 
c     this comes down to a simple linear numbering in the z-direct-
c     ion, and to the doubly-cyclic numbering in the other dirs     
c---------------------------------------------------------------------
      do c=2, p
         cell_coord(1,c) = mod(cell_coord(1,c-1)+1,p) 
         cell_coord(2,c) = mod(cell_coord(2,c-1)-1+p,p) 
         cell_coord(3,c) = c-1
      end do

c---------------------------------------------------------------------
c     offset all the coordinates by 1 to adjust for Fortran arrays
c---------------------------------------------------------------------
      do dir = 1, 3
         do c = 1, p
            cell_coord(dir,c) = cell_coord(dir,c) + 1
         end do
      end do
      
c---------------------------------------------------------------------
c     slice(dir,n) contains the sequence number of the cell that is in
c     coordinate plane n in the dir direction
c---------------------------------------------------------------------
      do dir = 1, 3
         do c = 1, p
            slice(dir,cell_coord(dir,c)) = c
         end do
      end do


c---------------------------------------------------------------------
c     fill the predecessor and successor entries, using the indices 
c     of the bottom cells (they are the same at each level of k 
c     anyway) acting as if full periodicity pertains; note that p is
c     added to those arguments to the mod functions that might
c     otherwise return wrong values when using the modulo function
c---------------------------------------------------------------------
      i = cell_coord(1,1)-1
      j = cell_coord(2,1)-1

      predecessor(1) = mod(i-1+p,p) + p*j
      predecessor(2) = i + p*mod(j-1+p,p)
      predecessor(3) = mod(i+1,p) + p*mod(j-1+p,p)
      successor(1)   = mod(i+1,p) + p*j
      successor(2)   = i + p*mod(j+1,p)
      successor(3)   = mod(i-1+p,p) + p*mod(j+1,p)

c---------------------------------------------------------------------
c     now compute the sizes of the cells                                    
c---------------------------------------------------------------------
      do dir= 1, 3
c---------------------------------------------------------------------
c     set cell_coord range for each direction                            
c---------------------------------------------------------------------
         size   = grid_points(dir)/p
         excess = mod(grid_points(dir),p)
         do c=1, ncells
            if (cell_coord(dir,c) .le. excess) then
               cell_size(dir,c) = size+1
               cell_low(dir,c) = (cell_coord(dir,c)-1)*(size+1)
               cell_high(dir,c) = cell_low(dir,c)+size
            else 
               cell_size(dir,c) = size
               cell_low(dir,c)  = excess*(size+1)+
     >              (cell_coord(dir,c)-excess-1)*size
               cell_high(dir,c) = cell_low(dir,c)+size-1
            endif
            if (cell_size(dir, c) .le. 2) then
               write(*,50)
 50            format(' Error: Cell size too small. Min size is 3')
               call MPI_Abort(mpi_comm_world,ierrcode,ierr)
               stop
            endif
         end do
      end do

      return
      end

c---------------------------------------------------------------------
c---------------------------------------------------------------------