diff --git a/examples/CoolingInput/check_output.py b/examples/CoolingInput/check_output.py
new file mode 100644
index 0000000000000000000000000000000000000000..a65a83b4c73dac741a227652b0e706246ce66bac
--- /dev/null
+++ b/examples/CoolingInput/check_output.py
@@ -0,0 +1,17 @@
+#!/usr/bin/env python3
+import h5py as h5
+import sys
+
+# First snapshot
+snap_filename = "coolingBox_0000.hdf5"
+
+# Read the initial state of the gas
+f = h5.File(snap_filename,'r')
+he = f["/PartType0/HeDensity"][:]
+
+if (he != 5.).any():
+ print("Test Failed")
+ print("Received", he)
+
+else:
+ print("Test Succeed")
diff --git a/examples/CoolingInput/coolingBox.yml b/examples/CoolingInput/coolingBox.yml
new file mode 100644
index 0000000000000000000000000000000000000000..3d34e1fd9e1b58c9b292378be673d7e176644b7a
--- /dev/null
+++ b/examples/CoolingInput/coolingBox.yml
@@ -0,0 +1,48 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 2.0e33 # Solar masses
+ UnitLength_in_cgs: 3.0857e21 # Kiloparsecs
+ UnitVelocity_in_cgs: 1.0e5 # Kilometers per second
+ 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: 1e-5 # The end time of the simulation (in internal units).
+ dt_min: 1e-5 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-5 # 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: 1e-2 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-3 # 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).
+ 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_cgs: 1.0e-22 # 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
+
+# Cooling with Grackle 2.0
+GrackleCooling:
+ GrackleCloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
+ UVbackground: 0 # Enable or not the UV background
+ GrackleRedshift: 0 # Redshift to use (-1 means time based redshift)
+ GrackleHSShieldingDensityThreshold: 1.1708e-26 # self shielding threshold in internal system of units
diff --git a/examples/CoolingInput/getCoolingTable.sh b/examples/CoolingInput/getCoolingTable.sh
new file mode 100644
index 0000000000000000000000000000000000000000..809958bfc236e5ab0f7be924c62789fa13b3ac29
--- /dev/null
+++ b/examples/CoolingInput/getCoolingTable.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5
diff --git a/examples/CoolingInput/getGlass.sh b/examples/CoolingInput/getGlass.sh
new file mode 100644
index 0000000000000000000000000000000000000000..ffd92e88deae6e91237059adac2a6c2067caee46
--- /dev/null
+++ b/examples/CoolingInput/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
diff --git a/examples/CoolingInput/makeIC.py b/examples/CoolingInput/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..7452c86838a88fc1fa806ada656a08d4a358d27e
--- /dev/null
+++ b/examples/CoolingInput/makeIC.py
@@ -0,0 +1,107 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Stefan Arridge (stefan.arridge@durhama.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 with a constant density and pressure
+
+# Parameters
+periodic= 1 # 1 For periodic box
+boxSize = 1 # 1 kiloparsec
+rho = 3.2e3 # Density in code units (3.2e6 is 0.1 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"
+
+#---------------------------------------------------
+
+# Read id, position and h from glass
+glass = h5py.File("glassCube_32.hdf5", "r")
+ids = glass["/PartType0/ParticleIDs"][:]
+pos = glass["/PartType0/Coordinates"][:,:] * boxSize
+h = glass["/PartType0/SmoothingLength"][:] * boxSize
+
+# Compute basic properties
+numPart = size(pos) / 3
+mass = boxSize**3 * rho / numPart
+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.0857e21
+grp.attrs["Unit mass in cgs (U_M)"] = 2.0e33
+grp.attrs["Unit time in cgs (U_t)"] = 3.0857e16
+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
+
+m = full((numPart, 1), mass)
+ds = grp.create_dataset('Masses', (numPart,1), 'f')
+ds[()] = m
+
+h = reshape(h, (numPart, 1))
+ds = grp.create_dataset('SmoothingLength', (numPart, 1), 'f')
+ds[()] = h
+
+u = full((numPart, 1), internalEnergy)
+ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
+ds[()] = u
+
+ids = reshape(ids, (numPart, 1))
+ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
+ds[()] = ids
+
+ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = pos
+
+he = 5. * ones((numPart, 1))
+ds = grp.create_dataset('HeDensity', (numPart, 1), 'f')
+ds[()] = he
+
+file.close()
+
+print numPart
diff --git a/examples/CoolingInput/run.sh b/examples/CoolingInput/run.sh
new file mode 100644
index 0000000000000000000000000000000000000000..d952144be7d7d29bd5caad1271a8e56c851edba5
--- /dev/null
+++ b/examples/CoolingInput/run.sh
@@ -0,0 +1,27 @@
+
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e glassCube_32.hdf5 ]
+then
+ echo "Fetching initial glass file for the cooling box example..."
+ ./getGlass.sh
+fi
+if [ ! -e coolingBox.hdf5 ]
+then
+ echo "Generating initial conditions for the cooling box example..."
+ python makeIC.py
+fi
+
+# Get the Grackle cooling table
+if [ ! -e CloudyData_UVB=HM2012.h5 ]
+then
+ echo "Fetching the Cloudy tables required by Grackle..."
+ ./getCoolingTable.sh
+fi
+
+# Run SWIFT
+../swift -s -C -t 1 coolingBox.yml
+
+# Check energy conservation and cooling rate
+python check_output.py