diff --git a/examples/Cooling/CoolingBox/coolingBox.yml b/examples/Cooling/CoolingBox/coolingBox.yml
index ca97793aa9b4de43872db38f55c4060a2b27ddc9..f776823d1ba2b4703288913d4bf67f3e4851bebd 100644
--- a/examples/Cooling/CoolingBox/coolingBox.yml
+++ b/examples/Cooling/CoolingBox/coolingBox.yml
@@ -58,7 +58,7 @@ EAGLECooling:
   dir_name:                ./coolingtables/
   H_reion_z:               11.5
   H_reion_eV_p_H:          2.0
-  He_reion_z_center:       3.5
+  He_reion_z_centre:       3.5
   He_reion_z_sigma:        0.5
   He_reion_eV_p_H:         2.0
 
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 425190bc1531f8fa4822a705d6ee5d8656e04586..11426ac89970917cf121351e266aaa67cfd99847 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -41,7 +41,7 @@ check_PROGRAMS = testGreetings testReading testTimeIntegration \
 		 testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
 		 testGravityDerivatives testPotentialSelf testPotentialPair testEOS testUtilities \
 		 testSelectOutput testCbrt testCosmology testOutputList test27cellsStars \
-		 test27cellsStars_subset testCooling testFeedback testHashmap testHydroMPIrules
+		 test27cellsStars_subset testCooling testComovingCooling testFeedback testHashmap testHydroMPIrules
 
 # Rebuild tests when SWIFT is updated.
 $(check_PROGRAMS): ../src/.libs/libswiftsim.a
@@ -130,6 +130,8 @@ testUtilities_SOURCES = testUtilities.c
 
 testCooling_SOURCES = testCooling.c
 
+testComovingCooling_SOURCES = testComovingCooling.c
+
 testFeedback_SOURCES = testFeedback.c
 
 testHashmap_SOURCES = testHashmap.c
diff --git a/tests/testComovingCooling.c b/tests/testComovingCooling.c
new file mode 100644
index 0000000000000000000000000000000000000000..0b3b405c562b13c656cb596b25148bfad0935a5f
--- /dev/null
+++ b/tests/testComovingCooling.c
@@ -0,0 +1,237 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#include "../config.h"
+
+/* Local headers. */
+#include "../src/cooling/EAGLE/cooling_rates.h"
+#include "../src/cooling/EAGLE/cooling_tables.h"
+#include "../src/cooling/EAGLE/interpolate.h"
+#include "swift.h"
+
+/*
+ * @brief Assign particle density and entropy corresponding to the
+ * hydrogen number density and internal energy specified.
+ *
+ * @param p Particle data structure
+ * @param xp extra particle structure
+ * @param us unit system struct
+ * @param cooling Cooling function data structure
+ * @param cosmo Cosmology data structure
+ * @param phys_const Physical constants data structure
+ * @param nh_cgs Hydrogen number density (cgs units)
+ * @param u Internal energy (cgs units)
+ * @param ti_current integertime to set cosmo quantities
+ */
+void set_quantities(struct part *restrict p, struct xpart *restrict xp,
+                    const struct unit_system *restrict us,
+                    const struct cooling_function_data *restrict cooling,
+                    struct cosmology *restrict cosmo,
+                    const struct phys_const *restrict phys_const, float nh_cgs,
+                    double u_cgs, integertime_t ti_current) {
+  /* calculate density */
+  double hydrogen_number_density = nh_cgs / cooling->number_density_to_cgs;
+  p->rho = hydrogen_number_density * phys_const->const_proton_mass /
+           p->chemistry_data.metal_mass_fraction[chemistry_element_H];
+
+  /* update entropy based on internal energy */
+  float pressure = (u_cgs)*cooling->internal_energy_from_cgs * p->rho *
+                   (hydro_gamma_minus_one);
+  p->entropy = pressure * (pow(p->rho, -hydro_gamma));
+  xp->entropy_full = p->entropy;
+
+  p->entropy_dt = 0.f;
+}
+
+/*
+ * @brief Tests cooling integration scheme by comparing EAGLE
+ * integration to subcycled explicit equation.
+ */
+int main(int argc, char **argv) {
+  // Declare relevant structs
+  struct swift_params *params = malloc(sizeof(struct swift_params));
+  struct unit_system us;
+  struct chemistry_global_data chem_data;
+  struct part p;
+  struct xpart xp;
+  struct phys_const phys_const;
+  struct cooling_function_data cooling;
+  struct cosmology cosmo;
+  char *parametersFileName = "./testCooling.yml";
+
+  float nh_cgs;  // hydrogen number density
+  double u_cgs;  // internal energy
+
+  const float seconds_per_year = 3.154e7;
+
+  /* Number of values to test for in redshift,
+   * hydrogen number density and internal energy */
+  const int n_z = 10;
+  const int n_nh = 10;
+  const int n_u = 10;
+
+  /* Number of subcycles and tolerance used to compare
+   * subcycled and implicit solution. Note, high value
+   * of tolerance due to mismatch between explicit and
+   * implicit solution for large timesteps */
+  const float integration_tolerance = 0.2;
+
+  /* Read the parameter file */
+  if (params == NULL) error("Error allocating memory for the parameter file.");
+  message("Reading runtime parameters from file '%s'", parametersFileName);
+  parser_read_file(parametersFileName, params);
+
+  /* Init units */
+  units_init_from_params(&us, params, "InternalUnitSystem");
+  phys_const_init(&us, params, &phys_const);
+
+  /* Init chemistry */
+  chemistry_init(params, &us, &phys_const, &chem_data);
+  chemistry_first_init_part(&phys_const, &us, &cosmo, &chem_data, &p, &xp);
+  chemistry_part_has_no_neighbours(&p, &xp, &chem_data, &cosmo);
+  chemistry_print(&chem_data);
+
+  /* Init cosmology */
+  cosmology_init(params, &us, &phys_const, &cosmo);
+  cosmology_print(&cosmo);
+
+  /* Set dt */
+  const int timebin = 38;
+  float dt_cool, dt_therm;
+
+  /* Init hydro_props */
+  struct hydro_props hydro_properties;
+  hydro_props_init(&hydro_properties, &phys_const, &us, params);
+
+  /* Init cooling */
+  cooling_init(params, &us, &phys_const, &hydro_properties, &cooling);
+  cooling_print(&cooling);
+  cooling_update(&cosmo, &cooling, 0);
+
+  /* Init entropy floor */
+  struct entropy_floor_properties floor_props;
+  entropy_floor_init(&floor_props, &phys_const, &us, &hydro_properties, params);
+
+  /* Cooling function needs to know the minimal energy. Set it to the lowest
+   * internal energy in the cooling table. */
+  hydro_properties.minimal_internal_energy =
+      exp(M_LN10 * cooling.Therm[0]) * cooling.internal_energy_from_cgs;
+
+  /* Calculate abundance ratios */
+  float *abundance_ratio;
+  abundance_ratio = malloc((chemistry_element_count + 2) * sizeof(float));
+  abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
+
+  /* extract mass fractions, calculate table indices and offsets */
+  float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
+  float HeFrac =
+      p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
+      (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
+  int He_i;
+  float d_He;
+  get_index_1d(cooling.HeFrac, eagle_cooling_N_He_frac, HeFrac, &He_i, &d_He);
+
+  /* calculate spacing in nh and u */
+  const float log_u_min_cgs = 11, log_u_max_cgs = 17;
+  const float log_nh_min_cgs = -6, log_nh_max_cgs = 3;
+  const float delta_log_nh_cgs = (log_nh_max_cgs - log_nh_min_cgs) / n_nh;
+  const float delta_log_u_cgs = (log_u_max_cgs - log_u_min_cgs) / n_u;
+
+  /* Declare variables we will be checking */
+  double du_dt_implicit, du_dt_check;
+  integertime_t ti_current;
+
+  /* Loop over values of nh and u */
+  for (int nh_i = 0; nh_i < n_nh; nh_i++) {
+    nh_cgs = exp(M_LN10 * log_nh_min_cgs + delta_log_nh_cgs * nh_i);
+    for (int u_i = 0; u_i < n_u; u_i++) {
+      u_cgs = exp(M_LN10 * log_u_min_cgs + delta_log_u_cgs * u_i);
+
+      /* Calculate cooling solution at redshift zero if we're doing the comoving
+       * check */
+      /* reset quantities to nh, u, and z that we want to test */
+      ti_current = max_nr_timesteps;
+      cosmology_update(&cosmo, &phys_const, ti_current);
+      set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh_cgs, u_cgs,
+                     ti_current);
+
+      /* Set dt */
+      const integertime_t ti_step = get_integer_timestep(timebin);
+      const integertime_t ti_begin =
+          get_integer_time_begin(ti_current - 1, timebin);
+      dt_cool = cosmology_get_delta_time(&cosmo, ti_begin, ti_begin + ti_step);
+      dt_therm =
+          cosmology_get_therm_kick_factor(&cosmo, ti_begin, ti_begin + ti_step);
+
+      cooling_init(params, &us, &phys_const, &hydro_properties, &cooling);
+      cooling_update(&cosmo, &cooling, 0);
+
+      /* compute implicit solution */
+      cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties,
+                        &floor_props, &cooling, &p, &xp, dt_cool, dt_therm);
+      du_dt_check = hydro_get_physical_internal_energy_dt(&p, &cosmo);
+
+      /* Now we can test the cooling at various redshifts and compare the result
+       * to the redshift zero solution */
+      for (int z_i = 0; z_i <= n_z; z_i++) {
+        ti_current = max_nr_timesteps / n_z * z_i + 1;
+
+        /* reset to get the comoving density */
+        cosmology_update(&cosmo, &phys_const, ti_current);
+        cosmo.z = 0.f;
+        set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const,
+                       nh_cgs * cosmo.a * cosmo.a * cosmo.a,
+                       u_cgs / cosmo.a2_inv, ti_current);
+
+        /* Load the appropriate tables */
+        cooling_init(params, &us, &phys_const, &hydro_properties, &cooling);
+        cooling_update(&cosmo, &cooling, 0);
+
+        /* compute implicit solution */
+        cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties,
+                          &floor_props, &cooling, &p, &xp, dt_cool, dt_therm);
+        du_dt_implicit = hydro_get_physical_internal_energy_dt(&p, &cosmo);
+
+        /* check if the two solutions are consistent */
+        if (fabs((du_dt_implicit - du_dt_check) / du_dt_check) >
+                integration_tolerance ||
+            (du_dt_check == 0.0 && du_dt_implicit != 0.0))
+          error(
+              "Solutions do not match. scale factor %.5e z %.5e nh_cgs %.5e "
+              "u_cgs %.5e dt (years) %.5e du cgs implicit %.5e reference %.5e "
+              "error %.5e",
+              cosmo.a, cosmo.z, nh_cgs, u_cgs,
+              dt_cool * units_cgs_conversion_factor(&us, UNIT_CONV_TIME) /
+                  seconds_per_year,
+              du_dt_implicit *
+                  units_cgs_conversion_factor(&us,
+                                              UNIT_CONV_ENERGY_PER_UNIT_MASS) *
+                  dt_therm,
+              du_dt_check *
+                  units_cgs_conversion_factor(&us,
+                                              UNIT_CONV_ENERGY_PER_UNIT_MASS) *
+                  dt_therm,
+              fabs((du_dt_implicit - du_dt_check) / du_dt_check));
+      }
+    }
+  }
+  message("done comoving cooling test");
+
+  free(params);
+  return 0;
+}
diff --git a/tests/testCooling.c b/tests/testCooling.c
index 727a9638b09b871e866fe787438a5707fd43ec6b..d0efcb9f3b212b09d0342ccea3040c00eb39f2af 100644
--- a/tests/testCooling.c
+++ b/tests/testCooling.c
@@ -19,10 +19,11 @@
 #include "../config.h"
 
 /* Local headers. */
+#include "../src/cooling/EAGLE/cooling_rates.h"
+#include "../src/cooling/EAGLE/cooling_tables.h"
+#include "../src/cooling/EAGLE/interpolate.h"
 #include "swift.h"
 
-#if 0
-
 /*
  * @brief Assign particle density and entropy corresponding to the
  * hydrogen number density and internal energy specified.
@@ -42,29 +43,24 @@ void set_quantities(struct part *restrict p, struct xpart *restrict xp,
                     const struct cooling_function_data *restrict cooling,
                     struct cosmology *restrict cosmo,
                     const struct phys_const *restrict phys_const, float nh_cgs,
-                    double u, integertime_t ti_current) {
-
-  /* Update cosmology quantities */
-  cosmology_update(cosmo, phys_const, ti_current);
-
+                    double u_cgs, integertime_t ti_current) {
   /* calculate density */
-  double hydrogen_number_density = nh_cgs / cooling->number_density_scale;
+  double hydrogen_number_density = nh_cgs / cooling->number_density_to_cgs;
   p->rho = hydrogen_number_density * phys_const->const_proton_mass /
-           p->chemistry_data.metal_mass_fraction[chemistry_element_H] *
-           (cosmo->a * cosmo->a * cosmo->a);
+           p->chemistry_data.metal_mass_fraction[chemistry_element_H];
 
   /* update entropy based on internal energy */
-  float pressure = (u * cosmo->a * cosmo->a) / cooling->internal_energy_scale *
-                   p->rho * (hydro_gamma_minus_one);
+  float pressure = (u_cgs)*cooling->internal_energy_from_cgs * p->rho *
+                   (hydro_gamma_minus_one);
   p->entropy = pressure * (pow(p->rho, -hydro_gamma));
   xp->entropy_full = p->entropy;
+
+  p->entropy_dt = 0.f;
 }
 
 /*
- * @brief Produces contributions to cooling rates for different
- * hydrogen number densities, from different metals,
- * tests 1d and 4d table interpolations produce
- * same results for cooling rate, dlambda/du and temperature.
+ * @brief Tests cooling integration scheme by comparing EAGLE
+ * integration to subcycled explicit equation.
  */
 int main(int argc, char **argv) {
   // Declare relevant structs
@@ -78,14 +74,16 @@ int main(int argc, char **argv) {
   struct cosmology cosmo;
   char *parametersFileName = "./testCooling.yml";
 
-  float nh;  // hydrogen number density
-  double u;  // internal energy
+  float nh_cgs;  // hydrogen number density
+  double u_cgs;  // internal energy
+
+  const float seconds_per_year = 3.154e7;
 
   /* Number of values to test for in redshift,
    * hydrogen number density and internal energy */
-  const int n_z = 50;
-  const int n_nh = 50;
-  const int n_u = 50;
+  const int n_z = 10;
+  const int n_nh = 10;
+  const int n_u = 10;
 
   /* Number of subcycles and tolerance used to compare
    * subcycled and implicit solution. Note, high value
@@ -94,10 +92,6 @@ int main(int argc, char **argv) {
   const int n_subcycle = 1000;
   const float integration_tolerance = 0.2;
 
-  /* Set dt */
-  const float dt_cool = 1.0e-5;
-  const float dt_therm = 1.0e-5;
-
   /* Read the parameter file */
   if (params == NULL) error("Error allocating memory for the parameter file.");
   message("Reading runtime parameters from file '%s'", parametersFileName);
@@ -110,17 +104,35 @@ int main(int argc, char **argv) {
   /* Init chemistry */
   chemistry_init(params, &us, &phys_const, &chem_data);
   chemistry_first_init_part(&phys_const, &us, &cosmo, &chem_data, &p, &xp);
+  chemistry_part_has_no_neighbours(&p, &xp, &chem_data, &cosmo);
   chemistry_print(&chem_data);
 
   /* Init cosmology */
   cosmology_init(params, &us, &phys_const, &cosmo);
   cosmology_print(&cosmo);
 
+  /* Set dt */
+  const int timebin = 38;
+  float dt_cool, dt_therm;
+
+  /* Init hydro_props */
+  struct hydro_props hydro_properties;
+  hydro_props_init(&hydro_properties, &phys_const, &us, params);
+
   /* Init cooling */
-  cooling_init(params, &us, &phys_const, &cooling);
+  cooling_init(params, &us, &phys_const, &hydro_properties, &cooling);
   cooling_print(&cooling);
   cooling_update(&cosmo, &cooling, 0);
 
+  /* Init entropy floor */
+  struct entropy_floor_properties floor_props;
+  entropy_floor_init(&floor_props, &phys_const, &us, &hydro_properties, params);
+
+  /* Cooling function needs to know the minimal energy. Set it to the lowest
+   * internal energy in the cooling table. */
+  hydro_properties.minimal_internal_energy =
+      exp(M_LN10 * cooling.Therm[0]) * cooling.internal_energy_from_cgs;
+
   /* Calculate abundance ratios */
   float *abundance_ratio;
   abundance_ratio = malloc((chemistry_element_count + 2) * sizeof(float));
@@ -133,72 +145,89 @@ int main(int argc, char **argv) {
       (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
   int He_i;
   float d_He;
-  get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
-
-  /* Cooling function needs to know the minimal energy. Set it to the lowest
-   * internal energy in the cooling table. */
-  struct hydro_props hydro_properties;
-  hydro_properties.minimal_internal_energy =
-      exp(M_LN10 * cooling.Therm[0]) / cooling.internal_energy_scale;
+  get_index_1d(cooling.HeFrac, eagle_cooling_N_He_frac, HeFrac, &He_i, &d_He);
 
   /* calculate spacing in nh and u */
-  const float delta_nh = (cooling.nH[cooling.N_nH - 1] - cooling.nH[0]) / n_nh;
-  const float delta_u =
-      (cooling.Therm[cooling.N_Temp - 1] - cooling.Therm[0]) / n_u;
+  const float log_u_min_cgs = 11, log_u_max_cgs = 17;
+  const float log_nh_min_cgs = -6, log_nh_max_cgs = 3;
+  const float delta_log_nh_cgs = (log_nh_max_cgs - log_nh_min_cgs) / n_nh;
+  const float delta_log_u_cgs = (log_u_max_cgs - log_u_min_cgs) / n_u;
+
+  /* Declare variables we will be checking */
+  double du_dt_implicit, du_dt_check;
+  integertime_t ti_current;
+
+  /* Loop over values of nh and u */
+  for (int nh_i = 0; nh_i < n_nh; nh_i++) {
+    nh_cgs = exp(M_LN10 * log_nh_min_cgs + delta_log_nh_cgs * nh_i);
+    for (int u_i = 0; u_i < n_u; u_i++) {
+      u_cgs = exp(M_LN10 * log_u_min_cgs + delta_log_u_cgs * u_i);
 
-  for (int z_i = 0; z_i < n_z; z_i++) {
-    integertime_t ti_current = max_nr_timesteps / n_z * z_i;
-    for (int nh_i = 0; nh_i < n_nh; nh_i++) {
-      nh = exp(M_LN10 * cooling.nH[0] + delta_nh * nh_i);
-      for (int u_i = 0; u_i < n_u; u_i++) {
-        u = exp(M_LN10 * cooling.Therm[0] + delta_u * u_i);
+      /* Loop over z */
+      for (int z_i = 0; z_i <= n_z; z_i++) {
+        ti_current = max_nr_timesteps / n_z * z_i + 1;
 
         /* update nh, u, z */
-        set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh, u,
-                       ti_current);
+        cosmology_update(&cosmo, &phys_const, ti_current);
+        cooling_init(params, &us, &phys_const, &hydro_properties, &cooling);
+        cooling_update(&cosmo, &cooling, 0);
+        set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh_cgs,
+                       u_cgs, ti_current);
+
+        /* Set dt */
+        const integertime_t ti_step = get_integer_timestep(timebin);
+        const integertime_t ti_begin =
+            get_integer_time_begin(ti_current - 1, timebin);
+        dt_cool =
+            cosmology_get_delta_time(&cosmo, ti_begin, ti_begin + ti_step);
+        dt_therm = cosmology_get_therm_kick_factor(&cosmo, ti_begin,
+                                                   ti_begin + ti_step);
 
         /* calculate subcycled solution */
         for (int t_subcycle = 0; t_subcycle < n_subcycle; t_subcycle++) {
           p.entropy_dt = 0;
           cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties,
-                            &cooling, &p, &xp, dt_cool / n_subcycle,
-                            dt_therm / n_subcycle);
+                            &floor_props, &cooling, &p, &xp,
+                            dt_cool / n_subcycle, dt_therm / n_subcycle);
           xp.entropy_full += p.entropy_dt * dt_therm / n_subcycle;
         }
-        double u_subcycled =
-            hydro_get_physical_internal_energy(&p, &xp, &cosmo) *
-            cooling.internal_energy_scale;
+        du_dt_check = hydro_get_physical_internal_energy_dt(&p, &cosmo);
 
         /* reset quantities to nh, u, and z that we want to test */
-        set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh, u,
-                       ti_current);
+        cosmology_update(&cosmo, &phys_const, ti_current);
+        set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh_cgs,
+                       u_cgs, ti_current);
 
         /* compute implicit solution */
-        cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties, &cooling,
-                          &p, &xp, dt_cool, dt_therm);
-        double u_implicit =
-            hydro_get_physical_internal_energy(&p, &xp, &cosmo) *
-            cooling.internal_energy_scale;
+        cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties,
+                          &floor_props, &cooling, &p, &xp, dt_cool, dt_therm);
+        du_dt_implicit = hydro_get_physical_internal_energy_dt(&p, &cosmo);
 
         /* check if the two solutions are consistent */
-        if (fabs((u_implicit - u_subcycled) / u_subcycled) >
-            integration_tolerance)
-          message(
-              "implicit and subcycled solutions do not match. z_i %d nh_i %d "
-              "u_i %d implicit %.5e subcycled %.5e error %.5e",
-              z_i, nh_i, u_i, u_implicit, u_subcycled,
-              fabs((u_implicit - u_subcycled) / u_subcycled));
+        if (fabs((du_dt_implicit - du_dt_check) / du_dt_check) >
+                integration_tolerance ||
+            (du_dt_check == 0.0 && du_dt_implicit != 0.0))
+          error(
+              "Solutions do not match. scale factor %.5e z %.5e nh_cgs %.5e "
+              "u_cgs %.5e dt (years) %.5e du cgs implicit %.5e reference %.5e "
+              "error %.5e",
+              cosmo.a, cosmo.z, nh_cgs, u_cgs,
+              dt_cool * units_cgs_conversion_factor(&us, UNIT_CONV_TIME) /
+                  seconds_per_year,
+              du_dt_implicit *
+                  units_cgs_conversion_factor(&us,
+                                              UNIT_CONV_ENERGY_PER_UNIT_MASS) *
+                  dt_therm,
+              du_dt_check *
+                  units_cgs_conversion_factor(&us,
+                                              UNIT_CONV_ENERGY_PER_UNIT_MASS) *
+                  dt_therm,
+              fabs((du_dt_implicit - du_dt_check) / du_dt_check));
       }
     }
   }
-  message("done test");
+  message("done explicit subcycling cooling test");
 
   free(params);
   return 0;
 }
-
-#else
-
-int main(int argc, char **argv) { return 0; }
-
-#endif
diff --git a/tests/testCooling.yml b/tests/testCooling.yml
index faec32cdfec20b48af7341889c79b60bd2f6bb5b..fab3ae8959112a078e92621a0031ee14dbc2bc0d 100644
--- a/tests/testCooling.yml
+++ b/tests/testCooling.yml
@@ -9,7 +9,7 @@ InternalUnitSystem:
 # Cosmological parameters
 Cosmology:
   h:              0.6777        # Reduced Hubble constant
-  a_begin:        0.1     # Initial scale-factor of the simulation
+  a_begin:        0.5           # Initial scale-factor of the simulation
   a_end:          1.0           # Final scale factor of the simulation
   Omega_m:        0.307         # Matter density parameter
   Omega_lambda:   0.693         # Dark-energy density parameter
@@ -81,27 +81,51 @@ GrackleCooling:
   MaxSteps: 1000
   ConvergenceLimit: 1e-2
 
-EagleCooling:
-  filename:                /cosma5/data/Eagle/BG_Tables/CoolingTables/
-  reionisation_redshift:   8.989
+EAGLECooling:
+  dir_name:                ./coolingtables/
+  H_reion_z:               11.5
+  H_reion_eV_p_H:          2.0
   He_reion_z_centre:       3.5
   He_reion_z_sigma:        0.5
-  He_reion_ev_pH:          2.0
+  He_reion_eV_p_H:         2.0
+
+#EAGLEChemistry:
+#  InitMetallicity:         0.014
+#  InitAbundance_Hydrogen:  0.70649785
+#  InitAbundance_Helium:    0.28055534
+#  InitAbundance_Carbon:    2.0665436e-3
+#  InitAbundance_Nitrogen:  8.3562563e-4
+#  InitAbundance_Oxygen:    5.4926244e-3
+#  InitAbundance_Neon:      1.4144605e-3
+#  InitAbundance_Magnesium: 5.907064e-4
+#  InitAbundance_Silicon:   6.825874e-4
+#  InitAbundance_Iron:      1.1032152e-3
+#  CalciumOverSilicon:      0.0941736
+#  SulphurOverSilicon:      0.6054160
+EAGLEChemistry:              # Solar abundances
+  init_abundance_metal:      0.014
+  init_abundance_Hydrogen:   0.70649785
+  init_abundance_Helium:     0.28055534
+  init_abundance_Carbon:     2.0665436e-3
+  init_abundance_Nitrogen:   8.3562563e-4
+  init_abundance_Oxygen:     5.4926244e-3
+  init_abundance_Neon:       1.4144605e-3
+  init_abundance_Magnesium:  5.907064e-4
+  init_abundance_Silicon:    6.825874e-4
+  init_abundance_Iron:       1.1032152e-3
 
-EAGLEChemistry:
-  InitMetallicity:         0.014
-  InitAbundance_Hydrogen:  0.70649785
-  InitAbundance_Helium:    0.28055534
-  InitAbundance_Carbon:    2.0665436e-3
-  InitAbundance_Nitrogen:  8.3562563e-4
-  InitAbundance_Oxygen:    5.4926244e-3
-  InitAbundance_Neon:      1.4144605e-3
-  InitAbundance_Magnesium: 5.907064e-4
-  InitAbundance_Silicon:   6.825874e-4
-  InitAbundance_Iron:      1.1032152e-3
-  CalciumOverSilicon:      0.0941736
-  SulphurOverSilicon:      0.6054160
 
 GearChemistry:
   InitialMetallicity: 0.01295
 
+# Parameters for the EAGLE "equation of state"
+EAGLEEntropyFloor:
+  Jeans_density_threshold_H_p_cm3: 0.1       # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+  Jeans_over_density_threshold:    10.       # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
+  Jeans_temperature_norm_K:        8000      # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
+  Jeans_gamma_effective:           1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
+  Cool_density_threshold_H_p_cm3: 1e-5       # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+  Cool_over_density_threshold:    10.        # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
+  Cool_temperature_norm_K:        8000       # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
+  Cool_gamma_effective:           1.         # Slope the of the EAGLE Cool limiter entropy floor
+