       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
+#ifdef INDSOFT
+      node_eps[tid]  = 0.01f;
+#endif
     }
     __syncthreads();
-  #pragma unroll
+#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
+#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
+#ifdef INDSOFT
+    node_eps[tid]  = temp;
+    __syncthreads();
+#endif
   } //if n_approx > 0
 
   if(n_direct > 0)

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

          //float darkMatterMass = octopole1.w;
          //float stellarMass = monopole.w - darkMatterMass;

          #ifdef INDSOFT
            float temp          = node_eps[tid]; //Backup value in the shmem into register
            node_eps[tid]       = octopole0.w;
          #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; //Restore original value in shmem
             __syncthreads();
          #endif
 	}
	__syncthreads();
#endif

#if 1
	/***********************************/
	/******       DIRECT          ******/
	/***********************************/

        int *sh_body = &approx[DIM];

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

        body_list[tid] = direct[tid];                                            //copy list of bodies from previous pass to body_list
        sh_body  [tid] = jbody;                                                  //store the leafs first body id into shared memory

        // step 1
        #ifdef OLDPREFIX
          calc_prefix<DIM2>(prefix0, tid, flag);
        #else
          inclusive_scan_block<ADDOP<int>, int>(prefix0, (int)flag, tid);       // inclusive scan on flags to construct array
        #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
          inclusive_scan_block<ADDOP<int>, int>(prefix0, nbody, tid);        // inclusive scan to compute memory offset for each body
          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

        #define NJMAX (DIM*4)
	while (n_bodies > 0) {
	  int nb    = min(n_bodies, NJMAX - n_direct);                    //Make sure number of bides to be extracted does not exceed
                                                                          //the amount of allocated shared memory

          // step 0                                                      //nullify part of the body_list that will be filled with bodies
	  for (int i = n_direct; i < n_direct + nb; i += DIM){           //from the leaves that are being processed
            body_list[i + tid] = 0;
          }
	  __syncthreads();

          //step 1:
	  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();

          //step 2:
         #ifdef 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
          #endif
	  nb = direct[prefix1[nl_pre + nl - 1]];                        // number of bodies stored in these leaves

	  // step 3:
	  for (int i = n_direct; i < n_direct + nb; i += DIM) {          //segmented fill of the body_list
	    int j = prefix1[nl_pre + body_list[i + tid] - 1];            // compute the first body in shared j-body array
	    body_list[i + tid] = (i + tid - n_direct) -                 //add to the index of the first j-body in a child
                                 (prefix0[j] - 1) + sh_body[j];         //the index of the first child in body_list array
	  }
	  __syncthreads();


         /**************************************************
          *  example of what is accomplished in steps 0-4   *
          *       ---------------------------               *
          * step 0: body_list = 000000000000000000000       *
          * step 1: body_list = 100010001000000100100       *
          * step 2: body_list = 111122223333333444555       *
          * step 3: body_list = 012301230123456012012       *
          *         assuming that sh_body[j] = 0            *
         ***************************************************/

	  n_bodies     -= nb;                                   //subtract from n_bodies number of bodies that have been extracted
   	  nl_pre       += nl;                                   //increase the number of leaves that where processed
	  direct [tid] -= nb;                                   //subtract the number of extracted bodies in this pass
 	  prefix0[tid] = max(prefix0[tid] - nb, 0);             //same here, but do not let the number be negative (GT200 bug!?)
	  n_direct     += nb;                                  //increase the number of bodies to be procssed

	  while(n_direct >= DIM) {
	    n_direct -= DIM;


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

            #ifdef INDSOFT
              float   temp = sh_eps [tid];        //Store the value from Shmem into a register
              sh_eps [tid] = body_vel[body_list[n_direct + tid]].w;  //Load the softening
            #endif
	    __syncthreads();
#pragma unroll
	    for (int j = 0; j < DIMx; j++)
            {
              int selfGrav = (body_i != sh_jid[offs + j]);
// 	      if (body_i != sh_jid[offs + j]) //If statement replaced by multiplication
              {
		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;
		}
	      }
	    }
            #ifdef INDSOFT
              //Restore the shmem value after all threads have used it
              __syncthreads();
              sh_eps[tid] = temp;
            #endif
	    __syncthreads();
	  }

	}
	direct[tid] = body_list[tid];
	__syncthreads();
#endif
      } //end lvl



  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];
      float4 octopole1 = multipole_data[address + 2];
#else
      float4 monopole  = tex1Dfetch(texMultipole, address);
      float4 octopole0 = tex1Dfetch(texMultipole, address + 1);
      float4 octopole1 = tex1Dfetch(texMultipole, address + 2);
#endif

      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
