       int     *MEM_BUF) 
 {
   const int blockDim2 = NTHREAD2;
-  const int shMemSize = 10 * (1 << blockDim2);
+  const int shMemSize = 14 * (1 << blockDim2);
   __shared__ int shmem_pool[shMemSize];
 
   const int nWarps2 = blockDim2 - WARP_SIZE2;

{


  /*********** set necessary thread constants **********/

  const int DIM2 = WARP_SIZE2;
  const int DIM  = WARP_SIZE;

  /*********** shared memory distribution **********/

  //  begin,    end,   size
  // -----------------------
  const int stack_sz = (LMEM_STACK_SIZE << SHIFT) << DIM2;  /* stack allocated per thread-block */
  const int nWarps2 = BLOCKDIM2 - WARP_SIZE2;
  int *approxL = lmem + stack_sz + (LMEM_EXTRA_SIZE >> nWarps2) * warpId;

  int *directS = shmem;                              //  0*DIM,  1*DIM,  1*DIM
  int *nodesS  = directS + DIM;                      //  1*DIM, 10*DIM,  9*DIM
  int *prefix  = nodesS  + DIM*8;                    //  9*DIM, 10*DIM,  1*DIM

  const int NJMAX = DIM*3;
  int    *body_list = (int*   )&nodesS   [  DIM]; //  2*DIM,   5*DIM,  2*DIM
  float  *sh_mass   = (float* )&body_list[NJMAX]; //  5*DIM,   6*DIM,  1*DIM
  float3 *sh_pos    = (float3*)&sh_mass  [  DIM]; //  6*DIM,   9*DIM   3*DIM

  int *approxM = approxL;
  int *directM = directS;

        int offset  = prefix[laneId];
        offset     += n_offset - nchild;                                  // convert inclusive into exclusive scan for referencing purpose

        for (int i = n_offset; i < n_offset + n_total; i += WARP_SIZE)         //nullify part of the array that will be filled with children
          nodesM[laneId + i] = 0;                                          //but do not touch those parts which has already been filled

#if 0  /* the following gives different result than then one in else */
