diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index 6c9087664fb1feba10b4ccf4e46a3e0b1e143dc1..eca6667ca4a9d218521af32a228f125b89978b90 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -52,6 +52,7 @@
 #include "proxy.h"
 #include "task_order.h"
 #include "timers.h"
+#include "star_formation.h"
 
 extern int engine_max_parts_per_ghost;
 extern int engine_max_sparts_per_ghost;
@@ -916,7 +917,7 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
 
       /* Particle recoupling */
       // ALEXEI: make sure we're only doing this for the types of star formation which actually decouple particles
-      if (with_star_formation) {
+      if (with_star_formation && swift_star_formation_model_decouples_particles) {
 	    c->part_recouple = scheduler_addtask(s, task_type_part_recouple,
 	    				     task_subtype_none, 0, 0, c, NULL);
 
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 101e1c095ad58d366ea2350d4e3ed03654ef5124..7417d52729c382813aa17b4b7aa26bfe364c0ba4 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -804,7 +804,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
     struct part *restrict p, const struct cosmology *cosmo) {
 
   p->force.h_dt *= p->h * hydro_dimension_inv;
-  if (p->id == SIMBA_DEBUG_ID) message("p->force.h_dt %.5e p->h %.5e hydro_dimension_inv %.5e", p->force.h_dt, p->h, hydro_dimension_inv);
 
   p->entropy_dt =
       0.5f * gas_entropy_from_internal_energy(p->rho, p->entropy_dt);
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 84fb07b97c2530c1f7ae9eab9755b3e68880fb60..d2996b827fd4c5f3339e3236bda3962fc2c2ac06 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -667,7 +667,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Get the time derivative for h. */
   pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
-  if (pi->id == SIMBA_DEBUG_ID && !isfinite(pi->force.h_dt)) message("pi id %llu pj id %llu pi->force.h_dt %e mi %e dvdr %e r_inv %e rhoi %e wj_dr %e %e", pi->id, pj->id, pi->force.h_dt, mj, dvdr, r_inv, rhoi, wj_dr,  mj * dvdr * r_inv / rhoi * wj_dr);
 
   /* Update the signal velocity. */
   pi->force.v_sig = max(pi->force.v_sig, v_sig);
diff --git a/src/star_formation.h b/src/star_formation.h
index a527db6903eaf4749f0bd404cc70e2ddf5a36e1d..e29dbe30024c7bc3f5fb53248c51eb7ecbca2923 100644
--- a/src/star_formation.h
+++ b/src/star_formation.h
@@ -33,18 +33,23 @@
 /* Import the right star formation law definition */
 #if defined(STAR_FORMATION_NONE)
 #define swift_star_formation_model_creates_stars 1
+#define swift_star_formation_model_decouples_particles 0
 #include "./star_formation/none/star_formation.h"
 #elif defined(STAR_FORMATION_QLA)
 #define swift_star_formation_model_creates_stars 0
+#define swift_star_formation_model_decouples_particles 0
 #include "./star_formation/QLA/star_formation.h"
 #elif defined(STAR_FORMATION_EAGLE)
 #define swift_star_formation_model_creates_stars 1
+#define swift_star_formation_model_decouples_particles 0
 #include "./star_formation/EAGLE/star_formation.h"
 #elif defined(STAR_FORMATION_GEAR)
 #define swift_star_formation_model_creates_stars 1
+#define swift_star_formation_model_decouples_particles 0
 #include "./star_formation/GEAR/star_formation.h"
 #elif defined(STAR_FORMATION_SIMBA)
 #define swift_star_formation_model_creates_stars 1
+#define swift_star_formation_model_decouples_particles 1
 #include "./star_formation/SIMBA/star_formation.h"
 #else
 #error "Invalid choice of star formation law"
diff --git a/src/timestep_limiter_iact.h b/src/timestep_limiter_iact.h
index 394654ad3057a44584d66ded7a244913bc4dfa65..8c2e846f36ed2fad2baa5c843774f7d2134b11ef 100644
--- a/src/timestep_limiter_iact.h
+++ b/src/timestep_limiter_iact.h
@@ -36,8 +36,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_timebin(
     struct part *restrict pi, struct part *restrict pj, const float a,
     const float H) {
 
-  //if (pi->id == SIMBA_DEBUG_ID) message("pi id %llu time_bin %d ngb_time_bin %d", pi->id, pi->time_bin, pi->limiter_data.min_ngb_time_bin);
-  //if (pj->id == SIMBA_DEBUG_ID) message("pj id %llu time_bin %d ngb_time_bin %d", pj->id, pj->time_bin, pj->limiter_data.min_ngb_time_bin);
   /* Update the minimal time-bin */
   if (pj->time_bin > 0)
     pi->limiter_data.min_ngb_time_bin =
@@ -65,8 +63,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_timebin(
     struct part *restrict pi, const struct part *restrict pj, const float a,
     const float H) {
 
-  //if (pi->id == SIMBA_DEBUG_ID) message("pi id %llu time_bin %d ngb_time_bin %d pj id %llu time_bin %d ngb_time_bin %d", pi->id, pi->time_bin, pi->limiter_data.min_ngb_time_bin, pj->id, pj->time_bin, pj->limiter_data.min_ngb_time_bin);
-  //if (pj->id == SIMBA_DEBUG_ID) message("pj id %llu time_bin %d ngb_time_bin %d pi id %llu time_bin %d ngb_time_bin %d", pj->id, pj->time_bin, pj->limiter_data.min_ngb_time_bin, pi->id, pi->time_bin, pi->limiter_data.min_ngb_time_bin);
   /* Update the minimal time-bin */
   if (pj->time_bin > 0)
     pi->limiter_data.min_ngb_time_bin =