diff --git a/src/hydro/GizmoMFM/hydro.h b/src/hydro/GizmoMFM/hydro.h
index dad33486edfb3c06063627db01f6f33993e65e58..b6f6dc2b5ee383bb23787aeaac0c6959b3ff278e 100644
--- a/src/hydro/GizmoMFM/hydro.h
+++ b/src/hydro/GizmoMFM/hydro.h
@@ -339,7 +339,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
      this is why we divide by the volume, and not by the density */
   p->P = hydro_gamma_minus_one * energy * volume_inv;
   p->A = p->conserved.entropy / m;
-  p->P = p->A * pow_gamma(p->rho);
 #endif
 
   /* sanity checks */
@@ -475,6 +474,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   p->flux.momentum[1] = 0.0f;
   p->flux.momentum[2] = 0.0f;
   p->flux.energy = 0.0f;
+
+  p->force.Ekinmax = 0.0f;
 }
 
 /**
@@ -654,7 +655,6 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
       p->conserved.mass * gas_internal_energy_from_entropy(0.0f, 0.0f);
 #else
   p->conserved.energy += p->flux.energy * dt_therm;
-  p->conserved.energy = hydro_one_over_gamma_minus_one * p->conserved.entropy * pow_gamma_minus_one(p->rho);
 #endif
 
 #ifndef HYDRO_GAMMA_5_3
@@ -704,6 +704,18 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     p->conserved.momentum[0] += dt_grav * p->conserved.mass * a_grav[0];
     p->conserved.momentum[1] += dt_grav * p->conserved.mass * a_grav[1];
     p->conserved.momentum[2] += dt_grav * p->conserved.mass * a_grav[2];
+
+    /* apply the entropy switch */
+    const float dEgrav = p->conserved.mass *
+                         sqrtf(a_grav[0] * a_grav[0] + a_grav[1] * a_grav[1] +
+                               a_grav[2] * a_grav[2]) *
+                         p->h;
+    if (p->conserved.energy <
+            0.001 * (p->force.Ekinmax + p->conserved.energy) ||
+        p->conserved.energy < 0.001 * dEgrav) {
+      p->conserved.energy = hydro_one_over_gamma_minus_one *
+                            p->conserved.entropy * pow_gamma_minus_one(p->rho);
+    }
   }
 
   /* Set the velocities: */
diff --git a/src/hydro/GizmoMFM/hydro_iact.h b/src/hydro/GizmoMFM/hydro_iact.h
index 08e3de538a5e2cd775a69abab9bcb11fb09e7a9b..f0a3b39b89fa3a8ee55d9c26a2cba115db827e23 100644
--- a/src/hydro/GizmoMFM/hydro_iact.h
+++ b/src/hydro/GizmoMFM/hydro_iact.h
@@ -260,6 +260,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   } else
     vmax = 0.0f;
 
+  /* calculate flux limiter quantities */
+  const float Wjminvi[3] = {Wj[1] - vi[0], Wj[2] - vi[1], Wj[3] - vi[2]};
+  const float Ekinngbj = 0.5f * pj->conserved.mass *
+                         (Wjminvi[0] * Wjminvi[0] + Wjminvi[1] * Wjminvi[1] +
+                          Wjminvi[2] * Wjminvi[2]);
+  pi->force.Ekinmax = max(pi->force.Ekinmax, Ekinngbj);
+  if (mode == 1) {
+    const float Wiminvj[3] = {Wi[1] - vj[0], Wi[2] - vj[1], Wi[3] - vj[2]};
+    const float Ekinngbi = 0.5f * pi->conserved.mass *
+                           (Wiminvj[0] * Wiminvj[0] + Wiminvj[1] * Wiminvj[1] +
+                            Wiminvj[2] * Wiminvj[2]);
+    pj->force.Ekinmax = max(pj->force.Ekinmax, Ekinngbi);
+  }
+
   /* Velocity on the axis linking the particles */
   float dvdr = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
                (Wi[3] - Wj[3]) * dx[2];
diff --git a/src/hydro/GizmoMFM/hydro_part.h b/src/hydro/GizmoMFM/hydro_part.h
index 0adfed995252e0c0bf8223bf25fd032ded8d377a..17bc03836c99530c2734efac86875a281d54c2c8 100644
--- a/src/hydro/GizmoMFM/hydro_part.h
+++ b/src/hydro/GizmoMFM/hydro_part.h
@@ -128,10 +128,8 @@ struct part {
 
         /* Maximum signal velocity among all the neighbours of the particle. The
          * signal velocity encodes information about the relative fluid
-         * velocities
-         * AND particle velocities of the neighbour and this particle, as well
-         * as
-         * the sound speed of both particles. */
+         * velocities AND particle velocities of the neighbour and this
+         * particle, as well as the sound speed of both particles. */
         float vmax;
 
       } timestepvars;
@@ -142,6 +140,10 @@ struct part {
         /* Needed to drift the primitive variables. */
         float h_dt;
 
+        /* Maximum kinetic energy of all neighbouring particles (in the rest
+         * frame of this particle). Used for the entropy switch. */
+        float Ekinmax;
+
       } force;
     };
   };