diff --git a/gprMax/cmds_geometry/cmds_geometry.py b/gprMax/cmds_geometry/cmds_geometry.py
index e102e8fd..55261b85 100644
--- a/gprMax/cmds_geometry/cmds_geometry.py
+++ b/gprMax/cmds_geometry/cmds_geometry.py
@@ -35,7 +35,7 @@ class UserObjectGeometry:
         self.autotranslate = True
 
     def __str__(self):
-        """Readble user string as per hash commands."""
+        """Readable string of parameters given to object."""
         s = ''
         for k, v in self.kwargs.items():
             if isinstance(v, tuple) or isinstance(v, list):
@@ -45,7 +45,7 @@ class UserObjectGeometry:
         return f'{self.hash}: {s[:-1]}'
 
     def params_str(self):
-        """Readble string of parameters given to object."""
+        """Readable user string as per hash commands."""
         return self.hash + ': ' + str(self.kwargs)
 
     def create(self, grid, uip):
diff --git a/gprMax/cmds_geometry/geometry_objects_read.py b/gprMax/cmds_geometry/geometry_objects_read.py
index 36b61d5a..a5a01597 100644
--- a/gprMax/cmds_geometry/geometry_objects_read.py
+++ b/gprMax/cmds_geometry/geometry_objects_read.py
@@ -25,8 +25,7 @@ import gprMax.config as config
 from .cmds_geometry import UserObjectGeometry
 from ..cython.geometry_primitives import build_voxels_from_array
 from ..exceptions import CmdInputError
-# from ..hash_cmds_file import check_cmd_names
-# from ..hash_cmds_multiuse import process_multicmds
+from ..hash_cmds_file import get_user_objects
 from ..utilities import round_value
 
 log = logging.getLogger(__name__)
@@ -69,11 +68,15 @@ class GeometryObjectsRead(UserObjectGeometry):
             # Strip out any newline characters and comments that must begin with double hashes
             materials = [line.rstrip() + '{' + matstr + '}\n' for line in f if(not line.startswith('##') and line.rstrip('\n'))]
 
-        # Check validity of command names
-        singlecmdsimport, multicmdsimport, geometryimport = check_cmd_names(materials, checkessential=False)
+        # build scene
+        # API for multiple scenes / model runs
+        scene = config.get_model_config().get_scene()
+        user_objs = get_user_objects(materials, check=False)
+        for user_obj in user_objs:
+            scene.add(user_obj)
 
-        # Process parameters for commands that can occur multiple times in the model
-        process_multicmds(multicmdsimport)
+        # Creates the internal simulation objects
+        # scene.create_internal_objects(self.G)
 
         # Update material type
         for material in G.materials:
diff --git a/gprMax/cmds_single_use.py b/gprMax/cmds_single_use.py
index cebb902b..5ede6f2a 100644
--- a/gprMax/cmds_single_use.py
+++ b/gprMax/cmds_single_use.py
@@ -272,9 +272,9 @@ class NumThreads(UserObjectSingle):
 
     def __str__(self):
         try:
-            return '#n_threads: {}'.format(self.kwargs['n'])
+            return f"#num_threads: {self.kwargs['n']}"
         except KeyError:
-            return '#n_threads:'
+            return '#num_threads:'
 
     def create(self, G, uip):
         try:
@@ -284,7 +284,7 @@ class NumThreads(UserObjectSingle):
         if n < 1:
             raise CmdInputError(self.__str__() + ' requires the value to be an integer not less than one')
 
-        config.model_configs[G.model_num].ompthreads = set_omp_threads(n)
+        config.get_model_config().ompthreads = set_omp_threads(n)
 
 
 class TimeStepStabilityFactor(UserObjectSingle):
@@ -386,9 +386,7 @@ class SrcSteps(UserObjectSingle):
         except KeyError:
             raise CmdInputError('#src_steps: requires exactly three parameters')
 
-        log.info(f'Simple sources will step {G.srcsteps[0] * G.dx:g}m, \
-                 {G.srcsteps[1] * G.dy:g}m, {G.srcsteps[2] * G.dz:g}m \
-                 for each model run.')
+        log.info(f'Simple sources will step {G.srcsteps[0] * G.dx:g}m, {G.srcsteps[1] * G.dy:g}m, {G.srcsteps[2] * G.dz:g}m for each model run.')
 
 
 class RxSteps(UserObjectSingle):
@@ -409,9 +407,7 @@ class RxSteps(UserObjectSingle):
         except KeyError:
             raise CmdInputError('#rx_steps: requires exactly three parameters')
 
-        log.info(f'All receivers will step {G.rxsteps[0] * G.dx:g}m, \
-                 {G.rxsteps[1] * G.dy:g}m, {G.rxsteps[2] * G.dz:g}m \
-                 for each model run.')
+        log.info(f'All receivers will step {G.rxsteps[0] * G.dx:g}m, {G.rxsteps[1] * G.dy:g}m, {G.rxsteps[2] * G.dz:g}m for each model run.')
 
 
 class ExcitationFile(UserObjectSingle):
@@ -509,7 +505,7 @@ class OutputDir(UserObjectSingle):
         self.order = 12
 
     def create(self, grid, uip):
-        config.model_configs[grid.model_num].set_output_file_path(self.kwargs['dir'])
+        config.get_model_config().set_output_file_path(self.kwargs['dir'])
 
 
 class NumberOfModelRuns(UserObjectSingle):
diff --git a/gprMax/config.py b/gprMax/config.py
index 7bdf6f4e..d890c664 100644
--- a/gprMax/config.py
+++ b/gprMax/config.py
@@ -43,18 +43,22 @@ sim_config = None
 # Instance of ModelConfig that hold model configuration parameters.
 model_configs = []
 
+# Each model in a simulation is given a unique number when the instance of
+# ModelConfig is created
+model_num = 0
+
+def get_model_config():
+    """Return ModelConfig instace for specific model."""
+    return model_configs[model_num]
+
+
 class ModelConfig:
     """Configuration parameters for a model.
         N.B. Multiple models can exist within a simulation
     """
 
-    def __init__(self, model_num):
-        """
-        Args:
-            model_num (int): Model number.
-        """
+    def __init__(self):
 
-        self.i = model_num # Indexed from 0
         self.mode = '3D'
         self.grids = []
         self.ompthreads = None # Number of OpenMP threads
@@ -75,11 +79,11 @@ class ModelConfig:
         self.reuse_geometry = False
 
         # String to print at start of each model run
-        inputfilestr = f'\n--- Model {self.i + 1}/{sim_config.model_end}, input file: {sim_config.input_file_path}'
+        inputfilestr = f'\n--- Model {model_num + 1}/{sim_config.model_end}, input file: {sim_config.input_file_path}'
         self.set_inputfilestr(inputfilestr)
 
         if not sim_config.single_model:
-            self.appendmodelnumber = str(self.i + 1) # Indexed from 1
+            self.appendmodelnumber = str(model_num + 1) # Indexed from 1
         else:
             self.appendmodelnumber = ''
 
@@ -108,7 +112,7 @@ class ModelConfig:
 
     def get_scene(self):
         if sim_config.scenes:
-            return sim_config.scenes[self.i]
+            return sim_config.scenes[model_num]
         else: return None
 
     def get_usernamespace(self):
@@ -116,8 +120,8 @@ class ModelConfig:
                 'e0': e0, # Permittivity of free space (F/m)
                 'm0': m0, # Permeability of free space (H/m)
                 'z0': np.sqrt(m0 / e0), # Impedance of free space (Ohms)
-                'number_model_runs': sim_config.model_end + 1,
-                'current_model_run': self.i + 1,
+                'number_model_runs': sim_config.model_end,
+                'current_model_run': model_num + 1,
                 'inputfile': sim_config.input_file_path.resolve()}
 
     def set_inputfilestr(self, inputfilestr):
@@ -169,7 +173,6 @@ class SimulationConfig:
         """
 
         self.args = args
-        log.debug('Fix parsing args')
 
         # General settings for the simulation
         #   inputfilepath: path to inputfile location
@@ -177,17 +180,15 @@ class SimulationConfig:
         #   messages: whether to print all messages as output to stdout or not
         #   progressbars: whether to show progress bars on stdoout or not
         #   cpu, cuda, opencl: solver type
+        #   subgrid: whether the simulation uses sub-grids
         #   precision: data type for electromagnetic field output (single/double)
-        #   autotranslate: 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.
         self.general = {'messages': True,
                         'progressbars': True,
                         'cpu': True,
                         'cuda': False,
                         'opencl': False,
-                        'precision': 'single'}
+                        'subgrid': False,
+                        'precision': 'double'}
 
         self.em_consts = {'c': c, # Speed of light in free space (m/s)
                           'e0': e0, # Permittivity of free space (F/m)
@@ -217,9 +218,9 @@ class SimulationConfig:
 
         # Subgrid parameter may not exist if user enters via CLI
         try:
-            self.subgrid = args.subgrid
+            self.general['subgrid'] = self.args.subgrid
         except AttributeError:
-            self.subgrid = False
+            self.general['subgrid'] = False
 
         # Scenes parameter may not exist if user enters via CLI
         try:
diff --git a/gprMax/config_parser.py b/gprMax/config_parser.py
index 23f90ea7..bb066b67 100755
--- a/gprMax/config_parser.py
+++ b/gprMax/config_parser.py
@@ -34,14 +34,11 @@ def write_simulation_config(args):
         config.sim_config = config.SimulationConfig(args)
 
 
-def write_model_config(model_num):
+def write_model_config():
     """Write model level configuration parameters to config module. As there can
         only be one instance of the config module objects are always found via
         'import gprMax.config'
-
-    Args:
-        model_num (int): Model number.
     """
 
-    model_config = config.ModelConfig(model_num)
+    model_config = config.ModelConfig()
     config.model_configs.append(model_config)
diff --git a/gprMax/contexts.py b/gprMax/contexts.py
index 04edd8d9..dfe08bd8 100644
--- a/gprMax/contexts.py
+++ b/gprMax/contexts.py
@@ -35,7 +35,7 @@ log = logging.getLogger(__name__)
 
 class Context:
     """Generic context for the model to run in. Sub-class with specific contexts
-    e.g. an MPI context.
+        e.g. an MPI context.
     """
 
     def __init__(self):
@@ -90,16 +90,17 @@ class NoMPIContext(Context):
         """Specialise how the models are farmed out."""
 
         for i in self.model_range:
-            write_model_config(i)
+            config.model_num = i
+            write_model_config()
             # Always create a solver for the first model.
             # The next model to run only gets a new solver if the
             # geometry is not re-used.
             if i != 0 and config.sim_config.args.geometry_fixed:
-                config.model_configs[i].reuse_geometry = True
+                config.get_model_config().reuse_geometry = True
                 # Ensure re-used G is associated correctly with model
-                G.model_num = i
+                # G.model_num = i
             else:
-                G = create_G(i)
+                G = create_G()
 
             model = ModelBuildRun(G)
             model.build()
diff --git a/gprMax/fields_outputs.py b/gprMax/fields_outputs.py
index 69212ed8..1a8ce898 100644
--- a/gprMax/fields_outputs.py
+++ b/gprMax/fields_outputs.py
@@ -98,9 +98,12 @@ __global__ void store_outputs(int NRX, int iteration, const int* __restrict__ rx
 """)
 
 
-def write_hdf5_outputfile(outputfile, G):
-    write_hdf5_main_grid_outputfile(outputfile, G)
-    write_hdf5_sub_grid_outputfile(outputfile, G)
+def write_hdf5_outputfiles(outputfile, G):
+    if G.rxs:
+        write_hdf5_main_grid_outputfile(outputfile, G)
+    sg_rxs = [True for sg in G.subgrids if sg.rxs]
+    if sg_rxs:
+        write_hdf5_sub_grid_outputfile(outputfile, G)
 
 
 def write_hdf5_main_grid_outputfile(outputfile, G):
diff --git a/gprMax/fractals.py b/gprMax/fractals.py
index 24b3d836..6df75a02 100644
--- a/gprMax/fractals.py
+++ b/gprMax/fractals.py
@@ -93,7 +93,7 @@ class FractalSurface:
         A = fftpack.fftshift(A)
 
         # Generate fractal
-        generate_fractal2D(surfacedims[0], surfacedims[1], config.model_configs[G.model_num].ompthreads,
+        generate_fractal2D(surfacedims[0], surfacedims[1], config.get_model_config().ompthreads,
                            self.b, self.weighting, v1, A, self.fractalsurface)
         # Shift the zero frequency component to start of the array
         self.fractalsurface = fftpack.ifftshift(self.fractalsurface)
@@ -180,7 +180,7 @@ class FractalVolume:
         A = fftpack.fftshift(A)
 
         # Generate fractal
-        generate_fractal3D(self.nx, self.ny, self.nz, config.model_configs[G.model_num].ompthreads,
+        generate_fractal3D(self.nx, self.ny, self.nz, config.get_model_config().ompthreads,
                            self.b, self.weighting, v1, A, self.fractalvolume)
 
         # Shift the zero frequency component to the start of the array
diff --git a/gprMax/geometry_outputs.py b/gprMax/geometry_outputs.py
index b72a1aa8..54025bce 100644
--- a/gprMax/geometry_outputs.py
+++ b/gprMax/geometry_outputs.py
@@ -113,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.parts
-        self.filename = Path(*parts[:-1], self.filename + config.model_configs[self.G.model_num].appendmodelnumber)
+        parts = config.get_model_config().output_file_path.parts
+        self.filename = Path(*parts[:-1], self.filename + config.get_model_config().appendmodelnumber)
         self.filename = self.filename.with_suffix(self.fileext)
 
     def write_vtk(self, G, pbar):
@@ -452,8 +452,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].output_file_path.parts
-        self.filename = Path(*parts[:-1], self.filename + config.model_configs[self.G.model_num].appendmodelnumber)
+        parts = config.get_model_config().output_file_path.parts
+        self.filename = Path(*parts[:-1], self.filename + config.get_model_config().appendmodelnumber)
         self.filename = self.filename.with_suffix(self.fileext)
 
     def write_vtk(self, *args):
diff --git a/gprMax/grid.py b/gprMax/grid.py
index 474ae3d4..1170c7bc 100644
--- a/gprMax/grid.py
+++ b/gprMax/grid.py
@@ -41,15 +41,9 @@ class FDTDGrid:
         accessing regularly used parameters.
     """
 
-    def __init__(self, model_num):
-        """
-        Args:
-            model_num (int): Model number.
-        """
-
+    def __init__(self):
         self.title = ''
         self.name = 'main_grid'
-        self.model_num = model_num
 
         self.nx = 0
         self.ny = 0
@@ -170,26 +164,26 @@ class FDTDGrid:
 
     def initialise_dispersive_arrays(self):
         """Initialise field arrays when there are dispersive materials present."""
-        self.Tx = np.zeros((config.model_configs[self.model_num].materials['maxpoles'],
-                            self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.model_configs[self.model_num].materials['dispersivedtype'])
-        self.Ty = np.zeros((config.model_configs[self.model_num].materials['maxpoles'],
-                            self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.model_configs[self.model_num].materials['dispersivedtype'])
-        self.Tz = np.zeros((config.model_configs[self.model_num].materials['maxpoles'],
-                            self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.model_configs[self.model_num].materials['dispersivedtype'])
+        self.Tx = np.zeros((config.get_model_config().materials['maxpoles'],
+                            self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.get_model_config().materials['dispersivedtype'])
+        self.Ty = np.zeros((config.get_model_config().materials['maxpoles'],
+                            self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.get_model_config().materials['dispersivedtype'])
+        self.Tz = np.zeros((config.get_model_config().materials['maxpoles'],
+                            self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.get_model_config().materials['dispersivedtype'])
 
     def initialise_dispersive_update_coeff_array(self):
         """Initialise array for storing update coefficients when there are dispersive
             materials present.
         """
         self.updatecoeffsdispersive = np.zeros((len(self.materials), 3 *
-                                                config.model_configs[self.model_num].materials['maxpoles']),
-                                                dtype=config.model_configs[self.model_num].materials['dispersivedtype'])
+                                                config.get_model_config().materials['maxpoles']),
+                                                dtype=config.get_model_config().materials['dispersivedtype'])
 
     def reset_fields(self):
         """Clear arrays for field components and PMLs."""
         # Clear arrays for field components
         self.initialise_field_arrays()
-        if config.model_configs[self.model_num].materials['maxpoles'] > 0:
+        if config.get_model_config().materials['maxpoles'] > 0:
             self.initialise_dispersive_arrays()
 
         # Clear arrays for fields in PML
@@ -242,9 +236,9 @@ class FDTDGrid:
             mem_use (int): Memory (bytes).
         """
 
-        mem_use = int(3 * config.model_configs[self.model_num].materials['maxpoles'] *
+        mem_use = int(3 * config.get_model_config().materials['maxpoles'] *
                        (self.nx + 1) * (self.ny + 1) * (self.nz + 1) *
-                       np.dtype(config.model_configs[self.model_num].materials['dispersivedtype']).itemsize)
+                       np.dtype(config.get_model_config().materials['dispersivedtype']).itemsize)
         return mem_use
 
     def tmx(self):
@@ -293,8 +287,7 @@ class FDTDGrid:
 class CUDAGrid(FDTDGrid):
     """Additional grid methods for solving on GPU using CUDA."""
 
-    def __init__(self, model_num):
-        super().__init__(model_num)
+    def __init__(self):
 
         # Threads per block - used for main electric/magnetic field updates
         self.tpb = (128, 1, 1)
@@ -345,7 +338,7 @@ class CUDAGrid(FDTDGrid):
         # Copy arrays geometry and field component arrays
         self.htod_geometry_arrays()
         self.htod_field_arrays()
-        if config.model_configs[self.model_num].materials['maxpoles'] > 0:
+        if config.get_model_config().materials['maxpoles'] > 0:
             self.htod_dispersive_arrays()
 
         # Copy arrays for fields in PML
@@ -404,7 +397,7 @@ def dispersion_analysis(G):
 
                     # Set maximum frequency to a threshold drop from maximum power, ignoring DC value
                     try:
-                        freqthres = np.where(power[freqmaxpower:] < -config.model_configs[G.model_num].numdispersion['highestfreqthres'])[0][0] + freqmaxpower
+                        freqthres = np.where(power[freqmaxpower:] < -config.get_model_config().numdispersion['highestfreqthres'])[0][0] + freqmaxpower
                         results['maxfreq'].append(freqs[freqthres])
                     except ValueError:
                         results['error'] = 'unable to calculate maximum power from waveform, most likely due to undersampling.'
@@ -445,9 +438,9 @@ def dispersion_analysis(G):
         minwavelength = minvelocity / results['maxfreq']
 
         # Maximum spatial step
-        if '3D' in config.model_configs[G.model_num].mode:
+        if '3D' in config.get_model_config().mode:
             delta = max(G.dx, G.dy, G.dz)
-        elif '2D' in config.model_configs[G.model_num].mode:
+        elif '2D' in config.get_model_config().mode:
             if G.nx == 1:
                 delta = max(G.dy, G.dz)
             elif G.ny == 1:
@@ -462,7 +455,7 @@ def dispersion_analysis(G):
         results['N'] = minwavelength / delta
 
         # Check grid sampling will result in physical wave propagation
-        if int(np.floor(results['N'])) >= config.model_configs[G.model_num].numdispersion['mingridsampling']:
+        if int(np.floor(results['N'])) >= config.get_model_config().numdispersion['mingridsampling']:
             # Numerical phase velocity
             vp = np.pi / (results['N'] * np.arcsin((1 / S) * np.sin((np.pi * S) / results['N'])))
 
diff --git a/gprMax/hash_cmds_file.py b/gprMax/hash_cmds_file.py
index 745bb4de..ad4b668c 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.')
+            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.\n')
             # String to hold Python code to be executed
             pythoncode = ''
             x += 1
@@ -168,7 +168,7 @@ def write_processed_file(processedlines, G):
         G (FDTDGrid): Parameters describing a grid in a model.
     """
 
-    parts = config.model_configs[G.model_num].output_file_path.parts
+    parts = config.get_model_config().output_file_path.parts
     processedfile = (Path(*parts[:-1], parts[-1] + '_processed.in'))
 
     with open(processedfile, 'w') as f:
@@ -295,7 +295,7 @@ def parse_hash_commands(scene, G):
 
     with open(config.sim_config.input_file_path) as inputfile:
 
-        usernamespace = config.model_configs[G.model_num].get_usernamespace()
+        usernamespace = config.get_model_config().get_usernamespace()
 
         # Read input file and process any Python and include file commands
         processedlines = process_python_include_code(inputfile, usernamespace)
diff --git a/gprMax/materials.py b/gprMax/materials.py
index 7dc8ac84..bd369bb5 100644
--- a/gprMax/materials.py
+++ b/gprMax/materials.py
@@ -131,18 +131,18 @@ class DispersiveMaterial(Material):
 
         # The implementation of the dispersive material modelling comes from the
         # derivation in: http://dx.doi.org/10.1109/TAP.2014.2308549
-        self.w = np.zeros(config.model_configs[G.model_num].materials['maxpoles'],
-                          dtype=config.model_configs[G.model_num].materials['dispersivedtype'])
-        self.q = np.zeros(config.model_configs[G.model_num].materials['maxpoles'],
-                          dtype=config.model_configs[G.model_num].materials['dispersivedtype'])
-        self.zt = np.zeros(config.model_configs[G.model_num].materials['maxpoles'],
-                           dtype=config.model_configs[G.model_num].materials['dispersivedtype'])
-        self.zt2 = np.zeros(config.model_configs[G.model_num].materials['maxpoles'],
-                            dtype=config.model_configs[G.model_num].materials['dispersivedtype'])
-        self.eqt = np.zeros(config.model_configs[G.model_num].materials['maxpoles'],
-                            dtype=config.model_configs[G.model_num].materials['dispersivedtype'])
-        self.eqt2 = np.zeros(config.model_configs[G.model_num].materials['maxpoles'],
-                            dtype=config.model_configs[G.model_num].materials['dispersivedtype'])
+        self.w = np.zeros(config.get_model_config().materials['maxpoles'],
+                          dtype=config.get_model_config().materials['dispersivedtype'])
+        self.q = np.zeros(config.get_model_config().materials['maxpoles'],
+                          dtype=config.get_model_config().materials['dispersivedtype'])
+        self.zt = np.zeros(config.get_model_config().materials['maxpoles'],
+                           dtype=config.get_model_config().materials['dispersivedtype'])
+        self.zt2 = np.zeros(config.get_model_config().materials['maxpoles'],
+                            dtype=config.get_model_config().materials['dispersivedtype'])
+        self.eqt = np.zeros(config.get_model_config().materials['maxpoles'],
+                            dtype=config.get_model_config().materials['dispersivedtype'])
+        self.eqt2 = np.zeros(config.get_model_config().materials['maxpoles'],
+                            dtype=config.get_model_config().materials['dispersivedtype'])
 
         for x in range(self.poles):
             if 'debye' in self.type:
@@ -223,7 +223,7 @@ def process_materials(G):
                                 print a table.
     """
 
-    if config.model_configs[G.model_num].materials['maxpoles'] == 0:
+    if config.get_model_config().materials['maxpoles'] == 0:
         materialsdata = [['\nID', '\nName', '\nType', '\neps_r', 'sigma\n[S/m]',
                           '\nmu_r', 'sigma*\n[Ohm/m]', 'Dielectric\nsmoothable']]
     else:
@@ -243,7 +243,7 @@ def process_materials(G):
         # Add update coefficients to overall storage for dispersive materials
         if hasattr(material, 'poles'):
             z = 0
-            for pole in range(config.model_configs[G.model_num].materials['maxpoles']):
+            for pole in range(config.get_model_config().materials['maxpoles']):
                 G.updatecoeffsdispersive[material.numID, z:z + 3] = (config.sim_config.em_consts['e0'] *
                                                                      material.eqt2[pole], material.eqt[pole], material.zt[pole])
                 z += 3
@@ -255,7 +255,7 @@ def process_materials(G):
         materialtext.append(material.type)
         materialtext.append(f'{material.er:g}')
         materialtext.append(f'{material.se:g}')
-        if config.model_configs[G.model_num].materials['maxpoles'] > 0:
+        if config.get_model_config().materials['maxpoles'] > 0:
             if 'debye' in material.type:
                 materialtext.append('\n'.join('{:g}'.format(deltaer) for deltaer in material.deltaer))
                 materialtext.append('\n'.join('{:g}'.format(tau) for tau in material.tau))
@@ -367,8 +367,8 @@ class PeplinskiSoil:
                 m.type = 'debye'
                 m.averagable = False
                 m.poles = 1
-                if m.poles > config.model_configs[G.model_num].materials['maxpoles']:
-                    config.model_configs[G.model_num].materials['maxpoles'] = m.poles
+                if m.poles > config.get_model_config().materials['maxpoles']:
+                    config.get_model_config().materials['maxpoles'] = m.poles
                 m.er = eri
                 m.se = sig
                 m.deltaer.append(er - eri)
diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py
index 44e69982..9a6fe6e0 100644
--- a/gprMax/model_build_run.py
+++ b/gprMax/model_build_run.py
@@ -35,7 +35,7 @@ import gprMax.config as config
 from .cython.yee_cell_build import build_electric_components
 from .cython.yee_cell_build import build_magnetic_components
 from .exceptions import GeneralError
-from .fields_outputs import write_hdf5_outputfile
+from .fields_outputs import write_hdf5_outputfiles
 from .grid import dispersion_analysis
 from .hash_cmds_file import parse_hash_commands
 from .materials import Material
@@ -65,7 +65,7 @@ class ModelBuildRun:
         # used with threaded model building methods, e.g. fractals. Can be
         # changed by #num_threads command in input file or via API later for
         # use with CPU solver.
-        config.model_configs[self.G.model_num].ompthreads = set_omp_threads(config.model_configs[self.G.model_num].ompthreads)
+        config.get_model_config().ompthreads = set_omp_threads(config.get_model_config().ompthreads)
 
     def build(self):
         """Builds the Yee cells for a model."""
@@ -76,14 +76,14 @@ class ModelBuildRun:
         self.p = psutil.Process()
 
         # Normal model reading/building process; bypassed if geometry information to be reused
-        self.reuse_geometry() if config.model_configs[G.model_num].reuse_geometry else self.build_geometry()
+        self.reuse_geometry() if config.get_model_config().reuse_geometry else self.build_geometry()
 
-        log.info(f'\nOutput directory: {config.model_configs[G.model_num].output_file_path.parent.resolve()}')
+        log.info(f'\nOutput directory: {config.get_model_config().output_file_path.parent.resolve()}')
 
         # Adjust position of simple sources and receivers if required
         if G.srcsteps[0] != 0 or G.srcsteps[1] != 0 or G.srcsteps[2] != 0:
             for source in itertools.chain(G.hertziandipoles, G.magneticdipoles):
-                if G.model_num == 0:
+                if config.model_num == 0:
                     if (source.xcoord + G.srcsteps[0] * config.sim_config.model_end < 0 or
                         source.xcoord + G.srcsteps[0] * config.sim_config.model_end > G.nx or
                         source.ycoord + G.srcsteps[1] * config.sim_config.model_end < 0 or
@@ -91,12 +91,12 @@ class ModelBuildRun:
                         source.zcoord + G.srcsteps[2] * config.sim_config.model_end < 0 or
                         source.zcoord + G.srcsteps[2] * config.sim_config.model_end > G.nz):
                         raise GeneralError('Source(s) will be stepped to a position outside the domain.')
-                source.xcoord = source.xcoordorigin + G.model_num * G.srcsteps[0]
-                source.ycoord = source.ycoordorigin + G.model_num * G.srcsteps[1]
-                source.zcoord = source.zcoordorigin + G.model_num * G.srcsteps[2]
+                source.xcoord = source.xcoordorigin + config.model_num * G.srcsteps[0]
+                source.ycoord = source.ycoordorigin + config.model_num * G.srcsteps[1]
+                source.zcoord = source.zcoordorigin + config.model_num * G.srcsteps[2]
         if G.rxsteps[0] != 0 or G.rxsteps[1] != 0 or G.rxsteps[2] != 0:
             for receiver in G.rxs:
-                if G.model_num == 0:
+                if config.model_num == 0:
                     if (receiver.xcoord + G.rxsteps[0] * config.sim_config.model_end < 0 or
                         receiver.xcoord + G.rxsteps[0] * config.sim_config.model_end > G.nx or
                         receiver.ycoord + G.rxsteps[1] * config.sim_config.model_end < 0 or
@@ -104,9 +104,9 @@ class ModelBuildRun:
                         receiver.zcoord + G.rxsteps[2] * config.sim_config.model_end < 0 or
                         receiver.zcoord + G.rxsteps[2] * config.sim_config.model_end > G.nz):
                         raise GeneralError('Receiver(s) will be stepped to a position outside the domain.')
-                receiver.xcoord = receiver.xcoordorigin + G.model_num * G.rxsteps[0]
-                receiver.ycoord = receiver.ycoordorigin + G.model_num * G.rxsteps[1]
-                receiver.zcoord = receiver.zcoordorigin + G.model_num * G.rxsteps[2]
+                receiver.xcoord = receiver.xcoordorigin + config.model_num * G.rxsteps[0]
+                receiver.ycoord = receiver.ycoordorigin + config.model_num * G.rxsteps[1]
+                receiver.zcoord = receiver.zcoordorigin + config.model_num * G.rxsteps[2]
 
         # Write files for any geometry views and geometry object outputs
         if not (G.geometryviews or G.geometryobjectswrite) and config.sim_config.args.geometry_only:
@@ -132,7 +132,7 @@ class ModelBuildRun:
     def build_geometry(self):
         G = self.G
 
-        log.info(config.model_configs[G.model_num].inputfilestr)
+        log.info(config.get_model_config().inputfilestr)
 
         scene = self.build_scene()
 
@@ -152,7 +152,7 @@ class ModelBuildRun:
             gb.update_voltage_source_materials()
             gb.grid.initialise_field_arrays()
             gb.grid.initialise_std_update_coeff_arrays()
-            if config.model_configs[gb.grid.model_num].materials['maxpoles'] > 0:
+            if config.get_model_config().materials['maxpoles'] > 0:
                 gb.grid.initialise_dispersive_arrays()
                 gb.grid.initialise_dispersive_update_coeff_array()
             gb.build_materials()
@@ -161,10 +161,10 @@ class ModelBuildRun:
             results = dispersion_analysis(gb.grid)
             if results['error']:
                 log.warning(Fore.RED + f"\nNumerical dispersion analysis ({gb.grid.name}) not carried out as {results['error']}" + Style.RESET_ALL)
-            elif results['N'] < config.model_configs[gb.grid.model_num].numdispersion['mingridsampling']:
+            elif results['N'] < config.get_model_config().numdispersion['mingridsampling']:
                 raise GeneralError(f"\nNon-physical wave propagation in {gb.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[gb.grid.model_num].numdispersion['maxnumericaldisp']):
+                  config.get_model_config().numdispersion['maxnumericaldisp']):
                 log.warning(Fore.RED + f"\n{gb.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 ({gb.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")
@@ -172,14 +172,14 @@ class ModelBuildRun:
     def reuse_geometry(self):
         # Reset iteration number
         self.G.iteration = 0
-        config.model_configs[self.G.model_num].set_inputfilestr(f'\n--- Model {config.model_configs[self.G.model_num].appendmodelnumber}/{config.sim_config.model_end}, input file (not re-processed, i.e. geometry fixed): {config.sim_config.input_file_path}')
-        log.info(config.model_configs[self.G.model_num].inputfilestr)
+        config.get_model_config().set_inputfilestr(f'\n--- Model {config.get_model_config().appendmodelnumber}/{config.sim_config.model_end}, input file (not re-processed, i.e. geometry fixed): {config.sim_config.input_file_path}')
+        log.info(config.get_model_config().inputfilestr)
         for grid in [self.G] + self.G.subgrids:
             grid.reset_fields()
 
     def build_scene(self):
         # API for multiple scenes / model runs
-        scene = config.model_configs[self.G.model_num].get_scene()
+        scene = config.get_model_config().get_scene()
 
         # If there is no scene, process the hash commands
         if not scene:
@@ -197,14 +197,14 @@ class ModelBuildRun:
             to file(s).
         """
 
-        # Write an output file in HDF5 format
-        write_hdf5_outputfile(config.model_configs[self.G.model_num].output_file_path_ext, self.G)
+        # Write an output file(s) in HDF5 format
+        write_hdf5_outputfiles(config.get_model_config().output_file_path_ext, self.G)
 
         # Write any snapshots to file
         if self.G.snapshots:
             # Create directory for snapshots
-            config.model_configs[self.G.model_num].set_snapshots_file_path()
-            snapshotdir = config.model_configs[self.G.model_num].snapshot_file_path
+            config.get_model_config().set_snapshots_file_path()
+            snapshotdir = config.get_model_config().snapshot_file_path
             snapshotdir.mkdir(exist_ok=True)
 
             log.info('')
@@ -244,16 +244,16 @@ class ModelBuildRun:
 
         # 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')
-            if config.model_configs[self.G.model_num].ompthreads > config.sim_config.hostinfo['physicalcores']:
-                log.warning(Fore.RED + f"You have specified more threads ({config.model_configs[self.G.model_num].ompthreads}) than available physical CPU cores ({config.sim_config.hostinfo['physicalcores']}). This may lead to degraded performance." + Style.RESET_ALL)
+            log.info(f'CPU (OpenMP) threads for solving: {config.get_model_config().ompthreads}\n')
+            if config.get_model_config().ompthreads > config.sim_config.hostinfo['physicalcores']:
+                log.warning(Fore.RED + f"You have specified more threads ({config.get_model_config().ompthreads}) than available physical CPU cores ({config.sim_config.hostinfo['physicalcores']}). This may lead to degraded performance." + Style.RESET_ALL)
         # Print information about any GPU in use
         elif config.sim_config.general['cuda']:
-            log.info(f"GPU for solving: {config.model_configs[self.G.model_num].cuda['gpu'].deviceID} - {config.model_configs[self.G.model_num].cuda['gpu'].name}\n")
+            log.info(f"GPU for solving: {config.get_model_config().cuda['gpu'].deviceID} - {config.get_model_config().cuda['gpu'].name}\n")
 
         # Prepare iterator
         if config.sim_config.is_messages():
-            iterator = tqdm(range(self.G.iterations), desc=f'Running model {self.G.model_num + 1}/{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 {config.model_num + 1}/{config.sim_config.model_end}', ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.sim_config.general['progressbars'])
         else:
             iterator = range(self.G.iterations)
 
@@ -298,11 +298,11 @@ class GridBuilder:
         pbar.close()
 
     def tm_grid_update(self):
-        if '2D TMx' == config.model_configs[self.grid.model_num].mode:
+        if '2D TMx' == config.get_model_config().mode:
             self.grid.tmx()
-        elif '2D TMy' == config.model_configs[self.grid.model_num].mode:
+        elif '2D TMy' == config.get_model_config().mode:
             self.grid.tmy()
-        elif '2D TMz' == config.model_configs[self.grid.model_num].mode:
+        elif '2D TMz' == config.get_model_config().mode:
             self.grid.tmz()
 
     def update_voltage_source_materials(self):
diff --git a/gprMax/solvers.py b/gprMax/solvers.py
index 078d830f..1de5a38b 100644
--- a/gprMax/solvers.py
+++ b/gprMax/solvers.py
@@ -24,22 +24,17 @@ from .updates import CPUUpdates
 from .updates import CUDAUpdates
 
 
-def create_G(model_num):
+def create_G():
     """Create grid object according to solver.
 
-    Args:
-        model_num (int): Model number.
-
     Returns:
         G (FDTDGrid): Holds essential parameters describing the model.
     """
 
     if config.sim_config.general['cpu']:
-        G = FDTDGrid(model_num)
+        G = FDTDGrid()
     elif config.sim_config.general['cuda']:
-        G = CUDAGrid(model_num)
-    elif config.sim_config.subgrid:
-        G = FDTDGrid(model_num)
+        G = CUDAGrid()
 
     return G
 
@@ -54,15 +49,7 @@ def create_solver(G):
         solver (Solver): solver object.
     """
 
-    if config.sim_config.general['cpu']:
-        updates = CPUUpdates(G)
-        solver = Solver(updates)
-        props = updates.adapt_dispersive_config()
-        updates.set_dispersive_updates(props)
-    elif config.sim_config.general['cuda']:
-        updates = CUDAUpdates(G)
-        solver = Solver(updates)
-    elif config.sim_config.subgrid:
+    if config.sim_config.general['subgrid']:
         updates = create_subgrid_updates(G)
         solver = Solver(updates, hsg=True)
         # A large range of different functions exist to advance the time step for
@@ -70,6 +57,15 @@ def create_solver(G):
         # the required numerical precision and dispersive material type.
         props = updates.adapt_dispersive_config()
         updates.set_dispersive_updates(props)
+    elif config.sim_config.general['cpu']:
+        updates = CPUUpdates(G)
+        solver = Solver(updates)
+        props = updates.adapt_dispersive_config()
+        updates.set_dispersive_updates(props)
+    elif config.sim_config.general['cuda']:
+        updates = CUDAUpdates(G)
+        solver = Solver(updates)
+
 
     return solver
 
diff --git a/gprMax/subgrids/subgrid_hsg.py b/gprMax/subgrids/subgrid_hsg.py
index ea944c4f..56a4239c 100644
--- a/gprMax/subgrids/subgrid_hsg.py
+++ b/gprMax/subgrids/subgrid_hsg.py
@@ -19,6 +19,7 @@
 from colorama import init, Fore, Style
 init()
 
+import gprMax.config as config
 from .base import SubGridBase
 from ..cython.fields_updates_hsg import cython_update_is
 from ..cython.fields_updates_hsg import cython_update_magnetic_os
@@ -45,16 +46,16 @@ class SubGridHSG(SubGridBase):
         # Hy = c0Hy - c3Ex + c1Ez
         # Hx = c0Hx - c2Ez + c3Ey
         # bottom and top
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx, self.nwy + 1, self.nwz, 1, self.Hy, precursors.ex_bottom, precursors.ex_top, self.IDlookup['Hy'], 1, -1, 3, self.nthreads)
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx + 1, self.nwy, self.nwz, 1, self.Hx, precursors.ey_bottom, precursors.ey_top, self.IDlookup['Hx'], -1, 1, 3, self.nthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx, self.nwy + 1, self.nwz, 1, self.Hy, precursors.ex_bottom, precursors.ex_top, self.IDlookup['Hy'], 1, -1, 3, config.get_model().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx + 1, self.nwy, self.nwz, 1, self.Hx, precursors.ey_bottom, precursors.ey_top, self.IDlookup['Hx'], -1, 1, 3, config.get_model().ompthreads)
 
         # left and right
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwy, self.nwz + 1, self.nwx, 2, self.Hz, precursors.ey_left, precursors.ey_right, self.IDlookup['Hz'], 1, -1, 1, self.nthreads)
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwy + 1, self.nwz, self.nwx, 2, self.Hy, precursors.ez_left, precursors.ez_right, self.IDlookup['Hy'], -1, 1, 1, self.nthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwy, self.nwz + 1, self.nwx, 2, self.Hz, precursors.ey_left, precursors.ey_right, self.IDlookup['Hz'], 1, -1, 1, config.get_model().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwy + 1, self.nwz, self.nwx, 2, self.Hy, precursors.ez_left, precursors.ez_right, self.IDlookup['Hy'], -1, 1, 1, config.get_model().ompthreads)
 
         # front and back
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx, self.nwz + 1, self.nwy, 3, self.Hz, precursors.ex_front, precursors.ex_back, self.IDlookup['Hz'], -1, 1, 2, self.nthreads)
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx + 1, self.nwz, self.nwy, 3, self.Hx, precursors.ez_front, precursors.ez_back, self.IDlookup['Hx'], 1, -1, 2, self.nthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx, self.nwz + 1, self.nwy, 3, self.Hz, precursors.ex_front, precursors.ex_back, self.IDlookup['Hz'], -1, 1, 2, config.get_model().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx + 1, self.nwz, self.nwy, 3, self.Hx, precursors.ez_front, precursors.ez_back, self.IDlookup['Hx'], 1, -1, 2, config.get_model().ompthreads)
 
     def update_electric_is(self, precursors):
         # Args: nwl, nwm, nwn, face, field, inc_field, lookup_id, sign, mod, co
@@ -64,16 +65,16 @@ class SubGridHSG(SubGridBase):
         # Ez = c0(Ez) + c1(dHy) - c2(dHx)
 
         # bottom and top
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx, self.nwy + 1, self.nwz, 1, self.Ex, precursors.hy_bottom, precursors.hy_top, self.IDlookup['Ex'], 1, -1, 3, self.nthreads)
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx + 1, self.nwy, self.nwz, 1, self.Ey, precursors.hx_bottom, precursors.hx_top, self.IDlookup['Ey'], -1, 1, 3, self.nthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx, self.nwy + 1, self.nwz, 1, self.Ex, precursors.hy_bottom, precursors.hy_top, self.IDlookup['Ex'], 1, -1, 3, config.get_model().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx + 1, self.nwy, self.nwz, 1, self.Ey, precursors.hx_bottom, precursors.hx_top, self.IDlookup['Ey'], -1, 1, 3, config.get_model().ompthreads)
 
         # left and right
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwy, self.nwz + 1, self.nwx, 2, self.Ey, precursors.hz_left, precursors.hz_right, self.IDlookup['Ey'], 1, -1, 1, self.nthreads)
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwy + 1, self.nwz, self.nwx, 2, self.Ez, precursors.hy_left, precursors.hy_right, self.IDlookup['Ez'], -1, 1, 1, self.nthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwy, self.nwz + 1, self.nwx, 2, self.Ey, precursors.hz_left, precursors.hz_right, self.IDlookup['Ey'], 1, -1, 1, config.get_model().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwy + 1, self.nwz, self.nwx, 2, self.Ez, precursors.hy_left, precursors.hy_right, self.IDlookup['Ez'], -1, 1, 1, config.get_model().ompthreads)
 
         # front and back
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx, self.nwz + 1, self.nwy, 3, self.Ex, precursors.hz_front, precursors.hz_back, self.IDlookup['Ex'], -1, 1, 2, self.nthreads)
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx + 1, self.nwz, self.nwy, 3, self.Ez, precursors.hx_front, precursors.hx_back, self.IDlookup['Ez'], 1, -1, 2, self.nthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx, self.nwz + 1, self.nwy, 3, self.Ex, precursors.hz_front, precursors.hz_back, self.IDlookup['Ex'], -1, 1, 2, config.get_model().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx + 1, self.nwz, self.nwy, 3, self.Ez, precursors.hx_front, precursors.hx_back, self.IDlookup['Ez'], 1, -1, 2, config.get_model().ompthreads)
 
     def update_electric_os(self, main_grid):
         i_l = self.i0 - self.is_os_sep
@@ -90,16 +91,16 @@ class SubGridHSG(SubGridBase):
         # Ey = c0(Ey) + c3(dHx) - c1(dHz)
         # Ez = c0(Ez) + c1(dHy) - c2(dHx)
 
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 3, i_l, i_u, k_l, k_u + 1, j_l, j_u, self.nwy, main_grid.IDlookup['Ex'], main_grid.Ex, self.Hz, 2, 1, -1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 3, i_l, i_u + 1, k_l, k_u, j_l, j_u, self.nwy, main_grid.IDlookup['Ez'], main_grid.Ez, self.Hx, 2, -1, 1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 3, i_l, i_u, k_l, k_u + 1, j_l, j_u, self.nwy, main_grid.IDlookup['Ex'], main_grid.Ex, self.Hz, 2, 1, -1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 3, i_l, i_u + 1, k_l, k_u, j_l, j_u, self.nwy, main_grid.IDlookup['Ez'], main_grid.Ez, self.Hx, 2, -1, 1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
 
         # Left and Right
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 2, j_l, j_u, k_l, k_u + 1, i_l, i_u, self.nwx, main_grid.IDlookup['Ey'], main_grid.Ey, self.Hz, 1, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 2, j_l, j_u + 1, k_l, k_u, i_l, i_u, self.nwx, main_grid.IDlookup['Ez'], main_grid.Ez, self.Hy, 1, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 2, j_l, j_u, k_l, k_u + 1, i_l, i_u, self.nwx, main_grid.IDlookup['Ey'], main_grid.Ey, self.Hz, 1, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 2, j_l, j_u + 1, k_l, k_u, i_l, i_u, self.nwx, main_grid.IDlookup['Ez'], main_grid.Ez, self.Hy, 1, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
 
         # Bottom and Top
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 1, i_l, i_u, j_l, j_u + 1, k_l, k_u, self.nwz, main_grid.IDlookup['Ex'], main_grid.Ex, self.Hy, 3, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 1, i_l, i_u + 1, j_l, j_u, k_l, k_u, self.nwz, main_grid.IDlookup['Ey'], main_grid.Ey, self.Hx, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 1, i_l, i_u, j_l, j_u + 1, k_l, k_u, self.nwz, main_grid.IDlookup['Ex'], main_grid.Ex, self.Hy, 3, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 1, i_l, i_u + 1, j_l, j_u, k_l, k_u, self.nwz, main_grid.IDlookup['Ey'], main_grid.Ey, self.Hx, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
 
     def update_magnetic_os(self, main_grid):
 
@@ -118,16 +119,16 @@ class SubGridHSG(SubGridBase):
         # Args: sub_grid, normal, l_l, l_u, m_l, m_u, n_l, n_u, nwn, lookup_id, field, inc_field, co, sign_n, sign_f):
 
         # Front and back
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 3, i_l, i_u, k_l, k_u + 1, j_l - 1, j_u, self.nwy, main_grid.IDlookup['Hz'], main_grid.Hz, self.Ex, 2, 1, -1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 3, i_l, i_u + 1, k_l, k_u, j_l - 1, j_u, self.nwy, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ez, 2, -1, 1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 3, i_l, i_u, k_l, k_u + 1, j_l - 1, j_u, self.nwy, main_grid.IDlookup['Hz'], main_grid.Hz, self.Ex, 2, 1, -1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 3, i_l, i_u + 1, k_l, k_u, j_l - 1, j_u, self.nwy, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ez, 2, -1, 1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
 
         # Left and Right
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 2, j_l, j_u, k_l, k_u + 1, i_l - 1, i_u, self.nwx, main_grid.IDlookup['Hz'], main_grid.Hz, self.Ey, 1, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 2, j_l, j_u + 1, k_l, k_u, i_l - 1, i_u, self.nwx, main_grid.IDlookup['Hy'], main_grid.Hy, self.Ez, 1, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 2, j_l, j_u, k_l, k_u + 1, i_l - 1, i_u, self.nwx, main_grid.IDlookup['Hz'], main_grid.Hz, self.Ey, 1, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 2, j_l, j_u + 1, k_l, k_u, i_l - 1, i_u, self.nwx, main_grid.IDlookup['Hy'], main_grid.Hy, self.Ez, 1, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
 
         # bottom and top
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, i_l, i_u, j_l, j_u + 1, k_l - 1, k_u, self.nwz, main_grid.IDlookup['Hy'], main_grid.Hy, self.Ex, 3, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, i_l, i_u + 1, j_l, j_u, k_l - 1, k_u, self.nwz, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ey, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, i_l, i_u, j_l, j_u + 1, k_l - 1, k_u, self.nwz, main_grid.IDlookup['Hy'], main_grid.Hy, self.Ex, 3, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, i_l, i_u + 1, j_l, j_u, k_l - 1, k_u, self.nwz, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ey, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
 
     def __str__(self):
 
@@ -141,7 +142,7 @@ class SubGridHSG(SubGridBase):
         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'Main Grid Indices: lower left({self.i0}, {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'
 
diff --git a/gprMax/subgrids/updates.py b/gprMax/subgrids/updates.py
index fb299771..5185d07d 100644
--- a/gprMax/subgrids/updates.py
+++ b/gprMax/subgrids/updates.py
@@ -127,7 +127,7 @@ class SubgridUpdater(CPUUpdates):
         """This is the first half of the subgrid update. Takes the time step
             up to the main grid electric update.
         """
-        
+
         G = self.G
         sub_grid = self.grid
         precursors = self.precursors
diff --git a/gprMax/updates.py b/gprMax/updates.py
index 726c4783..4d27736c 100644
--- a/gprMax/updates.py
+++ b/gprMax/updates.py
@@ -101,7 +101,7 @@ class CPUUpdates:
     def update_electric_a(self):
         """Update electric field components."""
         # All materials are non-dispersive so do standard update.
-        if config.model_configs[self.grid.model_num].materials['maxpoles'] == 0:
+        if config.get_model_config().materials['maxpoles'] == 0:
             update_electric_cpu(self.grid.nx,
                                 self.grid.ny,
                                 self.grid.nz,
@@ -122,7 +122,7 @@ class CPUUpdates:
                                      self.grid.ny,
                                      self.grid.nz,
                                      config.sim_config.hostinfo['ompthreads'],
-                                     config.model_configs[self.grid.model_num].materials['maxpoles'],
+                                     config.get_model_config().materials['maxpoles'],
                                      self.grid.updatecoeffsE,
                                      self.grid.updatecoeffsdispersive,
                                      self.grid.ID,
@@ -162,12 +162,12 @@ class CPUUpdates:
             updated after the electric field has been updated by the PML and
             source updates.
         """
-        if config.model_configs[self.grid.model_num].materials['maxpoles'] != 0:
+        if config.get_model_config().materials['maxpoles'] != 0:
             self.dispersive_update_b(self.grid.nx,
                                      self.grid.ny,
                                      self.grid.nz,
                                      config.sim_config.hostinfo['ompthreads'],
-                                     config.model_configs[self.grid.model_num].materials['maxpoles'],
+                                     config.get_model_config().materials['maxpoles'],
                                      self.grid.updatecoeffsdispersive,
                                      self.grid.ID,
                                      self.grid.Tx,
@@ -183,7 +183,7 @@ class CPUUpdates:
         Returns:
             props (Props): Dispersive material properties.
         """
-        if config.model_configs[self.grid.model_num].materials['maxpoles'] > 1:
+        if config.get_model_config().materials['maxpoles'] > 1:
             poles = 'multi'
         else:
             poles = '1'
@@ -193,7 +193,7 @@ class CPUUpdates:
         else:
             type = 'double'
 
-        if config.model_configs[self.grid.model_num].materials['dispersivedtype'] == config.sim_config.dtypes['complex']:
+        if config.get_model_config().materials['dispersivedtype'] == config.sim_config.dtypes['complex']:
             dispersion = 'complex'
         else:
             dispersion = 'real'
@@ -258,7 +258,7 @@ class CUDAUpdates:
         self.drv.init()
 
         # Create device handle and context on specifc GPU device (and make it current context)
-        self.dev = self.drv.Device(config.model_configs[self.grid.model_num].cuda['gpu'].deviceID)
+        self.dev = self.drv.Device(config.get_model_config().cuda['gpu'].deviceID)
         self.ctx = self.dev.make_context()
 
         # Initialise arrays on GPU, prepare kernels, and get kernel functions
@@ -276,7 +276,7 @@ class CUDAUpdates:
         """Electric and magnetic field updates - prepare kernels, and
             get kernel functions.
         """
-        if config.model_configs[self.grid.model_num].materials['maxpoles'] > 0:
+        if config.get_model_config().materials['maxpoles'] > 0:
             kernels_fields = self.source_module(kernel_template_fields.substitute(
                                                 REAL=config.sim_config.dtypes['C_float_or_double'],
                                                 COMPLEX=config.sim_config.dtypes['C_complex'],
@@ -320,7 +320,7 @@ class CUDAUpdates:
         # - get kernel functions and initialise array on GPU
         # If there are any dispersive materials (updates are split into two
         # parts as they require present and updated electric field values).
-        if config.model_configs[self.grid.model_num].materials['maxpoles'] > 0:
+        if config.get_model_config().materials['maxpoles'] > 0:
             self.dispersive_update_a = kernels_fields.get_function("update_electric_dispersive_A")
             self.dispersive_update_b = kernels_fields.get_function("update_electric_dispersive_B")
 
@@ -329,7 +329,7 @@ class CUDAUpdates:
         self.grid.set_blocks_per_grid()
         self.grid.htod_geometry_arrays()
         self.grid.htod_field_arrays()
-        if config.model_configs[self.grid.model_num].materials['maxpoles'] > 0:
+        if config.get_model_config().materials['maxpoles'] > 0:
             self.grid.htod_dispersive_arrays()
 
     def _set_pml_kernels(self):
@@ -445,8 +445,8 @@ class CUDAUpdates:
 
         # Check if coefficient arrays will fit on constant memory of GPU
         if (self.grid.updatecoeffsE.nbytes + self.grid.updatecoeffsH.nbytes
-            > config.model_configs[self.grid.model_num].cuda['gpu'].constmem):
-            raise GeneralError(f"Too many materials in the model to fit onto constant memory of size {human_size(config.model_configs[self.grid.model_num].cuda['gpu'].constmem)} on {config.model_configs[self.grid.model_num].cuda['gpu'].deviceID} - {config.model_configs[self.grid.model_num].cuda['gpu'].name} GPU")
+            > config.get_model_config().cuda['gpu'].constmem):
+            raise GeneralError(f"Too many materials in the model to fit onto constant memory of size {human_size(config.get_model_config().cuda['gpu'].constmem)} on {config.get_model_config().cuda['gpu'].deviceID} - {config.get_model_config().cuda['gpu'].name} GPU")
 
         updatecoeffsE = kernelE.get_global('updatecoeffsE')[0]
         updatecoeffsH = kernelH.get_global('updatecoeffsH')[0]
@@ -478,7 +478,7 @@ class CUDAUpdates:
 
         for i, snap in enumerate(self.grid.snapshots):
             if snap.time == iteration + 1:
-                snapno = 0 if config.model_configs[self.grid.model_num].cuda['snapsgpu2cpu'] else i
+                snapno = 0 if config.get_model_config().cuda['snapsgpu2cpu'] else i
                 self.store_snapshot_gpu(np.int32(snapno),
                                         np.int32(snap.xs),
                                         np.int32(snap.xf),
@@ -503,7 +503,7 @@ class CUDAUpdates:
                                         self.snapHz_gpu.gpudata,
                                         block=Snapshot.tpb,
                                         grid=Snapshot.bpg)
-                if config.model_configs[self.grid.model_num].cuda['snapsgpu2cpu']:
+                if config.get_model_config().cuda['snapsgpu2cpu']:
                     dtoh_snapshot_array(self.grid.snapEx_gpu.get(),
                                         self.grid.snapEy_gpu.get(),
                                         self.grid.snapEz_gpu.get(),
@@ -553,7 +553,7 @@ class CUDAUpdates:
     def update_electric_a(self):
         """Update electric field components."""
         # All materials are non-dispersive so do standard update.
-        if config.model_configs[self.grid.model_num].materials['maxpoles'] == 0:
+        if config.get_model_config().materials['maxpoles'] == 0:
             self.update_electric_gpu(np.int32(self.grid.nx),
                                      np.int32(self.grid.ny),
                                      np.int32(self.grid.nz),
@@ -573,7 +573,7 @@ class CUDAUpdates:
             self.dispersive_update_a(np.int32(self.grid.nx),
                                      np.int32(self.grid.ny),
                                      np.int32(self.grid.nz),
-                                     np.int32(config.model_configs[self.grid.model_num].materials['maxpoles']),
+                                     np.int32(config.get_model_config().materials['maxpoles']),
                                      self.grid.updatecoeffsdispersive_gpu.gpudata,
                                      self.grid.Tx_gpu.gpudata,
                                      self.grid.Ty_gpu.gpudata,
@@ -638,11 +638,11 @@ class CUDAUpdates:
             updated after the electric field has been updated by the PML and
             source updates.
         """
-        if config.model_configs[self.grid.model_num].materials['maxpoles'] > 0:
+        if config.get_model_config().materials['maxpoles'] > 0:
             self.dispersive_update_b(np.int32(self.grid.nx),
                                      np.int32(self.grid.ny),
                                      np.int32(self.grid.nz),
-                                     np.int32(config.model_configs[self.grid.model_num].materials['maxpoles']),
+                                     np.int32(config.get_model_config().materials['maxpoles']),
                                      self.grid.updatecoeffsdispersive_gpu.gpudata,
                                      self.grid.Tx_gpu.gpudata,
                                      self.grid.Ty_gpu.gpudata,
@@ -690,7 +690,7 @@ class CUDAUpdates:
                           self.grid)
 
         # Copy data from any snapshots back to correct snapshot objects
-        if self.grid.snapshots and not config.model_configs[self.grid.model_num].cuda['snapsgpu2cpu']:
+        if self.grid.snapshots and not config.get_model_config().cuda['snapsgpu2cpu']:
             for i, snap in enumerate(self.grid.snapshots):
                 dtoh_snapshot_arra(self.snapEx_gpu.get(),
                                    self.snapEy_gpu.get(),
diff --git a/gprMax/utilities.py b/gprMax/utilities.py
index 0cdf7485..e68892cb 100644
--- a/gprMax/utilities.py
+++ b/gprMax/utilities.py
@@ -296,12 +296,12 @@ def get_host_info():
 
         # CPU information
         try:
-            cpuIDinfo = subprocess.check_output("cat /proc/cpuinfo", shell=True, stderr=subprocess.STDOUT).decode('utf-8').strip()
+            # Locale to ensure English
+            myenv = {**os.environ, 'LANG': 'en_US.utf8'}
+            cpuIDinfo = subprocess.check_output("cat /proc/cpuinfo", shell=True, stderr=subprocess.STDOUT, env=myenv).decode('utf-8').strip()
             for line in cpuIDinfo.split('\n'):
                 if re.search('model name', line):
-                    cpuID = re.sub('.*model name.*:', '', line, 1).strip()
-            # Run lscpu with English locale
-            myenv = {**os.environ, 'LANG': 'en_US.utf8'}
+                    cpuID = re.sub('.*model name.*:', '', line, 1).strip()            
             allcpuinfo = subprocess.check_output("lscpu", shell=True, stderr=subprocess.STDOUT, env=myenv).decode('utf-8').strip()
             for line in allcpuinfo.split('\n'):
                 if 'Socket(s)' in line:
@@ -393,22 +393,21 @@ def mem_check_host(mem):
         raise GeneralError(f"Memory (RAM) required ~{human_size(mem)} exceeds {human_size(config.sim_config.hostinfo['ram'], a_kilobyte_is_1024_bytes=True)} detected!\n")
 
 
-def mem_check_gpu_snaps(model_num, total_mem, snaps_mem):
+def mem_check_gpu_snaps(total_mem, snaps_mem):
     """Check if the required amount of memory (RAM) for all snapshots can fit
         on specified GPU.
 
     Args:
-        model_num (int): Model number.
         total_mem (int): Total memory required for model (bytes).
         snaps_mem (int): Memory required for all snapshots (bytes).
     """
-    if total_mem - snaps_mem > config.model_configs[model_num].cuda['gpu'].totalmem:
-        raise GeneralError(f"Memory (RAM) required ~{human_size(total_mem)} exceeds {human_size(config.model_configs[model_num].cuda['gpu'].totalmem, a_kilobyte_is_1024_bytes=True)} detected on specified {config.model_configs[model_num].cuda['gpu'].deviceID} - {config.model_configs[model_num].cuda['gpu'].name} GPU!\n")
+    if total_mem - snaps_mem > config.get_model_config().cuda['gpu'].totalmem:
+        raise GeneralError(f"Memory (RAM) required ~{human_size(total_mem)} exceeds {human_size(config.get_model_config().cuda['gpu'].totalmem, a_kilobyte_is_1024_bytes=True)} detected on specified {config.get_model_config().cuda['gpu'].deviceID} - {config.get_model_config().cuda['gpu'].name} GPU!\n")
 
     # If the required memory without the snapshots will fit on the GPU then
     # transfer and store snaphots on host
-    if snaps_mem != 0 and total_mem - snaps_mem < config.model_configs[model_num].cuda['gpu'].totalmem:
-        config.model_configs[model_num].cuda['snapsgpu2cpu'] = True
+    if snaps_mem != 0 and total_mem - snaps_mem < config.get_model_config().cuda['gpu'].totalmem:
+        config.get_model_config().cuda['snapsgpu2cpu'] = True
 
 
 def mem_check_all(grids):
@@ -427,36 +426,36 @@ def mem_check_all(grids):
     total_snaps_mem = 0
 
     for grid in grids:
-        config.model_configs[grid.model_num].mem_use += grid.mem_est_basic()
+        config.get_model_config().mem_use += grid.mem_est_basic()
 
         # Set datatype for dispersive arrays if there are any dispersive materials.
-        if config.model_configs[grid.model_num].materials['maxpoles'] != 0:
+        if config.get_model_config().materials['maxpoles'] != 0:
             drudelorentz = any([m for m in grid.materials if 'drude' in m.type or 'lorentz' in m.type])
             if drudelorentz:
-                config.model_configs[grid.model_num].materials['dispersivedtype'] = config.sim_config.dtypes['complex']
-                config.model_configs[grid.model_num].materials['dispersiveCdtype'] = config.sim_config.dtypes['C_complex']
+                config.get_model_config().materials['dispersivedtype'] = config.sim_config.dtypes['complex']
+                config.get_model_config().materials['dispersiveCdtype'] = config.sim_config.dtypes['C_complex']
             else:
-                config.model_configs[grid.model_num].materials['dispersivedtype'] = config.sim_config.dtypes['float_or_double']
-                config.model_configs[grid.model_num].materials['dispersiveCdtype'] = config.sim_config.dtypes['C_float_or_double']
+                config.get_model_config().materials['dispersivedtype'] = config.sim_config.dtypes['float_or_double']
+                config.get_model_config().materials['dispersiveCdtype'] = config.sim_config.dtypes['C_float_or_double']
 
             # Update estimated memory (RAM) usage
-            config.model_configs[grid.model_num].mem_use += grid.mem_est_dispersive()
+            config.get_model_config().mem_use += grid.mem_est_dispersive()
 
         # Calculate snapshot memory
         if grid.snapshots:
             for snap in grid.snapshots:
                 # 2 x required to account for electric and magnetic fields
-                config.model_configs[grid.model_num].mem_use += int(2 * snap.datasizefield)
+                config.get_model_config().mem_use += int(2 * snap.datasizefield)
                 total_snaps_mem += int(2 * snap.datasizefield)
 
-        total_mem += config.model_configs[grid.model_num].mem_use
+        total_mem += config.get_model_config().mem_use
 
     # Check if there is sufficient memory on host
     mem_check_host(total_mem)
 
     # Check if there is sufficient memory for any snapshots on GPU
     if config.sim_config.general['cuda']:
-        mem_check_gpu_snaps(grid.model_num, total_mem, total_snaps_mem)
+        mem_check_gpu_snaps(total_mem, total_snaps_mem)
 
     return total_mem
 
diff --git a/setup.py b/setup.py
index 9e844fda..1fbce0a1 100644
--- a/setup.py
+++ b/setup.py
@@ -123,8 +123,8 @@ with open('README.rst','r') as fd:
     long_description = fd.read()
 
 # Python version
-if sys.version_info[:2] < (3, 6):
-    sys.exit('\nExited: Requires Python 3.6 or newer!\n')
+if sys.version_info[:2] < (3, 7):
+    sys.exit('\nExited: Requires Python 3.7 or newer!\n')
 
 # Process 'build' command line argument
 if 'build' in sys.argv: