diff --git a/phigrape.cpp b/phigrape.cpp
index c747a9e..33c243f 100644
--- a/phigrape.cpp
+++ b/phigrape.cpp
@@ -58,19 +58,19 @@ Last redaction : 2019.04.16 12:55
 
 //#define ADD_BH1        // add the Single BH
 
-// #define ADD_BH2            // add the Binary BH's
-// #define ADD_N_BH        // eps_BH = 0.0, but added only the Newtonian forces
-// #define ADD_PN_BH        // extra - added also the Post-Newton forces
-// #define ADD_SPIN_BH        // extra - added the SPIN for the BH's - DEFAULT !!!
+//#define ADD_BH2            // add the Binary BH's
+//#define ADD_N_BH        // eps_BH = 0.0, but added only the Newtonian forces
+//#define ADD_PN_BH        // extra - added also the Post-Newton forces
+//#define ADD_SPIN_BH        // extra - added the SPIN for the BH's - DEFAULT !!!
 
-// #define BH_OUT           // extra output for BH's (live)
-// #define BH_OUT_NB        // extra output for the BH's neighbours (live)
+//#define BH_OUT           // extra output for BH's (live)
+//#define BH_OUT_NB        // extra output for the BH's neighbours (live)
 
-// #define BBH_INF            // BBH influence sphere...
-// #define R_INF  10.0        // Factor for the influence sphere... if ( R < R_INF * DR_BBH )
-// #define R_INF2 (R_INF*R_INF)
+//#define BBH_INF            // BBH influence sphere...
+//#define R_INF  10.0        // Factor for the influence sphere... if ( R < R_INF * DR_BBH )
+//#define R_INF2 (R_INF*R_INF)
 
-// #define DT_MIN_WARNING         // dt < dt_min warning !!!
+//#define DT_MIN_WARNING         // dt < dt_min warning !!!
 
 //#define BH_OUT_NB_EXT        // extra output for the BH's neighbours (external)
 
@@ -137,7 +137,7 @@ Last redaction : 2019.04.16 12:55
 
 //#define AMD
 
-// #define ACT_DEF_LL
+//#define ACT_DEF_LL
 
 #if defined(ACT_DEF_LL) && defined(ACT_DEF_GRAPITE)
 #error "Contradicting preprocessor flags!"
@@ -310,7 +310,6 @@ FILE   *inp, *out, *tmp_file, *dbg;
 double CPU_time_real0, CPU_time_user0, CPU_time_syst0;
 double CPU_time_real,  CPU_time_user,  CPU_time_syst;
 
-
 #ifdef TIMING
 
 double CPU_tmp_real0, CPU_tmp_user0, CPU_tmp_syst0;
@@ -404,12 +403,10 @@ double r_sc[M_SC_DIM], m_sc[M_SC_DIM], p_sc[M_SC_DIM];
 double M_R, pot_out_R;
 #endif
 
-
 double pot_ext[N_MAX], pot_act_ext[N_MAX];    // for all EXTPOT
 
 #endif
 
-
 int    i_bh, i_bh1, i_bh2,
     num_bh = 0, num_bh1 = 0, num_bh2 = 0;
 
@@ -418,7 +415,6 @@ double m_bh, m_bh1, m_bh2, b_bh,
     x_ij, y_ij, z_ij,
     vx_ij, vy_ij, vz_ij, rv_ij;
 
-
 #ifdef BBH_INF
 int    inf_event[N_MAX];
 double x_bbhc[3], v_bbhc[3], DR2, tmp_r2;
@@ -581,7 +577,6 @@ struct rusage xxx;
 double sec_u, microsec_u, sec_s, microsec_s;
 struct timeval tv;
 
-
 getrusage(RUSAGE_SELF,&xxx);
 
 sec_u = xxx.ru_utime.tv_sec;
@@ -666,7 +661,6 @@ fclose(out);
 }
 /****************************************************************************/
 
-
 /****************************************************************************/
 void write_bh_data()
 {
@@ -764,7 +758,6 @@ fprintf(out,"%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .
 
 #endif        //    ADD_N_BH
 
-
 /* 2nd BH */
 
 //  i=N-1;
@@ -860,7 +853,6 @@ fclose(out);
 
 #else
 
-
 #ifdef ADD_BH1
 
 out = fopen("bh.dat","a");
@@ -1111,8 +1103,6 @@ fclose(out);
 #endif
 /****************************************************************************/
 
-
-
 /****************************************************************************/
 void calc_self_grav_zero()
 {
@@ -1198,7 +1188,6 @@ for (i=0; i<ni; i+=npipe)
 
     } /* i */
 
-
 /* Store the value of the local partial force etc... */
 
 for (i=0; i<n_act; i++)
@@ -1257,7 +1246,6 @@ for (i=0; i<M_SC; i++)
 #endif
 /****************************************************************************/
 
-
 /****************************************************************************/
 #ifdef EXTPOT_SC
 void def_SC_mass(double r, double *M_R, double *pot_out_R)
@@ -1272,14 +1260,12 @@ qqq = r*M_SC/r_sc[M_SC-1];
 ii1 = qqq - 1;
 ii2 = qqq + 1;
 
-
 //    if (myRank == rootRank)
 //      {
 //      printf("\t \t %06d %06d %.6E %.6E %.6E \n", ii1, ii2, r, M_R, pot_out_R);
 //      fflush(stdout);
 //      } /* if (myRank == rootRank) */
 
-
 for (ii=ii1; ii<ii2; ii++)
     {
     if ( (r_sc[ii] < r) && (r_sc[ii+1] > r) )
@@ -1327,7 +1313,6 @@ for (i=0; i<ni; i++)
     pot_ext[i] = 0.0;
     }
 
-
 #ifdef EXTPOT_GAL
 
 for (i=0; i<ni; i++)
@@ -1354,7 +1339,6 @@ for (i=0; i<ni; i++)
     y2_ij = SQR(y_ij);
     z2_ij = SQR(z_ij);
 
-
     /* for BULGE */
 
     r2 = SQR(b_bulge);
@@ -1377,7 +1361,6 @@ for (i=0; i<ni; i++)
     adot[i][1] -= tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
     adot[i][2] -= tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);
 
-
     /* for DISK */
 
     z2_tmp = z2_ij + SQR(b_disk);
@@ -1415,7 +1398,6 @@ for (i=0; i<ni; i++)
                         + 3.0*vx_ij*x_ij*z_ij*z2_tmp*(z_tmp + a_disk)
                         + 3.0*vy_ij*y_ij*z_ij*z2_tmp*(z_tmp + a_disk) ) / z2_tmp ;
 
-
     /* for HALO */
 
     r2 = SQR(b_halo);
@@ -1442,7 +1424,6 @@ for (i=0; i<ni; i++)
 
 #endif
 
-
 #ifdef EXTPOT_GAL_DEH
 
 for (i=0; i<ni; i++)
@@ -1456,7 +1437,6 @@ for (i=0; i<ni; i++)
     vy_ij = v_act[i][1];
     vz_ij = v_act[i][2];
 
-
     tmp_r = sqrt( SQR(x_ij) + SQR(y_ij) + SQR(z_ij) );
     tmp_r2 = SQR(tmp_r);
     tmp_r3 = POW3(tmp_r);
@@ -1466,17 +1446,14 @@ for (i=0; i<ni; i++)
     dum3  = POW3(dum);
     dum_g = pow( dum , -g_ext );
 
-
     pot_ext[i] -= ( (m_ext/r_ext) / (2.0-g_ext) ) * ( 1.0 - dum2 * dum_g );
 
-
     tmp = (m_ext/tmp_r3) * dum3 * dum_g;
 
     a[i][0] -= tmp * x_ij;
     a[i][1] -= tmp * y_ij;
     a[i][2] -= tmp * z_ij;
 
-
     rv_ij = ( vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij );
 
     tmp = ( m_ext/POW3(tmp_r+r_ext) ) * dum_g;
@@ -1488,11 +1465,8 @@ for (i=0; i<ni; i++)
 
     } /* i */
 
-
 #endif
 
-
-
 /* For simple Log potential... */
 
 #ifdef EXTPOT_GAL_LOG
@@ -1513,17 +1487,14 @@ for (i=0; i<ni; i++)
 
     r2_r2_halo = (r2 + r2_halo);
 
-
     pot_ext[i] -= 0.5*v2_halo * log(1.0 + r2/r2_halo);
 
-
     tmp = v2_halo/r2_r2_halo;
 
     a[i][0] -= tmp * x_ij;
     a[i][1] -= tmp * y_ij;
     a[i][2] -= tmp * z_ij;
 
-
     tmp /= (r2 + r2_halo);
 
     adot[i][0] += tmp * (rv_ij * x_ij - r2_r2_halo * vx_ij);
@@ -1534,8 +1505,6 @@ for (i=0; i<ni; i++)
 
 #endif
 
-
-
 /* For simple Plummer potential... */
 
 #ifdef EXTPOT_BH
@@ -1575,8 +1544,6 @@ for (i=0; i<ni; i++)
 
 #endif
 
-
-
 #ifdef EXTPOT_SC
 
 for (i=0; i<ni; i++)
@@ -1630,15 +1597,11 @@ for (i=0; i<ni; i++)
 
 #endif
 
-
-
 #endif        // EXTPOT
 
-
 }
 /****************************************************************************/
 
-
 /****************************************************************************/
 void calc_self_grav()
 {
@@ -1730,8 +1693,6 @@ for (i=0; i<ni; i+=npipe)
 
     } /* i */
 
-
-
 /* Store the new value of the local partial force etc... */
 
 for (i=0; i<n_act; i++)
@@ -1748,17 +1709,14 @@ for (i=0; i<n_act; i++)
 
     } /* i */
 
-
 #ifdef TIMING
 get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
 DT_ACT_GRAV += (CPU_tmp_user - CPU_tmp_user0);
 #endif
 
-
 }
 /****************************************************************************/
 
-
 /****************************************************************************/
 void calc_ext_grav()
 {
@@ -1774,7 +1732,6 @@ for (i=0; i<ni; i++)
     pot_act_ext[i] = 0.0;
     }
 
-
 #ifdef TIMING
 get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
 #endif
@@ -1827,7 +1784,6 @@ for (i=0; i<ni; i++)
     adot_act_new[i][1] -= tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
     adot_act_new[i][2] -= tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);
 
-
     /* for DISK */
 
     z2_tmp = z2_ij + SQR(b_disk);
@@ -1891,8 +1847,6 @@ for (i=0; i<ni; i++)
 
 #endif
 
-
-
 #ifdef EXTPOT_GAL_DEH
 
 for (i=0; i<ni; i++)
@@ -1906,7 +1860,6 @@ for (i=0; i<ni; i++)
     vy_ij = v_act_new[i][1];
     vz_ij = v_act_new[i][2];
 
-
     tmp_r = sqrt( SQR(x_ij) + SQR(y_ij) + SQR(z_ij) );
     tmp_r2 = SQR(tmp_r);
     tmp_r3 = POW3(tmp_r);
@@ -1916,17 +1869,14 @@ for (i=0; i<ni; i++)
     dum3  = POW3(dum);
     dum_g = pow( dum , -g_ext );
 
-
     pot_act_ext[i] -= ( (m_ext/r_ext) / (2.0-g_ext) ) * ( 1.0 - dum2 * dum_g );
 
-
     tmp = (m_ext/tmp_r3) * dum3 * dum_g;
 
     a_act_new[i][0] -= tmp * x_ij;
     a_act_new[i][1] -= tmp * y_ij;
     a_act_new[i][2] -= tmp * z_ij;
 
-
     rv_ij = ( vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij );
 
     tmp = ( m_ext/POW3(tmp_r+r_ext) ) * dum_g;
@@ -1938,11 +1888,8 @@ for (i=0; i<ni; i++)
 
     } /* i */
 
-
 #endif
 
-
-
 /* For simple Log potential... */
 
 #ifdef EXTPOT_GAL_LOG
@@ -1958,23 +1905,19 @@ for (i=0; i<ni; i++)
     vy_ij = v_act_new[i][1];
     vz_ij = v_act_new[i][2];
 
-
     r2 = SQR(x_ij) + SQR(y_ij) + SQR(z_ij);
     rv_ij = 2.0 * ( vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij );
 
     r2_r2_halo = (r2 + r2_halo);
 
-
     pot_act_ext[i] -= 0.5*v2_halo * log(1.0 + r2/r2_halo);
 
-
     tmp = v2_halo/r2_r2_halo;
 
     a_act_new[i][0] -= tmp * x_ij;
     a_act_new[i][1] -= tmp * y_ij;
     a_act_new[i][2] -= tmp * z_ij;
 
-
     tmp /= (r2 + r2_halo);
 
     adot_act_new[i][0] += tmp * (rv_ij * x_ij - r2_r2_halo * vx_ij);
@@ -1985,9 +1928,6 @@ for (i=0; i<ni; i++)
 
 #endif
 
-
-
-
     /* For simple Plummer potential... */
 
 #ifdef EXTPOT_BH
@@ -2027,8 +1967,6 @@ for (i=0; i<ni; i++)
 
 #endif
 
-
-
 #ifdef EXTPOT_SC
 
 for (i=0; i<ni; i++)
@@ -2076,7 +2014,6 @@ for (i=0; i<ni; i++)
 
 #endif
 
-
 #ifdef TIMING
 get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
 DT_EXT_GRAV += (CPU_tmp_user - CPU_tmp_user0);
@@ -2086,10 +2023,6 @@ DT_EXT_GRAV += (CPU_tmp_user - CPU_tmp_user0);
 }
 /****************************************************************************/
 
-
-
-
-
 /****************************************************************************/
 void energy_contr()
 {
@@ -2113,7 +2046,6 @@ if (Timesteps == 0.0)
     vcm_sum = 0.0;
     }
 
-
 mcm     = 0.0;
 rcm_mod = 0.0;
 vcm_mod = 0.0;
@@ -2142,14 +2074,12 @@ for (k=0;k<3;k++)
     vcm[k] /= mcm;
     } /* k */
 
-
 rcm_mod = sqrt(xcm[0]*xcm[0]+xcm[1]*xcm[1]+xcm[2]*xcm[2]);
 vcm_mod = sqrt(vcm[0]*vcm[0]+vcm[1]*vcm[1]+vcm[2]*vcm[2]);
 
 rcm_sum += rcm_mod;
 vcm_sum += vcm_mod;
 
-
 for (k=0;k<3;k++) mom[k] = 0.0;
 
 for (i=0; i<N; i++)
@@ -2159,7 +2089,6 @@ for (i=0; i<N; i++)
     mom[2] += m[i] * ( x[i][0] * v[i][1] - x[i][1] * v[i][0] );
     } /* i */
 
-
 get_CPU_time(&CPU_time_real, &CPU_time_user, &CPU_time_syst);
 
 E_tot = E_pot + E_kin + E_pot_ext;
@@ -2226,8 +2155,7 @@ printf("%.3E %.3E  % .4E %.4E % .4E  % .4E % .4E  %.2E\n",
         E_pot, E_kin, E_pot_ext, E_tot, DE_tot,
         CPU_time_user-CPU_time_user0);
 
-#endif        // #if defined(STARDESTR) || defined(STARDESTR_EXT)
-
+#endif        //#if defined(STARDESTR) || defined(STARDESTR_EXT)
 
 fflush(stdout);
 
@@ -2251,7 +2179,6 @@ fprintf(out,"%.8E \t %.8E %.8E %.8E \t % .8E % .8E % .8E % .8E % .8E \t % .8E %
         CPU_time_real-CPU_time_real0, CPU_time_user-CPU_time_user0, CPU_time_syst-CPU_time_syst0);
 fclose(out);
 
-
 E_tot_0 = E_tot;
 
 //  printf("3: %.6E %.6E \t % .8E % .8E \n", time_cur, Timesteps, E_tot, E_tot_0);
@@ -2305,19 +2232,15 @@ for (i=0; i<N; i++)
 }
 /****************************************************************************/
 
-
 /****************************************************************************/
 /****************************************************************************/
 /****************************************************************************/
 
 int main(int argc, char *argv[])
 {
-
-
     /* INIT the rand() !!! */
     srand(19640916);                 /* it is just my birthday :-) */
 
-
     /* Init MPI */
     MPI_Init(&argc, &argv);
 
@@ -2338,124 +2261,93 @@ int main(int argc, char *argv[])
 
     /* read the input parameters to the rootRank */
 
-    if (myRank == rootRank)
-        {
-
+    if (myRank == rootRank) {
         inp = fopen("phi-GRAPE.cfg","r");
-
         fscanf(inp,"%lE %lE %lE %lE %lE %lE %s", &eps, &t_end, &dt_disk, &dt_contr, &dt_bh, &eta, inp_fname);
-
         fclose(inp);
 
-
-
-    #ifdef BBH_INF
-        for (i=0; i<N; i++)
-        {
-        inf_event[i] = 0;
-        }
-    #endif
-
-
+#ifdef BBH_INF
+        for (i=0; i<N; i++) inf_event[i] = 0;
+#endif
 
         /*
-        eps         : Plummer softening parameter (can be even 0)
+        eps      : Plummer softening parameter (can be even 0)
         t_end    : end time of calculation
         dt_disk  : interval of snapshot files output (0xxx.dat)
         dt_contr : interval for the energy control output (contr.dat)
         dt_bh    : interval for BH output (bh.dat & bh_nb.dat)
-        eta         : parameter for timestep determination
+        eta      : parameter for timestep determination
         inp_data : name of the input file (data.inp)
         */
 
-
     /* Normalization... */
-    #ifdef NORM
-    inp = fopen("phi-GRAPE.norm","r");
-    fscanf(inp,"%lE %lE", &m_norm, &r_norm);
-    fclose(inp);
+#ifdef NORM
+        inp = fopen("phi-GRAPE.norm","r");
+        fscanf(inp,"%lE %lE", &m_norm, &r_norm);
+        fclose(inp);
 
-    m_norm *= Msol;
-    r_norm *= pc;
+        m_norm *= Msol;
+        r_norm *= pc;
 
-    v_norm = sqrt(G*m_norm/r_norm);
-    t_norm = (r_norm/v_norm);
-    #endif        // NORM
+        v_norm = sqrt(G*m_norm/r_norm);
+        t_norm = (r_norm/v_norm);
+#endif        // NORM
 
+        /* Read the data for EXTernal POTential */
+#ifdef EXTPOT
+        inp = fopen("phi-GRAPE.ext","r");
 
-    /* Read the data for EXTernal POTential */
-    #ifdef EXTPOT
+#ifdef EXTPOT_GAL
+        fscanf(inp,"%lE %lE %lE", &m_bulge, &a_bulge, &b_bulge);
+        fscanf(inp,"%lE %lE %lE", &m_disk, &a_disk, &b_disk);
+        fscanf(inp,"%lE %lE %lE", &m_halo, &a_halo, &b_halo);
 
-    inp = fopen("phi-GRAPE.ext","r");
+        m_bulge *= (Msol/m_norm);  a_bulge *= (kpc/r_norm);  b_bulge *= (kpc/r_norm);
+        m_disk  *= (Msol/m_norm);  a_disk  *= (kpc/r_norm);  b_disk  *= (kpc/r_norm);
+        m_halo  *= (Msol/m_norm);  a_halo  *= (kpc/r_norm);  b_halo  *= (kpc/r_norm);
+#endif
 
-    #ifdef EXTPOT_GAL
-    fscanf(inp,"%lE %lE %lE", &m_bulge, &a_bulge, &b_bulge);
-    fscanf(inp,"%lE %lE %lE", &m_disk, &a_disk, &b_disk);
-    fscanf(inp,"%lE %lE %lE", &m_halo, &a_halo, &b_halo);
+#ifdef EXTPOT_BH
+        fscanf(inp,"%lE %lE", &m_bh, &b_bh);
+#endif
 
-    m_bulge *= (Msol/m_norm);  a_bulge *= (kpc/r_norm);  b_bulge *= (kpc/r_norm);
-    m_disk  *= (Msol/m_norm);  a_disk  *= (kpc/r_norm);  b_disk  *= (kpc/r_norm);
-    m_halo  *= (Msol/m_norm);  a_halo  *= (kpc/r_norm);  b_halo  *= (kpc/r_norm);
-    #endif
+#ifdef EXTPOT_GAL_DEH
+        fscanf(inp,"%lE %lE %lE", &m_ext, &r_ext, &g_ext);
+#endif
 
-    #ifdef EXTPOT_BH
-    fscanf(inp,"%lE %lE", &m_bh, &b_bh);
-    #endif
+#ifdef EXTPOT_GAL_LOG
+        fscanf(inp,"%lE %lE", &v_halo, &r_halo);
+        v_halo *= (km/v_norm);  r_halo *= (pc/r_norm);
 
-    #ifdef EXTPOT_GAL_DEH
-    fscanf(inp,"%lE %lE %lE", &m_ext, &r_ext, &g_ext);
-    #endif
+        v2_halo = SQR(v_halo);
+        r2_halo = SQR(r_halo);
+#endif
 
-    #ifdef EXTPOT_GAL_LOG
-    fscanf(inp,"%lE %lE", &v_halo, &r_halo);
-    v_halo *= (km/v_norm);  r_halo *= (pc/r_norm);
-
-    v2_halo = SQR(v_halo);
-    r2_halo = SQR(r_halo);
-    #endif
-
-    fclose(inp);
-
-    #endif        // EXTPOT
+        fclose(inp);
+#endif        // EXTPOT
 
         /* read the global data for particles to the rootRank */
 
         read_data();
 
-    /* possible coordinate & velocity limits for ALL particles !!! */
-
-    #ifdef LIMITS_ALL_BEG
-
-        for (i=0; i<N; i++)
-        {
-
-        tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
-        tmp_rv = x[i][0]*v[i][0] + x[i][1]*v[i][1] + x[i][2]*v[i][2];
-
-        if ( (tmp_r > R_LIMITS_ALL) && (tmp_rv > 0.0) )
-            {
-
-            for (k=0;k<3;k++)
-            {
-            v[i][k] *= 1.0E-01;
-            } /* k */
+        /* possible coordinate & velocity limits for ALL particles !!! */
 
+#ifdef LIMITS_ALL_BEG
+        for (i=0; i<N; i++) {
+            tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
+            tmp_rv = x[i][0]*v[i][0] + x[i][1]*v[i][1] + x[i][2]*v[i][2];
+            if ( (tmp_r > R_LIMITS_ALL) && (tmp_rv > 0.0) ) {
+                for (k=0;k<3;k++) {
+                    v[i][k] *= 1.0E-01;
+                } /* k */
             } /* if */
-
         } /* i */
+#endif
 
-    #endif
-
-
-    #ifdef CMCORR
-
+#ifdef CMCORR
         /* possible CM correction of the initial datafile... */
-
         if ( (diskstep == 0) && (time_cur == 0.0) ) cm_corr();
-
-    #endif
-
-
+#endif
         printf("\n");
         printf("Begin the calculation of phi-GRAPE program on %03d processors\n", n_proc);
         printf("\n");
@@ -2466,212 +2358,189 @@ int main(int argc, char *argv[])
         printf("eta     = %.6E \n", eta);
         printf("\n");
 
-    #ifdef NORM
+#ifdef NORM
         printf("Normalization: \n");
         printf("\n");
         printf("m_norm = %.4E [Msol]  r_norm = %.4E [pc] \n", m_norm/Msol, r_norm/pc);
         printf("v_morm = %.4E [km/s]  t_morm = %.4E [Myr] \n", v_norm/km, t_norm/Myr);
         printf("\n");
-    #endif
+#endif
 
-    #ifdef EXTPOT
+#ifdef EXTPOT
         printf("External Potential: \n");
         printf("\n");
 
-    #ifdef EXTPOT_GAL
+#ifdef EXTPOT_GAL
         printf("m_bulge = %.4E [Msol]  a_bulge = %.4E  b_bulge = %.4E [kpc] \n", m_bulge*m_norm/Msol, a_bulge*r_norm/kpc, b_bulge*r_norm/kpc);
         printf("m_disk  = %.4E [Msol]  a_disk  = %.4E  b_disk  = %.4E [kpc] \n", m_disk*m_norm/Msol,  a_disk*r_norm/kpc,  b_disk*r_norm/kpc);
         printf("m_halo  = %.4E [Msol]  a_halo  = %.4E  b_halo  = %.4E [kpc] \n", m_halo*m_norm/Msol,  a_halo*r_norm/kpc,  b_halo*r_norm/kpc);
-    #endif
+#endif
 
-    #ifdef EXTPOT_BH
+#ifdef EXTPOT_BH
         printf("m_bh = %.4E  b_bh = %.4E \n", m_bh, b_bh);
-    #endif
+#endif
 
-    #ifdef EXTPOT_GAL_DEH
+#ifdef EXTPOT_GAL_DEH
         printf("m_ext  = %.6E r_ext  = %.6E \t g_ext = %.3E \n", m_ext, r_ext, g_ext);
-    #endif
+#endif
 
-    #ifdef EXTPOT_GAL_LOG
+#ifdef EXTPOT_GAL_LOG
         printf("v_halo = %.6E r_halo = %.6E \n", v_halo, r_halo);
-    #endif
+#endif
 
-    #ifdef EXTPOT_SC
+#ifdef EXTPOT_SC
         read_SC_mass();
 
-        printf("EXTPOT_SC: # of points %06d \t M_SC_TOT = %.6E \n", M_SC_DIM-1, m_sc[M_SC_DIM-1]);
-    #endif
+        printf("EXTPOT_SC:# of points %06d \t M_SC_TOT = %.6E \n", M_SC_DIM-1, m_sc[M_SC_DIM-1]);
+#endif
 
         printf("\n");
-    #endif
-
-
+#endif
 
         fflush(stdout);
 
+        if ( (diskstep == 0) && (time_cur == 0.0) ) {
+            out = fopen("contr.dat","w");
+            fclose(out);
 
-        if ( (diskstep == 0) && (time_cur == 0.0) )
-        {
+#ifdef TIMING
+            out = fopen("timing.dat","w");
+            fclose(out);
+#endif
 
-        out = fopen("contr.dat","w");
-        fclose(out);
+#ifdef ADD_BH2
 
-    #ifdef TIMING
-        out = fopen("timing.dat","w");
-        fclose(out);
-    #endif
+#ifdef BH_OUT
+            out = fopen("bh.dat","w");
+            fclose(out);
+#endif
 
-    #ifdef ADD_BH2
+#ifdef BH_OUT_NB
+            out = fopen("bh_nb.dat","w");
+            fclose(out);
+#endif
 
-    #ifdef BH_OUT
-        out = fopen("bh.dat","w");
-        fclose(out);
-    #endif
+#else
 
-    #ifdef BH_OUT_NB
-        out = fopen("bh_nb.dat","w");
-        fclose(out);
-    #endif
+#ifdef ADD_BH1
 
-    #else
+#ifdef BH_OUT
+            out = fopen("bh.dat","w");
+            fclose(out);
+#endif
 
-    #ifdef ADD_BH1
+#ifdef BH_OUT_NB
+            out = fopen("bh_nb.dat","w");
+            fclose(out);
+#endif
 
-    #ifdef BH_OUT
-        out = fopen("bh.dat","w");
-        fclose(out);
-    #endif
+#endif        // ADD_BH1
 
-    #ifdef BH_OUT_NB
-        out = fopen("bh_nb.dat","w");
-        fclose(out);
-    #endif
+#endif        // ADD_BH2
 
-    #endif        // ADD_BH1
+#ifdef STARDESTR
 
-    #endif        // ADD_BH2
+#ifdef ADD_BH2
+            out = fopen("mass-bh.dat","w");
+            m_bh1   = m[0];
+            m_bh2   = m[1];
+            num_bh1 = 0;
+            num_bh2 = 0;
+            fprintf(out,"%.8E \t %.8E  %06d \t %.8E  %06d \n",
+            time_cur, m_bh1, num_bh1, m_bh2, num_bh2);
+            fclose(out);
 
+            out = fopen("mass-bh.con","w");
+            m_bh1   = m[0];
+            m_bh2   = m[1];
+            num_bh1 = 0;
+            num_bh2 = 0;
+            fprintf(out,"%.8E \t %.8E  %06d \t %.8E  %06d \n",
+            time_cur, m_bh1, num_bh1, m_bh2, num_bh2);
+            fclose(out);
+#else
 
-    #ifdef STARDESTR
+#ifdef ADD_BH1
 
-    #ifdef ADD_BH2
+            out = fopen("mass-bh.dat","w");
+            m_bh1   = m[0];
+            num_bh1 = 0;
+            fprintf(out,"%.8E \t %.8E  %06d \n",
+            time_cur, m_bh1, num_bh1);
+            fclose(out);
 
-        out = fopen("mass-bh.dat","w");
-        m_bh1   = m[0];
-        m_bh2   = m[1];
-        num_bh1 = 0;
-        num_bh2 = 0;
-        fprintf(out,"%.8E \t %.8E  %06d \t %.8E  %06d \n",
-        time_cur, m_bh1, num_bh1, m_bh2, num_bh2);
-        fclose(out);
+            out = fopen("mass-bh.con","w");
+            m_bh1   = m[0];
+            num_bh1 = 0;
+            fprintf(out,"%.8E \t %.8E  %06d \n",
+            time_cur, m_bh1, num_bh1);
+            fclose(out);
 
-        out = fopen("mass-bh.con","w");
-        m_bh1   = m[0];
-        m_bh2   = m[1];
-        num_bh1 = 0;
-        num_bh2 = 0;
-        fprintf(out,"%.8E \t %.8E  %06d \t %.8E  %06d \n",
-        time_cur, m_bh1, num_bh1, m_bh2, num_bh2);
-        fclose(out);
+#endif        // ADD_BH1
 
-    #else
+#endif        // ADD_BH2
 
-    #ifdef ADD_BH1
+#endif        // STARDESTR
 
-        out = fopen("mass-bh.dat","w");
-        m_bh1   = m[0];
-        num_bh1 = 0;
-        fprintf(out,"%.8E \t %.8E  %06d \n",
-        time_cur, m_bh1, num_bh1);
-        fclose(out);
+#ifdef STARDESTR_EXT
+            num_bh = 0;
+            out = fopen("mass-bh.dat","w");
+            fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh);
+            fclose(out);
 
-        out = fopen("mass-bh.con","w");
-        m_bh1   = m[0];
-        num_bh1 = 0;
-        fprintf(out,"%.8E \t %.8E  %06d \n",
-        time_cur, m_bh1, num_bh1);
-        fclose(out);
+            num_bh = 0;
+            out = fopen("mass-bh.con","w");
+            fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh);
+            fclose(out);
 
-    #endif        // ADD_BH1
+#ifdef BH_OUT_NB_EXT
+            out = fopen("bh_nb.dat","w");
+            fclose(out);
+#endif
 
-    #endif        // ADD_BH2
+#endif        // STARDESTR_EXT
 
-    #endif        // STARDESTR
+#ifdef BBH_INF
+            out = fopen("bbh.inf","w");
+            fclose(out);
+#endif
+        } else {   // if (diskstep == 0)
 
+#ifdef STARDESTR
 
-    #ifdef STARDESTR_EXT
+#ifdef ADD_BH2
+            inp = fopen("mass-bh.con","r");
+            fscanf(inp,"%lE %lE %d %lE %d", &tmp, &m_bh1, &num_bh1, &m_bh2, &num_bh2);
+            fclose(inp);
+            printf("%.8E \t %.8E  %06d \t %.8E  %06d \n", tmp, m_bh1, num_bh1, m_bh2, num_bh2);
+            fflush(stdout);
+#else
 
-        num_bh = 0;
-        out = fopen("mass-bh.dat","w");
-        fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh);
-        fclose(out);
+#ifdef ADD_BH1
+            inp = fopen("mass-bh.con","r");
+            fscanf(inp,"%lE %lE %d", &tmp, &m_bh1, &num_bh1);
+            fclose(inp);
+            printf("%.8E \t %.8E \t %06d \n", tmp, m_bh1, num_bh1);
+            fflush(stdout);
+#endif        // ADD_BH1
 
-        num_bh = 0;
-        out = fopen("mass-bh.con","w");
-        fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh);
-        fclose(out);
+#endif        // ADD_BH2
 
-    #ifdef BH_OUT_NB_EXT
-        out = fopen("bh_nb.dat","w");
-        fclose(out);
-    #endif
-
-    #endif        // STARDESTR_EXT
-
-    #ifdef BBH_INF
-    out = fopen("bbh.inf","w");
-    fclose(out);
-    #endif
-
-
-        }
-        else    // if (diskstep == 0)
-        {
-
-    #ifdef STARDESTR
-
-    #ifdef ADD_BH2
-
-        inp = fopen("mass-bh.con","r");
-        fscanf(inp,"%lE %lE %d %lE %d", &tmp, &m_bh1, &num_bh1, &m_bh2, &num_bh2);
-        fclose(inp);
-        printf("%.8E \t %.8E  %06d \t %.8E  %06d \n", tmp, m_bh1, num_bh1, m_bh2, num_bh2);
-        fflush(stdout);
-
-    #else
-
-    #ifdef ADD_BH1
-
-        inp = fopen("mass-bh.con","r");
-        fscanf(inp,"%lE %lE %d", &tmp, &m_bh1, &num_bh1);
-        fclose(inp);
-        printf("%.8E \t %.8E \t %06d \n", tmp, m_bh1, num_bh1);
-        fflush(stdout);
-
-    #endif        // ADD_BH1
-
-    #endif        // ADD_BH2
-
-    #endif        // STARDESTR
-
-
-    #ifdef STARDESTR_EXT
-
-        inp = fopen("mass-bh.con","r");
-        fscanf(inp,"%lE %lE %d", &tmp, &m_bh, &num_bh);
-        fclose(inp);
-        printf("%.8E \t %.8E \t %06d \n", tmp, m_bh, num_bh);
-        fflush(stdout);
-
-    #endif        // STARDESTR_EXT
+#endif        // STARDESTR
 
+#ifdef STARDESTR_EXT
+            inp = fopen("mass-bh.con","r");
+            fscanf(inp,"%lE %lE %d", &tmp, &m_bh, &num_bh);
+            fclose(inp);
+            printf("%.8E \t %.8E \t %06d \n", tmp, m_bh, num_bh);
+            fflush(stdout);
+#endif        // STARDESTR_EXT
 
         }        // if (diskstep == 0)
 
-
         get_CPU_time(&CPU_time_real0, &CPU_time_user0, &CPU_time_syst0);
 
-        } /* if (myRank == rootRank) */
+    } /* if (myRank == rootRank) */
 
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
@@ -2686,16 +2555,16 @@ int main(int argc, char *argv[])
     MPI_Bcast(&dt_bh,    1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Bcast(&time_cur, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
 
-    #ifdef NORM
+#ifdef NORM
     MPI_Bcast(&m_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Bcast(&r_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Bcast(&v_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Bcast(&t_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-    #endif        // NORM
+#endif        // NORM
 
-    #ifdef EXTPOT
+#ifdef EXTPOT
 
-    #ifdef EXTPOT_GAL
+#ifdef EXTPOT_GAL
     MPI_Bcast(&m_bulge, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Bcast(&a_bulge, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Bcast(&b_bulge, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
@@ -2708,49 +2577,44 @@ int main(int argc, char *argv[])
     MPI_Bcast(&a_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Bcast(&b_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
 
-    #endif
+#endif
 
-    #ifdef EXTPOT_BH
+#ifdef EXTPOT_BH
     MPI_Bcast(&m_bh, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Bcast(&b_bh, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-    #endif
+#endif
 
-    #ifdef EXTPOT_GAL_DEH
+#ifdef EXTPOT_GAL_DEH
     MPI_Bcast(&m_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Bcast(&r_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Bcast(&g_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-    #endif
+#endif
 
-    #ifdef EXTPOT_GAL_LOG
+#ifdef EXTPOT_GAL_LOG
     MPI_Bcast(&v_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Bcast(&r_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
 
     MPI_Bcast(&v2_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Bcast(&r2_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-    #endif
+#endif
 
-    #endif        // EXTPOT
+#endif        // EXTPOT
 
-
-    #ifdef STARDESTR
+#ifdef STARDESTR
     MPI_Bcast(&m_bh1,   1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Bcast(&num_bh1, 1, MPI_INT,    rootRank, MPI_COMM_WORLD);
 
     MPI_Bcast(&m_bh2,   1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Bcast(&num_bh2, 1, MPI_INT,    rootRank, MPI_COMM_WORLD);
-    #endif        // STARDESTR
+#endif        // STARDESTR
 
-
-    #ifdef STARDESTR_EXT
+#ifdef STARDESTR_EXT
     MPI_Bcast(&num_bh, 1, MPI_INT,    rootRank, MPI_COMM_WORLD);
-    #endif        // STARDESTR_EXT
-
+#endif        // STARDESTR_EXT
 
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
 
-
-
     eta_s  = eta/ETA_S_CORR;
     eta_bh = eta/ETA_BH_CORR;
 
@@ -2770,24 +2634,19 @@ int main(int argc, char *argv[])
     t_contr = dt_contr*(1.0+floor(time_cur/dt_contr));
     t_bh    = dt_bh*(1.0+floor(time_cur/dt_bh));
 
-
-    if (myRank == rootRank)
-    {
-    printf("t_disk = %.6E   t_contr = %.6E   t_bh = %.6E \n", t_disk, t_contr, t_bh);
-    printf("\n");
-    fflush(stdout);
+    if (myRank == rootRank) {
+        printf("t_disk = %.6E   t_contr = %.6E   t_bh = %.6E \n", t_disk, t_contr, t_bh);
+        printf("\n");
+        fflush(stdout);
     } /* if (myRank == rootRank) */
 
-    for (i=0; i<N; i++)
-        {
+    for (i=0; i<N; i++) {
         t[i] = time_cur;
         dt[i] = dt_min;
-        }
-
+    }
 
     n_loc = N/n_proc;
 
-
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
 
@@ -2811,20 +2670,18 @@ int main(int argc, char *argv[])
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
 
-
-
     /* init the local GRAPE's */
 
-    #ifdef SAPPORO
+#ifdef SAPPORO
 
-    #ifdef LAOHUGPUDEVMGR
+#ifdef LAOHUGPUDEVMGR
     setRank(&myRank);
-    #endif
+#endif
     g6_open_(&clusterid);
     npipe = g6_npipes_();
     g6_set_tunit_(&new_tunit);
     g6_set_xunit_(&new_xunit);
-    #else
+#else
 
 #ifdef MPI_OVERRIDE
     numGPU = 1;
@@ -2843,27 +2700,25 @@ int main(int argc, char *argv[])
     g6_set_tunit(new_tunit);
     g6_set_xunit(new_xunit);
 
-    #ifndef YEBISU
+#ifndef YEBISU
     g6_set_overflow_flag_test_mode(aflag, jflag, pflag);
-    #endif
-    #endif
+#endif
+#endif
 
-    #ifdef AMD
+#ifdef AMD
     g6_initialize_jp_buffer(clusterid, 10000);
-    #endif
+#endif
 
-    #ifdef GUINNESS
+#ifdef GUINNESS
     g6_initialize_jp_buffer(clusterid, n_loc);
-    #endif
+#endif
 
-
-    #ifdef ETICS
+#ifdef ETICS
     grapite_read_particle_tags(N, "grapite.mask", myRank, n_loc);
     grapite_set_dt_exp(dt_exp);
     dt_exp = time_cur; // if we don't have a binary restart then we don't remember the coefficients, and we need to calculate them now.
     grapite_set_t_exp(t_exp);
-    #endif
-
+#endif
 
     /* initial load the particles to the local GRAPE's */
 
@@ -2871,34 +2726,31 @@ int main(int argc, char *argv[])
 
     /* load the nj particles to the G6 */
 
-    for (j=0; j<nj; j++)
-        {
-
-        for (k=0;k<3;k++)
-        {
-        a2by18[k] = 0.0; a1by6[k] = 0.0; aby2[k] = 0.0;
+    for (j=0; j<nj; j++) {
+        for (k=0;k<3;k++) {
+            a2by18[k] = 0.0; a1by6[k] = 0.0; aby2[k] = 0.0;
         } /* k */
 
-    #ifdef SAPPORO
+#ifdef SAPPORO
         g6_set_j_particle_(&clusterid, &j, &ind_loc[j], &t_loc[j], &dt_loc[j], &m_loc[j], a2by18, a1by6, aby2, v_loc[j], x_loc[j]);
-    #else
+#else
         g6_set_j_particle(clusterid, j, ind_loc[j], t_loc[j], dt_loc[j], m_loc[j], a2by18, a1by6, aby2, v_loc[j], x_loc[j]);
-    #endif
+#endif
 
-        } /* j */
+    } /* j */
 
-    #ifdef AMD
-        g6_flush_jp_buffer(clusterid);
-    #endif
+#ifdef AMD
+    g6_flush_jp_buffer(clusterid);
+#endif
 
-    #ifdef ETICS
-        double etics_length_scale;
-        if (myRank == rootRank) etics_length_scale = grapite_get_length_scale(N, m, x, v); // We don't want all ranks to do it, because they need to write a file and might confuse each other
-        MPI_Bcast(&etics_length_scale, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-        grapite_set_length_scale(etics_length_scale);
-    #endif
+#ifdef ETICS
+    double etics_length_scale;
+    if (myRank == rootRank) etics_length_scale = grapite_get_length_scale(N, m, x, v); // We don't want all ranks to do it, because they need to write a file and might confuse each other
+    MPI_Bcast(&etics_length_scale, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
+    grapite_set_length_scale(etics_length_scale);
+#endif
 
-    #ifdef ETICS_CEP
+#ifdef ETICS_CEP
     // First time only: get the CEP index
     grapite_cep_index = grapite_get_cep_index();
 
@@ -2917,22 +2769,19 @@ int main(int argc, char *argv[])
 
     double zeros[3] = {0., 0., 0.}; // We haven't calculated the force yet!
     grapite_update_cep(time_cur, xdc, vdc, zeros, zeros);
-    #endif
+#endif
 
     /* define the all particles as a active on all the processors for the first time grav calc. */
 
     n_act = N;
 
-    for (i=0; i<n_act; i++)
-        {
-
+    for (i=0; i<n_act; i++) {
         ind_act[i]  = ind[i];
         iii         = ind_act[i];
         m_act[i]    = m[iii];
 
-        for (k=0;k<3;k++)
-        {
-        x_act[i][k] = x[iii][k]; v_act[i][k] = v[iii][k];
+        for (k=0;k<3;k++) {
+            x_act[i][k] = x[iii][k]; v_act[i][k] = v[iii][k];
         } /* k */
 
         t_act[i]  = t[iii];
@@ -2940,10 +2789,8 @@ int main(int argc, char *argv[])
 
         } /* i */
 
-
     calc_self_grav_zero();
 
-
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
 
@@ -2955,26 +2802,23 @@ int main(int argc, char *argv[])
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
 
-
-
-    #ifdef ADD_BH2
+#ifdef ADD_BH2
 
     i_bh1 = 0;
     i_bh2 = 1;
 
-    #ifdef ADD_N_BH
+#ifdef ADD_N_BH
 
     m_bh1 = m[i_bh1];
     m_bh2 = m[i_bh2];
 
-    for (k=0;k<3;k++)
-        {
+    for (k=0;k<3;k++) {
         x_bh1[k] = x[i_bh1][k];
         v_bh1[k] = v[i_bh1][k];
 
         x_bh2[k] = x[i_bh2][k];
         v_bh2[k] = v[i_bh2][k];
-        }
+    }
 
     // calculate and "minus" the BH <-> BH _softened_ pot, acc & jerk
 
@@ -2987,14 +2831,13 @@ int main(int argc, char *argv[])
     pot[i_bh1] -= pot_bh1;
     pot[i_bh2] -= pot_bh2;
 
-    for (k=0;k<3;k++)
-        {
+    for (k=0;k<3;k++) {
         a[i_bh1][k] -= a_bh1[k];
         a[i_bh2][k] -= a_bh2[k];
 
         adot[i_bh1][k] -= adot_bh1[k];
         adot[i_bh2][k] -= adot_bh2[k];
-        }
+    }
 
     // calculate and "plus" the new BH <-> BH _unsoftened_ pot, acc, jerk
 
@@ -3007,58 +2850,51 @@ int main(int argc, char *argv[])
     pot[i_bh1] += pot_bh1;
     pot[i_bh2] += pot_bh2;
 
-    for (k=0;k<3;k++)
-        {
+    for (k=0;k<3;k++) {
         a[i_bh1][k] += a_bh1[k];
         a[i_bh2][k] += a_bh2[k];
 
         adot[i_bh1][k] += adot_bh1[k];
         adot[i_bh2][k] += adot_bh2[k];
-        }
+    }
 
-    #endif        // ADD_N_BH
+#endif        // ADD_N_BH
 
-
-    #ifdef ADD_PN_BH
+#ifdef ADD_PN_BH
 
     // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk
 
     dt_bh_tmp = dt[0];
 
     tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1,
-                            m_bh2, x_bh2, v_bh2, s_bh2,
-                            C_NB,  dt_bh_tmp, usedOrNot,
-                            a_pn1, adot_pn1,
-                            a_pn2, adot_pn2);
+                             m_bh2, x_bh2, v_bh2, s_bh2,
+                             C_NB,  dt_bh_tmp, usedOrNot,
+                             a_pn1, adot_pn1,
+                             a_pn2, adot_pn2);
 
-    for (k=0;k<3;k++)
-        {
+    for (k=0;k<3;k++) {
         a[i_bh1][k] += a_pn1[1][k] + a_pn1[2][k] + a_pn1[3][k] + a_pn1[4][k] + a_pn1[5][k] + a_pn1[6][k];
         a[i_bh2][k] += a_pn2[1][k] + a_pn2[2][k] + a_pn2[3][k] + a_pn2[4][k] + a_pn2[5][k] + a_pn2[6][k];
 
         adot[i_bh1][k] += adot_pn1[1][k] + adot_pn1[2][k] + adot_pn1[3][k] + adot_pn1[4][k] + adot_pn1[5][k] + adot_pn1[6][k];
         adot[i_bh2][k] += adot_pn2[1][k] + adot_pn2[2][k] + adot_pn2[3][k] + adot_pn2[4][k] + adot_pn2[5][k] + adot_pn2[6][k];
+    }
+
+    if (myRank == rootRank) {
+        if (tmp_i == 505) {
+            printf("PN RSDIST: %.8E \t %.8E \n", Timesteps, time_cur);
+            fflush(stdout);
+            exit(-1);
         }
+    }
 
-    if (myRank == rootRank)
-        {
-        if (tmp_i == 505)
-        {
-        printf("PN RSDIST: %.8E \t %.8E \n", Timesteps, time_cur);
-        fflush(stdout);
-        exit(-1);
-        }
-        }
+#endif        // ADD_PN_BH
 
-    #endif        // ADD_PN_BH
+#endif        // ADD_BH2
 
-    #endif        // ADD_BH2
-
-
-
-    #ifdef EXTPOT
+#ifdef EXTPOT
     calc_ext_grav_zero();
-    #endif
+#endif
 
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
@@ -3066,22 +2902,17 @@ int main(int argc, char *argv[])
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
 
-
     /* Energy control... */
-    if (myRank == rootRank)
-        {
-
+    if (myRank == rootRank) {
         energy_contr();
+    } /* if (myRank == rootRank) */
 
-        } /* if (myRank == rootRank) */
-
-    #ifdef ETICS_DUMP
-        if (diskstep==0) {
+#ifdef ETICS_DUMP
+    if (diskstep==0) {
         sprintf(out_fname, "coeffs.%06d.%02d.dat", 0, myRank);
         grapite_dump(out_fname, 2);
-        }
-    #endif
-
+    }
+#endif
 
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
@@ -3091,7 +2922,7 @@ int main(int argc, char *argv[])
     MPI_Scatter(a,    3*n_loc, MPI_DOUBLE, a_loc,    3*n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Scatter(adot, 3*n_loc, MPI_DOUBLE, adot_loc, 3*n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
 
-    #ifdef ETICS_CEP
+#ifdef ETICS_CEP
     // First calculate the DC
     grapite_calc_center(N, m, x, v, xcm, vcm, xdc, vdc);
 
@@ -3106,26 +2937,21 @@ int main(int argc, char *argv[])
     }
 
     grapite_update_cep(time_cur, xdc, vdc, a[grapite_cep_index], adot[grapite_cep_index]);
-    #endif
+#endif
 
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
 
-
-
-
-
     /* Define initial timestep for all particles on all nodes */
 
-    for (i=0; i<N; i++)
-        {
+    for (i=0; i<N; i++) {
         a2_mod = a[i][0]*a[i][0]+a[i][1]*a[i][1]+a[i][2]*a[i][2];
         adot2_mod = adot[i][0]*adot[i][0]+adot[i][1]*adot[i][1]+adot[i][2]*adot[i][2];
 
         if (adot2_mod == 0)
-        dt_tmp = eta_s;
+            dt_tmp = eta_s;
         else
-        dt_tmp = eta_s*sqrt(a2_mod/adot2_mod);
+            dt_tmp = eta_s*sqrt(a2_mod/adot2_mod);
 
         power = log(dt_tmp)/log(2.0) - 1;
 
@@ -3134,73 +2960,61 @@ int main(int argc, char *argv[])
         if (dt_tmp < dt_min) dt_tmp = dt_min;
         if (dt_tmp > dt_max) dt_tmp = dt_max;
 
-
-    #ifdef STARDESTR
+#ifdef STARDESTR
         if (m[i] == 0.0) dt_tmp = dt_max;
-    #endif
+#endif
 
-    #ifdef STARDESTR_EXT
+#ifdef STARDESTR_EXT
         if (m[i] == 0.0) dt_tmp = dt_max;
-    #endif
+#endif
 
         dt[i] = dt_tmp;
 
-
-    #ifdef DT_MIN_WARNING
-        if (myRank == 0)
-        {
-        if (dt[i] == dt_min)
-            {
-            printf("!!! Warning0: dt = dt_min = %.6E \t ind = %07d \n", dt[i], ind[i]);
-            fflush(stdout);
+#ifdef DT_MIN_WARNING
+        if (myRank == 0) {
+            if (dt[i] == dt_min) {
+                printf("!!! Warning0: dt = dt_min = %.6E \t ind = %07d \n", dt[i], ind[i]);
+                fflush(stdout);
             }
         }
-    #endif
+#endif
 
-        } /* i */
+    } /* i */
 
-
-
-
-    #ifdef ADD_BH2
+#ifdef ADD_BH2
 
     /* define the min. dt over all the part. and set it also for the BH... */
 
-        min_dt = dt[0];
+    min_dt = dt[0];
 
-        for (i=1; i<N; i++)
-        {
+    for (i=1; i<N; i++) {
         if ( dt[i] < min_dt ) min_dt = dt[i];
-        } /* i */
+    } /* i */
 
-        dt[0] = min_dt;
-        dt[1] = min_dt;
+    dt[0] = min_dt;
+    dt[1] = min_dt;
 
-    #else
+#else
 
-    #ifdef ADD_BH1
+#ifdef ADD_BH1
 
     /* define the min. dt over all the part. and set it also for the BH... */
 
-        min_dt = dt[0];
+    min_dt = dt[0];
 
-        for (i=1; i<N; i++)
-        {
+    for (i=1; i<N; i++) {
         if ( dt[i] < min_dt ) min_dt = dt[i];
-        } /* i */
+    } /* i */
 
-        dt[0] = min_dt;
+    dt[0] = min_dt;
 
-    #endif        // ADD_BH1
+#endif        // ADD_BH1
 
-    #endif        // ADD_BH2
+#endif        // ADD_BH2
 
-
-    #ifdef ACT_DEF_LL
+#ifdef ACT_DEF_LL
     CreateLinkList();
-    #endif
-
-
+#endif
 
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
@@ -3217,69 +3031,58 @@ int main(int argc, char *argv[])
 
     /* load the nj particles to the G6 */
 
-    for (j=0; j<nj; j++)
-        {
-        for (k=0;k<3;k++)
-        {
-        a2by18[k] = 0.0; a1by6[k] = adot_loc[j][k]*over6; aby2[k] = a_loc[j][k]*over2;
+    for (j=0; j<nj; j++) {
+        for (k=0;k<3;k++) {
+            a2by18[k] = 0.0; a1by6[k] = adot_loc[j][k]*over6; aby2[k] = a_loc[j][k]*over2;
         } /* k */
-
-    #ifdef SAPPORO
+#ifdef SAPPORO
         g6_set_j_particle_(&clusterid, &j, &ind_loc[j], &t_loc[j], &dt_loc[j], &m_loc[j], a2by18, a1by6, aby2, v_loc[j], x_loc[j]);
-    #else
+#else
         g6_set_j_particle(clusterid, j, ind_loc[j], t_loc[j], dt_loc[j], m_loc[j], a2by18, a1by6, aby2, v_loc[j], x_loc[j]);
-    #endif
+#endif
+    } /* j */
 
-        } /* j */
-
-    #ifdef AMD
+#ifdef AMD
         g6_flush_jp_buffer(clusterid);
-    #endif
+#endif
 
-    if (myRank == rootRank)
-        {
+    if (myRank == rootRank) {
 
-    #ifdef BH_OUT
+#ifdef BH_OUT
         /* Write BH data... */
         write_bh_data();
-    #endif
+#endif
 
-    #ifdef BH_OUT_NB
+#ifdef BH_OUT_NB
         /* Write BH NB data... */
         write_bh_nb_data();
-    #endif
+#endif
 
-    #ifdef BH_OUT_NB_EXT
+#ifdef BH_OUT_NB_EXT
         /* Write BH NB data... */
         write_bh_nb_data_ext();
-    #endif
+#endif
 
-        } /* if (myRank == rootRank) */
+    } /* if (myRank == rootRank) */
 
     /* Get the Starting time on rootRank */
 
-    if (myRank == rootRank)
-        {
+    if (myRank == rootRank) {
         get_CPU_time(&CPU_time_real0, &CPU_time_user0, &CPU_time_syst0);
         get_CPU_time(&CPU_time_real,  &CPU_time_user,  &CPU_time_syst);
         tmp_cpu = CPU_time_real-CPU_time_real0;
-        } /* if (myRank == rootRank) */
-
-
+    } /* if (myRank == rootRank) */
 
     Timesteps = 0.0;
     n_act_sum = 0.0;
 
-    for (i=1; i<N+1; i++)
-        {
+    for (i=1; i<N+1; i++) {
         n_act_distr[i-1] = 0.0;
-        }
+    }
 
     g6_calls = 0.0;
 
-
-
-    #ifdef TIMING
+#ifdef TIMING
     DT_TOT      = 0.0;
 
     DT_ACT_DEF1 = 0.0;
@@ -3298,1170 +3101,983 @@ int main(int argc, char *argv[])
     DT_STARDESTR = 0.0;
 
     DT_ACT_REDUCE = 0.0;
-    #endif
-
+#endif
 
     /* The main integration loop */
 
-    while (time_cur <= t_end)
-        {
+    while (time_cur <= t_end) {
 
-    /* Define the minimal time and the active particles on all the nodes (exclude the ZERO masses!!!) */
+        /* Define the minimal time and the active particles on all the nodes (exclude the ZERO masses!!!) */
 
+#ifdef ACT_DEF_LL
 
-    #ifdef ACT_DEF_LL
-
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
-    #endif
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
+#endif
 
         get_act_plist();
 
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
-    DT_ACT_DEF1 += (CPU_tmp_user - CPU_tmp_user0);
-    #endif
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
+        DT_ACT_DEF1 += (CPU_tmp_user - CPU_tmp_user0);
+#endif
 
+#else
 
-    #else
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
+#endif
 
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
-    #endif
-
-    #ifdef ACT_DEF_GRAPITE
+#ifdef ACT_DEF_GRAPITE
         min_t_loc = grapite_get_minimum_time();
-    #else
-    min_t_loc = t[0]+dt[0];
+#else
+        min_t_loc = t[0]+dt[0];
 
-    for (j=0; j<n_loc; j++)
-        {
-        jjj = j + myRank*n_loc;
-        tmp = t[jjj] + dt[jjj];
-        if ( tmp < min_t_loc ) min_t_loc = tmp;
+        for (j=0; j<n_loc; j++) {
+            jjj = j + myRank*n_loc;
+            tmp = t[jjj] + dt[jjj];
+            if ( tmp < min_t_loc ) min_t_loc = tmp;
         }
-    #endif
+#endif
 
-    /* Wait to all processors to finish his works... */
-    MPI_Barrier(MPI_COMM_WORLD);
+        /* Wait to all processors to finish his works... */
+        MPI_Barrier(MPI_COMM_WORLD);
 
-    /* Reduce the "global" min_t from min_t_loc "local" on all processors) */
-    MPI_Allreduce(&min_t_loc, &min_t, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+        /* Reduce the "global" min_t from min_t_loc "local" on all processors) */
+        MPI_Allreduce(&min_t_loc, &min_t, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
 
-    /* Wait to all processors to finish his works... */
-    MPI_Barrier(MPI_COMM_WORLD);
+        /* Wait to all processors to finish his works... */
+        MPI_Barrier(MPI_COMM_WORLD);
 
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
-    DT_ACT_DEF1 += (CPU_tmp_user - CPU_tmp_user0);
-    #endif
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
+        DT_ACT_DEF1 += (CPU_tmp_user - CPU_tmp_user0);
+#endif
 
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
+#endif
 
-
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
-    #endif
-
-    #ifdef ACT_DEF_GRAPITE
-    int ind_act_loc[N_MAX], n_act_loc;
-    grapite_active_search(min_t, ind_act_loc, &n_act_loc);
-    if (myRank > 0)
-        for (int i=0; i<n_act_loc; i++)
-        ind_act_loc[i] += myRank*n_loc;
-    int n_act_arr[256], displs[256]; // Assuming maximum of 256 processes... seems safe.
-    MPI_Allgather(&n_act_loc, 1, MPI_INT, n_act_arr, 1, MPI_INT, MPI_COMM_WORLD);
-    n_act = n_act_arr[0];
-    for (int i=1; i<n_proc; i++)
-        n_act += n_act_arr[i];
-    displs[0] = 0;
-    for (int i=1; i<n_proc; i++)
-        displs[i]=displs[i-1]+n_act_arr[i-1];
-    MPI_Allgatherv(ind_act_loc, n_act_loc, MPI_INT, ind_act, n_act_arr, displs, MPI_INT, MPI_COMM_WORLD);
-    #else
+#ifdef ACT_DEF_GRAPITE
+        int ind_act_loc[N_MAX], n_act_loc;
+        grapite_active_search(min_t, ind_act_loc, &n_act_loc);
+        if (myRank > 0)
+            for (int i=0; i<n_act_loc; i++)
+                ind_act_loc[i] += myRank*n_loc;
+        int n_act_arr[256], displs[256]; // Assuming maximum of 256 processes... seems safe.
+        MPI_Allgather(&n_act_loc, 1, MPI_INT, n_act_arr, 1, MPI_INT, MPI_COMM_WORLD);
+        n_act = n_act_arr[0];
+        for (int i=1; i<n_proc; i++)
+            n_act += n_act_arr[i];
+        displs[0] = 0;
+        for (int i=1; i<n_proc; i++)
+            displs[i]=displs[i-1]+n_act_arr[i-1];
+        MPI_Allgatherv(ind_act_loc, n_act_loc, MPI_INT, ind_act, n_act_arr, displs, MPI_INT, MPI_COMM_WORLD);
+#else
         n_act = 0;
 
-        for (i=0; i<N; i++)
-        {
-        if ( (t[i]+dt[i]) == min_t )
-            {
-            ind_act[n_act] = i;
-            n_act++;
+        for (i=0; i<N; i++) {
+            if ( (t[i]+dt[i]) == min_t ) {
+                ind_act[n_act] = i;
+                n_act++;
             }
         } /* i */
-    #endif      // ACT_DEF_GRAPITE
+#endif      // ACT_DEF_GRAPITE
 
-    #if defined(ACT_DEF_GRAPITE) && (defined(ADD_BH1) || defined(ADD_BH2))
-    #ifdef ADD_BH1
-    #define ACT_DEF_GRAPITE_NUMBH 1
-    #else
-    #define ACT_DEF_GRAPITE_NUMBH 2
-    #endif
+#if defined(ACT_DEF_GRAPITE) && (defined(ADD_BH1) || defined(ADD_BH2))
+#ifdef ADD_BH1
+#define ACT_DEF_GRAPITE_NUMBH 1
+#else
+#define ACT_DEF_GRAPITE_NUMBH 2
+#endif
         int acf_def_grapite_bh_count = 0;
         for (i=0; i<n_act; i++) {
-        if (ind_act[i]==0) {
-            i_bh1 = i;
-            acf_def_grapite_bh_count++;
+            if (ind_act[i]==0) {
+                i_bh1 = i;
+                acf_def_grapite_bh_count++;
         }
-    #ifdef ADD_BH2
-        else if (ind_act[i]==1) {
-            i_bh2 = i;
-            acf_def_grapite_bh_count++;
-        }
-    #endif
-        if (acf_def_grapite_bh_count==ACT_DEF_GRAPITE_NUMBH) break;
+#ifdef ADD_BH2
+            else if (ind_act[i]==1) {
+                i_bh2 = i;
+                acf_def_grapite_bh_count++;
+            }
+#endif
+            if (acf_def_grapite_bh_count==ACT_DEF_GRAPITE_NUMBH) break;
         }
         if (i==n_act) {
-        fprintf(stderr, "ERROR: black holes were not found in the active particle list");
-        exit(1);
+            fprintf(stderr, "ERROR: black holes were not found in the active particle list");
+            exit(1);
         }
-    #elif defined(ADD_BH1) || defined(ADD_BH2)
+#elif defined(ADD_BH1) || defined(ADD_BH2)
         i_bh1 = 0;
         i_bh2 = 1;
-    #endif
+#endif
 
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
-    DT_ACT_DEF2 += (CPU_tmp_user - CPU_tmp_user0);
-    #endif
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
+        DT_ACT_DEF2 += (CPU_tmp_user - CPU_tmp_user0);
+#endif
 
-    #endif      // ACT_DEF_LL
+#endif      // ACT_DEF_LL
 
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
+#endif
+        for (i=0; i<n_act; i++) {
+            iii = ind_act[i];
 
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
-    #endif
+            m_act[i] = m[iii];
 
-        for (i=0; i<n_act; i++)
-        {
-        iii = ind_act[i];
-
-        m_act[i] = m[iii];
-
-        for (k=0;k<3;k++)
-            {
-            x_act[i][k] = x[iii][k];
-            v_act[i][k] = v[iii][k];
+            for (k=0;k<3;k++) {
+                x_act[i][k] = x[iii][k];
+                v_act[i][k] = v[iii][k];
             } /* k */
 
-        t_act[i]  = t[iii];
-        dt_act[i] = dt[iii];
+            t_act[i]  = t[iii];
+            dt_act[i] = dt[iii];
 
-        pot_act[i] = pot[iii];
+            pot_act[i] = pot[iii];
 
-    #ifdef EXTPOT
-        pot_act_ext[i] = pot_ext[iii];
-    #endif
+#ifdef EXTPOT
+            pot_act_ext[i] = pot_ext[iii];
+#endif
 
-        for (k=0;k<3;k++)
-            {
-            a_act[i][k] = a[iii][k];
-            adot_act[i][k] = adot[iii][k];
+            for (k=0;k<3;k++) {
+                a_act[i][k] = a[iii][k];
+                adot_act[i][k] = adot[iii][k];
             } /* k */
 
         } /* i */
 
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
+        DT_ACT_DEF3 += (CPU_tmp_user - CPU_tmp_user0);
+#endif
 
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
-    DT_ACT_DEF3 += (CPU_tmp_user - CPU_tmp_user0);
-    #endif
+        /* predict the active particles positions etc... on all the nodes */
 
-    /* predict the active particles positions etc... on all the nodes */
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
+#endif
 
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
-    #endif
+        for (i=0; i<n_act; i++) {
+            dt_tmp = min_t - t_act[i];
+            dt2half = over2*dt_tmp*dt_tmp;
+            dt3over6 = over3*dt_tmp*dt2half;
 
-        for (i=0; i<n_act; i++)
-        {
-
-        dt_tmp = min_t - t_act[i];
-        dt2half = over2*dt_tmp*dt_tmp;
-        dt3over6 = over3*dt_tmp*dt2half;
-
-        for (k=0;k<3;k++)
-            {
-            x_act_new[i][k] = x_act[i][k] + v_act[i][k]*dt_tmp + a_act[i][k]*dt2half + adot_act[i][k]*dt3over6;
-            v_act_new[i][k] = v_act[i][k] + a_act[i][k]*dt_tmp + adot_act[i][k]*dt2half;
+            for (k=0;k<3;k++) {
+                x_act_new[i][k] = x_act[i][k] + v_act[i][k]*dt_tmp + a_act[i][k]*dt2half + adot_act[i][k]*dt3over6;
+                v_act_new[i][k] = v_act[i][k] + a_act[i][k]*dt_tmp + adot_act[i][k]*dt2half;
             } /* k */
 
         } /* i */
 
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
-    DT_ACT_PRED += (CPU_tmp_user - CPU_tmp_user0);
-    #endif
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
+        DT_ACT_PRED += (CPU_tmp_user - CPU_tmp_user0);
+#endif
 
-    calc_self_grav();
+        calc_self_grav();
 
+        /* Reduce the "global" vectors from "local" on all the nodes */
 
-    /* Reduce the "global" vectors from "local" on all the nodes */
-
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
-    #endif
-
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
+#endif
 
         MPI_Allreduce(pot_act_tmp,  pot_act_new,    n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 
-    MPI_Allreduce(a_act_tmp,    a_act_new,    3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-    MPI_Allreduce(adot_act_tmp, adot_act_new, 3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+        MPI_Allreduce(a_act_tmp,    a_act_new,    3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+        MPI_Allreduce(adot_act_tmp, adot_act_new, 3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 
-    /* Wait to all processors to finish his works... */
-    MPI_Barrier(MPI_COMM_WORLD);
+        /* Wait to all processors to finish his works... */
+        MPI_Barrier(MPI_COMM_WORLD);
 
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
-    DT_ACT_REDUCE += (CPU_tmp_user - CPU_tmp_user0);
-    #endif
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
+        DT_ACT_REDUCE += (CPU_tmp_user - CPU_tmp_user0);
+#endif
 
+#ifdef ADD_BH2
 
-
-    #ifdef ADD_BH2
-
-    #ifdef ADD_N_BH
-
-    m_bh1 = m_act[i_bh1];
-    m_bh2 = m_act[i_bh2];
-
-    for (k=0;k<3;k++)
-        {
-        x_bh1[k] = x_act_new[i_bh1][k];
-        v_bh1[k] = v_act_new[i_bh1][k];
-
-        x_bh2[k] = x_act_new[i_bh2][k];
-        v_bh2[k] = v_act_new[i_bh2][k];
-        }
-
-    // calculate and "minus" the BH <-> BH softened pot, acc & jerk
-
-    tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1,
-                            m_bh2, x_bh2, v_bh2,
-                            eps,
-                            &pot_bh1, a_bh1, adot_bh1,
-                            &pot_bh2, a_bh2, adot_bh2);
-
-    pot_act_new[i_bh1] -= pot_bh1;
-    pot_act_new[i_bh2] -= pot_bh2;
-
-    for (k=0;k<3;k++)
-        {
-        a_act_new[i_bh1][k] -= a_bh1[k];
-        a_act_new[i_bh2][k] -= a_bh2[k];
-
-        adot_act_new[i_bh1][k] -= adot_bh1[k];
-        adot_act_new[i_bh2][k] -= adot_bh2[k];
-        }
-
-    // calculate and "plus" the new BH <-> BH unsoftened pot, acc, jerk
-
-    tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1,
-                            m_bh2, x_bh2, v_bh2,
-                            eps_BH,
-                            &pot_bh1, a_bh1, adot_bh1,
-                            &pot_bh2, a_bh2, adot_bh2);
-
-    pot_act_new[i_bh1] += pot_bh1;
-    pot_act_new[i_bh2] += pot_bh2;
-
-    for (k=0;k<3;k++)
-        {
-        a_act_new[i_bh1][k] += a_bh1[k];
-        a_act_new[i_bh2][k] += a_bh2[k];
-
-        adot_act_new[i_bh1][k] += adot_bh1[k];
-        adot_act_new[i_bh2][k] += adot_bh2[k];
-        }
-
-    #endif        // ADD_N_BH
-
-
-    #ifdef ADD_PN_BH
-
-    // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk
-
-    dt_bh_tmp = dt[0];
-
-    tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1,
-                            m_bh2, x_bh2, v_bh2, s_bh2,
-                            C_NB,  dt_bh_tmp, usedOrNot,
-                            a_pn1, adot_pn1,
-                            a_pn2, adot_pn2);
-
-    for (k=0;k<3;k++)
-        {
-        a_act_new[i_bh1][k] += a_pn1[1][k] + a_pn1[2][k] + a_pn1[3][k] + a_pn1[4][k] + a_pn1[5][k] + a_pn1[6][k];
-        a_act_new[i_bh2][k] += a_pn2[1][k] + a_pn2[2][k] + a_pn2[3][k] + a_pn2[4][k] + a_pn2[5][k] + a_pn2[6][k];
-
-        adot_act_new[i_bh1][k] += adot_pn1[1][k] + adot_pn1[2][k] + adot_pn1[3][k] + adot_pn1[4][k] + adot_pn1[5][k] + adot_pn1[6][k];
-        adot_act_new[i_bh2][k] += adot_pn2[1][k] + adot_pn2[2][k] + adot_pn2[3][k] + adot_pn2[4][k] + adot_pn2[5][k] + adot_pn2[6][k];
-        }
-
-    if (myRank == rootRank)
-        {
-        if (tmp_i == 505)
-        {
-        printf("PN RSDIST: TS = %.8E \t t  = %.8E \n", Timesteps, time_cur);
-        fflush(stdout);
-        exit(-1);
-        }
-        }
-
-    #endif        // ADD_PN_BH
-
-    #endif        // ADD_BH2
-
-    #ifdef EXTPOT
-    calc_ext_grav();
-    #endif
-
-    /* correct the active particles positions etc... on all the nodes */
-
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
-    #endif
-
-    for (i=0; i<n_act; i++)
-        {
-
-        dt_tmp = min_t - t_act[i];
-
-        dt3over6 = dt_tmp*dt_tmp*dt_tmp/6.0;
-        dt4over24 = dt3over6*dt_tmp/4.0;
-        dt5over120 = dt4over24*dt_tmp/5.0;
-
-        dtinv = 1.0/dt_tmp;
-        dt2inv = dtinv*dtinv;
-        dt3inv = dt2inv*dtinv;
-
-        for (k=0; k<3; k++)
-        {
-        a0mia1 = a_act[i][k]-a_act_new[i][k];
-        ad04plad12 = 4.0*adot_act[i][k] + 2.0*adot_act_new[i][k];
-        ad0plad1 = adot_act[i][k] + adot_act_new[i][k];
-
-        a2[k] = -6.0*a0mia1*dt2inv - ad04plad12*dtinv;
-        a3[k] = 12.0*a0mia1*dt3inv + 6.0*ad0plad1*dt2inv;
-
-        x_act_new[i][k] += dt4over24*a2[k] + dt5over120*a3[k];
-        v_act_new[i][k] += dt3over6*a2[k] + dt4over24*a3[k];
-        } /* k */
-
-        a1abs = sqrt(a_act_new[i][0]*a_act_new[i][0]+a_act_new[i][1]*a_act_new[i][1]+a_act_new[i][2]*a_act_new[i][2]);
-        adot1abs = sqrt(adot_act_new[i][0]*adot_act_new[i][0]+adot_act_new[i][1]*adot_act_new[i][1]+adot_act_new[i][2]*adot_act_new[i][2]);
-
-        for (k=0; k<3; k++) a2dot1[k] = a2[k] + dt_tmp*a3[k];
-
-        a2dot1abs = sqrt(a2dot1[0]*a2dot1[0]+a2dot1[1]*a2dot1[1]+a2dot1[2]*a2dot1[2]);
-        a3dot1abs = sqrt(a3[0]*a3[0]+a3[1]*a3[1]+a3[2]*a3[2]);
-
-    #ifdef ADD_BH2
-
-        if ( (ind_act[i] == 0) || (ind_act[i] == 1) )
-        {
-        dt_new = sqrt(eta_bh*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
-        }
-        else
-        {
-        dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
-        }
-
-    #else
-
-    #ifdef ADD_BH1
-
-        if ( (ind_act[i] == 0) )
-        {
-        dt_new = sqrt(eta_bh*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
-        }
-        else
-        {
-        dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
-        }
-
-    #endif        // ADD_BH1
-
-        dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
-
-
-    #endif        // ADD_BH2
-
-
-    #ifdef STARDESTR
-        if (m_act[i] == 0.0) dt_new = dt_max;
-    #endif
-
-    #ifdef STARDESTR_EXT
-        if (m_act[i] == 0.0) dt_new = dt_max;
-    #endif
-
-        if (dt_new < dt_min) dt_tmp = dt_min;
-
-
-        if ( (dt_new < dt_tmp) && (dt_new > dt_min) )
-        {
-        power = log(dt_new)/log(2.0) - 1;
-
-        dt_tmp = pow(2.0, (double)power);
-        }
-
-        if ( (dt_new > 2.0*dt_tmp) &&
-            (fmod(min_t, 2.0*dt_tmp) == 0.0) &&
-            (2.0*dt_tmp <= dt_max) )
-        {
-        dt_tmp *= 2.0;
-        }
-
-        dt_act[i] = dt_tmp;
-
-        t_act[i] = min_t;
-
-        pot_act[i] = pot_act_new[i];
-
-        for (k=0; k<3; k++)
-        {
-            x_act[i][k] =    x_act_new[i][k];
-            v_act[i][k] =    v_act_new[i][k];
-            a_act[i][k] =    a_act_new[i][k];
-        adot_act[i][k] = adot_act_new[i][k];
-        } /* k */
-
-
-    #ifdef DT_MIN_WARNING
-        if (myRank == 0)
-        {
-        if (dt_act[i] == dt_min)
-            {
-            printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d \n", dt[i], ind_act[i]);
-            fflush(stdout);
-            }
-        }
-    #endif
-
-
-
-        } /* i */
-
-
-
-
-    /* define the min. dt over all the act. part. and set it also for the BH... */
-
-    #ifdef ADD_BH2
-
-        min_dt = dt_act[0];
-
-        for (i=1; i<n_act; i++)
-        {
-        if ( dt_act[i] < min_dt ) min_dt = dt_act[i];
-        } /* i */
-
-        dt_act[i_bh1] = min_dt;
-        dt_act[i_bh2] = min_dt;
-
-    #else
-
-    #ifdef ADD_BH1
-
-        min_dt = dt_act[0];
-
-        for (i=1; i<n_act; i++)
-        {
-        if ( dt_act[i] < min_dt ) min_dt = dt_act[i];
-        } /* i */
-
-        dt_act[i_bh1] = min_dt;
-
-    #endif        // ADD_BH1
-
-    #endif        // ADD_BH2
-
-
-
-    #ifdef BBH_INF
-    if (myRank == rootRank)
-        {
-
-        out = fopen("bbh.inf","a");
+#ifdef ADD_N_BH
 
         m_bh1 = m_act[i_bh1];
         m_bh2 = m_act[i_bh2];
 
-        for (k=0;k<3;k++)
-        {
-        x_bh1[k] = x_act[i_bh1][k];
-        x_bh2[k] = x_act[i_bh2][k];
+        for (k=0;k<3;k++) {
+            x_bh1[k] = x_act_new[i_bh1][k];
+            v_bh1[k] = v_act_new[i_bh1][k];
 
-        v_bh1[k] = v_act[i_bh1][k];
-        v_bh2[k] = v_act[i_bh2][k];
+            x_bh2[k] = x_act_new[i_bh2][k];
+            v_bh2[k] = v_act_new[i_bh2][k];
         }
 
-        for (k=0;k<3;k++)
-        {
-        x_bbhc[k] = (m_bh1*x_bh1[k] + m_bh2*x_bh2[k])/(m_bh1 + m_bh2);
-        v_bbhc[k] = (m_bh1*v_bh1[k] + m_bh2*v_bh2[k])/(m_bh1 + m_bh2);
+        // calculate and "minus" the BH <-> BH softened pot, acc & jerk
+
+        tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1,
+                                m_bh2, x_bh2, v_bh2,
+                                eps,
+                                &pot_bh1, a_bh1, adot_bh1,
+                                &pot_bh2, a_bh2, adot_bh2);
+
+        pot_act_new[i_bh1] -= pot_bh1;
+        pot_act_new[i_bh2] -= pot_bh2;
+
+        for (k=0;k<3;k++) {
+            a_act_new[i_bh1][k] -= a_bh1[k];
+            a_act_new[i_bh2][k] -= a_bh2[k];
+
+            adot_act_new[i_bh1][k] -= adot_bh1[k];
+            adot_act_new[i_bh2][k] -= adot_bh2[k];
         }
 
-        DR2 = SQR(x_bh1[0] - x_bh2[0]) + SQR(x_bh1[1] - x_bh2[1]) + SQR(x_bh1[2] - x_bh2[2]);
-        DV2 = SQR(v_bh1[0] - v_bh2[0]) + SQR(v_bh1[1] - v_bh2[1]) + SQR(v_bh1[2] - v_bh2[2]);
+    // calculate and "plus" the new BH <-> BH unsoftened pot, acc, jerk
 
-        EB = -(m_bh1 + m_bh2) / sqrt(DR2) + 0.5 * DV2;
+        tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1,
+                                m_bh2, x_bh2, v_bh2,
+                                eps_BH,
+                                &pot_bh1, a_bh1, adot_bh1,
+                                &pot_bh2, a_bh2, adot_bh2);
 
-        SEMI_a = -0.5 * (m_bh1 + m_bh2)/EB;
-        SEMI_a2 = SQR(SEMI_a);
+        pot_act_new[i_bh1] += pot_bh1;
+        pot_act_new[i_bh2] += pot_bh2;
 
-        for (i=2; i<n_act; i++)
-        {
+        for (k=0;k<3;k++) {
+            a_act_new[i_bh1][k] += a_bh1[k];
+            a_act_new[i_bh2][k] += a_bh2[k];
 
-        tmp_r2 = SQR(x_act[i][0] - x_bbhc[0]) + SQR(x_act[i][1] - x_bbhc[1]) + SQR(x_act[i][2] - x_bbhc[2]);
+            adot_act_new[i_bh1][k] += adot_bh1[k];
+            adot_act_new[i_bh2][k] += adot_bh2[k];
+        }
 
-        iii = ind_act[i];
+#endif        // ADD_N_BH
 
-        if ( (tmp_r2 < SEMI_a2*R_INF2) )
-            {
+#ifdef ADD_PN_BH
 
-        if ( (inf_event[iii] == 0) )
-            {
+        // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk
 
-            fprintf(out,"INF1 %.6E %.16E %07d %07d   %.6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E \n",
-                    Timesteps, time_cur, i, ind_act[i],
-                    sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2],
-                    m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[0],
-                    m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[1],
-                    sqrt(tmp_r2),
-                    m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i],
-                    dt_act[i]);
+        dt_bh_tmp = dt[0];
 
-            inf_event[iii] = 1;
+        tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1,
+                                 m_bh2, x_bh2, v_bh2, s_bh2,
+                                 C_NB,  dt_bh_tmp, usedOrNot,
+                                 a_pn1, adot_pn1,
+                                 a_pn2, adot_pn2);
+
+        for (k=0;k<3;k++) {
+            a_act_new[i_bh1][k] += a_pn1[1][k] + a_pn1[2][k] + a_pn1[3][k] + a_pn1[4][k] + a_pn1[5][k] + a_pn1[6][k];
+            a_act_new[i_bh2][k] += a_pn2[1][k] + a_pn2[2][k] + a_pn2[3][k] + a_pn2[4][k] + a_pn2[5][k] + a_pn2[6][k];
+
+            adot_act_new[i_bh1][k] += adot_pn1[1][k] + adot_pn1[2][k] + adot_pn1[3][k] + adot_pn1[4][k] + adot_pn1[5][k] + adot_pn1[6][k];
+            adot_act_new[i_bh2][k] += adot_pn2[1][k] + adot_pn2[2][k] + adot_pn2[3][k] + adot_pn2[4][k] + adot_pn2[5][k] + adot_pn2[6][k];
+        }
+
+        if (myRank == rootRank) {
+            if (tmp_i == 505) {
+                printf("PN RSDIST: TS = %.8E \t t  = %.8E \n", Timesteps, time_cur);
+                fflush(stdout);
+                exit(-1);
             }
-
-            }
-        else
-            {
-
-        if ( (inf_event[iii] == 1) )
-            {
-
-            fprintf(out,"INF2 %.6E %.16E %07d %07d   %.6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E \n",
-                    Timesteps, time_cur, i, ind_act[i],
-                    sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2],
-                    m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[0],
-                    m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[1],
-                    sqrt(tmp_r2),
-                    m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i],
-                    dt_act[i]);
-            }
-
-            inf_event[iii] = 0;
-
-            } /* if (tmp_r2 < DR2*R_INF2) */
-
-        } /* i */
-
-        fclose(out);
-
-        } /* if (myRank == rootRank) */
-    #endif
-
-
-
-
-    /* Return back the new coordinates + etc... of active part. to the global data... */
-
-    for (i=0; i<n_act; i++)
-        {
-
-        iii = ind_act[i];
-
-        m[iii] = m_act[i];
-
-        for (k=0;k<3;k++)
-        {
-        x[iii][k] = x_act[i][k];
-        v[iii][k] = v_act[i][k];
-        } /* k */
-
-        t[iii] = t_act[i];
-        dt[iii] = dt_act[i];
-
-        pot[iii] = pot_act[i];
-
-    #ifdef EXTPOT
-        pot_ext[iii] = pot_act_ext[i];
-    #endif
-
-        for (k=0;k<3;k++)
-        {
-        a[iii][k] = a_act[i][k];
-        adot[iii][k] = adot_act[i][k];
-        } /* k */
-
-        } /* i */
-
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
-    DT_ACT_CORR += (CPU_tmp_user - CPU_tmp_user0);
-    #endif
-
-
-
-
-    /* STARDESTR routine on all the nodes */
-
-    #ifdef STARDESTR
-
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
-    #endif
-
-
-    #ifdef ADD_BH2
-
-    star_destr(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act,
-                        N, ind, m, x, v, pot, a, adot, t, &m_bh1, &num_bh1, i_bh1);
-
-    m_act[i_bh1] = m_bh1;
-
-    star_destr(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act,
-                        N, ind, m, x, v, pot, a, adot, t, &m_bh2, &num_bh2, i_bh2);
-
-    m_act[i_bh2] = m_bh2;
-
-    #else
-
-    #ifdef ADD_BH1
-
-    star_destr(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act,
-                        N, ind, m, x, v, pot, a, adot, t, &m_bh1, &num_bh1, i_bh1);
-
-    m_act[i_bh1] = m_bh1;
-
-
-    #endif        // ADD_BH1
-
-    #endif        // ADD_BH2
-
-
-    #ifdef LIMITS
-    for (i=0; i<n_act; i++)
-        {
-        if (m_act[i] == 0.0)
-        {
-        x_act[i][0] = R_LIMITS + 1.0*my_rand2();
-        x_act[i][1] = R_LIMITS + 1.0*my_rand2();
-        x_act[i][2] = R_LIMITS + 1.0*my_rand2();
-
-        v_act[i][0] = 1.0E-06*my_rand2();
-        v_act[i][1] = 1.0E-06*my_rand2();
-        v_act[i][2] = 1.0E-06*my_rand2();
         }
-        }
-    #endif        // LIMITS
+#endif        // ADD_PN_BH
 
+#endif        // ADD_BH2
 
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
-    DT_STARDESTR += (CPU_tmp_user - CPU_tmp_user0);
-    #endif
+#ifdef EXTPOT
+        calc_ext_grav();
+#endif
 
+        /* correct the active particles positions etc... on all the nodes */
 
-    #endif        // STARDESTR
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
+#endif
 
+        for (i=0; i<n_act; i++) {
+            dt_tmp = min_t - t_act[i];
 
+            dt3over6 = dt_tmp*dt_tmp*dt_tmp/6.0;
+            dt4over24 = dt3over6*dt_tmp/4.0;
+            dt5over120 = dt4over24*dt_tmp/5.0;
 
-    #ifdef STARDESTR_EXT
+            dtinv = 1.0/dt_tmp;
+            dt2inv = dtinv*dtinv;
+            dt3inv = dt2inv*dtinv;
 
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
-    #endif
+            for (k=0; k<3; k++) {
+                a0mia1 = a_act[i][k]-a_act_new[i][k];
+                ad04plad12 = 4.0*adot_act[i][k] + 2.0*adot_act_new[i][k];
+                ad0plad1 = adot_act[i][k] + adot_act_new[i][k];
 
-    star_destr_ext(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act,
-                        N, m, x, v, a, adot, t, &m_bh, &num_bh);
+                a2[k] = -6.0*a0mia1*dt2inv - ad04plad12*dtinv;
+                a3[k] = 12.0*a0mia1*dt3inv + 6.0*ad0plad1*dt2inv;
 
-    #ifdef LIMITS
-    for (i=0; i<n_act; i++)
-        {
-        if (m_act[i] == 0.0)
-        {
-        x_act[i][0] = R_LIMITS + 1.0*my_rand2();
-        x_act[i][1] = R_LIMITS + 1.0*my_rand2();
-        x_act[i][2] = R_LIMITS + 1.0*my_rand2();
-
-        v_act[i][0] = 1.0E-06*my_rand2();
-        v_act[i][1] = 1.0E-06*my_rand2();
-        v_act[i][2] = 1.0E-06*my_rand2();
-        }
-        } /* i */
-    #endif        // LIMITS
-
-
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
-    DT_STARDESTR += (CPU_tmp_user - CPU_tmp_user0);
-    #endif
-
-    #endif        // STARDESTR_EXT
-
-
-
-    #if defined(STARDESTR) || defined(STARDESTR_EXT)
-
-
-    /* Return back the new coordinates + etc... of active part. to the global data... */
-
-    for (i=0; i<n_act; i++)
-        {
-
-        iii = ind_act[i];
-
-        m[iii] = m_act[i];
-
-        for (k=0;k<3;k++)
-        {
-        x[iii][k] = x_act[i][k];
-        v[iii][k] = v_act[i][k];
-        } /* k */
-
-        t[iii] = t_act[i];
-        dt[iii] = dt_act[i];
-
-        pot[iii] = pot_act[i];
-
-    #ifdef EXTPOT
-        pot_ext[iii] = pot_act_ext[i];
-    #endif
-
-        for (k=0;k<3;k++)
-        {
-        a[iii][k] = a_act[i][k];
-        adot[iii][k] = adot_act[i][k];
-        } /* k */
-
-        } /* i */
-
-    #endif        // #if defined(STARDESTR) || defined(STARDESTR_EXT)
-
-
-    /* load the new values for active particles to the local GRAPE's */
-
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
-    #endif
-
-    for (j=0; j<n_act; j++)
-        {
-    #ifdef ETICS_CEP
-        if (ind_act[j] == grapite_cep_index) grapite_update_cep(t_act[j], x_act[j], v_act[j], a_act[j], adot_act[j]); // All ranks should do it.
-    #endif
-        cur_rank = ind_act[j]/n_loc;
-
-        if (myRank == cur_rank)
-        {
-
-        jjj = ind_act[j] - myRank*n_loc;
-
-        for (k=0;k<3;k++)
-            {
-            a2by18[k] = 0.0; a1by6[k] = adot_act[j][k]*over6; aby2[k] = a_act[j][k]*over2;
+                x_act_new[i][k] += dt4over24*a2[k] + dt5over120*a3[k];
+                v_act_new[i][k] += dt3over6*a2[k] + dt4over24*a3[k];
             } /* k */
 
-    #ifdef SAPPORO
-        g6_set_j_particle_(&clusterid, &jjj, &ind_act[j], &t_act[j], &dt_act[j], &m_act[j], a2by18, a1by6, aby2, v_act[j], x_act[j]);
-    #else
-        g6_set_j_particle(clusterid, jjj, ind_act[j], t_act[j], dt_act[j], m_act[j], a2by18, a1by6, aby2, v_act[j], x_act[j]);
-    #endif
+            a1abs = sqrt(a_act_new[i][0]*a_act_new[i][0]+a_act_new[i][1]*a_act_new[i][1]+a_act_new[i][2]*a_act_new[i][2]);
+            adot1abs = sqrt(adot_act_new[i][0]*adot_act_new[i][0]+adot_act_new[i][1]*adot_act_new[i][1]+adot_act_new[i][2]*adot_act_new[i][2]);
 
-        } /* if (myRank == cur_rank) */
+            for (k=0; k<3; k++) a2dot1[k] = a2[k] + dt_tmp*a3[k];
 
-        } /* j */
+            a2dot1abs = sqrt(a2dot1[0]*a2dot1[0]+a2dot1[1]*a2dot1[1]+a2dot1[2]*a2dot1[2]);
+            a3dot1abs = sqrt(a3[0]*a3[0]+a3[1]*a3[1]+a3[2]*a3[2]);
 
-    #ifdef AMD
-        g6_flush_jp_buffer(clusterid);
-    #endif
+#ifdef ADD_BH2
+            if ( (ind_act[i] == 0) || (ind_act[i] == 1) ) {
+                dt_new = sqrt(eta_bh*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
+            } else {
+                dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
+            }
+#else
 
-    #ifdef TIMING
-    get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
-    DT_ACT_LOAD += (CPU_tmp_user - CPU_tmp_user0);
-    #endif
+#ifdef ADD_BH1
+            if ( (ind_act[i] == 0) ) {
+                dt_new = sqrt(eta_bh*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
+            } else {
+                dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
+            }
+#endif        // ADD_BH1
 
+            dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
 
+#endif        // ADD_BH2
 
-    /* Current time set to min_t */
+#ifdef STARDESTR
+            if (m_act[i] == 0.0) dt_new = dt_max;
+#endif
 
-    time_cur = min_t;
+#ifdef STARDESTR_EXT
+            if (m_act[i] == 0.0) dt_new = dt_max;
+#endif
 
-    Timesteps += 1.0;
-    n_act_sum += n_act;
-    n_act_distr[n_act-1]++;
+            if (dt_new < dt_min) dt_tmp = dt_min;
 
-    if (time_cur >= t_bh)
-        {
+            if ( (dt_new < dt_tmp) && (dt_new > dt_min) ) {
+                power = log(dt_new)/log(2.0) - 1;
 
-        if (myRank == rootRank)
-        {
+                dt_tmp = pow(2.0, (double)power);
+            }
 
-    #ifdef BH_OUT
-        /* Write BH data... */
-        write_bh_data();
-    #endif
+            if ((dt_new > 2.0*dt_tmp) && (fmod(min_t, 2.0*dt_tmp) == 0.0) && (2.0*dt_tmp <= dt_max)) {
+                dt_tmp *= 2.0;
+            }
 
-    #ifdef BH_OUT_NB
-        /* Write BH NB data... */
-        write_bh_nb_data();
-    #endif
+            dt_act[i] = dt_tmp;
 
-    #ifdef BH_OUT_NB_EXT
-        /* Write BH NB data... */
-        write_bh_nb_data_ext();
-    #endif
+            t_act[i] = min_t;
 
-        } /* if (myRank == rootRank) */
+            pot_act[i] = pot_act_new[i];
 
-        t_bh += dt_bh;
-        } /* if (time_cur >= t_bh) */
-
-
-
-
-
-        if (time_cur >= t_contr)
-        {
-
-        if (myRank == rootRank)
-        {
-
-    #ifdef STARDESTR
-
-    #ifdef ADD_BH2
-
-        out = fopen("mass-bh.dat","a");
-        fprintf(out,"%.8E \t %.8E  %06d \t %.8E  %06d \n",
-        time_cur, m_bh1, num_bh1, m_bh2, num_bh2);
-        fclose(out);
-
-        out = fopen("mass-bh.con","w");
-        fprintf(out,"%.8E \t %.8E  %06d \t %.8E  %06d \n",
-        time_cur, m_bh1, num_bh1, m_bh2, num_bh2);
-        fclose(out);
-
-    #else
-
-    #ifdef ADD_BH1
-
-        out = fopen("mass-bh.dat","a");
-        fprintf(out,"%.8E \t %.8E  %06d \n",
-        time_cur, m_bh1, num_bh1);
-        fclose(out);
-
-        out = fopen("mass-bh.con","w");
-        fprintf(out,"%.8E \t %.8E  %06d \n",
-        time_cur, m_bh1, num_bh1);
-        fclose(out);
-
-    #endif        // ADD_BH1
-
-    #endif        // ADD_BH2
-
-    #endif        // STARDESTR
-
-    #ifdef STARDESTR_EXT
-
-        out = fopen("phi-GRAPE.ext","w");
-        fprintf(out,"%.8E \t %.8E \n", m_bh, b_bh);
-        fclose(out);
-
-        out = fopen("mass-bh.dat","a");
-        fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh);
-        fclose(out);
-
-        out = fopen("mass-bh.con","w");
-        fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh);
-        fclose(out);
-
-    #endif        // STARDESTR_EXT
-
-        energy_contr();
-
-    #ifdef CMCORR111
-        for (i=0; i<N; i++)
-            {
-
-            for (k=0;k<3;k++)
-            {
-            x[i][k] -= xcm[k];
+            for (k=0; k<3; k++) {
+                x_act[i][k] =    x_act_new[i][k];
+                v_act[i][k] =    v_act_new[i][k];
+                a_act[i][k] =    a_act_new[i][k];
+                adot_act[i][k] = adot_act_new[i][k];
             } /* k */
 
+#ifdef DT_MIN_WARNING
+            if (myRank == 0) {
+                if (dt_act[i] == dt_min) {
+                    printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d \n", dt[i], ind_act[i]);
+                    fflush(stdout);
+                }
+            }
+#endif
+        } /* i */
+
+        /* define the min. dt over all the act. part. and set it also for the BH... */
+
+#ifdef ADD_BH2
+        min_dt = dt_act[0];
+
+        for (i=1; i<n_act; i++) {
+            if ( dt_act[i] < min_dt ) min_dt = dt_act[i];
+        } /* i */
+
+        dt_act[i_bh1] = min_dt;
+        dt_act[i_bh2] = min_dt;
+#else
+
+#ifdef ADD_BH1
+        min_dt = dt_act[0];
+
+        for (i=1; i<n_act; i++) {
+            if ( dt_act[i] < min_dt ) min_dt = dt_act[i];
+        } /* i */
+
+        dt_act[i_bh1] = min_dt;
+#endif        // ADD_BH1
+
+#endif        // ADD_BH2
+
+#ifdef BBH_INF
+        if (myRank == rootRank) {
+
+            out = fopen("bbh.inf","a");
+
+            m_bh1 = m_act[i_bh1];
+            m_bh2 = m_act[i_bh2];
+
+            for (k=0;k<3;k++) {
+                x_bh1[k] = x_act[i_bh1][k];
+                x_bh2[k] = x_act[i_bh2][k];
+
+                v_bh1[k] = v_act[i_bh1][k];
+                v_bh2[k] = v_act[i_bh2][k];
+            }
+
+            for (k=0;k<3;k++) {
+                x_bbhc[k] = (m_bh1*x_bh1[k] + m_bh2*x_bh2[k])/(m_bh1 + m_bh2);
+                v_bbhc[k] = (m_bh1*v_bh1[k] + m_bh2*v_bh2[k])/(m_bh1 + m_bh2);
+            }
+
+            DR2 = SQR(x_bh1[0] - x_bh2[0]) + SQR(x_bh1[1] - x_bh2[1]) + SQR(x_bh1[2] - x_bh2[2]);
+            DV2 = SQR(v_bh1[0] - v_bh2[0]) + SQR(v_bh1[1] - v_bh2[1]) + SQR(v_bh1[2] - v_bh2[2]);
+
+            EB = -(m_bh1 + m_bh2) / sqrt(DR2) + 0.5 * DV2;
+
+            SEMI_a = -0.5 * (m_bh1 + m_bh2)/EB;
+            SEMI_a2 = SQR(SEMI_a);
+
+            for (i=2; i<n_act; i++) {
+
+                tmp_r2 = SQR(x_act[i][0] - x_bbhc[0]) + SQR(x_act[i][1] - x_bbhc[1]) + SQR(x_act[i][2] - x_bbhc[2]);
+
+                iii = ind_act[i];
+
+                if ( (tmp_r2 < SEMI_a2*R_INF2) ) {
+
+                    if ( (inf_event[iii] == 0) ) {
+
+                        fprintf(out,"INF1 %.6E %.16E %07d %07d   %.6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E \n",
+                            Timesteps, time_cur, i, ind_act[i],
+                            sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2],
+                            m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[0],
+                            m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[1],
+                            sqrt(tmp_r2),
+                            m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i],
+                            dt_act[i]);
+
+                    inf_event[iii] = 1;
+                    }
+
+                } else {
+                    if ( (inf_event[iii] == 1) ) {
+                        fprintf(out,"INF2 %.6E %.16E %07d %07d   %.6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E \n",
+                            Timesteps, time_cur, i, ind_act[i],
+                            sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2],
+                            m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[0],
+                            m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[1],
+                            sqrt(tmp_r2),
+                            m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i],
+                            dt_act[i]);
+                    }
+
+                    inf_event[iii] = 0;
+
+                } /* if (tmp_r2 < DR2*R_INF2) */
+
             } /* i */
-    #endif
 
-
-
-        /* write cont data */
-
-        write_cont_data();
-
-        /* possible OUT for timing !!! */
-
-    #ifdef TIMING
-        out = fopen("timing.dat","a");
-
-        DT_TOT = DT_ACT_DEF1 + DT_ACT_DEF2 + DT_ACT_DEF3 + DT_ACT_PRED +
-                DT_ACT_GRAV + DT_EXT_GRAV + DT_GMC_GRAV +
-                DT_GMC_GMC_GRAV + DT_EXT_GMC_GRAV +
-                DT_ACT_CORR + DT_ACT_LOAD +
-                DT_STEVOL + DT_STARDISK + DT_STARDESTR +
-                DT_ACT_REDUCE;
-
-        fprintf(out,"%.8E \t %.6E \t %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f \t %.3f \t %.8E %.8E %.8E \t %.8E %.8E %.8E \n",
-                time_cur, DT_TOT,
-                100.0*DT_ACT_DEF1/DT_TOT, 100.0*DT_ACT_DEF2/DT_TOT, 100.0*DT_ACT_DEF3/DT_TOT, 100.0*DT_ACT_PRED/DT_TOT,
-                100.0*DT_ACT_GRAV/DT_TOT, 100.0*DT_EXT_GRAV/DT_TOT, 100.0*DT_GMC_GRAV/DT_TOT,
-                100.0*DT_GMC_GMC_GRAV/DT_TOT, 100.0*DT_EXT_GMC_GRAV/DT_TOT,
-                100.0*DT_ACT_CORR/DT_TOT, 100.0*DT_ACT_LOAD/DT_TOT,
-                100.0*DT_STEVOL/DT_TOT, 100.0*DT_STARDISK/DT_TOT, 100.0*DT_STARDESTR/DT_TOT,
-                100.0*DT_ACT_REDUCE/DT_TOT,
-                CPU_time_real-CPU_time_real0, CPU_time_user-CPU_time_user0, CPU_time_syst-CPU_time_syst0,
-                Timesteps, n_act_sum, 57.0*N*n_act_sum/(CPU_time_user-CPU_time_user0)/1.0E+09);
-
-        fclose(out);
-    #endif
+            fclose(out);
 
         } /* if (myRank == rootRank) */
+#endif
 
+        /* Return back the new coordinates + etc... of active part. to the global data... */
 
+        for (i=0; i<n_act; i++) {
 
+            iii = ind_act[i];
 
-    /* possible coordinate & velocity limits for ALL particles !!! */
+            m[iii] = m_act[i];
 
-    #ifdef LIMITS_ALL
-
-        for (i=0; i<N; i++)
-        {
-
-        tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
-        tmp_rv = x[i][0]*v[i][0] + x[i][1]*v[i][1] + x[i][2]*v[i][2];
-
-        if ( (tmp_r > R_LIMITS_ALL) && (tmp_rv > 0.0) )
-            {
-
-    #ifdef EXTPOT
-    #endif
-
-            for (k=0;k<3;k++)
-            {
-            v[i][k]    *= 1.0E-01;
+            for (k=0;k<3;k++) {
+                x[iii][k] = x_act[i][k];
+                v[iii][k] = v_act[i][k];
             } /* k */
 
-            } /* if */
+            t[iii] = t_act[i];
+            dt[iii] = dt_act[i];
 
+            pot[iii] = pot_act[i];
+
+#ifdef EXTPOT
+            pot_ext[iii] = pot_act_ext[i];
+#endif
+
+            for (k=0;k<3;k++) {
+                a[iii][k] = a_act[i][k];
+                adot[iii][k] = adot_act[i][k];
+            } /* k */
         } /* i */
 
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
+        DT_ACT_CORR += (CPU_tmp_user - CPU_tmp_user0);
+#endif
 
+        /* STARDESTR routine on all the nodes */
 
-    for (j=0; j<N; j++)
-        {
+#ifdef STARDESTR
 
-        cur_rank = ind[j]/n_loc;
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
+#endif
 
-        if (myRank == cur_rank)
-        {
+#ifdef ADD_BH2
 
-        jjj = ind[j] - myRank*n_loc;
+        star_destr(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act,
+                            N, ind, m, x, v, pot, a, adot, t, &m_bh1, &num_bh1, i_bh1);
 
-        for (k=0;k<3;k++)
+        m_act[i_bh1] = m_bh1;
+
+        star_destr(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act,
+                            N, ind, m, x, v, pot, a, adot, t, &m_bh2, &num_bh2, i_bh2);
+
+        m_act[i_bh2] = m_bh2;
+
+#else
+
+#ifdef ADD_BH1
+
+        star_destr(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act,
+                            N, ind, m, x, v, pot, a, adot, t, &m_bh1, &num_bh1, i_bh1);
+
+        m_act[i_bh1] = m_bh1;
+
+#endif        // ADD_BH1
+
+#endif        // ADD_BH2
+
+#ifdef LIMITS
+        for (i=0; i<n_act; i++)
             {
-            a2by18[k] = 0.0; a1by6[k] = adot[j][k]*over6; aby2[k] = a[j][k]*over2;
+            if (m_act[i] == 0.0)
+            {
+            x_act[i][0] = R_LIMITS + 1.0*my_rand2();
+            x_act[i][1] = R_LIMITS + 1.0*my_rand2();
+            x_act[i][2] = R_LIMITS + 1.0*my_rand2();
+
+            v_act[i][0] = 1.0E-06*my_rand2();
+            v_act[i][1] = 1.0E-06*my_rand2();
+            v_act[i][2] = 1.0E-06*my_rand2();
             }
+            }
+#endif        // LIMITS
 
-    #ifdef SAPPORO
-        g6_set_j_particle_(&clusterid, &jjj, &ind[j], &t[j], &dt[j], &m[j], a2by18, a1by6, aby2, v[j], x[j]);
-    #else
-        g6_set_j_particle(clusterid, jjj, ind[j], t[j], dt[j], m[j], a2by18, a1by6, aby2, v[j], x[j]);
-    #endif
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
+        DT_STARDESTR += (CPU_tmp_user - CPU_tmp_user0);
+#endif
 
-        }
+#endif        // STARDESTR
 
-        }
+#ifdef STARDESTR_EXT
 
-    #ifdef AMD
-        g6_flush_jp_buffer(clusterid);
-    #endif
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
+#endif
 
+        star_destr_ext(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act,
+                            N, m, x, v, a, adot, t, &m_bh, &num_bh);
 
-    #endif        // LIMITS_ALL
-
-
-
-
-    #ifdef CMCORR111
-
-    /* Wait to all processors to finish his works... */
-    MPI_Barrier(MPI_COMM_WORLD);
-
-    /* Broadcast the values of all particles to all processors... */
-    MPI_Bcast(x, 3*N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-    MPI_Bcast(v, 3*N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-
-    /* Scatter the "local" vectors from "global" */
-    MPI_Scatter(x, 3*n_loc, MPI_DOUBLE, x_loc, 3*n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-    MPI_Scatter(v, 3*n_loc, MPI_DOUBLE, v_loc, 3*n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-
-    /* Wait to all processors to finish his works... */
-    MPI_Barrier(MPI_COMM_WORLD);
-
-    #endif
-
-    #ifdef ETICS_CEP
-    // We are /inside/ a control step, so all particles must be synchronized; we can safely calculate their density centre. The acceleration and jerk currently in the memory are for the predicted position of the CEP, by calling grapite_calc_center we "correct" the position and velocity, but not the gravity at that point.
-
-    // First calculate the DC
-    grapite_calc_center(N, m, x, v, xcm, vcm, xdc, vdc);
-
-    // Now copy it to the global particle list
-    memcpy(x[grapite_cep_index], xdc, 3*sizeof(double));
-    memcpy(v[grapite_cep_index], vdc, 3*sizeof(double));
-
-    // Now copy it to the local particle list for tha appropriate rank
-    if (myRank == grapite_cep_index/n_loc) {
-        memcpy(x_loc[grapite_cep_index - myRank*n_loc], xdc, 3*sizeof(double));
-        memcpy(v_loc[grapite_cep_index - myRank*n_loc], vdc, 3*sizeof(double));
-    }
-    grapite_update_cep(time_cur, xdc, vdc, a[grapite_cep_index], adot[grapite_cep_index]);
-    #endif
-
-        t_contr += dt_contr;
-        } /* if (time_cur >= t_contr) */
-
-
-
-
-
-
-
-
-        if (time_cur >= t_disk)
-        {
-
-
-        if (myRank == rootRank)
-        {
-
-        diskstep++;
-
-        write_snap_data();
-
-        } /* if (myRank == rootRank) */
-
-    #ifdef ETICS_DUMP
-        sprintf(out_fname, "coeffs.%06d.%02d.dat", diskstep, myRank);
-        grapite_dump(out_fname, 2);
-    #endif
-
-
-
-
-    #ifdef LIMITS_NEW
-
-        for (i=0; i<N; i++)
-        {
-
-        tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
-
-        if ( tmp_r > 900.0 )
+#ifdef LIMITS
+        for (i=0; i<n_act; i++)
             {
-
-            m[i] = 0.0;
-
-            pot[i] = 0.0;
-
-    #ifdef EXTPOT
-            pot_ext[i] = 0.0;
-    #endif
-
-            for (k=0;k<3;k++)
+            if (m_act[i] == 0.0)
             {
-            x[i][k]    = 1000.0 + 1.0*my_rand2();
-            v[i][k]    = 1.0E-06*my_rand2();
-            a[i][k]    = 1.0E-06*my_rand2();
-            adot[i][k] = 1.0E-06*my_rand2();
+            x_act[i][0] = R_LIMITS + 1.0*my_rand2();
+            x_act[i][1] = R_LIMITS + 1.0*my_rand2();
+            x_act[i][2] = R_LIMITS + 1.0*my_rand2();
+
+            v_act[i][0] = 1.0E-06*my_rand2();
+            v_act[i][1] = 1.0E-06*my_rand2();
+            v_act[i][2] = 1.0E-06*my_rand2();
+            }
+            } /* i */
+#endif        // LIMITS
+
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
+        DT_STARDESTR += (CPU_tmp_user - CPU_tmp_user0);
+#endif
+
+#endif        // STARDESTR_EXT
+
+#if defined(STARDESTR) || defined(STARDESTR_EXT)
+
+        /* Return back the new coordinates + etc... of active part. to the global data... */
+
+        for (i=0; i<n_act; i++) {
+            iii = ind_act[i];
+
+            m[iii] = m_act[i];
+
+            for (k=0;k<3;k++) {
+                x[iii][k] = x_act[i][k];
+                v[iii][k] = v_act[i][k];
             } /* k */
 
-            t[i]  = 2.0*t_end;
+            t[iii] = t_act[i];
+            dt[iii] = dt_act[i];
 
-            dt[i] = 0.125;
+            pot[iii] = pot_act[i];
 
-            } /* if */
+#ifdef EXTPOT
+            pot_ext[iii] = pot_act_ext[i];
+#endif
+
+            for (k=0;k<3;k++) {
+                a[iii][k] = a_act[i][k];
+                adot[iii][k] = adot_act[i][k];
+            } /* k */
 
         } /* i */
 
+#endif        //#if defined(STARDESTR) || defined(STARDESTR_EXT)
 
+        /* load the new values for active particles to the local GRAPE's */
 
-    for (j=0; j<N; j++)
-        {
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
+#endif
 
-        cur_rank = ind[j]/n_loc;
+        for (j=0; j<n_act; j++) {
+#ifdef ETICS_CEP
+            if (ind_act[j] == grapite_cep_index) grapite_update_cep(t_act[j], x_act[j], v_act[j], a_act[j], adot_act[j]); // All ranks should do it.
+#endif
+            cur_rank = ind_act[j]/n_loc;
 
-        if (myRank == cur_rank)
-        {
+            if (myRank == cur_rank) {
 
-        jjj = ind[j] - myRank*n_loc;
+                jjj = ind_act[j] - myRank*n_loc;
 
-        for (k=0;k<3;k++)
-            {
-            a2by18[k] = 0.0; a1by6[k] = adot[j][k]*over6; aby2[k] = a[j][k]*over2;
-            } /* k */
+                for (k=0;k<3;k++) {
+                    a2by18[k] = 0.0; a1by6[k] = adot_act[j][k]*over6; aby2[k] = a_act[j][k]*over2;
+                } /* k */
 
-    #ifdef SAPPORO
-        g6_set_j_particle_(&clusterid, &jjj, &ind[j], &t[j], &dt[j], &m[j], a2by18, a1by6, aby2, v[j], x[j]);
-    #else
-        g6_set_j_particle(clusterid, jjj, ind[j], t[j], dt[j], m[j], a2by18, a1by6, aby2, v[j], x[j]);
-    #endif
+#ifdef SAPPORO
+                g6_set_j_particle_(&clusterid, &jjj, &ind_act[j], &t_act[j], &dt_act[j], &m_act[j], a2by18, a1by6, aby2, v_act[j], x_act[j]);
+#else
+                g6_set_j_particle(clusterid, jjj, ind_act[j], t_act[j], dt_act[j], m_act[j], a2by18, a1by6, aby2, v_act[j], x_act[j]);
+#endif
 
-        } /* if (myRank == cur_rank) */
+            } /* if (myRank == cur_rank) */
 
         } /* j */
 
-    #ifdef AMD
+#ifdef AMD
         g6_flush_jp_buffer(clusterid);
+#endif
+
+#ifdef TIMING
+        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
+        DT_ACT_LOAD += (CPU_tmp_user - CPU_tmp_user0);
+#endif
+
+        /* Current time set to min_t */
+
+        time_cur = min_t;
+
+        Timesteps += 1.0;
+        n_act_sum += n_act;
+        n_act_distr[n_act-1]++;
+
+        if (time_cur >= t_bh) {
+            if (myRank == rootRank) {
+#ifdef BH_OUT
+                /* Write BH data... */
+                write_bh_data();
+#endif
+
+#ifdef BH_OUT_NB
+                /* Write BH NB data... */
+                write_bh_nb_data();
+#endif
+
+#ifdef BH_OUT_NB_EXT
+                /* Write BH NB data... */
+                write_bh_nb_data_ext();
+#endif
+            } /* if (myRank == rootRank) */
+
+            t_bh += dt_bh;
+        } /* if (time_cur >= t_bh) */
+
+        if (time_cur >= t_contr) {
+            if (myRank == rootRank) {
+#ifdef STARDESTR
+
+#ifdef ADD_BH2
+
+            out = fopen("mass-bh.dat","a");
+            fprintf(out,"%.8E \t %.8E  %06d \t %.8E  %06d \n",
+            time_cur, m_bh1, num_bh1, m_bh2, num_bh2);
+            fclose(out);
+
+            out = fopen("mass-bh.con","w");
+            fprintf(out,"%.8E \t %.8E  %06d \t %.8E  %06d \n",
+            time_cur, m_bh1, num_bh1, m_bh2, num_bh2);
+            fclose(out);
+
+#else
+
+#ifdef ADD_BH1
+
+            out = fopen("mass-bh.dat","a");
+            fprintf(out,"%.8E \t %.8E  %06d \n",
+            time_cur, m_bh1, num_bh1);
+            fclose(out);
+
+            out = fopen("mass-bh.con","w");
+            fprintf(out,"%.8E \t %.8E  %06d \n",
+            time_cur, m_bh1, num_bh1);
+            fclose(out);
+
+#endif        // ADD_BH1
+
+#endif        // ADD_BH2
+
+#endif        // STARDESTR
+
+#ifdef STARDESTR_EXT
+
+            out = fopen("phi-GRAPE.ext","w");
+            fprintf(out,"%.8E \t %.8E \n", m_bh, b_bh);
+            fclose(out);
+
+            out = fopen("mass-bh.dat","a");
+            fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh);
+            fclose(out);
+
+            out = fopen("mass-bh.con","w");
+            fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh);
+            fclose(out);
+
+#endif        // STARDESTR_EXT
+
+                energy_contr();
+
+#ifdef CMCORR111
+                for (i=0; i<N; i++) {
+                    for (k=0;k<3;k++) {
+                        x[i][k] -= xcm[k];
+                    } /* k */
+                } /* i */
+#endif
+
+                /* write cont data */
+
+                write_cont_data();
+
+                /* possible OUT for timing !!! */
+
+#ifdef TIMING
+                out = fopen("timing.dat","a");
+
+                DT_TOT = DT_ACT_DEF1 + DT_ACT_DEF2 + DT_ACT_DEF3 + DT_ACT_PRED +
+                         DT_ACT_GRAV + DT_EXT_GRAV + DT_GMC_GRAV +
+                         DT_GMC_GMC_GRAV + DT_EXT_GMC_GRAV +
+                         DT_ACT_CORR + DT_ACT_LOAD +
+                         DT_STEVOL + DT_STARDISK + DT_STARDESTR +
+                         DT_ACT_REDUCE;
+
+                fprintf(out,"%.8E \t %.6E \t %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f \t %.3f \t %.8E %.8E %.8E \t %.8E %.8E %.8E \n",
+                        time_cur, DT_TOT,
+                        100.0*DT_ACT_DEF1/DT_TOT, 100.0*DT_ACT_DEF2/DT_TOT, 100.0*DT_ACT_DEF3/DT_TOT, 100.0*DT_ACT_PRED/DT_TOT,
+                        100.0*DT_ACT_GRAV/DT_TOT, 100.0*DT_EXT_GRAV/DT_TOT, 100.0*DT_GMC_GRAV/DT_TOT,
+                        100.0*DT_GMC_GMC_GRAV/DT_TOT, 100.0*DT_EXT_GMC_GRAV/DT_TOT,
+                        100.0*DT_ACT_CORR/DT_TOT, 100.0*DT_ACT_LOAD/DT_TOT,
+                        100.0*DT_STEVOL/DT_TOT, 100.0*DT_STARDISK/DT_TOT, 100.0*DT_STARDESTR/DT_TOT,
+                        100.0*DT_ACT_REDUCE/DT_TOT,
+                        CPU_time_real-CPU_time_real0, CPU_time_user-CPU_time_user0, CPU_time_syst-CPU_time_syst0,
+                        Timesteps, n_act_sum, 57.0*N*n_act_sum/(CPU_time_user-CPU_time_user0)/1.0E+09);
+
+                fclose(out);
+#endif
+
+            } /* if (myRank == rootRank) */
+
+            /* possible coordinate & velocity limits for ALL particles !!! */
+
+#ifdef LIMITS_ALL
+
+            for (i=0; i<N; i++) {
+
+                tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
+                tmp_rv = x[i][0]*v[i][0] + x[i][1]*v[i][1] + x[i][2]*v[i][2];
+
+                if ( (tmp_r > R_LIMITS_ALL) && (tmp_rv > 0.0) ) {
+                    for (k=0;k<3;k++) {
+                        v[i][k]    *= 1.0E-01;
+                    } /* k */
+                } /* if */
+            } /* i */
+
+            for (j=0; j<N; j++) {
+                cur_rank = ind[j]/n_loc;
+
+                if (myRank == cur_rank) {
+
+                    jjj = ind[j] - myRank*n_loc;
+
+                    for (k=0;k<3;k++) {
+                        a2by18[k] = 0.0; a1by6[k] = adot[j][k]*over6; aby2[k] = a[j][k]*over2;
+                    }
+
+    #ifdef SAPPORO
+                    g6_set_j_particle_(&clusterid, &jjj, &ind[j], &t[j], &dt[j], &m[j], a2by18, a1by6, aby2, v[j], x[j]);
+    #else
+                    g6_set_j_particle(clusterid, jjj, ind[j], t[j], dt[j], m[j], a2by18, a1by6, aby2, v[j], x[j]);
     #endif
+                }
+            }
 
+#ifdef AMD
+            g6_flush_jp_buffer(clusterid);
+#endif
 
+#endif        // LIMITS_ALL
 
-    #endif        // LIMITS_NEW
+#ifdef CMCORR111
+            /* Wait to all processors to finish his works... */
+            MPI_Barrier(MPI_COMM_WORLD);
 
+            /* Broadcast the values of all particles to all processors... */
+            MPI_Bcast(x, 3*N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
+            MPI_Bcast(v, 3*N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
 
+            /* Scatter the "local" vectors from "global" */
+            MPI_Scatter(x, 3*n_loc, MPI_DOUBLE, x_loc, 3*n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
+            MPI_Scatter(v, 3*n_loc, MPI_DOUBLE, v_loc, 3*n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
 
-        t_disk += dt_disk;
+            /* Wait to all processors to finish his works... */
+            MPI_Barrier(MPI_COMM_WORLD);
+#endif
+
+#ifdef ETICS_CEP
+            // We are /inside/ a control step, so all particles must be synchronized; we can safely calculate their density centre. The acceleration and jerk currently in the memory are for the predicted position of the CEP, by calling grapite_calc_center we "correct" the position and velocity, but not the gravity at that point.
+
+            // First calculate the DC
+            grapite_calc_center(N, m, x, v, xcm, vcm, xdc, vdc);
+
+            // Now copy it to the global particle list
+            memcpy(x[grapite_cep_index], xdc, 3*sizeof(double));
+            memcpy(v[grapite_cep_index], vdc, 3*sizeof(double));
+
+            // Now copy it to the local particle list for tha appropriate rank
+            if (myRank == grapite_cep_index/n_loc) {
+                memcpy(x_loc[grapite_cep_index - myRank*n_loc], xdc, 3*sizeof(double));
+                memcpy(v_loc[grapite_cep_index - myRank*n_loc], vdc, 3*sizeof(double));
+            }
+            grapite_update_cep(time_cur, xdc, vdc, a[grapite_cep_index], adot[grapite_cep_index]);
+#endif
+
+            t_contr += dt_contr;
+        } /* if (time_cur >= t_contr) */
+
+        if (time_cur >= t_disk) {
+            if (myRank == rootRank) {
+                diskstep++;
+                write_snap_data();
+            } /* if (myRank == rootRank) */
+
+#ifdef ETICS_DUMP
+            sprintf(out_fname, "coeffs.%06d.%02d.dat", diskstep, myRank);
+            grapite_dump(out_fname, 2);
+#endif
+
+#ifdef LIMITS_NEW
+
+            for (i=0; i<N; i++)
+            {
+
+            tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
+
+            if ( tmp_r > 900.0 )
+                {
+
+                m[i] = 0.0;
+
+                pot[i] = 0.0;
+
+#ifdef EXTPOT
+                pot_ext[i] = 0.0;
+#endif
+
+                for (k=0;k<3;k++)
+                {
+                x[i][k]    = 1000.0 + 1.0*my_rand2();
+                v[i][k]    = 1.0E-06*my_rand2();
+                a[i][k]    = 1.0E-06*my_rand2();
+                adot[i][k] = 1.0E-06*my_rand2();
+                } /* k */
+
+                t[i]  = 2.0*t_end;
+
+                dt[i] = 0.125;
+
+                } /* if */
+
+            } /* i */
+
+        for (j=0; j<N; j++)
+            {
+
+            cur_rank = ind[j]/n_loc;
+
+            if (myRank == cur_rank)
+            {
+
+            jjj = ind[j] - myRank*n_loc;
+
+            for (k=0;k<3;k++)
+                {
+                a2by18[k] = 0.0; a1by6[k] = adot[j][k]*over6; aby2[k] = a[j][k]*over2;
+                } /* k */
+
+#ifdef SAPPORO
+            g6_set_j_particle_(&clusterid, &jjj, &ind[j], &t[j], &dt[j], &m[j], a2by18, a1by6, aby2, v[j], x[j]);
+#else
+            g6_set_j_particle(clusterid, jjj, ind[j], t[j], dt[j], m[j], a2by18, a1by6, aby2, v[j], x[j]);
+#endif
+
+            } /* if (myRank == cur_rank) */
+
+            } /* j */
+
+#ifdef AMD
+            g6_flush_jp_buffer(clusterid);
+#endif
+
+#endif        // LIMITS_NEW
+
+            t_disk += dt_disk;
         } /* if (time_cur >= t_disk) */
 
-
-
-    #ifdef CPU_TIMELIMIT
-    if (myRank == rootRank)
-        {
-        get_CPU_time(&CPU_time_real, &CPU_time_user, &CPU_time_syst);
-        tmp_cpu = CPU_time_real-CPU_time_real0;
+#ifdef CPU_TIMELIMIT
+        if (myRank == rootRank) {
+            get_CPU_time(&CPU_time_real, &CPU_time_user, &CPU_time_syst);
+            tmp_cpu = CPU_time_real-CPU_time_real0;
         } /* if (myRank == rootRank) */
-    #endif
-
-
-    #ifdef ACT_DEF_LL
-    ModifyLinkList();
-
-    #endif
+#endif
 
+#ifdef ACT_DEF_LL
+        ModifyLinkList();
+#endif
     } /* while (time_cur < t_end) */
 
-
-
     /* close the local GRAPE's */
 
-    #ifdef SAPPORO
+#ifdef SAPPORO
     g6_close_(&clusterid);
-    #else
+#else
     g6_close(clusterid);
-    #endif
-
+#endif
 
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
@@ -4471,9 +4087,7 @@ int main(int argc, char *argv[])
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
 
-
-    if (myRank == rootRank)
-        {
+    if (myRank == rootRank) {
 
         /* Write some output for the timestep annalize... */
 
@@ -4485,18 +4099,12 @@ int main(int argc, char *argv[])
 
     } /* if (myRank == rootRank) */
 
-
-
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
 
     /* Finalize the MPI work */
     MPI_Finalize();
 
-    return(0);
+    return 0;
 
 }
-
-/****************************************************************************/
-/****************************************************************************/
-/****************************************************************************/