diff --git a/postprocessor.py b/postprocessor.py
index 5a3a049..ef1f2fa 100644
--- a/postprocessor.py
+++ b/postprocessor.py
@@ -944,26 +944,34 @@ class PostProcessor(HDF5Container):
         dmap_omega = rt_omega.process(self._cam[ax_los]).map.T
         return dmap_omega
 
-    def _toomreQ_disk(self, ax_los, omega_approx=True):
+    def _toomreQ_disk(self, ax_los, omega_approx=True, G1_units=True):
         """
         Compute the Toomre Q parameter
         """
 
         # Get maps
 
-        cs2 = self.get_value("/maps/T_mwavg_z")
-        coldens = self.get_value("/maps/coldens_z")
-        omega = self.get_value("/maps/omega_mwavg_z")
-
+        if G1_units:
+            G = 1.
+            cs2 = self.get_value("/maps/T_mwavg_z")
+            coldens = self.get_value("/maps/coldens_z")
+            omega = self.get_value("/maps/omega_mwavg_z")
+            space_coeff = self.lbox
+        else:
+            G = U.G.express(U.kg**-1 * U.m ** 3 * U.s ** -2)
+            cs2 = self.get_value("/maps/T_mwavg_z", unit=U.m ** 2 * U.s ** -2)
+            coldens = self.get_value("/maps/coldens_z", unit=U.kg * U.m ** -2)
+            omega = self.get_value("/maps/omega_mwavg_z", unit=U.s ** -1)
+            space_coeff = self.info["unit_length"].express(U.m)
         # Compute Q
         if omega_approx:
             # Use angular frequency as epiciclic frequency (true if the disk is Keplerian)
             kappa = omega
         else:
             # Get coordinates
-            im_extent = np.array(self.save.root.maps._v_attrs.im_extent) * self.lbox
+            im_extent = np.array(self.save.root.maps._v_attrs.im_extent) * space_coeff
             map_size = self.pp_params.pymses.map_size
-            center = np.array(self.pp_params.disk.center) * self.lbox
+            center = np.array(self.pp_params.disk.center) * space_coeff
 
             # Physical size of cells
             dx = (im_extent[1] - im_extent[0]) / map_size
@@ -982,8 +990,7 @@ class PostProcessor(HDF5Container):
             kappa_square = (2 * omega / rr) * Gr
             kappa = np.sqrt(kappa_square)
 
-        map_Q = (np.sqrt(cs2) * kappa) / (np.pi * self.G * coldens)
-
+        map_Q = (np.sqrt(cs2) * kappa) / (np.pi * G * coldens)
         return map_Q
 
     def _radial_bins(self, ax_los="z"):