       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


  if(n_approx > 0)
  {
#ifdef INDSOFT
    float temp = node_eps[tid];
#endif

    if (tid < n_approx) {
      int address      = (approx[tid] << 1) + approx[tid];
#ifndef TEXTURES
      float4 monopole  = multipole_data[address    ];
      float4 octopole0 = multipole_data[address + 1];


template<int DIM2, int SHIFT>
__device__ float4 approximate_gravity(int DIM2x, int DIM2y,
    int tid, int tx, int ty,
    int body_i, float4 pos_i,
    real4 group_pos,
    float eps2,
    uint2 node_begend,
    real4 *multipole_data,
    real4 *body_pos,
    int *shmem,
    int *lmem,
    int &ngb,
    int &apprCount, int &direCount,
    volatile float4 *boxSizeInfo,
    float4 groupSize,
    volatile float4 *boxCenterInfo,
    float group_eps,
    real4 *body_vel) {

  float4 acc_i = {0.0f, 0.0f, 0.0f, 0.0f};
  ngb = -1;
