diff --git a/configure.ac b/configure.ac
index 72c6fdafe0e2b6ca89050ae00bfbdc53e14f0bad..bc63e5143b691f355efd0b7cf074da8d40432e20 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1308,7 +1308,6 @@ if test "x$with_grackle" != "xno"; then
 
    have_grackle="yes"
 
-   echo $GRACKLE_LIBS
    AS_VAR_APPEND([GRACKLE_LIBS], ["$FCLIBS"])
 
    AC_CHECK_LIB(
@@ -1317,8 +1316,24 @@ if test "x$with_grackle" != "xno"; then
       [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! You need to have a version >= 3.1.1)],
+      [AC_MSG_ERROR(Cannot find grackle library! Please consult the documentation for specific version required.)],
       [$GRACKLE_LIBS])
+
+   AC_CHECK_LIB(
+      [grackle],
+      [set_velocity_units],
+      [ : ], # : means do nothing. Leaving this argument empty triggers AC default behaviour, which breaks stuff down the line.
+      [AC_MSG_ERROR(Wrong grackle library version. Please consult the documentation for specifics.)],
+      [$GRACKLE_LIBS])
+
+   AC_CHECK_LIB(
+      [grackle],
+      [get_grackle_version],
+      [ : ], # : means do nothing
+      [AC_MSG_ERROR(Wrong grackle library version. Please consult the documentation for specifics.)],
+      [$GRACKLE_LIBS])
+
+
 fi
 AC_SUBST([GRACKLE_LIBS])
 AC_SUBST([GRACKLE_INCS])
diff --git a/doc/RTD/source/RadiativeTransfer/GEAR_RT.rst b/doc/RTD/source/RadiativeTransfer/GEAR_RT.rst
index 7fde36f0968b4b6e3b1515b99d3c5078c3affc83..0fd66e35c882e5d1fe601baedc26f0ac4e8137c5 100644
--- a/doc/RTD/source/RadiativeTransfer/GEAR_RT.rst
+++ b/doc/RTD/source/RadiativeTransfer/GEAR_RT.rst
@@ -30,7 +30,7 @@ Compiling for GEAR RT
     you to select a hydro Riemann solver, e.g ``--with-riemann-solver=hllc``.
 
 -   The thermochemistry requires the `grackle <https://github.com/grackle-project/grackle>`_ 
-    library version 3.2. Grackle is a chemistry and cooling library presented in 
+    library version above 3.2.1. [#f4]_ Grackle is a chemistry and cooling library presented in 
     `B. Smith et al. 2017 <https://ui.adsabs.harvard.edu/abs/2017MNRAS.466.2217S>`_.
     Please note that the current implementation is not (yet) as
     advanced as the :ref:`GEAR subgrid model grackle cooling <gear_grackle_cooling>`, 
@@ -38,6 +38,17 @@ Compiling for GEAR RT
     grackle cooling in combination with GEAR RT. You can however follow the Grackle 
     installation instructions documented there.
 
+.. warning::
+    (State 2023) Grackle is experiencing current development, and the API is subject
+    to changes in the future. For convenience, a frozen version is hosted as a fork
+    on github here: https://github.com/mladenivkovic/grackle-swift .
+    The version available there will be tried and tested and ensured to work with
+    GEAR-RT. 
+
+    Additionally, that repository hosts files necessary to install that specific 
+    version of grackle with spack.
+
+
 
 
 Compulsory Runtime Parameters
@@ -468,3 +479,10 @@ useful:
    `I. Iliev et al. 2006 <https://ui.adsabs.harvard.edu/abs/2006MNRAS.369.1625I>`_ 
    paper, but that is very specialized and shouldn't have much use in real 
    applications.
+
+.. [#f4] Grackle version 3.2.1 still contained a bug related to the use of their
+   "threadsafe functions" that could lead to disastrous outcomes. That was fixed
+   in commit `a59489f`. So versions after 3.2.1 should work as expected.
+   To be safe, we recommend you use the forked grackle repository specifically
+   intended to freeze a stable version for use with SWIFT. You can find that fork
+   on github: https://github.com/mladenivkovic/grackle-swift .
diff --git a/examples/RadiativeTransferTests/CoolingTest/plotSolution.py b/examples/RadiativeTransferTests/CoolingTest/plotSolution.py
index b8622ff805835589115dac07804ea31095709c7e..5ea851990a6d0284839b723ccb4467da1e52a8db 100755
--- a/examples/RadiativeTransferTests/CoolingTest/plotSolution.py
+++ b/examples/RadiativeTransferTests/CoolingTest/plotSolution.py
@@ -48,7 +48,7 @@ mass_units = unyt.Msun
 
 def mean_molecular_weight(XH0, XHp, XHe0, XHep, XHepp):
     """
-    Determines the mean molecular weight for given 
+    Determines the mean molecular weight for given
     mass fractions of
         hydrogen:   XH0
         H+:         XHp
@@ -75,7 +75,7 @@ def mean_molecular_weight(XH0, XHp, XHe0, XHep, XHepp):
 
 def gas_temperature(u, mu, gamma):
     """
-    Compute the gas temperature given the specific internal 
+    Compute the gas temperature given the specific internal
     energy u and the mean molecular weight mu
     """
 
@@ -89,7 +89,7 @@ def gas_temperature(u, mu, gamma):
 
 def get_snapshot_list(snapshot_basename="output"):
     """
-    Find the snapshot(s) that are to be plotted 
+    Find the snapshot(s) that are to be plotted
     and return their names as list
     """
 
@@ -175,7 +175,7 @@ def get_snapshot_data(snaplist):
     Returns:
         numpy arrays of:
             time
-            temperatures 
+            temperatures
             mean molecular weights
             mass fractions
     """
@@ -189,7 +189,7 @@ def get_snapshot_data(snaplist):
         # allow to read in solutions with only cooling, without RT
         with_rt = False
 
-    times = np.zeros(nsnaps) * mass_units
+    times = np.zeros(nsnaps) * unyt.Myr
     temperatures = np.zeros(nsnaps) * unyt.K
     mean_molecular_weights = np.zeros(nsnaps)
     mass_fractions = np.zeros((nsnaps, 5))
diff --git a/src/rt/GEAR/rt.h b/src/rt/GEAR/rt.h
index 9d0b21ce119ec36317028f78a9d358dc9cc0207a..72a9e25b221c162063dbffbad0e8d84abe6ed37e 100644
--- a/src/rt/GEAR/rt.h
+++ b/src/rt/GEAR/rt.h
@@ -568,7 +568,7 @@ __attribute__((always_inline)) INLINE static void rt_kick_extra(
     /* Update the mass fraction changes due to interparticle fluxes */
     const float current_mass = p->conserved.mass;
 
-    if (current_mass <= 0.f || p->rho <= 0.f) {
+    if ((current_mass <= 0.f) || (p->rho <= 0.f)) {
       /* Deal with vacuum. Let hydro deal with actuall mass < 0, just do your
        * mass fractions thing. */
       p->rt_data.tchem.mass_fraction_HI = 0.f;
@@ -683,7 +683,8 @@ __attribute__((always_inline)) INLINE static void rt_clean(
    * segfaults since we didn't malloc the stuff */
   if (!restart) {
     /* Clean up grackle data. This is a call to a grackle function */
-    _free_chemistry_data(&props->grackle_chemistry_data, &grackle_rates);
+    local_free_chemistry_data(&props->grackle_chemistry_data,
+                              &props->grackle_chemistry_rates);
 
     for (int g = 0; g < RT_NGROUPS; g++) {
       free(props->energy_weighted_cross_sections[g]);
diff --git a/src/rt/GEAR/rt_grackle_utils.h b/src/rt/GEAR/rt_grackle_utils.h
index 815e36ca8777bacd54d83c9b71721d3bba08633f..599221737365832fc1f8df665fe67757006a93e3 100644
--- a/src/rt/GEAR/rt_grackle_utils.h
+++ b/src/rt/GEAR/rt_grackle_utils.h
@@ -19,6 +19,9 @@
 #ifndef SWIFT_RT_GRACKLE_UTILS_H
 #define SWIFT_RT_GRACKLE_UTILS_H
 
+/* skip deprecation warnings. I cleaned old API calls. */
+#define OMIT_LEGACY_INTERNAL_GRACKLE_FUNC
+
 /* need hydro gamma */
 #include "hydro.h"
 
@@ -45,6 +48,7 @@
  **/
 __attribute__((always_inline)) INLINE static void rt_init_grackle(
     code_units *grackle_units, chemistry_data *grackle_chemistry_data,
+    chemistry_data_storage *grackle_chemistry_rates,
     float hydrogen_mass_fraction, const int grackle_verb,
     const int case_B_recombination, const struct unit_system *us) {
 
@@ -66,23 +70,24 @@ __attribute__((always_inline)) INLINE static void rt_init_grackle(
 
   /* Chemistry Parameters */
   /* -------------------- */
-  /* More details on https://grackle.readthedocs.io/en/latest/Parameters.html */
+  /* More details on
+   * https://grackle.readthedocs.io/en/grackle-3.2.0/Integration.html#chemistry-data
+   */
 
-  /* TODO: cleanup later with all other grackle stuff */
-  /* rtp->grackle_chemistry_data = _set_default_chemistry_parameters(); */
-  if (set_default_chemistry_parameters(grackle_chemistry_data) == 0) {
+  if (local_initialize_chemistry_parameters(grackle_chemistry_data) == 0) {
     error("Error in set_default_chemistry_parameters.");
   }
+
   /* chemistry on */
   grackle_chemistry_data->use_grackle = 1;
   /* cooling on */
   /* NOTE: without cooling on, it also won't heat... */
   grackle_chemistry_data->with_radiative_cooling = 1;
-  /* 6 species atomic H and He-> */
+  /* 6 species atomic H and He */
   grackle_chemistry_data->primordial_chemistry = 1;
-  /* dust processes */
+  /* No dust processes */
   grackle_chemistry_data->dust_chemistry = 0;
-  /* H2 formation on dust */
+  /* No H2 formation on dust */
   grackle_chemistry_data->h2_on_dust = 0;
   /* metal cooling (uses Cloudy) off (for now) */
   grackle_chemistry_data->metal_cooling = 0;
@@ -90,27 +95,21 @@ __attribute__((always_inline)) INLINE static void rt_init_grackle(
   grackle_chemistry_data->cmb_temperature_floor = 1;
   /* UV background off */
   grackle_chemistry_data->UVbackground = 0;
-  /* data file - currently not used-> */
+  /* data file - currently not used */
   grackle_chemistry_data->grackle_data_file = "";
   /* adiabatic index */
   grackle_chemistry_data->Gamma = hydro_gamma;
   /* we'll provide grackle with ionization and heating rates from RT */
   grackle_chemistry_data->use_radiative_transfer = 1;
+
   /* fraction by mass of Hydrogen in the metal-free portion of the gas */
   grackle_chemistry_data->HydrogenFractionByMass = hydrogen_mass_fraction;
   /* Use case B recombination? (On-the-spot approximation) */
   grackle_chemistry_data->CaseBRecombination = case_B_recombination;
 
-  /* TODO: cleanup later with all other grackle stuff */
-  /* Initialize the chemistry_data_storage object to be
-   * able to use local functions */
-  /* chemistry_data_storage* chem_data_storage = */
-  /*     malloc(sizeof(chemistry_data_storage)); */
-  /* rtp->grackle_chemistry_rates = chem_data_storage; */
-  /* if (!_initialize_chemistry_data(&rtp->grackle_chemistry_data, */
-  /*                                 rtp->grackle_chemistry_rates, */
-  /*                                 &rtp->grackle_units)) { */
-  if (initialize_chemistry_data(grackle_units) == 0) {
+  if (local_initialize_chemistry_data(grackle_chemistry_data,
+                                      grackle_chemistry_rates,
+                                      grackle_units) == 0) {
     error("Error in initialize_chemistry_data");
   }
 }
@@ -177,6 +176,8 @@ rt_get_grackle_particle_fields(grackle_field_data *grackle_fields,
   grackle_fields->HDI_density = NULL;
   /* for metal_cooling = 1 */
   grackle_fields->metal_density = NULL;
+  /* for use_dust_density_field = 1 */
+  grackle_fields->dust_density = NULL;
 
   /* volumetric heating rate (provide in units [erg s^-1 cm^-3]) */
   grackle_fields->volumetric_heating_rate = NULL;
@@ -196,6 +197,10 @@ rt_get_grackle_particle_fields(grackle_field_data *grackle_fields,
    * (provide in units [erg s^-1 cm^-3] / nHI in cgs) */
   grackle_fields->RT_heating_rate = malloc(FIELD_SIZE * sizeof(gr_float));
 
+  grackle_fields->H2_self_shielding_length = NULL;
+  grackle_fields->H2_custom_shielding_factor = NULL;
+  grackle_fields->isrf_habing = NULL;
+
   for (int i = 0; i < FIELD_SIZE; i++) {
 
     grackle_fields->density[i] = density;
@@ -249,24 +254,37 @@ __attribute__((always_inline)) INLINE static void rt_clean_grackle_fields(
   free(grackle_fields->grid_start);
   free(grackle_fields->grid_end);
 
+  /* initial quantities */
   free(grackle_fields->density);
   free(grackle_fields->internal_energy);
   /* free(grackle_fields->x_velocity); */
   /* free(grackle_fields->y_velocity); */
   /* free(grackle_fields->z_velocity); */
+
+  /* for primordial_chemistry >= 1 */
   free(grackle_fields->HI_density);
   free(grackle_fields->HII_density);
   free(grackle_fields->HeI_density);
   free(grackle_fields->HeII_density);
   free(grackle_fields->HeIII_density);
   free(grackle_fields->e_density);
+
+  /* for primordial_chemistry >= 2 */
   /* free(grackle_fields->HM_density); */
   /* free(grackle_fields->H2I_density); */
   /* free(grackle_fields->H2II_density); */
+
+  /* for primordial_chemistry >= 3 */
   /* free(grackle_fields->DI_density); */
   /* free(grackle_fields->DII_density); */
   /* free(grackle_fields->HDI_density); */
+
+  /* for metal_cooling = 1 */
   /* free(grackle_fields->metal_density); */
+
+  /* for use_dust_density_field = 1 */
+  /* free(grackle_fields->dust_density); */
+
   /* free(grackle_fields->volumetric_heating_rate); */
   /* free(grackle_fields->specific_heating_rate); */
 
@@ -275,6 +293,10 @@ __attribute__((always_inline)) INLINE static void rt_clean_grackle_fields(
   free(grackle_fields->RT_HeII_ionization_rate);
   free(grackle_fields->RT_H2_dissociation_rate);
   free(grackle_fields->RT_heating_rate);
+
+  /* free(grackle_fields->H2_self_shielding_length); */
+  /* free(grackle_fields->H2_custom_shielding_factor); */
+  /* free(grackle_fields->isrf_habing); */
 }
 
 /**
diff --git a/src/rt/GEAR/rt_properties.h b/src/rt/GEAR/rt_properties.h
index a9226cc7568a0458d8cbe8998c9c2a216c333d13..49b0c9c60d50067f7927faa22bcc6df5f549db2d 100644
--- a/src/rt/GEAR/rt_properties.h
+++ b/src/rt/GEAR/rt_properties.h
@@ -73,7 +73,7 @@ struct rt_props {
   float mass_fraction_HeI_init;
   float mass_fraction_HeII_init;
   float mass_fraction_HeIII_init;
-  float number_density_electrons_init; /* todo: do we need this? */
+  /* float number_density_electrons_init; [> todo: do we need this? <] */
 
   /* Hydrogen and Helium mass fractions of the non-metal portion of the gas */
   float hydrogen_mass_fraction;
@@ -114,17 +114,16 @@ struct rt_props {
   /*! grackle chemistry data */
   chemistry_data grackle_chemistry_data;
 
-  /* use case B recombination? */
+  /*! grackle chemistry data storage
+   * (needed for local function calls) */
+  chemistry_data_storage grackle_chemistry_rates;
+
+  /*! use case B recombination? */
   int case_B_recombination;
 
-  /* make grackle talkative? */
+  /*! make grackle talkative? */
   int grackle_verbose;
 
-  /* TODO: cleanup later with all other grackle stuff */
-  /*! grackle chemistry data storage
-   * (needed for local function calls) */
-  /* chemistry_data_storage* grackle_chemistry_rates; */
-
 #ifdef SWIFT_RT_DEBUG_CHECKS
   /* radiation emitted by stars this step. This is not really a property,
    * but a placeholder to sum up a global variable. It's being reset
@@ -452,8 +451,8 @@ __attribute__((always_inline)) INLINE static void rt_props_init(
   rtp->case_B_recombination = parser_get_opt_param_int(
       params, "GEARRT:case_B_recombination", /*default=*/1);
   rt_init_grackle(&rtp->grackle_units, &rtp->grackle_chemistry_data,
-                  rtp->hydrogen_mass_fraction, rtp->grackle_verbose,
-                  rtp->case_B_recombination, us);
+                  &rtp->grackle_chemistry_rates, rtp->hydrogen_mass_fraction,
+                  rtp->grackle_verbose, rtp->case_B_recombination, us);
 
   /* Pre-compute interaction rates/cross sections */
   /* -------------------------------------------- */
@@ -501,6 +500,7 @@ __attribute__((always_inline)) INLINE static void rt_struct_restore(
                       "RT properties struct");
   /* Set up stuff that needs array allocation */
   rt_init_grackle(&props->grackle_units, &props->grackle_chemistry_data,
+                  &props->grackle_chemistry_rates,
                   props->hydrogen_mass_fraction, props->grackle_verbose,
                   props->case_B_recombination, us);
 
diff --git a/src/rt/GEAR/rt_thermochemistry.h b/src/rt/GEAR/rt_thermochemistry.h
index 917016a48adca7a29e254549daa03697721b6e4e..e59ee3abda1dc33f59c232005546b793958248a9 100644
--- a/src/rt/GEAR/rt_thermochemistry.h
+++ b/src/rt/GEAR/rt_thermochemistry.h
@@ -151,9 +151,9 @@ INLINE static void rt_do_thermochemistry(
   /* Note: `grackle_rates` is a global variable defined by grackle itself.
    * Using a manually allocd and initialized variable here fails with MPI
    * for some reason. */
-  if (local_solve_chemistry(&rt_props->grackle_chemistry_data, &grackle_rates,
-                            &rt_props->grackle_units, &particle_grackle_data,
-                            dt) == 0)
+  if (local_solve_chemistry(
+          &rt_props->grackle_chemistry_data, &rt_props->grackle_chemistry_rates,
+          &rt_props->grackle_units, &particle_grackle_data, dt) == 0)
     error("Error in solve_chemistry.");
 
   /* copy updated grackle data to particle */
@@ -282,9 +282,9 @@ __attribute__((always_inline)) INLINE static float rt_tchem_get_tchem_time(
    * Using a manually allocd and initialized variable here fails with MPI
    * for some reason. */
   gr_float tchem_time;
-  if (local_calculate_cooling_time(&rt_props->grackle_chemistry_data,
-                                   &grackle_rates, &rt_props->grackle_units,
-                                   &particle_grackle_data, &tchem_time) == 0)
+  if (local_calculate_cooling_time(
+          &rt_props->grackle_chemistry_data, &rt_props->grackle_chemistry_rates,
+          &rt_props->grackle_units, &particle_grackle_data, &tchem_time) == 0)
     error("Error in calculate_cooling_time.");
 
   /* Clean up after yourself. */
diff --git a/src/rt/GEAR/rt_unphysical.h b/src/rt/GEAR/rt_unphysical.h
index 83c31063d06bf95c0b009e53078b088db4b5f23d..79b60b160b5f3490dda9c5ea92753ee1bfa1fb08 100644
--- a/src/rt/GEAR/rt_unphysical.h
+++ b/src/rt/GEAR/rt_unphysical.h
@@ -183,15 +183,23 @@ rt_check_unphysical_mass_fractions(struct part* restrict p) {
    * particles, while they themselves remain inactive. The density of such
    * inactive particles however remains zero until the particle is active
    * again. See issue #833. */
-  if (p->conserved.mass * p->rho <= 0.f) {
-    /* Deal with unphysical situations and vacuum. */
-    p->rt_data.tchem.mass_fraction_HI = 0.f;
-    p->rt_data.tchem.mass_fraction_HII = 0.f;
-    p->rt_data.tchem.mass_fraction_HeI = 0.f;
-    p->rt_data.tchem.mass_fraction_HeII = 0.f;
-    p->rt_data.tchem.mass_fraction_HeIII = 0.f;
-    return;
-  }
+
+  /* Actually, turns out that the better solution is to exception-handle the
+   * rho <= 0 case outside of this unphysical check, and not call this check
+   * at all in that case. So assume we have mass > 0, rho > 0 at this point.
+   * Problems arise during start-up, where rho = 0, and stuff gets reset where
+   * it shouldn't. */
+
+  /* TODO: verify you don't need this any longer. */
+  /* if (check_rho && (p->conserved.mass * p->rho <= 0.f)) { */
+  /*   [> Deal with unphysical situations and vacuum. <] */
+  /*   p->rt_data.tchem.mass_fraction_HI = 0.f; */
+  /*   p->rt_data.tchem.mass_fraction_HII = 0.f; */
+  /*   p->rt_data.tchem.mass_fraction_HeI = 0.f; */
+  /*   p->rt_data.tchem.mass_fraction_HeII = 0.f; */
+  /*   p->rt_data.tchem.mass_fraction_HeIII = 0.f; */
+  /*   return; */
+  /* } */
 
 #ifdef SWIFT_RT_DEBUG_CHECKS
   if (p->rt_data.tchem.mass_fraction_HI < -1e4)