diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
index 1b81cbc213a0e12111d35e5efde31606d60fcd59..df72c9f56007513e33e814dc4bba94ca79e5c528 100644
--- a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
@@ -36,6 +36,9 @@ Statistics:
   time_first:          0.01 # Time of the first stat dump (non-cosmological run) (in internal units)
   delta_time:          1.05 # Time between statistics output
 
+Scheduler:
+  links_per_tasks: 25
+
 # Parameters for the self-gravity scheme
 Gravity:
   eta:                    0.025     # Constant dimensionless multiplier for time integration.
diff --git a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
index 569dbdd1487faea311318fa508d5a14832f074d4..57c366e79690d7a7111f81f38850b7e1be71d652 100644
--- a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
@@ -49,7 +49,7 @@ Gravity:
   theta:                  0.7     # Opening angle (Multipole acceptance criterion)
   comoving_softening:     0.0026994 # Comoving softening length (in internal units).
   max_physical_softening: 0.0007    # Physical softening length (in internal units).
-  mesh_side_length:       64
+  mesh_side_length:       256
 
 # Parameters for the hydrodynamics scheme
 SPH:
@@ -85,6 +85,9 @@ EAGLECooling:
   He_reion_z_sigma:        0.5
   He_reion_eV_p_H:         2.0
 
+Scheduler:
+  links_per_tasks:    25
+
 # EAGLE star formation parameters
 EAGLEStarFormation:
   EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
diff --git a/src/cell.c b/src/cell.c
index 99425f2987808f938772eb72eecc9e7d4e8382b1..d7513c8638f659dc30557c0f5b04788edfdc6f1b 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -160,6 +160,11 @@ struct cell_split_pair cell_split_pairs[13] = {
       {5, 6, 9},
       {7, 6, 12}}}};
 
+extern int cell_to_check;
+extern int parent_cell_to_check;
+extern int super_cell_to_check;
+int CHECK = 0;
+
 /**
  * @brief Get the size of the cell subtree.
  *
@@ -1072,7 +1077,8 @@ int cell_pack_sf_counts(struct cell *restrict c,
   /* Pack this cell's data. */
   pcells[0].stars.delta_from_rebuild = c->stars.parts - c->stars.parts_rebuild;
   pcells[0].stars.count = c->stars.count;
-
+  pcells[0].stars.dx_max_part = c->stars.dx_max_part;
+  
 #ifdef SWIFT_DEBUG_CHECKS
   if (c->stars.parts_rebuild == NULL)
     error("Star particles array at rebuild is NULL! c->depth=%d", c->depth);
@@ -1122,6 +1128,7 @@ int cell_unpack_sf_counts(struct cell *restrict c,
   /* Unpack this cell's data. */
   c->stars.count = pcells[0].stars.count;
   c->stars.parts = c->stars.parts_rebuild + pcells[0].stars.delta_from_rebuild;
+  c->stars.dx_max_part = pcells[0].stars.dx_max_part;
 
   /* Fill in the progeny, depth-first recursion. */
   int count = 1;
@@ -2541,7 +2548,11 @@ void cell_activate_hydro_sorts(struct cell *c, int sid, struct scheduler *s) {
  * @brief Activate the sorts up a cell hierarchy.
  */
 void cell_activate_stars_sorts_up(struct cell *c, struct scheduler *s) {
+
   if (c == c->hydro.super) {
+
+    if(CHECK) message("in if");
+
 #ifdef SWIFT_DEBUG_CHECKS
     if (c->stars.sorts == NULL)
       error("Trying to activate un-existing c->stars.sorts");
@@ -2551,11 +2562,31 @@ void cell_activate_stars_sorts_up(struct cell *c, struct scheduler *s) {
       cell_activate_drift_spart(c, s);
     }
   } else {
+    if(CHECK) message("in else");
+
+    if(CHECK) {
+
+    int sub_sort = cell_get_flag(c, cell_flag_do_stars_sub_sort);
+    int parent_sub_sort = cell_get_flag(c->parent, cell_flag_do_stars_sub_sort);
+
+      message("depth=%d cellID=%d sub_sort=%d parent->sub_sort=%d",
+	      c->depth, c->nodeID, sub_sort, parent_sub_sort);
+      message("super->stars.sorts=%p", c->hydro.super->stars.sorts);
+      message("super->stars.sorts->skip=%d", c->hydro.super->stars.sorts->skip);
+    } 
+
     for (struct cell *parent = c->parent;
          parent != NULL && !cell_get_flag(parent, cell_flag_do_stars_sub_sort);
          parent = parent->parent) {
+
       cell_set_flag(parent, cell_flag_do_stars_sub_sort);
+      if(CHECK) message("parent->depth=%d", parent->depth);
+
       if (parent == c->hydro.super) {
+
+	if(CHECK) message("in other if");
+
+
 #ifdef SWIFT_DEBUG_CHECKS
         if (parent->stars.sorts == NULL)
           error("Trying to activate un-existing parents->stars.sorts");
@@ -2572,8 +2603,16 @@ void cell_activate_stars_sorts_up(struct cell *c, struct scheduler *s) {
  * @brief Activate the sorts on a given cell, if needed.
  */
 void cell_activate_stars_sorts(struct cell *c, int sid, struct scheduler *s) {
+
+  if(CHECK)
+    message("Activating sorts for cell %d", c->cellID);
+
   /* Do we need to re-sort? */
   if (c->stars.dx_max_sort > space_maxreldx * c->dmin) {
+
+    if(c->cellID == cell_to_check)
+      message("In first if");
+
     /* Climb up the tree to active the sorts in that direction */
     for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
       if (finger->stars.requires_sorts) {
@@ -2586,6 +2625,10 @@ void cell_activate_stars_sorts(struct cell *c, int sid, struct scheduler *s) {
 
   /* Has this cell been sorted at all for the given sid? */
   if (!(c->stars.sorted & (1 << sid)) || c->nodeID != engine_rank) {
+
+    if(CHECK)
+      message("In second if");
+
     atomic_or(&c->stars.do_sort, (1 << sid));
     cell_activate_stars_sorts_up(c, s);
   }
@@ -2744,21 +2787,52 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
   /* Otherwise, pair interation */
   else {
 
-    /* Get the orientation of the pair. */
-    double shift[3];
-    const int sid = space_getsid(s->space, &ci, &cj, shift);
-
     const int ci_active = cell_is_active_stars(ci, e) ||
-                          (with_star_formation && cell_is_active_hydro(ci, e));
+      (with_star_formation && cell_is_active_hydro(ci, e));
     const int cj_active = cell_is_active_stars(cj, e) ||
-                          (with_star_formation && cell_is_active_hydro(cj, e));
+      (with_star_formation && cell_is_active_hydro(cj, e));
 
     /* Should we even bother? */
     if (!ci_active && !cj_active) return;
 
+    /* Get the orientation of the pair. */
+    double shift[3];
+    const int sid = space_getsid(s->space, &ci, &cj, shift);
+    
+    /* if(e->nodeID == 7 && ci->cellID == super_cell_to_check) */
+    /*   message("Found the super ci! ci_active=%d cj_active=%d depth=%d", ci_active, cj_active, ci->depth); */
+
+    /* if(e->nodeID == 7 && cj->cellID == super_cell_to_check) */
+    /*   message("Found the super cj! ci_active=%d cj_active=%d depth=%d", ci_active, cj_active, cj->depth); */
+
+    /* if(e->nodeID == 7 && ci->depth > 0 && ci->parent->cellID == super_cell_to_check) */
+    /*   message("Found the parent ci! ci_active=%d cj_active=%d depth=%d", ci_active, cj_active, ci->depth); */
+
+    /* if(e->nodeID == 7 && cj->depth > 0 &&cj->parent->cellID == super_cell_to_check) */
+    /*   message("Found the parent cj! ci_active=%d cj_active=%d depth=%d", ci_active, cj_active, cj->depth); */
+
+    /* if(e->nodeID == 7 && cj->hydro.super->cellID == super_cell_to_check) */
+    /*   message("Found a cell with super-cell= %d depth=%d cellID=%d ci_active=%d cj_active=%d cj->requires_sorts=%d cj->do_sort=%d sid=%d cj->dx_max_part=%e cj->dx_max_part_old=%e cj->dx_max_sort=%e cj->dx_max_sort_old=%e", */
+    /* 	      cj->hydro.super->cellID, cj->depth, cj->cellID, ci_active, cj_active, cj->stars.requires_sorts, cj->stars.do_sort, sid, */
+    /* 	      cj->stars.dx_max_part, cj->stars.dx_max_part_old, */
+    /* 	      cj->stars.dx_max_sort, cj->stars.dx_max_sort_old); */
+       
     /* recurse? */
     if (cell_can_recurse_in_pair_stars_task(ci, cj) &&
         cell_can_recurse_in_pair_stars_task(cj, ci)) {
+
+      /* if(e->nodeID == 7 && ci->cellID == super_cell_to_check) */
+      /* 	message("Found the super ci! Recursing!"); */
+      
+      /* if(e->nodeID == 7 && cj->cellID == super_cell_to_check) */
+      /* 	message("Found the super cj! Recursing!"); */
+
+      /* if(e->nodeID == 7 && ci->depth > 0 && ci->parent->cellID == super_cell_to_check) */
+      /* 	message("Found the parent ci! Recursing!"); */
+      
+      /* if(e->nodeID == 7 && cj->depth > 0 && cj->parent->cellID == super_cell_to_check) */
+      /* 	message("Found the parent cj! Recursing!"); */
+
       const struct cell_split_pair *csp = &cell_split_pairs[sid];
       for (int k = 0; k < csp->count; k++) {
         const int pid = csp->pairs[k].pid;
@@ -2772,7 +2846,16 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
     /* Otherwise, activate the sorts and drifts. */
     else {
 
-      if (ci_active) {
+      if (cell_is_active_stars(ci, e) ||
+	  (with_star_formation && cell_is_active_hydro(ci, e))) {
+
+    /* if(e->nodeID == 7 && cj->hydro.super->cellID == super_cell_to_check && cj->depth==3) */
+    /*   message("ACTIVATING  ci!!! Found a cell with super-cell= %d depth=%d cellID=%d ci_active=%d cj_active=%d cj->requires_sorts=%d cj->do_sort=%d sid=%d cj->dx_max_part=%e cj->dx_max_part_old=%e cj->dx_max_sort=%e cj->dx_max_sort_old=%e", */
+    /* 	      cj->hydro.super->cellID, cj->depth, cj->cellID, ci_active, cj_active, cj->stars.requires_sorts, cj->stars.do_sort, sid, */
+    /* 	      cj->stars.dx_max_part, cj->stars.dx_max_part_old, */
+    /* 	      cj->stars.dx_max_sort, cj->stars.dx_max_sort_old); */
+
+
         /* We are going to interact this pair, so store some values. */
         atomic_or(&cj->hydro.requires_sorts, 1 << sid);
         atomic_or(&ci->stars.requires_sorts, 1 << sid);
@@ -2789,7 +2872,19 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
         cell_activate_stars_sorts(ci, sid, s);
       }
 
-      if (cj_active) {
+      if (cell_is_active_stars(cj, e) ||
+	  (with_star_formation && cell_is_active_hydro(cj, e))) {
+
+      /* 	if(e->nodeID == 7 && cj->hydro.super->cellID == super_cell_to_check && cj->depth==3) { */
+      /* message("ACTIVATING  cj!!! Found a cell with super-cell= %d depth=%d cellID=%d ci_active=%d cj_active=%d cj->requires_sorts=%d cj->do_sort=%d sid=%d cj->dx_max_part=%e cj->dx_max_part_old=%e cj->dx_max_sort=%e cj->dx_max_sort_old=%e", */
+      /* 	      cj->hydro.super->cellID, cj->depth, cj->cellID, ci_active, cj_active, cj->stars.requires_sorts, cj->stars.do_sort, sid, */
+      /* 	      cj->stars.dx_max_part, cj->stars.dx_max_part_old, */
+      /* 	      cj->stars.dx_max_sort, cj->stars.dx_max_sort_old); */
+
+      /* if(sid == 3) CHECK =1; */
+      /* 	} */
+
+
         /* We are going to interact this pair, so store some values. */
         atomic_or(&cj->stars.requires_sorts, 1 << sid);
         atomic_or(&ci->hydro.requires_sorts, 1 << sid);
@@ -2804,6 +2899,8 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
         /* Do we need to sort the cells? */
         cell_activate_hydro_sorts(ci, sid, s);
         cell_activate_stars_sorts(cj, sid, s);
+
+	CHECK = 0;
       }
     }
   } /* Otherwise, pair interation */
@@ -3512,6 +3609,19 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s,
         (ci_nodeID == nodeID || cj_nodeID == nodeID)) {
       scheduler_activate(s, t);
 
+      if(ci->cellID == cell_to_check)
+	message("Activating task ci case t->type=%s/%s", taskID_names[t->type], subtaskID_names[t->subtype]);
+
+      if(cj != NULL && cj->cellID == cell_to_check)
+	message("Activating task cj case t->type=%s/%s", taskID_names[t->type], subtaskID_names[t->subtype]);
+
+
+      if(ci->cellID == super_cell_to_check)
+	message("Activating super task ci case t->type=%s/%s", taskID_names[t->type], subtaskID_names[t->subtype]);
+
+      if(cj != NULL && cj->cellID == super_cell_to_check)
+	message("Activating super task cj case t->type=%s/%s", taskID_names[t->type], subtaskID_names[t->subtype]);
+
       if (t->type == task_type_pair) {
         /* Do ci */
         if (ci_active) {
@@ -4801,6 +4911,35 @@ void cell_check_spart_pos(const struct cell *c,
 #endif
 }
 
+void cell_check_sort_flags(const struct cell* c) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  const int do_hydro_sub_sort = cell_get_flag(c, cell_flag_do_hydro_sub_sort);
+  const int do_stars_sub_sort = cell_get_flag(c, cell_flag_do_stars_sub_sort);
+
+  if(do_hydro_sub_sort)
+    error("cell %d has a hydro sub_sort flag set. Node=%d depth=%d maxdepth=%d",
+          c->cellID, c->nodeID, c->depth, c->maxdepth);
+
+  if(do_stars_sub_sort) {
+    message("cell %d has a stars sub_sort flag set. Node=%d depth=%d maxdepth=%d super=%p",
+             c->cellID, c->nodeID, c->depth, c->maxdepth, c->hydro.super);
+   message("c->stars.count=%d", c->stars.count);
+   message("super->cellID=%d super->sorts=%p super->depth=%d", c->hydro.super->cellID,
+c->hydro.super->stars.sorts, c->hydro.super->depth);
+   message("super->sorts->skip=%d",c->hydro.super->stars.sorts->skip);
+
+  error("oooo");
+}
+
+  if(c->split) {
+    for(int k = 0; k < 8; ++k) {
+      if(c->progeny[k] != NULL) cell_check_sort_flags(c->progeny[k]);
+    } 
+  }
+#endif
+}
+
 /**
  * @brief Recursively update the pointer and counter for #spart after the
  * addition of a new particle.
@@ -5190,6 +5329,11 @@ struct spart *cell_convert_part_to_spart(struct engine *e, struct cell *c,
   /* Did we run out of free spart slots? */
   if (sp == NULL) return NULL;
 
+  /* Copy over the distance since rebuild */
+  sp->x_diff[0] = xp->x_diff[0];
+  sp->x_diff[1] = xp->x_diff[1];
+  sp->x_diff[2] = xp->x_diff[2];
+
   /* Destroy the gas particle and get it's gpart friend */
   struct gpart *gp = cell_convert_part_to_gpart(e, c, p, xp);
 
diff --git a/src/cell.h b/src/cell.h
index 461ed45028e2363ec5c6143e913de78611cda8e5..c57b55cf9fef63abe4f5c9f74c1c6d8c9c734bb1 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -253,6 +253,26 @@ struct pcell_step_black_holes {
   float dx_max_part;
 };
 
+/**
+ * @brief Cell information to propagate the new counts of star particles.
+ */
+struct pcell_sf {
+
+  /*! Stars variables */
+  struct {
+
+    /* Distance by which the stars pointer has moved since the last rebuild */
+    ptrdiff_t delta_from_rebuild;
+
+    /* Number of particles in the cell */
+    int count;
+
+    /*! Maximum part movement in this cell since last construction. */
+    float dx_max_part;
+    
+  } stars;
+};
+
 /** Bitmasks for the cell flags. Beware when adding flags that you don't exceed
     the size of the flags variable in the struct cell. */
 enum cell_flags {
@@ -271,23 +291,6 @@ enum cell_flags {
   cell_flag_do_bh_sub_drift = (1UL << 12)
 };
 
-/**
- * @brief Cell information to propagate the new counts of star particles.
- */
-struct pcell_sf {
-
-  /*! Stars variables */
-  struct {
-
-    /* Distance by which the stars pointer has moved since the last rebuild */
-    ptrdiff_t delta_from_rebuild;
-
-    /* Number of particles in the cell */
-    int count;
-
-  } stars;
-};
-
 /**
  * @brief Cell within the tree structure.
  *
@@ -880,6 +883,7 @@ void cell_clear_limiter_flags(struct cell *c, void *data);
 void cell_set_super_mapper(void *map_data, int num_elements, void *extra_data);
 void cell_check_spart_pos(const struct cell *c,
                           const struct spart *global_sparts);
+void cell_check_sort_flags(const struct cell* c);
 void cell_clear_stars_sort_flags(struct cell *c);
 int cell_has_tasks(struct cell *c);
 void cell_remove_part(const struct engine *e, struct cell *c, struct part *p,
diff --git a/src/engine.c b/src/engine.c
index 36cf0709044c885ae48890bf9051a905f57839f1..0a21e3659fb4016f938f8fdffd4081e619065b62 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2561,6 +2561,8 @@ void engine_rebuild(struct engine *e, int repartitioned,
     for (int k = 0; k < e->s->nr_local_cells; k++)
       cell_check_foreign_multipole(&e->s->cells_top[e->s->local_cells_top[k]]);
   }
+
+  space_check_sort_flags(e->s);
 #endif
 
   /* Run through the tasks and mark as skip or not. */
@@ -3834,6 +3836,7 @@ void engine_step(struct engine *e) {
 #ifdef SWIFT_DEBUG_CHECKS
   /* Make sure all woken-up particles have been processed */
   space_check_limiter(e->s);
+  space_check_sort_flags(e->s);
 #endif
 
   /* Collect information about the next time-step */
@@ -4722,22 +4725,22 @@ void engine_dump_snapshot(struct engine *e) {
 #endif
 
 /* Dump... */
-#if defined(HAVE_HDF5)
-#if defined(WITH_MPI)
-#if defined(HAVE_PARALLEL_HDF5)
-  write_output_parallel(e, e->snapshot_base_name, e->internal_units,
-                        e->snapshot_units, e->nodeID, e->nr_nodes,
-                        MPI_COMM_WORLD, MPI_INFO_NULL);
-#else
-  write_output_serial(e, e->snapshot_base_name, e->internal_units,
-                      e->snapshot_units, e->nodeID, e->nr_nodes, MPI_COMM_WORLD,
-                      MPI_INFO_NULL);
-#endif
-#else
-  write_output_single(e, e->snapshot_base_name, e->internal_units,
-                      e->snapshot_units);
-#endif
-#endif
+/* #if defined(HAVE_HDF5) */
+/* #if defined(WITH_MPI) */
+/* #if defined(HAVE_PARALLEL_HDF5) */
+/*   write_output_parallel(e, e->snapshot_base_name, e->internal_units, */
+/*                         e->snapshot_units, e->nodeID, e->nr_nodes, */
+/*                         MPI_COMM_WORLD, MPI_INFO_NULL); */
+/* #else */
+/*   write_output_serial(e, e->snapshot_base_name, e->internal_units, */
+/*                       e->snapshot_units, e->nodeID, e->nr_nodes, MPI_COMM_WORLD, */
+/*                       MPI_INFO_NULL); */
+/* #endif */
+/* #else */
+/*   write_output_single(e, e->snapshot_base_name, e->internal_units, */
+/*                       e->snapshot_units); */
+/* #endif */
+/* #endif */
 
   /* Flag that we dumped a snapshot */
   e->step_props |= engine_step_prop_snapshot;
diff --git a/src/runner.c b/src/runner.c
index 858c024de1fb38103b9d99f8ce923288ed37a5d6..58105e9c9b8eedd91dd5ead8260354984ba52c59 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -142,6 +142,10 @@
 #undef FUNCTION_TASK_LOOP
 #undef FUNCTION
 
+int cell_to_check = -10000000;
+int parent_cell_to_check = -10000000;
+int super_cell_to_check = -10000000;
+
 /**
  * @brief Intermediate task after the density to check that the smoothing
  * lengths are correct.
@@ -1100,6 +1104,9 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
             /* Did we get a star? (Or did we run out of spare ones?) */
             if (sp != NULL) {
 
+              message("We formed a star id=%lld cellID=%d", sp->id,
+		      c->cellID);
+
               /* Copy the properties of the gas particle to the star particle */
               star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo,
                                              with_cosmology, phys_const,
@@ -1243,6 +1250,20 @@ void runner_do_sort_ascending(struct entry *sort, int N) {
 RUNNER_CHECK_SORTS(hydro)
 RUNNER_CHECK_SORTS(stars)
 
+void cell_clear_hydro_sort_flags2(struct cell *c)  {
+
+  c->hydro.do_sort = 0;
+  cell_clear_flag(c, cell_flag_do_hydro_sub_sort);
+  c->hydro.requires_sorts = 0;
+
+  if(c->split) {
+    for(int k = 0; k < 8; ++k){
+      if(c->progeny[k] != NULL)
+	cell_clear_hydro_sort_flags2(c->progeny[k]);
+    }
+  }
+}
+
 /**
  * @brief Sort the particles in the given cell along all cardinal directions.
  *
@@ -1308,15 +1329,20 @@ void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
     float dx_max_sort = 0.0f;
     float dx_max_sort_old = 0.0f;
     for (int k = 0; k < 8; k++) {
-      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
-        /* Only propagate cleanup if the progeny is stale. */
-        runner_do_hydro_sort(r, c->progeny[k], flags,
-                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
-                                         space_maxreldx * c->progeny[k]->dmin),
-                             0);
-        dx_max_sort = max(dx_max_sort, c->progeny[k]->hydro.dx_max_sort);
-        dx_max_sort_old =
+      if (c->progeny[k] != NULL) {
+	
+	if( c->progeny[k]->hydro.count > 0) {
+	  /* Only propagate cleanup if the progeny is stale. */
+	  runner_do_hydro_sort(r, c->progeny[k], flags,
+			       cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
+					   space_maxreldx * c->progeny[k]->dmin),
+			       0);
+	  dx_max_sort = max(dx_max_sort, c->progeny[k]->hydro.dx_max_sort);
+	  dx_max_sort_old =
             max(dx_max_sort_old, c->progeny[k]->hydro.dx_max_sort_old);
+	} else {
+	  cell_clear_hydro_sort_flags2(c->progeny[k]);
+	}
       }
     }
     c->hydro.dx_max_sort = dx_max_sort;
@@ -1465,6 +1491,21 @@ void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
   if (clock) TIMER_TOC(timer_dosort);
 }
 
+void cell_clear_stars_sort_flags2(struct cell *c)  {
+
+  c->stars.do_sort = 0;
+  cell_clear_flag(c, cell_flag_do_stars_sub_sort);
+  c->stars.requires_sorts = 0;
+
+  if(c->split) {
+    for(int k = 0; k < 8; ++k){
+      if(c->progeny[k] != NULL)
+	cell_clear_stars_sort_flags2(c->progeny[k]);
+    }
+  }
+}
+
+
 /**
  * @brief Sort the stars particles in the given cell along all cardinal
  * directions.
@@ -1491,6 +1532,16 @@ void runner_do_stars_sort(struct runner *r, struct cell *c, int flags,
   if (c->hydro.super == NULL) error("Task called above the super level!!!");
 #endif
 
+  if(c->cellID == cell_to_check) {
+    message("Sorting stars in cellID=%d count=%d is_super=%d super=%d", c->cellID,
+	    c->stars.count, c == c->hydro.super, c->hydro.super->cellID);
+  }
+
+  if(c->cellID == super_cell_to_check) {
+    message("Sorting stars in cellID=%d count=%d is_super=%d", c->cellID,
+	    c->stars.count, c == c->hydro.super);
+  }
+
   /* We need to do the local sorts plus whatever was requested further up. */
   flags |= c->stars.do_sort;
   if (cleanup) {
@@ -1531,7 +1582,8 @@ void runner_do_stars_sort(struct runner *r, struct cell *c, int flags,
     float dx_max_sort = 0.0f;
     float dx_max_sort_old = 0.0f;
     for (int k = 0; k < 8; k++) {
-      if (c->progeny[k] != NULL && c->progeny[k]->stars.count > 0) {
+      if (c->progeny[k] != NULL){ 
+	if(c->progeny[k]->stars.count > 0) {
         /* Only propagate cleanup if the progeny is stale. */
         const int cleanup_prog =
             cleanup && (c->progeny[k]->stars.dx_max_sort_old >
@@ -1540,6 +1592,9 @@ void runner_do_stars_sort(struct runner *r, struct cell *c, int flags,
         dx_max_sort = max(dx_max_sort, c->progeny[k]->stars.dx_max_sort);
         dx_max_sort_old =
             max(dx_max_sort_old, c->progeny[k]->stars.dx_max_sort_old);
+	} else {
+	  cell_clear_stars_sort_flags2(c->progeny[k]);
+	}
       }
     }
     c->stars.dx_max_sort = dx_max_sort;
@@ -1621,6 +1676,8 @@ void runner_do_stars_sort(struct runner *r, struct cell *c, int flags,
 
       /* And the individual sort distances if we are a local cell */
       for (int k = 0; k < count; k++) {
+	if(sparts[k].id == 155966626889L)
+	  message("Sorting star %lld", sparts[k].id);
         sparts[k].x_diff_sort[0] = 0.0f;
         sparts[k].x_diff_sort[1] = 0.0f;
         sparts[k].x_diff_sort[2] = 0.0f;
@@ -2349,8 +2406,11 @@ static void runner_do_unskip_hydro(struct cell *c, struct engine *e) {
 static void runner_do_unskip_stars(struct cell *c, struct engine *e,
                                    const int with_star_formation) {
 
+  const int non_empty = c->stars.count > 0 ||
+                        (with_star_formation && c->hydro.count > 0);
+
   /* Ignore empty cells. */
-  if (c->stars.count == 0) return;
+  if (!non_empty) return;
 
   const int ci_active = cell_is_active_stars(c, e) ||
                         (with_star_formation && cell_is_active_hydro(c, e));
@@ -2562,7 +2622,7 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
 
   /* Anything to do here? */
   if (!cell_is_starting_hydro(c, e) && !cell_is_starting_gravity(c, e) &&
-      !cell_is_starting_stars(c, e))
+      !cell_is_starting_stars(c, e) && !cell_is_starting_black_holes(c, e))
     return;
 
   /* Recurse? */
@@ -2746,7 +2806,7 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
 
   /* Anything to do here? */
   if (!cell_is_active_hydro(c, e) && !cell_is_active_gravity(c, e) &&
-      !cell_is_active_stars(c, e))
+      !cell_is_active_stars(c, e) && !cell_is_active_black_holes(c, e))
     return;
 
   /* Recurse? */
@@ -3117,6 +3177,9 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
       /* Get a handle on the part. */
       struct spart *restrict sp = &sparts[k];
 
+      if(sp->id == 155966626889L)
+	message("getting time-step for spart id=%lld", sp->id);
+
       /* need to be updated ? */
       if (spart_is_active(sp, e)) {
 
@@ -3135,6 +3198,10 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
         sp->time_bin = get_time_bin(ti_new_step);
         sp->gpart->time_bin = get_time_bin(ti_new_step);
 
+	if(sp->id == 155966626889L)
+	  message("new time-step for spart id=%lld is %d", sp->id,
+		  sp->time_bin);
+      
         /* Number of updated s-particles */
         s_updated++;
         g_updated++;
@@ -3818,6 +3885,18 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
       time_bin_min = min(time_bin_min, sparts[k].time_bin);
       time_bin_max = max(time_bin_max, sparts[k].time_bin);
       h_max = max(h_max, sparts[k].h);
+
+      if(sparts[k].id == 155966626889L) {
+	message("Received star %lld", sparts[k].id);
+	cell_to_check = c->cellID;
+	parent_cell_to_check = c->parent->cellID;
+	super_cell_to_check = c->hydro.super->cellID;
+	message("Cell to check: %d depth=%d", cell_to_check, c->depth);
+	message("Parent cell to check: %d", parent_cell_to_check);
+	message("Super to check: %d", super_cell_to_check);
+	message("Super to check sort task: %p", c->hydro.super->stars.sorts);
+	message("Super to check sort task skip: %d", c->hydro.super->stars.sorts->skip);
+      }
     }
 
     /* Convert into a time */
@@ -3833,7 +3912,7 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
         ti_stars_end_min =
             min(ti_stars_end_min, c->progeny[k]->stars.ti_end_min);
         ti_stars_end_max =
-            max(ti_stars_end_max, c->progeny[k]->grav.ti_end_max);
+            max(ti_stars_end_max, c->progeny[k]->stars.ti_end_max);
         h_max = max(h_max, c->progeny[k]->stars.h_max);
       }
     }
@@ -3841,7 +3920,7 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
 
 #ifdef SWIFT_DEBUG_CHECKS
   if (ti_stars_end_min < ti_current &&
-      !(e->policy & engine_policy_star_formation))
+      !(r->e->policy & engine_policy_star_formation))
     error(
         "Received a cell at an incorrect time c->ti_end_min=%lld, "
         "e->ti_current=%lld.",
@@ -3920,7 +3999,7 @@ void runner_do_recv_bpart(struct runner *r, struct cell *c, int clear_sorts,
         ti_black_holes_end_min =
             min(ti_black_holes_end_min, c->progeny[k]->black_holes.ti_end_min);
         ti_black_holes_end_max =
-            max(ti_black_holes_end_max, c->progeny[k]->grav.ti_end_max);
+            max(ti_black_holes_end_max, c->progeny[k]->black_holes.ti_end_max);
         h_max = max(h_max, c->progeny[k]->black_holes.h_max);
       }
     }
diff --git a/src/runner_doiact_stars.h b/src/runner_doiact_stars.h
index be4f83f8137a1908f0d3dc7f321b2d11ea28bd3d..822f57a1843fdc5a25e402dd44b461d00ec22419 100644
--- a/src/runner_doiact_stars.h
+++ b/src/runner_doiact_stars.h
@@ -1308,7 +1308,8 @@ void DOSUB_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
       /* Do any of the cells need to be sorted first? */
       if (!(ci->stars.sorted & (1 << sid)) ||
           ci->stars.dx_max_sort_old > ci->dmin * space_maxreldx) {
-        error("Interacting unsorted cell (sparts).");
+        error("Interacting unsorted cell (sparts). ci->nodeID=%d ci->stars.count=%d ci->stars[0].id=%lld ci->cellID=%d ci->super->cellID=%d", 
+	      ci->nodeID, ci->stars.count, ci->stars.parts[0].id, ci->cellID, ci->hydro.super->cellID);
       }
 
       if (!(cj->hydro.sorted & (1 << sid)) ||
@@ -1333,7 +1334,11 @@ void DOSUB_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
 
       if (!(cj->stars.sorted & (1 << sid)) ||
           cj->stars.dx_max_sort_old > cj->dmin * space_maxreldx) {
-        error("Interacting unsorted cell (sparts).");
+        error("Interacting unsorted cell (sparts). cj->nodeID=%d cj->stars.count=%d cj->stars[0].id=%lld cj->cellID=%d cj->parent->cellID=%d cj->super->cellID=%d stars.ti_end_min=%lld ti_current=%lld cj->depth=%d cj->super->depth=%d cj->requires_sorts=%d cj->do_sort=%d sid=%d dx_max_part=%e dx_max_part_old=%e dx_max_sort=%e dx_max_sort_old=%e", 
+	      cj->nodeID, cj->stars.count, cj->stars.parts[0].id, cj->cellID, cj->parent->cellID, cj->hydro.super->cellID, cj->stars.ti_end_min,
+	      e->ti_current, cj->depth, cj->hydro.super->depth, cj->stars.requires_sorts, cj->stars.do_sort, sid,
+	      cj->stars.dx_max_part, cj->stars.dx_max_part_old,
+	      cj->stars.dx_max_sort, cj->stars.dx_max_sort_old);
       }
     }
 
diff --git a/src/space.c b/src/space.c
index bdf71370c245af7946419cec4b61f5db4a191693..b5765a01e7ce0dbca86aadf9ea19b2b679c27e78 100644
--- a/src/space.c
+++ b/src/space.c
@@ -5192,6 +5192,42 @@ void space_check_limiter(struct space *s) {
 #endif
 }
 
+void space_check_sort_flags_mapper(void *map_data, int nr_cells,
+				   void *extra_data) {
+  
+#ifdef SWIFT_DEBUG_CHECKS
+  
+const struct space *s = (struct space*) extra_data;
+ int *local_cells_top = map_data;
+
+ for(int ind = 0; ind < nr_cells; ++ind) {
+  const struct cell *c = &s->cells_top[local_cells_top[ind]];
+
+  cell_check_sort_flags(c);
+}
+  
+#endif
+}
+
+/**
+ * @brief Checks that all cells have cleared their sort flags.
+ *
+ * Should only be used for debugging purposes.
+ *
+ * @param s The #space to check.
+ */
+void space_check_sort_flags(struct space *s) {
+#ifdef SWIFT_DEBUG_CHECKS
+
+  threadpool_map(&s->e->threadpool, space_check_sort_flags_mapper, 
+		 s->local_cells_with_tasks_top,
+                 s->nr_local_cells_with_tasks, 
+		 sizeof(int), 1, s);
+#else
+  error("Calling debugging code without debugging flag activated.");
+#endif
+}
+
 /**
  * @brief Resets all the individual cell task counters to 0.
  *
diff --git a/src/space.h b/src/space.h
index 4a2d5d8ce92d49ef129fc32c9332bc811e67f795..88f38209e7d8d82b021cd518dfa9806e4ec844eb 100644
--- a/src/space.h
+++ b/src/space.h
@@ -354,6 +354,7 @@ void space_check_top_multipoles_drift_point(struct space *s,
                                             integertime_t ti_drift);
 void space_check_timesteps(struct space *s);
 void space_check_limiter(struct space *s);
+void space_check_sort_flags(struct space *s);
 void space_replicate(struct space *s, int replicate, int verbose);
 void space_generate_gas(struct space *s, const struct cosmology *cosmo,
                         int periodic, const double dim[3], int verbose);
diff --git a/src/stars/EAGLE/stars.h b/src/stars/EAGLE/stars.h
index 09926d6a8bdd5af4eff120a9a6c57d15eb7e49c4..1a8564463cced21aaef79968f645153669425bf1 100644
--- a/src/stars/EAGLE/stars.h
+++ b/src/stars/EAGLE/stars.h
@@ -46,6 +46,9 @@ __attribute__((always_inline)) INLINE static void stars_init_spart(
   sp->num_ngb_density = 0;
 #endif
 
+  if(sp->id == 155966626889L)
+    message("Initialising spart id=%lld", sp->id);
+
   sp->density.wcount = 0.f;
   sp->density.wcount_dh = 0.f;
 }