diff --git a/boututils/analyse_equil_2.py b/boututils/analyse_equil_2.py
index 93b98fc..d315c6a 100644
--- a/boututils/analyse_equil_2.py
+++ b/boututils/analyse_equil_2.py
@@ -12,7 +12,6 @@
 from past.utils import old_div
 
 import numpy
-from bunch import Bunch
 from . import local_min_max
 from scipy.interpolate import RectBivariateSpline
 from matplotlib.pyplot import contour, gradient,  annotate, plot, draw
@@ -35,8 +34,8 @@ def analyse_equil(F, R, Z):
 
     Returns
     -------
-    bunch
-        A structure of critical points containing:
+    object
+        An object of critical points containing:
 
           n_opoint, n_xpoint   - Number of O- and X-points
           primary_opt          - Index of plasma centre O-point
@@ -169,8 +168,7 @@ def analyse_equil(F, R, Z):
     print("Number of O-points: "+str(n_opoint))
 
     if n_opoint == 0 :
-        print("No O-points! Giving up on this equilibrium")
-        return Bunch(n_opoint=0, n_xpoint=0, primary_opt=-1)
+        raise RuntimeError("No O-points! Giving up on this equilibrium")
 
 
   #;;;;;;;;;;;;;; Find plasma centre ;;;;;;;;;;;;;;;;;;;
@@ -259,8 +257,8 @@ def analyse_equil(F, R, Z):
     
 
 
-  #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-  # Put results into a structure
+    #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+    # Put results into a structure
   
     result = Bunch(n_opoint=n_opoint, n_xpoint=n_xpoint, # Number of O- and X-points
             primary_opt=primary_opt, # Which O-point is the plasma centre
diff --git a/boututils/bunch.py b/boututils/bunch.py
new file mode 100644
index 0000000..2bc1ca0
--- /dev/null
+++ b/boututils/bunch.py
@@ -0,0 +1,6 @@
+# what we need from bunch
+
+class Bunch:
+    def __init__(self, **dict):
+        for k in dict:
+            setattr(self, k, dict[k])
diff --git a/boututils/check_scaling.py b/boututils/check_scaling.py
index be22fc0..af59b0b 100644
--- a/boututils/check_scaling.py
+++ b/boututils/check_scaling.py
@@ -44,15 +44,14 @@ def check_order(error_list, expected_order, tolerance=2.e-1, spacing=None):
     if len(error_list) < 2:
         raise RuntimeError("Expected at least 2 data points to calculate error")
 
-    success=True
+    success = True
+
     for i in range(len(error_list)-1):
-        if spacing is None:
-            actual_order = log(errors[i] / errors[i+1]) / log(2)
-        else:
-            actual_order = log(errors[i] / errors[i+1]) / log(grid_spacing[i] / grid_spacing[i+1])
+        grid_spacing = 2 if spacing is None else spacing[i] / spacing[i+1]
+        actual_order = log(error_list[i] / error_list[i+1]) / log(grid_spacing)
 
         if not isclose(actual_order, expected_order, atol=tolerance, rtol=0):
-            success=False
+            success = False
     return success
 
 
diff --git a/boututils/datafile.py b/boututils/datafile.py
index e8d385b..00bb987 100644
--- a/boututils/datafile.py
+++ b/boututils/datafile.py
@@ -729,11 +729,12 @@ def read(self, name, ranges=None, asBoutArray=True):
             if var is None:
                 return None
 
-        if asBoutArray:
-            attributes = self.attributes(n)
+        attributes = self.attributes(n) if asBoutArray else {}
+
+        time_dependent = attributes.get("bout_type", "none").endswith("_t")
 
         ndims = len(var.shape)
-        if ndims == 1 and var.shape[0] == 1:
+        if ndims == 1 and var.shape[0] == 1 and not time_dependent:
             data = var
             if asBoutArray:
                 data = BoutArray(data, attributes=attributes)
diff --git a/boututils/efit_analyzer.py b/boututils/efit_analyzer.py
index 0b15f42..6577b78 100644
--- a/boututils/efit_analyzer.py
+++ b/boututils/efit_analyzer.py
@@ -4,7 +4,7 @@
 from builtins import range
 from past.utils import old_div
 import numpy as np
-from bunch import Bunch
+from boututils.bunch import Bunch
 from .radial_grid import radial_grid
 from .analyse_equil_2 import analyse_equil
 from pylab import (cm, clabel, contour, draw, legend, plot, setp, show,
diff --git a/boututils/geqdsk.py b/boututils/geqdsk.py
index 4335b03..d63caf1 100755
--- a/boututils/geqdsk.py
+++ b/boututils/geqdsk.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 
 from __future__ import print_function
 from __future__ import division
diff --git a/boututils/moment_xyzt.py b/boututils/moment_xyzt.py
index 8d7dc4c..675d1f1 100644
--- a/boututils/moment_xyzt.py
+++ b/boututils/moment_xyzt.py
@@ -3,8 +3,7 @@
 from builtins import range
 from past.utils import old_div
 import numpy as np
-from bunch import Bunch
-
+from boututils.bunch import Bunch
 
 def RMSvalue( vec1d):
 #;
@@ -16,9 +15,9 @@ def RMSvalue( vec1d):
     valrms=np.sqrt(old_div(np.sum((vec1d-valav)**2),nel))
     acvec=vec1d-valav
 
-
-    return Bunch(valrms=valrms, valav=valav, acvec=acvec)
- 
+    return Bunch(valrms=valrms,
+                 valav=valav,
+                 acvec=acvec)
 
 
 
@@ -30,8 +29,6 @@ def moment_xyzt( sig_xyzt, *args):#rms=None, dc=None, ac=None):
 #; -AC (DC subtracted out), i.e., a function of (x,y,z,t)
 #;-------------------------------------------------------------------
 
-    try:                     # return to caller
- 
         d = np.shape(sig_xyzt)
         if np.size(d) != 4 :
             print("Error: Variable must be 4D (x,y,z,t)")
@@ -73,8 +70,5 @@ def moment_xyzt( sig_xyzt, *args):#rms=None, dc=None, ac=None):
             res.ac = ac
 
         if 'RMS' not in args and 'DC' not in args and 'AC' not in args :
-            print('Wrong argument')
+            raise RuntimeError('Wrong argument')
         return res
-    except:
-        print('moment_xyz failed')
-        return
diff --git a/boututils/read_geqdsk.py b/boututils/read_geqdsk.py
index 755b3f1..f665f59 100644
--- a/boututils/read_geqdsk.py
+++ b/boututils/read_geqdsk.py
@@ -1,8 +1,8 @@
 from __future__ import print_function
 from builtins import range
 import numpy
-from bunch import Bunch
 from geqdsk import Geqdsk
+from boututils.bunch import Bunch
 
 def read_geqdsk (file):
     
@@ -67,8 +67,6 @@ def read_geqdsk (file):
         for j in range(0,nyefit):
             r[i,j] = rgrid1 + xdim*i/(nxefit-1)
             z[i,j] = (zmid-0.5*zdim) + zdim*j/(nyefit-1)
-  
-
 
     f=f.T
     
@@ -89,4 +87,4 @@ def read_geqdsk (file):
             pres=pres, # Plasma pressure in nt/m^2 on uniform flux grid
             qpsi=qpsi, # q values on uniform flux grid
             nbdry=nbdry, rbdry=rbdry, zbdry=zbdry, # Plasma boundary
-            nlim=nlim, xlim=xlim, ylim=ylim) # Wall boundary
\ No newline at end of file
+            nlim=nlim, xlim=xlim, ylim=ylim) # Wall boundary
diff --git a/boututils/run_wrapper.py b/boututils/run_wrapper.py
index c25e343..9a91d2c 100644
--- a/boututils/run_wrapper.py
+++ b/boututils/run_wrapper.py
@@ -2,21 +2,20 @@
 
 from builtins import str
 import os
+import pathlib
 import re
 import subprocess
+from subprocess import call, Popen, STDOUT, PIPE
 
-try:
-    # Python 2.4 onwards
-    from subprocess import call, Popen, STDOUT, PIPE
-    lib = "call"
-except ImportError:
-    # FIXME: drop support for python < 2.4!
-    # Use os.system (depreciated)
-    from os import popen4, system
-    lib = "system"
 
+if os.name == "nt":
+    # Default on Windows
+    DEFAULT_MPIRUN = "mpiexec.exe -n"
+else:
+    DEFAULT_MPIRUN = "mpirun -np"
 
-def getmpirun(default="mpirun -np"):
+
+def getmpirun(default=DEFAULT_MPIRUN):
   """Return environment variable named MPIRUN, if it exists else return
    a default mpirun command
 
@@ -53,27 +52,20 @@ def shell(command, pipe=False):
     """
     output = None
     status = 0
-    if lib == "system":
-        if pipe:
-            handle = popen4(command)
-            output = handle[1].read()
-        else:
-            status = system(command)
+    if pipe:
+        child = Popen(command, stderr=STDOUT, stdout=PIPE, shell=True)
+        # This returns a b'string' which is casted to string in
+        # python 2. However, as we want to use f.write() in our
+        # runtest, we cast this to utf-8 here
+        output = child.stdout.read().decode("utf-8", "ignore")
+        # Wait for the process to finish. Note that child.wait()
+        # would have deadlocked the system as stdout is PIPEd, we
+        # therefore use communicate, which in the end also waits for
+        # the process to finish
+        child.communicate()
+        status = child.returncode
     else:
-        if pipe:
-            child = Popen(command, stderr=STDOUT, stdout=PIPE, shell=True)
-            # This returns a b'string' which is casted to string in
-            # python 2. However, as we want to use f.write() in our
-            # runtest, we cast this to utf-8 here
-            output = child.stdout.read().decode("utf-8", "ignore")
-            # Wait for the process to finish. Note that child.wait()
-            # would have deadlocked the system as stdout is PIPEd, we
-            # therefore use communicate, which in the end also waits for
-            # the process to finish
-            child.communicate()
-            status = child.returncode
-        else:
-            status = call(command, shell=True)
+        status = call(command, shell=True)
 
     return status, output
 
@@ -93,6 +85,18 @@ def determineNumberOfCPUs():
         The number of CPUs
     """
 
+    # cpuset
+    # cpuset may restrict the number of *available* processors
+    try:
+        m = re.search(r'(?m)^Cpus_allowed:\s*(.*)$',
+                      open('/proc/self/status').read())
+        if m:
+            res = bin(int(m.group(1).replace(',', ''), 16)).count('1')
+            if res > 0:
+                return res
+    except IOError:
+        pass
+
     # Python 2.6+
     try:
         import multiprocessing
@@ -228,8 +232,12 @@ def launch(command, runcmd=None, nproc=None, mthread=None,
         cmd = cmd + " > "+output
 
     if mthread is not None:
-        cmd = "OMP_NUM_THREADS={j} ".format(j=mthread)+cmd
-
+        if os.name == "nt":
+            # We're on windows, so we have to do it a little different
+            cmd = 'cmd /C "set OMP_NUM_THREADS={} && {}"'.format(mthread, cmd)
+        else:
+            cmd = "OMP_NUM_THREADS={} {}".format(mthread, cmd)
+        
     if verbose == True:
          print(cmd)
 
@@ -277,3 +285,48 @@ def launch_safe(command, *args, **kwargs):
                            "Output was\n\n%s"%
                            (s,command,out))
     return s, out
+
+
+def build_and_log(test):
+    """Run make and redirect the output to a log file. Prints input
+
+    On Windows, does nothing because executable should have already
+    been built
+
+    """
+
+    if os.name == "nt":
+        return
+
+    print("Making {}".format(test))
+
+    if os.path.exists("makefile") or os.path.exists("Makefile"):
+        return shell_safe("make > make.log")
+
+    ctest_filename = "CTestTestfile.cmake"
+    if not os.path.exists(ctest_filename):
+        raise RuntimeError("Could not build: no makefile and no CMake files detected")
+
+    # We're using CMake, but we need to know the target name. If
+    # bout_add_integrated_test was used (which it should have been!),
+    # then the test name is the same as the target name
+    with open(ctest_filename, "r") as f:
+        contents = f.read()
+    match = re.search("add_test.(.*) ", contents)
+    if match is None:
+        raise RuntimeError("Using CMake, but could not determine test name")
+    test_name = match.group(1)
+
+    # Now we need to find the build directory. It'll be the first
+    # parent containing CMakeCache.txt
+    here = pathlib.Path(".").absolute()
+    for parent in here.parents:
+        if (parent / "CMakeCache.txt").exists():
+            return shell_safe(
+                "cmake --build {} --target {} > make.log".format(parent, test_name)
+            )
+
+    # We've just looked up the entire directory structure and not
+    # found the build directory, this could happen if CMakeCache was
+    # deleted, in which case we can't build anyway
+    raise RuntimeError("Using CMake, but could not find build directory")
diff --git a/boututils/surface_average.py b/boututils/surface_average.py
index d58e9a3..340d252 100644
--- a/boututils/surface_average.py
+++ b/boututils/surface_average.py
@@ -11,17 +11,16 @@
 from boututils.calculus import deriv
 from boututils.int_func import int_func
 from .idl_tabulate import idl_tabulate
-from bunch import bunchify
 
 
-def surface_average(var, g, area=None):
+def surface_average(var, grid, area=None):
     """Average a variable over a surface
 
     Parameters
     ----------
     var : array_like
         3-D or 4D variable to integrate (either [x,y,z] or [t,x,y,z])
-    g : dict
+    grid : dict
         A dictionary of various grid quantities
     area : bool
         Average by flux-surface area = (B/Bp)*dl * R*dz
@@ -35,8 +34,6 @@ def surface_average(var, g, area=None):
 
     s = np.ndim(var)
 
-
-
     if s == 4 :
         nx = np.shape(var)[1]
         ny = np.shape(var)[2]
@@ -44,28 +41,21 @@ def surface_average(var, g, area=None):
 
         result = np.zeros((nx,nt))
         for t in range (nt):
-
-            result[:,t] = surface_average(var[t,:,:,:], g, area=area)
+            result[:,t] = surface_average(var[t,:,:,:], grid, area=area)
 
         return result
     elif s != 3 :
-        print("ERROR: surface_average var must be 3 or 4D")
-        return 0
+        raise RuntimeError("ERROR: surface_average var must be 3D or 4D")
 
-
-  # 3D [x,y,z]
+    # 3D [x,y,z]
     nx = np.shape(var)[0]
     ny = np.shape(var)[1]
-#    nz = np.shape(var)[2]
-
-# Use bunch to create grid structure
-    grid=bunchify(g)
 
 
-  # Calculate poloidal angle from grid
+    # Calculate poloidal angle from grid
     theta = np.zeros((nx,ny))
 
-  #status = gen_surface(mesh=grid) ; Start generator
+    #status = gen_surface(mesh=grid) ; Start generator
     xi = -1
     yi = np.arange(0,ny,dtype=int)
     last = 0
@@ -76,17 +66,17 @@ def surface_average(var, g, area=None):
             last = 1
 
         dtheta = 2.*np.pi / np.float(ny)
-        r = grid.Rxy[xi,yi]
-        z = grid.Zxy[xi,yi]
+        r = grid['Rxy'][xi,yi]
+        z = grid['Zxy'][xi,yi]
         n = np.size(r)
 
         dl = old_div(np.sqrt( deriv(r)**2 + deriv(z)**2 ), dtheta)
         if area:
-            dA = (old_div(grid.Bxy[xi,yi],grid.Bpxy[xi,yi]))*r*dl
+            dA = (old_div(grid['Bxy'][xi,yi],grid['Bpxy'][xi,yi]))*r*dl
             A = int_func(np.arange(n),dA)
             theta[xi,yi] = 2.*np.pi*A/A[n-1]
         else:
-            nu = dl * (grid.Btxy[xi,yi]) / ((grid.Bpxy[xi,yi]) * r )
+            nu = dl * (grid['Btxy'][xi,yi]) / ((grid['Bpxy'][xi,yi]) * r )
             theta[xi,yi] = int_func(np.arange(n)*dtheta,nu)
             theta[xi,yi] = 2.*np.pi*theta[xi,yi] / theta[xi,yi[n-1]]
 
diff --git a/boututils/volume_integral.py b/boututils/volume_integral.py
index 9fdb083..7e5d4da 100644
--- a/boututils/volume_integral.py
+++ b/boututils/volume_integral.py
@@ -8,17 +8,15 @@
 from past.utils import old_div
 import numpy as np
 from boututils.calculus import deriv
-from bunch import bunchify
 
-
-def volume_integral(var, g, xr=False):
+def volume_integral(var, grid, xr=False):
     """Integrate a variable over a volume
 
     Parameters
     ----------
     var : array_like
         Variable to integrate
-    g : dict
+    grid : dict
         A dictionary of various grid quantities
     xr : (int, int), optional
         Range of x indices (default: all of x)
@@ -32,9 +30,6 @@ def volume_integral(var, g, xr=False):
 
     s = np.ndim(var)
 
-    grid=bunchify(g)
-
-
     if s == 4 :
         # 4D [t,x,y,z] - integrate for each t
         nx = np.shape(var)[1]
@@ -84,20 +79,20 @@ def volume_integral(var, g, xr=False):
 
         if (xi >= np.min(xr)) & (xi <= np.max(xr)) :
             dtheta = 2.*np.pi / np.float(ny)
-            r = grid.Rxy[xi,yi]
-            z = grid.Zxy[xi,yi]
+            r = grid['Rxy'][xi,yi]
+            z = grid['Zxy'][xi,yi]
             n = np.size(r)
             dl = old_div(np.sqrt( deriv(r)**2 + deriv(z)**2 ), dtheta)
 
       # Area of flux-surface
-            dA = (grid.Bxy[xi,yi]/grid.Bpxy[xi,yi]*dl) * (r*2.*np.pi)
+            dA = (grid['Bxy'][xi,yi]/grid['Bpxy'][xi,yi]*dl) * (r*2.*np.pi)
       # Volume
             if xi == nx-1 :
-                dpsi = (grid.psixy[xi,yi] - grid.psixy[xi-1,yi])
+                dpsi = (grid['psixy'][xi,yi] - grid['psixy'][xi-1,yi])
             else:
-                dpsi = (grid.psixy[xi+1,yi] - grid.psixy[xi,yi])
+                dpsi = (grid['psixy'][xi+1,yi] - grid['psixy'][xi,yi])
 
-            dV = dA * dpsi / (r*(grid.Bpxy[xi,yi])) # May need factor of 2pi
+            dV = dA * dpsi / (r*(grid['Bpxy'][xi,yi])) # May need factor of 2pi
             dV = np.abs(dV)
 
             result = result + np.sum(var[xi,yi] * dV)
diff --git a/requirements.txt b/requirements.txt
index 1dbb35b..7f6d297 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -3,7 +3,6 @@ matplotlib==3.2.1
 h5py==2.10.0
 scipy==1.4.1
 future==0.18.2
-bunch==1.0.1
 PyQt5==5.12.2 [mayavi]
 mayavi==4.6.2 [mayavi]
 netCDF4==1.5.3
diff --git a/setup.py b/setup.py
index 260d573..5de0663 100644
--- a/setup.py
+++ b/setup.py
@@ -50,7 +50,6 @@
                       'h5py',
                       'boututils',
                       'future',
-                      'bunch',
                       'netCDF4'],
     extras_require={
                     'mayavi': ['mayavi', 'PyQt5']},