     return calc_prefix1<DIM2>(prefix, tid, value);
   else
   {
-    const int WARP_SIZE2 = 5;
-    const int WARP_SIZE  = 1 << WARP_SIZE2;
     prefix[tid] = inclusive_scan_warp<WARP_SIZE2>(value);
     __syncthreads();
 


  //  begin,    end,   size
  // -----------------------
#ifdef _ORIG_SHMEM_

  int *approxS = (int*)&shmem  [     0];            //  0*DIM,  2*DIM,  2*DIM
  int *directS = (int*)&approxS[ 2*DIM];            //  2*DIM,  3*DIM,  1*DIM
  int *nodesS = (int*)&directS [   DIM];            //  3*DIM, 12*DIM,  9*DIM
  int *prefix = (int*)&nodesS  [9 *DIM];            // 12*DIM, 14*DIM,  2*DIM
  int *sh_body = &approxS[DIM];

  int *prefix0 = &prefix[  0];
  int *prefix1 = &prefix[DIM];

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

  int *approxM = approxS;
  int *directM = directS;
  int * nodesM =  nodesS;

#else   /* !_ORIG_SHMEM_ */

  const int stack_sz = (LMEM_STACK_SIZE << SHIFT) << DIM2;  /* stack allocated per thread-block */
  int *approxL = lmem + stack_sz; 

  int *directS = shmem;                              //  0*DIM,  1*DIM,  1*DIM
  int *nodesS  = directS + DIM;                      //  1*DIM, 10*DIM,  9*DIM
  int *prefixS = nodesS  + DIM*9;                    // 10*DIM, 12*DIM,  2*DIM

  int *prefix  = prefixS;
  int *prefix0 = prefixS;
  int *prefix1 = prefixS + 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    *sh_body   = nodesS + DIM*8;             //  9*DIM,  10*DIM,  1*DIM

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

#endif /* _ORIG_SHMEM_ */


  float  *node_mon0 = sh_mass;

        /***********************************/


        flag            = split && leaf && use_node;                                //flag = split + leaf + use_node
        const int jbody = node_data & BODYMASK;                                     //the first body in the leaf
        const int nbody = (((node_data & INVBMASK) >> LEAFBIT)+1) & BTEST(flag);    //number of bodies in the leaf masked with the flag

        body_list[tid] = directM[tid];                                            //copy list of bodies from previous pass to body_list

        // step 1
        calc_prefix<DIM2>(prefix, tid, flag);                                // inclusive scan on flags to construct array
        const int offset1 = prefix[tid];
        
        // step 2
        int n_bodies  = calc_prefix<DIM2>(prefix, tid, nbody);              // inclusive scan to compute memory offset for each body
        offset = prefix[tid];
        __syncthreads();

        if (flag) prefix[offset1 - 1] = tid;                             //with tidś whose leaves have to be opened
        __syncthreads();                                                      //thread barrier, make sure all warps completed the job

        directM[tid]  = offset;                                       //Store a copy of inclusive scan in direct
        offset       -= nbody;                                              //convert inclusive int oexclusive scan
        offset       += 1;                                                  //add unity, since later prefix0[tid] == 0 used to check barrier

        int nl_pre = 0;                                                     //Number of leaves that have already been processed

