diff --git a/configure.ac b/configure.ac
index 238e29468dd1a06d054b44d6d78b73afb3555407..ce3844381556ee99bc4a60915b1cdafffc6a1231 100644
--- a/configure.ac
+++ b/configure.ac
@@ -25,7 +25,7 @@ AC_DEFINE([PACKAGE_URL],["www.swiftsim.com"], [Package web pages])
 AC_COPYRIGHT
 AC_CONFIG_SRCDIR([src/space.c])
 AC_CONFIG_AUX_DIR([.])
-AM_INIT_AUTOMAKE
+AM_INIT_AUTOMAKE([subdir-objects])
 
 # Add local macro collection.
 AC_CONFIG_MACRO_DIR([m4])
@@ -421,6 +421,43 @@ AC_SUBST([METIS_LIBS])
 AC_SUBST([METIS_INCS])
 AM_CONDITIONAL([HAVEMETIS],[test -n "$METIS_LIBS"])
 
+# Check for grackle. 
+have_grackle="no"
+AC_ARG_WITH([grackle],
+    [AS_HELP_STRING([--with-grackle=PATH],
+       [root directory where grackle is installed @<:@yes/no@:>@]
+    )],
+    [],
+    [with_grackle="no"]
+)
+if test "x$with_grackle" != "xno"; then
+   AC_PROG_FC
+   AC_FC_LIBRARY_LDFLAGS
+   if test "x$with_grackle" != "xyes" -a "x$with_grackle" != "x"; then
+      GRACKLE_LIBS="-L$with_grackle/lib -lgrackle"
+      GRACKLE_INCS="-I$with_grackle/include"
+   else
+      GRACKLE_LIBS="-lgrackle"
+      GRACKLE_INCS=""
+   fi
+   
+   have_grackle="yes"
+      
+   AC_CHECK_LIB(
+      [grackle],
+      [initialize_grackle],
+      [AC_DEFINE([HAVE_GRACKLE],1,[The GRACKLE library appears to be present.]),
+        AC_DEFINE([CONFIG_BFLOAT_8],1,[Use doubles in grackle])
+      ],
+      [AC_MSG_ERROR(Cannot find grackle library!)],
+      [$GRACKLE_LIBS $GRACKLE_INCS $FCLIBS]
+   )
+fi
+AC_SUBST([GRACKLE_LIBS])
+AC_SUBST([GRACKLE_INCS])
+AM_CONDITIONAL([HAVEGRACKLE],[test -n "$GRACKLE_LIBS"])
+
+
 #  Check for tcmalloc a fast malloc that is part of the gperftools.
 have_tcmalloc="no"
 AC_ARG_WITH([tcmalloc],
@@ -948,6 +985,7 @@ AC_MSG_RESULT([
    Metis enabled    : $have_metis
    FFTW3 enabled    : $have_fftw3
    libNUMA enabled  : $have_numa
+   GRACKLE enabled  : $have_grackle
    Using tcmalloc   : $have_tcmalloc
    Using jemalloc   : $have_jemalloc
    CPU profiler     : $have_profiler
diff --git a/examples/CoolingBox/coolingBox.yml b/examples/CoolingBox/coolingBox.yml
index 1dbf06dec7c5ce2060b3dbb9be140fbb639a2a7a..ea84eeac27603d0c5d8185c44ecf1c0374ca7c76 100644
--- a/examples/CoolingBox/coolingBox.yml
+++ b/examples/CoolingBox/coolingBox.yml
@@ -39,3 +39,10 @@ LambdaCooling:
   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/CoolingBox/getCoolingTable.sh b/examples/CoolingBox/getCoolingTable.sh
new file mode 100755
index 0000000000000000000000000000000000000000..809958bfc236e5ab0f7be924c62789fa13b3ac29
--- /dev/null
+++ b/examples/CoolingBox/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/CoolingBox/run.sh b/examples/CoolingBox/run.sh
index 19e787df716145c1f5aa7744f4c7204c1c7f1064..61106d1f4819145ad4fb3a017918ca0e074a87c2 100755
--- a/examples/CoolingBox/run.sh
+++ b/examples/CoolingBox/run.sh
@@ -13,6 +13,13 @@ then
     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
 
diff --git a/examples/Makefile.am b/examples/Makefile.am
index 496d2c005004a1a68276fafb17e618c4b3ba9f8a..bddb32943781e9326e05cc5c55a87086fc96cf02 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -24,7 +24,7 @@ AM_CFLAGS = -I$(top_srcdir)/src $(HDF5_CPPFLAGS)
 AM_LDFLAGS = $(HDF5_LDFLAGS)
 
 # Extra libraries.
-EXTRA_LIBS = $(HDF5_LIBS) $(FFTW_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS)
+EXTRA_LIBS = $(HDF5_LIBS) $(FFTW_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(GRACKLE_LIBS)
 
 # MPI libraries.
 MPI_LIBS = $(METIS_LIBS) $(MPI_THREAD_LIBS)
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 948731d7879ec3e34f6e53f58dc72f6d7a6145ec..8975981fd0b3a4a1dab8e2daf278cdd40e2ddbf0 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -143,3 +143,9 @@ LambdaCooling:
   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: 1 # 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/src/Makefile.am b/src/Makefile.am
index b820e8060778fdf8e7313d1b47596e370927d274..c50c831b3e27ab07ade0c83ccb5cb3b4f78f6724 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -25,7 +25,7 @@ AM_LDFLAGS = $(HDF5_LDFLAGS) $(FFTW_LIBS) -version-info 0:0:0
 GIT_CMD = @GIT_CMD@
 
 # Additional dependencies for shared libraries.
-EXTRA_LIBS = $(HDF5_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS)
+EXTRA_LIBS = $(HDF5_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(GRACKLE_LIB)
 
 # MPI libraries.
 MPI_LIBS = $(METIS_LIBS) $(MPI_THREAD_LIBS)
@@ -48,6 +48,11 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
     gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h
 
+GRACKLE_SRC =
+if HAVEGRACKLE
+GRACKLE_SRC += cooling/grackle/grackle_wrapper.c
+endif
+
 # Common source files
 AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
     serial_io.c timers.c debug.c scheduler.c proxy.c parallel_io.c \
@@ -57,7 +62,9 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
     runner_doiact_fft.c threadpool.c cooling.c sourceterms.c \
     statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
     part_type.c xmf.c gravity_properties.c gravity.c \
-    collectgroup.c hydro_space.c equation_of_state.c
+    collectgroup.c hydro_space.c equation_of_state.c \
+    $(GRACKLE_SRC)
+
 
 # Include files for distribution, not installation.
 nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
@@ -115,6 +122,8 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
 		 cooling/none/cooling.h cooling/none/cooling_struct.h \
 	         cooling/const_du/cooling.h cooling/const_du/cooling_struct.h \
                  cooling/const_lambda/cooling.h cooling/const_lambda/cooling_struct.h \
+                 cooling/grackle/cooling.h cooling/grackle/cooling_struct.h \
+		 cooling/grackle/grackle_wrapper.h \
                  memswap.h dump.h logger.h
 
 
@@ -122,12 +131,14 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
 libswiftsim_la_SOURCES = $(AM_SOURCES)
 libswiftsim_la_CFLAGS = $(AM_CFLAGS)
 libswiftsim_la_LDFLAGS = $(AM_LDFLAGS) $(EXTRA_LIBS)
+libswiftsim_la_LIBADD = $(GRACKLE_LIBS)
 
 # Sources and flags for MPI library
 libswiftsim_mpi_la_SOURCES = $(AM_SOURCES)
 libswiftsim_mpi_la_CFLAGS = $(AM_CFLAGS) $(MPI_FLAGS)
 libswiftsim_mpi_la_LDFLAGS = $(AM_LDFLAGS) $(MPI_LIBS) $(EXTRA_LIBS)
 libswiftsim_mpi_la_SHORTNAME = mpi
+libswiftsim_mpi_la_LIBADD = $(GRACKLE_LIBS)
 
 
 # Versioning. If any sources change then update the version_string.h file with
diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h
new file mode 100644
index 0000000000000000000000000000000000000000..50700071904be662609e513b34a0215d0ef61d7f
--- /dev/null
+++ b/src/cooling/grackle/cooling.h
@@ -0,0 +1,240 @@
+/*******************************************************************************
+ * 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_GRACKLE_H
+#define SWIFT_COOLING_GRACKLE_H
+
+/**
+ * @file src/cooling/none/cooling.h
+ * @brief Empty infrastructure for the cases without cooling function
+ */
+
+/* Some standard headers. */
+#include <float.h>
+#include <math.h>
+
+/* Local includes. */
+#include "error.h"
+#include "hydro.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/* include the grackle wrapper */
+#include "grackle_wrapper.h"
+
+/**
+ * @brief Compute the cooling rate
+ *
+ * We do nothing.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param cooling The #cooling_function_data used in the run.
+ * @param p Pointer to the particle data.
+ * @param dt The time-step of this particle.
+ *
+ * @return du / dt
+ */
+__attribute__((always_inline)) INLINE static double cooling_rate(
+    const struct phys_const* restrict phys_const,
+    const struct unit_system* restrict us,
+    const struct cooling_function_data* restrict cooling,
+    struct part* restrict p, float dt) {
+
+  if (cooling->GrackleRedshift == -1) error("TODO time dependant redshift");
+
+  /* Get current internal energy (dt=0) */
+  double u_old = hydro_get_internal_energy(p);
+  double u_new = u_old;
+  /* Get current density */
+  const float rho = hydro_get_density(p);
+  /* Actual scaling fractor */
+  const float a_now = 1. / (1. + cooling->GrackleRedshift);
+
+  /* 0.02041 (= 1 Zsun in Grackle v2.0, but = 1.5761 Zsun in
+     Grackle v2.1) */
+  double Z = 0.02041;
+
+  if (wrap_do_cooling(rho, &u_new, dt, Z, a_now) == 0) {
+    error("Error in do_cooling.\n");
+    return 0;
+  }
+
+  return (u_new - u_old) / dt;
+}
+
+/**
+ * @brief Apply the cooling function to a particle.
+ *
+ * We do nothing.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param cooling The #cooling_function_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* restrict phys_const,
+    const struct unit_system* restrict us,
+    const struct cooling_function_data* restrict cooling,
+    struct part* restrict p, struct xpart* restrict xp, float dt) {
+
+  if (dt == 0.) return;
+
+  /* Current du_dt */
+  const float hydro_du_dt = hydro_get_internal_energy_dt(p);
+
+  float du_dt;
+  float delta_u;
+
+  du_dt = cooling_rate(phys_const, us, cooling, p, dt);
+
+  delta_u = du_dt * dt;
+
+  /* record energy lost */
+  xp->cooling_data.radiated_energy += -delta_u * hydro_get_mass(p);
+
+  /* Update the internal energy */
+  hydro_set_internal_energy_dt(p, hydro_du_dt + du_dt);
+}
+
+/**
+ * @brief Computes the cooling time-step.
+ *
+ * We return FLT_MAX so as to impose no limit on the time-step.
+ *
+ * @param cooling The #cooling_function_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_function_data* restrict cooling,
+    const struct phys_const* restrict phys_const,
+    const struct unit_system* restrict us, const struct part* restrict p) {
+
+  return FLT_MAX;
+}
+
+/**
+ * @brief Sets the cooling properties of the (x-)particles to a valid start
+ * state.
+ *
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void cooling_init_part(
+    const struct part* restrict p, struct xpart* restrict xp) {
+
+  xp->cooling_data.radiated_energy = 0.f;
+}
+
+/**
+ * @brief Returns the total radiated energy by this particle.
+ *
+ * @param xp The extended particle data
+ */
+__attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
+    const struct xpart* restrict xp) {
+
+  return xp->cooling_data.radiated_energy;
+}
+
+/**
+ * @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 unit_system* us,
+    const struct phys_const* phys_const,
+    struct cooling_function_data* cooling) {
+
+  double units_density, units_length, units_time;
+  int grackle_chemistry;
+  int UVbackground;
+
+  parser_get_param_string(parameter_file, "GrackleCooling:GrackleCloudyTable",
+                          cooling->GrackleCloudyTable);
+  cooling->UVbackground =
+      parser_get_param_int(parameter_file, "GrackleCooling:UVbackground");
+  cooling->GrackleRedshift =
+      parser_get_param_double(parameter_file, "GrackleCooling:GrackleRedshift");
+  cooling->GrackleHSShieldingDensityThreshold = parser_get_param_double(
+      parameter_file, "GrackleCooling:GrackleHSShieldingDensityThreshold");
+
+  UVbackground = cooling->UVbackground;
+  grackle_chemistry = 0; /* forced to be zero : read table */
+
+  units_density = us->UnitMass_in_cgs / pow(us->UnitLength_in_cgs, 3);
+  units_length = us->UnitLength_in_cgs;
+  units_time = us->UnitTime_in_cgs;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  float threshold = cooling->GrackleHSShieldingDensityThreshold;
+
+  threshold /= phys_const->const_proton_mass;
+  threshold /= pow(us->UnitLength_in_cgs, 3);
+
+  message("***************************************");
+  message("initializing grackle cooling function");
+  message("");
+  message("CloudyTable                        = %s",
+          cooling->GrackleCloudyTable);
+  message("UVbackground                       = %d", UVbackground);
+  message("GrackleRedshift                    = %g", cooling->GrackleRedshift);
+  message("GrackleHSShieldingDensityThreshold = %g atom/cm3", threshold);
+#endif
+
+  if (wrap_init_cooling(cooling->GrackleCloudyTable, UVbackground,
+                        units_density, units_length, units_time,
+                        grackle_chemistry) != 1) {
+    error("Error in initialize_chemistry_data.");
+  }
+
+#ifdef SWIFT_DEBUG_CHECKS
+  grackle_print_data();
+  message("");
+  message("***************************************");
+#endif
+}
+
+/**
+ * @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_function_data* cooling) {
+
+  message("Cooling function is 'Grackle'.");
+  message("CloudyTable                        = %s",
+          cooling->GrackleCloudyTable);
+  message("UVbackground                       = %d", cooling->UVbackground);
+  message("GrackleRedshift                    = %g", cooling->GrackleRedshift);
+  message("GrackleHSShieldingDensityThreshold = %g atom/cm3",
+          cooling->GrackleHSShieldingDensityThreshold);
+}
+
+#endif /* SWIFT_COOLING_GRACKLE_H */
diff --git a/src/cooling/grackle/cooling_struct.h b/src/cooling/grackle/cooling_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..b59af375cba0ed35fcb0159cb63605af8d94f030
--- /dev/null
+++ b/src/cooling/grackle/cooling_struct.h
@@ -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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_COOLING_STRUCT_NONE_H
+#define SWIFT_COOLING_STRUCT_NONE_H
+
+/**
+ * @file src/cooling/none/cooling_struct.h
+ * @brief Empty infrastructure for the cases without cooling function
+ */
+
+/**
+ * @brief Properties of the cooling function.
+ */
+struct cooling_function_data {
+
+  /* Filename of the Cloudy Table */
+  char GrackleCloudyTable[200];
+
+  /* Enable/Disable UV backgroud */
+  int UVbackground;
+
+  /* Redshift to use for the UV backgroud (-1 to use cosmological one) */
+  double GrackleRedshift;
+
+  /* Density Threshold for the shielding */
+  double GrackleHSShieldingDensityThreshold;
+};
+
+/**
+ * @brief Properties of the cooling stored in the particle data
+ */
+struct cooling_xpart_data {
+
+  /*! Energy radiated away by this particle since the start of the run */
+  float radiated_energy;
+};
+
+#endif /* SWIFT_COOLING_STRUCT_NONE_H */
diff --git a/src/cooling/grackle/grackle_wrapper.c b/src/cooling/grackle/grackle_wrapper.c
new file mode 100644
index 0000000000000000000000000000000000000000..9974717b670c937811cea54349124d9a38bed476
--- /dev/null
+++ b/src/cooling/grackle/grackle_wrapper.c
@@ -0,0 +1,227 @@
+/***********************************************************************
+/
+/ Grackle c wrapper
+/
+/
+/ Copyright (c) 2013, Enzo/Grackle Development Team.
+/
+/ Distributed under the terms of the Enzo Public Licence.
+/
+/ The full license is in the file LICENSE, distributed with this
+/ software.
+************************************************************************/
+
+#include "grackle_wrapper.h"
+
+#ifdef SWIFT_DEBUG_CHECKS
+#include <assert.h>
+#define GRACKLE_ASSERT(X) assert((X))
+#else
+#define GRACKLE_ASSERT(X)
+#endif
+
+code_units my_units;
+
+// arrays passed to grackle as input and to be filled
+#define FIELD_SIZE 1
+
+gr_float HDI_density[FIELD_SIZE];
+
+// Set grid dimension and size.
+// grid_start and grid_end are used to ignore ghost zones.
+const int grid_rank = 1;
+
+int wrap_init_cooling(char *CloudyTable, int UVbackground, double udensity,
+                      double ulength, double utime, int grackle_chemistry) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  grackle_verbose = 1;
+#endif
+  double velocity_units;
+
+  // First, set up the units system.
+  // These are conversions from code units to cgs.
+
+  my_units.a_units = 1.0;  // units for the expansion factor (1/1+zi)
+
+  my_units.comoving_coordinates =
+      0; /* so, according to the doc, we assume here all physical quantities to
+            be in proper coordinate (not comobile)  */
+  my_units.density_units = udensity;
+  my_units.length_units = ulength;
+  my_units.time_units = utime;
+  velocity_units =
+      my_units.a_units * my_units.length_units / my_units.time_units;
+  my_units.velocity_units = velocity_units;
+
+  // Second, create a chemistry object for parameters and rate data.
+  if (set_default_chemistry_parameters() == 0) {
+    error("Error in set_default_chemistry_parameters.");
+  }
+
+  // Set parameter values for chemistry.
+  grackle_data.use_grackle = 1;
+  grackle_data.with_radiative_cooling = 1;
+  grackle_data.primordial_chemistry =
+      grackle_chemistry;           // molecular network with H, He, D
+  grackle_data.metal_cooling = 1;  // metal cooling on
+  grackle_data.UVbackground = UVbackground;
+  grackle_data.grackle_data_file = CloudyTable;
+
+  // Finally, initialize the chemistry object.
+  // snl: a_value is not the true initial ae!! This should get set during
+  // update_UVbackground
+  double initial_redshift = 0.;
+  double a_value = 1. / (1. + initial_redshift);
+
+  // Finally, initialize the chemistry object.
+  if (initialize_chemistry_data(&my_units, a_value) == 0) {
+    error("Error in initialize_chemistry_data.");
+  }
+
+  return 1;
+}
+
+int wrap_set_UVbackground_on() {
+  // The UV background rates is enabled
+  grackle_data.UVbackground = 1;
+  return 1;
+}
+
+int wrap_set_UVbackground_off() {
+  // The UV background rates is disabled
+  grackle_data.UVbackground = 0;
+  return 1;
+}
+
+int wrap_get_cooling_time(double rho, double u, double Z, double a_now,
+                          double *coolingtime) {
+  gr_float den_factor = 1.0;
+  gr_float u_factor = 1.0;
+
+  gr_float x_velocity[FIELD_SIZE] = {0.0};
+  gr_float y_velocity[FIELD_SIZE] = {0.0};
+  gr_float z_velocity[FIELD_SIZE] = {0.0};
+
+  gr_float density[FIELD_SIZE] = {rho * den_factor};
+  gr_float metal_density[FIELD_SIZE] = {Z * density[0]};
+  gr_float energy[FIELD_SIZE] = {u * u_factor};
+
+  gr_float cooling_time[FIELD_SIZE] = {0.0};
+
+  int grid_dimension[3] = {1, 0, 0};
+  int grid_start[3] = {0, 0, 0};
+  int grid_end[3] = {0, 0, 0};
+
+  if (FIELD_SIZE != 1) {
+    error("field_size must currently be set to 1.");
+  }
+
+  if (calculate_cooling_time_table(&my_units, a_now, grid_rank, grid_dimension,
+                                   grid_start, grid_end, density, energy,
+                                   x_velocity, y_velocity, z_velocity,
+                                   metal_density, cooling_time) == 0) {
+    error("Error in calculate_cooling_time.");
+  }
+
+  // return updated chemistry and energy
+  for (int i = 0; i < FIELD_SIZE; i++) {
+    coolingtime[i] = cooling_time[i];
+  }
+
+  return 1;
+}
+
+int wrap_do_cooling(double rho, double *u, double dt, double Z, double a_now) {
+
+  GRACKLE_ASSERT(FIELD_SIZE == 1);
+
+  gr_float den_factor = 1.0;
+  gr_float u_factor = 1.0;
+
+  gr_float x_velocity[FIELD_SIZE] = {0.0};
+  gr_float y_velocity[FIELD_SIZE] = {0.0};
+  gr_float z_velocity[FIELD_SIZE] = {0.0};
+
+  gr_float density[FIELD_SIZE] = {rho * den_factor};
+  gr_float metal_density[FIELD_SIZE] = {Z * density[0]};
+  gr_float energy[FIELD_SIZE] = {(*u) * u_factor};
+
+  int grid_dimension[3] = {1, 0, 0};
+  int grid_start[3] = {0, 0, 0};
+  int grid_end[3] = {0, 0, 0};
+
+  if (solve_chemistry_table(&my_units, a_now, dt, grid_rank, grid_dimension,
+                            grid_start, grid_end, density, energy, x_velocity,
+                            y_velocity, z_velocity, metal_density) == 0) {
+    error("Error in solve_chemistry.");
+    return 0;
+  }
+  // return updated chemistry and energy
+  for (int i = 0; i < FIELD_SIZE; i++) {
+    u[i] = energy[i] / u_factor;
+  }
+
+  return 1;
+}
+
+void grackle_print_data() {
+  message("Grackle Data:");
+  message("\t Data file: %s", grackle_data.grackle_data_file);
+  message("\t With grackle: %i", grackle_data.use_grackle);
+  message("\t With radiative cooling: %i", grackle_data.with_radiative_cooling);
+  message("\t With UV background: %i", grackle_data.UVbackground);
+  message("\t With primordial chemistry: %i",
+          grackle_data.primordial_chemistry);
+  message("\t Number temperature bins: %i",
+          grackle_data.NumberOfTemperatureBins);
+  message("\t T = (%g, ..., %g)", grackle_data.TemperatureStart,
+          grackle_data.TemperatureEnd);
+
+  message("Primordial Cloudy");
+  cloudy_print_data(grackle_data.cloudy_primordial, 1);
+  if (grackle_data.metal_cooling) {
+    message("Metal Cooling");
+    cloudy_print_data(grackle_data.cloudy_metal, 0);
+  }
+
+  message("\t Gamma: %g", grackle_data.Gamma);
+
+  /* UVB */
+  if (grackle_data.UVbackground && grackle_data.primordial_chemistry != 0) {
+    struct UVBtable uvb = grackle_data.UVbackground_table;
+    long long N = uvb.Nz;
+    message("\t UV Background");
+    message("\t\t Redshift from %g to %g with %lli steps", uvb.zmin, uvb.zmax,
+            N);
+    message("\t\t z = (%g, ..., %g)", uvb.z[0], uvb.z[N - 1]);
+  }
+}
+
+void cloudy_print_data(const cloudy_data c, const int print_mmw) {
+  long long N = c.data_size;
+  message("\t Data size: %lli", N);
+  message("\t Grid rank: %lli", c.grid_rank);
+
+  char msg[200] = "\t Dimension: (";
+  for (long long i = 0; i < c.grid_rank; i++) {
+    char tmp[200] = "%lli%s";
+    if (i == c.grid_rank - 1)
+      sprintf(tmp, tmp, c.grid_dimension[i], ")");
+    else
+      sprintf(tmp, tmp, c.grid_dimension[i], ", ");
+
+    strcat(msg, tmp);
+  }
+  message("%s", msg);
+
+  if (c.heating_data)
+    message("\t Heating: (%g, ..., %g)", c.heating_data[0],
+            c.heating_data[N - 1]);
+  if (c.cooling_data)
+    message("\t Cooling: (%g, ..., %g)", c.cooling_data[0],
+            c.cooling_data[N - 1]);
+  if (c.mmw_data && print_mmw)
+    message("\t Mean molecular weigth: (%g, ..., %g)", c.mmw_data[0],
+            c.mmw_data[N - 1]);
+}
diff --git a/src/cooling/grackle/grackle_wrapper.h b/src/cooling/grackle/grackle_wrapper.h
new file mode 100644
index 0000000000000000000000000000000000000000..10c8b793ccd7785e899a97b7602cd3995e42bb20
--- /dev/null
+++ b/src/cooling/grackle/grackle_wrapper.h
@@ -0,0 +1,48 @@
+/***********************************************************************
+/
+/ Grackle c wrapper
+/
+/
+/ Copyright (c) 2013, Enzo/Grackle Development Team.
+/
+/ Distributed under the terms of the Enzo Public Licence.
+/
+/ The full license is in the file LICENSE, distributed with this
+/ software.
+************************************************************************/
+#ifndef SWIFT_COOLING_GRACKLE_WRAPPER_H
+#define SWIFT_COOLING_GRACKLE_WRAPPER_H
+
+#include "config.h"
+#include "error.h"
+
+#include <chemistry_data.h>
+#include <grackle.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "config.h"
+#include "error.h"
+
+int wrap_init_cooling(char *CloudyTable, int UVbackground, double udensity,
+                      double ulength, double utime, int grackle_chemistry);
+
+int wrap_init_cooling(char *CloudyTable, int UVbackground, double udensity,
+                      double ulength, double utime, int grackle_chemistry);
+
+int wrap_set_UVbackground_on();
+
+int wrap_set_UVbackground_off();
+
+int wrap_get_cooling_time(double rho, double u, double Z, double a_now,
+                          double *coolingtime);
+
+int wrap_do_cooling(double density, double *energy, double dtime, double Z,
+                    double a_now);
+
+void grackle_print_data();
+
+void cloudy_print_data(const cloudy_data c, const int print_mmw);
+
+#endif /* SWIFT_COOLING_GRACKLE_WRAPPER_H */
diff --git a/src/engine.c b/src/engine.c
index 497e453ec14767df9a422f810e6ec447243c525c..f21c7cfc8af721b9f8d895cffecb0c3605d34861 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -5311,7 +5311,7 @@ void engine_init(struct engine *e, struct space *s,
   if (e->policy & engine_policy_hydro)
     if (e->nodeID == 0) hydro_props_print(e->hydro_properties);
 
-  /* Print information about the hydro scheme */
+  /* Print information about the gravity scheme */
   if (e->policy & engine_policy_self_gravity)
     if (e->nodeID == 0) gravity_props_print(e->gravity_properties);
 
diff --git a/src/units.c b/src/units.c
index fef137971e58c51027da357c11007ccc0946e4b7..c9038924fa540f7df86f44db885cc2ec28672a4e 100644
--- a/src/units.c
+++ b/src/units.c
@@ -56,7 +56,7 @@ void units_init_cgs(struct unit_system* us) {
 
 /**
  * @brief Initialises the unit_system structure with the constants given in
- * rhe parameter file.
+ * the parameter file.
  *
  * @param us The unit_system to initialize.
  * @param params The parsed parameter file.
@@ -81,7 +81,7 @@ void units_init(struct unit_system* us, const struct swift_params* params,
 
 /**
  * @brief Initialises the unit_system structure with the constants given in
- * rhe parameter file. Uses a default if the values are not present in the file.
+ * the parameter file. Uses a default if the values are not present in the file.
  *
  * @param us The unit_system to initialize.
  * @param params The parsed parameter file.
@@ -594,6 +594,7 @@ double units_conversion_factor(const struct unit_system* from,
  * @param us The #unit_system
  */
 void units_print(const struct unit_system* us) {
+
   message("Units:");
   message("\tUnit Mass:        %g", us->UnitMass_in_cgs);
   message("\tUnit Length:      %g", us->UnitLength_in_cgs);