diff --git a/configure.ac b/configure.ac
index d8cd381b359335486f7e2ea517bb08b65103de47..a0aee165452d5d67a3a145c3f03785a22b74df19 100644
--- a/configure.ac
+++ b/configure.ac
@@ -478,8 +478,8 @@ if test "x$with_grackle" != "xno"; then
       
    AC_CHECK_LIB(
       [grackle],
-      [initialize_grackle],
-      [AC_DEFINE([HAVE_GRACKLE],1,[The GRACKLE library appears to be present.]),
+      [initialize_chemistry_data],
+      [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!)],
diff --git a/src/Makefile.am b/src/Makefile.am
index c50c831b3e27ab07ade0c83ccb5cb3b4f78f6724..df1ed0a670892ecd2a41b229a8707ffb993a7cc3 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -48,11 +48,6 @@ 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 \
@@ -62,8 +57,7 @@ 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 \
-    $(GRACKLE_SRC)
+    collectgroup.c hydro_space.c equation_of_state.c
 
 
 # Include files for distribution, not installation.
@@ -123,7 +117,6 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.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
 
 
diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h
index 50700071904be662609e513b34a0215d0ef61d7f..f24105d70fbd681ea29e897bb74e43a09968ac28 100644
--- a/src/cooling/grackle/cooling.h
+++ b/src/cooling/grackle/cooling.h
@@ -27,6 +27,7 @@
 /* Some standard headers. */
 #include <float.h>
 #include <math.h>
+#include <grackle.h>
 
 /* Local includes. */
 #include "error.h"
@@ -36,13 +37,67 @@
 #include "physical_constants.h"
 #include "units.h"
 
-/* include the grackle wrapper */
-#include "grackle_wrapper.h"
+/* need to rework (and check) code if changed */
+#define GRACKLE_NPART 1
+#define GRACKLE_RANK 3
 
 /**
- * @brief Compute the cooling rate
+ * @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 Prints the properties of the cooling model to stdout.
  *
- * We do nothing.
+ * @param cooling The properties of the cooling function.
+ */
+__attribute__((always_inline))INLINE static void cooling_print_backend(
+    const struct cooling_function_data* cooling) {
+
+  message("Cooling function is 'Grackle'.");
+  message("Using Grackle           = %i", cooling->chemistry.use_grackle);
+  message("Chemical network        = %i", cooling->chemistry.primordial_chemistry);
+  message("Radiative cooling       = %i", cooling->chemistry.with_radiative_cooling);
+  message("Metal cooling           = %i", cooling->chemistry.metal_cooling);
+  
+  message("CloudyTable             = %s",
+          cooling->cloudy_table);
+  message("UVbackground            = %d", cooling->uv_background);
+  message("Redshift                = %g", cooling->redshift);
+  message("Density Self Shielding  = %g",
+          cooling->density_self_shielding);
+  message("Units:");
+  message("\tComoving     = %i", cooling->units.comoving_coordinates);
+  message("\tLength       = %g", cooling->units.length_units);
+  message("\tDensity      = %g", cooling->units.density_units);
+  message("\tTime         = %g", cooling->units.time_units);
+  message("\tScale Factor = %g", cooling->units.a_units);
+
+}
+
+
+/**
+ * @brief Compute the cooling rate
  *
  * @param phys_const The physical constants in internal units.
  * @param us The internal system of units.
@@ -58,33 +113,103 @@ __attribute__((always_inline)) INLINE static double cooling_rate(
     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;
+  if (cooling->chemistry.primordial_chemistry > 1)
+    error("Not implemented");
+
+  /* set current time */
+  code_units units = cooling->units;
+  if (cooling->redshift == -1)
+    error("TODO time dependant redshift");
+  else
+    units.a_value = 1. / (1. + cooling->redshift);
+
+  /* initialize data */
+  grackle_field_data data;
+
+  /* set values */
+  /* grid */
+  int grid_dimension[GRACKLE_RANK] = {GRACKLE_NPART, 1, 1};
+  int grid_start[GRACKLE_RANK] = {0, 0, 0};
+  int grid_end[GRACKLE_RANK] = {GRACKLE_NPART - 1, 0, 0};
+  
+  data.grid_rank = GRACKLE_RANK;
+  data.grid_dimension = grid_dimension;
+  data.grid_start = grid_start;
+  data.grid_end = grid_end;
+
+  /* general particle data */
+  gr_float density = hydro_get_density(p);
+  const double energy_before = hydro_get_internal_energy(p);
+  gr_float energy = energy_before;
+  gr_float vx = 0;
+  gr_float vy = 0;
+  gr_float vz = 0;
+
+  data.density = &density;
+  data.internal_energy = &energy;
+  data.x_velocity = &vx;
+  data.y_velocity = &vy;
+  data.z_velocity = &vz;
+
+  /* /\* primordial chemistry >= 1 *\/ */
+  /* gr_float HI_density = density; */
+  /* gr_float HII_density = 0.; */
+  /* gr_float HeI_density = 0.; */
+  /* gr_float HeII_density = 0.; */
+  /* gr_float HeIII_density = 0.; */
+  /* gr_float e_density = 0.; */
+  
+  /* data.HI_density = &HI_density; */
+  /* data.HII_density = &HII_density; */
+  /* data.HeI_density = &HeI_density; */
+  /* data.HeII_density = &HeII_density; */
+  /* data.HeIII_density = &HeIII_density; */
+  /* data.e_density = &e_density; */
+
+  /* /\* primordial chemistry >= 2 *\/ */
+  /* gr_float HM_density = 0.; */
+  /* gr_float H2I_density = 0.; */
+  /* gr_float H2II_density = 0.; */
+  
+  /* data.HM_density = &HM_density; */
+  /* data.H2I_density = &H2I_density; */
+  /* data.H2II_density = &H2II_density; */
+  
+  /* /\* primordial chemistry >= 3 *\/ */
+  /* gr_float DI_density = 0.; */
+  /* gr_float DII_density = 0.; */
+  /* gr_float HDI_density = 0.; */
+  
+  /* data.DI_density = &DI_density; */
+  /* data.DII_density = &DII_density; */
+  /* data.HDI_density = &HDI_density; */
+
+  /* metal cooling = 1 */
+  gr_float metal_density = density * grackle_data->SolarMetalFractionByMass;
+
+  data.metal_density = &metal_density;
+
+  /* /\* volumetric heating rate *\/ */
+  /* gr_float volumetric_heating_rate = 0.; */
+
+  /* data.volumetric_heating_rate = &volumetric_heating_rate; */
+  
+  /* /\* specific heating rate *\/ */
+  /* gr_float specific_heating_rate = 0.; */
+
+  /* data.specific_heating_rate = &specific_heating_rate; */
+
+  /* solve chemistry with table */
+  if (solve_chemistry(&units, &data, dt) == 0) {
+    error("Error in solve_chemistry.");
   }
-
-  return (u_new - u_old) / dt;
+  
+  return (energy - energy_before) / 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.
@@ -102,15 +227,11 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
   /* 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;
+  /* compute cooling rate */
+  const float du_dt = cooling_rate(phys_const, us, cooling, p, dt);
 
   /* record energy lost */
-  xp->cooling_data.radiated_energy += -delta_u * hydro_get_mass(p);
+  xp->cooling_data.radiated_energy += - du_dt * dt * hydro_get_mass(p);
 
   /* Update the internal energy */
   hydro_set_internal_energy_dt(p, hydro_du_dt + du_dt);
@@ -134,30 +255,6 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
   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.
  *
@@ -166,33 +263,76 @@ __attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
  * @param phys_const The physical constants in internal units.
  * @param cooling The cooling properties to initialize
  */
-static INLINE void cooling_init_backend(
+__attribute__((always_inline))INLINE static 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;
-
+    /* read parameters */
   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(
+                          cooling->cloudy_table);
+  cooling->uv_background =
+    parser_get_param_int(parameter_file, "GrackleCooling:UVbackground");
+
+  cooling->redshift =
+    parser_get_param_double(parameter_file, "GrackleCooling:GrackleRedshift");
+
+  cooling->density_self_shielding = parser_get_param_double(
       parameter_file, "GrackleCooling:GrackleHSShieldingDensityThreshold");
+  
+#ifdef SWIFT_DEBUG_CHECKS
+  /* enable verbose for grackle */
+  grackle_verbose = 1;
+#endif
 
-  UVbackground = cooling->UVbackground;
-  grackle_chemistry = 0; /* forced to be zero : read table */
+  /* Set up the units system.
+     These are conversions from code units to cgs. */
+
+  /* first cosmo */
+  cooling->units.a_units = 1.0;  // units for the expansion factor (1/1+zi)
+  cooling->units.a_value = 1.0;
+
+  /* We assume here all physical quantities to
+     be in proper coordinate (not comobile)  */
+  cooling->units.comoving_coordinates = 0;
+
+  /* then units */
+  cooling->units.density_units = us->UnitMass_in_cgs / pow(us->UnitLength_in_cgs, 3);
+  cooling->units.length_units = us->UnitLength_in_cgs;
+  cooling->units.time_units = us->UnitTime_in_cgs;
+  cooling->units.velocity_units =
+    cooling->units.a_units * cooling->units.length_units / cooling->units.time_units;
+
+  chemistry_data *chemistry = &cooling->chemistry;
+  
+  /* Create a chemistry object for parameters and rate data. */
+  if (set_default_chemistry_parameters(chemistry) == 0) {
+    error("Error in set_default_chemistry_parameters.");
+  }
+
+  // Set parameter values for chemistry.
+  chemistry->use_grackle = 1;
+  chemistry->with_radiative_cooling = 1;
+  /* molecular network with H, He, D
+   From Cloudy table */
+  chemistry->primordial_chemistry = 0;
+  chemistry->metal_cooling = 1;  // metal cooling on
+  chemistry->UVbackground = cooling->uv_background;
+  chemistry->grackle_data_file = cooling->cloudy_table;
+  chemistry->use_radiative_transfer = 0;
+  chemistry->use_volumetric_heating_rate = 0;
+  chemistry->use_specific_heating_rate = 0;
+
+  /* Initialize the chemistry object. */
+  if (initialize_chemistry_data(&cooling->units) == 0) {
+    error("Error in initialize_chemistry_data.");
+  }
 
-  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;
+  if (GRACKLE_NPART != 1)
+    error("Grackle with multiple particles not implemented");
+  float threshold = cooling->density_self_shielding;
 
   threshold /= phys_const->const_proton_mass;
   threshold /= pow(us->UnitLength_in_cgs, 3);
@@ -200,41 +340,14 @@ static INLINE void cooling_init_backend(
   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
+  cooling_print_backend(cooling);
+  message("Density Self Shielding = %g atom/cm3", threshold);
 
-  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
index b59af375cba0ed35fcb0159cb63605af8d94f030..ddc07ed70e36e175b3535e1a48be3ff8349d363c 100644
--- a/src/cooling/grackle/cooling_struct.h
+++ b/src/cooling/grackle/cooling_struct.h
@@ -19,6 +19,9 @@
 #ifndef SWIFT_COOLING_STRUCT_NONE_H
 #define SWIFT_COOLING_STRUCT_NONE_H
 
+/* include grackle */
+#include <grackle.h>
+
 /**
  * @file src/cooling/none/cooling_struct.h
  * @brief Empty infrastructure for the cases without cooling function
@@ -30,16 +33,22 @@
 struct cooling_function_data {
 
   /* Filename of the Cloudy Table */
-  char GrackleCloudyTable[200];
+  char cloudy_table[200];
 
   /* Enable/Disable UV backgroud */
-  int UVbackground;
+  int uv_background;
 
   /* Redshift to use for the UV backgroud (-1 to use cosmological one) */
-  double GrackleRedshift;
+  double redshift;
 
   /* Density Threshold for the shielding */
-  double GrackleHSShieldingDensityThreshold;
+  double density_self_shielding;
+
+  /* unit system */
+  code_units units;
+
+  /* grackle chemistry data */
+  chemistry_data chemistry;
 };
 
 /**
diff --git a/src/cooling/grackle/grackle_wrapper.c b/src/cooling/grackle/grackle_wrapper.c
deleted file mode 100644
index 9974717b670c937811cea54349124d9a38bed476..0000000000000000000000000000000000000000
--- a/src/cooling/grackle/grackle_wrapper.c
+++ /dev/null
@@ -1,227 +0,0 @@
-/***********************************************************************
-/
-/ 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
deleted file mode 100644
index 10c8b793ccd7785e899a97b7602cd3995e42bb20..0000000000000000000000000000000000000000
--- a/src/cooling/grackle/grackle_wrapper.h
+++ /dev/null
@@ -1,48 +0,0 @@
-/***********************************************************************
-/
-/ 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 */