diff --git a/gprMax/cmds_geometry/fractal_box_builder.py b/gprMax/cmds_geometry/fractal_box_builder.py
index ce72cc78..f120ca05 100644
--- a/gprMax/cmds_geometry/fractal_box_builder.py
+++ b/gprMax/cmds_geometry/fractal_box_builder.py
@@ -73,7 +73,8 @@ class FractalBoxBuilder(UserObjectGeometry):
 
                 # If there is only 1 bin then a normal material is being used, otherwise a mixing model
                 if volume.nbins == 1:
-                    volume.fractalvolume = np.ones((volume.nx, volume.ny, volume.nz), dtype=config.sim_config.dtype['float_or_double'])
+                    volume.fractalvolume = np.ones((volume.nx, volume.ny, volume.nz),
+                                                    dtype=config.sim_config.dtypes['float_or_double'])
                     materialnumID = next(x.numID for x in grid.materials if x.ID == volume.operatingonID)
                     volume.fractalvolume *= materialnumID
                 else:
@@ -115,7 +116,8 @@ class FractalBoxBuilder(UserObjectGeometry):
                                     if surface.fractalsurface[j - surface.ys, k - surface.zs] > 0:
                                         height = 0
                                         for i in range(volume.xs, surface.fractalrange[1]):
-                                            if i < surface.fractalsurface[j - surface.ys, k - surface.zs] and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] != 1:
+                                            if (i < surface.fractalsurface[j - surface.ys, k - surface.zs] and
+                                                volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] != 1):
                                                 y, z = g.calculate_blade_geometry(blade, height)
                                                 # Add y, z coordinates to existing location
                                                 yy = int(j - volume.ys + y)
@@ -136,7 +138,8 @@ class FractalBoxBuilder(UserObjectGeometry):
                                         depth = 0
                                         i = volume.xf - 1
                                         while i > volume.xs:
-                                            if i > volume.originalxf - (surface.fractalsurface[j - surface.ys, k - surface.zs] - volume.originalxf) and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] == 1:
+                                            if (i > volume.originalxf - (surface.fractalsurface[j - surface.ys, k - surface.zs] -
+                                                volume.originalxf) and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] == 1):
                                                 y, z = g.calculate_root_geometry(root, depth)
                                                 # Add y, z coordinates to existing location
                                                 yy = int(j - volume.ys + y)
@@ -181,7 +184,8 @@ class FractalBoxBuilder(UserObjectGeometry):
                                     if surface.fractalsurface[i - surface.xs, k - surface.zs] > 0:
                                         height = 0
                                         for j in range(volume.ys, surface.fractalrange[1]):
-                                            if j < surface.fractalsurface[i - surface.xs, k - surface.zs] and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] != 1:
+                                            if (j < surface.fractalsurface[i - surface.xs, k - surface.zs] and
+                                                volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] != 1):
                                                 x, z = g.calculate_blade_geometry(blade, height)
                                                 # Add x, z coordinates to existing location
                                                 xx = int(i - volume.xs + x)
@@ -202,7 +206,8 @@ class FractalBoxBuilder(UserObjectGeometry):
                                         depth = 0
                                         j = volume.yf - 1
                                         while j > volume.ys:
-                                            if j > volume.originalyf - (surface.fractalsurface[i - surface.xs, k - surface.zs] - volume.originalyf) and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] == 1:
+                                            if (j > volume.originalyf - (surface.fractalsurface[i - surface.xs, k - surface.zs] -
+                                                volume.originalyf) and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] == 1):
                                                 x, z = g.calculate_root_geometry(root, depth)
                                                 # Add x, z coordinates to existing location
                                                 xx = int(i - volume.xs + x)
@@ -247,7 +252,8 @@ class FractalBoxBuilder(UserObjectGeometry):
                                     if surface.fractalsurface[i - surface.xs, j - surface.ys] > 0:
                                         height = 0
                                         for k in range(volume.zs, surface.fractalrange[1]):
-                                            if k < surface.fractalsurface[i - surface.xs, j - surface.ys] and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] != 1:
+                                            if (k < surface.fractalsurface[i - surface.xs, j - surface.ys] and
+                                                volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] != 1):
                                                 x, y = g.calculate_blade_geometry(blade, height)
                                                 # Add x, y coordinates to existing location
                                                 xx = int(i - volume.xs + x)
@@ -268,7 +274,8 @@ class FractalBoxBuilder(UserObjectGeometry):
                                         depth = 0
                                         k = volume.zf - 1
                                         while k > volume.zs:
-                                            if k > volume.originalzf - (surface.fractalsurface[i - surface.xs, j - surface.ys] - volume.originalzf) and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] == 1:
+                                            if (k > volume.originalzf - (surface.fractalsurface[i - surface.xs, j - surface.ys] -
+                                                volume.originalzf) and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] == 1):
                                                 x, y = g.calculate_root_geometry(root, depth)
                                                 # Add x, y coordinates to existing location
                                                 xx = int(i - volume.xs + x)
@@ -287,7 +294,10 @@ class FractalBoxBuilder(UserObjectGeometry):
                 grassnumID = next((x.numID for x in grid.materials if x.ID == 'grass'), 0)
                 data = volume.fractalvolume.astype('int16', order='C')
                 mask = volume.mask.copy(order='C')
-                build_voxels_from_array_mask(volume.xs, volume.ys, volume.zs, waternumID, grassnumID, volume.averaging, mask, data, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
+                build_voxels_from_array_mask(volume.xs, volume.ys, volume.zs,
+                                             waternumID, grassnumID, volume.averaging,
+                                             mask, data, grid.solid, grid.rigidE,
+                                             grid.rigidH, grid.ID)
 
             else:
                 if volume.nbins == 1:
@@ -297,4 +307,6 @@ class FractalBoxBuilder(UserObjectGeometry):
                     volume.fractalvolume += volume.mixingmodel.startmaterialnum
 
                 data = volume.fractalvolume.astype('int16', order='C')
-                build_voxels_from_array(volume.xs, volume.ys, volume.zs, 0, volume.averaging, data, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
+                build_voxels_from_array(volume.xs, volume.ys, volume.zs, 0,
+                                        volume.averaging, data, grid.solid,
+                                        grid.rigidE, grid.rigidH, grid.ID)
diff --git a/gprMax/cmds_multiple.py b/gprMax/cmds_multiple.py
index 90fff1e3..5663c147 100644
--- a/gprMax/cmds_multiple.py
+++ b/gprMax/cmds_multiple.py
@@ -26,6 +26,7 @@ from .cmds_geometry.cmds_geometry import UserObjectGeometry
 from .exceptions import CmdInputError
 from .geometry_outputs import GeometryObjects as GeometryObjectsUser
 from .materials import Material as MaterialUser
+from .materials import DispersiveMaterial as DispersiveMaterialUser
 from .materials import PeplinskiSoil as PeplinskiSoilUser
 from .pml import CFSParameter
 from .pml import CFS
@@ -508,7 +509,7 @@ class Rx(UserObjectMulti):
 
         try:
             r.ID = self.kwargs['id']
-            outputs = self.kwargs['outputs']
+            outputs = [self.kwargs['outputs']]
             # Get allowable outputs
             allowableoutputs = RxUser.allowableoutputs_gpu if config.sim_config.general['cuda'] else RxUser.allowableoutputs
             # Check and add field output names
@@ -767,20 +768,28 @@ class AddDebyeDispersion(UserObjectMulti):
             raise CmdInputError(f"'{self.params_str()}' material(s) {notfound} do not exist")
 
         for material in materials:
-            material.type = 'debye'
-            material.poles = poles
-            material.averagable = False
+            disp_material = DispersiveMaterialUser(material.numID, material.ID)
+            disp_material.er = material.er
+            disp_material.se = material.se
+            disp_material.mr = material.mr
+            disp_material.sm = material.sm
+            disp_material.type = 'debye'
+            disp_material.poles = poles
+            disp_material.averagable = False
             for i in range(0, poles):
-                # N.B Not checking if relaxation times are greater than time-step
                 if tau[i] > 0:
-                    material.deltaer.append(er_delta[i])
-                    material.tau.append(tau[i])
+                    log.debug('Not checking if relaxation times are greater than time-step')
+                    disp_material.deltaer.append(er_delta[i])
+                    disp_material.tau.append(tau[i])
                 else:
                     raise CmdInputError(f"'{self.params_str()}' requires positive values for the permittivity difference.")
-            if material.poles > MaterialUser.maxpoles:
-                MaterialUser.maxpoles = material.poles
+            if disp_material.poles > config.model_configs[grid.model_num].materials['maxpoles']:
+                config.model_configs[grid.model_num].materials['maxpoles'] = disp_material.poles
 
-            log.info(f"Debye disperion added to {material.ID} with delta_eps_r={', '.join('%4.2f' % deltaer for deltaer in material.deltaer)}, and tau={', '.join('%4.3e' % tau for tau in material.tau)} secs created.")
+            # Replace original material with newly created DispersiveMaterial
+            grid.materials = [disp_material if mat.numID==material.numID else mat for mat in grid.materials]
+
+            log.info(f"Debye disperion added to {disp_material.ID} with delta_eps_r={', '.join('%4.2f' % deltaer for deltaer in disp_material.deltaer)}, and tau={', '.join('%4.3e' % tau for tau in disp_material.tau)} secs created.")
 
 
 class AddLorentzDispersion(UserObjectMulti):
@@ -1028,7 +1037,7 @@ class GeometryView(UserObjectMulti):
 
         g = GeometryViewUser(xs, ys, zs, xf, yf, zf, dx, dy, dz, filename, fileext, grid)
 
-        log.info(f'Geometry view from {xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, {yf * grid.dy:g}m, {zf * grid.dz:g}m, discretisation {dx * grid.dx:g}m, {dy * grid.dy:g}m, {dz * grid.dz:g}m, multi_grid={self.multi_grid}, grid={grid.name}, with filename base {g.basefilename} created.')
+        log.info(f'Geometry view from {xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, {yf * grid.dy:g}m, {zf * grid.dz:g}m, discretisation {dx * grid.dx:g}m, {dy * grid.dy:g}m, {dz * grid.dz:g}m, multi_grid={self.multi_grid}, grid={grid.name}, with filename base {g.filename} created.')
 
         # Append the new GeometryView object to the geometry views list
         grid.geometryviews.append(g)
diff --git a/gprMax/config.py b/gprMax/config.py
index eab33301..a0606e93 100644
--- a/gprMax/config.py
+++ b/gprMax/config.py
@@ -187,7 +187,7 @@ class SimulationConfig:
                         'cuda': False,
                         'opencl': False,
                         'precision': 'single',
-                        'autotranslate': False}
+                        'autotranslate': self.args.autotranslate}
         log.debug('Should autotranslate be a ModelConfig parameter?')
 
         self.em_consts = {'c': c, # Speed of light in free space (m/s)
diff --git a/gprMax/geometry_outputs.py b/gprMax/geometry_outputs.py
index dc5c7561..1d4aba60 100644
--- a/gprMax/geometry_outputs.py
+++ b/gprMax/geometry_outputs.py
@@ -19,7 +19,6 @@
 import logging
 import os
 from pathlib import Path
-import sys
 
 import h5py
 import numpy as np
@@ -61,7 +60,7 @@ class GeometryView:
         self.dx = dx
         self.dy = dy
         self.dz = dz
-        self.basefilename = filename
+        self.filename = filename
         self.fileext = fileext
         self.G = G
 
@@ -89,19 +88,19 @@ class GeometryView:
             self.vtk_nylines = self.ny * (self.nx + 1) * (self.nz + 1)
             self.vtk_nzlines = self.nz * (self.nx + 1) * (self.ny + 1)
             self.vtk_numlines = self.vtk_nxlines + self.vtk_nylines + self.vtk_nzlines
-            self.vtk_connectivity_offset = round_value((self.vtk_numpoints *
-                                                       self.vtk_numpoint_components *
-                                                       np.dtype(np.float32).itemsize) +
-                                                       np.dtype(np.uint32).itemsize)
-            self.vtk_offsets_offset = round_value(self.vtk_connectivity_offset +
-                                                 (self.vtk_numlines *
-                                                 self.vtk_numline_components *
-                                                 np.dtype(np.uint32).itemsize) +
-                                                 np.dtype(np.uint32).itemsize)
-            self.vtk_materials_offset = round_value(self.vtk_offsets_offset +
-                                                   (self.vtk_numlines *
-                                                   np.dtype(np.uint32).itemsize) +
-                                                   np.dtype(np.uint32).itemsize)
+            self.vtk_connectivity_offset = ((self.vtk_numpoints *
+                                           self.vtk_numpoint_components *
+                                           np.dtype(np.float32).itemsize) +
+                                           np.dtype(np.uint32).itemsize)
+            self.vtk_offsets_offset = (self.vtk_connectivity_offset +
+                                      (self.vtk_numlines *
+                                      self.vtk_numline_components *
+                                      np.dtype(np.uint32).itemsize) +
+                                      np.dtype(np.uint32).itemsize)
+            self.vtk_materials_offset = (self.vtk_offsets_offset +
+                                        (self.vtk_numlines *
+                                        np.dtype(np.uint32).itemsize) +
+                                        np.dtype(np.uint32).itemsize)
             vtk_cell_offsets = (((self.vtk_numline_components * self.vtk_numlines) +
                                self.vtk_numline_components - self.vtk_numline_components - 1) //
                                self.vtk_numline_components + 1)
@@ -114,8 +113,8 @@ class GeometryView:
 
     def set_filename(self):
         """Construct filename from user-supplied name and model run number."""
-        parts = config.model_configs[self.G.model_num].output_file_path.with_suffix('').parts
-        self.filename = Path(*parts[:-1], parts[-1])
+        parts = config.model_configs[self.G.model_num].output_file_path.parts
+        self.filename = Path(*parts[:-1], self.filename + config.model_configs[self.G.model_num].appendmodelnumber)
         self.filename = self.filename.with_suffix(self.fileext)
 
     def write_vtk(self, G, pbar):
@@ -142,19 +141,15 @@ class GeometryView:
             for index, rx in enumerate(G.rxs):
                 self.rxs[rx.xcoord, rx.ycoord, rx.zcoord] = index + 1
 
-            vtk_srcs_pml_offset = round_value((np.dtype(np.uint32).itemsize *
-                                               self.vtk_nxcells *
-                                               self.vtk_nycells *
-                                               self.vtk_nzcells) +
-                                               np.dtype(np.uint32).itemsize)
-            vtk_rxs_offset = round_value((np.dtype(np.uint32).itemsize *
-                                          self.vtk_nxcells * self.vtk_nycells *
-                                          self.vtk_nzcells) +
-                                          np.dtype(np.uint32).itemsize +
-                                          (np.dtype(np.int8).itemsize *
-                                          self.vtk_nxcells * self.vtk_nycells *
-                                          self.vtk_nzcells) +
-                                          np.dtype(np.uint32).itemsize)
+            vtk_srcs_pml_offset = ((np.dtype(np.uint32).itemsize * self.vtk_nxcells *
+                                  self.vtk_nycells * self.vtk_nzcells) +
+                                  np.dtype(np.uint32).itemsize)
+            vtk_rxs_offset = ((np.dtype(np.uint32).itemsize * self.vtk_nxcells *
+                             self.vtk_nycells * self.vtk_nzcells) +
+                             np.dtype(np.uint32).itemsize +
+                             (np.dtype(np.int8).itemsize * self.vtk_nxcells *
+                             self.vtk_nycells * self.vtk_nzcells) +
+                             np.dtype(np.uint32).itemsize)
 
             with open(self.filename, 'wb') as f:
                 f.write('<?xml version="1.0"?>\n'.encode('utf-8'))
@@ -428,7 +423,7 @@ class GeometryViewFineMultiGrid:
         self.nx = G.nx
         self.ny = G.ny
         self.nz = G.nz
-        self.basefilename = filename
+        self.filename = filename
         self.fileext = '.vtp'
         self.sg_views = []
 
@@ -455,8 +450,8 @@ class GeometryViewFineMultiGrid:
 
     def set_filename(self):
         """Construct filename from user-supplied name and model run number."""
-        parts = config.model_configs[self.G.model_num].input_file_path.parts
-        self.filename = Path(*parts[:-1], parts[-1] + config.model_configs[self.G.model_num].appendmodelnumber)
+        parts = config.model_configs[self.G.model_num].output_file_path.parts
+        self.filename = Path(*parts[:-1], self.filename + config.model_configs[self.G.model_num].appendmodelnumber)
         self.filename = self.filename.with_suffix(self.fileext)
 
     def write_vtk(self, *args):
diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py
index c0b439ec..fe279e12 100644
--- a/gprMax/gprMax.py
+++ b/gprMax/gprMax.py
@@ -37,6 +37,7 @@ def run(
     mpicomm=None,
     gpu=None,
     subgrid=None,
+    autotranslate=False,
     geometry_only=False,
     geometry_fixed=False,
     write_processed=False,
@@ -45,10 +46,13 @@ def run(
         programming interface (API). Run the simulation for the
         given list of scenes.
 
-    :param scenes: list of the scenes to run the model. Multiple scene objects can given in order to run multiple simulation runs. Each scene must contain the essential simulation objects
+    :param scenes: list of the scenes to run the model. Multiple scene objects
+                    can given in order to run multiple simulation runs. Each
+                    scene must contain the essential simulation objects
     :type scenes: list, optional
 
-    :param inputfile:  input file path. Can also run simulation by providing an input file.
+    :param inputfile:  input file path. Can also run simulation by providing an
+                        input file.
     :type inputfile: str, optional
 
     :param outputfile: file path to the output data file.
@@ -57,32 +61,54 @@ def run(
     :param n: number of required simulation runs.
     :type n: int, non-optional
 
-    :param task: task identifier (model number) when running simulation as a job array on open grid scheduler/grid engine. for further details see the parallel performance section of the user guide
+    :param task: task identifier (model number) when running simulation as a
+                    job array on open grid scheduler/grid engine. For further
+                    details see the parallel performance section of the User Guide.
     :type task: int, optional
 
-    :param restart: model number to start/restart simulation from. It would typically be used to restart a series of models from a specific model number, with the n argument, e.g. to restart from A-scan 45 when creating a B-scan with 60 traces
+    :param restart: model number to start/restart simulation from. It would
+                    typically be used to restart a series of models from a
+                    specific model number, with the n argument, e.g. to restart
+                    from A-scan 45 when creating a B-scan with 60 traces.
     :type restart: int, optional
 
-    :param mpi: number of Message Passing Interface (MPI) tasks, i.e. master + workers, for MPI task farm. This option is most usefully combined with n to allow individual models to be farmed out using a MPI task farm, e.g. to create a B-scan with 60 traces and use MPI to farm out each trace1. For further details see the parallel performance section of the User Guide
+    :param mpi: number of Message Passing Interface (MPI) tasks,
+                i.e. master + workers, for MPI task farm. This option is most
+                usefully combined with n to allow individual models to be farmed
+                out using a MPI task farm, e.g. to create a B-scan with 60 traces
+                and use MPI to farm out each trace1. For further details see the
+                parallel performance section of the User Guide.
     :type mpi: int, optional
 
-    :param mpi_no_spawn: flag to use MPI task farm without spawn mechanism. For further details see the parallel performance section of the User Guide.
+    :param mpi_no_spawn: flag to use MPI task farm without spawn mechanism.
+                            For further details see the parallel performance
+                            section of the User Guide.
     :type mpi_no_spawn: bool, optional
 
-    :param gpu: flag to use NVIDIA GPU or list of NVIDIA GPU device ID(s) for specific GPU card(s)
+    :param gpu: flag to use NVIDIA GPU or list of NVIDIA GPU device ID(s) for
+                specific GPU card(s).
     :type gpu: list or bool, optional
 
     :param subgrid: flag to use sub-gridding.
     :type subgrid: bool, optional
 
-    :param geometry_only: build a model and produce any geometry views but do not run the simulation.
+    :param autotranslate: for sub-gridding - auto translate objects with main grid
+                            coordinates to their equivalent local grid coordinate
+                            within the subgrid. If this option is off users must
+                            specify sub-grid object point within the global
+                            subgrid space.
+    :type autotranslate: bool, optional
+
+    :param geometry_only: build a model and produce any geometry views but do
+                            not run the simulation.
     :type geometry_only: bool, optional
 
-    :param geometry_fixed: build a model and produce any geometry views but do not run the simulation.
+    :param geometry_fixed: run a series of models where the geometry does not
+                            change between models.
     :type geometry_fixed: bool, optional
 
     :param write_processed: write another input file after any Python code and
-    in the original input file has been processed.
+                            in the original input file has been processed.
     :type write_processed: bool, optional
     """
 
@@ -101,7 +127,8 @@ def run(
     args.mpi_no_spawn = mpi_no_spawn
     args.mpicomm = mpicomm
     args.gpu = gpu
-    args.subgrid=subgrid
+    args.subgrid = subgrid
+    args.autotranslate = autotranslate
     args.geometry_only = geometry_only
     args.geometry_fixed = geometry_fixed
     args.write_processed = write_processed
diff --git a/gprMax/grid.py b/gprMax/grid.py
index 53fe8e19..f2321cdc 100644
--- a/gprMax/grid.py
+++ b/gprMax/grid.py
@@ -48,7 +48,7 @@ class FDTDGrid:
         """
 
         self.title = ''
-        self.name = 'Main'
+        self.name = 'main_grid'
         self.model_num = model_num
 
         self.nx = 0
@@ -346,7 +346,7 @@ class CUDAGrid(FDTDGrid):
 
     def reset_fields(self):
         """Clear arrays for field components and PMLs."""
-        
+
         super().reset_fields()
 
         # Clear arrays for field components
diff --git a/gprMax/hash_cmds_file.py b/gprMax/hash_cmds_file.py
index 3ed0d650..745bb4de 100644
--- a/gprMax/hash_cmds_file.py
+++ b/gprMax/hash_cmds_file.py
@@ -62,7 +62,7 @@ def process_python_include_code(inputfile, usernamespace):
 
         # Process any Python code
         if(inputlines[x].startswith('#python:')):
-
+            log.warning('#python blocks are deprecated and will be removed in the next release of gprMax. Please convert your model to use our Python API instead.')
             # String to hold Python code to be executed
             pythoncode = ''
             x += 1
diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py
index b9a5a147..ba8b61e6 100644
--- a/gprMax/model_build_run.py
+++ b/gprMax/model_build_run.py
@@ -209,16 +209,17 @@ class ModelBuildRun:
             gb.build_materials()
 
         # Check to see if numerical dispersion might be a problem
-        results = dispersion_analysis(G)
-        if results['error']:
-            log.warning(Fore.RED + f"\nNumerical dispersion analysis not carried out as {results['error']}" + Style.RESET_ALL)
-        elif results['N'] < config.model_configs[G.model_num].numdispersion['mingridsampling']:
-            raise GeneralError(f"Non-physical wave propagation: Material '{results['material'].ID}' has wavelength sampled by {results['N']} cells, less than required minimum for physical wave propagation. Maximum significant frequency estimated as {results['maxfreq']:g}Hz")
-        elif (results['deltavp'] and np.abs(results['deltavp']) >
-              config.model_configs[G.model_num].numdispersion['maxnumericaldisp']):
-            log.warning(Fore.RED + f"\nPotentially significant numerical dispersion. Estimated largest physical phase-velocity error is {results['deltavp']:.2f}% in material '{results['material'].ID}' whose wavelength sampled by {results['N']} cells. Maximum significant frequency estimated as {results['maxfreq']:g}Hz" + Style.RESET_ALL)
-        elif results['deltavp']:
-            log.info(f"\nNumerical dispersion analysis: estimated largest physical phase-velocity error is {results['deltavp']:.2f}% in material '{results['material'].ID}' whose wavelength sampled by {results['N']} cells. Maximum significant frequency estimated as {results['maxfreq']:g}Hz")
+        for grid in grids:
+            results = dispersion_analysis(grid)
+            if results['error']:
+                log.warning(Fore.RED + f"\nNumerical dispersion analysis ({grid.name}) not carried out as {results['error']}" + Style.RESET_ALL)
+            elif results['N'] < config.model_configs[G.model_num].numdispersion['mingridsampling']:
+                raise GeneralError(f"\nNon-physical wave propagation in {grid.name} detected. Material '{results['material'].ID}' has wavelength sampled by {results['N']} cells, less than required minimum for physical wave propagation. Maximum significant frequency estimated as {results['maxfreq']:g}Hz")
+            elif (results['deltavp'] and np.abs(results['deltavp']) >
+                  config.model_configs[G.model_num].numdispersion['maxnumericaldisp']):
+                log.warning(Fore.RED + f"\n{grid.name} has potentially significant numerical dispersion. Estimated largest physical phase-velocity error is {results['deltavp']:.2f}% in material '{results['material'].ID}' whose wavelength sampled by {results['N']} cells. Maximum significant frequency estimated as {results['maxfreq']:g}Hz" + Style.RESET_ALL)
+            elif results['deltavp']:
+                log.info(f"\nNumerical dispersion analysis ({grid.name}): estimated largest physical phase-velocity error is {results['deltavp']:.2f}% in material '{results['material'].ID}' whose wavelength sampled by {results['N']} cells. Maximum significant frequency estimated as {results['maxfreq']:g}Hz")
 
     def reuse_geometry(self):
         # Reset iteration number
@@ -293,8 +294,6 @@ class ModelBuildRun:
             tsolve (float): time taken to execute solving (seconds).
         """
 
-        log.info(f'\nOutput file: {config.model_configs[self.G.model_num].output_file_path_ext.name}')
-
         # Check number of OpenMP threads
         if config.sim_config.general['cpu']:
             log.info(f'CPU (OpenMP) threads for solving: {config.model_configs[self.G.model_num].ompthreads}\n')
@@ -306,7 +305,7 @@ class ModelBuildRun:
 
         # Prepare iterator
         if config.sim_config.is_messages():
-            iterator = tqdm(range(self.G.iterations), desc='Running simulation, model ' + str(self.G.model_num + 1) + '/' + str(config.sim_config.model_end), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.sim_config.general['progressbars'])
+            iterator = tqdm(range(self.G.iterations), desc=f'Running model {self.G.model_num + 1}/{config.sim_config.model_end}, {config.model_configs[self.G.model_num].output_file_path_ext.name}', ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.sim_config.general['progressbars'])
         else:
             iterator = range(self.G.iterations)
 
@@ -329,7 +328,7 @@ class GridBuilder:
     def build_pmls(self):
         log.info('')
         pbar = tqdm(total=sum(1 for value in self.grid.pmlthickness.values() if value > 0),
-                    desc=f'Building {self.grid.name} Grid PML boundaries',
+                    desc=f'Building PML boundaries ({self.grid.name})',
                     ncols=get_terminal_width() - 1, file=sys.stdout,
                     disable=not config.sim_config.general['progressbars'])
         for pml_id, thickness in self.grid.pmlthickness.items():
@@ -342,7 +341,7 @@ class GridBuilder:
         # Build the model, i.e. set the material properties (ID) for every edge
         # of every Yee cell
         log.info('')
-        pbar = tqdm(total=2, desc=f'Building {self.grid.name} Grid',
+        pbar = tqdm(total=2, desc=f'Building Yee cells ({self.grid.name})',
                     ncols=get_terminal_width() - 1, file=sys.stdout,
                     disable=not config.sim_config.general['progressbars'])
         build_electric_components(self.grid.solid, self.grid.rigidE, self.grid.ID, self.grid)
@@ -373,5 +372,5 @@ class GridBuilder:
         materialstable.outer_border = False
         materialstable.justify_columns[0] = 'right'
 
-        log.info(f'\n{self.grid.name} Grid Materials:')
+        log.info(f'\nMaterials ({self.grid.name}):')
         log.info(materialstable.table)
diff --git a/gprMax/scene.py b/gprMax/scene.py
index 2b306670..f8438847 100644
--- a/gprMax/scene.py
+++ b/gprMax/scene.py
@@ -133,11 +133,6 @@ class Scene:
         # Process main grid multiple commands
         self.process_cmds(self.multiple_cmds, G)
 
-        # Estimate and check memory (RAM) usage
-        log.debug('Fix memory checks')
-        # G.memory_check()
-        #snapshot_memory_check(G)
-
         # Initialise an array for volumetric material IDs (solid), boolean
         # arrays for specifying materials not to be averaged (rigid),
         # an array for cell edge IDs (ID)
diff --git a/gprMax/subgrids/base.py b/gprMax/subgrids/base.py
index 86176d63..172af766 100644
--- a/gprMax/subgrids/base.py
+++ b/gprMax/subgrids/base.py
@@ -24,8 +24,8 @@ from ..grid import FDTDGrid
 
 class SubGridBase(FDTDGrid):
 
-    def __init__(self, **kwargs):
-        super().__init__()
+    def __init__(self, *args, **kwargs):
+        super().__init__(args[0])
 
         self.mode = '3D'
         self.ratio = kwargs['ratio']
diff --git a/gprMax/subgrids/subgrid_hsg.py b/gprMax/subgrids/subgrid_hsg.py
index 9c615e48..ea944c4f 100644
--- a/gprMax/subgrids/subgrid_hsg.py
+++ b/gprMax/subgrids/subgrid_hsg.py
@@ -30,8 +30,8 @@ class SubGridHSG(SubGridBase):
 
     gridtype = '3DSUBGRID'
 
-    def __init__(self, **kwargs):
-        super().__init__(**kwargs)
+    def __init__(self, model_num, **kwargs):
+        super().__init__(model_num, **kwargs)
         self.gridtype = SubGridHSG.gridtype
 
     def update_magnetic_is(self, precursors):
@@ -40,7 +40,7 @@ class SubGridHSG(SubGridBase):
             Args:
                 nwl, nwm, nwn, face, field, inc_field, lookup_id, sign, mod, co
         """
-        
+
         # Hz = c0Hz - c1Ey + c2Ex
         # Hy = c0Hy - c3Ex + c1Ez
         # Hx = c0Hx - c2Ez + c3Ey
@@ -131,37 +131,29 @@ class SubGridHSG(SubGridBase):
 
     def __str__(self):
 
-        self.memory_estimate_basic()
+        mem_use = self.mem_est_basic()
 
         s = '\n'
         s += Fore.CYAN
         s += 'Sub Grid HSG\n'
-        s += 'Name: {}\n'.format(self.name)
-        s += 'dx, dy, dz: {}m {}m {}m\n'.format(self.dx, self.dy, self.dz)
-        s += 'dt: {}s\n'.format(self.dt)
-        s += 'Memory Estimate: {}\n'.format(human_size(self.memoryusage))
-        s += 'Position: ({}m, {}m, {}m), ({}m, {}m, {}m)\n'.format(self.x1,
-                                                                   self.y1,
-                                                                   self.z1,
-                                                                   self.x2,
-                                                                   self.y2,
-                                                                   self.z2)
-        s += 'Main Grid Indices: lower left({}, {}, {}), upper right({}, {}, {})\n'.format(self.i0, self.j0, self.k0, self.i1, self.j1, self.k1)
-        s += 'Total Cells: {} {} {}\n'.format(self.nx, self.ny, self.nz)
-        s += 'Working Region Cells: {} {} {}\n'.format(self.nwx,
-                                                       self.nwy,
-                                                       self.nwz)
+        s += f'Name: {self.name}\n'
+        s += f'dx, dy, dz: {self.dx}m {self.dy}m {self.dz}m\n'
+        s += f'dt: {self.dt}s\n'
+        s += f'Memory Estimate: {human_size(mem_use)}\n'
+        s += f'Position: ({self.x1}m, {self.y1}m, {self.z1}m), ({self.x2}m, {self.y2}m, {self.z2}m)\n'
+        s += f'Main Grid Indices: lower left({self.io}, {self.j0}, {self.k0}), upper right({self.i1}, {self.j1}, {self.k1})\n'
+        s += f'Total Cells: {self.nx} {self.ny} {self.nz}\n'
+        s += f'Working Region Cells: {self.nwx} {self.nwy} {self.nwz}\n'
+
         for h in self.hertziandipoles:
-            s += 'Hertizian dipole: {} {} {}\n'.format(h.xcoord,
-                                                       h.ycoord,
-                                                       h.zcoord)
+            s += f'Hertizian dipole: {h.xcoord} {h.ycoord} {h.zcoord}\n'
             s += str([x for x in self.waveforms
                       if x.ID == h.waveformID][0]) + '\n'
         for r in self.rxs:
-            s += 'Receiver: {} {} {}\n'.format(r.xcoord, r.ycoord, r.zcoord)
+            s += f'Receiver: {r.xcoord} {r.ycoord} {r.zcoord}\n'
 
         for tl in self.transmissionlines:
-            s += 'Transmission Line: {} {} {}\n'.format(tl.xcoord, tl.ycoord, tl.zcoord)
+            s += f'Transmission Line: {tl.xcoord} {tl.ycoord} {tl.zcoord}\n'
             s += str([x for x in self.waveforms
                       if x.ID == tl.waveformID][0]) + '\n'
         s += Style.RESET_ALL
diff --git a/gprMax/subgrids/user_objects.py b/gprMax/subgrids/user_objects.py
index ca5b7657..7421f959 100644
--- a/gprMax/subgrids/user_objects.py
+++ b/gprMax/subgrids/user_objects.py
@@ -96,9 +96,6 @@ class SubGridBase(UserObjectMulti):
         # Set the temporal discretisation
         sg.calculate_dt()
 
-        # ensure stability
-        sg.round_time_step()
-
         # set the indices related to the subgrids main grid placement
         self.set_main_grid_indices(sg, grid, uip, p1, p2)
 
@@ -128,8 +125,6 @@ class SubGridBase(UserObjectMulti):
 
         # Copy over built in materials
         sg.materials = [copy(m) for m in grid.materials if m.numID in range(0, grid.n_built_in_materials + 1)]
-        # use same number of threads
-        sg.nthreads = grid.nthreads
 
         # Dont mix and match different subgrids
         for sg_made in grid.subgrids:
@@ -194,10 +189,8 @@ class SubGridHSG(SubGridBase):
         self.hash = '#subgrid_hsg'
 
     def create(self, grid, uip):
-        sg = SubGridHSGUser(**self.kwargs)
+        sg = SubGridHSGUser(grid.model_num, **self.kwargs)
         self.setup(sg, grid, uip)
-        if config.is_messages():
-            print(sg)
         return sg
 
 
@@ -219,4 +212,4 @@ class ReferenceRx(Rx):
             r.offset = ratio // 2
 
         except KeyError:
-            raise CmdInputError("'{}' has an no ratio parameter".format(self.__str__()))
+            raise CmdInputError(f"'{self.__str__()}' has an no ratio parameter")
diff --git a/gprMax/user_inputs.py b/gprMax/user_inputs.py
index 85b624bd..86862221 100644
--- a/gprMax/user_inputs.py
+++ b/gprMax/user_inputs.py
@@ -46,10 +46,10 @@ def create_user_input_points(grid, user_obj):
     if isinstance(grid, SubGridBase):
         # Local object configuration trumps. User can turn of autotranslate for
         # specfic objects.
-        if not user_obj.autotranslate and config.general['autotranslate']:
+        if not user_obj.autotranslate and config.sim_config.general['autotranslate']:
             return MainGridUserInput(grid)
 
-        if config.general['autotranslate']:
+        if config.sim_config.general['autotranslate']:
             return SubgridUserInput(grid)
         else:
             return MainGridUserInput(grid)
diff --git a/gprMax/utilities.py b/gprMax/utilities.py
index dbf043e6..af47e686 100644
--- a/gprMax/utilities.py
+++ b/gprMax/utilities.py
@@ -28,7 +28,7 @@ import subprocess
 from shutil import get_terminal_size
 import sys
 import textwrap
-from time import process_time
+from time import perf_counter
 
 from colorama import init
 from colorama import Fore
@@ -485,4 +485,4 @@ def detect_check_gpus(deviceIDs):
 
 def timer():
     """Function to return the current process wide time in fractional seconds."""
-    return process_time()
+    return perf_counter()
diff --git a/todo b/todo
deleted file mode 100644
index f9902507..00000000
--- a/todo
+++ /dev/null
@@ -1 +0,0 @@
-memory update for all grid (model_build_run.py)
diff --git a/tools/plot_Ascan.py b/tools/plot_Ascan.py
index 6764d019..ad53231e 100644
--- a/tools/plot_Ascan.py
+++ b/tools/plot_Ascan.py
@@ -18,7 +18,6 @@
 
 import argparse
 import os
-import sys
 
 import h5py
 import numpy as np
@@ -51,7 +50,7 @@ def mpl_plot(filename, outputs=Rx.defaultoutputs, fft=False):
 
     # Check there are any receivers
     if nrx == 0:
-        raise CmdInputError('No receivers found in {}'.format(filename))
+        raise CmdInputError(f'No receivers found in {filename}')
 
     # Check for single output component when doing a FFT
     if fft:
@@ -77,7 +76,7 @@ def mpl_plot(filename, outputs=Rx.defaultoutputs, fft=False):
                 output = outputs[0]
 
             if output not in availableoutputs:
-                raise CmdInputError('{} output requested to plot, but the available output for receiver 1 is {}'.format(output, ', '.join(availableoutputs)))
+                raise CmdInputError(f"{output} output requested to plot, but the available output for receiver 1 is {', '.join(availableoutputs)}")
 
             outputdata = f[path + output][:] * polarity
 
@@ -165,7 +164,7 @@ def mpl_plot(filename, outputs=Rx.defaultoutputs, fft=False):
 
                 # Check if requested output is in file
                 if output not in availableoutputs:
-                    raise CmdInputError('Output(s) requested to plot: {}, but available output(s) for receiver {} in the file: {}'.format(', '.join(outputs), rx, ', '.join(availableoutputs)))
+                    raise CmdInputError(f"Output(s) requested to plot: {', '.join(outputs)}, but available output(s) for receiver {rx} in the file: {', '.join(availableoutputs)}")
 
                 outputdata = f[path + output][:] * polarity
 
diff --git a/user_models/example_subgrid.py b/user_models/example_subgrid.py
deleted file mode 100644
index bfef3e7e..00000000
--- a/user_models/example_subgrid.py
+++ /dev/null
@@ -1,64 +0,0 @@
-import gprMax
-
-scene = gprMax.Scene()
-
-# global configuration
-# objects within the inner surface can be positioned automatically
-# in this configuration coordinates in the main grid Space
-# are converted to subgrid local space automatically
-# if objects or required in the region between the inner Surface
-# and the outer surface autotranslate must be set to False
-# and the local position must be set explicitly. Addionally,
-# in this instance, objects which occur within the IS OS gap
-# must be positioned in the subgrid and the main grid.
-gprMax.config.general['autotranslate'] = True
-
-# Points are specified as iterables
-p1 = (1, 1, 1)
-
-# single objects
-dxdydz = gprMax.Discretisation(p1=(1e-2, 1e-2, 1e-2))
-messages = gprMax.Messages(yn='y')
-tw = gprMax.TimeWindow(iterations=250)
-domain = gprMax.Domain(p1=p1)
-title = gprMax.Title(name='mysimulation')
-waveform = gprMax.Waveform(wave_type='gaussiandotnorm', amp=1, freq=400e6, id='mywaveform')
-
-# position in subgrid automatically
-dipole = gprMax.HertzianDipole(p1=(0.5, 0.5, 0.5), polarisation='y', waveform_id='mywaveform')
-
-# equivalent position when autotranslate is false
-#rx = gprMax.Rx(p1=(50 * 1e-2 / 5, 50 * 1e-2 / 5, 50 * 1e-2 / 5))
-
-# rx in the subgrid
-# subgrid data is stored in an additional output file called
-# outputfile_<subgrid_id>.out
-rx1 = gprMax.Rx(p1=(0.52, 0.5, 0.5))
-
-# rx in the main grid
-rx2 = gprMax.Rx(p1=(0.2, 0.5, 0.5))
-
-# set the position of the Inner Surface
-subgrid = gprMax.SubGridHSG(p1=(45e-2, 45e-2, 45e-2),
-                            p2=(55e-2, 55e-2, 55e-2),
-                            ratio=3,
-                            id='subgrid1',
-                            filter=False)
-                            #pml_separation=6,
-                            #subgrid_pml_thickness=10)
-
-scene = gprMax.Scene()
-scene.add(dxdydz)
-scene.add(messages)
-scene.add(tw)
-scene.add(domain)
-scene.add(title)
-scene.add(subgrid)
-scene.add(rx2)
-
-subgrid.add(waveform)
-subgrid.add(dipole)
-subgrid.add(rx1)
-
-#gprMax.run(scenes=[scene], n=1, geometry_only=False, outputfile='mysimulation', geometry_fixed=False, subgrid=True)
-gprMax.run(scenes=[scene], n=1, geometry_only=False, outputfile='mysimulation', geometry_fixed=False, subgrid=True)
diff --git a/user_models/sub-gridding/gssi_400_over_fractal_subsurface.py b/user_models/sub-gridding/gssi_400_over_fractal_subsurface.py
new file mode 100644
index 00000000..7ab9b6e6
--- /dev/null
+++ b/user_models/sub-gridding/gssi_400_over_fractal_subsurface.py
@@ -0,0 +1,137 @@
+# Copyright (C) 2015-2019, John Hartley
+#
+# This module is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License.
+# To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/.
+
+from pathlib import Path
+
+import gprMax
+from user_libs.antennas.GSSI import antenna_like_GSSI_400
+
+import numpy as np
+
+# file path step
+fn = Path(__file__).with_suffix('')
+parts = fn.parts
+
+# subgrid Discretisation is 1 mm in x, y, z directions. This allows us
+# to model the geometry of the antenna
+dl_s = 1e-3
+# subgridding ratio. This must always be an odd integer multiple. In this case
+# the main grid discrestisation is 9e-3 m.
+ratio = 9
+dl = dl_s * 9
+
+# estimated return time for signal to propagate 1 metre and back
+tw = 2 / 3e8 * (np.sqrt(3.2) + np.sqrt(9))
+
+# domain extent
+x = 3
+y = 1
+z = 2
+
+scene = gprMax.Scene()
+
+title_gpr = gprMax.Title(name=fn.name)
+dxdydz = gprMax.Discretisation(p1=(dl, dl, dl))
+domain = gprMax.Domain(p1=(x, y, z))
+time_window = gprMax.TimeWindow(time=tw)
+
+scene.add(domain)
+scene.add(title_gpr)
+scene.add(dxdydz)
+scene.add(time_window)
+
+# half space material
+halfspace_m = gprMax.Material(er=3.2, se=0.397e-3, mr=1, sm=0, id='soil')
+scene.add(halfspace_m)
+
+antenna_case = (0.3, 0.3, 0.178)
+
+bounding_box = 2 * dl
+
+# pml + boundary_cells + is_os + subgrid+boundary + half antenna
+x0 = (10 + 15 + 5 + 2) * dl + antenna_case[0] / 2
+
+#antenna_p = (x / 2, y / 2, z - 30 * dl - bounding_box - antenna_case[2])
+
+# Position of antenna base
+antenna_p = (x / 2, y / 2, 170 * dl)
+
+sg_x0 = antenna_p[0] - antenna_case[0] / 2 - bounding_box
+sg_y0 = antenna_p[1] - antenna_case[1] / 2 - bounding_box
+sg_z0 = antenna_p[2] - bounding_box
+
+sg_x1 = antenna_p[0] + antenna_case[0] / 2 + bounding_box
+sg_y1 = antenna_p[1] + antenna_case[1] / 2 + bounding_box
+sg_z1 = antenna_p[2] + antenna_case[2] + bounding_box
+
+sg = gprMax.SubGridHSG(p1=[sg_x0, sg_y0, sg_z0], p2=[sg_x1, sg_y1, sg_z1], ratio=ratio, id='subgrid1')
+scene.add(sg)
+
+# half space box in main grid
+halfspace = gprMax.Box(p1=(0, 0, 0), p2=(x, y, antenna_p[2]), material_id='soil')
+scene.add(halfspace)
+
+# position half space box in the subgrid. The  halfspace has to be positioned
+# manually because it traverses the grid. Grid traversal requries that objects extend
+# beyond the OS. Turning off autotranslate allows you to place objects beyond the OS.
+
+# PML seperation from the OS
+ps = ratio // 2 + 2
+# number of pml cells in the subgrid
+pc = 6
+# is os seperation
+isos = 3 * ratio
+
+h = antenna_p[2] - sg_z0 + (ps + pc + isos) * dl_s
+
+# half space box
+halfspace = gprMax.Box(p1=(0, 0, 0), p2=(411 * dl_s, 411 * dl_s, h), material_id='soil')
+# position the box using local coordinates3e8 / 400e6
+halfspace.autotranslate = False
+sg.add(halfspace)
+
+# Import the antenna model and add components to subgrid
+gssi_objects = antenna_like_GSSI_400(*antenna_p, resolution=dl_s)
+for obj in gssi_objects:
+    sg.add(obj)
+
+# half space box
+halfspace = gprMax.Box(p1=(0, 0, 0), p2=(x, y, antenna_p[2]), material_id='soil')
+scene.add(halfspace)
+
+for i in range(1, 51):
+    snap = gprMax.Snapshot(p1=(0, y / 2, 0), p2=(x, y / 2 + dl, z), dl=(dl, dl, dl),
+                           filename=Path(*parts[:-1], parts[-1] + '_' + str(i)).name,
+                           time=i * tw / 50)
+    scene.add(snap)
+
+# create a geometry view of the main grid and the sub grid stitched together
+gv = gprMax.GeometryView(p1=(0, 0, 0),
+                         p2=domain.props.p1,
+                         dl=dl,
+                         filename=fn.name,
+                         output_type='f',
+                         multi_grid=True)
+
+# create a geometry view of the main grid and the sub grid stitched together
+gv_normal = gprMax.GeometryView(p1=(0, 0, 0),
+                         p2=domain.props.p1,
+                         dl=dl,
+                         filename=Path(*parts[:-1], parts[-1] + '_voxels').name,
+                         output_type='n')
+# scene.add(gv)
+scene.add(gv_normal)
+
+# half space material
+layer_m = gprMax.Material(er=9, se=0.397e-3, mr=1, sm=0, id='soil_2')
+scene.add(layer_m)
+
+fb = gprMax.FractalBox(p1=(0, 0, 0), p2=(3, 1, 1), frac_dim=1.5, weighting=(1, 1, 1), n_materials=1, mixing_model_id='soil_2', id='fbox', seed=1)
+scene.add(fb)
+
+sr = gprMax.AddSurfaceRoughness(p1=(0, 0, 1), p2=(3, 1, 1), frac_dim=1.5, weighting=(1, 1), limits=(0.4, 1.2), fractal_box_id='fbox', seed=1)
+scene.add(sr)
+
+gprMax.run(scenes=[scene], n=1, geometry_only=False, outputfile=fn.parent, subgrid=True, autotranslate=True)