           __syncthreads();
 
           //step 1:
+#if 0
           if (flag && (direct[tid] <= nb) && (prefix0[tid] > 0))        //make sure that the thread indeed carries a leaf
             body_list[n_direct + prefix0[tid] - 1] = 1;                 //whose bodies will be extracted
           __syncthreads();
+#else
+          if (flag && (direct[tid] <= nb) && (offset1 > 0))        //make sure that the thread indeed carries a leaf
+            body_list[n_direct + offset1 - 1] = 1;                 //whose bodies will be extracted
+          __syncthreads();
+#endif
 
           //step 2:
-#if 1 //def OLDPREFIX
+#if 0 //def OLDPREFIX
           int nl = calc_prefix<DIM2>(nb, &body_list[n_direct], tid);
 #else
-          int nl = inclusive_scan_array<ADDOP<int>, int>              // inclusive scan to compute number of leaves to process
-            (&body_list[n_direct], nb, tid);            // to make sure that there is enough shared memory for bodies
+          int nl = inclusive_scan_array<DIM2>              // inclusive scan to compute number of leaves to process
+            (&body_list[n_direct], prefix_shmem, nb, tid);            // to make sure that there is enough shared memory for bodies
 #endif
           nb = direct[prefix1[nl_pre + nl - 1]];                        // number of bodies stored in these leaves
 

        sh_body  [tid] = jbody;                                                  //store the leafs first body id into shared memory

        // step 1
#ifdef OLDPREFIX
        calc_prefix<DIM2>(prefix0, tid, flag);
#else
#if 0
        inclusive_scan_block<ADDOP<int>, int, DIM2>(prefix0, (int)flag, tid);       // inclusive scan on flags to construct array
#else
        v0 = inclusive_scan_block<ADDOP<int>, int, DIM2>((int)flag, tid);       // inclusive scan on flags to construct array
        prefix0[tid] = v0;
        __syncthreads();
#endif
#endif

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

        // step 2
#ifdef OLDPREFIX
        int n_bodies  = calc_prefix<DIM2>(prefix0, tid, nbody);
#else
#if 0
        inclusive_scan_block<ADDOP<int>, int, DIM2>(prefix0, nbody, tid);        // inclusive scan to compute memory offset for each body
#else
        v0 = inclusive_scan_block<ADDOP<int>, int, DIM2>(nbody, tid);        // inclusive scan to compute memory offset for each body
        prefix[tid] = v0;
        __syncthreads();
#endif
        int n_bodies = prefix0[blockDim.x - 1];                            //Total number of bides extract from the leaves
        __syncthreads();                                                   // thread barrier to make sure that warps completed their jobs
#endif

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

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


        if (flag && nodes[prefix[tid] + 5] == 0) nodes[prefix[tid] + 5] = child + 5; __syncthreads();
        if (flag && nodes[prefix[tid] + 6] == 0) nodes[prefix[tid] + 6] = child + 6; __syncthreads();
        if (flag && nodes[prefix[tid] + 7] == 0) nodes[prefix[tid] + 7] = child + 7; __syncthreads();
#else
        /* Thread with the node that is about to be split */
        if (flag                          ) nodes[offset    ] = child;    
        /* in the following 8 lines, we calculate indexes of all the children that have to be walked from the index of the first child */
        if (flag && nodes[offset + 1] == 0) nodes[offset + 1] = child + 1;
        if (flag && nodes[offset + 2] == 0) nodes[offset + 2] = child + 2;
        if (flag && nodes[offset + 3] == 0) nodes[offset + 3] = child + 3;
        if (flag && nodes[offset + 4] == 0) nodes[offset + 4] = child + 4;
        if (flag && nodes[offset + 5] == 0) nodes[offset + 5] = child + 5;
        if (flag && nodes[offset + 6] == 0) nodes[offset + 6] = child + 6;
        if (flag && nodes[offset + 7] == 0) nodes[offset + 7] = child + 7;
        __syncthreads();
#endif

        n_offset += n_total;    //Increase the offset in the array by the number of newly added nodes

