diff --git a/AUTHORS b/AUTHORS
index 4d43745609fc9d0ac2f149256a76ed6c581c9144..c822300c22885a05b42d58a51cc86af9da410429 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -8,3 +8,5 @@ John A. Regan 		john.a.regan@durham.ac.uk
 Angus Lepper		angus.lepper@ed.ac.uk
 Tom Theuns 		tom.theuns@durham.ac.uk
 Richard G. Bower	r.g.bower@durham.ac.uk
+Stefan Arridge		stefan.arridge@durham.ac.uk
+Massimiliano Culpo	massimiliano.culpo@googlemail.com
diff --git a/README b/README
index 2fc16c28f5b14904e566cf734f2c0e34235f93ae..562a54f3104884b1bf4c5e607700279161fad4c9 100644
--- a/README
+++ b/README
@@ -19,6 +19,7 @@ Usage: swift [OPTION]... PARAMFILE
 Valid options are:
   -a          Pin runners using processor affinity
   -c          Run with cosmological time integration
+  -C          Run with cooling
   -d          Dry run. Read the parameter file, allocate memory but does not read 
               the particles from ICs and exit before the start of time integration.
               Allows user to check validy of parameter and IC files as well as memory limits.
@@ -30,8 +31,9 @@ Valid options are:
   -n    {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop. 
   -s          Run with SPH
   -t    {int} The number of threads to use on each MPI rank. Defaults to 1 if not specified.
-  -v     [12] Increase the level of verbosity 1: MPI-rank 0 writes 
-              2: All MPI-ranks write
+  -v     [12] Increase the level of verbosity
+  	      1: MPI-rank 0 writes
+	      2: All MPI-ranks write
   -y    {int} Time-step frequency at which task graphs are dumped
   -h          Print this help message and exit
 
diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index 4703c7091550c8c496952d9e96b623e180c78a69..0658dd27e043b11c7480eaf8b01c0c40809e215c 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -760,8 +760,10 @@ WARN_LOGFILE           =
 # Note: If this tag is empty the current directory is searched.
 
 INPUT                  =  @top_srcdir@ @top_srcdir@/src @top_srcdir@/tests @top_srcdir@/examples
-INPUT		       += @top_srcdir@/src/hydro/Minimal @top_srcdir@/src/gravity/Default
-INPUT		       += @top_srcdir@/src/riemann 
+INPUT		       += @top_srcdir@/src/hydro/Minimal
+INPUT		       += @top_srcdir@/src/gravity/Default
+INPUT		       += @top_srcdir@/src/riemann
+INPUT		       += @top_srcdir@/src/cooling/const_lambda
 
 # This tag can be used to specify the character encoding of the source files
 # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
diff --git a/examples/CoolingBox/coolingBox.yml b/examples/CoolingBox/coolingBox.yml
new file mode 100644
index 0000000000000000000000000000000000000000..e13de6095066836853d9e9068330938f6260f38e
--- /dev/null
+++ b/examples/CoolingBox/coolingBox.yml
@@ -0,0 +1,44 @@
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     2.0e33   # Solar masses
+  UnitLength_in_cgs:   3.01e21   # Kilparsecs
+  UnitVelocity_in_cgs: 1.0e5   # Time unit is cooling time
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1.0    # The end time of the simulation (in internal units).
+  dt_min:     1e-6  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            coolingBox # Common part of the name of output files
+  time_first:          0.         # Time of the first output (in internal units)
+  delta_time:          1.0e-1       # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-2 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      0.1      # The tolerance for the targetted number of neighbours.
+  max_smoothing_length:  0.1      # Maximal smoothing length allowed (in internal units).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./coolingBox.hdf5     # The file to read
+
+# Dimensionless pre-factor for the time-step condition
+LambdaCooling:
+  lambda:                      0.0    # Cooling rate (in cgs units)
+  minimum_temperature:         1.0e4  # Minimal temperature (Kelvin)
+  mean_molecular_weight:       0.59   # Mean molecular weight
+  hydrogen_mass_abundance:     0.75   # Hydrogen mass abundance (dimensionless)
+  cooling_tstep_mult:          1.0    # Dimensionless pre-factor for the time-step condition
+  
diff --git a/examples/CoolingBox/energy_plot.py b/examples/CoolingBox/energy_plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..28cf9ab64decb5b56e98118a407221fac2bd4f16
--- /dev/null
+++ b/examples/CoolingBox/energy_plot.py
@@ -0,0 +1,84 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import h5py as h5
+import sys
+
+stats_filename = "./energy.txt"
+snap_filename = "coolingBox_000.hdf5"
+#plot_dir = "./"
+
+#some constants in cgs units
+k_b = 1.38E-16 #boltzmann
+m_p = 1.67e-24 #proton mass
+#initial conditions set in makeIC.py
+rho = 3.2e3
+P = 4.5e6
+n_H_cgs = 0.0001
+gamma = 5./3.
+T_init = 1.0e5
+
+#Read the units parameters from the snapshot
+f = h5.File(snap_filename,'r')
+units = f["InternalCodeUnits"]
+unit_mass = units.attrs["Unit mass in cgs (U_M)"]
+unit_length = units.attrs["Unit length in cgs (U_L)"]
+unit_time = units.attrs["Unit time in cgs (U_t)"]
+parameters = f["Parameters"]
+cooling_lambda = float(parameters.attrs["LambdaCooling:lambda"])
+min_T = float(parameters.attrs["LambdaCooling:minimum_temperature"])
+mu = float(parameters.attrs["LambdaCooling:mean_molecular_weight"])
+X_H = float(parameters.attrs["LambdaCooling:hydrogen_mass_abundance"])
+
+#get number of particles
+header = f["Header"]
+n_particles = header.attrs["NumPart_ThisFile"][0]
+#read energy and time arrays
+array = np.genfromtxt(stats_filename,skip_header = 1)
+time = array[:,0]
+total_energy = array[:,2]
+total_mass = array[:,1]
+
+time = time[1:]
+total_energy = total_energy[1:]
+total_mass = total_mass[1:]
+
+#conversions to cgs
+rho_cgs = rho * unit_mass / (unit_length)**3
+time_cgs = time * unit_time
+u_init_cgs = total_energy[0]/(total_mass[0]) * unit_length**2 / (unit_time)**2 
+
+#find the energy floor
+print min_T
+u_floor_cgs = k_b * min_T / (mu * m_p * (gamma - 1.))
+#find analytic solution
+analytic_time = np.linspace(time_cgs[0],time_cgs[-1],1000)
+print time_cgs[1]
+print analytic_time[1]
+du_dt_cgs = -cooling_lambda * n_H_cgs**2 / rho_cgs
+u_analytic = du_dt_cgs*(analytic_time - analytic_time[0]) + u_init_cgs
+cooling_time = u_init_cgs/(-du_dt_cgs)
+#rescale energy to initial energy
+total_energy /= total_energy[0]
+u_analytic /= u_init_cgs
+u_floor_cgs /= u_init_cgs
+# plot_title = r"$\Lambda \, = \, %1.1g \mathrm{erg}\mathrm{cm^3}\mathrm{s^{-1}} \, \, T_{init} = %1.1g\mathrm{K} \, \, T_{floor} = %1.1g\mathrm{K} \, \, n_H = %1.1g\mathrm{cm^{-3}}$" %(cooling_lambda,T_init,T_floor,n_H)
+# plot_filename = "energy_plot_creasey_no_cooling_T_init_1p0e5_n_H_0p1.png"
+#analytic_solution = np.zeros(n_snaps-1)
+for i in range(u_analytic.size):
+    if u_analytic[i]<u_floor_cgs:
+        u_analytic[i] = u_floor_cgs
+plt.plot(time_cgs,total_energy,'k',label = "Numerical solution")
+plt.plot(analytic_time,u_analytic,'--r',lw = 2.0,label = "Analytic Solution")
+plt.plot((cooling_time,cooling_time),(0,1),'b',label = "Cooling time")
+plt.plot((time_cgs[0],time_cgs[0]),(0,1),'m',label = "First output")
+plt.title(r"$n_H = %1.1e \, \mathrm{cm}^{-3}$" %n_H_cgs)
+plt.xlabel("Time (seconds)")
+plt.ylabel("Energy/Initial energy")
+plt.ylim(0.999,1.001)
+plt.xlim(0,min(10.0*cooling_time,time_cgs[-1]))
+plt.legend(loc = "upper right")    
+if (int(sys.argv[1])==0):
+    plt.show()
+else:
+    plt.savefig(full_plot_filename,format = "png")
+    plt.close()
diff --git a/examples/CoolingBox/makeIC.py b/examples/CoolingBox/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..f35c9243d4fa71f872fd27520de14a23073c4b9d
--- /dev/null
+++ b/examples/CoolingBox/makeIC.py
@@ -0,0 +1,109 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # 
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ # 
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ # GNU General Public License for more details.
+ # 
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ # 
+ ##############################################################################
+
+import h5py
+import sys
+from numpy import *
+
+# Generates a swift IC file containing a cartesian distribution of particles
+# at a constant density and pressure in a cubic box
+
+# Parameters
+periodic= 1           # 1 For periodic box
+boxSize = 1           #1 kiloparsec    
+L = int(sys.argv[1])  # Number of particles along one axis
+rho = 3.2e3          # Density in code units (0.01 hydrogen atoms per cm^3)
+P = 4.5e6          # Pressure in code units (at 10^5K)
+gamma = 5./3.         # Gas adiabatic index
+eta = 1.2349          # 48 ngbs with cubic spline kernel
+fileName = "coolingBox.hdf5" 
+
+#---------------------------------------------------
+numPart = L**3
+mass = boxSize**3 * rho / numPart
+print mass
+internalEnergy = P / ((gamma - 1.)*rho)
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 3.08e21 
+grp.attrs["Unit mass in cgs (U_M)"] = 2.0e33 
+grp.attrs["Unit time in cgs (U_t)"] = 3.08e16 
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+#Particle group
+grp = file.create_group("/PartType0")
+
+v  = zeros((numPart, 3))
+ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+ds[()] = v
+v = zeros(1)
+
+m = full((numPart, 1), mass)
+ds = grp.create_dataset('Masses', (numPart,1), 'f')
+ds[()] = m
+m = zeros(1)
+
+h = full((numPart, 1), eta * boxSize / L)
+ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
+ds[()] = h
+h = zeros(1)
+
+u = full((numPart, 1), internalEnergy)
+ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
+ds[()] = u
+u = zeros(1)
+
+
+ids = linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1))
+ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
+ds[()] = ids + 1
+x      = ids % L;
+y      = ((ids - x) / L) % L;
+z      = (ids - x - L * y) / L**2;
+coords = zeros((numPart, 3))
+coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L)
+coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
+coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
+ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = coords
+
+file.close()
diff --git a/examples/CoolingBox/run.sh b/examples/CoolingBox/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..c78eec9da6c486bc31a60ab7a8521ce6a6a63165
--- /dev/null
+++ b/examples/CoolingBox/run.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+echo "Generating initial conditions for the cooling box example..."
+
+python makeIC.py 10
+
+../swift -s -t 1 coolingBox.yml -C 2>&1 | tee output.log
+
+python energy_plot.py 0
diff --git a/examples/CosmoVolume/cosmoVolume.yml b/examples/CosmoVolume/cosmoVolume.yml
index 548b85078952954f2bf97280be6feb25eb6ef444..13cea318144d296183d630a53d78c69d050c1abe 100644
--- a/examples/CosmoVolume/cosmoVolume.yml
+++ b/examples/CosmoVolume/cosmoVolume.yml
@@ -6,11 +6,6 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1   # Amperes
   UnitTemp_in_cgs:     1   # Kelvin
 
-# Parameters for the task scheduling
-Scheduler:
-  cell_sub_size:    6000     # Value used for the original scaling tests
-  cell_split_size:  300      # Value used for the original scaling tests
-
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
diff --git a/examples/EAGLE_12/eagle_12.yml b/examples/EAGLE_12/eagle_12.yml
index d9d4e73626deed0d13527a0ba3c4ece0e7a33624..80714d87f4afa7d7e4d41ce7bf56faed856208ef 100644
--- a/examples/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_12/eagle_12.yml
@@ -6,11 +6,6 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
 
-# Parameters for the task scheduling
-#Scheduler:
-  cell_sub_size:    64000000     # Value used for the original scaling tests
-#  cell_split_size:  300      # Value used for the original scaling tests
-
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
diff --git a/examples/EAGLE_25/eagle_25.yml b/examples/EAGLE_25/eagle_25.yml
index cd5a1211b68c4ec2bebcd565713cc52a7b8de6df..6afb737677040ba0605d4e3116800f079a059be4 100644
--- a/examples/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_25/eagle_25.yml
@@ -6,11 +6,6 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
 
-# Parameters for the task scheduling
-Scheduler:
-  cell_sub_size:    6000     # Value used for the original scaling tests
-  cell_split_size:  300      # Value used for the original scaling tests
-
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
diff --git a/examples/EAGLE_50/eagle_50.yml b/examples/EAGLE_50/eagle_50.yml
index 38a9c165796359dfdabef423154299ee04ae8c76..d9e5d46326780fe5b2abde025b50a7ec667b19b1 100644
--- a/examples/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_50/eagle_50.yml
@@ -6,11 +6,6 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
 
-# Parameters for the task scheduling
-Scheduler:
-  cell_sub_size:    6000     # Value used for the original scaling tests
-  cell_split_size:  300      # Value used for the original scaling tests
-
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
diff --git a/examples/Makefile.am b/examples/Makefile.am
index 38ddcf61935b2099a785fa9a7f427807177b7cc0..d9e1f2fe741098fe2051155fc1ff2d66d4751cee 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -68,6 +68,7 @@ swift_fixdt_mpi_LDADD =  ../src/.libs/libswiftsim_mpi.a $(MPI_LIBS) $(EXTRA_LIBS
 EXTRA_DIST = BigCosmoVolume/makeIC.py \
 	     BigPerturbedBox/makeIC_fcc.py \
 	     CosmoVolume/cosmoVolume.yml CosmoVolume/getIC.sh CosmoVolume/run.sh \
+	     CoolingBox/coolingBox.yml CoolingBox/energy_plot.py CoolingBox/makeIC.py CoolingBox/run.sh \
 	     EAGLE_12/eagle_12.yml EAGLE_12/getIC.sh EAGLE_12/README EAGLE_12/run.sh \
 	     EAGLE_25/eagle_25.yml EAGLE_25/getIC.sh EAGLE_25/README EAGLE_25/run.sh \
 	     EAGLE_50/eagle_50.yml EAGLE_50/getIC.sh EAGLE_50/README EAGLE_50/run.sh \
diff --git a/examples/main.c b/examples/main.c
index c2d3e3b8b926b72b0e24519da87e4139657b0bfb..084a34723a7b1c3d6e075a67936c14b2380e69ad 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -58,6 +58,7 @@ void print_help_message() {
   printf("Valid options are:\n");
   printf("  %2s %8s %s\n", "-a", "", "Pin runners using processor affinity");
   printf("  %2s %8s %s\n", "-c", "", "Run with cosmological time integration");
+  printf("  %2s %8s %s\n", "-C", "", "Run with cooling");
   printf(
       "  %2s %8s %s\n", "-d", "",
       "Dry run. Read the parameter file, allocate memory but does not read ");
@@ -83,8 +84,8 @@ void print_help_message() {
   printf("  %2s %8s %s\n", "-t", "{int}",
          "The number of threads to use on each MPI rank. Defaults to 1 if not "
          "specified.");
-  printf("  %2s %8s %s\n", "-v", "[12]",
-         "Increase the level of verbosity 1: MPI-rank 0 writes ");
+  printf("  %2s %8s %s\n", "-v", "[12]", "Increase the level of verbosity");
+  printf("  %2s %8s %s\n", "", "", "1: MPI-rank 0 writes ");
   printf("  %2s %8s %s\n", "", "", "2: All MPI-ranks write");
   printf("  %2s %8s %s\n", "-y", "{int}",
          "Time-step frequency at which task graphs are dumped");
@@ -146,6 +147,7 @@ int main(int argc, char *argv[]) {
   int nsteps = -2;
   int with_cosmology = 0;
   int with_external_gravity = 0;
+  int with_cooling = 0;
   int with_self_gravity = 0;
   int with_hydro = 0;
   int with_fp_exceptions = 0;
@@ -157,13 +159,16 @@ int main(int argc, char *argv[]) {
 
   /* Parse the parameters */
   int c;
-  while ((c = getopt(argc, argv, "acdDef:gGhn:st:v:y:")) != -1) switch (c) {
+  while ((c = getopt(argc, argv, "acCdDef:gGhn:st:v:y:")) != -1) switch (c) {
       case 'a':
         with_aff = 1;
         break;
       case 'c':
         with_cosmology = 1;
         break;
+      case 'C':
+        with_cooling = 1;
+        break;
       case 'd':
         dry_run = 1;
         break;
@@ -325,12 +330,6 @@ int main(int argc, char *argv[]) {
   struct hydro_props hydro_properties;
   hydro_props_init(&hydro_properties, params);
 
-  /* Initialise the external potential properties */
-  struct external_potential potential;
-  if (with_external_gravity)
-    potential_init(params, &prog_const, &us, &potential);
-  if (with_external_gravity && myrank == 0) potential_print(&potential);
-
   /* Read particles and space information from (GADGET) ICs */
   char ICfileName[200] = "";
   parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
@@ -433,6 +432,17 @@ int main(int argc, char *argv[]) {
     message("nr of cells at depth %i is %i.", data[0], data[1]);
   }
 
+  /* Initialise the external potential properties */
+  struct external_potential potential;
+  if (with_external_gravity)
+    potential_init(params, &prog_const, &us, &potential);
+  if (with_external_gravity && myrank == 0) potential_print(&potential);
+
+  /* Initialise the cooling function properties */
+  struct cooling_data cooling;
+  if (with_cooling) cooling_init(params, &us, &prog_const, &cooling);
+  if (with_cooling && myrank == 0) cooling_print(&cooling);
+
   /* Construct the engine policy */
   int engine_policies = ENGINE_POLICY | engine_policy_steal;
   if (with_drift_all) engine_policies |= engine_policy_drift_all;
@@ -440,13 +450,14 @@ int main(int argc, char *argv[]) {
   if (with_self_gravity) engine_policies |= engine_policy_self_gravity;
   if (with_external_gravity) engine_policies |= engine_policy_external_gravity;
   if (with_cosmology) engine_policies |= engine_policy_cosmology;
+  if (with_cooling) engine_policies |= engine_policy_cooling;
 
   /* Initialize the engine with the space and policies. */
   if (myrank == 0) clocks_gettime(&tic);
   struct engine e;
   engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, with_aff,
               engine_policies, talking, &us, &prog_const, &hydro_properties,
-              &potential);
+              &potential, &cooling);
   if (myrank == 0) {
     clocks_gettime(&toc);
     message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
@@ -519,7 +530,7 @@ int main(int argc, char *argv[]) {
 
       /* Make sure output file is empty, only on one rank. */
       char dumpfile[30];
-      snprintf(dumpfile, 30, "thread_info_MPI-step%d.dat", j);
+      snprintf(dumpfile, 30, "thread_info_MPI-step%d.dat", j + 1);
       FILE *file_thread;
       if (myrank == 0) {
         file_thread = fopen(dumpfile, "w");
@@ -571,7 +582,7 @@ int main(int argc, char *argv[]) {
 
 #else
       char dumpfile[30];
-      snprintf(dumpfile, 30, "thread_info-step%d.dat", j);
+      snprintf(dumpfile, 30, "thread_info-step%d.dat", j + 1);
       FILE *file_thread;
       file_thread = fopen(dumpfile, "w");
       /* Add some information to help with the plots */
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index b1e106b94d4b42e51836ce1c36d733062746be6d..10bd09a18f7e8099e4db559e76957d19d8f90164 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -9,8 +9,8 @@ InternalUnitSystem:
 # Parameters for the task scheduling
 Scheduler:
   nr_queues:        0        # (Optional) The number of task queues to use. Use 0  to let the system decide.
-  cell_max_size:    8000000  # (Optional) Maximal number of interactions per task (this is the default value).
-  cell_sub_size:    8000000  # (Optional) Maximal number of interactions per sub-task  (this is the default value).
+  cell_max_size:    8000000  # (Optional) Maximal number of interactions per task if we force the split (this is the default value).
+  cell_sub_size:    64000000 # (Optional) Maximal number of interactions per sub-task  (this is the default value).
   cell_split_size:  400      # (Optional) Maximal number of particles per cell (this is the default value).
 
 # Parameters governing the time integration
@@ -63,7 +63,7 @@ DomainDecomposition:
   initial_grid_z:    10
   repartition_type:   b     # (Optional) The re-decomposition strategy ("n", "b", "v", "e" or "x").
  
-# Parameters related to external potentials
+# Parameters related to external potentials --------------------------------------------
   
 # Point mass external potentials
 PointMass:
@@ -73,6 +73,7 @@ PointMass:
   mass:            1e10     # mass of external point mass (internal units)
   timestep_mult:   0.03     # Dimensionless pre-factor for the time-step condition
 
+# Isothermal potential parameters
 IsothermalPotential:
   position_x:      100.     # Location of centre of isothermal potential (internal units)
   position_y:      100.
@@ -80,9 +81,27 @@ IsothermalPotential:
   vrot:            200.     # Rotation speed of isothermal potential (internal units)
   timestep_mult:   0.03     # Dimensionless pre-factor for the time-step condition
 
+# Disk-patch potential parameters
 Disk-PatchPotential:
   surface_density: 10.      # Surface density of the disk (internal units)
   scale_height:    100.     # Scale height of the disk (internal units)
   z_disk:          200.     # Disk height (internal units)
   timestep_mult:   0.03     # Dimensionless pre-factor for the time-step condition
   growth_time:     5.       # (Optional) Time for the disk to grow to its final size (multiple of the dynamical time)
+
+# Parameters related to cooling function  ----------------------------------------------
+
+# Constant du/dt cooling function
+ConstCooling:
+  cooling_rate: 1.          # Cooling rate (du/dt) (internal units)
+  min_energy:   1.          # Minimal internal energy per unit mass (internal units)
+  cooling_tstep_mult: 1.    # Dimensionless pre-factor for the time-step condition
+
+# Constant lambda cooling function
+LambdaCooling:
+  lambda:                      2.0   # Cooling rate (in cgs units)
+  minimum_temperature:         1.0e4 # Minimal temperature (Kelvin)
+  mean_molecular_weight:       0.59  # Mean molecular weight
+  hydrogen_mass_abundance:     0.75  # Hydrogen mass abundance (dimensionless)
+  cooling_tstep_mult:          1.0   # Dimensionless pre-factor for the time-step condition
+
diff --git a/examples/plot_scaling_results.py b/examples/plot_scaling_results.py
index cf9c96d41fe8f0897acb6be1f9af2269abee6174..5a76e9870bd3ec55807c7b79c475c62b14119e5c 100755
--- a/examples/plot_scaling_results.py
+++ b/examples/plot_scaling_results.py
@@ -48,7 +48,7 @@ threadList = []
 hexcols = ['#332288', '#88CCEE', '#44AA99', '#117733', '#999933', '#DDCC77',
            '#CC6677', '#882255', '#AA4499', '#661100', '#6699CC', '#AA4466',
            '#4477AA']
-linestyle = (hexcols[0],hexcols[1],hexcols[3],hexcols[5],hexcols[6])
+linestyle = (hexcols[0],hexcols[1],hexcols[3],hexcols[5],hexcols[6],hexcols[8])
 #cmdLine = './swift_fixdt -s -t 16 cosmoVolume.yml'
 #platform = 'KNL'
 
@@ -68,6 +68,9 @@ elif len(sys.argv) == 5:
 elif len(sys.argv) == 6:
   inputFileNames = (sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4],sys.argv[5])
   numOfSeries = 5
+elif len(sys.argv) == 7:
+  inputFileNames = (sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4],sys.argv[5],sys.argv[6])
+  numOfSeries = 6
 
 # Get the names of the branch, Git revision, hydro scheme and hydro kernel
 def parse_header(inputFile):
@@ -76,7 +79,8 @@ def parse_header(inputFile):
     for line in f:
       if 'Branch:' in line:
         s = line.split()
-        branch.append(s[2])
+        line = s[2:]
+        branch.append(" ".join(line))
       elif 'Revision:' in line:
         s = line.split() 
         revision.append(s[2])
@@ -129,8 +133,8 @@ def parse_files():
     parse_header(file_list[0])
     
     version.append(branch[i] + " " + revision[i] + "\n" + hydro_scheme[i] + 
-                   "\n" + hydro_kernel[i] + r", $N_{ngb}$=" + hydro_neighbours[i] + 
-                   r", $\eta$=" + hydro_eta[i])                  
+                   "\n" + hydro_kernel[i] + r", $N_{ngb}=%d$"%float(hydro_neighbours[i]) + 
+                   r", $\eta=%.3f$"%float(hydro_eta[i]))
     times.append([])
     totalTime.append([])
     speedUp.append([])
@@ -176,7 +180,7 @@ def print_results(times,totalTime,parallelEff,version):
 
 def plot_results(times,totalTime,speedUp,parallelEff):
   
-  fig, axarr = plt.subplots(2, 2, figsize=(10,10), frameon=False)
+  fig, axarr = plt.subplots(2, 2, figsize=(10,10), frameon=True)
   speedUpPlot = axarr[0, 0]
   parallelEffPlot = axarr[0, 1]
   totalTimePlot = axarr[1, 0]
@@ -218,7 +222,7 @@ def plot_results(times,totalTime,speedUp,parallelEff):
   totalTimePlot.set_xlim([0.9, 10**(np.floor(np.log10(threadList[i][-1]))+0.5)])
   totalTimePlot.set_ylim([10**np.floor(np.log10(np.min(totalTime)*0.6)), 1.2*10**np.floor(np.log10(np.max(totalTime) * 1.5)+1)])
   
-  totalTimePlot.legend(bbox_to_anchor=(1.14, 1), loc=2, borderaxespad=0.,prop={'size':12})
+  totalTimePlot.legend(bbox_to_anchor=(1.14, 0.97), loc=2, borderaxespad=0.,prop={'size':12}, frameon=False)
   emptyPlot.axis('off')
   
   for i, txt in enumerate(threadList[0]):
diff --git a/src/Makefile.am b/src/Makefile.am
index 430514528dcc919956b161d0e3879a1db1d8693c..2343ab99ffd90a27e588344c1fae4f1491b4625e 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -43,7 +43,8 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     engine.h swift.h serial_io.h timers.h debug.h scheduler.h proxy.h parallel_io.h \
     common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h parser.h \
     physical_constants.h physical_constants_cgs.h potentials.h version.h \
-    hydro_properties.h threadpool.h
+    hydro_properties.h threadpool.h cooling.h
+
 
 # Common source files
 AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
@@ -51,7 +52,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
     units.c common_io.c single_io.c multipole.c version.c map.c \
     kernel_hydro.c tools.c part.c partition.c clocks.c parser.c \
     physical_constants.c potentials.c hydro_properties.c \
-    runner_doiact_fft.c threadpool.c
+    runner_doiact_fft.c threadpool.c cooling.c
 
 # Include files for distribution, not installation.
 nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
@@ -71,7 +72,9 @@ nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_h
 		 hydro/Gizmo/hydro.h hydro/Gizmo/hydro_iact.h hydro/Gizmo/hydro_io.h \
                  hydro/Gizmo/hydro_debug.h hydro/Gizmo/hydro_part.h \
 	         riemann.h riemann/riemann_hllc.h riemann/riemann_trrs.h \
-		 riemann/riemann_exact.h riemann/riemann_vacuum.h 
+		 riemann/riemann_exact.h riemann/riemann_vacuum.h \
+	         cooling/const_du/cooling.h cooling/const_lambda/cooling.h
+
 
 # Sources and flags for regular library
 libswiftsim_la_SOURCES = $(AM_SOURCES)
diff --git a/src/cell.h b/src/cell.h
index 9f3971dd59c943ebd49a9ac1c812b612ab68bdaf..fd836206241c55cd6b9da2e29157c11a14c5a892 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -163,6 +163,9 @@ struct cell {
   /* Task for external gravity */
   struct task *grav_external;
 
+  /* Task for cooling */
+  struct task *cooling;
+
   /* Number of tasks that are associated with this cell. */
   int nr_tasks;
 
diff --git a/src/const.h b/src/const.h
index 5f7bf0d26431887c8e37d96b20b0eba460de55a6..552d49c8f1e4dd9f4fa2855e5fed548acd5c0b3f 100644
--- a/src/const.h
+++ b/src/const.h
@@ -94,6 +94,11 @@
 //#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
 //#define EXTERNAL_POTENTIAL_DISK_PATCH
 
+/* Cooling properties */
+//#define COOLING_CONST_DU
+#define COOLING_CONST_LAMBDA
+//#define COOLING_GRACKLE
+
 /* Are we debugging ? */
 //#define SWIFT_DEBUG_CHECKS
 
diff --git a/src/cooling.c b/src/cooling.c
new file mode 100644
index 0000000000000000000000000000000000000000..102adba9d521dbe53caeff3c4409292bb64a29b2
--- /dev/null
+++ b/src/cooling.c
@@ -0,0 +1,54 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* This object's header. */
+#include "cooling.h"
+
+/**
+ * @brief Initialises the cooling properties.
+ *
+ * Calls cooling_init_backend for the chosen cooling function.
+ *
+ * @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 cooling The cooling properties to initialize
+ */
+void cooling_init(const struct swift_params* parameter_file,
+                  const struct UnitSystem* us,
+                  const struct phys_const* phys_const,
+                  struct cooling_data* cooling) {
+
+  cooling_init_backend(parameter_file, us, phys_const, cooling);
+}
+
+/**
+ * @brief Prints the properties of the cooling model to stdout.
+ *
+ * Calls cooling_print_backend for the chosen cooling function.
+ *
+ * @param cooling The properties of the cooling function.
+ */
+void cooling_print(const struct cooling_data* cooling) {
+
+  cooling_print_backend(cooling);
+}
diff --git a/src/cooling.h b/src/cooling.h
new file mode 100644
index 0000000000000000000000000000000000000000..034ee2329d91b17b875932bc4eb40d0d80a05111
--- /dev/null
+++ b/src/cooling.h
@@ -0,0 +1,52 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_COOLING_H
+#define SWIFT_COOLING_H
+
+/**
+ * @file src/cooling.h
+ * @brief Branches between the different cooling functions.
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers. */
+#include "const.h"
+
+/* Import the right cooling definition */
+#if defined(COOLING_CONST_DU)
+#include "./cooling/const_du/cooling.h"
+#elif defined(COOLING_CONST_LAMBDA)
+#include "./cooling/const_lambda/cooling.h"
+#elif defined(COOLING_GRACKLE)
+#include "./cooling/grackle/cooling.h"
+#else
+#error "Invalid choice of cooling function."
+#endif
+
+/* Common functions */
+void cooling_init(const struct swift_params* parameter_file,
+                  const struct UnitSystem* us,
+                  const struct phys_const* phys_const,
+                  struct cooling_data* cooling);
+
+void cooling_print(const struct cooling_data* cooling);
+
+#endif /* SWIFT_COOLING_H */
diff --git a/src/cooling/const_du/cooling.h b/src/cooling/const_du/cooling.h
new file mode 100644
index 0000000000000000000000000000000000000000..634723f5a5a9af56aefc5fee1d409f583b243eba
--- /dev/null
+++ b/src/cooling/const_du/cooling.h
@@ -0,0 +1,139 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *                    Richard Bower (r.g.bower@durham.ac.uk)
+ *                    Stefan Arridge  (stefan.arridge@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_COOLING_CONST_DU_H
+#define SWIFT_COOLING_CONST_DU_H
+
+/**
+ * @file src/cooling/const/cooling.h
+ * @brief Routines related to the "constant cooling" cooling function.
+ *
+ * This is the simplest possible cooling function. A constant cooling rate with
+ * a minimal energy floor is applied.
+ */
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local includes. */
+#include "const.h"
+#include "error.h"
+#include "hydro.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/**
+ * @brief Properties of the cooling function.
+ */
+struct cooling_data {
+
+  /*! Cooling rate in internal units. du_dt = -cooling_rate */
+  float cooling_rate;
+
+  /*! Minimally allowed internal energy of the particles */
+  float min_energy;
+
+  /*! Constant multiplication factor for time-step criterion */
+  float cooling_tstep_mult;
+};
+
+/**
+ * @brief Apply the cooling function to a particle.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param cooling The #cooling_data used in the run.
+ * @param p Pointer to the particle data.
+ * @param dt The time-step of this particle.
+ */
+__attribute__((always_inline)) INLINE static void cooling_cool_part(
+    const struct phys_const* const phys_const, const struct UnitSystem* us,
+    const struct cooling_data* cooling, struct part* p, float dt) {
+
+  /* Get current internal energy (dt=0) */
+  const float u_old = hydro_get_internal_energy(p, 0.f);
+
+  /* Get cooling function properties */
+  const float du_dt = -cooling->cooling_rate;
+  const float u_floor = cooling->min_energy;
+
+  /* Constant cooling with a minimal floor */
+  float u_new;
+  if (u_old - du_dt * dt > u_floor) {
+    u_new = u_old + du_dt * dt;
+  } else {
+    u_new = u_floor;
+  }
+
+  /* Update the internal energy */
+  hydro_set_internal_energy(p, u_new);
+}
+
+/**
+ * @brief Computes the cooling time-step.
+ *
+ * @param cooling The #cooling_data used in the run.
+ * @param phys_const The physical constants in internal units.
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static double cooling_timestep(
+    const struct cooling_data* cooling,
+    const struct phys_const* const phys_const, const struct part* const p) {
+
+  const float cooling_rate = cooling->cooling_rate;
+  const float internal_energy =
+      hydro_get_internal_energy(p, 0);  // dt = 0 because using current entropy
+  return cooling->cooling_tstep_mult * internal_energy / cooling_rate;
+}
+
+/**
+ * @brief Initialises the cooling properties.
+ *
+ * @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 cooling The cooling properties to initialize
+ */
+static INLINE void cooling_init_backend(
+    const struct swift_params* parameter_file, const struct UnitSystem* us,
+    const struct phys_const* phys_const, struct cooling_data* cooling) {
+
+  cooling->cooling_rate =
+      parser_get_param_double(parameter_file, "ConstCooling:cooling_rate");
+  cooling->min_energy =
+      parser_get_param_double(parameter_file, "ConstCooling:min_energy");
+  cooling->cooling_tstep_mult = parser_get_param_double(
+      parameter_file, "ConstCooling:cooling_tstep_mult");
+}
+
+/**
+ * @brief Prints the properties of the cooling model to stdout.
+ *
+ * @param cooling The properties of the cooling function.
+ */
+static INLINE void cooling_print_backend(const struct cooling_data* cooling) {
+
+  message("Cooling function is 'Constant cooling' with rate %f and floor %f",
+          cooling->cooling_rate, cooling->min_energy);
+}
+
+#endif /* SWIFT_COOLING_CONST_DU_H */
diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
new file mode 100644
index 0000000000000000000000000000000000000000..3fd583b55a31534996e23e4dce71bcae9394d7de
--- /dev/null
+++ b/src/cooling/const_lambda/cooling.h
@@ -0,0 +1,212 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *                    Richard Bower (r.g.bower@durham.ac.uk)
+ *                    Stefan Arridge  (stefan.arridge@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_COOLING_CONST_LAMBDA_H
+#define SWIFT_COOLING_CONST_LAMBDA_H
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local includes. */
+#include "const.h"
+#include "error.h"
+#include "hydro.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/* Cooling Properties */
+struct cooling_data {
+
+  /*! Cooling rate in cgs units. Defined by 'rho * du/dt = -lambda * n_H^2'*/
+  float lambda;
+
+  /*! Minimum temperature (in Kelvin) for all gas particles*/
+  float min_temperature;
+
+  /*! Fraction of gas mass that is Hydrogen. Used to calculate n_H*/
+  float hydrogen_mass_abundance;
+
+  /* 'mu', used to convert min_temperature to min_internal energy*/
+  float mean_molecular_weight;
+
+  /*! Minimally allowed internal energy of the particles */
+  float min_energy;
+  float min_energy_cgs;
+
+  /*! Constant multiplication factor for time-step criterion */
+  float cooling_tstep_mult;
+};
+
+/**
+ * @brief Calculates du/dt in code units for a particle.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param cooling The #cooling_data used in the run.
+ * @param p Pointer to the particle data..
+ */
+__attribute__((always_inline)) INLINE static float cooling_rate(
+    const struct phys_const* const phys_const, const struct UnitSystem* us,
+    const struct cooling_data* cooling, const struct part* p) {
+
+  /* Get particle properties */
+  /* Density */
+  const float rho = hydro_get_density(p);
+  /* Get cooling function properties */
+  const float X_H = cooling->hydrogen_mass_abundance;
+  /* lambda should always be set in cgs units */
+  const float lambda_cgs = cooling->lambda;
+
+  /*convert from internal code units to cgs*/
+  const float rho_cgs =
+      rho * units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
+  const float m_p_cgs = phys_const->const_proton_mass *
+                        units_cgs_conversion_factor(us, UNIT_CONV_MASS);
+  const float n_H_cgs = X_H * rho_cgs / m_p_cgs;
+
+  /* Calculate du_dt */
+  const float du_dt_cgs = -lambda_cgs * n_H_cgs * n_H_cgs / rho_cgs;
+
+  /* Convert du/dt back to internal code units */
+  const float du_dt =
+      du_dt_cgs * units_cgs_conversion_factor(us, UNIT_CONV_TIME) /
+      units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+
+  return du_dt;
+}
+
+/**
+ * @brief Apply the cooling function to a particle.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param cooling The #cooling_data used in the run.
+ * @param p Pointer to the particle data.
+ * @param dt The time-step of this particle.
+ */
+__attribute__((always_inline)) INLINE static void cooling_cool_part(
+    const struct phys_const* const phys_const, const struct UnitSystem* us,
+    const struct cooling_data* cooling, struct part* p, float dt) {
+
+  /* Get current internal energy (dt=0) */
+  const float u_old = hydro_get_internal_energy(p, 0.f);
+
+  /* Internal energy floor */
+  const float u_floor = cooling->min_energy;
+
+  /* Calculate du_dt */
+  const float du_dt = cooling_rate(phys_const, us, cooling, p);
+
+  /* Intergrate cooling equation, but enforce energy floor */
+  float u_new;
+  if (u_old + du_dt * dt > u_floor) {
+    u_new = u_old + du_dt * dt;
+  } else {
+    u_new = u_floor;
+  }
+
+  /* Update the internal energy */
+  hydro_set_internal_energy(p, u_new);
+
+  /* if (-(u_new_test - u_old) / u_old > 1.0e-6) */
+  /*   error( */
+  /*       "Particle has not successfully cooled: u_old = %g , du_dt = %g , dt =
+   * " */
+  /*       "%g, du_dt*dt = %g, u_old + du_dt*dt = %g, u_new = %g\n", */
+  /*       u_old, du_dt, dt, du_dt * dt, u_new, u_new_test); */
+}
+
+/**
+ * @brief Computes the time-step due to cooling
+ *
+ * @param cooling The #cooling_data used in the run.
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float cooling_timestep(
+    const struct cooling_data* cooling,
+    const struct phys_const* const phys_const, const struct UnitSystem* us,
+    const struct part* const p) {
+
+  /* Get du_dt */
+  const float du_dt = cooling_rate(phys_const, us, cooling, p);
+
+  /* Get current internal energy (dt=0) */
+  const float u = hydro_get_internal_energy(p, 0.f);
+
+  return u / du_dt;
+}
+
+/**
+ * @brief Initialises the cooling properties.
+ *
+ * @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 cooling The cooling properties to initialize
+ */
+static INLINE void cooling_init_backend(
+    const struct swift_params* parameter_file, const struct UnitSystem* us,
+    const struct phys_const* phys_const, struct cooling_data* cooling) {
+
+  cooling->lambda =
+      parser_get_param_double(parameter_file, "LambdaCooling:lambda");
+  cooling->min_temperature = parser_get_param_double(
+      parameter_file, "LambdaCooling:minimum_temperature");
+  cooling->hydrogen_mass_abundance = parser_get_param_double(
+      parameter_file, "LambdaCooling:hydrogen_mass_abundance");
+  cooling->mean_molecular_weight = parser_get_param_double(
+      parameter_file, "LambdaCooling:mean_molecular_weight");
+  cooling->cooling_tstep_mult = parser_get_param_double(
+      parameter_file, "LambdaCooling:cooling_tstep_mult");
+
+  /*convert minimum temperature into minimum internal energy*/
+  const float u_floor =
+      phys_const->const_boltzmann_k * cooling->min_temperature /
+      (hydro_gamma_minus_one * cooling->mean_molecular_weight *
+       phys_const->const_proton_mass);
+  const float u_floor_cgs =
+      u_floor * units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+
+  cooling->min_energy = u_floor;
+  cooling->min_energy_cgs = u_floor_cgs;
+}
+
+/**
+ * @brief Prints the properties of the cooling model to stdout.
+ *
+ * @param cooling The properties of the cooling function.
+ */
+static INLINE void cooling_print_backend(const struct cooling_data* cooling) {
+
+  message(
+      "Cooling function is 'Constant lambda' with "
+      "(lambda,min_temperature,hydrogen_mass_abundance,mean_molecular_weight) "
+      "=  (%g,%g,%g,%g)",
+      cooling->lambda, cooling->min_temperature,
+      cooling->hydrogen_mass_abundance, cooling->mean_molecular_weight);
+}
+
+#endif /* SWIFT_COOLING_CONST_LAMBDA_H */
diff --git a/src/engine.c b/src/engine.c
index 9c0fc5d86ed19902dad21484f37e7af4fdeede90..b9f903237a36286c976cad23c6ffd98e39c13d70 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -68,7 +68,7 @@
 #include "units.h"
 #include "version.h"
 
-const char *engine_policy_names[14] = {"none",
+const char *engine_policy_names[15] = {"none",
                                        "rand",
                                        "steal",
                                        "keep",
@@ -81,7 +81,8 @@ const char *engine_policy_names[14] = {"none",
                                        "self_gravity",
                                        "external_gravity",
                                        "cosmology_integration",
-                                       "drift_all"};
+                                       "drift_all",
+                                       "cooling"};
 
 /** The rank of the engine as a global variable (for messages). */
 int engine_rank;
@@ -95,7 +96,6 @@ int engine_rank;
  *
  * @return The new #link pointer.
  */
-
 void engine_addlink(struct engine *e, struct link **l, struct task *t) {
 
   /* Get the next free link. */
@@ -185,6 +185,8 @@ void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c,
 
   struct scheduler *s = &e->sched;
   const int is_fixdt = (e->policy & engine_policy_fixdt) == engine_policy_fixdt;
+  const int is_with_cooling =
+      (e->policy & engine_policy_cooling) == engine_policy_cooling;
 
   /* Is this the super-cell? */
   if (super == NULL && (c->density != NULL || (c->count > 0 && !c->split))) {
@@ -220,6 +222,10 @@ void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c,
       c->extra_ghost = scheduler_addtask(s, task_type_extra_ghost,
                                          task_subtype_none, 0, 0, c, NULL, 0);
 #endif
+
+      if (is_with_cooling)
+        c->cooling = scheduler_addtask(s, task_type_cooling, task_subtype_none,
+                                       0, 0, c, NULL, 0);
     }
   }
 
@@ -773,7 +779,6 @@ void engine_addtasks_send(struct engine *e, struct cell *ci, struct cell *cj,
  * @param t_gradient The recv_gradient #task, if it has already been created.
  * @param t_ti The recv_ti #task, if required and has already been created.
  */
-
 void engine_addtasks_recv(struct engine *e, struct cell *c, struct task *t_xv,
                           struct task *t_rho, struct task *t_gradient,
                           struct task *t_ti) {
@@ -1584,7 +1589,7 @@ static inline void engine_make_hydro_loops_dependencies(struct scheduler *sched,
 void engine_make_extra_hydroloop_tasks(struct engine *e) {
 
   struct scheduler *sched = &e->sched;
-  int nr_tasks = sched->nr_tasks;
+  const int nr_tasks = sched->nr_tasks;
   const int nodeID = e->nodeID;
 
   for (int ind = 0; ind < nr_tasks; ind++) {
@@ -1780,11 +1785,18 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
       }
 #endif
     }
+
     /* External gravity tasks should depend on init and unlock the kick */
     else if (t->type == task_type_grav_external) {
       scheduler_addunlock(sched, t->ci->init, t);
       scheduler_addunlock(sched, t, t->ci->kick);
     }
+
+    /* Cooling tasks should depend on kick and does not unlock anything since
+     it is the last task*/
+    else if (t->type == task_type_cooling) {
+      scheduler_addunlock(sched, t->ci->kick, t);
+    }
   }
 }
 
@@ -1871,7 +1883,7 @@ void engine_maketasks(struct engine *e) {
   /* Count the number of tasks associated with each cell and
      store the density tasks in each cell, and make each sort
      depend on the sorts of its progeny. */
-  engine_count_and_link_tasks(e);
+  if (e->policy & engine_policy_hydro) engine_count_and_link_tasks(e);
 
   /* Append hierarchical tasks to each cells */
   if (e->policy & engine_policy_hydro)
@@ -1886,7 +1898,7 @@ void engine_maketasks(struct engine *e) {
   /* Run through the tasks and make force tasks for each density task.
      Each force task depends on the cell ghosts and unlocks the kick task
      of its super-cell. */
-  engine_make_extra_hydroloop_tasks(e);
+  if (e->policy & engine_policy_hydro) engine_make_extra_hydroloop_tasks(e);
 
   /* Add the dependencies for the self-gravity stuff */
   if (e->policy & engine_policy_self_gravity) engine_link_gravity_tasks(e);
@@ -2221,7 +2233,6 @@ void engine_print_task_counts(struct engine *e) {
  *
  * @param e The #engine.
  */
-
 void engine_rebuild(struct engine *e) {
 
   const ticks tic = getticks();
@@ -2790,6 +2801,11 @@ void engine_step(struct engine *e) {
     mask |= 1 << task_type_grav_external;
   }
 
+  /* Add the tasks corresponding to cooling to the masks */
+  if (e->policy & engine_policy_cooling) {
+    mask |= 1 << task_type_cooling;
+  }
+
   /* Add MPI tasks if need be */
   if (e->policy & engine_policy_mpi) {
 
@@ -3130,6 +3146,7 @@ void engine_unpin() {
  * @param physical_constants The #phys_const used for this run.
  * @param hydro The #hydro_props used for this run.
  * @param potential The properties of the external potential.
+ * @param cooling The properties of the cooling function.
  */
 void engine_init(struct engine *e, struct space *s,
                  const struct swift_params *params, int nr_nodes, int nodeID,
@@ -3137,7 +3154,8 @@ void engine_init(struct engine *e, struct space *s,
                  const struct UnitSystem *internal_units,
                  const struct phys_const *physical_constants,
                  const struct hydro_props *hydro,
-                 const struct external_potential *potential) {
+                 const struct external_potential *potential,
+                 const struct cooling_data *cooling) {
 
   /* Clean-up everything */
   bzero(e, sizeof(struct engine));
@@ -3189,6 +3207,7 @@ void engine_init(struct engine *e, struct space *s,
   e->physical_constants = physical_constants;
   e->hydro_properties = hydro;
   e->external_potential = potential;
+  e->cooling_data = cooling;
   e->parameter_file = params;
   engine_rank = nodeID;
 
diff --git a/src/engine.h b/src/engine.h
index 9f544b5b7b87731842a9945adb316ccb48d72542..c35163a44e1c6fa4e399fb421bfbfd0224b9486c 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -38,6 +38,7 @@
 
 /* Includes. */
 #include "clocks.h"
+#include "cooling.h"
 #include "parser.h"
 #include "partition.h"
 #include "potentials.h"
@@ -62,7 +63,8 @@ enum engine_policy {
   engine_policy_self_gravity = (1 << 9),
   engine_policy_external_gravity = (1 << 10),
   engine_policy_cosmology = (1 << 11),
-  engine_policy_drift_all = (1 << 12)
+  engine_policy_drift_all = (1 << 12),
+  engine_policy_cooling = (1 << 13),
 };
 
 extern const char *engine_policy_names[];
@@ -202,6 +204,9 @@ struct engine {
   /* Properties of external gravitational potential */
   const struct external_potential *external_potential;
 
+  /* Properties of the cooling scheme */
+  const struct cooling_data *cooling_data;
+
   /* The (parsed) parameter file */
   const struct swift_params *parameter_file;
 };
@@ -217,7 +222,8 @@ void engine_init(struct engine *e, struct space *s,
                  const struct UnitSystem *internal_units,
                  const struct phys_const *physical_constants,
                  const struct hydro_props *hydro,
-                 const struct external_potential *potential);
+                 const struct external_potential *potential,
+                 const struct cooling_data *cooling);
 void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
                    unsigned int submask);
 void engine_prepare(struct engine *e, int nodrift);
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 9dab5d7fd96a833ba9ea56889139c04000634645..37ce34e5910a033fab82dbc69839888a28c5ab12 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -20,6 +20,7 @@
 #include <float.h>
 #include "adiabatic_index.h"
 #include "approx_math.h"
+#include "equation_of_state.h"
 #include "hydro_gradients.h"
 #include "minmax.h"
 
@@ -502,3 +503,37 @@ __attribute__((always_inline)) INLINE static float hydro_get_density(
 
   return p->primitives.rho;
 }
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed internal
+ * energy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param u The new internal energy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
+    struct part* restrict p, float u) {
+
+  /* conserved.energy is NOT the specific energy (u), but the total thermal
+     energy (u*m) */
+  p->conserved.energy = u * p->conserved.mass;
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed entropy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param S The new entropy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_entropy(
+    struct part* restrict p, float S) {
+
+  p->conserved.energy = gas_internal_energy_from_entropy(p->primitives.rho, S) *
+                        p->conserved.mass;
+}
diff --git a/src/runner.c b/src/runner.c
index c0c6ea325b0122fa970fb7734333852078c69ffa..352a64664460343dd9004d60a9102b4875fce600 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -42,6 +42,7 @@
 #include "atomic.h"
 #include "cell.h"
 #include "const.h"
+#include "cooling.h"
 #include "debug.h"
 #include "drift.h"
 #include "engine.h"
@@ -158,6 +159,55 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
   if (timer) TIMER_TOC(timer_dograv_external);
 }
 
+/**
+ * @brief Calculate change in thermal state of particles induced
+ * by radiative cooling and heating.
+ *
+ * @param r runner task
+ * @param c cell
+ * @param timer 1 if the time is to be recorded.
+ */
+void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
+
+  struct part *restrict parts = c->parts;
+  const int count = c->count;
+  const int ti_current = r->e->ti_current;
+  const struct cooling_data *cooling = r->e->cooling_data;
+  const struct phys_const *constants = r->e->physical_constants;
+  const struct UnitSystem *us = r->e->internalUnits;
+  const double timeBase = r->e->timeBase;
+
+  TIMER_TIC;
+
+  /* Recurse? */
+  if (c->split) {
+    for (int k = 0; k < 8; k++)
+      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
+    return;
+  }
+
+#ifdef TASK_VERBOSE
+  OUT;
+#endif
+
+  /* Loop over the parts in this cell. */
+  for (int i = 0; i < count; i++) {
+
+    /* Get a direct pointer on the part. */
+    struct part *restrict p = &parts[i];
+
+    /* Kick has already updated ti_end, so need to check ti_begin */
+    if (p->ti_begin == ti_current) {
+
+      const double dt = (p->ti_end - p->ti_begin) * timeBase;
+
+      cooling_cool_part(constants, us, cooling, p, dt);
+    }
+  }
+
+  if (timer) TIMER_TOC(timer_do_cooling);
+}
+
 /**
  * @brief Sort the entries in ascending order using QuickSort.
  *
@@ -1281,6 +1331,9 @@ void *runner_main(void *data) {
         case task_type_grav_external:
           runner_do_grav_external(r, t->ci, 1);
           break;
+        case task_type_cooling:
+          runner_do_cooling(r, t->ci, 1);
+          break;
         default:
           error("Unknown task type.");
       }
diff --git a/src/runner.h b/src/runner.h
index bc9e202b6a4bc62b385db5f233b3fb5f2ef08c18..be19ab61b997f5730a04fa8c01c0787e8b99a8b2 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -53,6 +53,7 @@ void runner_do_sort(struct runner *r, struct cell *c, int flag, int clock);
 void runner_do_kick(struct runner *r, struct cell *c, int timer);
 void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer);
 void runner_do_init(struct runner *r, struct cell *c, int timer);
+void runner_do_cooling(struct runner *r, struct cell *c, int timer);
 void *runner_main(void *data);
 void runner_do_drift_mapper(void *map_data, int num_elements, void *extra_data);
 
diff --git a/src/space.h b/src/space.h
index 66d82c6f78a08447851a8dfdaea1231e5778693b..72b17405f13766ad2ccc9d53712068f28172067b 100644
--- a/src/space.h
+++ b/src/space.h
@@ -40,7 +40,7 @@
 #define space_cellallocchunk 1000
 #define space_splitsize_default 400
 #define space_maxsize_default 8000000
-#define space_subsize_default 8000000
+#define space_subsize_default 64000000
 #define space_stretch 1.10f
 #define space_maxreldx 0.25f
 
diff --git a/src/swift.h b/src/swift.h
index 7e3116c1de8bc8e6cc2f89d0d5cbe9771ffbf33a..80a4b1a1f792ceab3184d671f9aace277a36f419 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -27,6 +27,7 @@
 #include "cell.h"
 #include "clocks.h"
 #include "const.h"
+#include "cooling.h"
 #include "cycle.h"
 #include "debug.h"
 #include "engine.h"
diff --git a/src/task.c b/src/task.c
index 2cecce9f066727c9059b08764df647bd8f0f3901..4bc50cd0b3daac48edea7ea3040c04e8c42dbb02 100644
--- a/src/task.c
+++ b/src/task.c
@@ -46,11 +46,12 @@
 #include "inline.h"
 #include "lock.h"
 
+/* Task type names. */
 const char *taskID_names[task_type_count] = {
-    "none",       "sort",    "self",         "pair",          "sub_self",
-    "sub_pair",   "init",    "ghost",        "extra_ghost",   "kick",
-    "kick_fixdt", "send",    "recv",         "grav_gather_m", "grav_fft",
-    "grav_mm",    "grav_up", "grav_external"};
+    "none",       "sort",    "self",          "pair",          "sub_self",
+    "sub_pair",   "init",    "ghost",         "extra_ghost",   "kick",
+    "kick_fixdt", "send",    "recv",          "grav_gather_m", "grav_fft",
+    "grav_mm",    "grav_up", "grav_external", "cooling"};
 
 const char *subtaskID_names[task_subtype_count] = {
     "none", "density", "gradient", "force", "grav", "tend"};
@@ -116,6 +117,7 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
     case task_type_sort:
     case task_type_ghost:
     case task_type_extra_ghost:
+    case task_type_cooling:
       return task_action_part;
       break;
 
diff --git a/src/task.h b/src/task.h
index 5f26ec05b026518de1234162d1f5f1d195672507..f070451fe4e79e0c16dc3dcca1ce145c08841c47 100644
--- a/src/task.h
+++ b/src/task.h
@@ -53,6 +53,7 @@ enum task_types {
   task_type_grav_mm,
   task_type_grav_up,
   task_type_grav_external,
+  task_type_cooling,
   task_type_count
 } __attribute__((packed));
 
diff --git a/src/timers.h b/src/timers.h
index 0bf1f3f413d123c84cc30edff73ca9dfce4ea159..b93a34df4d90251d4686afaef0be51b36dc9a25e 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -61,6 +61,7 @@ enum {
   timer_qsteal,
   timer_runners,
   timer_step,
+  timer_do_cooling,
   timer_count,
 };