diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index 3fd357a2d8778f5ca8b014935d538350eccb99c6..bfb5cd1ce39a9908573c66406f41b56561a870d6 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -33,7 +33,7 @@
  * @param dt Time since the last kick
  */
 __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
-    const struct part *restrict p, float dt) {
+    const struct part *restrict p) {
 
   return p->u;
 }
@@ -45,7 +45,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
  * @param dt Time since the last kick
  */
 __attribute__((always_inline)) INLINE static float hydro_get_pressure(
-    const struct part *restrict p, float dt) {
+    const struct part *restrict p) {
 
   return gas_pressure_from_internal_energy(p->rho, p->u);
 }
@@ -57,7 +57,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_pressure(
  * @param dt Time since the last kick
  */
 __attribute__((always_inline)) INLINE static float hydro_get_entropy(
-    const struct part *restrict p, float dt) {
+    const struct part *restrict p) {
 
   return gas_entropy_from_internal_energy(p->rho, p->u);
 }
@@ -69,7 +69,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
  * @param dt Time since the last kick
  */
 __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
-    const struct part *restrict p, float dt) {
+    const struct part *restrict p) {
 
   return p->force.soundspeed;
 }
@@ -97,34 +97,30 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
 }
 
 /**
- * @brief Modifies the thermal state of a particle to the imposed internal
- * energy
+ * @brief Returns the time derivative of internal energy of a particle
  *
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * We assume a constant density.
  *
- * @param p The particle
- * @param u The new internal energy
+ * @param p The particle of interest
  */
-__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
-    struct part *restrict p, float u) {
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
+    const struct part *restrict p) {
 
-  p->u = u;
+  return p->force.u_dt;
 }
 
 /**
- * @brief Modifies the thermal state of a particle to the imposed entropy
+ * @brief Returns the time derivative of internal energy of a particle
  *
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * We assume a constant density.
  *
- * @param p The particle
- * @param S The new entropy
+ * @param p The particle of interest.
+ * @param du_dt The new time derivative of the internal energy.
  */
-__attribute__((always_inline)) INLINE static void hydro_set_entropy(
-    struct part *restrict p, float S) {
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
+    struct part *restrict p, float du_dt) {
 
-  p->u = gas_internal_energy_from_entropy(p->rho, S);
+  p->force.u_dt = du_dt;
 }
 
 /**
@@ -152,26 +148,6 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
   return min(dt_cfl, dt_u_change);
 }
 
-/**
- * @brief Initialises the particles for the first time
- *
- * This function is called only once just after the ICs have been
- * read in to do some conversions.
- *
- * @param p The particle to act upon
- * @param xp The extended particle data to act upon
- */
-__attribute__((always_inline)) INLINE static void hydro_first_init_part(
-    struct part *restrict p, struct xpart *restrict xp) {
-
-  p->ti_begin = 0;
-  p->ti_end = 0;
-  xp->v_full[0] = p->v[0];
-  xp->v_full[1] = p->v[1];
-  xp->v_full[2] = p->v[2];
-  xp->u_full = p->u;
-}
-
 /**
  * @brief Prepares a particle for the density calculation.
  *
@@ -244,8 +220,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
  * @param time The current time
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
-    struct part *restrict p, struct xpart *restrict xp, int ti_current,
-    double timeBase) {
+    struct part *restrict p, struct xpart *restrict xp) {
 
   /* Some smoothing length multiples. */
   const float h = p->h;
@@ -270,17 +245,18 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   p->force.balsara = normDiv_v / (normDiv_v + normRot_v + 0.0001f * fc * h_inv);
 
   /* Viscosity parameter decay time */
-  const float tau = h / (2.f * const_viscosity_length * p->force.soundspeed);
+  /* const float tau = h / (2.f * const_viscosity_length * p->force.soundspeed);
+   */
 
   /* Viscosity source term */
-  const float S = max(-normDiv_v, 0.f);
+  /* const float S = max(-normDiv_v, 0.f); */
 
   /* Compute the particle's viscosity parameter time derivative */
-  const float alpha_dot = (const_viscosity_alpha_min - p->alpha) / tau +
-                          (const_viscosity_alpha_max - p->alpha) * S;
+  /* const float alpha_dot = (const_viscosity_alpha_min - p->alpha) / tau + */
+  /*                         (const_viscosity_alpha_max - p->alpha) * S; */
 
   /* Update particle's viscosity paramter */
-  p->alpha += alpha_dot * (p->ti_end - p->ti_begin) * timeBase;
+  /* p->alpha += alpha_dot * (p->ti_end - p->ti_begin) * timeBase; */  // MATTHIEU
 }
 
 /**
@@ -305,6 +281,22 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
   p->force.v_sig = 0.0f;
 }
 
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param p The particle.
+ * @param xp The extended data of this particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
+    struct part *restrict p, const struct xpart *restrict xp) {
+
+  /* Re-set the predicted velocities */
+  p->v[0] = xp->v_full[0];
+  p->v[1] = xp->v_full[1];
+  p->v[2] = xp->v_full[2];
+}
+
 /**
  * @brief Predict additional particle fields forward in time when drifting
  *
@@ -316,8 +308,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
  * @param timeBase The minimal time-step size
  */
 __attribute__((always_inline)) INLINE static void hydro_predict_extra(
-    struct part *restrict p, struct xpart *restrict xp, float dt, int t0,
-    int t1, double timeBase) {
+    struct part *restrict p, struct xpart *restrict xp, float dt) {
   float u, w;
 
   const float h_inv = 1.f / p->h;
@@ -368,8 +359,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
  * @param half_dt The half time-step for this kick
  */
 __attribute__((always_inline)) INLINE static void hydro_kick_extra(
-    struct part *restrict p, struct xpart *restrict xp, float dt,
-    float half_dt) {}
+    struct part *restrict p, struct xpart *restrict xp, float dt) {}
 
 /**
  *  @brief Converts hydro quantity of a particle at the start of a run
@@ -379,6 +369,28 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
  * @param p The particle to act upon
  */
 __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
-    struct part *restrict p) {}
+    struct part *restrict p, struct xpart *restrict xp) {}
+
+/**
+ * @brief Initialises the particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_first_init_part(
+    struct part *restrict p, struct xpart *restrict xp) {
+
+  p->time_bin = 0;
+  xp->v_full[0] = p->v[0];
+  xp->v_full[1] = p->v[1];
+  xp->v_full[2] = p->v[2];
+  xp->u_full = p->u;
+
+  hydro_reset_acceleration(p);
+  hydro_init_part(p);
+}
 
 #endif /* SWIFT_DEFAULT_HYDRO_H */
diff --git a/src/hydro/Default/hydro_debug.h b/src/hydro/Default/hydro_debug.h
index d02d3ef82c1b3d751731f49850c06df4b146b164..3be9c9e1760591423edbd218d19b46ddf9aad01e 100644
--- a/src/hydro/Default/hydro_debug.h
+++ b/src/hydro/Default/hydro_debug.h
@@ -25,11 +25,10 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       "x=[%.3e,%.3e,%.3e], "
       "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
       "h=%.3e, "
-      "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, t_begin=%d, t_end=%d\n",
+      "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, time_bin=%d\n",
       p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
       xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
-      p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->ti_begin,
-      p->ti_end);
+      p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->time_bin);
 }
 
 #endif /* SWIFT_DEFAULT_HYDRO_DEBUG_H */
diff --git a/src/hydro/Default/hydro_part.h b/src/hydro/Default/hydro_part.h
index c7464bcf338b1c5b81ffa91d92264c2bd35e9313..332eecb27fb65a6b4da48cbb595450a432c44615 100644
--- a/src/hydro/Default/hydro_part.h
+++ b/src/hydro/Default/hydro_part.h
@@ -55,12 +55,6 @@ struct part {
   /* Particle cutoff radius. */
   float h;
 
-  /* Particle time of beginning of time-step. */
-  int ti_begin;
-
-  /* Particle time of end of time-step. */
-  int ti_end;
-
   /* Particle internal energy. */
   float u;
 
@@ -125,6 +119,9 @@ struct part {
   /* Pointer to corresponding gravity part. */
   struct gpart* gpart;
 
+  /* Particle time-bin */
+  timebin_t time_bin;
+
 } SWIFT_STRUCT_ALIGN;
 
 #endif /* SWIFT_DEFAULT_HYDRO_PART_H */
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 66c9203e39e56d520eeace8858b0c618b45e6a22..ce46d8b2dde76e28be468ac89ba0e0bf74edf786 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -645,13 +645,13 @@ void write_output_parallel(struct engine* e, const char* baseName,
   size_t N[NUM_PARTICLE_TYPES] = {Ngas, Ndm, 0};
   long long N_total[NUM_PARTICLE_TYPES] = {0};
   long long offset[NUM_PARTICLE_TYPES] = {0};
-  MPI_Exscan(&N, &offset, NUM_PARTICLE_TYPES, MPI_LONG_LONG, MPI_SUM, comm);
+  MPI_Exscan(&N, &offset, NUM_PARTICLE_TYPES, MPI_LONG_LONG_INT, MPI_SUM, comm);
   for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype)
     N_total[ptype] = offset[ptype] + N[ptype];
 
   /* The last rank now has the correct N_total. Let's
    * broadcast from there */
-  MPI_Bcast(&N_total, 6, MPI_LONG_LONG, mpi_size - 1, comm);
+  MPI_Bcast(&N_total, 6, MPI_LONG_LONG_INT, mpi_size - 1, comm);
 
   /* Now everybody konws its offset and the total number of
    * particles of each
diff --git a/src/serial_io.c b/src/serial_io.c
index b9ad0fbaa856a889d3f84bb42013282f3640fd5e..49dd945006d5907c43d58babc28f72914484a742 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -536,7 +536,7 @@ void read_ic_serial(char* fileName, const struct UnitSystem* internal_units,
   /* Now need to broadcast that information to all ranks. */
   MPI_Bcast(flag_entropy, 1, MPI_INT, 0, comm);
   MPI_Bcast(periodic, 1, MPI_INT, 0, comm);
-  MPI_Bcast(&N_total, NUM_PARTICLE_TYPES, MPI_LONG_LONG, 0, comm);
+  MPI_Bcast(&N_total, NUM_PARTICLE_TYPES, MPI_LONG_LONG_INT, 0, comm);
   MPI_Bcast(dim, 3, MPI_DOUBLE, 0, comm);
   MPI_Bcast(ic_units, sizeof(struct UnitSystem), MPI_BYTE, 0, comm);
 
@@ -694,12 +694,12 @@ void write_output_serial(struct engine* e, const char* baseName,
   size_t N[NUM_PARTICLE_TYPES] = {Ngas, Ndm, 0};
   long long N_total[NUM_PARTICLE_TYPES] = {0};
   long long offset[NUM_PARTICLE_TYPES] = {0};
-  MPI_Exscan(&N, &offset, NUM_PARTICLE_TYPES, MPI_LONG_LONG, MPI_SUM, comm);
+  MPI_Exscan(&N, &offset, NUM_PARTICLE_TYPES, MPI_LONG_LONG_INT, MPI_SUM, comm);
   for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype)
     N_total[ptype] = offset[ptype] + N[ptype];
 
   /* The last rank now has the correct N_total. Let's broadcast from there */
-  MPI_Bcast(&N_total, 6, MPI_LONG_LONG, mpi_size - 1, comm);
+  MPI_Bcast(&N_total, 6, MPI_LONG_LONG_INT, mpi_size - 1, comm);
 
   /* Now everybody konws its offset and the total number of particles of each
    * type */