       int c_stack0 = n_stack_pre;
 
       /*********** walk a level **********/
-      while(c_stack0 < n_stack0) {
-
-	/***
-	**** --> fetch the list of nodes rom LMEM
-	***/
-	bool use_node = tid <  n_nodes0;
- 	{ prefix[tid] = nstack[ACCS<SHIFT>(c_stack0)];   c_stack0++; }
-	__syncthreads();
-	int node  = prefix[min(tid, n_nodes0 - 1)];
+      while(c_stack0 < n_stack0) 
+      {
+
+        /***
+         **** --> fetch the list of nodes rom LMEM
+         ***/
+        bool use_node = tid <  n_nodes0;
+        { prefix[tid] = nstack[ACCS<SHIFT>(c_stack0)];   c_stack0++; }
+        __syncthreads();
+        int node  = prefix[min(tid, n_nodes0 - 1)];
 
         if(n_nodes0 > 0){       //Work around pre 4.1 compiler bug
           n_nodes0 -= DIM;
         }
 
-	/***
-	**** --> process each of the nodes in the list in parallel
-	***/
-
-        #ifndef TEXTURES
-          float4 nodeSize = boxSizeInfo[node];                   //Fetch the size of the box. Size.w = child info
-          float4 node_pos = boxCenterInfo[node];                 //Fetch the center of the box. center.w = opening info
-        #else
-          float4 nodeSize =  tex1Dfetch(texNodeSize, node);
-          float4 node_pos =  tex1Dfetch(texNodeCenter, node);
-        #endif
+        /***
+         **** --> process each of the nodes in the list in parallel
+         ***/
+
+#ifndef TEXTURES
+        float4 nodeSize = boxSizeInfo[node];                   //Fetch the size of the box. Size.w = child info
+        float4 node_pos = boxCenterInfo[node];                 //Fetch the center of the box. center.w = opening info
+#else
+        float4 nodeSize =  tex1Dfetch(texNodeSize, node);
+        float4 node_pos =  tex1Dfetch(texNodeCenter, node);
+#endif
 
         int node_data = __float_as_int(nodeSize.w);
 
-        #ifdef INDSOFT
-          //Very inefficient this but for testing I have to live with it...
-          float node_eps_val = multipole_data[node*3 + 1].w;
-        #endif
-
-	//Check if a cell has to be opened
-	#ifdef IMPBH
-	  //Improved barnes hut method
-          #ifndef TEXTURES
-            float4 nodeCOM = multipole_data[node*3];
-          #else
-            float4 nodeCOM = tex1Dfetch(texMultipole,node*3);
-          #endif
-
-	  nodeCOM.w      = node_pos.w;
-          #ifdef INDSOFT
-            bool   split   = split_node_grav_impbh(nodeCOM, group_pos, groupSize, group_eps, node_eps_val);
-          #else
-            bool   split   = split_node_grav_impbh(nodeCOM, group_pos, groupSize);
-          #endif
-
-	#else
-          //Minimum distance method
-          #ifdef INDSOFT
-            bool   split   = split_node_grav_md(node_pos, nodeSize, group_pos, groupSize, group_eps, node_eps_val);  //Check if node should be split
-          #else
-            bool   split   = split_node_grav_md(node_pos, nodeSize, group_pos, groupSize);  //Check if node should be split
-          #endif
-	#endif
+        //Check if a cell has to be opened
+#ifndef TEXTURES
+        float4 nodeCOM = multipole_data[node*3];
+#else
+        float4 nodeCOM = tex1Dfetch(texMultipole,node*3);
+#endif
+
+        nodeCOM.w      = node_pos.w;
+        bool   split   = split_node_grav_impbh(nodeCOM, group_pos, groupSize);
+
 
         bool leaf       = node_pos.w <= 0;  //Small AND equal incase of a 1 particle cell       //Check if it is a leaf
-//         split = true;
-
-
-	uint mask    = BTEST((split && !leaf) && use_node);               // mask = #FFFFFFFF if use_node+split+not_a_leaf==true, otherwise zero
+        //         split = true;
+
+
+        uint mask    = BTEST((split && !leaf) && use_node);               // mask = #FFFFFFFF if use_node+split+not_a_leaf==true, otherwise zero
         int child    =    node_data & 0x0FFFFFFF;                         //Index to the first child of the node
         int nchild   = (((node_data & 0xF0000000) >> 28)) & mask;         //The number of children this node has
 
 
-  	/***
-	**** --> calculate prefix
-	***/
-
-	int *prefix0 = &prefix[  0];
-	int *prefix1 = &prefix[DIM];
-
-        #ifdef OLDPREFIX
-          int n_total = calc_prefix<DIM2>(prefix, tid,  nchild);
-          prefix[tid] += n_offset - nchild;
-          __syncthreads();
-        #else
-          inclusive_scan_block<ADDOP<int>, int>(prefix, nchild, tid);        // inclusive scan to compute memory offset of each child
-          int n_total = prefix[blockDim.x - 1];                              // fetch total number of children, i.e. offset of the last child -1
-          __syncthreads();                                                   // thread barrier to make sure that warps completed their jobs
-          prefix[tid] += n_offset - nchild;                                  // convert inclusive into exclusive scan for referencing purpose
-          __syncthreads();                                                   // thread barrier
-        #endif
+        /***
+         **** --> calculate prefix
+         ***/
+
+        int *prefix0 = &prefix[  0];
+        int *prefix1 = &prefix[DIM];
+
+#ifdef OLDPREFIX
+        int n_total = calc_prefix<DIM2>(prefix, tid,  nchild);
+        prefix[tid] += n_offset - nchild;
+        __syncthreads();
+#else
+        inclusive_scan_block<ADDOP<int>, int>(prefix, nchild, tid);        // inclusive scan to compute memory offset of each child
+        int n_total = prefix[blockDim.x - 1];                              // fetch total number of children, i.e. offset of the last child -1
+        __syncthreads();                                                   // thread barrier to make sure that warps completed their jobs
+        prefix[tid] += n_offset - nchild;                                  // convert inclusive into exclusive scan for referencing purpose
+        __syncthreads();                                                   // thread barrier
+#endif
 
         for (int i = n_offset; i < n_offset + n_total; i += DIM)         //nullify part of the array that will be filled with children
-            nodes[tid + i] = 0;                                          //but do not touch those parts which has already been filled
+          nodes[tid + i] = 0;                                          //but do not touch those parts which has already been filled
         __syncthreads();                                                 //Thread barrier to make sure all warps finished writing data
 
         bool flag = (split && !leaf) && use_node;                        //Flag = use_node + split + not_a_leaf;Use only non_leaf nodes that are to be split


      node_mon0[tid] = monopole.w;
      node_mon1[tid] = make_float3(monopole.x,  monopole.y,  monopole.z);
      node_oct0[tid] = make_float3(octopole0.x, octopole0.y, octopole0.z);
      node_oct1[tid] = make_float3(octopole1.x, octopole1.y, octopole1.z);

#ifdef INDSOFT
      node_eps[tid]  = octopole0.w;
#endif

    } else {
      //Set non-active memory locations to zero
      node_mon0[tid] = 0.0f;
      node_mon1[tid] = make_float3(1.0e10f, 1.0e10f, 1.0e10f);
      node_oct0[tid] = make_float3(0.0f, 0.0f, 0.0f);
      node_oct1[tid] = make_float3(0.0f, 0.0f, 0.0f);

#ifdef INDSOFT
      node_eps[tid]  = 0.01f;
#endif
    }
    __syncthreads();
#pragma unroll
    for (int i = 0; i < DIMx; i++)
    {
      apprCount++;
#ifdef DOGRAV
#ifdef INDSOFT
      acc_i = add_acc(acc_i, pos_i,
          node_mon0[offs + i], node_mon1[offs + i],
          node_oct0[offs + i], node_oct1[offs + i], node_eps[offs + i], eps2);

#else
      acc_i = add_acc(acc_i, pos_i,
          node_mon0[offs + i], node_mon1[offs + i],
          node_oct0[offs + i], node_oct1[offs + i], eps2);
#endif
#endif
    }

    __syncthreads();
#ifdef INDSOFT
    node_eps[tid]  = temp;
    __syncthreads();
#endif
  } //if n_approx > 0

  if(n_direct > 0)
  {
    if (tid < n_direct) {

      float4 posj = body_pos[direct[tid]];

      //     float4 posj  = tex1Dfetch(texBody, direct[tid]);
      sh_mass[tid] = posj.w;
      sh_pos [tid] = make_float3(posj.x, posj.y, posj.z);
      sh_jid [tid] = direct[tid];

#ifdef INDSOFT
      sh_eps [tid] = body_vel[direct[tid]].w;
#endif
    } else {
      sh_mass[tid] = 0.0f;
      sh_pos [tid] = make_float3(1.0e10f, 1.0e10f, 1.0e10f);
      sh_jid [tid] = -1;

#ifdef INDSOFT
      sh_eps [tid] = 0.01f;
#endif
    }

    __syncthreads();
#pragma unroll
    for (int j = 0; j < DIMx; j++) {
      if ((sh_jid[offs + j] >= 0)) {
        int selfGrav = (body_i != sh_jid[offs + j]);
        float ds2 =  1.0e10f;
        direCount++;
#ifdef DOGRAV
#ifdef INDSOFT
        acc_i = add_acc(acc_i, pos_i, sh_mass[offs + j], sh_pos[offs + j], sh_eps[offs + j], ds2, eps2, selfGrav);
#else
        acc_i = add_acc(acc_i, pos_i, sh_mass[offs + j], sh_pos[offs + j], ds2, eps2, selfGrav);
#endif
#endif
        ds2 += selfGrav*1.0e10f;
        if (ds2 < ds2_min) {
          ngb     = sh_jid[offs + j];
          ds2_min = ds2;
        }
      }
    }
    __syncthreads();
  }


    }
  }
  __syncthreads();
#endif

  return acc_i;
}
