diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index 890380b23ac1496854dc5359f9249993c166d2c3..1228aa69abfe812e5b2e73f066bb13be3292aa20 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -21,71 +21,134 @@
 
 #include "adiabatic_index.h"
 #include "approx_math.h"
+#include "cosmology.h"
+#include "dimension.h"
 #include "equation_of_state.h"
 #include "hydro_space.h"
+#include "kernel_hydro.h"
 #include "minmax.h"
 
 #include <float.h>
 
 /**
- * @brief Returns the internal energy of a particle
+ * @brief Returns the comoving internal energy of a particle
  *
  * @param p The particle of interest
- * @param dt Time since the last kick
  */
-__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
-    const struct part *restrict p) {
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_internal_energy(const struct part *restrict p) {
 
   return p->u;
 }
 
 /**
- * @brief Returns the pressure of a particle
+ * @brief Returns the physical internal energy of a particle
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_internal_energy(const struct part *restrict p,
+                                   const struct cosmology *cosmo) {
+
+  return cosmo->a_factor_internal_energy * p->u;
+}
+
+/**
+ * @brief Returns the comoving pressure of a particle
  *
  * @param p The particle of interest
- * @param dt Time since the last kick
  */
-__attribute__((always_inline)) INLINE static float hydro_get_pressure(
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
     const struct part *restrict p) {
 
   return gas_pressure_from_internal_energy(p->rho, p->u);
 }
 
 /**
- * @brief Returns the entropy of a particle
+ * @brief Returns the physical pressure of a particle
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
+    const struct part *restrict p, const struct cosmology *cosmo) {
+
+  return cosmo->a_factor_pressure *
+         gas_pressure_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the comoving entropy of a particle
  *
  * @param p The particle of interest
- * @param dt Time since the last kick
  */
-__attribute__((always_inline)) INLINE static float hydro_get_entropy(
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
     const struct part *restrict p) {
 
   return gas_entropy_from_internal_energy(p->rho, p->u);
 }
 
 /**
- * @brief Returns the sound speed of a particle
+ * @brief Returns the physical entropy of a particle
  *
  * @param p The particle of interest
- * @param dt Time since the last kick
+ * @param cosmo The cosmological model.
  */
-__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
-    const struct part *restrict p) {
+__attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
+    const struct part *restrict p, const struct cosmology *cosmo) {
+
+  /* Note: no cosmological conversion required here with our choice of
+   * coordinates. */
+  return gas_entropy_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the comoving sound speed of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_soundspeed(const struct part *restrict p) {
 
   return p->force.soundspeed;
 }
 
 /**
- * @brief Returns the density of a particle
+ * @brief Returns the physical sound speed of a particle
  *
  * @param p The particle of interest
+ * @param cosmo The cosmological model.
  */
-__attribute__((always_inline)) INLINE static float hydro_get_density(
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_soundspeed(const struct part *restrict p,
+                              const struct cosmology *cosmo) {
+
+  return cosmo->a_factor_sound_speed * p->force.soundspeed;
+}
+
+/**
+ * @brief Returns the comoving density of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
     const struct part *restrict p) {
 
   return p->rho;
 }
 
+/**
+ * @brief Returns the physical density of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_density(
+    const struct part *restrict p, const struct cosmology *cosmo) {
+
+  return p->rho * cosmo->a3_inv;
+}
+
 /**
  * @brief Returns the mass of a particle
  *
@@ -102,16 +165,20 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
  *
  * @param p The particle of interest
  * @param xp The extended data of the particle.
- * @param dt The time since the last kick.
+ * @param dt_kick_hydro The time (for hydro accelerations) since the last kick.
+ * @param dt_kick_grav The time (for gravity accelerations) since the last kick.
  * @param v (return) The velocities at the current time.
  */
 __attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
-    const struct part *restrict p, const struct xpart *xp, float dt,
-    float v[3]) {
-
-  v[0] = xp->v_full[0] + p->a_hydro[0] * dt;
-  v[1] = xp->v_full[1] + p->a_hydro[1] * dt;
-  v[2] = xp->v_full[2] + p->a_hydro[2] * dt;
+    const struct part *restrict p, const struct xpart *xp, float dt_kick_hydro,
+    float dt_kick_grav, float v[3]) {
+
+  v[0] = xp->v_full[0] + p->a_hydro[0] * dt_kick_hydro +
+         xp->a_grav[0] * dt_kick_grav;
+  v[1] = xp->v_full[1] + p->a_hydro[1] * dt_kick_hydro +
+         xp->a_grav[1] * dt_kick_grav;
+  v[2] = xp->v_full[2] + p->a_hydro[2] * dt_kick_hydro +
+         xp->a_grav[2] * dt_kick_grav;
 }
 
 /**
@@ -146,17 +213,19 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
  *
  * @param p Pointer to the particle data
  * @param xp Pointer to the extended particle data
- *
+ * @param hydro_properties The constants used in the scheme
+ * @param cosmo The cosmological model.
  */
 __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
     const struct part *restrict p, const struct xpart *restrict xp,
-    const struct hydro_props *restrict hydro_properties) {
+    const struct hydro_props *restrict hydro_properties,
+    const struct cosmology *restrict cosmo) {
 
-  const float CFL_condition = hydro_properties->CFL_condition;
+  const float CFL = hydro_properties->CFL_condition;
 
   /* CFL condition */
-  const float dt_cfl =
-      2.f * kernel_gamma * CFL_condition * p->h / p->force.v_sig;
+  const float dt_cfl = 2.f * kernel_gamma * CFL * cosmo->a * p->h /
+                       (cosmo->a_factor_sound_speed * p->force.v_sig);
 
   /* Limit change in u */
   const float dt_u_change =
@@ -187,6 +256,7 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
  */
 __attribute__((always_inline)) INLINE static void hydro_init_part(
     struct part *restrict p, const struct hydro_space *hs) {
+
   p->density.wcount = 0.f;
   p->density.wcount_dh = 0.f;
   p->rho = 0.f;
@@ -204,9 +274,10 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
  * and add the self-contribution term.
  *
  * @param p The particle to act upon
+ * @param cosmo The current cosmological model.
  */
 __attribute__((always_inline)) INLINE static void hydro_end_density(
-    struct part *restrict p) {
+    struct part *restrict p, const struct cosmology *cosmo) {
 
   /* Some smoothing length multiples. */
   const float h = p->h;
@@ -226,14 +297,15 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
   p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
 
   const float irho = 1.f / p->rho;
+  const float a_inv2 = cosmo->a2_inv;
 
   /* Finish calculation of the velocity curl components */
-  p->density.rot_v[0] *= h_inv_dim_plus_one * irho;
-  p->density.rot_v[1] *= h_inv_dim_plus_one * irho;
-  p->density.rot_v[2] *= h_inv_dim_plus_one * irho;
+  p->density.rot_v[0] *= h_inv_dim_plus_one * a_inv2 * irho;
+  p->density.rot_v[1] *= h_inv_dim_plus_one * a_inv2 * irho;
+  p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * irho;
 
   /* Finish calculation of the velocity divergence */
-  p->density.div_v *= h_inv_dim_plus_one * irho;
+  p->density.div_v *= h_inv_dim_plus_one * a_inv2 * irho;
 }
 
 /**
@@ -241,9 +313,11 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
  *
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
+ * @param cosmo The current cosmological model.
  */
 __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
-    struct part *restrict p, struct xpart *restrict xp) {
+    struct part *restrict p, struct xpart *restrict xp,
+    const struct cosmology *cosmo) {
 
   /* Some smoothing length multiples. */
   const float h = p->h;
@@ -268,10 +342,13 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
  *
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
- * @param time The current time
+ * @param cosmo The current cosmological model.
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
-    struct part *restrict p, struct xpart *restrict xp) {
+    struct part *restrict p, struct xpart *restrict xp,
+    const struct cosmology *cosmo) {
+
+  const float fac_mu = cosmo->a_factor_mu;
 
   /* Some smoothing length multiples. */
   const float h = p->h;
@@ -296,7 +373,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   p->force.P_over_rho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega);
 
   /* Balsara switch */
-  p->force.balsara = normDiv_v / (normDiv_v + normRot_v + 0.0001f * fc * h_inv);
+  p->force.balsara =
+      normDiv_v / (normDiv_v + normRot_v + 0.0001f * fac_mu * fc * h_inv);
 
   /* Viscosity parameter decay time */
   /* const float tau = h / (2.f * const_viscosity_length * p->force.soundspeed);
@@ -349,6 +427,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
   p->v[0] = xp->v_full[0];
   p->v[1] = xp->v_full[1];
   p->v[2] = xp->v_full[2];
+
+  p->u = xp->u_full;
 }
 
 /**
@@ -356,19 +436,18 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
  *
  * @param p The particle
  * @param xp The extended data of the particle
- * @param dt The drift time-step.
- * @param t0 The time at the start of the drift
- * @param t1 The time at the end of the drift
- * @param timeBase The minimal time-step size
+ * @param dt_drift The drift time-step for positions.
+ * @param dt_therm The drift time-step for thermal quantities.
  */
 __attribute__((always_inline)) INLINE static void hydro_predict_extra(
-    struct part *restrict p, struct xpart *restrict xp, float dt) {
+    struct part *restrict p, struct xpart *restrict xp, float dt_drift,
+    float dt_therm) {
   float u, w;
 
   const float h_inv = 1.f / p->h;
 
   /* Predict smoothing length */
-  const float w1 = p->force.h_dt * h_inv * dt;
+  const float w1 = p->force.h_dt * h_inv * dt_drift;
   if (fabsf(w1) < 0.2f)
     p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
   else
@@ -382,7 +461,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
     p->rho *= expf(w2);
 
   /* Predict internal energy */
-  w = p->force.u_dt / p->u * dt;
+  w = p->force.u_dt / p->u * dt_therm;
   if (fabsf(w) < 0.2f)
     u = p->u *= approx_expf(w);
   else
@@ -398,9 +477,10 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
  * Multiplies the forces and accelerationsby the appropiate constants
  *
  * @param p The particle to act upon
+ * @param cosmo The current cosmological model.
  */
 __attribute__((always_inline)) INLINE static void hydro_end_force(
-    struct part *restrict p) {
+    struct part *restrict p, const struct cosmology *cosmo) {
   p->force.h_dt *= p->h * hydro_dimension_inv;
 }
 
@@ -409,11 +489,10 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
  *
  * @param p The particle to act upon
  * @param xp The particle extended data to act upon
- * @param dt The time-step for this kick
- * @param half_dt The half time-step for this kick
+ * @param dt_therm The time-step for this kick (for thermodynamic quantities)
  */
 __attribute__((always_inline)) INLINE static void hydro_kick_extra(
-    struct part *restrict p, struct xpart *restrict xp, float dt) {}
+    struct part *restrict p, struct xpart *restrict xp, float dt_therm) {}
 
 /**
  *  @brief Converts hydro quantity of a particle at the start of a run
@@ -441,6 +520,9 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
   xp->v_full[0] = p->v[0];
   xp->v_full[1] = p->v[1];
   xp->v_full[2] = p->v[2];
+  xp->a_grav[0] = 0.f;
+  xp->a_grav[1] = 0.f;
+  xp->a_grav[2] = 0.f;
   xp->u_full = p->u;
 
   hydro_reset_acceleration(p);
diff --git a/src/hydro/Default/hydro_iact.h b/src/hydro/Default/hydro_iact.h
index 89b90be22b831066c8b9211bfe50208b4c4af67e..658b4aba83085610a49bb9d2579d4f20c70d6c5b 100644
--- a/src/hydro/Default/hydro_iact.h
+++ b/src/hydro/Default/hydro_iact.h
@@ -26,22 +26,32 @@
  * @brief SPH interaction functions following the Gadget-2 version of SPH.
  *
  * The interactions computed here are the ones presented in the Gadget-2 paper
- *and use the same
+ * and use the same
  * numerical coefficients as the Gadget-2 code. When used with the Spline-3
- *kernel, the results
+ * kernel, the results
  * should be equivalent to the ones obtained with Gadget-2 up to the rounding
- *errors and interactions
+ * errors and interactions
  * missed by the Gadget-2 tree-code neighbours search.
  *
  * The code uses internal energy instead of entropy as a thermodynamical
- *variable.
+ * variable.
  */
 
 /**
- * @brief Density loop
+ * @brief Density interaction between two particles.
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_density(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H) {
 
   float r = sqrtf(r2), ri = 1.0f / r;
   float xi, xj;
@@ -97,10 +107,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
 }
 
 /**
- * @brief Density loop (non-symmetric version)
+ * @brief Density interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    const struct part *restrict pj, float a, float H) {
 
   float r, ri;
   float xi;
@@ -145,10 +165,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
 }
 
 /**
- * @brief Force loop
+ * @brief Force interaction between two particles.
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_force(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H) {
 
   float r = sqrtf(r2), ri = 1.0f / r;
   float xi, xj;
@@ -160,6 +190,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   float f;
   int k;
 
+  /* Cosmological factors entering the EoMs */
+  const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+  const float a2_Hubble = a * a * H;
+
   /* Get some values in local variables. */
   mi = pi->mass;
   mj = pj->mass;
@@ -184,12 +218,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
 
   /* Compute dv dot r. */
   dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
-         (pi->v[2] - pj->v[2]) * dx[2];
+         (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
   dvdr *= ri;
 
   /* Compute the relative velocity. (This is 0 if the particles move away from
    * each other and negative otherwise) */
-  omega_ij = min(dvdr, 0.f);
+  omega_ij = min(fac_mu * dvdr, 0.f);
 
   /* Compute signal velocity */
   v_sig = pi->force.soundspeed + pj->force.soundspeed - 2.0f * omega_ij;
@@ -240,10 +274,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
 }
 
 /**
- * @brief Force loop (non-symmetric version)
+ * @brief Force interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    const struct part *restrict pj, float a, float H) {
 
   float r = sqrtf(r2), ri = 1.0f / r;
   float xi, xj;
@@ -255,6 +299,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   float f;
   int k;
 
+  /* Cosmological factors entering the EoMs */
+  const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+  const float a2_Hubble = a * a * H;
+
   /* Get some values in local variables. */
   // mi = pi->mass;
   mj = pj->mass;
@@ -279,12 +327,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Compute dv dot r. */
   dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
-         (pi->v[2] - pj->v[2]) * dx[2];
+         (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
   dvdr *= ri;
 
   /* Compute the relative velocity. (This is 0 if the particles move away from
    * each other and negative otherwise) */
-  omega_ij = min(dvdr, 0.f);
+  omega_ij = min(fac_mu * dvdr, 0.f);
 
   /* Compute signal velocity */
   v_sig = pi->force.soundspeed + pj->force.soundspeed - 2.0f * omega_ij;
diff --git a/src/hydro/Default/hydro_part.h b/src/hydro/Default/hydro_part.h
index 6290ee5ba37a567b3819a06b8b824d83071ee4c2..2a18e03cb533ca860f227a31152ef2058e0dd37d 100644
--- a/src/hydro/Default/hydro_part.h
+++ b/src/hydro/Default/hydro_part.h
@@ -34,6 +34,9 @@ struct xpart {
   /* Velocity at the last full step. */
   float v_full[3];
 
+  /* Gravitational acceleration at the last full step. */
+  float a_grav[3];
+
   /* Additional data used to record cooling information */
   struct cooling_xpart_data cooling_data;