diff --git a/lib/iris/coords.py b/lib/iris/coords.py
index 6402e73263..5519630c5d 100644
--- a/lib/iris/coords.py
+++ b/lib/iris/coords.py
@@ -1013,7 +1013,7 @@ def _sanity_check_bounds(self):
                              'are not defined for coordinates with more than '
                              '2 dimensions.'.format(self.name()))
 
-    def _discontiguity_in_bounds(self, rtol=1e-05, atol=1e-08):
+    def _discontiguity_in_bounds(self, rtol=1e-5, atol=1e-8):
         """
         Checks that the bounds of the coordinate are contiguous.
 
@@ -1064,9 +1064,14 @@ def mod360_adjust(compare_axis):
                         modification = (index.astype(int) * 360) * sign
                         upper_bounds -= modification
 
-                diffs_along_axis = np.abs(upper_bounds - lower_bounds)
-                contiguous_along_axis = np.allclose(upper_bounds, lower_bounds,
-                                                    rtol=rtol, atol=atol)
+                diffs_between_cells = np.abs(upper_bounds - lower_bounds)
+                cell_size = lower_bounds - upper_bounds
+                diffs_along_axis = diffs_between_cells > (atol +
+                                                          rtol * cell_size)
+
+                points_close_enough = diffs_along_axis <= (atol +
+                                                           rtol * cell_size)
+                contiguous_along_axis = np.all(points_close_enough)
                 return diffs_along_axis, contiguous_along_axis
 
             diffs_along_x, match_cell_x1 = mod360_adjust(compare_axis='x')
diff --git a/lib/iris/plot.py b/lib/iris/plot.py
index e3193e96a4..88f89cc79a 100644
--- a/lib/iris/plot.py
+++ b/lib/iris/plot.py
@@ -104,7 +104,7 @@ def get_span(coord):
             raise ValueError('The coordinate {!r} has {} dimensions.'
                              'Cell-based plotting is only supported for'
                              'coordinates with one or two dimensions.'
-                             .format(coord.name()), len(span))
+                             .format(coord.name(), len(span)))
 
     # Check the combination of coordinates spans enough (ndims) data
     # dimensions.
@@ -271,7 +271,7 @@ def _invert_yaxis(v_coord, axes=None):
             axes.invert_yaxis()
 
 
-def _check_bounds_contiguity_and_mask(coord, data, atol=None):
+def _check_bounds_contiguity_and_mask(coord, data, atol=None, rtol=None):
     """
     Checks that any discontiguities in the bounds of the given coordinate only
     occur where the data is masked.
@@ -297,6 +297,7 @@ def _check_bounds_contiguity_and_mask(coord, data, atol=None):
             tolerance.
 
     """
+    kwargs = {}
     data_is_masked = hasattr(data, 'mask')
     if data_is_masked:
         # When checking the location of the discontiguities, we check against
@@ -317,15 +318,16 @@ def _check_bounds_contiguity_and_mask(coord, data, atol=None):
 
     elif coord.ndim == 2:
         if atol:
-            kwargs = {'atol': atol}
-        else:
-            kwargs = {}
+            kwargs['atol'] = atol
+        if rtol:
+            kwargs['rtol'] = rtol
         contiguous, diffs = coord._discontiguity_in_bounds(**kwargs)
 
         if not contiguous and data_is_masked:
             diffs_along_x, diffs_along_y = diffs
 
-            # Check along both dimensions.
+            # Check along both dimensions that any discontiguous
+            # points are correctly masked.
             not_masked_at_discontiguity_along_x = np.any(
                 np.logical_and(mask_invert[:, :-1], diffs_along_x))
 
@@ -341,12 +343,20 @@ def _check_bounds_contiguity_and_mask(coord, data, atol=None):
         if not data_is_masked:
             raise ValueError('The bounds of the {} coordinate are not '
                              'contiguous. Not able to create a suitable grid'
-                             'to plot.'.format(coord.name()))
+                             'to plot. You can use '
+                             'iris.util.find_discontiguities() to identify '
+                             'discontiguities in your x and y coordinate '
+                             'bounds arrays.'.format(coord.name()))
         if not_masked_at_discontiguity:
             raise ValueError('The bounds of the {} coordinate are not '
                              'contiguous and data is not masked where the '
                              'discontiguity occurs. Not able to create a '
-                             'suitable grid to plot.'.format(
+                             'suitable grid to plot. You can use '
+                             'iris.util.find_discontiguities() to identify '
+                             'discontiguities in your x and y coordinate '
+                             'bounds arrays, and then mask them with '
+                             'iris.util.mask_cube()'
+                             ''.format(
                               coord.name()))
 
 
diff --git a/lib/iris/tests/stock/_stock_2d_latlons.py b/lib/iris/tests/stock/_stock_2d_latlons.py
index 395a50d8ea..4fb3979a06 100644
--- a/lib/iris/tests/stock/_stock_2d_latlons.py
+++ b/lib/iris/tests/stock/_stock_2d_latlons.py
@@ -302,17 +302,17 @@ def sample_cube(xargs, yargs):
     return cube
 
 
-def make_bounds_discontiguous_at_point(cube, at_iy, at_ix):
+def make_bounds_discontiguous_at_point(cube, at_iy, at_ix, in_y=False):
     """
-    Meddle with the XY grid bounds of a cube to make the grid discontiguous.
+    Meddle with the XY grid bounds of a 2D cube to make the grid discontiguous.
 
     Changes the points and bounds of a single gridcell, so that it becomes
-    discontinuous with the next gridcell to its right.
-    Also masks the cube data at the given point.
+    discontinuous with an adjacent gridcell : either the one to its right, or
+    the one above it (if 'in_y' is True).
 
-    The cube must have bounded 2d 'x' and 'y' coordinates.
+    Also masks the cube data at the given point.
 
-    TODO: add a switch to make a discontiguity in the *y*-direction instead ?
+    The cube must be 2-dimensional and have bounded 2d 'x' and 'y' coordinates.
 
     """
     x_coord = cube.coord(axis='x')
@@ -326,11 +326,21 @@ def adjust_coord(coord):
         pts, bds = coord.points, coord.bounds
         # Fetch the 4 bounds (bottom-left, bottom-right, top-right, top-left)
         bds_bl, bds_br, bds_tr, bds_tl = bds[at_iy, at_ix]
-        # Make a discontinuity "at" (iy, ix), by moving the right-hand edge of
-        # the cell to the midpoint of the existing left+right bounds.
-        new_bds_br = 0.5 * (bds_bl + bds_br)
-        new_bds_tr = 0.5 * (bds_tl + bds_tr)
-        bds_br, bds_tr = new_bds_br, new_bds_tr
+        if not in_y:
+            # Make a discontinuity "at" (iy, ix), by moving the right-hand edge
+            # of the cell to the midpoint of the existing left+right bounds.
+            new_bds_br = 0.5 * (bds_bl + bds_br)
+            new_bds_tr = 0.5 * (bds_tl + bds_tr)
+            bds_br, bds_tr = new_bds_br, new_bds_tr
+        else:
+            # Same but in the 'grid y direction' :
+            # Make a discontinuity "at" (iy, ix), by moving the **top** edge of
+            # the cell to the midpoint of the existing **top+bottom** bounds.
+            new_bds_tl = 0.5 * (bds_bl + bds_tl)
+            new_bds_tr = 0.5 * (bds_br + bds_tr)
+            bds_tl, bds_tr = new_bds_tl, new_bds_tr
+
+        # Write in the new bounds (all 4 corners).
         bds[at_iy, at_ix] = [bds_bl, bds_br, bds_tr, bds_tl]
         # Also reset the cell midpoint to the middle of the 4 new corners,
         # in case having a midpoint outside the corners might cause a problem.
@@ -341,10 +351,15 @@ def adjust_coord(coord):
 
     adjust_coord(x_coord)
     adjust_coord(y_coord)
-    # Also mask the relevant data point.
-    data = cube.data  # N.B. fetch all the data.
-    if not ma.isMaskedArray(data):
-        # Promote to masked array, to avoid converting mask to NaN.
-        data = ma.masked_array(data)
-    data[at_iy, at_ix] = ma.masked
-    cube.data = data
+
+    # Check which dimensions are spanned by each coordinate.
+    for coord in (x_coord, y_coord):
+        span = set(cube.coord_dims(coord))
+        if not span:
+            msg = 'The coordinate {!r} doesn\'t span a data dimension.'
+            raise ValueError(msg.format(coord.name()))
+
+        masked_data = ma.masked_array(cube.data)
+        masked_data[at_iy, at_ix] = ma.masked
+
+    cube.data = masked_data
diff --git a/lib/iris/tests/unit/coords/test_Coord.py b/lib/iris/tests/unit/coords/test_Coord.py
index 27b468044d..d42448d609 100644
--- a/lib/iris/tests/unit/coords/test_Coord.py
+++ b/lib/iris/tests/unit/coords/test_Coord.py
@@ -480,7 +480,8 @@ def test_2d_discontiguous_along_x(self):
         contiguous, diffs = coord._discontiguity_in_bounds()
         diffs_along_x, diffs_along_y = diffs
         self.assertFalse(contiguous)
-        self.assertArrayEqual(diffs_along_x, np.array([2, 2, 2]).reshape(3, 1))
+        self.assertArrayEqual(diffs_along_x,
+                              np.array([True, True, True]).reshape(3, 1))
         self.assertTrue(not diffs_along_y.any())
 
     def test_2d_discontiguous_along_y(self):
@@ -490,7 +491,7 @@ def test_2d_discontiguous_along_y(self):
         diffs_along_x, diffs_along_y = diffs
         self.assertFalse(contiguous)
         self.assertTrue(not diffs_along_x.any())
-        self.assertArrayEqual(diffs_along_y, np.array([[2, 2, 2]]))
+        self.assertArrayEqual(diffs_along_y, np.array([[True, True, True]]))
 
     def test_2d_discontiguous_along_x_and_y(self):
         coord = AuxCoord(np.array([[1, 5], [3, 5]]),
@@ -498,8 +499,8 @@ def test_2d_discontiguous_along_x_and_y(self):
                                           [[2, 4, 4, 2], [4, 6, 6, 4]]]))
         contiguous, diffs = coord._discontiguity_in_bounds()
         diffs_along_x, diffs_along_y = diffs
-        exp_x_diffs = np.array([2, 0]).reshape(2, 1)
-        exp_y_diffs = np.array([2, 0]).reshape(1, 2)
+        exp_x_diffs = np.array([True, False]).reshape(2, 1)
+        exp_y_diffs = np.array([True, False]).reshape(1, 2)
         self.assertFalse(contiguous)
         self.assertArrayEqual(diffs_along_x, exp_x_diffs)
         self.assertArrayEqual(diffs_along_y, exp_y_diffs)
@@ -508,10 +509,11 @@ def test_2d_contiguous_along_x_atol(self):
         coord = AuxCoord(self.points_3by3[:, ::2],
                          bounds=self.lon_bounds_3by3[:, ::2, :])
         # Set a high atol that allows small discontiguities.
-        contiguous, diffs = coord._discontiguity_in_bounds(atol=2)
+        contiguous, diffs = coord._discontiguity_in_bounds(atol=5)
         diffs_along_x, diffs_along_y = diffs
         self.assertTrue(contiguous)
-        self.assertArrayEqual(diffs_along_x, np.array([2, 2, 2]).reshape(3, 1))
+        self.assertArrayEqual(diffs_along_x,
+                              np.array([False, False, False]).reshape(3, 1))
         self.assertTrue(not diffs_along_y.any())
 
     def test_2d_one_cell(self):
@@ -571,7 +573,7 @@ def test_2d_discontiguous_mod_360(self):
         contiguous, diffs = coord._discontiguity_in_bounds()
         diffs_along_x, diffs_along_y = diffs
         self.assertFalse(contiguous)
-        self.assertArrayEqual(diffs_along_x, np.array([[170], [170]]))
+        self.assertArrayEqual(diffs_along_x, np.array([[True], [True]]))
         self.assertTrue(not diffs_along_y.any())
 
     def test_2d_contiguous_mod_360_not_longitude(self):
@@ -599,7 +601,7 @@ def test_2d_discontiguous_mod_360_not_longitude(self):
         contiguous, diffs = coord._discontiguity_in_bounds()
         diffs_along_x, diffs_along_y = diffs
         self.assertFalse(contiguous)
-        self.assertArrayEqual(diffs_along_x, np.array([[100], [100]]))
+        self.assertArrayEqual(diffs_along_x, np.array([[True], [True]]))
         self.assertTrue(not diffs_along_y.any())
 
 
diff --git a/lib/iris/tests/unit/plot/test__check_bounds_contiguity_and_mask.py b/lib/iris/tests/unit/plot/test__check_bounds_contiguity_and_mask.py
index 225aa2890e..75c2b5be2f 100644
--- a/lib/iris/tests/unit/plot/test__check_bounds_contiguity_and_mask.py
+++ b/lib/iris/tests/unit/plot/test__check_bounds_contiguity_and_mask.py
@@ -33,6 +33,9 @@
 from iris.tests.stock import (sample_2d_latlons,
                               make_bounds_discontiguous_at_point)
 
+from iris.plot import _check_bounds_contiguity_and_mask
+
+
 if tests.MPL_AVAILABLE:
     import iris.plot as iplt
 
@@ -43,20 +46,21 @@ def test_1d_not_checked(self):
         # Test a 1D coordinate, which is not checked as atol is not set.
         coord = DimCoord([1, 3, 5], bounds=[[0, 2], [2, 4], [5, 6]])
         data = np.array([278, 300, 282])
-        iplt._check_bounds_contiguity_and_mask(coord, data)
+        # Make sure contiguity checking doesn't throw an error
+        _check_bounds_contiguity_and_mask(coord, data)
 
     def test_1d_contiguous(self):
-        # Test a 1D coordinate which is contiguous.
+        # Test that a 1D coordinate which is contiguous does not fail.
         coord = DimCoord([1, 3, 5], bounds=[[0, 2], [2, 4], [4, 6]])
         data = np.array([278, 300, 282])
-        iplt._check_bounds_contiguity_and_mask(coord, data, atol=1e-3)
+        _check_bounds_contiguity_and_mask(coord, data, atol=1e-3)
 
     def test_1d_discontigous_masked(self):
         # Test a 1D coordinate which is discontiguous but masked at
         # discontiguities.
         coord = DimCoord([1, 3, 5], bounds=[[0, 2], [2, 4], [5, 6]])
         data = ma.array(np.array([278, 300, 282]), mask=[0, 1, 0])
-        iplt._check_bounds_contiguity_and_mask(coord, data, atol=1e-3)
+        _check_bounds_contiguity_and_mask(coord, data, atol=1e-3)
 
     def test_1d_discontigous_unmasked(self):
         # Test a 1D coordinate which is discontiguous and unmasked at
@@ -66,13 +70,13 @@ def test_1d_discontigous_unmasked(self):
         msg = 'coordinate are not contiguous and data is not masked where ' \
               'the discontiguity occurs'
         with self.assertRaisesRegexp(ValueError, msg):
-            iplt._check_bounds_contiguity_and_mask(coord, data, atol=1e-3)
+            _check_bounds_contiguity_and_mask(coord, data, atol=1e-3)
 
     def test_2d_contiguous(self):
-        # Test a 2D coordinate which is contiguous.
+        # Test that a 2D coordinate which is contiguous does not throw
+        # an error.
         cube = sample_2d_latlons()
-        iplt._check_bounds_contiguity_and_mask(cube.coord('longitude'),
-                                               cube.data)
+        _check_bounds_contiguity_and_mask(cube.coord('longitude'), cube.data)
 
     def test_2d_contiguous_atol(self):
         # Check the atol is passed correctly.
@@ -82,18 +86,16 @@ def test_2d_contiguous_atol(self):
             # Discontiguity returns two objects that are unpacked in
             # `_check_bounds_contiguity_and_mask`.
             discontiguity_check.return_value = [True, None]
-            iplt._check_bounds_contiguity_and_mask(cube.coord('longitude'),
-                                                   cube.data,
-                                                   atol=1e-3)
+            _check_bounds_contiguity_and_mask(cube.coord('longitude'),
+                                              cube.data, atol=1e-3)
         discontiguity_check.assert_called_with(atol=1e-3)
 
     def test_2d_discontigous_masked(self):
-        # Test a 2D coordinate which is discontiguous but masked at
-        # discontiguities.
+        # Test that a 2D coordinate which is discontiguous but masked at
+        # discontiguities doesn't error.
         cube = sample_2d_latlons()
         make_bounds_discontiguous_at_point(cube, 3, 4)
-        iplt._check_bounds_contiguity_and_mask(cube.coord('longitude'),
-                                               cube.data)
+        _check_bounds_contiguity_and_mask(cube.coord('longitude'), cube.data)
 
     def test_2d_discontigous_unmasked(self):
         # Test a 2D coordinate which is discontiguous and unmasked at
@@ -103,8 +105,8 @@ def test_2d_discontigous_unmasked(self):
         msg = 'coordinate are not contiguous'
         cube.data[3, 4] = ma.nomask
         with self.assertRaisesRegexp(ValueError, msg):
-            iplt._check_bounds_contiguity_and_mask(cube.coord('longitude'),
-                                                   cube.data)
+            _check_bounds_contiguity_and_mask(cube.coord('longitude'),
+                                              cube.data)
 
 
 if __name__ == "__main__":
diff --git a/lib/iris/tests/unit/util/test_find_discontiguities.py b/lib/iris/tests/unit/util/test_find_discontiguities.py
new file mode 100644
index 0000000000..eeaae1596c
--- /dev/null
+++ b/lib/iris/tests/unit/util/test_find_discontiguities.py
@@ -0,0 +1,99 @@
+# (C) British Crown Copyright 2018, Met Office
+#
+# This file is part of Iris.
+#
+# Iris is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by the
+# Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# Iris is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with Iris.  If not, see <http://www.gnu.org/licenses/>.
+"""Test function :func:`iris.util.find_discontiguities"""
+
+from __future__ import (absolute_import, division, print_function)
+from six.moves import (filter, input, map, range, zip)  # noqa
+
+# Import iris.tests first so that some things can be initialised before
+# importing anything else.
+import iris.tests as tests
+
+import numpy as np
+
+from iris.tests.stock import simple_3d
+from iris.tests.stock import sample_2d_latlons
+from iris.tests.stock import make_bounds_discontiguous_at_point
+
+from iris.util import find_discontiguities
+
+
+def full2d_global():
+    return sample_2d_latlons(transformed=True)
+
+
+@tests.skip_data
+class Test(tests.IrisTest):
+    def setUp(self):
+        # Set up a 2d lat-lon cube with 2d coordinates that have been
+        # transformed so they are not in a regular lat-lon grid.
+        # Then generate a discontiguity at a single lat-lon point.
+        self.testcube_discontig = full2d_global()
+        make_bounds_discontiguous_at_point(self.testcube_discontig, 3, 3)
+        # Repeat that for a discontiguity in the grid 'Y' direction.
+        self.testcube_discontig_along_y = full2d_global()
+        make_bounds_discontiguous_at_point(self.testcube_discontig_along_y,
+                                           2, 4, in_y=True)
+
+    def test_find_discontiguities(self):
+        # Check that the mask we generate when making the discontiguity
+        # matches that generated by find_discontiguities
+        cube = self.testcube_discontig
+        expected = cube.data.mask
+        returned = find_discontiguities(cube)
+        self.assertTrue(np.all(expected == returned))
+
+    def test_find_discontiguities_in_y(self):
+        # Check that the mask we generate when making the discontiguity
+        # matches that generated by find_discontiguities
+        cube = self.testcube_discontig_along_y
+        expected = cube.data.mask
+        returned = find_discontiguities(cube)
+        self.assertTrue(np.all(expected == returned))
+
+    def test_find_discontiguities_1d_coord(self):
+        # Check that an error is raised when we try and use
+        # find_discontiguities on 1D coordinates:
+        cube = simple_3d()
+        with self.assertRaises(NotImplementedError):
+            find_discontiguities(cube)
+
+    def test_find_discontiguities_with_atol(self):
+        cube = self.testcube_discontig
+        # Choose a very large absolute tolerance which will result in fine
+        # discontiguities being disregarded
+        atol = 100
+        # Construct an array the size of the points array filled with 'False'
+        # to represent a mask showing no discontiguities
+        expected = np.zeros(cube.shape, dtype=bool)
+        returned = find_discontiguities(cube, abs_tol=atol)
+        self.assertTrue(np.all(expected == returned))
+
+    def test_find_discontiguities_with_rtol(self):
+        cube = self.testcube_discontig
+        # Choose a very large relative tolerance which will result in fine
+        # discontiguities being disregarded
+        rtol = 1000
+        # Construct an array the size of the points array filled with 'False'
+        # to represent a mask showing no discontiguities
+        expected = np.zeros(cube.shape, dtype=bool)
+        returned = find_discontiguities(cube, rel_tol=rtol)
+        self.assertTrue(np.all(expected == returned))
+
+
+if __name__ == '__main__':
+    tests.main()
diff --git a/lib/iris/tests/unit/util/test_mask_cube.py b/lib/iris/tests/unit/util/test_mask_cube.py
new file mode 100644
index 0000000000..ce14662174
--- /dev/null
+++ b/lib/iris/tests/unit/util/test_mask_cube.py
@@ -0,0 +1,61 @@
+# (C) British Crown Copyright 2018, Met Office
+#
+# This file is part of Iris.
+#
+# Iris is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by the
+# Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# Iris is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with Iris.  If not, see <http://www.gnu.org/licenses/>.
+"""Test function :func:`iris.util.mask_cube"""
+
+from __future__ import (absolute_import, division, print_function)
+from six.moves import (filter, input, map, range, zip)  # noqa
+
+# Import iris.tests first so that some things can be initialised before
+# importing anything else.
+import iris.tests as tests
+
+import numpy as np
+import numpy.ma as ma
+
+from iris.util import mask_cube
+from iris.tests.stock import sample_2d_latlons
+from iris.tests.stock import make_bounds_discontiguous_at_point
+
+
+def full2d_global():
+    return sample_2d_latlons(transformed=True)
+
+
+@tests.skip_data
+class Test(tests.IrisTest):
+    def setUp(self):
+        # Set up a 2d cube with a masked discontiguity to test masking
+        # of 2-dimensional cubes
+        self.cube_2d = full2d_global()
+        make_bounds_discontiguous_at_point(self.cube_2d, 3, 3)
+
+    def test_mask_cube_2d(self):
+        # This tests the masking of a 2d data array
+        cube = self.cube_2d
+        discontiguity_array = ma.getmaskarray(cube.data).copy()
+        expected = cube.copy()
+
+        # Remove mask so that we can pass an unmasked data set to
+        # mask_discontiguities, and check that it masks the correct point by
+        # comparing with masked data
+        cube.data.mask = ma.nomask
+        returned = mask_cube(cube, discontiguity_array)
+        self.assertTrue(np.all(expected.data.mask == returned.data.mask))
+
+
+if __name__ == '__main__':
+    tests.main()
diff --git a/lib/iris/util.py b/lib/iris/util.py
index 4f4c42dc4b..a12a5e0102 100644
--- a/lib/iris/util.py
+++ b/lib/iris/util.py
@@ -33,13 +33,13 @@
 import os.path
 import sys
 import tempfile
-import time
 
 import cf_units
 import numpy as np
 import numpy.ma as ma
 
 import iris
+import iris.coords
 import iris.exceptions
 
 
@@ -1554,3 +1554,111 @@ def _meshgrid(*xi, **kwargs):
         if mxii.dtype != xii.dtype:
             mxi[i] = mxii.astype(xii.dtype)
     return mxi
+
+
+def find_discontiguities(cube, rel_tol=1e-5, abs_tol=1e-8):
+    """
+    Searches coord for discontiguities in the bounds array, returned as a
+    boolean array (True where discontiguities are present).
+
+    Args:
+
+    * cube (`iris.cube.Cube`):
+        The cube to be checked for discontinuities in its 'x' and 'y'
+        coordinates.
+
+    Kwargs:
+
+    * rel_tol (float):
+        The relative equality tolerance to apply in coordinate bounds
+        checking.
+
+    * abs_tol (float):
+        The absolute value tolerance to apply in coordinate bounds
+        checking.
+
+    Returns:
+
+    * result (`numpy.ndarray` of bool) :
+        true/false map of which cells in the cube XY grid have
+        discontiguities in the coordinate points array.
+
+        This can be used as the input array for
+        :func:`iris.util.mask_discontiguities`.
+
+    Examples::
+
+        # Find any unknown discontiguities in your cube's x and y arrays:
+        discontiguities = iris.util.find_discontiguities(cube)
+
+        # Pass the resultant boolean array to `iris.util.mask_discontiguities`
+        # with a cube slice; this will use the boolean array to mask
+        # any discontiguous data points before plotting:
+        masked_cube_slice = iris.util.mask_discontiguities(cube[0],
+                                                           discontiguities)
+
+        # Plot the masked cube slice:
+        iplt.pcolormesh(masked_cube_slice)
+
+    """
+    lats_and_lons = ['latitude', 'grid_latitude',
+                     'longitude', 'grid_longitude']
+    spatial_coords = [coord for coord in cube.aux_coords if
+                      coord.name() in lats_and_lons]
+    dim_err_msg = 'Discontiguity searches are currently only supported for ' \
+                  '2-dimensional coordinates.'
+    if len(spatial_coords) != 2:
+        raise NotImplementedError(dim_err_msg)
+
+    # Check which dimensions are spanned by each coordinate.
+    for coord in spatial_coords:
+        if coord.ndim != 2:
+            raise NotImplementedError(dim_err_msg)
+        else:
+            span = set(cube.coord_dims(coord))
+        if not span:
+            msg = 'The coordinate {!r} doesn\'t span a data dimension.'
+            raise ValueError(msg.format(coord.name()))
+
+    # Check that the 2d coordinate arrays are the same shape as each other
+    if len(spatial_coords) == 2:
+        assert spatial_coords[0].points.shape == spatial_coords[1].points.shape
+
+    # Set up unmasked boolean array the same size as the coord points array:
+    bad_points_boolean = np.zeros(spatial_coords[0].points.shape, dtype=bool)
+
+    for coord in spatial_coords:
+        _, (diffs_x, diffs_y) = coord._discontiguity_in_bounds(rtol=rel_tol,
+                                                               atol=abs_tol)
+
+        bad_points_boolean[:, :-1] = np.logical_or(bad_points_boolean[:, :-1],
+                                                   diffs_x)
+        # apply mask for y-direction discontiguities:
+        bad_points_boolean[:-1, :] = np.logical_or(bad_points_boolean[:-1, :],
+                                                   diffs_y)
+    return bad_points_boolean
+
+
+def mask_cube(cube, points_to_mask):
+    """
+    Masks any cells in the data array which correspond to cells marked `True`
+    in the `points_to_mask` array.
+
+    Args:
+
+    * cube (`iris.cube.Cube`):
+        A 2-dimensional instance of :class:`iris.cube.Cube`.
+
+    * points_to_mask (`numpy.ndarray` of bool):
+        A 2d boolean array of Truth values representing points to mask in the
+        x and y arrays of the cube.
+
+    Returns:
+
+    * result (`iris.cube.Cube`):
+        A cube whose data array is masked at points specified by input array.
+
+    """
+    cube.data = ma.masked_array(cube.data)
+    cube.data[points_to_mask] = ma.masked
+    return cube