diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 00b38de1..04073d30 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -6,6 +6,20 @@ recent versions
 """""""""""""""
 
 
+*latest (dev, will become v0.14.0)*
+-----------------------------------
+
+- ``models.Model``: Removed deprecated parameters ``res_{x;y;z}``.
+
+- ``io.save``:
+
+  - Removed deprecated parameter ``backend`` in ``save``.
+  - File extension has to be provided in ``save``, no default any longer.
+
+- ``meshes``: Deprecated the old meshing routines ``get_stretched_h``,
+  ``get_domain``, ``get_hx``; will be removed in the future.
+
+
 *v0.13.0* : CLI
 ---------------
 
diff --git a/emg3d/io.py b/emg3d/io.py
index 5ec4ce0f..b0128f2d 100644
--- a/emg3d/io.py
+++ b/emg3d/io.py
@@ -19,7 +19,6 @@
 
 import os
 import json
-import warnings
 from datetime import datetime
 
 import numpy as np
@@ -46,7 +45,7 @@
 }
 
 
-def save(fname, backend=None, compression="gzip", **kwargs):
+def save(fname, **kwargs):
     """Save surveys, meshes, models, fields, and more to disk.
 
     Serialize and save data to disk in different formats (see parameter
@@ -66,18 +65,14 @@ def save(fname, backend=None, compression="gzip", **kwargs):
         File name inclusive ending, which defines the used data format.
         Implemented are currently:
 
-        - `.h5` (default): Uses `h5py` to store inputs to a hierarchical,
-          compressed binary hdf5 file. Recommended file format, but requires
-          the module `h5py`. Default format if ending is not provided or not
-          recognized.
+        - `.h5`: Uses `h5py` to store inputs to a hierarchical, compressed
+          binary hdf5 file. Recommended file format, but requires the module
+          `h5py`.
         - `.npz`: Uses `numpy` to store inputs to a flat, compressed binary
-          file. Default format if `h5py` is not installed.
+          file.
         - `.json`: Uses `json` to store inputs to a hierarchical, plain text
           file.
 
-    backend : deprecated
-        Set the appropriate file-ending in `fname` instead.
-
     compression : int or str, optional
         Passed through to h5py, default is 'gzip'.
 
@@ -104,6 +99,7 @@ def save(fname, backend=None, compression="gzip", **kwargs):
 
     """
     # Get and remove optional kwargs.
+    compression = kwargs.pop('compression', 'gzip')
     json_indent = kwargs.pop('json_indent', 2)
     collect_classes = kwargs.pop('collect_classes', False)
     verb = kwargs.pop('verb', 1)
@@ -120,37 +116,8 @@ def save(fname, backend=None, compression="gzip", **kwargs):
     # sorted TensorMesh, Field, and Model instances.
     data = _dict_serialize(kwargs, collect_classes=collect_classes)
 
-    # Deprecated backends.
-    if backend is not None:
-        mesg = ("\n    The use of `backend` is deprecated and will be removed."
-                "\n    Use the file-endings [`.h5`, `.npz`, `.json`] instead.")
-        warnings.warn(mesg, DeprecationWarning)
-
-        if backend == 'numpy':
-            backend = 'npz'
-        elif backend == 'h5py':
-            backend = 'h5'
-        else:
-            backend = backend
-
-        # Add file-ending if necessary.
-        if not full_path.endswith('.'+backend):
-            full_path += '.'+backend
-
-    else:
-
-        # Get backend from `fname`.
-        if full_path.split('.')[-1] in ['npz', 'h5', 'json']:
-            backend = full_path.split('.')[-1]
-        else:  # Fallback to default.
-            if isinstance(h5py, str):
-                backend = 'npz'
-            else:
-                backend = 'h5'
-            full_path += '.'+backend
-
-    # Save data depending on the backend.
-    if backend == "npz":
+    # Save data depending on the extension.
+    if full_path.endswith('.npz'):
 
         # Convert hierarchical dict to a flat dict.
         data = _dict_flatten(data)
@@ -158,7 +125,7 @@ def save(fname, backend=None, compression="gzip", **kwargs):
         # Store flattened data.
         np.savez_compressed(full_path, **data)
 
-    elif backend == "h5":
+    elif full_path.endswith('.h5'):
 
         # Check if h5py is installed.
         if isinstance(h5py, str):
@@ -168,7 +135,7 @@ def save(fname, backend=None, compression="gzip", **kwargs):
         with h5py.File(full_path, "w") as h5file:
             _hdf5_add_to(data, h5file, compression)
 
-    elif backend == "json":
+    elif full_path.endswith('.json'):
 
         # Move arrays to lists and decompose complex data.
         data = _dict_dearray_decomp(data)
@@ -178,7 +145,8 @@ def save(fname, backend=None, compression="gzip", **kwargs):
             json.dump(data, f, indent=json_indent)
 
     else:
-        raise ValueError(f"Unknown backend '{backend}'.")
+        ext = full_path.split('.')[-1]
+        raise ValueError(f"Unknown extension '.{ext}'.")
 
     # Print file info.
     if verb > 0:
@@ -197,8 +165,7 @@ def load(fname, **kwargs):
     Parameters
     ----------
     fname : str
-        File name including extension. Used backend depends on the file
-        extensions:
+        File name including extension. Possibilities:
 
         - '.npz': numpy-binary
         - '.h5': h5py-binary (needs `h5py`)
@@ -696,7 +663,7 @@ def _hdf5_get_from(h5file):
 def _compare_dicts(dict1, dict2, verb=False, **kwargs):
     """Return True if the two dicts `dict1` and `dict2` are the same.
 
-    Private method, not foolproof. Useful for developing new backends.
+    Private method, not foolproof. Useful for developing new extensions.
 
     If `verb=True`, it prints it key starting with the following legend:
 
diff --git a/emg3d/meshes.py b/emg3d/meshes.py
index 5483441a..0caad753 100644
--- a/emg3d/meshes.py
+++ b/emg3d/meshes.py
@@ -17,6 +17,7 @@
 # License for the specific language governing permissions and limitations under
 # the License.
 
+import warnings
 from copy import deepcopy
 
 import numpy as np
@@ -31,6 +32,11 @@ class dTensorMesh:
 __all__ = ['TensorMesh', 'get_hx_h0', 'get_cell_numbers', 'get_stretched_h',
            'get_domain', 'get_hx']
 
+MESHWARNING = (
+    "\n    `get_stretched_h`, `get_domain`, and `get_hx` are"
+    "\n    deprecated and will be removed. Use `get_hx_h0` instead."
+)
+
 
 class _TensorMesh:
     """Minimal TensorMesh for internal multigrid computation.
@@ -599,7 +605,7 @@ def get_stretched_h(min_width, domain, nx, x0=0, x1=None, resp_domain=False):
     Returns `nx` cell widths within `domain`, where the minimum cell width is
     `min_width`. The cells are not stretched within `x0` and `x1`, and outside
     uses a power-law stretching. The actual stretching factor and the number of
-    cells left and right of `x0` and `x1` are find in a minimization process.
+    cells left and right of `x0` and `x1` are found in a minimization process.
 
     The domain is not completely respected. The starting point of the domain
     is, but the endpoint of the domain might slightly shift (this is more
@@ -656,6 +662,7 @@ def get_stretched_h(min_width, domain, nx, x0=0, x1=None, resp_domain=False):
         Cell widths of mesh.
 
     """
+    warnings.warn(MESHWARNING, DeprecationWarning)
 
     # Cast to arrays
     domain = np.array(domain, dtype=np.float64)
@@ -806,6 +813,7 @@ def get_domain(x0=0, freq=1, res=0.3, limits=None, min_width=None,
         Start- and end-points of computation domain.
 
     """
+    warnings.warn(MESHWARNING, DeprecationWarning)
 
     # Set fact_pos to fact_neg if not provided.
     if fact_pos is None:
@@ -835,11 +843,11 @@ def get_domain(x0=0, freq=1, res=0.3, limits=None, min_width=None,
 def get_hx(alpha, domain, nx, x0, resp_domain=True):
     r"""Return cell widths for given input.
 
-    Find the number of cells left and right of `x0`, `nl` and `nr`
-    respectively, for the provided alpha. For this, we solve
+    Find the number of cells left (``nl``) and right (``nr``) of the center
+    ``x0`` for the provided alpha. For this, we solve
 
     .. math::   \frac{x_\text{max}-x_0}{x_0-x_\text{min}} =
-                \frac{a^{nr}-1}{a^{nl}-1}
+                \frac{a^\text{nr}-1}{a^\text{nl}-1}
 
     where :math:`a = 1+\alpha`.
 
@@ -851,17 +859,17 @@ def get_hx(alpha, domain, nx, x0, resp_domain=True):
         Stretching factor `a` is given by ``a=1+alpha``.
 
     domain : list
-        [start, end] of model domain.
+        ``[x_min, x_max]`` of model domain.
 
     nx : int
         Number of cells.
 
     x0 : float
-        Center of the grid. `x0` is restricted to `domain`.
+        Center of the grid. ``x0`` is restricted to ``domain``.
 
     resp_domain : bool
-        If False (default), then the domain-end might shift slightly to assure
-        that the same stretching factor is applied throughout. If set to True,
+        If False, then the domain-end might shift slightly to assure that the
+        same stretching factor is applied throughout. If set to True (default),
         however, the domain is respected absolutely. This will introduce one
         stretch-factor which is different from the other stretch factors, to
         accommodate the restriction. This one-off factor is between the left-
@@ -871,9 +879,12 @@ def get_hx(alpha, domain, nx, x0, resp_domain=True):
     Returns
     -------
     hx : ndarray
-        Cell widths of mesh.
+        Cell widths of mesh. All points are given by
+        ``np.r_[xmin, xmin+np.cumsum(hx)]``
 
     """
+    warnings.warn(MESHWARNING, DeprecationWarning)
+
     if alpha <= 0.:  # If alpha <= 0: equal spacing (no stretching at all)
         hx = np.ones(nx)*np.diff(np.squeeze(domain))/nx
 
@@ -882,7 +893,7 @@ def get_hx(alpha, domain, nx, x0, resp_domain=True):
 
         # Get hx depending if x0 is on the domain boundary or not.
         if np.isclose(x0, domain[0]) or np.isclose(x0, domain[1]):
-            # Get al a's
+            # Get all a's
             alr = np.diff(domain)*alpha/(a**nx-1)*a**np.arange(nx)
             if x0 == domain[1]:
                 alr = alr[::-1]
diff --git a/emg3d/models.py b/emg3d/models.py
index 467b5fea..92659357 100644
--- a/emg3d/models.py
+++ b/emg3d/models.py
@@ -17,7 +17,6 @@
 # License for the specific language governing permissions and limitations under
 # the License.
 
-import warnings
 from copy import deepcopy
 
 import numpy as np
@@ -89,29 +88,10 @@ class Model:
 
     """
 
-    _res_warning = (
-        "\n    The keywords `res_{x;y;z}` are deprecated and will be removed."
-        "\n    Use the keywords `property_{x;y;z}` instead."
-    )
-
     def __init__(self, grid, property_x=1., property_y=None, property_z=None,
                  mu_r=None, epsilon_r=None, mapping='Resistivity', **kwargs):
         """Initiate a new model."""
 
-        # Check kwargs; purely for backwards compatibility.
-        if any([key.startswith('res') for key in kwargs.keys()]):
-            warnings.warn(self._res_warning, DeprecationWarning)
-
-            res_x = kwargs.pop('res_x', None)
-            res_y = kwargs.pop('res_y', None)
-            res_z = kwargs.pop('res_z', None)
-
-            property_x = property_x if res_x is None else res_x
-            property_y = property_y if res_y is None else res_y
-            property_z = property_z if res_z is None else res_z
-
-            mapping = 'Resistivity'
-
         # Ensure no kwargs left.
         if kwargs:
             raise TypeError(f"Unexpected **kwargs: {list(kwargs.keys())}")
@@ -295,15 +275,6 @@ def from_dict(cls, inp):
 
         """
         try:
-            # Check `res`; purely for backwards compatibility.
-            if any([key.startswith('res') for key in inp]):
-                warnings.warn(cls._res_warning, DeprecationWarning)
-
-                inp['property_x'] = inp.pop('res_x', None)
-                inp['property_y'] = inp.pop('res_y', None)
-                inp['property_z'] = inp.pop('res_z', None)
-                inp['mapping'] = 'Resistivity'
-
             return cls(grid=inp['vnC'], property_x=inp['property_x'],
                        property_y=inp['property_y'],
                        property_z=inp['property_z'], mu_r=inp['mu_r'],
@@ -367,25 +338,6 @@ def property_z(self, property_z):
         self._property_z = self._check_parameter(
                 property_z, 'property_z', True)
 
-    # Backwards compatibility for deprecated attributes.
-    @property
-    def res_x(self):
-        warnings.warn(self._res_warning, DeprecationWarning)
-        fmap = maps.MapResistivity()
-        return fmap.forward(self.map.backward(self.property_x))
-
-    @property
-    def res_y(self):
-        warnings.warn(self._res_warning, DeprecationWarning)
-        fmap = maps.MapResistivity()
-        return fmap.forward(self.map.backward(self.property_y))
-
-    @property
-    def res_z(self):
-        warnings.warn(self._res_warning, DeprecationWarning)
-        fmap = maps.MapResistivity()
-        return fmap.forward(self.map.backward(self.property_z))
-
     # MAGNETIC PERMEABILITIES
     @property
     def mu_r(self):
diff --git a/tests/create_data/regression.py b/tests/create_data/regression.py
index ea3514c0..0a527ade 100644
--- a/tests/create_data/regression.py
+++ b/tests/create_data/regression.py
@@ -8,14 +8,14 @@
 # # # # # # # # # # 1. Homogeneous VTI fullspace # # # # # # # # # #
 
 freq = 1.
-hx_min, xdomain = meshes.get_domain(x0=0, freq=freq)
-hy_min, ydomain = meshes.get_domain(x0=0, freq=freq)
-hz_min, zdomain = meshes.get_domain(x0=250, freq=freq)
-nx = 2**3
-hx = meshes.get_stretched_h(hx_min, xdomain, nx, 0)
-hy = meshes.get_stretched_h(hy_min, ydomain, nx, 0)
-hz = meshes.get_stretched_h(hz_min, zdomain, nx, 250)
-input_grid = {'h': [hx, hy, hz], 'x0': (xdomain[0], ydomain[0], zdomain[0])}
+hxy, hxy0 = meshes.get_hx_h0(
+        freq=freq, res=2, domain=[-50, 50], max_domain=1500,
+        possible_nx=meshes.get_cell_numbers(100, 2, 1), verb=0)
+hz, hz0 = meshes.get_hx_h0(
+        freq=freq, res=3.3, domain=[200, 300], max_domain=1500,
+        possible_nx=meshes.get_cell_numbers(100, 2, 1), verb=0)
+
+input_grid = {'h': [hxy, hxy, hz], 'x0': (hxy0, hxy0, hz0)}
 grid = meshes.TensorMesh(**input_grid)
 
 input_model = {
@@ -67,14 +67,15 @@
 src = [50., 110., 250., 25, 15]
 freq = 0.375
 
-hx_min, xdomain = meshes.get_domain(x0=0, freq=.1)
-hy_min, ydomain = meshes.get_domain(x0=0, freq=.1)
-hz_min, zdomain = meshes.get_domain(x0=250, freq=.1)
-hx = meshes.get_stretched_h(hx_min, xdomain, 8, 0)
-hy = meshes.get_stretched_h(hy_min, ydomain, 4, 0)
-hz = meshes.get_stretched_h(hz_min, zdomain, 16, 250)
-grid = meshes.TensorMesh([hx, hy, hz], x0=(xdomain[0], ydomain[0], zdomain[0]))
+hxy, hxy0 = meshes.get_hx_h0(
+        freq=freq, res=50, domain=[src[0], src[1]], max_domain=5000,
+        possible_nx=meshes.get_cell_numbers(100, 2, 1), verb=0)
+hz, hz0 = meshes.get_hx_h0(
+        freq=freq, res=3.3, domain=[src[2]-20, src[2]+20], max_domain=5000,
+        possible_nx=meshes.get_cell_numbers(100, 2, 1), verb=0)
 
+input_grid = {'h': [hxy, hxy, hz], 'x0': (hxy0, hxy0, hz0)}
+grid = meshes.TensorMesh(**input_grid)
 
 # Initialize model
 # Create a model with random resistivities between [0, 50)
@@ -99,10 +100,9 @@
 clevel = 10  # Way to high
 
 efield = solver.solve(
-        grid, model, sfield,
-        semicoarsening=semicoarsening, linerelaxation=linerelaxation,
-        verb=verb, tol=tol, maxit=maxit, nu_init=nu_init, nu_pre=nu_pre,
-        nu_coarse=nu_coarse, nu_post=nu_post, clevel=clevel)
+        grid, model, sfield, semicoarsening=semicoarsening,
+        linerelaxation=linerelaxation, tol=tol, maxit=maxit, nu_init=nu_init,
+        nu_pre=nu_pre, nu_coarse=nu_coarse, nu_post=nu_post, clevel=clevel)
 
 hfield = fields.get_h_field(grid, model, efield)
 
@@ -114,7 +114,7 @@
     'inp': {
         'semicoarsening': semicoarsening,
         'linerelaxation': linerelaxation,
-        'verb': verb,
+        'verb': 4,
         'tol': tol,
         'maxit': maxit,
         'nu_init': nu_init,
@@ -152,14 +152,14 @@
 # # # # # # # # # # 4. Homogeneous VTI fullspace LAPLACE # # # # # # # # # #
 
 freq = -2*np.pi
-hx_min, xdomain = meshes.get_domain(x0=0, freq=freq)
-hy_min, ydomain = meshes.get_domain(x0=0, freq=freq)
-hz_min, zdomain = meshes.get_domain(x0=250, freq=freq)
-nx = 2**3
-hx = meshes.get_stretched_h(hx_min, xdomain, nx, 0)
-hy = meshes.get_stretched_h(hy_min, ydomain, nx, 0)
-hz = meshes.get_stretched_h(hz_min, zdomain, nx, 250)
-input_grid_l = {'h': [hx, hy, hz], 'x0': (xdomain[0], ydomain[0], zdomain[0])}
+hxy, hxy0 = meshes.get_hx_h0(
+        freq=freq, res=2, domain=[-1, 1], max_domain=1500,
+        possible_nx=meshes.get_cell_numbers(100, 2, 1), verb=0)
+hz, hz0 = meshes.get_hx_h0(
+        freq=freq, res=3.3, domain=[240, 260], max_domain=1500,
+        possible_nx=meshes.get_cell_numbers(100, 2, 1), verb=0)
+
+input_grid_l = {'h': [hxy, hxy, hz], 'x0': (hxy0, hxy0, hz0)}
 grid_l = meshes.TensorMesh(**input_grid_l)
 
 input_model_l = {
diff --git a/tests/data/regression.npz b/tests/data/regression.npz
index 6770da64..cba67e19 100644
Binary files a/tests/data/regression.npz and b/tests/data/regression.npz differ
diff --git a/tests/test_core.py b/tests/test_core.py
index efb15ab1..0e794be3 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -29,9 +29,9 @@ def test_amat_x(njit):
     x = np.arange(1, grid.nCx+1)*2
     y = 1/np.arange(1, grid.nCy+1)
     z = np.arange(1, grid.nCz+1)[::-1]/10
-    res_x = np.outer(np.outer(x, y), z).ravel()
+    property_x = np.outer(np.outer(x, y), z).ravel()
     freq = 0.319
-    model = models.Model(grid, res_x, 0.8*res_x, 2*res_x)
+    model = models.Model(grid, property_x, 0.8*property_x, 2*property_x)
 
     # Create a source field
     sfield = fields.get_source_field(grid=grid, src=src, freq=freq)
@@ -119,11 +119,11 @@ def test_gauss_seidel(njit):
             [hx, hy, hz], np.array([-hx.sum()/2, -hy.sum()/2, -hz.sum()/2]))
 
         # Initialize model with some resistivities.
-        res_x = np.arange(grid.nC)+1
-        res_y = 0.5*np.arange(grid.nC)+1
-        res_z = 2*np.arange(grid.nC)+1
+        property_x = np.arange(grid.nC)+1
+        property_y = 0.5*np.arange(grid.nC)+1
+        property_z = 2*np.arange(grid.nC)+1
 
-        model = models.Model(grid, res_x, res_y, res_z)
+        model = models.Model(grid, property_x, property_y, property_z)
 
         # Initialize source field.
         sfield = fields.get_source_field(grid, src, freq)
diff --git a/tests/test_fields.py b/tests/test_fields.py
index 90a39660..af434522 100644
--- a/tests/test_fields.py
+++ b/tests/test_fields.py
@@ -301,7 +301,7 @@ def test_get_h_field():
     # Mainly regression tests, not ideal.
 
     # Check it does still the same (pure regression).
-    dat = REGRES['Data']['reg_2']
+    dat = REGRES['reg_2']
     grid = dat['grid']
     model = dat['model']
     efield = dat['result']
@@ -313,7 +313,7 @@ def test_get_h_field():
     assert hout.is_electric is False
 
     # Add some mu_r - Just 1, to trigger, and compare.
-    dat = REGRES['Data']['res']
+    dat = REGRES['res']
     grid = dat['grid']
     efield = dat['Fresult']
     model1 = models.Model(**dat['input_model'])
@@ -396,7 +396,7 @@ def test_get_receiver():
     # Check it works with model parameters.
     model = models.Model(grid, np.ones(grid.vnC))
     out10 = fields.get_receiver(
-            grid, model.res_x, (-10, -10, -10), 'linear', True)
+            grid, model.property_x, (-10, -10, -10), 'linear', True)
     assert_allclose(out10, 1.)
     assert out10.__class__ != utils.EMArray
 
diff --git a/tests/test_io.py b/tests/test_io.py
index 5de40e29..3cc8c9cc 100644
--- a/tests/test_io.py
+++ b/tests/test_io.py
@@ -52,11 +52,11 @@ class TensorMesh:
     field.ensure_pec
 
     # Some model.
-    res_x = create_dummy(*grid.vnC, False)
-    res_y = res_x/2.0
-    res_z = res_x*1.4
-    mu_r = res_x*1.11
-    model = models.Model(grid, res_x, res_y, res_z, mu_r=mu_r)
+    property_x = create_dummy(*grid.vnC, False)
+    property_y = property_x/2.0
+    property_z = property_x*1.4
+    mu_r = property_x*1.11
+    model = models.Model(grid, property_x, property_y, property_z, mu_r=mu_r)
 
     # Save it.
     io.save(tmpdir+'/test.npz', emg3d=grid, discretize=grid2, model=model,
@@ -97,25 +97,15 @@ class TensorMesh:
     with pytest.raises(TypeError, match="Unexpected "):
         io.load('ttt.npz', stupidkeyword='a')
 
-    # Unknown backend/extension.
-    with pytest.raises(ValueError, match="Unknown backend 'what"):
-        io.save(tmpdir+'/testwrongbackend', something=1, backend='what?')
-    io.save(tmpdir+'/testwrongbackend.abc', something=1)
+    # Unknown extension.
     with pytest.raises(ValueError, match="Unknown extension '.abc'"):
-        io.load(tmpdir+'/testwrongbackend.abc')
-    if h5py:
-        io.load(tmpdir+'/testwrongbackend.abc.h5')
-    else:
-        io.load(tmpdir+'/testwrongbackend.abc.npz')
-
-    # Ensure deprecated backend/extension still work.
-    if h5py:
-        io.save(tmpdir+'/ttt', backend='h5py')
-    io.save(tmpdir+'/ttt', backend='numpy')
+        io.save(tmpdir+'/testwrongextension.abc', something=1)
+    with pytest.raises(ValueError, match="Unknown extension '.abc'"):
+        io.load(tmpdir+'/testwrongextension.abc')
 
     # Test h5py.
     if h5py:
-        io.save(tmpdir+'/test', emg3d=grid, discretize=grid2,
+        io.save(tmpdir+'/test.h5', emg3d=grid, discretize=grid2,
                 a=1.0, b=1+1j, c=True,
                 model=model, field=field, what={'f': field.fx},
                 collect_classes=True)
@@ -131,18 +121,15 @@ class TensorMesh:
 
         assert io._compare_dicts(out_h5, out_npz) is True
     else:
-        with pytest.raises(ImportError):
-            # Ensure deprecated backend/extension still work.
-            io.save(tmpdir+'/ttt', backend='h5py')
         with pytest.raises(ImportError):
             io.save(tmpdir+'/test.h5', grid=grid)
         with pytest.raises(ImportError):
             io.load(str(tmpdir+'/test-h5.h5'))
 
     # Test json.
-    io.save(tmpdir+'/test', emg3d=grid, discretize=grid2,
+    io.save(tmpdir+'/test.json', emg3d=grid, discretize=grid2,
             a=1.0, b=1+1j,
-            model=model, field=field, what={'f': field.fx}, backend='json',
+            model=model, field=field, what={'f': field.fx},
             collect_classes=True)
     out_json = io.load(str(tmpdir+'/test.json'))
     assert out_json['Model']['model'] == model
diff --git a/tests/test_meshes.py b/tests/test_meshes.py
index a8e3a9c9..b1e4951d 100644
--- a/tests/test_meshes.py
+++ b/tests/test_meshes.py
@@ -223,7 +223,7 @@ def test_get_hx():
 
 def test_TensorMesh():
     # Load mesh created with discretize.TensorMesh.
-    grid = REGRES['Data']['grid']
+    grid = REGRES['grid']
 
     # Use this grid instance to create emg3d equivalent.
     emg3dgrid = meshes.TensorMesh(
diff --git a/tests/test_models.py b/tests/test_models.py
index ad1f3a28..aff460e4 100644
--- a/tests/test_models.py
+++ b/tests/test_models.py
@@ -29,10 +29,10 @@ def test_regression(self, capsys):
                 [np.array([2, 2]), np.array([3, 4]), np.array([0.5, 2])],
                 np.zeros(3))
 
-        res_x = create_dummy(*grid.vnC, False)
-        res_y = res_x/2.0
-        res_z = res_x*1.4
-        mu_r = res_x*1.11
+        property_x = create_dummy(*grid.vnC, False)
+        property_y = property_x/2.0
+        property_z = property_x*1.4
+        mu_r = property_x*1.11
 
         _, _ = capsys.readouterr()  # Clean-up
         # Using defaults; check backwards compatibility for freq.
@@ -112,19 +112,21 @@ def test_regression(self, capsys):
 
         # Check with all inputs
         gridvol = grid.vol.reshape(grid.vnC, order='F')
-        model3 = models.Model(grid, res_x, res_y, res_z, mu_r=mu_r)
+        model3 = models.Model(
+            grid, property_x, property_y, property_z, mu_r=mu_r)
         vmodel3 = models.VolumeModel(grid, model3, sfield)
-        assert_allclose(model3.res_x, model3.res_y*2)
-        assert_allclose(model3.res_x.shape, grid.vnC)
-        assert_allclose(model3.res_x, model3.res_z/1.4)
+        assert_allclose(model3.property_x, model3.property_y*2)
+        assert_allclose(model3.property_x.shape, grid.vnC)
+        assert_allclose(model3.property_x, model3.property_z/1.4)
         assert_allclose(gridvol/mu_r, vmodel3.zeta)
         # Check with all inputs
-        model3b = models.Model(grid, res_x.ravel('F'), res_y.ravel('F'),
-                               res_z.ravel('F'), mu_r=mu_r.ravel('F'))
+        model3b = models.Model(
+            grid, property_x.ravel('F'), property_y.ravel('F'),
+            property_z.ravel('F'), mu_r=mu_r.ravel('F'))
         vmodel3b = models.VolumeModel(grid, model3b, sfield)
-        assert_allclose(model3b.res_x, model3b.res_y*2)
-        assert_allclose(model3b.res_x.shape, grid.vnC)
-        assert_allclose(model3b.res_x, model3b.res_z/1.4)
+        assert_allclose(model3b.property_x, model3b.property_y*2)
+        assert_allclose(model3b.property_x.shape, grid.vnC)
+        assert_allclose(model3b.property_x, model3b.property_z/1.4)
         assert_allclose(gridvol/mu_r, vmodel3b.zeta)
 
         # Check setters vnC
@@ -135,16 +137,16 @@ def test_regression(self, capsys):
         model3.property_z = tres*4.0
         model3.mu_r = tres*5.0
         model3.epsilon_r = tres*6.0
-        assert_allclose(tres*2., model3.res_x)
-        assert_allclose(tres*3., model3.res_y)
-        assert_allclose(tres*4., model3.res_z)
+        assert_allclose(tres*2., model3.property_x)
+        assert_allclose(tres*3., model3.property_y)
+        assert_allclose(tres*4., model3.property_z)
         assert_allclose(tres*6., model3.epsilon_r)
 
         # Check eta
         iomep = sfield.sval*models.epsilon_0
-        eta_x = sfield.smu0*(1./model3.res_x + iomep)*gridvol
-        eta_y = sfield.smu0*(1./model3.res_y + iomep)*gridvol
-        eta_z = sfield.smu0*(1./model3.res_z + iomep)*gridvol
+        eta_x = sfield.smu0*(1./model3.property_x + iomep)*gridvol
+        eta_y = sfield.smu0*(1./model3.property_y + iomep)*gridvol
+        eta_z = sfield.smu0*(1./model3.property_z + iomep)*gridvol
         vmodel3 = models.VolumeModel(grid, model3, sfield)
         assert_allclose(vmodel3.eta_x, eta_x)
         assert_allclose(vmodel3.eta_y, eta_y)
@@ -158,16 +160,16 @@ def test_regression(self, capsys):
 
         # Check a couple of out-of-range failures
         with pytest.raises(ValueError, match='`property_x` must be all'):
-            _ = models.Model(grid, res_x=res_x*0)
-        Model = models.Model(grid, res_x=res_x)
+            _ = models.Model(grid, property_x=property_x*0)
+        Model = models.Model(grid, property_x=property_x)
         with pytest.raises(ValueError, match='`property_x` must be all'):
-            Model._check_parameter(res_x*0, 'property_x')
+            Model._check_parameter(property_x*0, 'property_x')
         with pytest.raises(ValueError, match='`property_x` must be all'):
             Model._check_parameter(-1.0, 'property_x')
         with pytest.raises(ValueError, match='`property_y` must be all'):
-            _ = models.Model(grid, res_y=np.inf)
+            _ = models.Model(grid, property_y=np.inf)
         with pytest.raises(ValueError, match='`property_z` must be all'):
-            _ = models.Model(grid, res_z=res_z*np.inf)
+            _ = models.Model(grid, property_z=property_z*np.inf)
         with pytest.raises(ValueError, match='`mu_r` must be all'):
             _ = models.Model(grid, mu_r=-1)
 
@@ -241,24 +243,6 @@ def test_equal_mapping(self):
 
         assert check is False
 
-    def test_old_dict(self):
-        grid = meshes.TensorMesh(
-                [np.array([2, 2]), np.array([3, 4]), np.array([0.5, 2])],
-                np.zeros(3))
-
-        model1 = models.Model(grid, 1., 2., 3.)
-        mydict = model1.to_dict()
-        mydict['res_x'] = mydict['property_x']
-        mydict['res_y'] = mydict['property_y']
-        mydict['res_z'] = mydict['property_z']
-        del mydict['property_x']
-        del mydict['property_y']
-        del mydict['property_z']
-
-        model2 = models.Model.from_dict(mydict)
-
-        assert model1 == model2
-
     def test_negative_values(self):
         # Create some dummy data
         grid = meshes.TensorMesh(
@@ -290,8 +274,8 @@ class TestModelOperators:
     model_int = models.Model(mesh_base, 1.)
     model_1_a = models.Model(mesh_base, 1., 2.)
     model_1_b = models.Model(mesh_base, 2., 4.)
-    model_2_a = models.Model(mesh_base, 1., res_z=3.)
-    model_2_b = models.Model(mesh_base, 2., res_z=6.)
+    model_2_a = models.Model(mesh_base, 1., property_z=3.)
+    model_2_b = models.Model(mesh_base, 2., property_z=6.)
     model_3_a = models.Model(mesh_base, 1., 2., 3.)
     model_3_b = models.Model(mesh_base, 2., 4., 6.)
     model_mu_a = models.Model(mesh_base, 1., mu_r=1.)
@@ -333,19 +317,19 @@ def test_add(self):
 
         # All different cases
         a = self.model_1_a + self.model_1_a
-        assert a.res_x == self.model_1_b.res_x
-        assert a.res_y == self.model_1_b.res_y
-        assert a.res_z.base is a.res_x.base
+        assert a.property_x == self.model_1_b.property_x
+        assert a.property_y == self.model_1_b.property_y
+        assert a.property_z.base is a.property_x.base
 
         a = self.model_2_a + self.model_2_a
-        assert a.res_x == self.model_2_b.res_x
-        assert a.res_y.base is a.res_x.base
-        assert a.res_z == self.model_2_b.res_z
+        assert a.property_x == self.model_2_b.property_x
+        assert a.property_y.base is a.property_x.base
+        assert a.property_z == self.model_2_b.property_z
 
         a = self.model_3_a + self.model_3_a
-        assert a.res_x == self.model_3_b.res_x
-        assert a.res_y == self.model_3_b.res_y
-        assert a.res_z == self.model_3_b.res_z
+        assert a.property_x == self.model_3_b.property_x
+        assert a.property_y == self.model_3_b.property_y
+        assert a.property_z == self.model_3_b.property_z
 
         # mu_r and epsilon_r
         a = self.model_mu_a + self.model_mu_a
@@ -360,19 +344,19 @@ def test_sub(self):
 
         # All different cases
         a = self.model_1_b - self.model_1_a
-        assert a.res_x == self.model_1_a.res_x
-        assert a.res_y == self.model_1_a.res_y
-        assert a.res_z.base is a.res_x.base
+        assert a.property_x == self.model_1_a.property_x
+        assert a.property_y == self.model_1_a.property_y
+        assert a.property_z.base is a.property_x.base
 
         a = self.model_2_b - self.model_2_a
-        assert a.res_x == self.model_2_a.res_x
-        assert a.res_y.base is a.res_x.base
-        assert a.res_z == self.model_2_a.res_z
+        assert a.property_x == self.model_2_a.property_x
+        assert a.property_y.base is a.property_x.base
+        assert a.property_z == self.model_2_a.property_z
 
         a = self.model_3_b - self.model_3_a
-        assert a.res_x == self.model_3_a.res_x
-        assert a.res_y == self.model_3_a.res_y
-        assert a.res_z == self.model_3_a.res_z
+        assert a.property_x == self.model_3_a.property_x
+        assert a.property_y == self.model_3_a.property_y
+        assert a.property_z == self.model_3_a.property_z
 
         # mu_r and epsilon_r
         a = self.model_mu_b - self.model_mu_a
diff --git a/tests/test_solver.py b/tests/test_solver.py
index 01157b93..74a13f95 100644
--- a/tests/test_solver.py
+++ b/tests/test_solver.py
@@ -27,7 +27,7 @@ def create_dummy(nx, ny, nz):
 def test_solver_homogeneous(capsys):
     # Regression test for homogeneous halfspace.
     # Not very sophisticated; replace/extend by more detailed tests.
-    dat = REGRES['Data']['res']
+    dat = REGRES['res']
 
     grid = meshes.TensorMesh(**dat['input_grid'])
     model = models.Model(**dat['input_model'])
@@ -45,8 +45,8 @@ def test_solver_homogeneous(capsys):
 
     # Experimental:
     # Check if norms are also the same, at least for first two cycles.
-    assert "1.509e-01  after   1 F-cycles   [9.161e-07, 0.151]   0 0" in out
-    assert "1.002e-01  after   2 F-cycles   [6.082e-07, 0.664]   0 0" in out
+    assert "4.091e-02  after   1 F-cycles   [1.937e-07, 0.041]   0 0" in out
+    assert "4.529e-03  after   2 F-cycles   [2.144e-08, 0.111]   0 0" in out
 
     # Check all fields (ex, ey, and ez)
     assert_allclose(dat['Fresult'], efield)
@@ -158,7 +158,7 @@ def test_solver_homogeneous(capsys):
 
 def test_solver_heterogeneous(capsys):
     # Regression test for heterogeneous case.
-    dat = REGRES['Data']['reg_2']
+    dat = REGRES['reg_2']
     grid = dat['grid']
     model = dat['model']
     sfield = dat['sfield']
@@ -211,7 +211,7 @@ def test_solver_heterogeneous(capsys):
 def test_one_liner(capsys):
     grid = meshes.TensorMesh(
             [np.ones(8), np.ones(8), np.ones(8)], x0=np.array([0, 0, 0]))
-    model = models.Model(grid, res_x=1.5, res_y=1.8, res_z=3.3)
+    model = models.Model(grid, property_x=1.5, property_y=1.8, property_z=3.3)
     sfield = fields.get_source_field(grid, src=[4, 4, 4, 0, 0], freq=10.0)
 
     # Dynamic one-liner.
@@ -238,7 +238,7 @@ def test_one_liner(capsys):
 def test_solver_homogeneous_laplace():
     # Regression test for homogeneous halfspace in Laplace domain.
     # Not very sophisticated; replace/extend by more detailed tests.
-    dat = REGRES['Data']['lap']
+    dat = REGRES['lap']
 
     grid = meshes.TensorMesh(**dat['input_grid'])
     model = models.Model(**dat['input_model'])
@@ -248,13 +248,13 @@ def test_solver_homogeneous_laplace():
     efield = solver.solve(grid, model, sfield, verb=1)
 
     # Check all fields (ex, ey, and ez)
-    assert_allclose(dat['Fresult'], efield)
+    assert_allclose(dat['Fresult'], efield, rtol=5e-6)
 
     # BiCGSTAB with some print checking.
     efield = solver.solve(grid, model, sfield, verb=1, sslsolver=True)
 
     # Check all fields (ex, ey, and ez)
-    assert_allclose(dat['bicresult'], efield)
+    assert_allclose(dat['bicresult'], efield, rtol=5e-6)
 
     # If efield is complex, assert it fails.
     efield = fields.Field(grid, dtype=np.complex_)
@@ -294,9 +294,9 @@ def test_smoothing():
         x = np.arange(1, grid.nCx+1)*2
         y = 1/np.arange(1, grid.nCy+1)
         z = np.arange(1, grid.nCz+1)[::-1]/10
-        res_x = np.outer(np.outer(x, y), z).ravel()
+        property_x = np.outer(np.outer(x, y), z).ravel()
         freq = 0.319
-        model = models.Model(grid, res_x, 0.8*res_x, 2*res_x)
+        model = models.Model(grid, property_x, 0.8*property_x, 2*property_x)
 
         # Create a source field
         sfield = fields.get_source_field(grid=grid, src=src, freq=freq)
@@ -403,9 +403,9 @@ def test_residual():
     x = np.arange(1, grid.nCx+1)*2
     y = 1/np.arange(1, grid.nCy+1)
     z = np.arange(1, grid.nCz+1)[::-1]/10
-    res_x = np.outer(np.outer(x, y), z).ravel()
+    property_x = np.outer(np.outer(x, y), z).ravel()
     freq = 0.319
-    model = models.Model(grid, res_x, 0.8*res_x, 2*res_x)
+    model = models.Model(grid, property_x, 0.8*property_x, 2*property_x)
 
     # Create a source field
     sfield = fields.get_source_field(grid=grid, src=src, freq=freq)
@@ -438,9 +438,11 @@ def test_krylov(capsys):
     # Just check here for bicgstab-error.
 
     # Load any case.
-    dat = REGRES['Data']['res']
+    dat = REGRES['res']
     grid = meshes.TensorMesh(**dat['input_grid'])
     model = models.Model(**dat['input_model'])
+    model.property_x /= 100000  # Set stupid input to make bicgstab fail.
+    model.property_y *= 100000  # Set stupid input to make bicgstab fail.
     sfield = fields.get_source_field(**dat['input_source'])
     vmodel = models.VolumeModel(grid, model, sfield)
     efield = fields.Field(grid)  # Initiate e-field.
@@ -448,8 +450,7 @@ def test_krylov(capsys):
     # Get var-instance
     var = solver.MGParameters(
             cycle=None, sslsolver=True, semicoarsening=False,
-            linerelaxation=False, vnC=grid.vnC, verb=4,
-            maxit=-1,  # Set stupid input to make bicgstab fail.
+            linerelaxation=False, vnC=grid.vnC, verb=4, maxit=-1,
     )
     var.l2_refe = sl.norm(sfield, check_finite=False)