diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt
index 4a5772f0..a471ea19 100644
--- a/src/clib/CMakeLists.txt
+++ b/src/clib/CMakeLists.txt
@@ -129,6 +129,7 @@ add_library(Grackle_Grackle
   set_default_chemistry_parameters.c
   solve_chemistry.c
   update_UVbackground_rates.c
+  unit_handling.h
   utils.c
 
   # auto-generated C source files
diff --git a/src/clib/calculate_cooling_time.c b/src/clib/calculate_cooling_time.c
index 7c76d75c..f897c124 100644
--- a/src/clib/calculate_cooling_time.c
+++ b/src/clib/calculate_cooling_time.c
@@ -14,10 +14,10 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
+#include "grackle.h"
 #include "grackle_macros.h"
-#include "grackle_types.h"
-#include "grackle_chemistry_data.h"
 #include "phys_constants.h"
+#include "unit_handling.h"
 #include "utils.h"
 
 extern chemistry_data *grackle_data;
@@ -25,8 +25,6 @@ extern chemistry_data_storage grackle_rates;
 
 /* function prototypes */
 
-double get_temperature_units(code_units *my_units);
-
 int update_UVbackground_rates(chemistry_data *my_chemistry,
                               chemistry_data_storage *my_rates,
                               photo_rate_storage *my_uvb_rates,
@@ -90,6 +88,17 @@ int local_calculate_cooling_time(chemistry_data *my_chemistry,
   if (!my_chemistry->use_grackle)
     return SUCCESS;
 
+  /* do unit-handling */
+  code_units units = determine_code_units(my_units, my_rates,
+                                          my_fields->current_a_value,
+                                          my_chemistry->unit_handling,
+                                          "local_calculate_cooling_time");
+  if (units.a_units < 0) {
+    return FAIL;
+  } else {
+    my_units = &units;
+  }
+
   /* Update UV background rates. */
   photo_rate_storage my_uvb_rates;
 
diff --git a/src/clib/calculate_dust_temperature.c b/src/clib/calculate_dust_temperature.c
index 025da5a6..096fe73b 100644
--- a/src/clib/calculate_dust_temperature.c
+++ b/src/clib/calculate_dust_temperature.c
@@ -14,10 +14,10 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
+#include "grackle.h"
 #include "grackle_macros.h"
-#include "grackle_types.h"
-#include "grackle_chemistry_data.h"
 #include "phys_constants.h"
+#include "unit_handling.h"
 #ifdef _OPENMP
 #include <omp.h>
 #endif
@@ -27,8 +27,6 @@ extern chemistry_data_storage grackle_rates;
 
 /* function prototypes */
 
-double get_temperature_units(code_units *my_units);
-
 int local_calculate_temperature(chemistry_data *my_chemistry,
                                 chemistry_data_storage *my_rates,
                                 code_units *my_units,
@@ -65,6 +63,17 @@ int local_calculate_dust_temperature(chemistry_data *my_chemistry,
   if (my_chemistry->dust_chemistry < 1 && my_chemistry->h2_on_dust < 1)
     return SUCCESS;
 
+  /* do unit-handling */
+  code_units units = determine_code_units(my_units, my_rates,
+                                          my_fields->current_a_value,
+                                          my_chemistry->unit_handling,
+                                          "local_calculate_dust_temperature");
+  if (units.a_units < 0) {
+    return FAIL;
+  } else {
+    my_units = &units;
+  }
+
   double co_length_units, co_density_units;
   if (my_units->comoving_coordinates == TRUE) {
     co_length_units = my_units->length_units;
diff --git a/src/clib/calculate_gamma.c b/src/clib/calculate_gamma.c
index 1b177904..5bda5bd0 100644
--- a/src/clib/calculate_gamma.c
+++ b/src/clib/calculate_gamma.c
@@ -14,11 +14,11 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
+#include "grackle.h"
 #include "grackle_macros.h"
-#include "grackle_types.h"
-#include "grackle_chemistry_data.h"
 #include "phys_constants.h"
 #include "index_helper.h"
+#include "unit_handling.h"
 #ifdef _OPENMP
 #include <omp.h>
 #endif
@@ -43,6 +43,18 @@ int local_calculate_gamma(chemistry_data *my_chemistry,
 
   if (!my_chemistry->use_grackle)
     return SUCCESS;
+
+  /* do unit-handling */
+  code_units units = determine_code_units(my_units, my_rates,
+                                          my_fields->current_a_value,
+                                          my_chemistry->unit_handling,
+                                          "local_calculate_gamma");
+  if (units.a_units < 0) {
+    return FAIL;
+  } else {
+    my_units = &units;
+  }
+
  
   const grackle_index_helper ind_helper = _build_index_helper(my_fields);
   int outer_ind, index;
diff --git a/src/clib/calculate_pressure.c b/src/clib/calculate_pressure.c
index 40f14274..9a2ddd60 100644
--- a/src/clib/calculate_pressure.c
+++ b/src/clib/calculate_pressure.c
@@ -14,11 +14,11 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
+#include "grackle.h"
 #include "grackle_macros.h"
-#include "grackle_types.h"
-#include "grackle_chemistry_data.h"
 #include "phys_constants.h"
 #include "index_helper.h"
+#include "unit_handling.h"
 #ifdef _OPENMP
 #include <omp.h>
 #endif
@@ -26,8 +26,6 @@
 extern chemistry_data *grackle_data;
 extern chemistry_data_storage grackle_rates;
 
-double get_temperature_units(code_units *my_units);
-
 int local_calculate_pressure(chemistry_data *my_chemistry,
                              chemistry_data_storage *my_rates,
                              code_units *my_units,
@@ -38,6 +36,18 @@ int local_calculate_pressure(chemistry_data *my_chemistry,
   if (!my_chemistry->use_grackle)
     return SUCCESS;
 
+  /* do unit-handling */
+  code_units units = determine_code_units(my_units, my_rates,
+                                          my_fields->current_a_value,
+                                          my_chemistry->unit_handling,
+                                          "calculate_pressure");
+  if (units.a_units < 0) {
+    return FAIL;
+  } else {
+    my_units = &units;
+  }
+
+
   double tiny_number = 1.e-20;
   const grackle_index_helper ind_helper = _build_index_helper(my_fields);
   int outer_ind, index;
diff --git a/src/clib/calculate_temperature.c b/src/clib/calculate_temperature.c
index 4466db60..61c4abb4 100644
--- a/src/clib/calculate_temperature.c
+++ b/src/clib/calculate_temperature.c
@@ -14,11 +14,11 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
+#include "grackle.h"
 #include "grackle_macros.h"
-#include "grackle_types.h"
-#include "grackle_chemistry_data.h"
 #include "phys_constants.h"
 #include "index_helper.h"
+#include "unit_handling.h"
 #ifdef _OPENMP
 #include <omp.h>
 #endif
@@ -47,8 +47,6 @@ extern void FORTRAN_NAME(calc_temp_cloudy_g)(
         double *priPar1, double *priPar2, double *priPar3, 
  	long long *priDataSize, double *priMMW);
 
-double get_temperature_units(code_units *my_units);
-
 int local_calculate_pressure(chemistry_data *my_chemistry,
                              chemistry_data_storage *my_rates,
                              code_units *my_units,
@@ -71,6 +69,17 @@ int local_calculate_temperature(chemistry_data *my_chemistry,
   if (!my_chemistry->use_grackle)
     return SUCCESS;
 
+  /* do unit-handling */
+  code_units units = determine_code_units(my_units, my_rates,
+                                          my_fields->current_a_value,
+                                          my_chemistry->unit_handling,
+                                          "local_solve_chemistry");
+  if (units.a_units < 0) {
+    return FAIL;
+  } else {
+    my_units = &units;
+  }
+
   /* Compute the pressure first. */
  
   if (my_chemistry->primordial_chemistry > 0) {
diff --git a/src/clib/grackle_chemistry_data_fields.def b/src/clib/grackle_chemistry_data_fields.def
index bcc8aa4a..9124b3df 100644
--- a/src/clib/grackle_chemistry_data_fields.def
+++ b/src/clib/grackle_chemistry_data_fields.def
@@ -208,6 +208,10 @@ ENTRY(collisional_ionisation_rates, INT, 1) //Collisional ionisation
 ENTRY(recombination_cooling_rates, INT, 1) //Recombination cooling
 ENTRY(bremsstrahlung_cooling_rates, INT, 1) //Bremsstrahlung cooling
 
+/* Flag to specify unit-handling. Use the GR_UNIT_HANDLING_LEGACY and
+   GR_UNIT_HANDLING_AUTOMATIC to set the flag values */
+ENTRY(unit_handling, INT, -3)
+
 /* maximum number of subcycle iterations for solve_chemistry */
 ENTRY(max_iterations, INT, 10000)
 
diff --git a/src/clib/initialize_chemistry_data.c b/src/clib/initialize_chemistry_data.c
index ed51108b..f045f8d7 100644
--- a/src/clib/initialize_chemistry_data.c
+++ b/src/clib/initialize_chemistry_data.c
@@ -236,6 +236,21 @@ int local_initialize_chemistry_data(chemistry_data *my_chemistry,
 //omp_set_schedule( omp_sched_auto,    chunk_size );
 # endif
 
+  if (my_chemistry->unit_handling == -3) {
+    if (grackle_verbose) {
+      fprintf(stderr, ("WARNING: unit_handling is unset. Defaulting to legacy "
+                       "handling.\n"));
+    }
+    my_chemistry->unit_handling = GR_UNIT_HANDLING_LEGACY;
+  } else if ((my_chemistry->unit_handling != GR_UNIT_HANDLING_LEGACY) &&
+             (my_chemistry->unit_handling != GR_UNIT_HANDLING_AUTOMATIC)) {
+    fprintf(stderr, "unit_handling has an invalid value\n");
+    return GR_FAIL;
+  }
+
+  /* store a copy of the initial units */
+  my_rates->initial_units = *my_units;
+
   /* Only allow a units to be one with proper coordinates. */
   if (my_units->comoving_coordinates == FALSE &&
       my_units->a_units != 1.0) {
@@ -342,8 +357,6 @@ int local_initialize_chemistry_data(chemistry_data *my_chemistry,
     return GR_FAIL;
   }
 
-  /* store a copy of the initial units */
-  my_rates->initial_units = *my_units;
 
   if (grackle_verbose) {
     time_t timer;
diff --git a/src/clib/set_default_chemistry_parameters.c b/src/clib/set_default_chemistry_parameters.c
index 964338c0..d9940394 100644
--- a/src/clib/set_default_chemistry_parameters.c
+++ b/src/clib/set_default_chemistry_parameters.c
@@ -16,8 +16,7 @@
 #include <math.h>
 
 #include "grackle_macros.h"
-#include "grackle_types.h"
-#include "grackle_chemistry_data.h"
+#include "grackle.h"
 
 int grackle_verbose = FALSE;
 
@@ -53,6 +52,7 @@ int gr_initialize_field_data(grackle_field_data *my_fields)
   my_fields->grid_dimension = NULL;
   my_fields->grid_start = NULL;
   my_fields->grid_end = NULL;
+  my_fields->current_a_value = GR_SPECIFY_INITIAL_A_VALUE;
   my_fields->grid_dx = -1.0;
 
   // now, modify all members holding datafields to have values of NULL
diff --git a/src/clib/solve_chemistry.c b/src/clib/solve_chemistry.c
index d4398af0..23d0dbe5 100644
--- a/src/clib/solve_chemistry.c
+++ b/src/clib/solve_chemistry.c
@@ -13,10 +13,10 @@
 
 #include <stdio.h>
 #include <math.h>
+#include "grackle.h"
 #include "grackle_macros.h"
-#include "grackle_types.h"
-#include "grackle_chemistry_data.h"
 #include "phys_constants.h"
+#include "unit_handling.h"
 #include "utils.h"
 
 extern chemistry_data *grackle_data;
@@ -24,8 +24,6 @@ extern chemistry_data_storage grackle_rates;
 
 /* function prototypes */
 
-double get_temperature_units(code_units *my_units);
-
 int update_UVbackground_rates(chemistry_data *my_chemistry,
                               chemistry_data_storage *my_rates,
                               photo_rate_storage *my_uvb_rates,
@@ -102,6 +100,17 @@ int local_solve_chemistry(chemistry_data *my_chemistry,
   if (!my_chemistry->use_grackle)
     return SUCCESS;
 
+  /* do unit-handling */
+  code_units units = determine_code_units(my_units, my_rates,
+                                          my_fields->current_a_value,
+                                          my_chemistry->unit_handling,
+                                          "local_solve_chemistry");
+  if (units.a_units < 0) {
+    return FAIL;
+  } else {
+    my_units = &units;
+  }
+
   /* Update UV background rates. */
   photo_rate_storage my_uvb_rates;
 
diff --git a/src/clib/unit_handling.h b/src/clib/unit_handling.h
new file mode 100644
index 00000000..05fa532d
--- /dev/null
+++ b/src/clib/unit_handling.h
@@ -0,0 +1,113 @@
+/***********************************************************************
+/
+/ Implement and declare automatic unit-handling
+/
+/ 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 UNIT_HANDLING_H
+#define UNIT_HANDLING_H
+
+#include <stdio.h>  // provides fprintf
+#include "grackle.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif /* __cplusplus */
+
+
+
+/// Automatically compute a code_units instance when Grackle is configured
+/// with automatic code-units.
+///
+/// @param specified_units Pointer to the specified units
+/// @param my_rates Pointer to the chemistry_data_storage object that contains
+///     the initial rates that were specified when the grackle_solver was
+///     initialized
+/// @param current_a_value The specified current scale-factor
+/// @param unit_handling The value of Grackle's unit handling flag
+/// @param func_name Name of the function that is calling the error check
+///
+/// @returns a fully initialized code_units struct. If there was any kind of
+///     issue, then the a_units member will be negative.
+///
+/// @note
+/// This is declared inline to reduce as much overhead as possible
+///
+/// @note
+/// This is currently a proof-of-concept, some improvements are definitely
+/// possible. For example, we could:
+///   * directly pass in my_rates->initial_units rather than my_rates
+///   * refactor gr_query_units so that we could directly call into its
+///     internals (the error-checking and string comparisons introduce a bunch
+///     of unnecessary overhead)
+static inline code_units determine_code_units (
+    const code_units* specified_units,
+    const chemistry_data_storage* my_rates,
+    double current_a_value,
+    int unit_handling,
+    const char* func_name)
+{
+  code_units out;
+  out.a_units = -1.0;
+
+  if (unit_handling == GR_UNIT_HANDLING_LEGACY) {
+    // handle error-cases
+    if (current_a_value != GR_SPECIFY_INITIAL_A_VALUE) {
+      fprintf(stderr,
+          "ERROR in %s: current_a_value shouldn't be specified when using "
+          "Grackle's legacy unit handling\n", func_name);
+      out.a_units = -1.0;
+      return out;
+    } else if (specified_units == NULL) {
+      fprintf(stderr,
+          "ERROR in %s: code_units argument can't be NULL when using "
+          "Grackle's legacy unit handling\n", func_name);
+      return out;
+    }
+
+    // this is the happy path
+    out = *specified_units;
+    return out;
+
+  } else if (unit_handling == GR_UNIT_HANDLING_AUTOMATIC) {
+    // handle error-case
+    if (specified_units != NULL) {
+      fprintf(stderr,
+          "ERROR in %s: code_units argument must be NULL when using Grackle's "
+          "automatic unit handling\n", func_name);
+      return out;
+    }
+
+    // this is the happy path
+    out = my_rates->initial_units;
+    out.a_value = current_a_value;
+    out.density_units = gr_query_units(my_rates, "density_units",
+                                       current_a_value);
+    out.length_units = gr_query_units(my_rates, "length_units",
+                                      current_a_value);
+    out.time_units = gr_query_units(my_rates, "time_units",
+                                    current_a_value);
+    set_velocity_units(&out);
+    return out;
+
+  } else {
+    fprintf(stderr,
+        "ERROR in %s: Grackle's unit_handling parameter has an unrecognized "
+        "value\n", func_name);
+    return out;
+  }
+
+}
+
+#ifdef __cplusplus
+} /* extern "C" */
+#endif /* __cplusplus */
+
+#endif /* UNIT_HANDLING_H */
+
diff --git a/src/example/cxx_example.C b/src/example/cxx_example.C
index c91f4f20..b7589e8f 100644
--- a/src/example/cxx_example.C
+++ b/src/example/cxx_example.C
@@ -77,6 +77,7 @@ int main(int argc, char *argv[])
 
   // Create struct for storing grackle field data
   grackle_field_data my_fields;
+  gr_initialize_field_data(&my_fields);
 
   // Set grid dimension and size.
   // grid_start and grid_end are used to ignore ghost zones.
diff --git a/src/include/grackle.h b/src/include/grackle.h
index 3796e787..14e1e3d5 100644
--- a/src/include/grackle.h
+++ b/src/include/grackle.h
@@ -26,6 +26,9 @@ extern "C" {
 #define GR_FAIL 0
 
 #define GR_SPECIFY_INITIAL_A_VALUE -1
+#define GR_UNIT_HANDLING_LEGACY -2
+#define GR_UNIT_HANDLING_AUTOMATIC -1
+
 
 extern int grackle_verbose;
 
diff --git a/src/include/grackle_chemistry_data.h b/src/include/grackle_chemistry_data.h
index 22ce5007..9a2cc223 100644
--- a/src/include/grackle_chemistry_data.h
+++ b/src/include/grackle_chemistry_data.h
@@ -176,6 +176,9 @@ typedef struct
   int recombination_cooling_rates; //Recombination cooling
   int bremsstrahlung_cooling_rates; //Bremsstrahlung cooling
 
+  /* Flag to specify unit-handling */
+  int unit_handling;
+
   /* maximum number of subcycle iterations for solve_chemistry */
   int max_iterations;
 
diff --git a/src/include/grackle_fortran_interface.def b/src/include/grackle_fortran_interface.def
index 6629c69e..6f166d24 100644
--- a/src/include/grackle_fortran_interface.def
+++ b/src/include/grackle_fortran_interface.def
@@ -40,7 +40,9 @@ c     This is the fortran definition of grackle_field_data
          TYPE(C_PTR) :: grid_dimension
          TYPE(C_PTR) :: grid_start
          TYPE(C_PTR) :: grid_end
+         REAL(C_DOUBLE) :: current_a_value
          REAL(C_DOUBLE) :: grid_dx
+
          
          TYPE(C_PTR) :: density
          TYPE(C_PTR) :: HI_density
@@ -252,4 +254,4 @@ c     The following define the fortran interfaces to the C routines
             IMPORT
             TYPE(grackle_field_data), INTENT(INOUT) :: my_fields
          END FUNCTION gr_initialize_field_data
-      END INTERFACE
\ No newline at end of file
+      END INTERFACE
diff --git a/src/include/grackle_types.h b/src/include/grackle_types.h
index 3a5f2d86..f24b71a1 100644
--- a/src/include/grackle_types.h
+++ b/src/include/grackle_types.h
@@ -50,6 +50,7 @@ typedef struct
   int *grid_start;
   int *grid_end;
 
+  double current_a_value;
   gr_float grid_dx;
 
   gr_float *density;