diff --git a/configure.ac b/configure.ac
index eff47c85f76c19176619476fc3f01783156a6765..d770dd06e99e1083a2a02a7ac9a1b4a4dac38646 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1443,7 +1443,7 @@ AC_SUBST([NUMA_INCS])
 # As an example for this, see the call to AC_ARG_WITH for cooling.
 AC_ARG_WITH([subgrid],
 	[AS_HELP_STRING([--with-subgrid=<subgrid>],
-		[Master switch for subgrid methods. Inexperienced user should start from here @<:@none, GEAR, EAGLE default: none@:>@]
+		[Master switch for subgrid methods. Inexperienced user should start here. Options are: @<:@none, GEAR, QLA, EAGLE default: none@:>@]
 	)],
 	[with_subgrid="$withval"],
 	[with_subgrid=none]
@@ -1478,6 +1478,18 @@ case "$with_subgrid" in
 	with_subgrid_task_order=none
 	enable_fof=no
    ;;
+   QLA)
+	with_subgrid_cooling=QLA
+	with_subgrid_chemistry=QLA
+	with_subgrid_tracers=none
+	with_subgrid_entropy_floor=QLA
+	with_subgrid_stars=none
+	with_subgrid_star_formation=QLA
+	with_subgrid_feedback=none
+	with_subgrid_black_holes=none
+	with_subgrid_task_order=default
+	enable_fof=no
+   ;;
    EAGLE)
 	with_subgrid_cooling=EAGLE
 	with_subgrid_chemistry=EAGLE
@@ -1810,6 +1822,9 @@ case "$with_cooling" in
       primordial_chemistry=${with_cooling:8}
       AC_DEFINE_UNQUOTED([COOLING_GRACKLE_MODE], [$primordial_chemistry], [Grackle chemistry network])
    ;;
+   QLA)
+      AC_DEFINE([COOLING_QLA], [1], [Cooling following the Quick-Lyman-alpha model])
+   ;;
    EAGLE)
       AC_DEFINE([COOLING_EAGLE], [1], [Cooling following the EAGLE model])
    ;;
@@ -1845,6 +1860,9 @@ case "$with_chemistry" in
       number_element=${with_chemistry:5}
       AC_DEFINE_UNQUOTED([GEAR_CHEMISTRY_ELEMENT_COUNT], [$number_element], [Number of element to follow])
    ;;
+   QLA)
+      AC_DEFINE([CHEMISTRY_QLA], [1], [Chemistry taken from the Quick-Lyman-alpha model])
+   ;;
    EAGLE)
       AC_DEFINE([CHEMISTRY_EAGLE], [1], [Chemistry taken from the EAGLE model])
    ;;
@@ -2083,6 +2101,9 @@ case "$with_entropy_floor" in
    none)
       AC_DEFINE([ENTROPY_FLOOR_NONE], [1], [No entropy floor])
    ;;
+   QLA)
+      AC_DEFINE([ENTROPY_FLOOR_QLA], [1], [Quick Lyman-alpha entropy floor])
+   ;;
    EAGLE)
       AC_DEFINE([ENTROPY_FLOOR_EAGLE], [1], [EAGLE entropy floor])
    ;;
@@ -2140,6 +2161,9 @@ case "$with_star_formation" in
    none)
       AC_DEFINE([STAR_FORMATION_NONE], [1], [No star formation])
    ;;
+   QLA)
+      AC_DEFINE([STAR_FORMATION_QLA], [1], [Quick Lyman-alpha star formation model)])
+   ;;
    EAGLE)
       AC_DEFINE([STAR_FORMATION_EAGLE], [1], [EAGLE star formation model (Schaye and Dalla Vecchia (2008))])
    ;;
@@ -2180,6 +2204,9 @@ AC_SUBST([GIT_CMD])
 DX_INIT_DOXYGEN(libswift,doc/Doxyfile,doc/)
 AM_CONDITIONAL([HAVE_DOXYGEN], [test "$ac_cv_path_ac_pt_DX_DOXYGEN" != ""])
 
+# Check if using QLA cooling
+AM_CONDITIONAL([HAVEQLACOOLING], [test $with_cooling = "QLA"])
+
 # Check if using EAGLE cooling
 AM_CONDITIONAL([HAVEEAGLECOOLING], [test $with_cooling = "EAGLE"])
 
diff --git a/src/Makefile.am b/src/Makefile.am
index 7678220bbf0a05fab3af60248245493c3be9886b..dc9e48f51c92b68db7d0f2a24ae10e69afdd805b 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -59,6 +59,12 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     black_holes.h black_holes_io.h black_holes_properties.h black_holes_struct.h \
     feedback.h feedback_struct.h feedback_properties.h task_order.h
 
+# source files for EAGLE cooling
+QLA_COOLING_SOURCES =
+if HAVEQLACOOLING
+QLA_COOLING_SOURCES += cooling/QLA/cooling.c cooling/QLA/cooling_tables.c
+endif
+
 # source files for EAGLE cooling
 EAGLE_COOLING_SOURCES =
 if HAVEEAGLECOOLING
@@ -102,6 +108,7 @@ AM_SOURCES = space.c runner_main.c runner_doiact_hydro.c runner_doiact_limiter.c
     chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.c \
     outputlist.c velociraptor_dummy.c logger_io.c memuse.c mpiuse.c memuse_rnodes.c fof.c \
     hashmap.c pressure_floor.c \
+    $(QLA_COOLING_SOURCES) \
     $(EAGLE_COOLING_SOURCES) $(EAGLE_FEEDBACK_SOURCES) \
     $(GRACKLE_COOLING_SOURCES) $(GEAR_FEEDBACK_SOURCES)
 
@@ -222,6 +229,8 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
                  cooling/grackle/cooling_io.h \
 		 cooling/EAGLE/cooling.h cooling/EAGLE/cooling_struct.h cooling/EAGLE/cooling_tables.h \
                  cooling/EAGLE/cooling_io.h cooling/EAGLE/interpolate.h cooling/EAGLE/cooling_rates.h \
+		 cooling/QLA/cooling.h cooling/QLA/cooling_struct.h cooling/QLA/cooling_tables.h \
+                 cooling/QLA/cooling_io.h cooling/QLA/interpolate.h cooling/QLA/cooling_rates.h \
                  chemistry/none/chemistry.h \
 		 chemistry/none/chemistry_io.h \
 		 chemistry/none/chemistry_struct.h \
@@ -234,6 +243,10 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
 		 chemistry/EAGLE/chemistry_io.h \
 		 chemistry/EAGLE/chemistry_struct.h\
 		 chemistry/EAGLE/chemistry_iact.h \
+                 chemistry/QLA/chemistry.h \
+		 chemistry/QLA/chemistry_io.h \
+		 chemistry/QLA/chemistry_struct.h\
+		 chemistry/QLA/chemistry_iact.h \
 	         entropy_floor/none/entropy_floor.h \
                  entropy_floor/EAGLE/entropy_floor.h \
 		 tracers/none/tracers.h tracers/none/tracers_struct.h \
diff --git a/src/chemistry.h b/src/chemistry.h
index f9daa41db22a69a09f06be3fb560a68edac2f078..c6113cbbc64e3d1dcbd4c0cea0fdef25b1aaa75f 100644
--- a/src/chemistry.h
+++ b/src/chemistry.h
@@ -35,6 +35,9 @@
 #elif defined(CHEMISTRY_GEAR)
 #include "./chemistry/GEAR/chemistry.h"
 #include "./chemistry/GEAR/chemistry_iact.h"
+#elif defined(CHEMISTRY_QLA)
+#include "./chemistry/QLA/chemistry.h"
+#include "./chemistry/QLA/chemistry_iact.h"
 #elif defined(CHEMISTRY_EAGLE)
 #include "./chemistry/EAGLE/chemistry.h"
 #include "./chemistry/EAGLE/chemistry_iact.h"
diff --git a/src/chemistry/QLA/chemistry.h b/src/chemistry/QLA/chemistry.h
new file mode 100644
index 0000000000000000000000000000000000000000..69139ea6320b228df5b9990e93fb0377490839b9
--- /dev/null
+++ b/src/chemistry/QLA/chemistry.h
@@ -0,0 +1,321 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_CHEMISTRY_QLA_H
+#define SWIFT_CHEMISTRY_QLA_H
+
+/**
+ * @file src/chemistry/none/chemistry.h
+ * @brief Empty infrastructure for the cases without chemistry function
+ */
+
+/* Some standard headers. */
+#include <float.h>
+#include <math.h>
+
+/* Local includes. */
+#include "chemistry_struct.h"
+#include "cosmology.h"
+#include "error.h"
+#include "hydro.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/**
+ * @brief Return a string containing the name of a given #chemistry_element.
+ */
+__attribute__((always_inline)) INLINE static const char*
+chemistry_get_element_name(enum chemistry_element elem) {
+
+  static const char* chemistry_element_names[chemistry_element_count] = {};
+
+  return chemistry_element_names[elem];
+}
+
+/**
+ * @brief Initialises the chemistry properties.
+ *
+ * Nothing to do here.
+ *
+ * @param parameter_file The parsed parameter file.
+ * @param us The current internal system of units.
+ * @param phys_const The physical constants in internal units.
+ * @param data The global chemistry information (to be filled).
+ */
+static INLINE void chemistry_init_backend(struct swift_params* parameter_file,
+                                          const struct unit_system* us,
+                                          const struct phys_const* phys_const,
+                                          struct chemistry_global_data* data) {}
+
+/**
+ * @brief Prints the properties of the chemistry model to stdout.
+ *
+ * Nothing to do here.
+ *
+ * @brief The #chemistry_global_data containing information about the current
+ * model.
+ */
+static INLINE void chemistry_print_backend(
+    const struct chemistry_global_data* data) {
+
+  message("Chemistry function is 'Quick Lyman-alpha'.");
+}
+
+/**
+ * @brief Finishes the density calculation.
+ *
+ * Nothing to do here.
+ *
+ * @param p The particle to act upon
+ * @param cd The global chemistry information.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_end_density(
+    struct part* restrict p, const struct chemistry_global_data* cd,
+    const struct cosmology* cosmo) {}
+
+/**
+ * @brief Updates to the chemistry data after the hydro force loop.
+ *
+ * @param p The particle to act upon.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_end_force(
+    struct part* restrict p, const struct cosmology* cosmo) {}
+
+/**
+ * @brief Computes the chemistry-related time-step constraint.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param cosmo The current cosmological model.
+ * @param us The internal system of units.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param cd The global properties of the chemistry scheme.
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float chemistry_timestep(
+    const struct phys_const* restrict phys_const,
+    const struct cosmology* restrict cosmo,
+    const struct unit_system* restrict us,
+    const struct hydro_props* hydro_props,
+    const struct chemistry_global_data* cd, const struct part* restrict p) {
+  return FLT_MAX;
+}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ *
+ * Nothing to do here.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cd #chemistry_global_data containing chemistry informations.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void
+chemistry_part_has_no_neighbours(struct part* restrict p,
+                                 struct xpart* restrict xp,
+                                 const struct chemistry_global_data* cd,
+                                 const struct cosmology* cosmo) {}
+
+/**
+ * @brief Sets the chemistry properties of the (x-)particles to a valid start
+ * state.
+ *
+ * Nothing to do here.
+ *
+ * @param phys_const The physical constant in internal units.
+ * @param us The unit system.
+ * @param cosmo The current cosmological model.
+ * @param data The global chemistry information used for this run.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_first_init_part(
+    const struct phys_const* restrict phys_const,
+    const struct unit_system* restrict us,
+    const struct cosmology* restrict cosmo,
+    const struct chemistry_global_data* data, const struct part* restrict p,
+    struct xpart* restrict xp) {}
+
+/**
+ * @brief Sets the chemistry properties of the (x-)particles to a valid start
+ * state.
+ *
+ * Nothing to do here.
+ *
+ * @param p Pointer to the particle data.
+ * @param data The global chemistry information.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_init_part(
+    struct part* restrict p, const struct chemistry_global_data* data) {}
+
+/**
+ * @brief Sets the chemistry properties of the sparticles to a valid start
+ * state.
+ *
+ * Nothing to do here.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param cosmo The current cosmological model.
+ * @param data The global chemistry information.
+ * @param sp Pointer to the sparticle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_first_init_spart(
+    const struct chemistry_global_data* data, struct spart* restrict sp) {}
+
+/**
+ * @brief Initialise the chemistry properties of a black hole with
+ * the chemistry properties of the gas it is born from.
+ *
+ * Nothing to do here.
+ *
+ * @param bp_data The black hole data to initialise.
+ * @param p_data The gas data to use.
+ * @param gas_mass The mass of the gas particle.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_bpart_from_part(
+    struct chemistry_bpart_data* bp_data,
+    const struct chemistry_part_data* p_data, const double gas_mass) {}
+
+/**
+ * @brief Add the chemistry data of a gas particle to a black hole.
+ *
+ * Nothing to do here.
+ *
+ * @param bp_data The black hole data to add to.
+ * @param p_data The gas data to use.
+ * @param gas_mass The mass of the gas particle.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_add_part_to_bpart(
+    struct chemistry_bpart_data* bp_data,
+    const struct chemistry_part_data* p_data, const double gas_mass) {}
+
+/**
+ * @brief Add the chemistry data of a black hole to another one.
+ *
+ * Nothing to do here.
+ *
+ * @param bp_data The black hole data to add to.
+ * @param swallowed_data The black hole data to use.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_add_bpart_to_bpart(
+    struct chemistry_bpart_data* bp_data,
+    const struct chemistry_bpart_data* swallowed_data) {}
+
+/**
+ * @brief Split the metal content of a particle into n pieces
+ *
+ * Nothing to do here.
+ *
+ * @param p The #part.
+ * @param n The number of pieces to split into.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_split_part(
+    struct part* p, const double n) {}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * star particle to be used in feedback/enrichment related routines.
+ *
+ * No metallicity treatment here -> return 0.
+ *
+ * @param sp Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_total_metal_mass_fraction_for_feedback(const struct spart* sp) {
+
+  return 0.f;
+}
+
+/**
+ * @brief Returns the abundance array (metal mass fractions) of the
+ * star particle to be used in feedback/enrichment related routines.
+ *
+ * No metallicity treatment here -> return NULL array.
+ *
+ * @param sp Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float const*
+chemistry_get_metal_mass_fraction_for_feedback(const struct spart* sp) {
+
+  return NULL;
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * gas particle to be used in cooling related routines.
+ *
+ * No metallicity treatment here -> return 0.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_total_metal_mass_fraction_for_cooling(const struct part* p) {
+
+  return 0.f;
+}
+
+/**
+ * @brief Returns the abundance array (metal mass fractions) of the
+ * gas particle to be used in cooling related routines.
+ *
+ * No metallicity treatment here -> return NULL array.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float const*
+chemistry_get_metal_mass_fraction_for_cooling(const struct part* p) {
+
+  return NULL;
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * gas particle to be used in star formation related routines.
+ *
+ * No metallicity treatment here -> return 0.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_total_metal_mass_fraction_for_star_formation(
+    const struct part* p) {
+
+  return 0.f;
+}
+
+/**
+ * @brief Returns the abundance array (metal mass fractions) of the
+ * gas particle to be used in star formation related routines.
+ *
+ * No metallicity treatment here -> return NULL array.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float const*
+chemistry_get_metal_mass_fraction_for_star_formation(const struct part* p) {
+
+  return NULL;
+}
+
+#endif /* SWIFT_CHEMISTRY_QLA_H */
diff --git a/src/chemistry/QLA/chemistry_iact.h b/src/chemistry/QLA/chemistry_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..89f0d4a3c8a2dd4fbde40b4525196d439ba66f85
--- /dev/null
+++ b/src/chemistry/QLA/chemistry_iact.h
@@ -0,0 +1,61 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_QLA_CHEMISTRY_IACT_H
+#define SWIFT_QLA_CHEMISTRY_IACT_H
+
+/**
+ * @file none/chemistry_iact.h
+ * @brief Density computation
+ */
+
+/**
+ * @brief do chemistry computation after the runner_iact_density (symmetric
+ * version)
+ *
+ * @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_chemistry(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H) {}
+
+/**
+ * @brief do chemistry computation after the runner_iact_density (non symmetric
+ * version)
+ *
+ * @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_chemistry(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    const struct part *restrict pj, float a, float H) {}
+
+#endif /* SWIFT_QLA_CHEMISTRY_IACT_H */
diff --git a/src/chemistry/QLA/chemistry_io.h b/src/chemistry/QLA/chemistry_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..8a7b78cbfc5795ce93cc71f0a6772131af090f05
--- /dev/null
+++ b/src/chemistry/QLA/chemistry_io.h
@@ -0,0 +1,105 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_CHEMISTRY_IO_QLA_H
+#define SWIFT_CHEMISTRY_IO_QLA_H
+
+#include "io_properties.h"
+
+/**
+ * @brief Specifies which particle fields to read from a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to read.
+ *
+ * @return Returns the number of fields to write.
+ */
+INLINE static int chemistry_read_particles(struct part* parts,
+                                           struct io_props* list) {
+
+  /* update list according to hydro_io */
+
+  /* Return the number of fields to read */
+  return 0;
+}
+
+/**
+ * @brief Specifies which particle fields to write to a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return Returns the number of fields to write.
+ */
+INLINE static int chemistry_write_particles(const struct part* parts,
+                                            struct io_props* list) {
+
+  /* update list according to hydro_io */
+
+  /* Return the number of fields to write */
+  return 0;
+}
+
+/**
+ * @brief Specifies which sparticle fields to write to a dataset
+ *
+ * @param sparts The sparticle array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return Returns the number of fields to write.
+ */
+INLINE static int chemistry_write_sparticles(const struct spart* sparts,
+                                             struct io_props* list) {
+
+  /* update list according to hydro_io */
+
+  /* Return the number of fields to write */
+  return 0;
+}
+
+/**
+ * @brief Specifies which bparticle fields to write to a dataset
+ *
+ * @param bparts The bparticle array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return Returns the number of fields to write.
+ */
+INLINE static int chemistry_write_bparticles(const struct bpart* bparts,
+                                             struct io_props* list) {
+
+  /* update list according to hydro_io */
+
+  /* Return the number of fields to write */
+  return 0;
+}
+
+#ifdef HAVE_HDF5
+
+/**
+ * @brief Writes the current model of chemistry to the file
+ * @param h_grp The HDF5 group in which to write
+ */
+INLINE static void chemistry_write_flavour(hid_t h_grp) {
+
+  io_write_attribute_s(h_grp, "Chemistry Model", "Quick Lyman-alpha");
+  io_write_attribute_d(h_grp, "Chemistry element count", 0);
+}
+#endif
+
+#endif /* SWIFT_CHEMISTRY_IO_QLA_H */
diff --git a/src/chemistry/QLA/chemistry_struct.h b/src/chemistry/QLA/chemistry_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..97194433b1d1e33a5ae66d67f344ac48d58f71ce
--- /dev/null
+++ b/src/chemistry/QLA/chemistry_struct.h
@@ -0,0 +1,53 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_CHEMISTRY_STRUCT_QLA_H
+#define SWIFT_CHEMISTRY_STRUCT_QLA_H
+
+/**
+ * @file src/chemistry/none/chemistry_struct.h
+ * @brief Empty infrastructure for the cases without chemistry function
+ */
+
+/**
+ * @brief The individual elements traced in the model.
+ */
+enum chemistry_element { chemistry_element_count = 0 };
+
+/**
+ * @brief Global chemical abundance information.
+ *
+ * Nothing here.
+ */
+struct chemistry_global_data {};
+
+/**
+ * @brief Chemistry properties carried by the #part.
+ *
+ * Nothing here.
+ */
+struct chemistry_part_data {};
+
+/**
+ * @brief Chemistry properties carried by the #bpart.
+ *
+ * Nothing here.
+ */
+struct chemistry_bpart_data {};
+
+#endif /* SWIFT_CHEMISTRY_STRUCT_QLA_H */
diff --git a/src/chemistry_io.h b/src/chemistry_io.h
index 877f96b0488164f8b06ceffc45aa3e1fa5809937..c94461c84bdf243990fa8f28f0049178e180d9ec 100644
--- a/src/chemistry_io.h
+++ b/src/chemistry_io.h
@@ -27,6 +27,8 @@
 #include "./chemistry/none/chemistry_io.h"
 #elif defined(CHEMISTRY_GEAR)
 #include "./chemistry/GEAR/chemistry_io.h"
+#elif defined(CHEMISTRY_QLA)
+#include "./chemistry/QLA/chemistry_io.h"
 #elif defined(CHEMISTRY_EAGLE)
 #include "./chemistry/EAGLE/chemistry_io.h"
 #else
diff --git a/src/chemistry_struct.h b/src/chemistry_struct.h
index ffeecb3cef172416075a7413c94de9384bd0e7f7..05074d8759be660876544aaa20cffec6cddb8bce 100644
--- a/src/chemistry_struct.h
+++ b/src/chemistry_struct.h
@@ -32,6 +32,8 @@
 #include "./chemistry/none/chemistry_struct.h"
 #elif defined(CHEMISTRY_GEAR)
 #include "./chemistry/GEAR/chemistry_struct.h"
+#elif defined(CHEMISTRY_QLA)
+#include "./chemistry/QLA/chemistry_struct.h"
 #elif defined(CHEMISTRY_EAGLE)
 #include "./chemistry/EAGLE/chemistry_struct.h"
 #else
diff --git a/src/cooling.h b/src/cooling.h
index 8592025234b28cb4ff74524cffd407ff2c6f708b..bb0fbc58e67358a4ccac7fea2370b9f802575985 100644
--- a/src/cooling.h
+++ b/src/cooling.h
@@ -45,6 +45,8 @@
 #include "./cooling/Compton/cooling.h"
 #elif defined(COOLING_GRACKLE)
 #include "./cooling/grackle/cooling.h"
+#elif defined(COOLING_QLA)
+#include "./cooling/QLA/cooling.h"
 #elif defined(COOLING_EAGLE)
 #include "./cooling/EAGLE/cooling.h"
 #else
diff --git a/src/cooling/QLA/cooling.c b/src/cooling/QLA/cooling.c
new file mode 100644
index 0000000000000000000000000000000000000000..9805261124644bfa8096f22162b9cca1b73380eb
--- /dev/null
+++ b/src/cooling/QLA/cooling.c
@@ -0,0 +1,911 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+/**
+ * @file src/cooling/QLA/cooling.c
+ * @brief Quick Lyman-alpha cooling functions
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <float.h>
+#include <hdf5.h>
+#include <math.h>
+#include <time.h>
+
+/* Local includes. */
+#include "chemistry.h"
+#include "cooling.h"
+#include "cooling_rates.h"
+#include "cooling_struct.h"
+#include "cooling_tables.h"
+#include "entropy_floor.h"
+#include "error.h"
+#include "exp10.h"
+#include "hydro.h"
+#include "interpolate.h"
+#include "io_properties.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "space.h"
+#include "units.h"
+
+/* Maximum number of iterations for bisection scheme */
+static const int bisection_max_iterations = 150;
+
+/* Tolerances for termination criteria. */
+static const float explicit_tolerance = 0.05;
+static const float bisection_tolerance = 1.0e-6;
+static const double bracket_factor = 1.5;
+
+/**
+ * @brief Find the index of the current redshift along the redshift dimension
+ * of the cooling tables.
+ *
+ * Since the redshift table is not evenly spaced, compare z with each
+ * table value in decreasing order starting with the previous redshift index
+ *
+ * The returned difference is expressed in units of the table separation. This
+ * means dx = (x - table[i]) / (table[i+1] - table[i]). It is always between
+ * 0 and 1.
+ *
+ * @param z Redshift we are searching for.
+ * @param z_index (return) Index of the redshift in the table.
+ * @param dz (return) Difference in redshift between z and table[z_index].
+ * @param cooling #cooling_function_data structure containing redshift table.
+ */
+__attribute__((always_inline)) INLINE void get_redshift_index(
+    const float z, int *z_index, float *dz,
+    struct cooling_function_data *restrict cooling) {
+
+  /* Before the earliest redshift or before hydrogen reionization, flag for
+   * collisional cooling */
+  if (z > cooling->H_reion_z) {
+    *z_index = qla_cooling_N_redshifts;
+    *dz = 0.0;
+  }
+
+  /* From reionization use the cooling tables */
+  else if (z > cooling->Redshifts[qla_cooling_N_redshifts - 1] &&
+           z <= cooling->H_reion_z) {
+    *z_index = qla_cooling_N_redshifts + 1;
+    *dz = 0.0;
+  }
+
+  /* At the end, just use the last value */
+  else if (z <= cooling->Redshifts[0]) {
+    *z_index = 0;
+    *dz = 0.0;
+  }
+
+  /* Normal case: search... */
+  else {
+
+    /* start at the previous index and search */
+    for (int i = cooling->previous_z_index; i >= 0; i--) {
+      if (z > cooling->Redshifts[i]) {
+
+        *z_index = i;
+        cooling->previous_z_index = i;
+
+        *dz = (z - cooling->Redshifts[i]) /
+              (cooling->Redshifts[i + 1] - cooling->Redshifts[i]);
+        break;
+      }
+    }
+  }
+}
+
+/**
+ * @brief Common operations performed on the cooling function at a
+ * given time-step or redshift. Predominantly used to read cooling tables
+ * above and below the current redshift, if not already read in.
+ *
+ * Also calls the additional H reionisation energy injection if need be.
+ *
+ * @param cosmo The current cosmological model.
+ * @param cooling The #cooling_function_data used in the run.
+ * @param s The space data, including a pointer to array of particles
+ */
+void cooling_update(const struct cosmology *cosmo,
+                    struct cooling_function_data *cooling, struct space *s) {
+
+  /* Current redshift */
+  const float redshift = cosmo->z;
+
+  /* What is the current table index along the redshift axis? */
+  int z_index = -1;
+  float dz = 0.f;
+  get_redshift_index(redshift, &z_index, &dz, cooling);
+  cooling->dz = dz;
+
+  /* Extra energy for reionization? */
+  if (!cooling->H_reion_done) {
+
+    /* Does this timestep straddle Hydrogen reionization? If so, we need to
+     * input extra heat */
+    if (cosmo->z <= cooling->H_reion_z && cosmo->z_old > cooling->H_reion_z) {
+
+      if (s == NULL) error("Trying to do H reionization on an empty space!");
+
+      /* Inject energy to all particles */
+      cooling_Hydrogen_reionization(cooling, cosmo, s);
+
+      /* Flag that reionization happened */
+      cooling->H_reion_done = 1;
+    }
+  }
+
+  /* Do we already have the correct tables loaded? */
+  if (cooling->z_index == z_index) return;
+
+  /* Which table should we load ? */
+  if (z_index >= qla_cooling_N_redshifts) {
+
+    if (z_index == qla_cooling_N_redshifts + 1) {
+
+      /* Bewtween re-ionization and first table */
+      get_redshift_invariant_table(cooling, /* photodis=*/0);
+
+    } else {
+
+      /* Above re-ionization */
+      get_redshift_invariant_table(cooling, /* photodis=*/1);
+    }
+
+  } else {
+
+    /* Normal case: two tables bracketing the current z */
+    const int low_z_index = z_index;
+    const int high_z_index = z_index + 1;
+
+    get_cooling_table(cooling, low_z_index, high_z_index);
+  }
+
+  /* Store the currently loaded index */
+  cooling->z_index = z_index;
+}
+
+/**
+ * @brief Bisection integration scheme
+ *
+ * @param u_ini_cgs Internal energy at beginning of hydro step in CGS.
+ * @param n_H_cgs Hydrogen number density in CGS.
+ * @param redshift Current redshift.
+ * @param n_H_index Particle hydrogen number density index.
+ * @param d_n_H Particle hydrogen number density offset.
+ * @param He_index Particle helium fraction index.
+ * @param d_He Particle helium fraction offset.
+ * @param Lambda_He_reion_cgs Cooling rate coming from He reionization.
+ * @param ratefact_cgs Multiplication factor to get a cooling rate.
+ * @param cooling #cooling_function_data structure.
+ * @param abundance_ratio Array of ratios of metal abundance to solar.
+ * @param dt_cgs timestep in CGS.
+ * @param ID ID of the particle (for debugging).
+ */
+INLINE static double bisection_iter(
+    const double u_ini_cgs, const double n_H_cgs, const double redshift,
+    const int n_H_index, const float d_n_H, const int He_index,
+    const float d_He, const double Lambda_He_reion_cgs,
+    const double ratefact_cgs,
+    const struct cooling_function_data *restrict cooling,
+    const float abundance_ratio[qla_cooling_N_abundances], const double dt_cgs,
+    const long long ID) {
+
+  /* Bracketing */
+  double u_lower_cgs = u_ini_cgs;
+  double u_upper_cgs = u_ini_cgs;
+
+  /*************************************/
+  /* Let's get a first guess           */
+  /*************************************/
+
+  double LambdaNet_cgs =
+      Lambda_He_reion_cgs +
+      qla_cooling_rate(log10(u_ini_cgs), redshift, n_H_cgs, abundance_ratio,
+                       n_H_index, d_n_H, He_index, d_He, cooling);
+
+  /*************************************/
+  /* Let's try to bracket the solution */
+  /*************************************/
+
+  if (LambdaNet_cgs < 0) {
+
+    /* we're cooling! */
+    u_lower_cgs /= bracket_factor;
+    u_upper_cgs *= bracket_factor;
+
+    /* Compute a new rate */
+    LambdaNet_cgs =
+        Lambda_He_reion_cgs +
+        qla_cooling_rate(log10(u_lower_cgs), redshift, n_H_cgs, abundance_ratio,
+                         n_H_index, d_n_H, He_index, d_He, cooling);
+
+    int i = 0;
+    while (u_lower_cgs - u_ini_cgs - LambdaNet_cgs * ratefact_cgs * dt_cgs >
+               0 &&
+           i < bisection_max_iterations) {
+
+      u_lower_cgs /= bracket_factor;
+      u_upper_cgs /= bracket_factor;
+
+      /* Compute a new rate */
+      LambdaNet_cgs = Lambda_He_reion_cgs +
+                      qla_cooling_rate(log10(u_lower_cgs), redshift, n_H_cgs,
+                                       abundance_ratio, n_H_index, d_n_H,
+                                       He_index, d_He, cooling);
+      i++;
+    }
+
+    if (i >= bisection_max_iterations) {
+      error(
+          "particle %llu exceeded max iterations searching for bounds when "
+          "cooling, u_ini_cgs %.5e n_H_cgs %.5e",
+          ID, u_ini_cgs, n_H_cgs);
+    }
+  } else {
+
+    /* we are heating! */
+    u_lower_cgs /= bracket_factor;
+    u_upper_cgs *= bracket_factor;
+
+    /* Compute a new rate */
+    LambdaNet_cgs =
+        Lambda_He_reion_cgs +
+        qla_cooling_rate(log10(u_upper_cgs), redshift, n_H_cgs, abundance_ratio,
+                         n_H_index, d_n_H, He_index, d_He, cooling);
+
+    int i = 0;
+    while (u_upper_cgs - u_ini_cgs - LambdaNet_cgs * ratefact_cgs * dt_cgs <
+               0 &&
+           i < bisection_max_iterations) {
+
+      u_lower_cgs *= bracket_factor;
+      u_upper_cgs *= bracket_factor;
+
+      /* Compute a new rate */
+      LambdaNet_cgs = Lambda_He_reion_cgs +
+                      qla_cooling_rate(log10(u_upper_cgs), redshift, n_H_cgs,
+                                       abundance_ratio, n_H_index, d_n_H,
+                                       He_index, d_He, cooling);
+      i++;
+    }
+
+    if (i >= bisection_max_iterations) {
+      error(
+          "particle %llu exceeded max iterations searching for bounds when "
+          "heating, u_ini_cgs %.5e n_H_cgs %.5e",
+          ID, u_ini_cgs, n_H_cgs);
+    }
+  }
+
+  /********************************************/
+  /* We now have an upper and lower bound.    */
+  /* Let's iterate by reducing the bracketing */
+  /********************************************/
+
+  /* bisection iteration */
+  int i = 0;
+  double u_next_cgs;
+
+  do {
+
+    /* New guess */
+    u_next_cgs = 0.5 * (u_lower_cgs + u_upper_cgs);
+
+    /* New rate */
+    LambdaNet_cgs =
+        Lambda_He_reion_cgs +
+        qla_cooling_rate(log10(u_next_cgs), redshift, n_H_cgs, abundance_ratio,
+                         n_H_index, d_n_H, He_index, d_He, cooling);
+#ifdef SWIFT_DEBUG_CHECKS
+    if (u_next_cgs <= 0)
+      error(
+          "Got negative energy! u_next_cgs=%.5e u_upper=%.5e u_lower=%.5e "
+          "Lambda=%.5e",
+          u_next_cgs, u_upper_cgs, u_lower_cgs, LambdaNet_cgs);
+#endif
+
+    /* Where do we go next? */
+    if (u_next_cgs - u_ini_cgs - LambdaNet_cgs * ratefact_cgs * dt_cgs > 0.0) {
+      u_upper_cgs = u_next_cgs;
+    } else {
+      u_lower_cgs = u_next_cgs;
+    }
+
+    i++;
+  } while (fabs(u_upper_cgs - u_lower_cgs) / u_next_cgs > bisection_tolerance &&
+           i < bisection_max_iterations);
+
+  if (i >= bisection_max_iterations)
+    error("Particle id %llu failed to converge", ID);
+
+  return u_upper_cgs;
+}
+
+/**
+ * @brief Apply the cooling function to a particle.
+ *
+ * We want to compute u_new such that u_new = u_old + dt * du/dt(u_new, X),
+ * where X stands for the metallicity, density and redshift. These are
+ * kept constant.
+ *
+ * We first compute du/dt(u_old). If dt * du/dt(u_old) is small enough, we
+ * use an explicit integration and use this as our solution.
+ *
+ * Otherwise, we try to find a solution to the implicit time-integration
+ * problem. This leads to the root-finding problem:
+ *
+ * f(u_new) = u_new - u_old - dt * du/dt(u_new, X) = 0
+ *
+ * We first try a few Newton-Raphson iteration if it does not converge, we
+ * revert to a bisection scheme.
+ *
+ * This is done by first bracketing the solution and then iterating
+ * towards the solution by reducing the window down to a certain tolerance.
+ * Note there is always at least one solution since
+ * f(+inf) is < 0 and f(-inf) is > 0.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param cosmo The current cosmological model.
+ * @param hydro_properties the hydro_props struct
+ * @param floor_props Properties of the entropy floor.
+ * @param cooling The #cooling_function_data used in the run.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ * @param dt The cooling time-step of this particle.
+ * @param dt_therm The hydro time-step of this particle.
+ * @param time The current time (since the Big Bang or start of the run) in
+ * internal units.
+ */
+void cooling_cool_part(const struct phys_const *phys_const,
+                       const struct unit_system *us,
+                       const struct cosmology *cosmo,
+                       const struct hydro_props *hydro_properties,
+                       const struct entropy_floor_properties *floor_props,
+                       const struct cooling_function_data *cooling,
+                       struct part *restrict p, struct xpart *restrict xp,
+                       const float dt, const float dt_therm,
+                       const double time) {
+
+  /* No cooling happens over zero time */
+  if (dt == 0.) return;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (cooling->Redshifts == NULL)
+    error(
+        "Cooling function has not been initialised. Did you forget the "
+        "--cooling runtime flag?");
+#endif
+
+  /* Get internal energy at the last kick step */
+  const float u_start = hydro_get_physical_internal_energy(p, xp, cosmo);
+
+  /* Get the change in internal energy due to hydro forces */
+  const float hydro_du_dt = hydro_get_physical_internal_energy_dt(p, cosmo);
+
+  /* Get internal energy at the end of the step (assuming dt does not
+   * increase) */
+  double u_0 = (u_start + hydro_du_dt * dt_therm);
+
+  /* Check for minimal energy */
+  u_0 = max(u_0, hydro_properties->minimal_internal_energy);
+
+  /* Convert to CGS units */
+  const double u_0_cgs = u_0 * cooling->internal_energy_to_cgs;
+  const double dt_cgs = dt * units_cgs_conversion_factor(us, UNIT_CONV_TIME);
+
+  /* Change in redshift over the course of this time-step
+     (See cosmology theory document for the derivation) */
+  const double delta_redshift = -dt * cosmo->H * cosmo->a_inv;
+
+  /* Get this particle's abundance ratios compared to solar
+   * Note that we need to add S and Ca that are in the tables but not tracked
+   * by the particles themselves.
+   * The order is [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe] */
+  float abundance_ratio[qla_cooling_N_abundances];
+  abundance_ratio_to_solar(p, cooling, phys_const, abundance_ratio);
+
+  /* Get the Hydrogen and Helium mass fractions */
+  const float XH = 1. - phys_const->const_primordial_He_fraction;
+  const float XHe = phys_const->const_primordial_He_fraction;
+
+  /* Get the Helium mass fraction. Note that this is He / (H + He), i.e. a
+   * metal-free Helium mass fraction as per the Wiersma+08 definition */
+  const float HeFrac = XHe / (XH + XHe);
+
+  /* convert Hydrogen mass fraction into physical Hydrogen number density */
+  const double n_H =
+      hydro_get_physical_density(p, cosmo) * XH / phys_const->const_proton_mass;
+  const double n_H_cgs = n_H * cooling->number_density_to_cgs;
+
+  /* ratefact = n_H * n_H / rho; Might lead to round-off error: replaced by
+   * equivalent expression  below */
+  const double ratefact_cgs = n_H_cgs * (XH * cooling->inv_proton_mass_cgs);
+
+  /* compute hydrogen number density and helium fraction table indices and
+   * offsets (These are fixed for any value of u, so no need to recompute them)
+   */
+  int He_index, n_H_index;
+  float d_He, d_n_H;
+  get_index_1d(cooling->HeFrac, qla_cooling_N_He_frac, HeFrac, &He_index,
+               &d_He);
+  get_index_1d(cooling->nH, qla_cooling_N_density, log10(n_H_cgs), &n_H_index,
+               &d_n_H);
+
+  /* Start by computing the cooling (heating actually) rate from Helium
+     re-ionization as this needs to be added on no matter what */
+
+  /* Get helium and hydrogen reheating term */
+  const double Helium_reion_heat_cgs =
+      qla_helium_reionization_extraheat(cosmo->z, delta_redshift, cooling);
+
+  /* Convert this into a rate */
+  const double Lambda_He_reion_cgs =
+      Helium_reion_heat_cgs / (dt_cgs * ratefact_cgs);
+
+  /* Let's compute the internal energy at the end of the step */
+  /* Initialise to the initial energy to appease compiler; this will never not
+     be overwritten. */
+  double u_final_cgs = u_0_cgs;
+
+  /* First try an explicit integration (note we ignore the derivative) */
+  const double LambdaNet_cgs =
+      Lambda_He_reion_cgs + qla_cooling_rate(log10(u_0_cgs), cosmo->z, n_H_cgs,
+                                             abundance_ratio, n_H_index, d_n_H,
+                                             He_index, d_He, cooling);
+
+  /* if cooling rate is small, take the explicit solution */
+  if (fabs(ratefact_cgs * LambdaNet_cgs * dt_cgs) <
+      explicit_tolerance * u_0_cgs) {
+
+    u_final_cgs = u_0_cgs + ratefact_cgs * LambdaNet_cgs * dt_cgs;
+
+  } else {
+
+    /* Otherwise, go the bisection route. */
+    u_final_cgs =
+        bisection_iter(u_0_cgs, n_H_cgs, cosmo->z, n_H_index, d_n_H, He_index,
+                       d_He, Lambda_He_reion_cgs, ratefact_cgs, cooling,
+                       abundance_ratio, dt_cgs, p->id);
+  }
+
+  /* Convert back to internal units */
+  double u_final = u_final_cgs * cooling->internal_energy_from_cgs;
+
+  /* We now need to check that we are not going to go below any of the limits */
+
+  /* Absolute minimum */
+  const double u_minimal = hydro_properties->minimal_internal_energy;
+  u_final = max(u_final, u_minimal);
+
+  /* Limit imposed by the entropy floor */
+  const double A_floor = entropy_floor(p, cosmo, floor_props);
+  const double rho_physical = hydro_get_physical_density(p, cosmo);
+  const double u_floor =
+      gas_internal_energy_from_entropy(rho_physical, A_floor);
+  u_final = max(u_final, u_floor);
+
+  /* Expected change in energy over the next kick step
+     (assuming no change in dt) */
+  const double delta_u = u_final - max(u_start, u_floor);
+
+  /* Turn this into a rate of change (including cosmology term) */
+  const float cooling_du_dt = delta_u / dt_therm;
+
+  /* Update the internal energy time derivative */
+  hydro_set_physical_internal_energy_dt(p, cosmo, cooling_du_dt);
+}
+
+/**
+ * @brief Computes the cooling time-step.
+ *
+ * The time-step is not set by the properties of cooling.
+ *
+ * @param cooling The #cooling_function_data used in the run.
+ * @param phys_const #phys_const data struct.
+ * @param us The internal system of units.
+ * @param cosmo #cosmology struct.
+ * @param hydro_props the properties of the hydro scheme.
+ * @param p #part data.
+ * @param xp extended particle data.
+ */
+__attribute__((always_inline)) INLINE float cooling_timestep(
+    const struct cooling_function_data *restrict cooling,
+    const struct phys_const *restrict phys_const,
+    const struct cosmology *restrict cosmo,
+    const struct unit_system *restrict us,
+    const struct hydro_props *hydro_props, const struct part *restrict p,
+    const struct xpart *restrict xp) {
+
+  return FLT_MAX;
+}
+
+/**
+ * @brief Sets the cooling properties of the (x-)particles to a valid start
+ * state.
+ *
+ * @param phys_const #phys_const data structure.
+ * @param us The internal system of units.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param cosmo #cosmology data structure.
+ * @param cooling #cooling_function_data struct.
+ * @param p #part data.
+ * @param xp Pointer to the #xpart data.
+ */
+__attribute__((always_inline)) INLINE void cooling_first_init_part(
+    const struct phys_const *restrict phys_const,
+    const struct unit_system *restrict us,
+    const struct hydro_props *hydro_props,
+    const struct cosmology *restrict cosmo,
+    const struct cooling_function_data *restrict cooling,
+    const struct part *restrict p, struct xpart *restrict xp) {}
+
+/**
+ * @brief Compute the temperature of a #part based on the cooling function.
+ *
+ * We use the Temperature table of the Wiersma+08 set. This computes the
+ * equilibirum temperature of a gas for a given redshift, Hydrogen density,
+ * internal energy per unit mass and Helium fraction.
+ *
+ * The temperature returned is consistent with the cooling rates.
+ *
+ * @param phys_const #phys_const data structure.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param us The internal system of units.
+ * @param cosmo #cosmology data structure.
+ * @param cooling #cooling_function_data struct.
+ * @param p #part data.
+ * @param xp Pointer to the #xpart data.
+ */
+float cooling_get_temperature(
+    const struct phys_const *restrict phys_const,
+    const struct hydro_props *restrict hydro_props,
+    const struct unit_system *restrict us,
+    const struct cosmology *restrict cosmo,
+    const struct cooling_function_data *restrict cooling,
+    const struct part *restrict p, const struct xpart *restrict xp) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (cooling->Redshifts == NULL)
+    error(
+        "Cooling function has not been initialised. Did you forget the "
+        "--temperature runtime flag?");
+#endif
+
+  /* Get physical internal energy */
+  const float u = hydro_get_physical_internal_energy(p, xp, cosmo);
+  const double u_cgs = u * cooling->internal_energy_to_cgs;
+
+  /* Get the Hydrogen and Helium mass fractions */
+  const float XH = 1. - phys_const->const_primordial_He_fraction;
+  const float XHe = phys_const->const_primordial_He_fraction;
+
+  /* Get the Helium mass fraction. Note that this is He / (H + He), i.e. a
+   * metal-free Helium mass fraction as per the Wiersma+08 definition */
+  const float HeFrac = XHe / (XH + XHe);
+
+  /* Convert Hydrogen mass fraction into Hydrogen number density */
+  const float rho = hydro_get_physical_density(p, cosmo);
+  const double n_H = rho * XH / phys_const->const_proton_mass;
+  const double n_H_cgs = n_H * cooling->number_density_to_cgs;
+
+  /* compute hydrogen number density and helium fraction table indices and
+   * offsets */
+  int He_index, n_H_index;
+  float d_He, d_n_H;
+  get_index_1d(cooling->HeFrac, qla_cooling_N_He_frac, HeFrac, &He_index,
+               &d_He);
+  get_index_1d(cooling->nH, qla_cooling_N_density, log10(n_H_cgs), &n_H_index,
+               &d_n_H);
+
+  /* Compute the log10 of the temperature by interpolating the table */
+  const double log_10_T = qla_convert_u_to_temp(
+      log10(u_cgs), cosmo->z, n_H_index, He_index, d_n_H, d_He, cooling);
+
+  /* Undo the log! */
+  return exp10(log_10_T);
+}
+
+/**
+ * @brief Returns the total radiated energy by this particle.
+ *
+ * @param xp #xpart data struct
+ */
+__attribute__((always_inline)) INLINE float cooling_get_radiated_energy(
+    const struct xpart *restrict xp) {
+
+  return 0.f;
+}
+
+/**
+ * @brief Split the coolong content of a particle into n pieces
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param n The number of pieces to split into.
+ */
+void cooling_split_part(struct part *p, struct xpart *xp, double n) {}
+
+/**
+ * @brief Inject a fixed amount of energy to each particle in the simulation
+ * to mimic Hydrogen reionization.
+ *
+ * @param cooling The properties of the cooling model.
+ * @param cosmo The cosmological model.
+ * @param s The #space containing the particles.
+ */
+void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling,
+                                   const struct cosmology *cosmo,
+                                   struct space *s) {
+
+  struct part *parts = s->parts;
+  struct xpart *xparts = s->xparts;
+
+  /* Energy to inject in internal units */
+  const float extra_heat =
+      cooling->H_reion_heat_cgs * cooling->internal_energy_from_cgs;
+
+  message("Applying extra energy for H reionization!");
+
+  /* Loop through particles and set new heat */
+  for (size_t i = 0; i < s->nr_parts; i++) {
+
+    struct part *p = &parts[i];
+    struct xpart *xp = &xparts[i];
+
+    const float old_u = hydro_get_physical_internal_energy(p, xp, cosmo);
+    const float new_u = old_u + extra_heat;
+
+    hydro_set_physical_internal_energy(p, xp, cosmo, new_u);
+    hydro_set_drifted_physical_internal_energy(p, cosmo, new_u);
+  }
+}
+
+/**
+ * @brief Initialises properties stored in the cooling_function_data struct
+ *
+ * @param parameter_file The parsed parameter file
+ * @param us Internal system of units data structure
+ * @param phys_const #phys_const data structure
+ * @param hydro_props The properties of the hydro scheme.
+ * @param cooling #cooling_function_data struct to initialize
+ */
+void cooling_init_backend(struct swift_params *parameter_file,
+                          const struct unit_system *us,
+                          const struct phys_const *phys_const,
+                          const struct hydro_props *hydro_props,
+                          struct cooling_function_data *cooling) {
+
+  /* Read model parameters */
+
+  /* Directory for cooling tables */
+  parser_get_param_string(parameter_file, "QLACooling:dir_name",
+                          cooling->cooling_table_path);
+
+  /* Despite the names, the values of H_reion_heat_cgs and He_reion_heat_cgs
+   * that are read in are actually in units of electron volts per proton mass.
+   * We later convert to units just below */
+
+  cooling->H_reion_done = 0;
+  cooling->H_reion_z =
+      parser_get_param_float(parameter_file, "QLACooling:H_reion_z");
+  cooling->H_reion_heat_cgs =
+      parser_get_param_float(parameter_file, "QLACooling:H_reion_eV_p_H");
+  cooling->He_reion_z_centre =
+      parser_get_param_float(parameter_file, "QLACooling:He_reion_z_centre");
+  cooling->He_reion_z_sigma =
+      parser_get_param_float(parameter_file, "QLACooling:He_reion_z_sigma");
+  cooling->He_reion_heat_cgs =
+      parser_get_param_float(parameter_file, "QLACooling:He_reion_eV_p_H");
+
+  /* Optional parameters to correct the abundances */
+  cooling->Ca_over_Si_ratio_in_solar = parser_get_opt_param_float(
+      parameter_file, "QLACooling:Ca_over_Si_in_solar", 1.f);
+  cooling->S_over_Si_ratio_in_solar = parser_get_opt_param_float(
+      parameter_file, "QLACooling:S_over_Si_in_solar", 1.f);
+
+  /* Convert H_reion_heat_cgs and He_reion_heat_cgs to cgs
+   * (units used internally by the cooling routines). This is done by
+   * multiplying by 'eV/m_H' in internal units, then converting to cgs units.
+   * Note that the dimensions of these quantities are energy/mass = velocity^2
+   */
+
+  cooling->H_reion_heat_cgs *=
+      phys_const->const_electron_volt / phys_const->const_proton_mass *
+      units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+
+  cooling->He_reion_heat_cgs *=
+      phys_const->const_electron_volt / phys_const->const_proton_mass *
+      units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+
+  /* Read in the list of redshifts */
+  get_cooling_redshifts(cooling);
+
+  /* Read in cooling table header */
+  char fname[qla_table_path_name_length + 12];
+  sprintf(fname, "%sz_0.000.hdf5", cooling->cooling_table_path);
+  read_cooling_header(fname, cooling);
+
+  /* Allocate space for cooling tables */
+  allocate_cooling_tables(cooling);
+
+  /* Compute conversion factors */
+  cooling->internal_energy_to_cgs =
+      units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+  cooling->internal_energy_from_cgs = 1. / cooling->internal_energy_to_cgs;
+  cooling->number_density_to_cgs =
+      units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY);
+
+  /* Store some constants in CGS units */
+  const double proton_mass_cgs =
+      phys_const->const_proton_mass *
+      units_cgs_conversion_factor(us, UNIT_CONV_MASS);
+  cooling->inv_proton_mass_cgs = 1. / proton_mass_cgs;
+  cooling->T_CMB_0 = phys_const->const_T_CMB_0 *
+                     units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
+
+  /* Compute the coefficient at the front of the Compton cooling expression */
+  const double radiation_constant =
+      4. * phys_const->const_stefan_boltzmann / phys_const->const_speed_light_c;
+  const double compton_coefficient =
+      4. * radiation_constant * phys_const->const_thomson_cross_section *
+      phys_const->const_boltzmann_k /
+      (phys_const->const_electron_mass * phys_const->const_speed_light_c);
+  const float dimension_coefficient[5] = {1, 2, -3, 0, -5};
+
+  /* This should be ~1.0178085e-37 g cm^2 s^-3 K^-5 */
+  const double compton_coefficient_cgs =
+      compton_coefficient *
+      units_general_cgs_conversion_factor(us, dimension_coefficient);
+
+#ifdef SWIFT_DEBUG_CHECKS
+  const double expected_compton_coefficient_cgs = 1.0178085e-37;
+  if (fabs(compton_coefficient_cgs - expected_compton_coefficient_cgs) /
+          expected_compton_coefficient_cgs >
+      0.01)
+    error("compton coefficient incorrect.");
+#endif
+
+  /* And now the Compton rate */
+  cooling->compton_rate_cgs = compton_coefficient_cgs * cooling->T_CMB_0 *
+                              cooling->T_CMB_0 * cooling->T_CMB_0 *
+                              cooling->T_CMB_0;
+
+  /* Set the redshift indices to invalid values */
+  cooling->z_index = -10;
+
+  /* set previous_z_index and to last value of redshift table*/
+  cooling->previous_z_index = qla_cooling_N_redshifts - 2;
+}
+
+/**
+ * @brief Restore cooling tables (if applicable) after
+ * restart
+ *
+ * @param cooling the #cooling_function_data structure
+ * @param cosmo #cosmology structure
+ */
+void cooling_restore_tables(struct cooling_function_data *cooling,
+                            const struct cosmology *cosmo) {
+
+  /* Read redshifts */
+  get_cooling_redshifts(cooling);
+
+  /* Read cooling header */
+  char fname[qla_table_path_name_length + 12];
+  sprintf(fname, "%sz_0.000.hdf5", cooling->cooling_table_path);
+  read_cooling_header(fname, cooling);
+
+  /* Allocate memory for the tables */
+  allocate_cooling_tables(cooling);
+
+  /* Force a re-read of the cooling tables */
+  cooling->z_index = -10;
+  cooling->previous_z_index = qla_cooling_N_redshifts - 2;
+  cooling_update(cosmo, cooling, /*space=*/NULL);
+}
+
+/**
+ * @brief Prints the properties of the cooling model to stdout.
+ *
+ * @param cooling #cooling_function_data struct.
+ */
+void cooling_print_backend(const struct cooling_function_data *cooling) {
+
+  message("Cooling function is 'Quick Lyman-alpha (EAGLE with primordial Z)'.");
+}
+
+/**
+ * @brief Clean-up the memory allocated for the cooling routines
+ *
+ * We simply free all the arrays.
+ *
+ * @param cooling the cooling data structure.
+ */
+void cooling_clean(struct cooling_function_data *cooling) {
+
+  /* Free the side arrays */
+  swift_free("cooling", cooling->Redshifts);
+  swift_free("cooling", cooling->nH);
+  swift_free("cooling", cooling->Temp);
+  swift_free("cooling", cooling->HeFrac);
+  swift_free("cooling", cooling->Therm);
+  swift_free("cooling", cooling->SolarAbundances);
+  swift_free("cooling", cooling->SolarAbundances_inv);
+
+  /* Free the tables */
+  swift_free("cooling-tables", cooling->table.metal_heating);
+  swift_free("cooling-tables", cooling->table.electron_abundance);
+  swift_free("cooling-tables", cooling->table.temperature);
+  swift_free("cooling-tables", cooling->table.H_plus_He_heating);
+  swift_free("cooling-tables", cooling->table.H_plus_He_electron_abundance);
+}
+
+/**
+ * @brief Write a cooling struct to the given FILE as a stream of bytes.
+ *
+ * @param cooling the struct
+ * @param stream the file stream
+ */
+void cooling_struct_dump(const struct cooling_function_data *cooling,
+                         FILE *stream) {
+
+  /* To make sure everything is restored correctly, we zero all the pointers to
+     tables. If they are not restored correctly, we would crash after restart on
+     the first call to the cooling routines. Helps debugging. */
+  struct cooling_function_data cooling_copy = *cooling;
+  cooling_copy.Redshifts = NULL;
+  cooling_copy.nH = NULL;
+  cooling_copy.Temp = NULL;
+  cooling_copy.Therm = NULL;
+  cooling_copy.SolarAbundances = NULL;
+  cooling_copy.SolarAbundances_inv = NULL;
+  cooling_copy.table.metal_heating = NULL;
+  cooling_copy.table.H_plus_He_heating = NULL;
+  cooling_copy.table.H_plus_He_electron_abundance = NULL;
+  cooling_copy.table.temperature = NULL;
+  cooling_copy.table.electron_abundance = NULL;
+
+  restart_write_blocks((void *)&cooling_copy,
+                       sizeof(struct cooling_function_data), 1, stream,
+                       "cooling", "cooling function");
+}
+
+/**
+ * @brief Restore a hydro_props struct from the given FILE as a stream of
+ * bytes.
+ *
+ * Read the structure from the stream and restore the cooling tables by
+ * re-reading them.
+ *
+ * @param cooling the struct
+ * @param stream the file stream
+ * @param cosmo #cosmology structure
+ */
+void cooling_struct_restore(struct cooling_function_data *cooling, FILE *stream,
+                            const struct cosmology *cosmo) {
+  restart_read_blocks((void *)cooling, sizeof(struct cooling_function_data), 1,
+                      stream, NULL, "cooling function");
+
+  cooling_restore_tables(cooling, cosmo);
+}
diff --git a/src/cooling/QLA/cooling.h b/src/cooling/QLA/cooling.h
new file mode 100644
index 0000000000000000000000000000000000000000..89680e3473b2acbf03267a213593c8925aa04a9e
--- /dev/null
+++ b/src/cooling/QLA/cooling.h
@@ -0,0 +1,91 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_COOLING_QLA_H
+#define SWIFT_COOLING_QLA_H
+
+/**
+ * @file src/cooling/EAGLE/cooling.h
+ * @brief EAGLE cooling function declarations
+ */
+
+/* Local includes. */
+#include "cooling_struct.h"
+
+struct part;
+struct xpart;
+struct cosmology;
+struct hydro_props;
+struct entropy_floor_properties;
+struct space;
+
+void cooling_update(const struct cosmology *cosmo,
+                    struct cooling_function_data *cooling, struct space *s);
+
+void cooling_cool_part(const struct phys_const *phys_const,
+                       const struct unit_system *us,
+                       const struct cosmology *cosmo,
+                       const struct hydro_props *hydro_properties,
+                       const struct entropy_floor_properties *floor_props,
+                       const struct cooling_function_data *cooling,
+                       struct part *restrict p, struct xpart *restrict xp,
+                       const float dt, const float dt_therm, const double time);
+
+float cooling_timestep(const struct cooling_function_data *restrict cooling,
+                       const struct phys_const *restrict phys_const,
+                       const struct cosmology *restrict cosmo,
+                       const struct unit_system *restrict us,
+                       const struct hydro_props *hydro_props,
+                       const struct part *restrict p,
+                       const struct xpart *restrict xp);
+
+void cooling_first_init_part(
+    const struct phys_const *restrict phys_const,
+    const struct unit_system *restrict us,
+    const struct hydro_props *hydro_props,
+    const struct cosmology *restrict cosmo,
+    const struct cooling_function_data *restrict cooling,
+    const struct part *restrict p, struct xpart *restrict xp);
+
+float cooling_get_temperature(
+    const struct phys_const *restrict phys_const,
+    const struct hydro_props *restrict hydro_props,
+    const struct unit_system *restrict us,
+    const struct cosmology *restrict cosmo,
+    const struct cooling_function_data *restrict cooling,
+    const struct part *restrict p, const struct xpart *restrict xp);
+
+float cooling_get_radiated_energy(const struct xpart *restrict xp);
+
+void cooling_split_part(struct part *p, struct xpart *xp, double n);
+
+void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling,
+                                   const struct cosmology *cosmo,
+                                   struct space *s);
+
+void cooling_init_backend(struct swift_params *parameter_file,
+                          const struct unit_system *us,
+                          const struct phys_const *phys_const,
+                          const struct hydro_props *hydro_props,
+                          struct cooling_function_data *cooling);
+
+void cooling_print_backend(const struct cooling_function_data *cooling);
+
+void cooling_clean(struct cooling_function_data *data);
+
+#endif /* SWIFT_COOLING_QLA_H */
diff --git a/src/cooling/QLA/cooling_io.h b/src/cooling/QLA/cooling_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..71303e5e4277159afd69037ac11954dbc7f7d0c9
--- /dev/null
+++ b/src/cooling/QLA/cooling_io.h
@@ -0,0 +1,73 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_COOLING_QLA_IO_H
+#define SWIFT_COOLING_QLA_IO_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local includes */
+#include "cooling.h"
+#include "io_properties.h"
+
+#ifdef HAVE_HDF5
+
+/**
+ * @brief Writes the current model of cooling to the file.
+ *
+ * @param h_grp The HDF5 group in which to write
+ * @param cooling The #cooling_function_data
+ */
+__attribute__((always_inline)) INLINE static void cooling_write_flavour(
+    hid_t h_grp, const struct cooling_function_data* cooling) {
+
+  io_write_attribute_s(h_grp, "Cooling Model",
+                       "Quick Lyman-alpha (EAGLE with primordial Z only)");
+}
+#endif
+
+INLINE static void convert_part_T(const struct engine* e, const struct part* p,
+                                  const struct xpart* xp, float* ret) {
+
+  ret[0] = cooling_get_temperature(e->physical_constants, e->hydro_properties,
+                                   e->internal_units, e->cosmology,
+                                   e->cooling_func, p, xp);
+}
+
+/**
+ * @brief Specifies which particle fields to write to a dataset
+ *
+ * @param parts The particle array.
+ * @param xparts The extended data particle array.
+ * @param list The list of i/o properties to write.
+ * @param cooling The #cooling_function_data
+ *
+ * @return Returns the number of fields to write.
+ */
+__attribute__((always_inline)) INLINE static int cooling_write_particles(
+    const struct part* parts, const struct xpart* xparts, struct io_props* list,
+    const struct cooling_function_data* cooling) {
+
+  list[0] = io_make_output_field_convert_part(
+      "Temperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, parts, xparts,
+      convert_part_T, "Temperatures of the gas particles");
+  return 1;
+}
+
+#endif /* SWIFT_COOLING_QLA_IO_H */
diff --git a/src/cooling/QLA/cooling_rates.h b/src/cooling/QLA/cooling_rates.h
new file mode 100644
index 0000000000000000000000000000000000000000..c193fca5b8e803f55b8bbc1a0fbd42ecf2f0911f
--- /dev/null
+++ b/src/cooling/QLA/cooling_rates.h
@@ -0,0 +1,519 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#ifndef SWIFT_QLA_COOLING_RATES_H
+#define SWIFT_QLA_COOLING_RATES_H
+
+#include "../config.h"
+
+/* Local includes. */
+#include "cooling_tables.h"
+#include "exp10.h"
+#include "interpolate.h"
+
+/**
+ * @brief Compute ratio of mass fraction to solar mass fraction
+ * for each element carried by a given particle.
+ *
+ * The solar abundances are taken from the tables themselves.
+ *
+ * The quick Lyman-alpha chemistry model does not track any elements.
+ * We use the primordial H and He to fill in the element abundance
+ * array. 
+ *
+ * We also re-order the elements such that they match the order of the
+ * tables. This is [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe].
+ * 
+ * The solar abundances table (from the cooling struct) is arranged as
+ * [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe].
+ *
+ * @param p Pointer to #part struct.
+ * @param cooling #cooling_function_data struct.
+ * @param phys_const the physical constants in internal units 
+ * @param ratio_solar (return) Array of ratios to solar abundances.
+ */
+__attribute__((always_inline)) INLINE static void abundance_ratio_to_solar(
+    const struct part *p, const struct cooling_function_data *cooling,
+    const struct phys_const *phys_const,
+    float ratio_solar[qla_cooling_N_abundances]) {
+
+  /* We just construct an array of primoridal abundances as fractions
+     of solar */
+
+  ratio_solar[0] = (1. - phys_const->const_primordial_He_fraction) *
+                   cooling->SolarAbundances_inv[0 /* H */];
+
+  ratio_solar[1] = phys_const->const_primordial_He_fraction *
+                   cooling->SolarAbundances_inv[1 /* He */];
+
+  ratio_solar[2] = 0.f;
+  ratio_solar[3] = 0.f;
+  ratio_solar[4] = 0.f;
+  ratio_solar[5] = 0.f;
+  ratio_solar[6] = 0.f;
+  ratio_solar[7] = 0.f;
+  ratio_solar[8] = 0.f;
+  ratio_solar[9] = 0.f;
+  ratio_solar[10] = 0.f;
+}
+
+/**
+ * @brief Computes the extra heat from Helium reionisation at a given redshift.
+ *
+ * We follow the implementation of Wiersma et al. 2009, MNRAS, 399, 574-600,
+ * section. 2. The calculation returns energy in CGS.
+ *
+ * Note that delta_z is negative.
+ *
+ * @param z The current redshift.
+ * @param delta_z The change in redhsift over the course of this time-step.
+ * @param cooling The #cooling_function_data used in the run.
+ * @return Helium reionization energy in CGS units.
+ */
+__attribute__((always_inline)) INLINE double qla_helium_reionization_extraheat(
+    const double z, const double delta_z,
+    const struct cooling_function_data *cooling) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (delta_z > 0.f) error("Invalid value for delta_z. Should be negative.");
+#endif
+
+  /* Recover the values we need */
+  const double z_centre = cooling->He_reion_z_centre;
+  const double z_sigma = cooling->He_reion_z_sigma;
+  const double heat_cgs = cooling->He_reion_heat_cgs;
+
+  double extra_heat = 0.;
+
+  /* Integral of the Gaussian between z and z - delta_z */
+  extra_heat += erf((z - delta_z - z_centre) / (M_SQRT2 * z_sigma));
+  extra_heat -= erf((z - z_centre) / (M_SQRT2 * z_sigma));
+
+  /* Multiply by the normalisation factor */
+  extra_heat *= heat_cgs * 0.5;
+
+  return extra_heat;
+}
+
+/**
+ * @brief Computes the log_10 of the temperature corresponding to a given
+ * internal energy, hydrogen number density, Helium fraction and redshift.
+ *
+ * Note that the redshift is implicitly passed in via the currently loaded
+ * tables in the #cooling_function_data.
+ *
+ * For the low-z case, we interpolate the flattened 4D table 'u_to_temp' that
+ * is arranged in the following way:
+ * - 1st dim: redshift, length = eagle_cooling_N_loaded_redshifts
+ * - 2nd dim: Hydrogen density, length = eagle_cooling_N_density
+ * - 3rd dim: Helium fraction, length = eagle_cooling_N_He_frac
+ * - 4th dim: Internal energy, length = eagle_cooling_N_temperature
+ *
+ * For the high-z case, we interpolate the flattened 3D table 'u_to_temp' that
+ * is arranged in the following way:
+ * - 1st dim: Hydrogen density, length = eagle_cooling_N_density
+ * - 2nd dim: Helium fraction, length = eagle_cooling_N_He_frac
+ * - 3rd dim: Internal energy, length = eagle_cooling_N_temperature
+ *
+ * @param log_10_u_cgs Log base 10 of internal energy in cgs.
+ * @param redshift Current redshift.
+ * @param n_H_index Index along the Hydrogen density dimension.
+ * @param He_index Index along the Helium fraction dimension.
+ * @param d_n_H Offset between Hydrogen density and table[n_H_index].
+ * @param d_He Offset between helium fraction and table[He_index].
+ * @param cooling #cooling_function_data structure.
+ *
+ * @return log_10 of the temperature.
+ */
+__attribute__((always_inline)) INLINE double qla_convert_u_to_temp(
+    const double log_10_u_cgs, const float redshift, const int n_H_index,
+    const int He_index, const float d_n_H, const float d_He,
+    const struct cooling_function_data *cooling) {
+
+  /* Get index of u along the internal energy axis */
+  int u_index;
+  float d_u;
+  get_index_1d(cooling->Therm, qla_cooling_N_temperature, log_10_u_cgs,
+               &u_index, &d_u);
+
+  /* Interpolate temperature table to return temperature for current
+   * internal energy (use 3D interpolation for high redshift table,
+   * otherwise 4D) */
+  float log_10_T;
+  if (redshift > cooling->Redshifts[qla_cooling_N_redshifts - 1]) {
+
+    log_10_T = interpolation_3d(cooling->table.temperature,   /* */
+                                n_H_index, He_index, u_index, /* */
+                                d_n_H, d_He, d_u,             /* */
+                                qla_cooling_N_density,        /* */
+                                qla_cooling_N_He_frac,        /* */
+                                qla_cooling_N_temperature);   /* */
+  } else {
+
+    log_10_T =
+        interpolation_4d(cooling->table.temperature,                  /* */
+                         /*z_index=*/0, n_H_index, He_index, u_index, /* */
+                         cooling->dz, d_n_H, d_He, d_u,               /* */
+                         qla_cooling_N_loaded_redshifts,              /* */
+                         qla_cooling_N_density,                       /* */
+                         qla_cooling_N_He_frac,                       /* */
+                         qla_cooling_N_temperature);                  /* */
+  }
+
+  /* Special case for temperatures below the start of the table */
+  if (u_index == 0 && d_u == 0.f) {
+
+    /* The temperature is multiplied by u / 10^T[0]
+     * where T[0] is the first entry in the table */
+    log_10_T += log_10_u_cgs - cooling->Temp[0];
+  }
+
+  return log_10_T;
+}
+
+/**
+ * @brief Compute the Compton cooling rate from the CMB at a given
+ * redshift, electron abundance, temperature and Hydrogen density.
+ *
+ * Uses an analytic formula.
+ *
+ * @param cooling The #cooling_function_data used in the run.
+ * @param redshift The current redshift.
+ * @param n_H_cgs The Hydrogen number density in CGS units.
+ * @param temperature The temperature.
+ * @param electron_abundance The electron abundance.
+ */
+__attribute__((always_inline)) INLINE double qla_Compton_cooling_rate(
+    const struct cooling_function_data *cooling, const double redshift,
+    const double n_H_cgs, const double temperature,
+    const double electron_abundance) {
+
+  const double zp1 = 1. + redshift;
+  const double zp1p2 = zp1 * zp1;
+  const double zp1p4 = zp1p2 * zp1p2;
+
+  /* CMB temperature at this redshift */
+  const double T_CMB = cooling->T_CMB_0 * zp1;
+
+  /* Compton cooling rate */
+  return cooling->compton_rate_cgs * (temperature - T_CMB) * zp1p4 *
+         electron_abundance / n_H_cgs;
+}
+
+/**
+ * @brief Computes the cooling rate corresponding to a given internal energy,
+ * hydrogen number density, Helium fraction, redshift and metallicity from
+ * all the possible channels.
+ *
+ * 1) Metal-free cooling:
+ * We interpolate the flattened 4D table 'H_and_He_net_heating' that is
+ * arranged in the following way:
+ * - 1st dim: redshift, length = eagle_cooling_N_loaded_redshifts
+ * - 2nd dim: Hydrogen density, length = eagle_cooling_N_density
+ * - 3rd dim: Helium fraction, length = eagle_cooling_N_He_frac
+ * - 4th dim: Internal energy, length = eagle_cooling_N_temperature
+ *
+ * 2) Electron abundance
+ * We compute the electron abundance by interpolating the flattened 4d table
+ * 'H_and_He_electron_abundance' that is arranged in the following way:
+ * - 1st dim: redshift, length = eagle_cooling_N_loaded_redshifts
+ * - 2nd dim: Hydrogen density, length = eagle_cooling_N_density
+ * - 3rd dim: Helium fraction, length = eagle_cooling_N_He_frac
+ * - 4th dim: Internal energy, length = eagle_cooling_N_temperature
+ *
+ * 3) Compton cooling is applied via the analytic formula.
+ *
+ * 4) Solar electron abudance
+ * We compute the solar electron abundance by interpolating the flattened 3d
+ * table 'solar_electron_abundance' that is arranged in the following way:
+ * - 1st dim: redshift, length = eagle_cooling_N_loaded_redshifts
+ * - 2nd dim: Hydrogen density, length = eagle_cooling_N_density
+ * - 3rd dim: Internal energy, length = eagle_cooling_N_temperature
+ *
+ * 5) Metal-line cooling
+ * For each tracked element we interpolate the flattened 4D table
+ * 'table_metals_net_heating' that is arrange in the following way:
+ * - 1st dim: element, length = eagle_cooling_N_metal
+ * - 2nd dim: redshift, length = eagle_cooling_N_loaded_redshifts
+ * - 3rd dim: Hydrogen density, length = eagle_cooling_N_density
+ * - 4th dim: Internal energy, length = eagle_cooling_N_temperature
+ *
+ * Note that this is a fake 4D interpolation as we do not interpolate
+ * along the 1st dimension. We just do this once per element.
+ *
+ * Since only the temperature changes when cooling a given particle,
+ * the redshift, hydrogen number density and helium fraction indices
+ * and offsets passed in.
+ *
+ * If the argument element_lambda is non-NULL, the routine also
+ * returns the cooling rate per element in the array.
+ *
+ * @param log10_u_cgs Log base 10 of internal energy per unit mass in CGS units.
+ * @param redshift The current redshift
+ * @param n_H_cgs The Hydrogen number density in CGS units.
+ * @param solar_ratio Array of ratios of particle metal abundances
+ * to solar metal abundances
+ *
+ * @param n_H_index Particle hydrogen number density index
+ * @param d_n_H Particle hydrogen number density offset
+ * @param He_index Particle helium fraction index
+ * @param d_He Particle helium fraction offset
+ * @param cooling Cooling data structure
+ *
+ * @param element_lambda (return) Cooling rate from each element
+ *
+ * @return The cooling rate
+ */
+INLINE static double qla_metal_cooling_rate(
+    const double log10_u_cgs, const double redshift, const double n_H_cgs,
+    const float solar_ratio[qla_cooling_N_abundances], const int n_H_index,
+    const float d_n_H, const int He_index, const float d_He,
+    const struct cooling_function_data *cooling, double *element_lambda) {
+
+  /* Temperature */
+  const double log_10_T = qla_convert_u_to_temp(
+      log10_u_cgs, redshift, n_H_index, He_index, d_n_H, d_He, cooling);
+
+  /* Get index along temperature dimension of the tables */
+  int T_index;
+  float d_T;
+  get_index_1d(cooling->Temp, qla_cooling_N_temperature, log_10_T, &T_index,
+               &d_T);
+
+  /**********************/
+  /* Metal-free cooling */
+  /**********************/
+
+  double Lambda_free;
+
+  if (redshift > cooling->Redshifts[qla_cooling_N_redshifts - 1]) {
+
+    /* If we're using the high redshift tables then we don't interpolate
+     * in redshift */
+    Lambda_free = interpolation_3d(cooling->table.H_plus_He_heating, /* */
+                                   n_H_index, He_index, T_index,     /* */
+                                   d_n_H, d_He, d_T,                 /* */
+                                   qla_cooling_N_density,            /* */
+                                   qla_cooling_N_He_frac,            /* */
+                                   qla_cooling_N_temperature);       /* */
+  } else {
+
+    /* Using normal tables, have to interpolate in redshift */
+    Lambda_free =
+        interpolation_4d(cooling->table.H_plus_He_heating,            /* */
+                         /*z_index=*/0, n_H_index, He_index, T_index, /* */
+                         cooling->dz, d_n_H, d_He, d_T,               /* */
+                         qla_cooling_N_loaded_redshifts,              /* */
+                         qla_cooling_N_density,                       /* */
+                         qla_cooling_N_He_frac,                       /* */
+                         qla_cooling_N_temperature);                  /* */
+  }
+
+  /* If we're testing cooling rate contributions write to array */
+  if (element_lambda != NULL) {
+    element_lambda[0] = Lambda_free;
+  }
+
+  /**********************/
+  /* Electron abundance */
+  /**********************/
+
+  double H_plus_He_electron_abundance;
+
+  if (redshift > cooling->Redshifts[qla_cooling_N_redshifts - 1]) {
+
+    H_plus_He_electron_abundance =
+        interpolation_3d(cooling->table.H_plus_He_electron_abundance, /* */
+                         n_H_index, He_index, T_index,                /* */
+                         d_n_H, d_He, d_T,                            /* */
+                         qla_cooling_N_density,                       /* */
+                         qla_cooling_N_He_frac,                       /* */
+                         qla_cooling_N_temperature);                  /* */
+
+  } else {
+
+    H_plus_He_electron_abundance =
+        interpolation_4d(cooling->table.H_plus_He_electron_abundance, /* */
+                         /*z_index=*/0, n_H_index, He_index, T_index, /* */
+                         cooling->dz, d_n_H, d_He, d_T,               /* */
+                         qla_cooling_N_loaded_redshifts,              /* */
+                         qla_cooling_N_density,                       /* */
+                         qla_cooling_N_He_frac,                       /* */
+                         qla_cooling_N_temperature);                  /* */
+  }
+
+  /**********************/
+  /* Compton cooling    */
+  /**********************/
+
+  double Lambda_Compton = 0.;
+
+  /* Do we need to add the inverse Compton cooling? */
+  /* It is *not* stored in the tables before re-ionisation */
+  if ((redshift > cooling->Redshifts[qla_cooling_N_redshifts - 1]) ||
+      (redshift > cooling->H_reion_z)) {
+
+    const double T = exp10(log_10_T);
+
+    /* Note the minus sign */
+    Lambda_Compton -= qla_Compton_cooling_rate(cooling, redshift, n_H_cgs, T,
+                                               H_plus_He_electron_abundance);
+  }
+
+  /* If we're testing cooling rate contributions write to array */
+  if (element_lambda != NULL) {
+    element_lambda[1] = Lambda_Compton;
+  }
+
+  /*******************************/
+  /* Solar electron abundance    */
+  /*******************************/
+
+  double solar_electron_abundance;
+
+  if (redshift > cooling->Redshifts[qla_cooling_N_redshifts - 1]) {
+
+    /* If we're using the high redshift tables then we don't interpolate
+     * in redshift */
+    solar_electron_abundance =
+        interpolation_2d(cooling->table.electron_abundance, /* */
+                         n_H_index, T_index,                /* */
+                         d_n_H, d_T,                        /* */
+                         qla_cooling_N_density,             /* */
+                         qla_cooling_N_temperature);        /* */
+
+  } else {
+
+    /* Using normal tables, have to interpolate in redshift */
+    solar_electron_abundance =
+        interpolation_3d(cooling->table.electron_abundance, /* */
+                         /*z_index=*/0, n_H_index, T_index, /* */
+                         cooling->dz, d_n_H, d_T,           /* */
+                         qla_cooling_N_loaded_redshifts,    /* */
+                         qla_cooling_N_density,             /* */
+                         qla_cooling_N_temperature);        /* */
+  }
+
+  const double electron_abundance_ratio =
+      H_plus_He_electron_abundance / solar_electron_abundance;
+
+  /**********************/
+  /* Metal-line cooling */
+  /**********************/
+
+  /* for each element the cooling rate is multiplied by the ratio of H, He
+   * electron abundance to solar electron abundance then by the ratio of the
+   * particle metal abundance to solar metal abundance. */
+
+  double lambda_metal[qla_cooling_N_metal + 2] = {0.};
+
+  if (redshift > cooling->Redshifts[qla_cooling_N_redshifts - 1]) {
+
+    /* Loop over the metals (ignore H and He) */
+    for (int elem = 2; elem < qla_cooling_N_metal + 2; elem++) {
+
+      if (solar_ratio[elem] > 0.) {
+
+        /* Note that we do not interpolate along the x-axis
+         * (element dimension) */
+        lambda_metal[elem] =
+            interpolation_3d_no_x(cooling->table.metal_heating,   /* */
+                                  elem - 2, n_H_index, T_index,   /* */
+                                  /*delta_elem=*/0.f, d_n_H, d_T, /* */
+                                  qla_cooling_N_metal,            /* */
+                                  qla_cooling_N_density,          /* */
+                                  qla_cooling_N_temperature);     /* */
+
+        lambda_metal[elem] *= electron_abundance_ratio;
+        lambda_metal[elem] *= solar_ratio[elem];
+      }
+    }
+
+  } else {
+
+    /* Loop over the metals (ignore H and He) */
+    for (int elem = 2; elem < qla_cooling_N_metal + 2; elem++) {
+
+      if (solar_ratio[elem] > 0.) {
+
+        /* Note that we do not interpolate along the x-axis
+         * (element dimension) */
+        lambda_metal[elem] = interpolation_4d_no_x(
+            cooling->table.metal_heating,                /* */
+            elem - 2, /*z_index=*/0, n_H_index, T_index, /* */
+            /*delta_elem=*/0.f, cooling->dz, d_n_H, d_T, /* */
+            qla_cooling_N_metal,                         /* */
+            qla_cooling_N_loaded_redshifts,              /* */
+            qla_cooling_N_density,                       /* */
+            qla_cooling_N_temperature);                  /* */
+
+        lambda_metal[elem] *= electron_abundance_ratio;
+        lambda_metal[elem] *= solar_ratio[elem];
+
+        // message("lambda[%d]=%e", elem, lambda_metal[elem]);
+      }
+    }
+  }
+
+  if (element_lambda != NULL) {
+    for (int elem = 2; elem < qla_cooling_N_metal + 2; ++elem) {
+      element_lambda[elem] = lambda_metal[elem];
+    }
+  }
+
+  /* Sum up all the contributions */
+  double Lambda_net = Lambda_free + Lambda_Compton;
+  for (int elem = 2; elem < qla_cooling_N_metal + 2; ++elem) {
+    Lambda_net += lambda_metal[elem];
+  }
+
+  return Lambda_net;
+}
+
+/**
+ * @brief Wrapper function used to calculate cooling rate.
+ * Table indices and offsets for redshift, hydrogen number density and
+ * helium fraction are passed it so as to compute them only once per particle.
+ *
+ * @param log10_u_cgs Log base 10 of internal energy per unit mass in CGS units.
+ * @param redshift The current redshift.
+ * @param n_H_cgs Hydrogen number density in CGS units.
+ * @param abundance_ratio Ratio of element abundance to solar.
+ *
+ * @param n_H_index Particle hydrogen number density index
+ * @param d_n_H Particle hydrogen number density offset
+ * @param He_index Particle helium fraction index
+ * @param d_He Particle helium fraction offset
+ * @param cooling #cooling_function_data structure
+ *
+ * @return The cooling rate
+ */
+INLINE static double qla_cooling_rate(
+    const double log10_u_cgs, const double redshift, const double n_H_cgs,
+    const float abundance_ratio[qla_cooling_N_abundances], const int n_H_index,
+    const float d_n_H, const int He_index, const float d_He,
+    const struct cooling_function_data *cooling) {
+
+  return qla_metal_cooling_rate(log10_u_cgs, redshift, n_H_cgs, abundance_ratio,
+                                n_H_index, d_n_H, He_index, d_He, cooling,
+                                /* element_lambda=*/NULL);
+}
+
+#endif /* SWIFT_QLA_COOLING_RATES_H */
diff --git a/src/cooling/QLA/cooling_struct.h b/src/cooling/QLA/cooling_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..eab843b7f1ed416b3c1d9eb7b5738a15332d42c8
--- /dev/null
+++ b/src/cooling/QLA/cooling_struct.h
@@ -0,0 +1,136 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_COOLING_STRUCT_QLA_H
+#define SWIFT_COOLING_STRUCT_QLA_H
+
+#define qla_table_path_name_length 500
+
+/**
+ * @brief struct containing cooling tables
+ */
+struct cooling_tables {
+
+  /* array of heating rates due to metals */
+  float *metal_heating;
+
+  /* array of heating rates due to hydrogen and helium */
+  float *H_plus_He_heating;
+
+  /* array of electron abundances due to hydrogen and helium */
+  float *H_plus_He_electron_abundance;
+
+  /* array of temperatures */
+  float *temperature;
+
+  /* array of electron abundances due to metals */
+  float *electron_abundance;
+};
+
+/**
+ * @brief Properties of the cooling function.
+ */
+struct cooling_function_data {
+
+  /*! Cooling tables */
+  struct cooling_tables table;
+
+  /*! Redshift bins */
+  float *Redshifts;
+
+  /*! Hydrogen number density bins */
+  float *nH;
+
+  /*! Temperature bins */
+  float *Temp;
+
+  /*! Helium fraction bins */
+  float *HeFrac;
+
+  /*! Internal energy bins */
+  float *Therm;
+
+  /*! Mass fractions of elements for solar abundances (from the tables) */
+  float *SolarAbundances;
+
+  /*! Inverse of the solar mass fractions */
+  float *SolarAbundances_inv;
+
+  /*! Filepath to the directory containing the HDF5 cooling tables */
+  char cooling_table_path[qla_table_path_name_length];
+
+  /*! Redshift of H reionization */
+  float H_reion_z;
+
+  /*! H reionization energy in CGS units */
+  float H_reion_heat_cgs;
+
+  /*! Have we already done H reioisation? */
+  int H_reion_done;
+
+  /*! Ca over Si abundance divided by the solar ratio for these elements */
+  float Ca_over_Si_ratio_in_solar;
+
+  /*! S over Si abundance divided by the solar ratio for these elements */
+  float S_over_Si_ratio_in_solar;
+
+  /*! Redshift of He reionization */
+  float He_reion_z_centre;
+
+  /*! Spread of the He reionization */
+  float He_reion_z_sigma;
+
+  /*! He reionization energy in CGS units */
+  float He_reion_heat_cgs;
+
+  /*! Internal energy conversion from internal units to CGS (for quick access)
+   */
+  double internal_energy_to_cgs;
+
+  /*! Internal energy conversion from CGS to internal units (for quick access)
+   */
+  double internal_energy_from_cgs;
+
+  /*! Number density conversion from internal units to CGS (for quick access) */
+  double number_density_to_cgs;
+
+  /*! Inverse of proton mass in cgs (for quick access) */
+  double inv_proton_mass_cgs;
+
+  /*! Temperature of the CMB at present day (for quick access) */
+  double T_CMB_0;
+
+  /*! Compton rate in cgs units */
+  double compton_rate_cgs;
+
+  /*! Index of the current redshift along the redshift index of the tables */
+  int z_index;
+
+  /*! Distance between the current redshift and table[z_index] */
+  float dz;
+
+  /*! Index of the previous tables along the redshift index of the tables */
+  int previous_z_index;
+};
+
+/**
+ * @brief Properties of the cooling stored in the extended particle data.
+ */
+struct cooling_xpart_data {};
+
+#endif /* SWIFT_COOLING_STRUCT_QLA_H */
diff --git a/src/cooling/QLA/cooling_tables.c b/src/cooling/QLA/cooling_tables.c
new file mode 100644
index 0000000000000000000000000000000000000000..f2159e9540c11ab3520a6779a03ef588718b3c63
--- /dev/null
+++ b/src/cooling/QLA/cooling_tables.c
@@ -0,0 +1,768 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/**
+ * @file src/cooling/QLA/cooling_tables.c
+ * @brief Functions to read QLA tables
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+#include <hdf5.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+/* Local includes. */
+#include "chemistry_struct.h"
+#include "cooling_struct.h"
+#include "cooling_tables.h"
+#include "error.h"
+#include "interpolate.h"
+
+/**
+ * @brief Names of the elements in the order they are stored in the files
+ */
+static const char *qla_tables_element_names[qla_cooling_N_metal] = {
+    "Carbon",  "Nitrogen", "Oxygen",  "Neon", "Magnesium",
+    "Silicon", "Sulphur",  "Calcium", "Iron"};
+
+/*! Number of elements in a z-slice of the H+He cooling rate tables */
+static const size_t num_elements_cooling_rate =
+    qla_cooling_N_temperature * qla_cooling_N_density;
+
+/*! Number of elements in a z-slice of the metal cooling rate tables */
+static const size_t num_elements_metal_heating =
+    qla_cooling_N_metal * qla_cooling_N_temperature * qla_cooling_N_density;
+
+/*! Number of elements in a z-slice of the metal electron abundance tables */
+static const size_t num_elements_electron_abundance =
+    qla_cooling_N_temperature * qla_cooling_N_density;
+
+/*! Number of elements in a z-slice of the temperature tables */
+static const size_t num_elements_temperature =
+    qla_cooling_N_He_frac * qla_cooling_N_temperature * qla_cooling_N_density;
+
+/*! Number of elements in a z-slice of the H+He cooling rate tables */
+static const size_t num_elements_HpHe_heating =
+    qla_cooling_N_He_frac * qla_cooling_N_temperature * qla_cooling_N_density;
+
+/*! Number of elements in a z-slice of the H+He electron abundance tables */
+static const size_t num_elements_HpHe_electron_abundance =
+    qla_cooling_N_He_frac * qla_cooling_N_temperature * qla_cooling_N_density;
+
+/**
+ * @brief Reads in EAGLE table of redshift values
+ *
+ * @param cooling #cooling_function_data structure
+ */
+void get_cooling_redshifts(struct cooling_function_data *cooling) {
+
+  /* Read the list of table redshifts */
+  char redshift_filename[qla_table_path_name_length + 16];
+  sprintf(redshift_filename, "%s/redshifts.dat", cooling->cooling_table_path);
+
+  FILE *infile = fopen(redshift_filename, "r");
+  if (infile == NULL) {
+    error("Cannot open the list of cooling table redshifts (%s)",
+          redshift_filename);
+  }
+
+  int N_Redshifts = -1;
+
+  /* Read the file */
+  if (!feof(infile)) {
+
+    char buffer[50];
+
+    /* Read the number of redshifts (1st line in the file) */
+    if (fgets(buffer, 50, infile) != NULL)
+      N_Redshifts = atoi(buffer);
+    else
+      error("Impossible to read the number of redshifts");
+
+    /* Be verbose about it */
+    message("Found cooling tables at %d redhsifts", N_Redshifts);
+
+    /* Check value */
+    if (N_Redshifts != qla_cooling_N_redshifts)
+      error("Invalid redshift length array.");
+
+    /* Allocate the list of redshifts */
+    if (swift_memalign("cooling", (void **)&cooling->Redshifts,
+                       SWIFT_STRUCT_ALIGNMENT,
+                       qla_cooling_N_redshifts * sizeof(float)) != 0)
+      error("Failed to allocate redshift table");
+
+    /* Read all the redshift values */
+    int count = 0;
+    while (!feof(infile)) {
+      if (fgets(buffer, 50, infile) != NULL) {
+        cooling->Redshifts[count] = atof(buffer);
+        count++;
+      }
+    }
+
+    /* Verify that the file was self-consistent */
+    if (count != N_Redshifts) {
+      error(
+          "Redshift file (%s) does not contain the correct number of redshifts "
+          "(%d vs. %d)",
+          redshift_filename, count, N_Redshifts);
+    }
+  } else {
+    error("Redshift file (%s) is empty!", redshift_filename);
+  }
+
+  /* We are done with this file */
+  fclose(infile);
+
+  /* QLA cooling assumes cooling->Redshifts table is in increasing order. Test
+   * this. */
+  for (int i = 0; i < N_Redshifts - 1; i++) {
+    if (cooling->Redshifts[i + 1] < cooling->Redshifts[i]) {
+      error("table should be in increasing order\n");
+    }
+  }
+}
+
+/**
+ * @brief Reads in QLA cooling table header. Consists of tables
+ * of values for temperature, hydrogen number density, helium fraction
+ * solar element abundances, and elements used to index the cooling tables.
+ *
+ * @param fname Filepath for cooling table from which to read header
+ * @param cooling Cooling data structure
+ */
+void read_cooling_header(const char *fname,
+                         struct cooling_function_data *cooling) {
+
+#ifdef HAVE_HDF5
+
+  int N_Temp, N_nH, N_He, N_SolarAbundances, N_Elements;
+
+  /* read sizes of array dimensions */
+  hid_t tempfile_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
+  if (tempfile_id < 0) error("unable to open file %s\n", fname);
+
+  /* read size of each table of values */
+  hid_t dataset =
+      H5Dopen(tempfile_id, "/Header/Number_of_temperature_bins", H5P_DEFAULT);
+  herr_t status =
+      H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &N_Temp);
+  if (status < 0) error("error reading number of temperature bins");
+  status = H5Dclose(dataset);
+  if (status < 0) error("error closing cooling dataset");
+
+  /* Check value */
+  if (N_Temp != qla_cooling_N_temperature)
+    error("Invalid temperature array length.");
+
+  dataset = H5Dopen(tempfile_id, "/Header/Number_of_density_bins", H5P_DEFAULT);
+  status =
+      H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &N_nH);
+  if (status < 0) error("error reading number of density bins");
+  status = H5Dclose(dataset);
+  if (status < 0) error("error closing cooling dataset");
+
+  /* Check value */
+  if (N_nH != qla_cooling_N_density) error("Invalid density array length.");
+
+  dataset =
+      H5Dopen(tempfile_id, "/Header/Number_of_helium_fractions", H5P_DEFAULT);
+  status =
+      H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &N_He);
+  if (status < 0) error("error reading number of He fraction bins");
+  status = H5Dclose(dataset);
+  if (status < 0) error("error closing cooling dataset");
+
+  /* Check value */
+  if (N_He != qla_cooling_N_He_frac)
+    error("Invalid Helium fraction array length.");
+
+  dataset = H5Dopen(tempfile_id, "/Header/Abundances/Number_of_abundances",
+                    H5P_DEFAULT);
+  status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                   &N_SolarAbundances);
+  if (status < 0) error("error reading number of solar abundance bins");
+  status = H5Dclose(dataset);
+  if (status < 0) error("error closing cooling dataset");
+
+  /* Check value */
+  if (N_SolarAbundances != qla_cooling_N_abundances)
+    error("Invalid solar abundances array length.");
+
+  dataset = H5Dopen(tempfile_id, "/Header/Number_of_metals", H5P_DEFAULT);
+  status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                   &N_Elements);
+  if (status < 0) error("error reading number of metal bins");
+  status = H5Dclose(dataset);
+  if (status < 0) error("error closing cooling dataset");
+
+  /* Check value */
+  if (N_Elements != qla_cooling_N_metal) error("Invalid metal array length.");
+
+  /* allocate arrays of values for each of the above quantities */
+  if (swift_memalign("cooling", (void **)&cooling->Temp, SWIFT_STRUCT_ALIGNMENT,
+                     N_Temp * sizeof(float)) != 0)
+    error("Failed to allocate temperature table");
+  if (swift_memalign("cooling", (void **)&cooling->Therm,
+                     SWIFT_STRUCT_ALIGNMENT, N_Temp * sizeof(float)) != 0)
+    error("Failed to allocate internal energy table");
+  if (swift_memalign("cooling", (void **)&cooling->nH, SWIFT_STRUCT_ALIGNMENT,
+                     N_nH * sizeof(float)) != 0)
+    error("Failed to allocate nH table");
+  if (swift_memalign("cooling", (void **)&cooling->HeFrac,
+                     SWIFT_STRUCT_ALIGNMENT, N_He * sizeof(float)) != 0)
+    error("Failed to allocate HeFrac table");
+  if (swift_memalign("cooling", (void **)&cooling->SolarAbundances,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     N_SolarAbundances * sizeof(float)) != 0)
+    error("Failed to allocate Solar abundances table");
+  if (swift_memalign("cooling", (void **)&cooling->SolarAbundances_inv,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     N_SolarAbundances * sizeof(float)) != 0)
+    error("Failed to allocate Solar abundances inverses table");
+
+  /* read in values for each of the arrays */
+  dataset = H5Dopen(tempfile_id, "/Solar/Temperature_bins", H5P_DEFAULT);
+  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                   cooling->Temp);
+  if (status < 0) error("error reading temperature bins");
+  status = H5Dclose(dataset);
+  if (status < 0) error("error closing cooling dataset");
+
+  dataset = H5Dopen(tempfile_id, "/Solar/Hydrogen_density_bins", H5P_DEFAULT);
+  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                   cooling->nH);
+  if (status < 0) error("error reading H density bins");
+  status = H5Dclose(dataset);
+  if (status < 0) error("error closing cooling dataset");
+
+  dataset = H5Dopen(tempfile_id, "/Metal_free/Helium_mass_fraction_bins",
+                    H5P_DEFAULT);
+  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                   cooling->HeFrac);
+  if (status < 0) error("error reading He fraction bins");
+  status = H5Dclose(dataset);
+  if (status < 0) error("error closing cooling dataset");
+
+  dataset = H5Dopen(tempfile_id, "/Header/Abundances/Solar_mass_fractions",
+                    H5P_DEFAULT);
+  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                   cooling->SolarAbundances);
+  if (status < 0) error("error reading solar mass fraction bins");
+  status = H5Dclose(dataset);
+  if (status < 0) error("error closing cooling dataset");
+
+  dataset = H5Dopen(tempfile_id, "/Metal_free/Temperature/Energy_density_bins",
+                    H5P_DEFAULT);
+  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                   cooling->Therm);
+  if (status < 0) error("error reading internal energy bins");
+  status = H5Dclose(dataset);
+  if (status < 0) error("error closing cooling dataset");
+
+  /* Convert to temperature, density and internal energy arrays to log10 */
+  for (int i = 0; i < N_Temp; i++) {
+    cooling->Temp[i] = log10(cooling->Temp[i]);
+    cooling->Therm[i] = log10(cooling->Therm[i]);
+  }
+  for (int i = 0; i < N_nH; i++) {
+    cooling->nH[i] = log10(cooling->nH[i]);
+  }
+
+    /* Compute inverse of solar mass fractions */
+#if defined(__ICC)
+#pragma novector
+#endif
+  for (int i = 0; i < N_SolarAbundances; ++i) {
+    cooling->SolarAbundances_inv[i] = 1.f / cooling->SolarAbundances[i];
+  }
+
+#else
+  error("Need HDF5 to read cooling tables");
+#endif
+}
+
+/**
+ * @brief Allocate space for cooling tables.
+ *
+ * @param cooling #cooling_function_data structure
+ */
+void allocate_cooling_tables(struct cooling_function_data *restrict cooling) {
+
+  /* Allocate arrays to store cooling tables. Arrays contain two tables of
+   * cooling rates with one table being for the redshift above current redshift
+   * and one below. */
+
+  if (swift_memalign("cooling-tables", (void **)&cooling->table.metal_heating,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     qla_cooling_N_loaded_redshifts *
+                         num_elements_metal_heating * sizeof(float)) != 0)
+    error("Failed to allocate metal_heating array");
+
+  if (swift_memalign("cooling-tables",
+                     (void **)&cooling->table.electron_abundance,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     qla_cooling_N_loaded_redshifts *
+                         num_elements_electron_abundance * sizeof(float)) != 0)
+    error("Failed to allocate electron_abundance array");
+
+  if (swift_memalign("cooling-tables", (void **)&cooling->table.temperature,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     qla_cooling_N_loaded_redshifts * num_elements_temperature *
+                         sizeof(float)) != 0)
+    error("Failed to allocate temperature array");
+
+  if (swift_memalign("cooling-tables",
+                     (void **)&cooling->table.H_plus_He_heating,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     qla_cooling_N_loaded_redshifts *
+                         num_elements_HpHe_heating * sizeof(float)) != 0)
+    error("Failed to allocate H_plus_He_heating array");
+
+  if (swift_memalign("cooling-tables",
+                     (void **)&cooling->table.H_plus_He_electron_abundance,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     qla_cooling_N_loaded_redshifts *
+                         num_elements_HpHe_electron_abundance *
+                         sizeof(float)) != 0)
+    error("Failed to allocate H_plus_He_electron_abundance array");
+}
+
+/**
+ * @brief Get the redshift invariant table of cooling rates (before reionization
+ * at redshift ~9) Reads in table of cooling rates and electron abundances due
+ * to metals (depending on temperature, hydrogen number density), cooling rates
+ * and electron abundances due to hydrogen and helium (depending on temperature,
+ * hydrogen number density and helium fraction), and temperatures (depending on
+ * internal energy, hydrogen number density and helium fraction; note: this is
+ * distinct from table of temperatures read in ReadCoolingHeader, as that table
+ * is used to index the cooling, electron abundance tables, whereas this one is
+ * used to obtain temperature of particle)
+ *
+ * @param cooling #cooling_function_data structure
+ * @param photodis Are we loading the photo-dissociation table?
+ */
+void get_redshift_invariant_table(
+    struct cooling_function_data *restrict cooling, const int photodis) {
+#ifdef HAVE_HDF5
+
+  /* Temporary tables */
+  float *net_cooling_rate = NULL;
+  float *electron_abundance = NULL;
+  float *temperature = NULL;
+  float *he_net_cooling_rate = NULL;
+  float *he_electron_abundance = NULL;
+
+  /* Allocate arrays for reading in cooling tables.  */
+  if (swift_memalign("cooling-temp", (void **)&net_cooling_rate,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     num_elements_cooling_rate * sizeof(float)) != 0)
+    error("Failed to allocate net_cooling_rate array");
+  if (swift_memalign("cooling-temp", (void **)&electron_abundance,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     num_elements_electron_abundance * sizeof(float)) != 0)
+    error("Failed to allocate electron_abundance array");
+  if (swift_memalign("cooling-temp", (void **)&temperature,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     num_elements_temperature * sizeof(float)) != 0)
+    error("Failed to allocate temperature array");
+  if (swift_memalign("cooling-temp", (void **)&he_net_cooling_rate,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     num_elements_HpHe_heating * sizeof(float)) != 0)
+    error("Failed to allocate he_net_cooling_rate array");
+  if (swift_memalign("cooling-temp", (void **)&he_electron_abundance,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     num_elements_HpHe_electron_abundance * sizeof(float)) != 0)
+    error("Failed to allocate he_electron_abundance array");
+
+  /* Decide which high redshift table to read. Indices set in cooling_update */
+  char filename[qla_table_path_name_length + 21];
+  if (photodis) {
+    sprintf(filename, "%sz_photodis.hdf5", cooling->cooling_table_path);
+    message("Reading cooling table 'z_photodis.hdf5'");
+  } else {
+    sprintf(filename, "%sz_8.989nocompton.hdf5", cooling->cooling_table_path);
+    message("Reading cooling table 'z_8.989nocompton.hdf5'");
+  }
+
+  hid_t file_id = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
+  if (file_id < 0) error("unable to open file %s\n", filename);
+
+  char set_name[64];
+
+  /* read in cooling rates due to metals */
+  for (int specs = 0; specs < qla_cooling_N_metal; specs++) {
+
+    /* Read in the cooling rate for this metal */
+    sprintf(set_name, "/%s/Net_Cooling", qla_tables_element_names[specs]);
+    hid_t dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
+    herr_t status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
+                            H5P_DEFAULT, net_cooling_rate);
+    if (status < 0) error("error reading metal cooling rate table");
+    status = H5Dclose(dataset);
+    if (status < 0) error("error closing cooling dataset");
+
+    /* Transpose from order tables are stored in (temperature, nH)
+     * to (metal species, nH, temperature) where fastest
+     * varying index is on right. Tables contain cooling rates but we
+     * want rate of change of internal energy, hence minus sign. */
+    for (int j = 0; j < qla_cooling_N_temperature; j++) {
+      for (int k = 0; k < qla_cooling_N_density; k++) {
+
+        /* Index in the HDF5 table */
+        const int hdf5_index = row_major_index_2d(
+            j, k, qla_cooling_N_temperature, qla_cooling_N_density);
+
+        /* Index in the internal table */
+        const int internal_index = row_major_index_3d(
+            specs, k, j, qla_cooling_N_metal, qla_cooling_N_density,
+            qla_cooling_N_temperature);
+
+        /* Change the sign and transpose */
+        cooling->table.metal_heating[internal_index] =
+            -net_cooling_rate[hdf5_index];
+      }
+    }
+  }
+
+  /* read in cooling rates due to H + He */
+  strcpy(set_name, "/Metal_free/Net_Cooling");
+  hid_t dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
+  herr_t status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
+                          H5P_DEFAULT, he_net_cooling_rate);
+  if (status < 0) error("error reading metal free cooling rate table");
+  status = H5Dclose(dataset);
+  if (status < 0) error("error closing cooling dataset");
+
+  /* read in Temperatures */
+  strcpy(set_name, "/Metal_free/Temperature/Temperature");
+  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
+  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                   temperature);
+  if (status < 0) error("error reading temperature table");
+  status = H5Dclose(dataset);
+  if (status < 0) error("error closing cooling dataset");
+
+  /* read in H + He electron abundances */
+  strcpy(set_name, "/Metal_free/Electron_density_over_n_h");
+  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
+  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                   he_electron_abundance);
+  if (status < 0) error("error reading electron density table");
+  status = H5Dclose(dataset);
+  if (status < 0) error("error closing cooling dataset");
+
+  /* Transpose from order tables are stored in (helium fraction, temperature,
+   * nH) to (nH, helium fraction, temperature) where fastest
+   * varying index is on right. Tables contain cooling rates but we
+   * want rate of change of internal energy, hence minus sign. */
+  for (int i = 0; i < qla_cooling_N_He_frac; i++) {
+    for (int j = 0; j < qla_cooling_N_temperature; j++) {
+      for (int k = 0; k < qla_cooling_N_density; k++) {
+
+        /* Index in the HDF5 table */
+        const int hdf5_index = row_major_index_3d(
+            i, j, k, qla_cooling_N_He_frac, qla_cooling_N_temperature,
+            qla_cooling_N_density);
+
+        /* Index in the internal table */
+        const int internal_index = row_major_index_3d(
+            k, i, j, qla_cooling_N_density, qla_cooling_N_He_frac,
+            qla_cooling_N_temperature);
+
+        /* Change the sign and transpose */
+        cooling->table.H_plus_He_heating[internal_index] =
+            -he_net_cooling_rate[hdf5_index];
+
+        /* Convert to log T and transpose */
+        cooling->table.temperature[internal_index] =
+            log10(temperature[hdf5_index]);
+
+        /* Just transpose */
+        cooling->table.H_plus_He_electron_abundance[internal_index] =
+            he_electron_abundance[hdf5_index];
+      }
+    }
+  }
+
+  /* read in electron densities due to metals */
+  strcpy(set_name, "/Solar/Electron_density_over_n_h");
+  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
+  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                   electron_abundance);
+  if (status < 0) error("error reading solar electron density table");
+  status = H5Dclose(dataset);
+  if (status < 0) error("error closing cooling dataset");
+
+  /* Transpose from order tables are stored in (temperature, nH) to
+   * (nH, temperature) where fastest varying index is on right. */
+  for (int i = 0; i < qla_cooling_N_temperature; i++) {
+    for (int j = 0; j < qla_cooling_N_density; j++) {
+
+      /* Index in the HDF5 table */
+      const int hdf5_index = row_major_index_2d(i, j, qla_cooling_N_temperature,
+                                                qla_cooling_N_density);
+
+      /* Index in the internal table */
+      const int internal_index = row_major_index_2d(j, i, qla_cooling_N_density,
+                                                    qla_cooling_N_temperature);
+
+      /* Just transpose */
+      cooling->table.electron_abundance[internal_index] =
+          electron_abundance[hdf5_index];
+    }
+  }
+
+  status = H5Fclose(file_id);
+  if (status < 0) error("error closing file");
+
+  swift_free("cooling-temp", net_cooling_rate);
+  swift_free("cooling-temp", electron_abundance);
+  swift_free("cooling-temp", temperature);
+  swift_free("cooling-temp", he_net_cooling_rate);
+  swift_free("cooling-temp", he_electron_abundance);
+
+#ifdef SWIFT_DEBUG_CHECKS
+  message("done reading in redshift invariant table");
+#endif
+
+#else
+  error("Need HDF5 to read cooling tables");
+#endif
+}
+
+/**
+ * @brief Get redshift dependent table of cooling rates.
+ * Reads in table of cooling rates and electron abundances due to
+ * metals (depending on temperature, hydrogen number density), cooling rates and
+ * electron abundances due to hydrogen and helium (depending on temperature,
+ * hydrogen number density and helium fraction), and temperatures (depending on
+ * internal energy, hydrogen number density and helium fraction; note: this is
+ * distinct from table of temperatures read in ReadCoolingHeader, as that table
+ * is used to index the cooling, electron abundance tables, whereas this one is
+ * used to obtain temperature of particle)
+ *
+ * @param cooling #cooling_function_data structure
+ * @param low_z_index Index of the lowest redshift table to load.
+ * @param high_z_index Index of the highest redshift table to load.
+ */
+void get_cooling_table(struct cooling_function_data *restrict cooling,
+                       const int low_z_index, const int high_z_index) {
+
+#ifdef HAVE_HDF5
+
+  /* Temporary tables */
+  float *net_cooling_rate = NULL;
+  float *electron_abundance = NULL;
+  float *temperature = NULL;
+  float *he_net_cooling_rate = NULL;
+  float *he_electron_abundance = NULL;
+
+  /* Allocate arrays for reading in cooling tables.  */
+  if (swift_memalign("cooling-temp", (void **)&net_cooling_rate,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     num_elements_cooling_rate * sizeof(float)) != 0)
+    error("Failed to allocate net_cooling_rate array");
+  if (swift_memalign("cooling-temp", (void **)&electron_abundance,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     num_elements_electron_abundance * sizeof(float)) != 0)
+    error("Failed to allocate electron_abundance array");
+  if (swift_memalign("cooling-temp", (void **)&temperature,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     num_elements_temperature * sizeof(float)) != 0)
+    error("Failed to allocate temperature array");
+  if (swift_memalign("cooling-temp", (void **)&he_net_cooling_rate,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     num_elements_HpHe_heating * sizeof(float)) != 0)
+    error("Failed to allocate he_net_cooling_rate array");
+  if (swift_memalign("cooling-temp", (void **)&he_electron_abundance,
+                     SWIFT_STRUCT_ALIGNMENT,
+                     num_elements_HpHe_electron_abundance * sizeof(float)) != 0)
+    error("Failed to allocate he_electron_abundance array");
+
+  /* Read in tables, transpose so that values for indices which vary most are
+   * adjacent. Repeat for redshift above and redshift below current value.  */
+  for (int z_index = low_z_index; z_index <= high_z_index; z_index++) {
+
+    /* Index along redhsift dimension for the subset of tables we read */
+    const int local_z_index = z_index - low_z_index;
+
+#ifdef SWIFT_DEBUG_CHECKS
+    if (local_z_index >= eagle_cooling_N_loaded_redshifts)
+      error("Reading invalid number of tables along z axis.");
+#endif
+
+    /* Open table for this redshift index */
+    char fname[qla_table_path_name_length + 12];
+    sprintf(fname, "%sz_%1.3f.hdf5", cooling->cooling_table_path,
+            cooling->Redshifts[z_index]);
+    message("Reading cooling table 'z_%1.3f.hdf5'",
+            cooling->Redshifts[z_index]);
+
+    hid_t file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
+    if (file_id < 0) error("unable to open file %s", fname);
+
+    char set_name[64];
+
+    /* read in cooling rates due to metals */
+    for (int specs = 0; specs < qla_cooling_N_metal; specs++) {
+
+      sprintf(set_name, "/%s/Net_Cooling", qla_tables_element_names[specs]);
+      hid_t dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
+      herr_t status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
+                              H5P_DEFAULT, net_cooling_rate);
+      if (status < 0) error("error reading metal cooling rate table");
+      status = H5Dclose(dataset);
+      if (status < 0) error("error closing cooling dataset");
+
+      /* Transpose from order tables are stored in (temperature, nH)
+       * to (metal species, redshift, nH, temperature) where fastest
+       * varying index is on right. Tables contain cooling rates but we
+       * want rate of change of internal energy, hence minus sign. */
+      for (int i = 0; i < qla_cooling_N_density; i++) {
+        for (int j = 0; j < qla_cooling_N_temperature; j++) {
+
+          /* Index in the HDF5 table */
+          const int hdf5_index = row_major_index_2d(
+              j, i, qla_cooling_N_temperature, qla_cooling_N_density);
+
+          /* Index in the internal table */
+          const int internal_index = row_major_index_4d(
+              specs, local_z_index, i, j, qla_cooling_N_metal,
+              qla_cooling_N_loaded_redshifts, qla_cooling_N_density,
+              qla_cooling_N_temperature);
+
+          /* Change the sign and transpose */
+          cooling->table.metal_heating[internal_index] =
+              -net_cooling_rate[hdf5_index];
+        }
+      }
+    }
+
+    /* read in cooling rates due to H + He */
+    strcpy(set_name, "/Metal_free/Net_Cooling");
+    hid_t dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
+    herr_t status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
+                            H5P_DEFAULT, he_net_cooling_rate);
+    if (status < 0) error("error reading metal free cooling rate table");
+    status = H5Dclose(dataset);
+    if (status < 0) error("error closing cooling dataset");
+
+    /* read in Temperature */
+    strcpy(set_name, "/Metal_free/Temperature/Temperature");
+    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
+    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                     temperature);
+    if (status < 0) error("error reading temperature table");
+    status = H5Dclose(dataset);
+    if (status < 0) error("error closing cooling dataset");
+
+    /* Read in H + He electron abundance */
+    strcpy(set_name, "/Metal_free/Electron_density_over_n_h");
+    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
+    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                     he_electron_abundance);
+    if (status < 0) error("error reading electron density table");
+    status = H5Dclose(dataset);
+    if (status < 0) error("error closing cooling dataset");
+
+    /* Transpose from order tables are stored in (helium fraction, temperature,
+     * nH) to (redshift, nH, helium fraction, temperature) where fastest
+     * varying index is on right. */
+    for (int i = 0; i < qla_cooling_N_He_frac; i++) {
+      for (int j = 0; j < qla_cooling_N_temperature; j++) {
+        for (int k = 0; k < qla_cooling_N_density; k++) {
+
+          /* Index in the HDF5 table */
+          const int hdf5_index = row_major_index_3d(
+              i, j, k, qla_cooling_N_He_frac, qla_cooling_N_temperature,
+              qla_cooling_N_density);
+
+          /* Index in the internal table */
+          const int internal_index = row_major_index_4d(
+              local_z_index, k, i, j, qla_cooling_N_loaded_redshifts,
+              qla_cooling_N_density, qla_cooling_N_He_frac,
+              qla_cooling_N_temperature);
+
+          /* Change the sign and transpose */
+          cooling->table.H_plus_He_heating[internal_index] =
+              -he_net_cooling_rate[hdf5_index];
+
+          /* Convert to log T and transpose */
+          cooling->table.temperature[internal_index] =
+              log10(temperature[hdf5_index]);
+
+          /* Just transpose */
+          cooling->table.H_plus_He_electron_abundance[internal_index] =
+              he_electron_abundance[hdf5_index];
+        }
+      }
+    }
+
+    /* read in electron densities due to metals */
+    strcpy(set_name, "/Solar/Electron_density_over_n_h");
+    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
+    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                     electron_abundance);
+    if (status < 0) error("error reading solar electron density table");
+    status = H5Dclose(dataset);
+    if (status < 0) error("error closing cooling dataset");
+
+    /* Transpose from order tables are stored in (temperature, nH) to
+     * (redshift, nH, temperature) where fastest varying index is on right. */
+    for (int i = 0; i < qla_cooling_N_temperature; i++) {
+      for (int j = 0; j < qla_cooling_N_density; j++) {
+
+        /* Index in the HDF5 table */
+        const int hdf5_index = row_major_index_2d(
+            i, j, qla_cooling_N_temperature, qla_cooling_N_density);
+
+        /* Index in the internal table */
+        const int internal_index = row_major_index_3d(
+            local_z_index, j, i, qla_cooling_N_loaded_redshifts,
+            qla_cooling_N_density, qla_cooling_N_temperature);
+
+        /* Just transpose */
+        cooling->table.electron_abundance[internal_index] =
+            electron_abundance[hdf5_index];
+      }
+    }
+
+    status = H5Fclose(file_id);
+    if (status < 0) error("error closing file");
+  }
+
+  swift_free("cooling-temp", net_cooling_rate);
+  swift_free("cooling-temp", electron_abundance);
+  swift_free("cooling-temp", temperature);
+  swift_free("cooling-temp", he_net_cooling_rate);
+  swift_free("cooling-temp", he_electron_abundance);
+
+#ifdef SWIFT_DEBUG_CHECKS
+  message("Done reading in general cooling table");
+#endif
+
+#else
+  error("Need HDF5 to read cooling tables");
+#endif
+}
diff --git a/src/cooling/QLA/cooling_tables.h b/src/cooling/QLA/cooling_tables.h
new file mode 100644
index 0000000000000000000000000000000000000000..67aa437fd9823f07a8543a0ef3cfc7ef2bab1b4e
--- /dev/null
+++ b/src/cooling/QLA/cooling_tables.h
@@ -0,0 +1,65 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_QLA_COOL_TABLES_H
+#define SWIFT_QLA_COOL_TABLES_H
+
+/**
+ * @file src/cooling/QLA/cooling.h
+ * @brief Quick Lyman-alpha cooling function (EAGLE with primordial Z only)
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+#include "cooling_struct.h"
+
+/*! Number of different bins along the redhsift axis of the tables */
+#define qla_cooling_N_redshifts 49
+
+/*! Number of redshift bins loaded at any given point int time */
+#define qla_cooling_N_loaded_redshifts 2
+
+/*! Number of different bins along the temperature axis of the tables */
+#define qla_cooling_N_temperature 176
+
+/*! Number of different bins along the density axis of the tables */
+#define qla_cooling_N_density 41
+
+/*! Number of different bins along the metal axis of the tables */
+#define qla_cooling_N_metal 9
+
+/*! Number of different bins along the metal axis of the tables */
+#define qla_cooling_N_He_frac 7
+
+/*! Number of different bins along the abundances axis of the tables */
+#define qla_cooling_N_abundances 11
+
+void get_cooling_redshifts(struct cooling_function_data *cooling);
+
+void read_cooling_header(const char *fname,
+                         struct cooling_function_data *cooling);
+
+void allocate_cooling_tables(struct cooling_function_data *restrict cooling);
+
+void get_redshift_invariant_table(
+    struct cooling_function_data *restrict cooling, const int photodis);
+void get_cooling_table(struct cooling_function_data *restrict cooling,
+                       const int low_z_index, const int high_z_index);
+
+#endif /* SWIFT_QLA_COOL_TABLES_H */
diff --git a/src/cooling/QLA/interpolate.h b/src/cooling/QLA/interpolate.h
new file mode 100644
index 0000000000000000000000000000000000000000..d49a482b1e0d031128052a0a79a905e5dcc56a48
--- /dev/null
+++ b/src/cooling/QLA/interpolate.h
@@ -0,0 +1,499 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_INTERPOL_QLA_H
+#define SWIFT_INTERPOL_QLA_H
+
+/**
+ * @file src/cooling/QLA/interpolate.h
+ * @brief Interpolation functions for Quick Lyman-alpha cooling tables
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local includes. */
+#include "align.h"
+#include "error.h"
+#include "inline.h"
+
+/**
+ * @brief Returns the 1d index of element with 2d indices x,y
+ * from a flattened 2d array in row major order
+ *
+ * @param x, y Indices of element of interest
+ * @param Nx, Ny Sizes of array dimensions
+ */
+__attribute__((always_inline)) INLINE int row_major_index_2d(const int x,
+                                                             const int y,
+                                                             const int Nx,
+                                                             const int Ny) {
+#ifdef SWIFT_DEBUG_CHECKS
+  assert(x < Nx);
+  assert(y < Ny);
+#endif
+  return x * Ny + y;
+}
+
+/**
+ * @brief Returns the 1d index of element with 3d indices x,y,z
+ * from a flattened 3d array in row major order
+ *
+ * @param x, y, z Indices of element of interest
+ * @param Nx, Ny, Nz Sizes of array dimensions
+ */
+__attribute__((always_inline)) INLINE int row_major_index_3d(
+    const int x, const int y, const int z, const int Nx, const int Ny,
+    const int Nz) {
+#ifdef SWIFT_DEBUG_CHECKS
+  assert(x < Nx);
+  assert(y < Ny);
+  assert(z < Nz);
+#endif
+  return x * Ny * Nz + y * Nz + z;
+}
+
+/**
+ * @brief Returns the 1d index of element with 4d indices x,y,z,w
+ * from a flattened 4d array in row major order
+ *
+ * @param x, y, z, w Indices of element of interest
+ * @param Nx, Ny, Nz, Nw Sizes of array dimensions
+ */
+__attribute__((always_inline)) INLINE int row_major_index_4d(
+    const int x, const int y, const int z, const int w, const int Nx,
+    const int Ny, const int Nz, const int Nw) {
+#ifdef SWIFT_DEBUG_CHECKS
+  assert(x < Nx);
+  assert(y < Ny);
+  assert(z < Nz);
+  assert(w < Nw);
+#endif
+  return x * Ny * Nz * Nw + y * Nz * Nw + z * Nw + w;
+}
+
+/**
+ * @brief Finds the index of a value in a table and compute delta to nearest
+ * element.
+ *
+ * This function assumes the table is monotonically increasing with a constant
+ * difference between adjacent values.
+ *
+ * The returned difference is expressed in units of the table separation. This
+ * means dx = (x - table[i]) / (table[i+1] - table[i]). It is always between
+ * 0 and 1.
+ *
+ * We use a small epsilon of 1e-4 to avoid out-of-range accesses due to
+ * rounding errors.
+ *
+ * @param table The table to search in.
+ * @param size The number of elements in the table.
+ * @param x The value to search for.
+ * @param i (return) The index in the table of the element.
+ * @param *dx (return) The difference between x and table[i]
+ */
+__attribute__((always_inline)) INLINE void get_index_1d(
+    const float *restrict table, const int size, const float x, int *i,
+    float *restrict dx) {
+
+  /* Small epsilon to avoid rounding issues leading to out-of-bound
+   * access when using the indices later to read data from the tables. */
+  const float epsilon = 1e-4f;
+
+  /* Indicate that the whole array is aligned on boundaries */
+  swift_align_information(float, table, SWIFT_STRUCT_ALIGNMENT);
+
+  /* Distance between elements in the array */
+  const float delta = (size - 1) / (table[size - 1] - table[0]);
+
+  if (x < table[0] + epsilon) {
+    /* We are below the first element */
+    *i = 0;
+    *dx = 0.f;
+  } else if (x < table[size - 1] - epsilon) {
+    /* Normal case */
+    *i = (x - table[0]) * delta;
+
+#ifdef SWIFT_DEBUG_CHECKS
+    if (*i > size || *i < 0) {
+      error(
+          "trying to get index for value outside table range. Table size: %d, "
+          "calculated index: %d, value: %.5e, table[0]: %.5e, grid size: %.5e",
+          size, *i, x, table[0], delta);
+    }
+#endif
+
+    *dx = (x - table[*i]) * delta;
+  } else {
+    /* We are after the last element */
+    *i = size - 2;
+    *dx = 1.f;
+  }
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (*dx < -0.001f || *dx > 1.001f) error("Invalid distance found dx=%e", *dx);
+#endif
+}
+
+/**
+ * @brief Interpolate a flattened 2D table at a given position.
+ *
+ * This function uses linear interpolation along each axis. It also
+ * assumes that the table is aligned on SWIFT_STRUCT_ALIGNMENT.
+ *
+ * @param table The 2D table to interpolate.
+ * @param xi, yi Indices of element of interest.
+ * @param Nx, Ny Sizes of array dimensions.
+ * @param dx, dy Distance between the point and the index in units of
+ * the grid spacing.
+ */
+__attribute__((always_inline)) INLINE float interpolation_2d(
+    const float *table, const int xi, const int yi, const float dx,
+    const float dy, const int Nx, const int Ny) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (dx < -0.001f || dx > 1.001f) error("Invalid dx=%e", dx);
+  if (dy < -0.001f || dy > 1.001f) error("Invalid dy=%e", dy);
+#endif
+
+  const float tx = 1.f - dx;
+  const float ty = 1.f - dy;
+
+  /* Indicate that the whole array is aligned on boundaries */
+  swift_align_information(float, table, SWIFT_STRUCT_ALIGNMENT);
+
+  /* Linear interpolation along each axis. We read the table 2^2=4 times */
+  float result = tx * ty * table[row_major_index_2d(xi + 0, yi + 0, Nx, Ny)];
+
+  result += tx * dy * table[row_major_index_2d(xi + 0, yi + 1, Nx, Ny)];
+  result += dx * ty * table[row_major_index_2d(xi + 1, yi + 0, Nx, Ny)];
+
+  result += dx * dy * table[row_major_index_2d(xi + 1, yi + 1, Nx, Ny)];
+
+  return result;
+}
+
+/**
+ * @brief Interpolate a flattened 3D table at a given position.
+ *
+ * This function uses linear interpolation along each axis. It also
+ * assumes that the table is aligned on SWIFT_STRUCT_ALIGNMENT.
+ *
+ * @param table The 3D table to interpolate.
+ * @param xi, yi, zi Indices of element of interest.
+ * @param Nx, Ny, Nz Sizes of array dimensions.
+ * @param dx, dy, dz Distance between the point and the index in units of
+ * the grid spacing.
+ */
+__attribute__((always_inline)) INLINE float interpolation_3d(
+    const float *table, const int xi, const int yi, const int zi,
+    const float dx, const float dy, const float dz, const int Nx, const int Ny,
+    const int Nz) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (dx < -0.001f || dx > 1.001f) error("Invalid dx=%e", dx);
+  if (dy < -0.001f || dy > 1.001f) error("Invalid dy=%e", dy);
+  if (dz < -0.001f || dz > 1.001f) error("Invalid dz=%e", dz);
+#endif
+
+  const float tx = 1.f - dx;
+  const float ty = 1.f - dy;
+  const float tz = 1.f - dz;
+
+  /* Indicate that the whole array is aligned on page boundaries */
+  swift_align_information(float, table, SWIFT_STRUCT_ALIGNMENT);
+
+  /* Linear interpolation along each axis. We read the table 2^3=8 times */
+  float result = tx * ty * tz *
+                 table[row_major_index_3d(xi + 0, yi + 0, zi + 0, Nx, Ny, Nz)];
+
+  result += tx * ty * dz *
+            table[row_major_index_3d(xi + 0, yi + 0, zi + 1, Nx, Ny, Nz)];
+  result += tx * dy * tz *
+            table[row_major_index_3d(xi + 0, yi + 1, zi + 0, Nx, Ny, Nz)];
+  result += dx * ty * tz *
+            table[row_major_index_3d(xi + 1, yi + 0, zi + 0, Nx, Ny, Nz)];
+
+  result += tx * dy * dz *
+            table[row_major_index_3d(xi + 0, yi + 1, zi + 1, Nx, Ny, Nz)];
+  result += dx * ty * dz *
+            table[row_major_index_3d(xi + 1, yi + 0, zi + 1, Nx, Ny, Nz)];
+  result += dx * dy * tz *
+            table[row_major_index_3d(xi + 1, yi + 1, zi + 0, Nx, Ny, Nz)];
+
+  result += dx * dy * dz *
+            table[row_major_index_3d(xi + 1, yi + 1, zi + 1, Nx, Ny, Nz)];
+
+  return result;
+}
+
+/**
+ * @brief Interpolate a flattened 3D table at a given position but avoid the
+ * x-dimension.
+ *
+ * This function uses linear interpolation along each axis.
+ * We look at the xi coordoniate but do not interpolate around it. We just
+ * interpolate the remaining 2 dimensions.
+ * The function also assumes that the table is aligned on
+ * SWIFT_STRUCT_ALIGNMENT.
+ *
+ * @param table The 3D table to interpolate.
+ * @param xi, yi, zi Indices of element of interest.
+ * @param Nx, Ny, Nz Sizes of array dimensions.
+ * @param dx, dy, dz Distance between the point and the index in units of
+ * the grid spacing.
+ */
+__attribute__((always_inline)) INLINE float interpolation_3d_no_x(
+    const float *table, const int xi, const int yi, const int zi,
+    const float dx, const float dy, const float dz, const int Nx, const int Ny,
+    const int Nz) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (dx != 0.f) error("Attempting to interpolate along x!");
+  if (dy < -0.001f || dy > 1.001f) error("Invalid dy=%e", dy);
+  if (dz < -0.001f || dz > 1.001f) error("Invalid dz=%e", dz);
+#endif
+
+  const float tx = 1.f;
+  const float ty = 1.f - dy;
+  const float tz = 1.f - dz;
+
+  /* Indicate that the whole array is aligned on page boundaries */
+  swift_align_information(float, table, SWIFT_STRUCT_ALIGNMENT);
+
+  /* Linear interpolation along each axis. We read the table 2^2=4 times */
+  /* Note that we intentionally kept the table access along the axis where */
+  /* we do not interpolate as comments in the code to allow readers to */
+  /* understand what is going on. */
+  float result = tx * ty * tz *
+                 table[row_major_index_3d(xi + 0, yi + 0, zi + 0, Nx, Ny, Nz)];
+
+  result += tx * ty * dz *
+            table[row_major_index_3d(xi + 0, yi + 0, zi + 1, Nx, Ny, Nz)];
+  result += tx * dy * tz *
+            table[row_major_index_3d(xi + 0, yi + 1, zi + 0, Nx, Ny, Nz)];
+  /* result += dx * ty * tz * */
+  /*           table[row_major_index_3d(xi + 1, yi + 0, zi + 0, Nx, Ny, Nz)]; */
+
+  result += tx * dy * dz *
+            table[row_major_index_3d(xi + 0, yi + 1, zi + 1, Nx, Ny, Nz)];
+  /* result += dx * ty * dz * */
+  /*           table[row_major_index_3d(xi + 1, yi + 0, zi + 1, Nx, Ny, Nz)]; */
+  /* result += dx * dy * tz * */
+  /*           table[row_major_index_3d(xi + 1, yi + 1, zi + 0, Nx, Ny, Nz)]; */
+
+  /* result += dx * dy * dz * */
+  /*           table[row_major_index_3d(xi + 1, yi + 1, zi + 1, Nx, Ny, Nz)]; */
+
+  return result;
+}
+
+/**
+ * @brief Interpolate a flattened 4D table at a given position.
+ *
+ * This function uses linear interpolation along each axis. It also
+ * assumes that the table is aligned on SWIFT_STRUCT_ALIGNMENT.
+ *
+ * @param table The 4D table to interpolate.
+ * @param xi, yi, zi, wi Indices of element of interest.
+ * @param Nx, Ny, Nz, Nw Sizes of array dimensions.
+ * @param dx, dy, dz, dw Distance between the point and the index in units of
+ * the grid spacing.
+ */
+__attribute__((always_inline)) INLINE float interpolation_4d(
+    const float *table, const int xi, const int yi, const int zi, const int wi,
+    const float dx, const float dy, const float dz, const float dw,
+    const int Nx, const int Ny, const int Nz, const int Nw) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (dx < -0.001f || dx > 1.001f) error("Invalid dx=%e", dx);
+  if (dy < -0.001f || dy > 1.001f) error("Invalid dy=%e", dy);
+  if (dz < -0.001f || dz > 1.001f) error("Invalid dz=%e", dz);
+  if (dw < -0.001f || dw > 1.001f) error("Invalid dw=%e", dw);
+#endif
+
+  const float tx = 1.f - dx;
+  const float ty = 1.f - dy;
+  const float tz = 1.f - dz;
+  const float tw = 1.f - dw;
+
+  /* Indicate that the whole array is aligned on page boundaries */
+  swift_align_information(float, table, SWIFT_STRUCT_ALIGNMENT);
+
+  /* Linear interpolation along each axis. We read the table 2^4=16 times */
+  float result =
+      tx * ty * tz * tw *
+      table[row_major_index_4d(xi + 0, yi + 0, zi + 0, wi + 0, Nx, Ny, Nz, Nw)];
+
+  result +=
+      tx * ty * tz * dw *
+      table[row_major_index_4d(xi + 0, yi + 0, zi + 0, wi + 1, Nx, Ny, Nz, Nw)];
+  result +=
+      tx * ty * dz * tw *
+      table[row_major_index_4d(xi + 0, yi + 0, zi + 1, wi + 0, Nx, Ny, Nz, Nw)];
+  result +=
+      tx * dy * tz * tw *
+      table[row_major_index_4d(xi + 0, yi + 1, zi + 0, wi + 0, Nx, Ny, Nz, Nw)];
+  result +=
+      dx * ty * tz * tw *
+      table[row_major_index_4d(xi + 1, yi + 0, zi + 0, wi + 0, Nx, Ny, Nz, Nw)];
+
+  result +=
+      tx * ty * dz * dw *
+      table[row_major_index_4d(xi + 0, yi + 0, zi + 1, wi + 1, Nx, Ny, Nz, Nw)];
+  result +=
+      tx * dy * tz * dw *
+      table[row_major_index_4d(xi + 0, yi + 1, zi + 0, wi + 1, Nx, Ny, Nz, Nw)];
+  result +=
+      dx * ty * tz * dw *
+      table[row_major_index_4d(xi + 1, yi + 0, zi + 0, wi + 1, Nx, Ny, Nz, Nw)];
+  result +=
+      tx * dy * dz * tw *
+      table[row_major_index_4d(xi + 0, yi + 1, zi + 1, wi + 0, Nx, Ny, Nz, Nw)];
+  result +=
+      dx * ty * dz * tw *
+      table[row_major_index_4d(xi + 1, yi + 0, zi + 1, wi + 0, Nx, Ny, Nz, Nw)];
+  result +=
+      dx * dy * tz * tw *
+      table[row_major_index_4d(xi + 1, yi + 1, zi + 0, wi + 0, Nx, Ny, Nz, Nw)];
+
+  result +=
+      dx * dy * dz * tw *
+      table[row_major_index_4d(xi + 1, yi + 1, zi + 1, wi + 0, Nx, Ny, Nz, Nw)];
+  result +=
+      dx * dy * tz * dw *
+      table[row_major_index_4d(xi + 1, yi + 1, zi + 0, wi + 1, Nx, Ny, Nz, Nw)];
+  result +=
+      dx * ty * dz * dw *
+      table[row_major_index_4d(xi + 1, yi + 0, zi + 1, wi + 1, Nx, Ny, Nz, Nw)];
+  result +=
+      tx * dy * dz * dw *
+      table[row_major_index_4d(xi + 0, yi + 1, zi + 1, wi + 1, Nx, Ny, Nz, Nw)];
+
+  result +=
+      dx * dy * dz * dw *
+      table[row_major_index_4d(xi + 1, yi + 1, zi + 1, wi + 1, Nx, Ny, Nz, Nw)];
+
+  return result;
+}
+
+/**
+ * @brief Interpolate a flattened 4D table at a given position but avoid the
+ * x-dimension.
+ *
+ * This function uses linear interpolation along each axis.
+ * We look at the xi coordoniate but do not interpolate around it. We just
+ * interpolate the remaining 3 dimensions.
+ * The function also assumes that the table is aligned on
+ * SWIFT_STRUCT_ALIGNMENT.
+ *
+ * @param table The 4D table to interpolate.
+ * @param xi, yi, zi, wi Indices of element of interest.
+ * @param Nx, Ny, Nz, Nw Sizes of array dimensions.
+ * @param dx, dy, dz, dw Distance between the point and the index in units of
+ * the grid spacing.
+ */
+__attribute__((always_inline)) INLINE float interpolation_4d_no_x(
+    const float *table, const int xi, const int yi, const int zi, const int wi,
+    const float dx, const float dy, const float dz, const float dw,
+    const int Nx, const int Ny, const int Nz, const int Nw) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (dx != 0.f) error("Attempting to interpolate along x!");
+  if (dy < -0.001f || dy > 1.001f) error("Invalid dy=%e", dy);
+  if (dz < -0.001f || dz > 1.001f) error("Invalid dz=%e", dz);
+  if (dw < -0.001f || dw > 1.001f) error("Invalid dw=%e", dw);
+#endif
+
+  const float tx = 1.f;
+  const float ty = 1.f - dy;
+  const float tz = 1.f - dz;
+  const float tw = 1.f - dw;
+
+  /* Indicate that the whole array is aligned on boundaries */
+  swift_align_information(float, table, SWIFT_STRUCT_ALIGNMENT);
+
+  /* Linear interpolation along each axis. We read the table 2^3=8 times */
+  /* Note that we intentionally kept the table access along the axis where */
+  /* we do not interpolate as comments in the code to allow readers to */
+  /* understand what is going on. */
+  float result =
+      tx * ty * tz * tw *
+      table[row_major_index_4d(xi + 0, yi + 0, zi + 0, wi + 0, Nx, Ny, Nz, Nw)];
+
+  result +=
+      tx * ty * tz * dw *
+      table[row_major_index_4d(xi + 0, yi + 0, zi + 0, wi + 1, Nx, Ny, Nz, Nw)];
+  result +=
+      tx * ty * dz * tw *
+      table[row_major_index_4d(xi + 0, yi + 0, zi + 1, wi + 0, Nx, Ny, Nz, Nw)];
+  result +=
+      tx * dy * tz * tw *
+      table[row_major_index_4d(xi + 0, yi + 1, zi + 0, wi + 0, Nx, Ny, Nz, Nw)];
+  /* result += */
+  /*     dx * ty * tz * tw * */
+  /*     table[row_major_index_4d(xi + 1, yi + 0, zi + 0, wi + 0, Nx, Ny, Nz,
+   * Nw)]; */
+
+  result +=
+      tx * ty * dz * dw *
+      table[row_major_index_4d(xi + 0, yi + 0, zi + 1, wi + 1, Nx, Ny, Nz, Nw)];
+  result +=
+      tx * dy * tz * dw *
+      table[row_major_index_4d(xi + 0, yi + 1, zi + 0, wi + 1, Nx, Ny, Nz, Nw)];
+  /* result += */
+  /*     dx * ty * tz * dw * */
+  /*     table[row_major_index_4d(xi + 1, yi + 0, zi + 0, wi + 1, Nx, Ny, Nz,
+   * Nw)]; */
+  result +=
+      tx * dy * dz * tw *
+      table[row_major_index_4d(xi + 0, yi + 1, zi + 1, wi + 0, Nx, Ny, Nz, Nw)];
+  /* result += */
+  /*     dx * ty * dz * tw * */
+  /*     table[row_major_index_4d(xi + 1, yi + 0, zi + 1, wi + 0, Nx, Ny, Nz,
+   * Nw)]; */
+  /* result += */
+  /*     dx * dy * tz * tw * */
+  /*     table[row_major_index_4d(xi + 1, yi + 1, zi + 0, wi + 0, Nx, Ny, Nz, */
+  /* Nw)]; */
+
+  /* result += */
+  /*     dx * dy * dz * tw * */
+  /*     table[row_major_index_4d(xi + 1, yi + 1, zi + 1, wi + 0, Nx, Ny, Nz, */
+  /* Nw)]; */
+  /* result += */
+  /*     dx * dy * tz * dw * */
+  /*     table[row_major_index_4d(xi + 1, yi + 1, zi + 0, wi + 1, Nx, Ny, Nz, */
+  /* Nw)]; */
+  /* result += */
+  /*     dx * ty * dz * dw * */
+  /*     table[row_major_index_4d(xi + 1, yi + 0, zi + 1, wi + 1, Nx, Ny, Nz,
+   * Nw)]; */
+  result +=
+      tx * dy * dz * dw *
+      table[row_major_index_4d(xi + 0, yi + 1, zi + 1, wi + 1, Nx, Ny, Nz, Nw)];
+
+  /* result += */
+  /*     dx * dy * dz * dw * */
+  /*     table[row_major_index_4d(xi + 1, yi + 1, zi + 1, wi + 1, Nx, Ny, Nz, */
+  /* Nw)]; */
+
+  return result;
+}
+
+#endif
diff --git a/src/cooling_io.h b/src/cooling_io.h
index 1ced353d7ff8320a48731545300274c654a20744..e56aef8ba05dbc35b6f78c13d2039541df12c0df 100644
--- a/src/cooling_io.h
+++ b/src/cooling_io.h
@@ -33,6 +33,8 @@
 #include "./cooling/Compton/cooling_io.h"
 #elif defined(COOLING_GRACKLE)
 #include "./cooling/grackle/cooling_io.h"
+#elif defined(COOLING_QLA)
+#include "./cooling/QLA/cooling_io.h"
 #elif defined(COOLING_EAGLE)
 #include "./cooling/EAGLE/cooling_io.h"
 #else
diff --git a/src/cooling_struct.h b/src/cooling_struct.h
index 93de8d1b7a0bcfd56d2b1a503aea1e8339bc8016..48a57f1fa3640f1d6ca64355c463c5533349265e 100644
--- a/src/cooling_struct.h
+++ b/src/cooling_struct.h
@@ -38,6 +38,8 @@
 #include "./cooling/Compton/cooling_struct.h"
 #elif defined(COOLING_GRACKLE)
 #include "./cooling/grackle/cooling_struct.h"
+#elif defined(COOLING_QLA)
+#include "./cooling/QLA/cooling_struct.h"
 #elif defined(COOLING_EAGLE)
 #include "./cooling/EAGLE/cooling_struct.h"
 #else