Repository: SciTools/iris
Pull request #6537 (2025-06-27 14:57:50)

Merge `FEATURE_wkt` into `main`
# Merge via merge-commit

(Not via squash-commit)


Patch:
diff --git a/docs/src/further_topics/netcdf_io.rst b/docs/src/further_topics/netcdf_io.rst
index 682918d5f..47d85ffea 100644
--- a/docs/src/further_topics/netcdf_io.rst
+++ b/docs/src/further_topics/netcdf_io.rst
@@ -235,3 +235,160 @@ Worked example:
     >>> my_coord.ignore_axis = True
     >>> print(guess_coord_axis(my_coord))
     None
+
+Multiple Coordinate Systems and Ordered Axes
+--------------------------------------------
+
+In a CF compliant NetCDF file, the coordinate variables associated with a
+data variable can specify a specific *coordinate system* that defines how
+the coordinate values relate to physical locations on the globe. For example,
+a coordinate might have values with units of metres that should be referenced
+against a *Transverse Mercator* projection with a specific origin. This
+information is not stored on the coordinate itself, but in a separate
+*grid mapping* variable. Furthermore, the grid mapping for a set of
+coordinates is associated with the data variable (not the coordinates
+variables) via the ``grid_mapping`` attribute.
+
+For example, a temperature variable defined on a *rotated pole* grid might
+look like this in a NetCDF file (extract of relevant variables):
+
+.. code-block:: text
+
+  float T(rlat,rlon) ;
+    T:long_name = "temperature" ;
+    T:units = "K" ;
+    T:grid_mapping = "rotated_pole" ;
+
+  char rotated_pole ;
+    rotated_pole:grid_mapping_name = "rotated_latitude_longitude" ;
+    rotated_pole:grid_north_pole_latitude = 32.5 ;
+    rotated_pole:grid_north_pole_longitude = 170. ;
+
+  float rlon(rlon) ;
+    rlon:long_name = "longitude in rotated pole grid" ;
+    rlon:units = "degrees" ;
+    rlon:standard_name = "grid_longitude";
+
+  float rlat(rlat) ;
+    rlat:long_name = "latitude in rotated pole grid" ;
+    rlat:units = "degrees" ;
+    rlat:standard_name = "grid_latitude";
+
+
+Note how the ``rotated pole`` grid mapping (coordinate system) is referenced
+from the data variable ``T:grid_mapping = "rotated_pole"`` and is implicitly
+associated with the dimension coordinate variables ``rlat`` and ``rlon``.
+
+
+Since version `1.8 of the CF Conventions
+<https://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html#grid-mappings-and-projections>`_
+, there has been support for a more explicit version of the ``grid_mapping``
+attribute. This allows for **multiple coordinate systems** to be defined for
+a data variable and individual coordinates to be explicitly associated with
+a coordinate system. This is achieved by use of an **extended syntax** in the
+``grid_mapping`` variable of a data variable:
+
+
+.. code-block:: text
+
+  <grid_mapping_var>: <coord_var> [<coord_var>] [<grid_mapping_var>: <coord_var> ...]
+
+where each ``grid_mapping_var`` identifies a grid mapping variable followed by
+the list of associated coordinate variables (``coord_var``). Note that with
+this syntax it is possible to specify multiple coordinate systems for a
+data variable.
+
+For example, consider the following *air pressure* variable that is
+defined on an *OSGB Transverse Mercator grid*:
+
+.. code-block:: text
+
+    float pres(y, x) ;
+        pres:standard_name = "air_pressure" ;
+        pres:units = "Pa" ;
+        pres:coordinates = "lat lon" ;
+        pres:grid_mapping = "crsOSGB: x y crsWGS84: lat lon" ;
+
+    double x(x) ;
+        x:standard_name = "projection_x_coordinate" ;
+        x:units = "m" ;
+
+    double y(y) ;
+        y:standard_name = "projection_y_coordinate" ;
+        y:units = "m" ;
+
+    double lat(y, x) ;
+        lat:standard_name = "latitude" ;
+        lat:units = "degrees_north" ;
+
+    double lon(y, x) ;
+        lon:standard_name = "longitude" ;
+        lon:units = "degrees_east" ;
+
+    int crsOSGB ;
+        crsOSGB:grid_mapping_name = "transverse_mercator" ;
+        crsOSGB:semi_major_axis = 6377563.396 ;
+        crsOSGB:inverse_flattening = 299.3249646 ;
+        <snip>
+
+    int crsWGS84 ;
+        crsWGS84:grid_mapping_name = "latitude_longitude" ;
+        crsWGS84:longitude_of_prime_meridian = 0. ;
+        <snip>
+
+
+The dimension coordinates ``x`` and ``y`` are explicitly defined on
+an a *transverse mercator* grid via the ``crsOSGB`` variable.
+
+However, with the extended grid syntax, it is also possible to define
+a second coordinate system on a standard **latitude_longitude** grid
+and associate it with the auxiliary ``lat`` and ``lon`` coordinates:
+
+::
+
+    pres:grid_mapping = "crsOSGB: x y crsWGS84: lat lon" ;
+
+
+Note, the *order* of the axes in the extended grid mapping specification is
+significant, but only when used in conjunction with a
+`CRS Well Known Text (WKT)`_ representation of the coordinate system where it
+should be consistent with the ``AXES ORDER`` specified in the ``crs_wkt``
+attribute.
+
+
+Effect on loading
+^^^^^^^^^^^^^^^^^
+
+When Iris loads a NetCDF file that uses the extended grid mapping syntax
+it will generate an :class:`iris.coord_systems.CoordSystem` for each
+coordinate system listed and attempt to attach it to the associated
+:class:`iris.coords.Coord` instances on the cube. Currently, Iris considers
+the ``crs_wkt`` supplementary and builds coordinate systems exclusively
+from the ``grid_mapping`` attribute.
+
+The :attr:`iris.cube.Cube.extended_grid_mapping` property will be set to
+``True`` for cubes loaded from NetCDF data variables utilising the extended
+``grid_mapping`` syntax.
+
+Effect on saving
+^^^^^^^^^^^^^^^^
+
+To maintain existing behaviour, saving an :class:`iris.cube.Cube` to
+a netCDF file will default to the "simple" grid mapping syntax, unless
+the cube was loaded from a file using the extended grid mapping syntax.
+If the cube contains multiple coordinate systems, only the coordinate
+system of the dimension coordinate(s) will be specified.
+
+To enable saving of multiple coordinate systems with ordered axes,
+set the :attr:`iris.cube.Cube.extended_grid_mapping` to ``True``.
+This will generate a ``grid_mapping`` attribute using the extended syntax
+to specify all coordinate systems on the cube. The axes ordering of the
+associated coordinate variables will be consistent with that of the
+generated ``crs_wkt`` attribute.
+
+Note, the ``crs_wkt`` attribute will only be generated when the
+extended grid mapping is also written, i.e. when
+``Cube.extended_grid_mapping=True``.
+
+
+.. _CRS Well Known Text (WKT): https://cfconventions.org/Data/cf-conventions/cf-conventions-1.12/cf-conventions.html#use-of-the-crs-well-known-text-format
\ No newline at end of file
diff --git a/lib/iris/cube.py b/lib/iris/cube.py
index d8006c1f2..04060160b 100644
--- a/lib/iris/cube.py
+++ b/lib/iris/cube.py
@@ -2461,6 +2461,32 @@ class Cube(CFVariableMixin):
 
         return result
 
+    def coord_systems(self) -> list[iris.coord_systems.CoordSystem]:
+        """Return a list of all coordinate systems used in cube coordinates."""
+        # Gather list of our unique CoordSystems on cube:
+        coord_systems = ClassDict(iris.coord_systems.CoordSystem)
+        for coord in self.coords():
+            if coord.coord_system:
+                coord_systems.add(coord.coord_system, replace=True)
+
+        return list(coord_systems.values())
+
+    @property
+    def extended_grid_mapping(self) -> bool:
+        """Return True if a cube will use extended grid mapping syntax to write axes order in grid_mapping.
+
+        Only relevant when saving a cube to NetCDF file format.
+
+        For more details see "Grid Mappings and Projections" in the CF Conventions document:
+        https://cfconventions.org/cf-conventions/conformance.html
+        """
+        return self.attributes.get("iris_extended_grid_mapping", False)
+
+    @extended_grid_mapping.setter
+    def extended_grid_mapping(self, ordered: bool) -> None:
+        """Set to True to enable extended grid mapping syntax."""
+        self.attributes["iris_extended_grid_mapping"] = ordered
+
     def _any_meshcoord(self) -> MeshCoord | None:
         """Return a MeshCoord if there are any, else None."""
         mesh_coords = self.coords(mesh_coords=True)
diff --git a/lib/iris/exceptions.py b/lib/iris/exceptions.py
index 450afce3a..870a64901 100644
--- a/lib/iris/exceptions.py
+++ b/lib/iris/exceptions.py
@@ -172,3 +172,9 @@ class DatalessError(ValueError):
             "operations are not currently supported."
         )
         return msg.format(super().__str__())
+
+
+class CFParseError(IrisError):
+    """Raised when a string associated with a CF defined syntax could not be parsed."""
+
+    pass
diff --git a/lib/iris/fileformats/_nc_load_rules/actions.py b/lib/iris/fileformats/_nc_load_rules/actions.py
index ee19b54b4..af84a84fa 100644
--- a/lib/iris/fileformats/_nc_load_rules/actions.py
+++ b/lib/iris/fileformats/_nc_load_rules/actions.py
@@ -203,7 +203,10 @@ def action_provides_grid_mapping(engine, gridmapping_fact):
 
         def build_outer(engine_, cf_var_):
             coordinate_system = builder(engine_, cf_var_)
-            engine_.cube_parts["coordinate_system"] = coordinate_system
+            # We can now handle more than one coordinate_system, so store as dictionary:
+            engine_.cube_parts["coordinate_systems"][cf_var_.cf_name] = (
+                coordinate_system
+            )
 
         # Part 1 - only building - adding takes place downstream in
         #  helpers.build_and_add_dimension/auxiliary_coordinate().
@@ -218,11 +221,8 @@ def action_provides_grid_mapping(engine, gridmapping_fact):
             ),
         )
 
-        # Check there is not an existing one.
-        # ATM this is guaranteed by the caller, "run_actions".
-        assert engine.fact_list("grid-type") == []
-
-        engine.add_fact("grid-type", (grid_mapping_type,))
+        # Store grid-mapping name along with grid-type to match them later on
+        engine.add_fact("grid-type", (var_name, grid_mapping_type))
 
     else:
         message = "Coordinate system not created. Debug info:\n"
@@ -343,7 +343,35 @@ def action_build_dimension_coordinate(engine, providescoord_fact):
         # Non-conforming lon/lat/projection coords will be classed as
         # dim-coords by cf.py, but 'action_provides_coordinate' will give them
         # a coord-type of 'miscellaneous' : hence, they have no coord-system.
-        coord_system = engine.cube_parts.get("coordinate_system")
+        #
+        # At this point, we need to match any "coordinate_system" entries in
+        # the engine to the coord we are building. There are a couple of cases here:
+        #  1. Simple `grid_mapping = crs` is used, in which case
+        #     we should just apply that mapping to all dim coords.
+        #  2. Extended `grid_mapping = crs: coord1 coord2 crs: coord3 coord4`
+        #     is used in which case we need to match the crs to the coord here.
+
+        # We can have multiple coordinate_system, so now stored as a list (note plural key)
+        coord_systems = engine.cube_parts.get("coordinate_systems")
+        coord_system = None
+
+        if len(coord_systems):
+            # Find which coord system applies to this coordinate.
+            cs_mappings = engine.cube_parts.get("coordinate_system_mappings")
+            if cs_mappings and coord_systems:
+                if len(coord_systems) == 1 and None in cs_mappings:
+                    # Simple grid mapping (a single coord_system with no explicit coords)
+                    # Applies to spatial DimCoord(s) only. In this case only one
+                    # coordinate_system will have been built, so just use it.
+                    (coord_system,) = coord_systems.values()
+                    (cs_name,) = cs_mappings.values()
+                else:
+                    # Extended grid mapping, e.g.
+                    #  `grid_mapping = "crs: coord1 coord2 crs: coord3 coord4"`
+                    # We need to search for coord system that references our coordinate.
+                    if cs_name := cs_mappings.get(cf_var.cf_name):
+                        coord_system = coord_systems.get(cs_name, None)
+
         # Translate the specific grid-mapping type to a grid-class
         if coord_system is None:
             succeed = True
@@ -352,8 +380,13 @@ def action_build_dimension_coordinate(engine, providescoord_fact):
             # Get a grid-class from the grid-type
             # i.e. one of latlon/rotated/projected, as for coord_grid_class.
             gridtypes_factlist = engine.fact_list("grid-type")
-            (gridtypes_fact,) = gridtypes_factlist  # only 1 fact
-            (cs_gridtype,) = gridtypes_fact  # fact contains 1 term
+
+            # potentially multiple grid-type facts; find one for CRS varname
+            cs_gridtype = None
+            for fact_cs_name, fact_cs_type in gridtypes_factlist:
+                if fact_cs_name == cs_name:
+                    cs_gridtype = fact_cs_type
+
             if cs_gridtype == "latitude_longitude":
                 cs_gridclass = "latlon"
             elif cs_gridtype == "rotated_latitude_longitude":
@@ -446,6 +479,7 @@ def action_build_auxiliary_coordinate(engine, auxcoord_fact):
     """Convert a CFAuxiliaryCoordinateVariable into a cube aux-coord."""
     (var_name,) = auxcoord_fact
     rule_name = "fc_build_auxiliary_coordinate"
+    cf_var = engine.cf_var.cf_group[var_name]
 
     # Identify any known coord "type" : latitude/longitude/time/time_period
     # If latitude/longitude, this sets the standard_name of the built AuxCoord
@@ -473,8 +507,27 @@ def action_build_auxiliary_coordinate(engine, auxcoord_fact):
     if coord_type:
         rule_name += f"_{coord_type}"
 
+    # Check if we have a coord_system specified for this coordinate.
+    # (Only possible via extended grid_mapping attribute)
+    coord_systems = engine.cube_parts.get("coordinate_systems")
+    coord_system = None
+
+    cs_mappings = engine.cube_parts.get("coordinate_system_mappings", None)
+    if cs_mappings and coord_systems:
+        if len(coord_systems) == 1 and None in cs_mappings:
+            # Simple grid_mapping - doesn't apply to AuxCoords (we need an explicit mapping)
+            pass
+        else:
+            # Extended grid mapping, e.g.
+            #  `grid_mapping = "crs: coord1 coord2 crs: coord3 coord4"`
+            # We need to search for coord system that references our coordinate.
+            if cs_name := cs_mappings.get(cf_var.cf_name):
+                coord_system = coord_systems.get(cs_name, None)
+
     cf_var = engine.cf_var.cf_group.auxiliary_coordinates[var_name]
-    hh.build_and_add_auxiliary_coordinate(engine, cf_var, coord_name=coord_name)
+    hh.build_and_add_auxiliary_coordinate(
+        engine, cf_var, coord_name=coord_name, coord_system=coord_system
+    )
 
     return rule_name
 
@@ -615,10 +668,9 @@ def run_actions(engine):
     # default (all cubes) action, always runs
     action_default(engine)  # This should run the default rules.
 
-    # deal with grid-mappings
+    # deal with grid-mappings; potentially multiple mappings if extended grid_mapping used.
     grid_mapping_facts = engine.fact_list("grid_mapping")
-    # For now, there should be at most *one* of these.
-    assert len(grid_mapping_facts) in (0, 1)
+
     for grid_mapping_fact in grid_mapping_facts:
         action_provides_grid_mapping(engine, grid_mapping_fact)
 
diff --git a/lib/iris/fileformats/_nc_load_rules/helpers.py b/lib/iris/fileformats/_nc_load_rules/helpers.py
index fc4762594..35c2e9692 100644
--- a/lib/iris/fileformats/_nc_load_rules/helpers.py
+++ b/lib/iris/fileformats/_nc_load_rules/helpers.py
@@ -143,6 +143,46 @@ CF_GRID_MAPPING_GEOSTATIONARY = "geostationary"
 CF_GRID_MAPPING_OBLIQUE = "oblique_mercator"
 CF_GRID_MAPPING_ROTATED_MERCATOR = "rotated_mercator"
 
+#
+# Regex for parsing grid_mapping (extended format)
+# Link to online regex101 playground: https://regex101.com/r/NcKzkQ/1
+#
+#   (\w+):                # Matches '<word>:' and stores in CAPTURE GROUP 1
+#   (                     # CAPTURE GROUP 2 for capturing multiple coords
+#       (?:               #  Non-capturing group for composing match
+#           \s+           #   Matches one or more blank characters
+#           (?!\w+:)      #   Negative look-ahead: don't match <word> followed by colon
+#           \w+           #   Matches a <word>
+#       )+                #  Repeats non-capturing group at least once.
+#   )                     # End of CAPTURE GROUP 2
+_GRID_MAPPING_PARSE_EXTENDED = re.compile(
+    r"""
+        (\w+):
+        (
+            (?:
+                \s+
+                (?!\w+:)
+                \w+
+            )+
+        )+
+    """,
+    re.VERBOSE,
+)
+_GRID_MAPPING_PARSE_SIMPLE = re.compile(r"^\w+$")
+_GRID_MAPPING_VALIDATORS = (  # fmt: skip
+    (
+        re.compile(r"\w+: +\w+:"),
+        "`<coord_system>:` identifier followed immediately by another `<coord_system>:` identifier",
+    ),
+    (
+        re.compile(r"\w+: *$"),
+        "`<coord_system>:` is empty - missing coordinate list",
+    ),
+    (
+        re.compile(r"^\w+ +\w+"),
+        "Multiple coordinates found without `<coord_system>:` identifier",
+    ),
+)
 #
 # CF Attribute Names.
 #
@@ -1936,6 +1976,41 @@ def is_grid_mapping(engine, cf_name, grid_mapping):
     return is_valid
 
 
+################################################################################
+def _parse_extended_grid_mapping(grid_mapping: str) -> dict[None | str, str]:
+    """Parse `grid_mapping` attribute and return list of coordinate system variables and associated coords."""
+    # Handles extended grid_mapping too. Possibilities:
+    #  grid_mapping = "crs"  : simple mapping; a single variable name with no coords
+    #  grid_mapping = "crs: lat lon"  : extended mapping; a variable name and list of coords
+    #  grid_mapping = "crs: lat lon other: var1 var2"  : multiple extended mappings
+
+    mappings: dict[None | str, str]
+
+    # try simple mapping first
+    if _GRID_MAPPING_PARSE_SIMPLE.match(grid_mapping):
+        mappings = {None: grid_mapping}  # simple single grid mapping variable
+    else:
+        # Try extended mapping:
+        # 1. Run validators to check for invalid expressions:
+        for v_re, v_msg in _GRID_MAPPING_VALIDATORS:
+            if len(match := v_re.findall(grid_mapping)):
+                msg = f"Invalid syntax in extended grid_mapping: {grid_mapping!r}\n{v_msg} : {match}"
+                raise iris.exceptions.CFParseError(msg)
+
+        # 2. Parse grid_mapping into list of [cs, (coords, ...)]:
+        result = _GRID_MAPPING_PARSE_EXTENDED.findall(grid_mapping)
+        if len(result) == 0:
+            msg = f"Failed to parse grid_mapping: {grid_mapping!r}"
+            raise iris.exceptions.CFParseError(msg)
+
+        # split second match group into list of coordinates:
+        mappings = {}
+        for cs, coords in result:
+            mappings.update({coord: cs for coord in coords.split()})
+
+    return mappings
+
+
 ################################################################################
 def _is_rotated(engine, cf_name, cf_attr_value):
     """Determine whether the CF coordinate variable is rotated."""
diff --git a/lib/iris/fileformats/cf.py b/lib/iris/fileformats/cf.py
index e69d686c0..23186072a 100644
--- a/lib/iris/fileformats/cf.py
+++ b/lib/iris/fileformats/cf.py
@@ -24,6 +24,8 @@ import warnings
 import numpy as np
 import numpy.ma as ma
 
+import iris.exceptions
+import iris.fileformats._nc_load_rules.helpers as hh
 from iris.fileformats.netcdf import _thread_safe_nc
 from iris.mesh.components import Connectivity
 import iris.util
@@ -654,7 +656,9 @@ class CFGridMappingVariable(CFVariable):
     cf_identity = "grid_mapping"
 
     @classmethod
-    def identify(cls, variables, ignore=None, target=None, warn=True):
+    def identify(
+        cls, variables, ignore=None, target=None, warn=True, coord_system_mappings=None
+    ):
         result = {}
         ignore, target = cls._identify_common(variables, ignore, target)
 
@@ -664,19 +668,61 @@ class CFGridMappingVariable(CFVariable):
             nc_var_att = getattr(nc_var, cls.cf_identity, None)
 
             if nc_var_att is not None:
-                name = nc_var_att.strip()
-
-                if name not in ignore:
-                    if name not in variables:
-                        if warn:
-                            message = "Missing CF-netCDF grid mapping variable %r, referenced by netCDF variable %r"
-                            warnings.warn(
-                                message % (name, nc_var_name),
-                                category=iris.warnings.IrisCfMissingVarWarning,
-                            )
-                    else:
-                        result[name] = CFGridMappingVariable(name, variables[name])
+                # All `grid_mapping` attributes will already have been parsed prior
+                # to `identify` being called and passed in as an argument. We can
+                # ignore the attribute here (it's just used to identify that a grid
+                # mapping exists for this data variable) and get the pre-parsed
+                # mapping from the `coord_mapping_systems` keyword:
+                cs_mappings = None
+                if coord_system_mappings:
+                    cs_mappings = coord_system_mappings.get(nc_var_name, None)
+
+                if not cs_mappings:
+                    # If cs_mappings is None, some parse error must have occurred and the
+                    # user will have already been warned by `_parse_extended_grid_mappings`
+                    continue
 
+                # group the cs_mappings by coordinate system, as we want to iterate over coord systems:
+                uniq_cs = set(cs_mappings.values())
+                cs_coord_mappings = {
+                    cs: [
+                        coord
+                        for coord, coord_cs in cs_mappings.items()
+                        if cs == coord_cs
+                    ]
+                    for cs in uniq_cs
+                }
+
+                for name, coords in cs_coord_mappings.items():
+                    if name not in ignore:
+                        if name not in variables:
+                            if warn:
+                                message = "Missing CF-netCDF grid mapping variable %r, referenced by netCDF variable %r"
+                                warnings.warn(
+                                    message % (name, nc_var_name),
+                                    category=iris.warnings.IrisCfMissingVarWarning,
+                                )
+                        else:
+                            # For extended grid_mapping, also check coord references exist:
+                            has_a_valid_coord = False
+                            if coords:
+                                for coord_name in coords:
+                                    # coord_name could be None if simple grid_mapping is used.
+                                    if coord_name is None or (
+                                        coord_name and coord_name in variables
+                                    ):
+                                        has_a_valid_coord = True
+                                    else:
+                                        message = "Missing CF-netCDF coordinate variable %r (associated with grid mapping variable %r), referenced by netCDF variable %r"
+                                        warnings.warn(
+                                            message % (coord_name, name, nc_var_name),
+                                            category=iris.warnings.IrisCfMissingVarWarning,
+                                        )
+                            #  Only add as a CFGridMappingVariable if at least one of its referenced coords exists:
+                            if has_a_valid_coord:
+                                result[name] = CFGridMappingVariable(
+                                    name, variables[name]
+                                )
         return result
 
 
@@ -1327,6 +1373,9 @@ class CFReader:
         #: Collection of CF-netCDF variables associated with this netCDF file
         self.cf_group = self.CFGroup()
 
+        # Result of parsing "grid_mapping" attribute; mapping of coordinate_system => coordinates
+        self._coord_system_mappings = {}
+
         # Issue load optimisation warning.
         if warn and self._dataset.file_format in [
             "NETCDF3_CLASSIC",
@@ -1395,6 +1444,18 @@ class CFReader:
         """Classify the netCDF variables into CF-netCDF variables."""
         netcdf_variable_names = list(variables.keys())
 
+        # Parse all instances of "grid_mapping" attributes and store in CFReader
+        # This avoids re-parsing the grid_mappings each time they are needed.
+        for nc_var in variables.values():
+            if grid_mapping_attr := getattr(nc_var, "grid_mapping", None):
+                try:
+                    cs_mappings = hh._parse_extended_grid_mapping(grid_mapping_attr)
+                    self._coord_system_mappings[nc_var.name] = cs_mappings
+                except iris.exceptions.CFParseError as e:
+                    msg = f"Error parsing `grid_mapping` attribute for {nc_var.name}: {str(e)}"
+                    warnings.warn(msg, category=iris.warnings.IrisCfWarning)
+                    continue
+
         # Identify all CF coordinate variables first. This must be done
         # first as, by CF convention, the definition of a CF auxiliary
         # coordinate variable may include a scalar CF coordinate variable,
@@ -1413,7 +1474,15 @@ class CFReader:
                 if issubclass(variable_type, CFGridMappingVariable)
                 else coordinate_names
             )
-            self.cf_group.update(variable_type.identify(variables, ignore=ignore))
+            kwargs = (
+                {"coord_system_mappings": self._coord_system_mappings}
+                if issubclass(variable_type, CFGridMappingVariable)
+                else {}
+            )
+
+            self.cf_group.update(
+                variable_type.identify(variables, ignore=ignore, **kwargs)
+            )
 
         # Identify global netCDF attributes.
         attr_dict = {
@@ -1555,12 +1624,18 @@ class CFReader:
             # Build CF variable relationships.
             for variable_type in self._variable_types:
                 ignore = []
+                kwargs = {}
                 # Avoid UGridAuxiliaryCoordinateVariables also being
                 # processed as CFAuxiliaryCoordinateVariables.
                 if not is_mesh_var:
                     ignore += ugrid_coord_names
                 # Prevent grid mapping variables being mis-identified as CF coordinate variables.
-                if not issubclass(variable_type, CFGridMappingVariable):
+                if issubclass(variable_type, CFGridMappingVariable):
+                    # pass parsed grid_mappings to CFGridMappingVariable types
+                    kwargs.update(
+                        {"coord_system_mappings": self._coord_system_mappings}
+                    )
+                else:
                     ignore += coordinate_names
 
                 match = variable_type.identify(
@@ -1568,6 +1643,7 @@ class CFReader:
                     ignore=ignore,
                     target=cf_variable.cf_name,
                     warn=False,
+                    **kwargs,
                 )
                 # Sanity check dimensionality coverage.
                 for cf_name in match:
diff --git a/lib/iris/fileformats/netcdf/loader.py b/lib/iris/fileformats/netcdf/loader.py
index b7895c4dd..216df6759 100644
--- a/lib/iris/fileformats/netcdf/loader.py
+++ b/lib/iris/fileformats/netcdf/loader.py
@@ -79,6 +79,12 @@ def _assert_case_specific_facts(engine, cf, cf_group):
     engine.cube_parts["coordinates"] = []
     engine.cube_parts["cell_measures"] = []
     engine.cube_parts["ancillary_variables"] = []
+    engine.cube_parts["coordinate_systems"] = {}
+
+    # Add the parsed coordinate reference system mappings
+    engine.cube_parts["coordinate_system_mappings"] = cf._coord_system_mappings.get(
+        engine.cf_var.cf_name, None
+    )
 
     # Assert facts for CF coordinates.
     for cf_name in cf_group.coordinates.keys():
@@ -445,6 +451,13 @@ def _load_cube_inner(engine, cf, cf_var, filename):
         for method in cube.cell_methods
     ]
 
+    # Set extended_grid_mapping property ONLY if extended grid_mapping was used.
+    # This avoids having an unnecessary `iris_extended_grid_mapping` attribute entry.
+    if cs_mappings := engine.cube_parts.get("coordinate_system_mappings", None):
+        # `None` as a mapping key implies simple mapping syntax (single coord system)
+        if None not in cs_mappings:
+            cube.extended_grid_mapping = True
+
     if DEBUG:
         # Show activation statistics for this data-var (i.e. cube).
         _actions_activation_stats(engine, cf_var.cf_name)
diff --git a/lib/iris/fileformats/netcdf/saver.py b/lib/iris/fileformats/netcdf/saver.py
index 95bd92d8a..553bb1e19 100644
--- a/lib/iris/fileformats/netcdf/saver.py
+++ b/lib/iris/fileformats/netcdf/saver.py
@@ -104,7 +104,12 @@ _CF_DATA_ATTRS = [
 _CF_GLOBAL_ATTRS = ["conventions", "featureType", "history", "title"]
 
 # UKMO specific attributes that should not be global.
-_UKMO_DATA_ATTRS = ["STASH", "um_stash_source", "ukmo__process_flags"]
+_UKMO_DATA_ATTRS = [
+    "STASH",
+    "um_stash_source",
+    "ukmo__process_flags",
+    "iris_extended_grid_mapping",
+]
 
 # TODO: whenever we advance to CF-1.11 we should then discuss a completion date
 #  for the deprecation of Rotated Mercator in coord_systems.py and
@@ -1934,6 +1939,214 @@ class Saver:
 
         return " ".join(cell_methods)
 
+    def _add_grid_mapping_to_dataset(self, cs, extended_grid_mapping=False):
+        """Create a CF-netCDF grid mapping variable and add to the dataset.
+
+        Parameters
+        ----------
+        cs : :class:`iris.coord_system.CoordSystem`
+            The :class:`iris.coord_system.CoordSystem` used to generate the
+            new netCDF CF grid mapping variable.
+
+        Returns
+        -------
+        None
+        """
+        cf_var_grid = self._dataset.createVariable(cs.grid_mapping_name, np.int32)
+        _setncattr(cf_var_grid, "grid_mapping_name", cs.grid_mapping_name)
+
+        def add_ellipsoid(ellipsoid):
+            cf_var_grid.longitude_of_prime_meridian = (
+                ellipsoid.longitude_of_prime_meridian
+            )
+            semi_major = ellipsoid.semi_major_axis
+            semi_minor = ellipsoid.semi_minor_axis
+            if semi_minor == semi_major:
+                cf_var_grid.earth_radius = semi_major
+            else:
+                cf_var_grid.semi_major_axis = semi_major
+                cf_var_grid.semi_minor_axis = semi_minor
+            if ellipsoid.datum is not None:
+                cf_var_grid.horizontal_datum_name = ellipsoid.datum
+
+        # latlon
+        if isinstance(cs, iris.coord_systems.GeogCS):
+            add_ellipsoid(cs)
+
+        # rotated latlon
+        elif isinstance(cs, iris.coord_systems.RotatedGeogCS):
+            if cs.ellipsoid:
+                add_ellipsoid(cs.ellipsoid)
+            cf_var_grid.grid_north_pole_latitude = cs.grid_north_pole_latitude
+            cf_var_grid.grid_north_pole_longitude = cs.grid_north_pole_longitude
+            cf_var_grid.north_pole_grid_longitude = cs.north_pole_grid_longitude
+
+        # tmerc
+        elif isinstance(cs, iris.coord_systems.TransverseMercator):
+            if cs.ellipsoid:
+                add_ellipsoid(cs.ellipsoid)
+            cf_var_grid.longitude_of_central_meridian = cs.longitude_of_central_meridian
+            cf_var_grid.latitude_of_projection_origin = cs.latitude_of_projection_origin
+            cf_var_grid.false_easting = cs.false_easting
+            cf_var_grid.false_northing = cs.false_northing
+            cf_var_grid.scale_factor_at_central_meridian = (
+                cs.scale_factor_at_central_meridian
+            )
+
+        # merc
+        elif isinstance(cs, iris.coord_systems.Mercator):
+            if cs.ellipsoid:
+                add_ellipsoid(cs.ellipsoid)
+            cf_var_grid.longitude_of_projection_origin = (
+                cs.longitude_of_projection_origin
+            )
+            cf_var_grid.false_easting = cs.false_easting
+            cf_var_grid.false_northing = cs.false_northing
+            # Only one of these should be set
+            if cs.standard_parallel is not None:
+                cf_var_grid.standard_parallel = cs.standard_parallel
+            elif cs.scale_factor_at_projection_origin is not None:
+                cf_var_grid.scale_factor_at_projection_origin = (
+                    cs.scale_factor_at_projection_origin
+                )
+
+        # lcc
+        elif isinstance(cs, iris.coord_systems.LambertConformal):
+            if cs.ellipsoid:
+                add_ellipsoid(cs.ellipsoid)
+            cf_var_grid.standard_parallel = cs.secant_latitudes
+            cf_var_grid.latitude_of_projection_origin = cs.central_lat
+            cf_var_grid.longitude_of_central_meridian = cs.central_lon
+            cf_var_grid.false_easting = cs.false_easting
+            cf_var_grid.false_northing = cs.false_northing
+
+        # polar stereo (have to do this before Stereographic because it subclasses it)
+        elif isinstance(cs, iris.coord_systems.PolarStereographic):
+            if cs.ellipsoid:
+                add_ellipsoid(cs.ellipsoid)
+            cf_var_grid.latitude_of_projection_origin = cs.central_lat
+            cf_var_grid.straight_vertical_longitude_from_pole = cs.central_lon
+            cf_var_grid.false_easting = cs.false_easting
+            cf_var_grid.false_northing = cs.false_northing
+            # Only one of these should be set
+            if cs.true_scale_lat is not None:
+                cf_var_grid.true_scale_lat = cs.true_scale_lat
+            elif cs.scale_factor_at_projection_origin is not None:
+                cf_var_grid.scale_factor_at_projection_origin = (
+                    cs.scale_factor_at_projection_origin
+                )
+            else:
+                cf_var_grid.scale_factor_at_projection_origin = 1.0
+
+        # stereo
+        elif isinstance(cs, iris.coord_systems.Stereographic):
+            if cs.ellipsoid:
+                add_ellipsoid(cs.ellipsoid)
+            cf_var_grid.longitude_of_projection_origin = cs.central_lon
+            cf_var_grid.latitude_of_projection_origin = cs.central_lat
+            cf_var_grid.false_easting = cs.false_easting
+            cf_var_grid.false_northing = cs.false_northing
+            # Only one of these should be set
+            if cs.true_scale_lat is not None:
+                msg = (
+                    "It is not valid CF to save a true_scale_lat for "
+                    "a Stereographic grid mapping."
+                )
+                raise ValueError(msg)
+            elif cs.scale_factor_at_projection_origin is not None:
+                cf_var_grid.scale_factor_at_projection_origin = (
+                    cs.scale_factor_at_projection_origin
+                )
+            else:
+                cf_var_grid.scale_factor_at_projection_origin = 1.0
+
+        # osgb (a specific tmerc)
+        elif isinstance(cs, iris.coord_systems.OSGB):
+            warnings.warn(
+                "OSGB coordinate system not yet handled",
+                category=iris.warnings.IrisSaveWarning,
+            )
+
+        # lambert azimuthal equal area
+        elif isinstance(cs, iris.coord_systems.LambertAzimuthalEqualArea):
+            if cs.ellipsoid:
+                add_ellipsoid(cs.ellipsoid)
+            cf_var_grid.longitude_of_projection_origin = (
+                cs.longitude_of_projection_origin
+            )
+            cf_var_grid.latitude_of_projection_origin = cs.latitude_of_projection_origin
+            cf_var_grid.false_easting = cs.false_easting
+            cf_var_grid.false_northing = cs.false_northing
+
+        # albers conical equal area
+        elif isinstance(cs, iris.coord_systems.AlbersEqualArea):
+            if cs.ellipsoid:
+                add_ellipsoid(cs.ellipsoid)
+            cf_var_grid.longitude_of_central_meridian = cs.longitude_of_central_meridian
+            cf_var_grid.latitude_of_projection_origin = cs.latitude_of_projection_origin
+            cf_var_grid.false_easting = cs.false_easting
+            cf_var_grid.false_northing = cs.false_northing
+            cf_var_grid.standard_parallel = cs.standard_parallels
+
+        # vertical perspective
+        elif isinstance(cs, iris.coord_systems.VerticalPerspective):
+            if cs.ellipsoid:
+                add_ellipsoid(cs.ellipsoid)
+            cf_var_grid.longitude_of_projection_origin = (
+                cs.longitude_of_projection_origin
+            )
+            cf_var_grid.latitude_of_projection_origin = cs.latitude_of_projection_origin
+            cf_var_grid.false_easting = cs.false_easting
+            cf_var_grid.false_northing = cs.false_northing
+            cf_var_grid.perspective_point_height = cs.perspective_point_height
+
+        # geostationary
+        elif isinstance(cs, iris.coord_systems.Geostationary):
+            if cs.ellipsoid:
+                add_ellipsoid(cs.ellipsoid)
+            cf_var_grid.longitude_of_projection_origin = (
+                cs.longitude_of_projection_origin
+            )
+            cf_var_grid.latitude_of_projection_origin = cs.latitude_of_projection_origin
+            cf_var_grid.false_easting = cs.false_easting
+            cf_var_grid.false_northing = cs.false_northing
+            cf_var_grid.perspective_point_height = cs.perspective_point_height
+            cf_var_grid.sweep_angle_axis = cs.sweep_angle_axis
+
+        # oblique mercator (and rotated variant)
+        # Use duck-typing over isinstance() - subclasses (i.e.
+        #  RotatedMercator) upset mock tests.
+        elif getattr(cs, "grid_mapping_name", None) == "oblique_mercator":
+            # RotatedMercator subclasses ObliqueMercator, and RM
+            #  instances are implicitly saved as OM due to inherited
+            #  properties. This is correct because CF 1.11 is removing
+            #  all mention of RM.
+            if cs.ellipsoid:
+                add_ellipsoid(cs.ellipsoid)
+            cf_var_grid.azimuth_of_central_line = cs.azimuth_of_central_line
+            cf_var_grid.latitude_of_projection_origin = cs.latitude_of_projection_origin
+            cf_var_grid.longitude_of_projection_origin = (
+                cs.longitude_of_projection_origin
+            )
+            cf_var_grid.false_easting = cs.false_easting
+            cf_var_grid.false_northing = cs.false_northing
+            cf_var_grid.scale_factor_at_projection_origin = (
+                cs.scale_factor_at_projection_origin
+            )
+
+        # other
+        else:
+            warnings.warn(
+                "Unable to represent the horizontal "
+                "coordinate system. The coordinate system "
+                "type %r is not yet implemented." % type(cs),
+                category=iris.warnings.IrisSaveWarning,
+            )
+
+        # add WKT string
+        if extended_grid_mapping:
+            cf_var_grid.crs_wkt = cs.as_cartopy_crs().to_wkt()
+
     def _create_cf_grid_mapping(self, cube, cf_var_cube):
         """Create CF-netCDF grid mapping and associated CF-netCDF variable.
 
@@ -1953,227 +2166,108 @@ class Saver:
         None
 
         """
-        cs = cube.coord_system("CoordSystem")
-        if cs is not None:
+        coord_systems = []
+        if cube.extended_grid_mapping:
+            # get unique list of all coord_systems on cube coords:
+            for coord in cube.coords():
+                if coord.coord_system and coord.coord_system not in coord_systems:
+                    coord_systems.append(coord.coord_system)
+        else:
+            # will only return a single coord system (if one exists):
+            if cs := cube.coord_system("CoordSystem"):
+                coord_systems.append(cs)
+
+        grid_mappings = {}
+        matched_all_coords = True
+
+        for cs in coord_systems:
             # Grid var not yet created?
             if cs not in self._coord_systems:
+                # handle potential duplicate netCDF variable names:
                 while cs.grid_mapping_name in self._dataset.variables:
                     aname = self._increment_name(cs.grid_mapping_name)
                     cs.grid_mapping_name = aname
 
-                cf_var_grid = self._dataset.createVariable(
-                    cs.grid_mapping_name, np.int32
+                # create grid mapping variable on dataset for this coordinate system:
+                self._add_grid_mapping_to_dataset(
+                    cs, extended_grid_mapping=cube.extended_grid_mapping
                 )
-                _setncattr(cf_var_grid, "grid_mapping_name", cs.grid_mapping_name)
-
-                def add_ellipsoid(ellipsoid):
-                    cf_var_grid.longitude_of_prime_meridian = (
-                        ellipsoid.longitude_of_prime_meridian
-                    )
-                    semi_major = ellipsoid.semi_major_axis
-                    semi_minor = ellipsoid.semi_minor_axis
-                    if semi_minor == semi_major:
-                        cf_var_grid.earth_radius = semi_major
-                    else:
-                        cf_var_grid.semi_major_axis = semi_major
-                        cf_var_grid.semi_minor_axis = semi_minor
-                    if ellipsoid.datum is not None:
-                        cf_var_grid.horizontal_datum_name = ellipsoid.datum
-
-                # latlon
-                if isinstance(cs, iris.coord_systems.GeogCS):
-                    add_ellipsoid(cs)
-
-                # rotated latlon
-                elif isinstance(cs, iris.coord_systems.RotatedGeogCS):
-                    if cs.ellipsoid:
-                        add_ellipsoid(cs.ellipsoid)
-                    cf_var_grid.grid_north_pole_latitude = cs.grid_north_pole_latitude
-                    cf_var_grid.grid_north_pole_longitude = cs.grid_north_pole_longitude
-                    cf_var_grid.north_pole_grid_longitude = cs.north_pole_grid_longitude
-
-                # tmerc
-                elif isinstance(cs, iris.coord_systems.TransverseMercator):
-                    if cs.ellipsoid:
-                        add_ellipsoid(cs.ellipsoid)
-                    cf_var_grid.longitude_of_central_meridian = (
-                        cs.longitude_of_central_meridian
-                    )
-                    cf_var_grid.latitude_of_projection_origin = (
-                        cs.latitude_of_projection_origin
-                    )
-                    cf_var_grid.false_easting = cs.false_easting
-                    cf_var_grid.false_northing = cs.false_northing
-                    cf_var_grid.scale_factor_at_central_meridian = (
-                        cs.scale_factor_at_central_meridian
-                    )
-
-                # merc
-                elif isinstance(cs, iris.coord_systems.Mercator):
-                    if cs.ellipsoid:
-                        add_ellipsoid(cs.ellipsoid)
-                    cf_var_grid.longitude_of_projection_origin = (
-                        cs.longitude_of_projection_origin
-                    )
-                    cf_var_grid.false_easting = cs.false_easting
-                    cf_var_grid.false_northing = cs.false_northing
-                    # Only one of these should be set
-                    if cs.standard_parallel is not None:
-                        cf_var_grid.standard_parallel = cs.standard_parallel
-                    elif cs.scale_factor_at_projection_origin is not None:
-                        cf_var_grid.scale_factor_at_projection_origin = (
-                            cs.scale_factor_at_projection_origin
-                        )
+                self._coord_systems.append(cs)
 
-                # lcc
-                elif isinstance(cs, iris.coord_systems.LambertConformal):
-                    if cs.ellipsoid:
-                        add_ellipsoid(cs.ellipsoid)
-                    cf_var_grid.standard_parallel = cs.secant_latitudes
-                    cf_var_grid.latitude_of_projection_origin = cs.central_lat
-                    cf_var_grid.longitude_of_central_meridian = cs.central_lon
-                    cf_var_grid.false_easting = cs.false_easting
-                    cf_var_grid.false_northing = cs.false_northing
-
-                # polar stereo (have to do this before Stereographic because it subclasses it)
-                elif isinstance(cs, iris.coord_systems.PolarStereographic):
-                    if cs.ellipsoid:
-                        add_ellipsoid(cs.ellipsoid)
-                    cf_var_grid.latitude_of_projection_origin = cs.central_lat
-                    cf_var_grid.straight_vertical_longitude_from_pole = cs.central_lon
-                    cf_var_grid.false_easting = cs.false_easting
-                    cf_var_grid.false_northing = cs.false_northing
-                    # Only one of these should be set
-                    if cs.true_scale_lat is not None:
-                        cf_var_grid.true_scale_lat = cs.true_scale_lat
-                    elif cs.scale_factor_at_projection_origin is not None:
-                        cf_var_grid.scale_factor_at_projection_origin = (
-                            cs.scale_factor_at_projection_origin
-                        )
-                    else:
-                        cf_var_grid.scale_factor_at_projection_origin = 1.0
-
-                # stereo
-                elif isinstance(cs, iris.coord_systems.Stereographic):
-                    if cs.ellipsoid:
-                        add_ellipsoid(cs.ellipsoid)
-                    cf_var_grid.longitude_of_projection_origin = cs.central_lon
-                    cf_var_grid.latitude_of_projection_origin = cs.central_lat
-                    cf_var_grid.false_easting = cs.false_easting
-                    cf_var_grid.false_northing = cs.false_northing
-                    # Only one of these should be set
-                    if cs.true_scale_lat is not None:
+            # create the `grid_mapping` attribute for the data variable:
+            if cube.extended_grid_mapping:
+                # Order the coordinates as per the order in the CRS/WKT string.
+                # (We should only ever have a coordinate system for horizontal
+                # spatial coords, so check for east/north directions)
+                ordered_coords = []
+                for ax_info in cs.as_cartopy_crs().axis_info:
+                    try:
+                        match ax_info.direction.lower():
+                            case "east":
+                                ordered_coords.append(
+                                    cube.coord(axis="X", coord_system=cs)
+                                )
+                            case "north":
+                                ordered_coords.append(
+                                    cube.coord(axis="Y", coord_system=cs)
+                                )
+                            case _:
+                                msg = (
+                                    f"Can't handle axis direction {ax_info.direction!r}"
+                                )
+                                raise NotImplementedError(msg)
+                    except iris.exceptions.CoordinateNotFoundError as e:
                         msg = (
-                            "It is not valid CF to save a true_scale_lat for "
-                            "a Stereographic grid mapping."
+                            f"Failed to assign coordinate for {ax_info.name} axis of "
+                            f"coordinate system {cs.grid_mapping_name}: {str(e)}"
                         )
-                        raise ValueError(msg)
-                    elif cs.scale_factor_at_projection_origin is not None:
-                        cf_var_grid.scale_factor_at_projection_origin = (
-                            cs.scale_factor_at_projection_origin
+                        warnings.warn(msg, category=iris.warnings.IrisSaveWarning)
+
+                        # fall back to simple grid mapping (single crs entry for DimCoord only)
+                        matched_all_coords = False
+
+                # Get variable names, and store them if necessary:
+                cfvar_names = []
+                for coord in ordered_coords:
+                    cfvar = self._name_coord_map.name(coord)
+                    if not cfvar:
+                        # not found - create and store it:
+                        cfvar = self._get_coord_variable_name(cube, coord)
+                        self._name_coord_map.append(
+                            cfvar, self._get_coord_variable_name(cube, coord)
                         )
-                    else:
-                        cf_var_grid.scale_factor_at_projection_origin = 1.0
-
-                # osgb (a specific tmerc)
-                elif isinstance(cs, iris.coord_systems.OSGB):
-                    warnings.warn(
-                        "OSGB coordinate system not yet handled",
-                        category=iris.warnings.IrisSaveWarning,
-                    )
-
-                # lambert azimuthal equal area
-                elif isinstance(cs, iris.coord_systems.LambertAzimuthalEqualArea):
-                    if cs.ellipsoid:
-                        add_ellipsoid(cs.ellipsoid)
-                    cf_var_grid.longitude_of_projection_origin = (
-                        cs.longitude_of_projection_origin
-                    )
-                    cf_var_grid.latitude_of_projection_origin = (
-                        cs.latitude_of_projection_origin
-                    )
-                    cf_var_grid.false_easting = cs.false_easting
-                    cf_var_grid.false_northing = cs.false_northing
-
-                # albers conical equal area
-                elif isinstance(cs, iris.coord_systems.AlbersEqualArea):
-                    if cs.ellipsoid:
-                        add_ellipsoid(cs.ellipsoid)
-                    cf_var_grid.longitude_of_central_meridian = (
-                        cs.longitude_of_central_meridian
-                    )
-                    cf_var_grid.latitude_of_projection_origin = (
-                        cs.latitude_of_projection_origin
-                    )
-                    cf_var_grid.false_easting = cs.false_easting
-                    cf_var_grid.false_northing = cs.false_northing
-                    cf_var_grid.standard_parallel = cs.standard_parallels
-
-                # vertical perspective
-                elif isinstance(cs, iris.coord_systems.VerticalPerspective):
-                    if cs.ellipsoid:
-                        add_ellipsoid(cs.ellipsoid)
-                    cf_var_grid.longitude_of_projection_origin = (
-                        cs.longitude_of_projection_origin
-                    )
-                    cf_var_grid.latitude_of_projection_origin = (
-                        cs.latitude_of_projection_origin
-                    )
-                    cf_var_grid.false_easting = cs.false_easting
-                    cf_var_grid.false_northing = cs.false_northing
-                    cf_var_grid.perspective_point_height = cs.perspective_point_height
-
-                # geostationary
-                elif isinstance(cs, iris.coord_systems.Geostationary):
-                    if cs.ellipsoid:
-                        add_ellipsoid(cs.ellipsoid)
-                    cf_var_grid.longitude_of_projection_origin = (
-                        cs.longitude_of_projection_origin
+                    cfvar_names.append(cfvar)
+
+                coord_string = " ".join(cfvar_names)
+                grid_mappings[cs.grid_mapping_name] = coord_string
+
+        # Refer to grid var
+        if len(coord_systems):
+            grid_mapping = None
+            if cube.extended_grid_mapping:
+                if matched_all_coords:
+                    grid_mapping = " ".join(
+                        f"{cs_name}: {cs_coords}"
+                        for cs_name, cs_coords in grid_mappings.items()
                     )
-                    cf_var_grid.latitude_of_projection_origin = (
-                        cs.latitude_of_projection_origin
-                    )
-                    cf_var_grid.false_easting = cs.false_easting
-                    cf_var_grid.false_northing = cs.false_northing
-                    cf_var_grid.perspective_point_height = cs.perspective_point_height
-                    cf_var_grid.sweep_angle_axis = cs.sweep_angle_axis
-
-                # oblique mercator (and rotated variant)
-                # Use duck-typing over isinstance() - subclasses (i.e.
-                #  RotatedMercator) upset mock tests.
-                elif getattr(cs, "grid_mapping_name", None) == "oblique_mercator":
-                    # RotatedMercator subclasses ObliqueMercator, and RM
-                    #  instances are implicitly saved as OM due to inherited
-                    #  properties. This is correct because CF 1.11 is removing
-                    #  all mention of RM.
-                    if cs.ellipsoid:
-                        add_ellipsoid(cs.ellipsoid)
-                    cf_var_grid.azimuth_of_central_line = cs.azimuth_of_central_line
-                    cf_var_grid.latitude_of_projection_origin = (
-                        cs.latitude_of_projection_origin
-                    )
-                    cf_var_grid.longitude_of_projection_origin = (
-                        cs.longitude_of_projection_origin
-                    )
-                    cf_var_grid.false_easting = cs.false_easting
-                    cf_var_grid.false_northing = cs.false_northing
-                    cf_var_grid.scale_factor_at_projection_origin = (
-                        cs.scale_factor_at_projection_origin
-                    )
-
-                # other
                 else:
-                    warnings.warn(
-                        "Unable to represent the horizontal "
-                        "coordinate system. The coordinate system "
-                        "type %r is not yet implemented." % type(cs),
-                        category=iris.warnings.IrisSaveWarning,
-                    )
-
-                self._coord_systems.append(cs)
+                    # We didn't match all coords required for grid mapping, default to
+                    # single coordinate system for DimCoords, if set.
+                    for ax in ["x", "y"]:
+                        try:
+                            if dim_cs := cube.coord(
+                                axis=ax, dim_coords=True
+                            ).coord_system:
+                                grid_mapping = dim_cs.grid_mapping_name
+                                break
+                        except iris.exceptions.CoordinateNotFoundError:
+                            pass
+            else:
+                # Not using extended_grid_mapping: only ever be on cs in this case:
+                grid_mapping = coord_systems[0].grid_mapping_name
 
-            # Refer to grid var
-            _setncattr(cf_var_cube, "grid_mapping", cs.grid_mapping_name)
+            if grid_mapping:
+                _setncattr(cf_var_cube, "grid_mapping", grid_mapping)
 
     def _create_cf_data_variable(
         self,
@@ -2311,8 +2405,8 @@ class Saver:
             attr_names = set(cube.attributes).intersection(local_keys)
 
         for attr_name in sorted(attr_names):
-            # Do not output 'conventions' attribute.
-            if attr_name.lower() == "conventions":
+            # Do not output 'conventions' or extended grid mapping attribute.
+            if attr_name.lower() in ["conventions", "iris_extended_grid_mapping"]:
                 continue
 
             value = cube.attributes[attr_name]


Test patch:
diff --git a/docs/src/whatsnew/latest.rst b/docs/src/whatsnew/latest.rst
index 592b99b48..55daabbad 100644
--- a/docs/src/whatsnew/latest.rst
+++ b/docs/src/whatsnew/latest.rst
@@ -40,6 +40,31 @@ This document explains the changes made to Iris for this release
    now preserved when it was not previously. See also: :ref:`load-problems`.
    (:pull:`6465`, :pull:`6529`)
 
+#. `@wjbenfold`_ and `@trexfeathers`_ added ``crs_wkt`` to the attributes when
+   saving a :class:`~iris.coord_systems.CoordSystem` to a NetCDF file. Note that
+   ``crs_wkt`` is considered *supplementary* by the CF conventions, with
+   ``grid_mapping`` being the primary source of information, and ``crs_wkt`` not
+   expected to contain conflicting information. Because of this, Iris generates
+   :class:`~iris.coord_systems.CoordSystem` exclusively from ``grid_mapping``
+   when loading, and writes a fresh ``crs_wkt`` whenever a
+   :class:`~iris.coord_systems.CoordSystem` is saved. If your use case goes
+   beyond the CF conventions, you can modify the save and load process for your
+   needs by using the `Ncdata`_ package.
+   See `CRS WKT in the CF Conventions`_ for more. (:issue:`3796`, :pull:`6519`)
+
+#. `@ukmo-ccbunney`_ and `@trexfeathers`_ added support for
+   **multiple coordinate systems** and **ordered coordinates** when loading
+   and saving NetCDF files.
+   This allows for coordinates to be explicitly associated with a coordinate
+   system via an extended syntax in the ``grid_mapping`` attribute of a NetCDF
+   data variable. This extended syntax also supports specification of multiple
+   coordinate systems per data variable. Setting the property
+   ``cube.extended_grid_mapping = True`` will enable extended grid mapping
+   syntax when saving a NetCDF file and also generate an associated **well known
+   text** attribute (``crs_wkt``; as described in :issue:`3796`).
+   See `CRS Grid Mappings and Projections`_ for more information.
+   (:issue:`3388`:, :pull:`6536`:)
+
 #. `@ESadek-MO`_ made MeshCoords immutable. :class:`iris.MeshCoord`s are now updated automatically when
    changing the attached mesh. All changes to the :class:`iris.MeshCoord` should instead be done to
    the relevant :class:`iris.Coord` located on the attached :class:`iris.MeshXY`. This change also affects 
@@ -166,4 +191,7 @@ This document explains the changes made to Iris for this release
 .. comment
     Whatsnew resources in alphabetical order:
 
+.. _CRS WKT in the CF Conventions: https://cfconventions.org/Data/cf-conventions/cf-conventions-1.12/cf-conventions.html#use-of-the-crs-well-known-text-format
+.. _CRS Grid Mappings and Projections: https://cfconventions.org/Data/cf-conventions/cf-conventions-1.12/cf-conventions.html#grid-mappings-and-projections
+.. _Ncdata: https://github.com/pp-mo/ncdata
 .. _Trusted Publishing: https://docs.pypi.org/trusted-publishers/
diff --git a/lib/iris/tests/integration/netcdf/test_coord_systems.py b/lib/iris/tests/integration/netcdf/test_coord_systems.py
index 6c101d902..9f4f272dd 100644
--- a/lib/iris/tests/integration/netcdf/test_coord_systems.py
+++ b/lib/iris/tests/integration/netcdf/test_coord_systems.py
@@ -8,38 +8,23 @@
 # importing anything else.
 import iris.tests as tests  # isort:skip
 
-from os.path import join as path_join
-import shutil
-import tempfile
 import warnings
 
+import numpy as np
 import pytest
 
 import iris
 from iris.coords import DimCoord
 from iris.cube import Cube
 from iris.tests import stock as stock
+from iris.tests._shared_utils import assert_CML
 from iris.tests.stock.netcdf import ncgen_from_cdl
 from iris.tests.unit.fileformats.netcdf.loader import test_load_cubes as tlc
 
 
-@tests.skip_data
-class TestCoordSystem(tests.IrisTest):
-    def setUp(self):
-        tlc.setUpModule()
-
-    def tearDown(self):
-        tlc.tearDownModule()
-
-    def test_load_laea_grid(self):
-        cube = iris.load_cube(
-            tests.get_data_path(
-                ("NetCDF", "lambert_azimuthal_equal_area", "euro_air_temp.nc")
-            )
-        )
-        self.assertCML(cube, ("netcdf", "netcdf_laea.cml"))
-
-    datum_cf_var_cdl = """
+@pytest.fixture
+def datum_cf_var_cdl():
+    return """
         netcdf output {
         dimensions:
             y = 4 ;
@@ -84,7 +69,10 @@ class TestCoordSystem(tests.IrisTest):
         }
     """
 
-    datum_wkt_cdl = """
+
+@pytest.fixture
+def datum_wkt_cdl():
+    return """
 netcdf output5 {
 dimensions:
     y = 4 ;
@@ -93,7 +81,7 @@ variables:
     float data(y, x) ;
         data :standard_name = "toa_brightness_temperature" ;
         data :units = "K" ;
-        data :grid_mapping = "mercator" ;
+        data :grid_mapping = "mercator: x y" ;
     int mercator ;
         mercator:grid_mapping_name = "mercator" ;
         mercator:longitude_of_prime_meridian = 0. ;
@@ -131,25 +119,94 @@ data:
 }
     """
 
-    def test_load_datum_wkt(self):
+
+@pytest.fixture
+def multi_cs_osgb_wkt():
+    return """
+netcdf osgb {
+dimensions:
+    y = 5 ;
+    x = 4 ;
+variables:
+    double x(x) ;
+        x:standard_name = "projection_x_coordinate" ;
+        x:long_name = "Easting" ;
+        x:units = "m" ;
+    double y(y) ;
+        y:standard_name = "projection_y_coordinate" ;
+        y:long_name = "Northing" ;
+        y:units = "m" ;
+    double lat(y, x) ;
+        lat:standard_name = "latitude" ;
+        lat:units = "degrees_north" ;
+    double lon(y, x) ;
+        lon:standard_name = "longitude" ;
+        lon:units = "degrees_east" ;
+    float temp(y, x) ;
+        temp:standard_name = "air_temperature" ;
+        temp:units = "K" ;
+        temp:coordinates = "lat lon" ;
+        temp:grid_mapping = "crsOSGB: x y crsWGS84: lat lon" ;
+    int crsOSGB ;
+        crsOSGB:grid_mapping_name = "transverse_mercator" ;
+        crsOSGB:semi_major_axis = 6377563.396 ;
+        crsOSGB:inverse_flattening = 299.3249646 ;
+        crsOSGB:longitude_of_prime_meridian = 0. ;
+        crsOSGB:latitude_of_projection_origin = 49. ;
+        crsOSGB:longitude_of_central_meridian = -2. ;
+        crsOSGB:scale_factor_at_central_meridian = 0.9996012717 ;
+        crsOSGB:false_easting = 400000. ;
+        crsOSGB:false_northing = -100000. ;
+        crsOSGB:unit = "metre" ;
+        crsOSGB:crs_wkt = "PROJCRS[\\"unknown\\",BASEGEOGCRS[\\"unknown\\",DATUM[\\"Unknown based on Airy 1830 ellipsoid\\",ELLIPSOID[\\"Airy 1830\\",6377563.396,299.324964600004,LENGTHUNIT[\\"metre\\",1,ID[\\"EPSG\\",9001]]]],PRIMEM[\\"Greenwich\\",0,ANGLEUNIT[\\"degree\\",0.0174532925199433],ID[\\"EPSG\\",8901]]],CONVERSION[\\"unknown\\",METHOD[\\"Transverse Mercator\\",ID[\\"EPSG\\",9807]],PARAMETER[\\"Latitude of natural origin\\",49,ANGLEUNIT[\\"degree\\",0.0174532925199433],ID[\\"EPSG\\",8801]],PARAMETER[\\"Longitude of natural origin\\",-2,ANGLEUNIT[\\"degree\\",0.0174532925199433],ID[\\"EPSG\\",8802]],PARAMETER[\\"Scale factor at natural origin\\",0.9996012717,SCALEUNIT[\\"unity\\",1],ID[\\"EPSG\\",8805]],PARAMETER[\\"False easting\\",400000,LENGTHUNIT[\\"metre\\",1],ID[\\"EPSG\\",8806]],PARAMETER[\\"False northing\\",-100000,LENGTHUNIT[\\"metre\\",1],ID[\\"EPSG\\",8807]]],CS[Cartesian,2],AXIS[\\"(E)\\",east,ORDER[1],LENGTHUNIT[\\"metre\\",1,ID[\\"EPSG\\",9001]]],AXIS[\\"(N)\\",north,ORDER[2],LENGTHUNIT[\\"metre\\",1,ID[\\"EPSG\\",9001]]]]" ;
+    int crsWGS84 ;
+        crsWGS84:grid_mapping_name = "latitude_longitude" ;
+        crsWGS84:longitude_of_prime_meridian = 0. ;
+        crsWGS84:semi_major_axis = 6378137. ;
+        crsWGS84:inverse_flattening = 298.257223563 ;
+        crsWGS84: crs_wkt = "GEOGCRS[\\"unknown\\",DATUM[\\"Unknown based on WGS 84 ellipsoid\\",ELLIPSOID[\\"WGS 84\\",6378137,298.257223562997,LENGTHUNIT[\\"metre\\",1,ID[\\"EPSG\\",9001]]]],PRIMEM[\\"Greenwich\\",0,ANGLEUNIT[\\"degree\\",0.0174532925199433],ID[\\"EPSG\\",8901]],CS[ellipsoidal,2],AXIS[\\"longitude\\",east,ORDER[1],ANGLEUNIT[\\"degree\\",0.0174532925199433,ID[\\"EPSG\\",9122]]],AXIS[\\"latitude\\",north,ORDER[2],ANGLEUNIT[\\"degree\\",0.0174532925199433,ID[\\"EPSG\\",9122]]]]" ;
+data:
+    x = 1,2,3,4 ;
+    y = 1,2,3,4,5 ;
+}
+    """
+
+
+@tests.skip_data
+class TestCoordSystem:
+    @pytest.fixture(autouse=True)
+    def _setup(self):
+        tlc.setUpModule()
+        yield
+        tlc.tearDownModule()
+
+    def test_load_laea_grid(self, request):
+        cube = iris.load_cube(
+            tests.get_data_path(
+                ("NetCDF", "lambert_azimuthal_equal_area", "euro_air_temp.nc")
+            )
+        )
+        assert_CML(request, cube, ("netcdf", "netcdf_laea.cml"))
+
+    def test_load_datum_wkt(self, datum_wkt_cdl):
         expected = "OSGB 1936"
-        nc_path = tlc.cdl_to_nc(self.datum_wkt_cdl)
+        nc_path = tlc.cdl_to_nc(datum_wkt_cdl)
         with iris.FUTURE.context(datum_support=True):
             cube = iris.load_cube(nc_path)
         test_crs = cube.coord("projection_y_coordinate").coord_system
         actual = str(test_crs.as_cartopy_crs().datum)
         assert actual == expected
 
-    def test_no_load_datum_wkt(self):
-        nc_path = tlc.cdl_to_nc(self.datum_wkt_cdl)
+    def test_no_load_datum_wkt(self, datum_wkt_cdl):
+        nc_path = tlc.cdl_to_nc(datum_wkt_cdl)
         with pytest.warns(FutureWarning, match="iris.FUTURE.datum_support"):
             cube = iris.load_cube(nc_path)
         test_crs = cube.coord("projection_y_coordinate").coord_system
         actual = str(test_crs.as_cartopy_crs().datum)
         assert actual == "unknown"
 
-    def test_no_datum_no_warn(self):
-        new_cdl = self.datum_wkt_cdl.splitlines()
+    def test_no_datum_no_warn(self, datum_wkt_cdl):
+        new_cdl = datum_wkt_cdl.splitlines()
         new_cdl = [line for line in new_cdl if "DATUM" not in line]
         new_cdl = "\n".join(new_cdl)
         nc_path = tlc.cdl_to_nc(new_cdl)
@@ -158,24 +215,40 @@ data:
             warnings.simplefilter("error", FutureWarning)
             _ = iris.load_cube(nc_path)
 
-    def test_load_datum_cf_var(self):
+    def test_load_datum_cf_var(self, datum_cf_var_cdl):
         expected = "OSGB 1936"
-        nc_path = tlc.cdl_to_nc(self.datum_cf_var_cdl)
+        nc_path = tlc.cdl_to_nc(datum_cf_var_cdl)
         with iris.FUTURE.context(datum_support=True):
             cube = iris.load_cube(nc_path)
         test_crs = cube.coord("projection_y_coordinate").coord_system
         actual = str(test_crs.as_cartopy_crs().datum)
         assert actual == expected
 
-    def test_no_load_datum_cf_var(self):
-        nc_path = tlc.cdl_to_nc(self.datum_cf_var_cdl)
+    def test_no_load_datum_cf_var(self, datum_cf_var_cdl):
+        nc_path = tlc.cdl_to_nc(datum_cf_var_cdl)
         with pytest.warns(FutureWarning, match="iris.FUTURE.datum_support"):
             cube = iris.load_cube(nc_path)
         test_crs = cube.coord("projection_y_coordinate").coord_system
         actual = str(test_crs.as_cartopy_crs().datum)
         assert actual == "unknown"
 
-    def test_save_datum(self):
+    def test_load_multi_cs_wkt(self, multi_cs_osgb_wkt):
+        nc_path = tlc.cdl_to_nc(multi_cs_osgb_wkt)
+        with iris.FUTURE.context(datum_support=True):
+            cube = iris.load_cube(nc_path)
+
+        assert len(cube.coord_systems()) == 2
+        for name in ["projection_y_coordinate", "projection_y_coordinate"]:
+            assert (
+                cube.coord(name).coord_system.grid_mapping_name == "transverse_mercator"
+            )
+        for name in ["latitude", "longitude"]:
+            assert (
+                cube.coord(name).coord_system.grid_mapping_name == "latitude_longitude"
+            )
+        assert cube.extended_grid_mapping is True
+
+    def test_save_datum(self, tmp_path):
         expected = "OSGB 1936"
         saved_crs = iris.coord_systems.Mercator(
             ellipsoid=iris.coord_systems.GeogCS.from_datum("OSGB36")
@@ -204,24 +277,73 @@ data:
                 (test_lon_coord, 2),
             ),
         )
-
-        with self.temp_filename(suffix=".nc") as filename:
-            iris.save(test_cube, filename)
-            with iris.FUTURE.context(datum_support=True):
-                cube = iris.load_cube(filename)
+        filename = tmp_path / "output.nc"
+        iris.save(test_cube, filename)
+        with iris.FUTURE.context(datum_support=True):
+            cube = iris.load_cube(filename)
 
         test_crs = cube.coord("projection_y_coordinate").coord_system
         actual = str(test_crs.as_cartopy_crs().datum)
         assert actual == expected
 
+    def test_save_multi_cs_wkt(self, tmp_path):
+        crsOSGB = iris.coord_systems.OSGB()
+        crsLatLon = iris.coord_systems.GeogCS(6e6)
 
-class TestLoadMinimalGeostationary(tests.IrisTest):
-    """Check we can load data with a geostationary grid-mapping, even when the
-    'false-easting' and 'false_northing' properties are missing.
+        dimx_coord = iris.coords.DimCoord(
+            np.arange(4), "projection_x_coordinate", coord_system=crsOSGB
+        )
+        dimy_coord = iris.coords.DimCoord(
+            np.arange(5), "projection_y_coordinate", coord_system=crsOSGB
+        )
 
-    """
+        auxlon_coord = iris.coords.AuxCoord(
+            np.arange(20).reshape((5, 4)),
+            standard_name="longitude",
+            coord_system=crsLatLon,
+        )
+        auxlat_coord = iris.coords.AuxCoord(
+            np.arange(20).reshape((5, 4)),
+            standard_name="latitude",
+            coord_system=crsLatLon,
+        )
+
+        test_cube = Cube(
+            np.ones(20).reshape((5, 4)),
+            standard_name="air_pressure",
+            units="Pa",
+            dim_coords_and_dims=(
+                (dimy_coord, 0),
+                (dimx_coord, 1),
+            ),
+            aux_coords_and_dims=(
+                (auxlat_coord, (0, 1)),
+                (auxlon_coord, (0, 1)),
+            ),
+        )
 
-    _geostationary_problem_cdl = """
+        test_cube.extended_grid_mapping = True
+
+        filename = tmp_path / "output.nc"
+        iris.save(test_cube, filename)
+        with iris.FUTURE.context(datum_support=True):
+            cube = iris.load_cube(filename)
+
+        assert len(cube.coord_systems()) == 2
+        for name in ["projection_y_coordinate", "projection_y_coordinate"]:
+            assert (
+                cube.coord(name).coord_system.grid_mapping_name == "transverse_mercator"
+            )
+        for name in ["latitude", "longitude"]:
+            assert (
+                cube.coord(name).coord_system.grid_mapping_name == "latitude_longitude"
+            )
+        assert cube.extended_grid_mapping is True
+
+
+@pytest.fixture(scope="module")
+def geostationary_problem_cdl():
+    return """
 netcdf geostationary_problem_case {
 dimensions:
     y = 2 ;
@@ -258,35 +380,34 @@ data:
  x = 0, 1, 2 ;
 
 }
-"""
-
-    @classmethod
-    def setUpClass(cls):
-        # Create a temp directory for transient test files.
-        cls.temp_dir = tempfile.mkdtemp()
-        cls.path_test_cdl = path_join(cls.temp_dir, "geos_problem.cdl")
-        cls.path_test_nc = path_join(cls.temp_dir, "geos_problem.nc")
-        # Create reference CDL and netcdf files from the CDL text.
+    """
+
+
+class TestLoadMinimalGeostationary:
+    """Check we can load data with a geostationary grid-mapping, even when the
+    'false-easting' and 'false_northing' properties are missing.
+
+    """
+
+    @pytest.fixture(scope="class")
+    def geostationary_problem_ncfile(self, tmp_path_factory, geostationary_problem_cdl):
+        tmp_path = tmp_path_factory.mktemp("geos")
+        cdl_path = tmp_path / "geos_problem.cdl"
+        nc_path = tmp_path / "geos_problem.nc"
         ncgen_from_cdl(
-            cdl_str=cls._geostationary_problem_cdl,
-            cdl_path=cls.path_test_cdl,
-            nc_path=cls.path_test_nc,
+            cdl_str=geostationary_problem_cdl,
+            cdl_path=cdl_path,
+            nc_path=nc_path,
         )
+        return nc_path
 
-    @classmethod
-    def tearDownClass(cls):
-        # Destroy the temp directory.
-        shutil.rmtree(cls.temp_dir)
-
-    def test_geostationary_no_false_offsets(self):
+    def test_geostationary_no_false_offsets(
+        self, tmp_path, geostationary_problem_ncfile
+    ):
         # Check we can load the test data and coordinate system properties are correct.
-        cube = iris.load_cube(self.path_test_nc)
+        cube = iris.load_cube(geostationary_problem_ncfile)
         # Check the coordinate system properties has the correct default properties.
         cs = cube.coord_system()
         assert isinstance(cs, iris.coord_systems.Geostationary)
         assert cs.false_easting == 0.0
         assert cs.false_northing == 0.0
-
-
-if __name__ == "__main__":
-    tests.main()
diff --git a/lib/iris/tests/integration/test_netcdf__loadsaveattrs.py b/lib/iris/tests/integration/test_netcdf__loadsaveattrs.py
index b89477bfe..a2d6e9c75 100644
--- a/lib/iris/tests/integration/test_netcdf__loadsaveattrs.py
+++ b/lib/iris/tests/integration/test_netcdf__loadsaveattrs.py
@@ -69,6 +69,9 @@ _LOCAL_TEST_ATTRS = (
     iris.fileformats.netcdf.saver._CF_DATA_ATTRS
     + iris.fileformats.netcdf.saver._UKMO_DATA_ATTRS
 )
+# Don't test iris_extended_grid_mapping, as it is a special attribute that
+# is not expected to always roundtrip.
+_LOCAL_TEST_ATTRS = [a for a in _LOCAL_TEST_ATTRS if a != "iris_extended_grid_mapping"]
 
 
 # Define a fixture to parametrise over the 'local-style' test attributes.
@@ -513,6 +516,7 @@ _SPECIAL_ATTRS = [
     "standard_error_multiplier",
     "STASH",
     "um_stash_source",
+    "iris_extended_grid_mapping",
 ]
 _MATRIX_ATTRNAMES = [attr for attr in _MATRIX_ATTRNAMES if attr not in _SPECIAL_ATTRS]
 
diff --git a/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/multi_cs.cdl b/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/multi_cs.cdl
new file mode 100644
index 000000000..547813fdf
--- /dev/null
+++ b/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/multi_cs.cdl
@@ -0,0 +1,40 @@
+dimensions:
+	projection_x_coordinate = 4 ;
+	projection_y_coordinate = 3 ;
+variables:
+	int64 air_pressure_anomaly(projection_y_coordinate, projection_x_coordinate) ;
+		air_pressure_anomaly:standard_name = "air_pressure_anomaly" ;
+		air_pressure_anomaly:grid_mapping = "transverse_mercator: projection_x_coordinate projection_y_coordinate latitude_longitude: longitude latitude" ;
+		air_pressure_anomaly:coordinates = "latitude longitude" ;
+	int transverse_mercator ;
+		transverse_mercator:grid_mapping_name = "transverse_mercator" ;
+		transverse_mercator:longitude_of_prime_meridian = 0. ;
+		transverse_mercator:semi_major_axis = 6377563.396 ;
+		transverse_mercator:semi_minor_axis = 6356256.909 ;
+		transverse_mercator:longitude_of_central_meridian = -2. ;
+		transverse_mercator:latitude_of_projection_origin = 49. ;
+		transverse_mercator:false_easting = -400000. ;
+		transverse_mercator:false_northing = 100000. ;
+		transverse_mercator:scale_factor_at_central_meridian = 0.9996012717 ;
+		transverse_mercator:crs_wkt = "PROJCRS[\"unknown\",BASEGEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6377563.396,299.324961266495,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]]],CONVERSION[\"unknown\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",49,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-2,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996012717,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",-400000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",100000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]]" ;
+	int latitude_longitude ;
+		latitude_longitude:grid_mapping_name = "latitude_longitude" ;
+		latitude_longitude:longitude_of_prime_meridian = 0. ;
+		latitude_longitude:semi_major_axis = 6377563.396 ;
+		latitude_longitude:semi_minor_axis = 6356256.909 ;
+		latitude_longitude:crs_wkt = "GEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6377563.396,299.324961266495,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]],CS[ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]]" ;
+	int64 projection_y_coordinate(projection_y_coordinate) ;
+		projection_y_coordinate:axis = "Y" ;
+		projection_y_coordinate:units = "m" ;
+		projection_y_coordinate:standard_name = "projection_y_coordinate" ;
+	int64 projection_x_coordinate(projection_x_coordinate) ;
+		projection_x_coordinate:axis = "X" ;
+		projection_x_coordinate:units = "m" ;
+		projection_x_coordinate:standard_name = "projection_x_coordinate" ;
+	int64 latitude(projection_y_coordinate, projection_x_coordinate) ;
+		latitude:units = "degrees_north" ;
+		latitude:standard_name = "latitude" ;
+	int64 longitude(projection_y_coordinate, projection_x_coordinate) ;
+		longitude:units = "degrees_east" ;
+		longitude:standard_name = "longitude" ;
+}
diff --git a/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/multi_cs_missing_coord.cdl b/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/multi_cs_missing_coord.cdl
new file mode 100644
index 000000000..34c7061c9
--- /dev/null
+++ b/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/multi_cs_missing_coord.cdl
@@ -0,0 +1,37 @@
+dimensions:
+	projection_x_coordinate = 4 ;
+	projection_y_coordinate = 3 ;
+variables:
+	int64 air_pressure_anomaly(projection_y_coordinate, projection_x_coordinate) ;
+		air_pressure_anomaly:standard_name = "air_pressure_anomaly" ;
+		air_pressure_anomaly:grid_mapping = "transverse_mercator" ;
+		air_pressure_anomaly:coordinates = "longitude" ;
+	int transverse_mercator ;
+		transverse_mercator:grid_mapping_name = "transverse_mercator" ;
+		transverse_mercator:longitude_of_prime_meridian = 0. ;
+		transverse_mercator:semi_major_axis = 6377563.396 ;
+		transverse_mercator:semi_minor_axis = 6356256.909 ;
+		transverse_mercator:longitude_of_central_meridian = -2. ;
+		transverse_mercator:latitude_of_projection_origin = 49. ;
+		transverse_mercator:false_easting = -400000. ;
+		transverse_mercator:false_northing = 100000. ;
+		transverse_mercator:scale_factor_at_central_meridian = 0.9996012717 ;
+		transverse_mercator:crs_wkt = "PROJCRS[\"unknown\",BASEGEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6377563.396,299.324961266495,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]]],CONVERSION[\"unknown\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",49,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-2,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996012717,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",-400000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",100000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]]" ;
+	int latitude_longitude ;
+		latitude_longitude:grid_mapping_name = "latitude_longitude" ;
+		latitude_longitude:longitude_of_prime_meridian = 0. ;
+		latitude_longitude:semi_major_axis = 6377563.396 ;
+		latitude_longitude:semi_minor_axis = 6356256.909 ;
+		latitude_longitude:crs_wkt = "GEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6377563.396,299.324961266495,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]],CS[ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]]" ;
+	int64 projection_y_coordinate(projection_y_coordinate) ;
+		projection_y_coordinate:axis = "Y" ;
+		projection_y_coordinate:units = "m" ;
+		projection_y_coordinate:standard_name = "projection_y_coordinate" ;
+	int64 projection_x_coordinate(projection_x_coordinate) ;
+		projection_x_coordinate:axis = "X" ;
+		projection_x_coordinate:units = "m" ;
+		projection_x_coordinate:standard_name = "projection_x_coordinate" ;
+	int64 longitude(projection_y_coordinate, projection_x_coordinate) ;
+		longitude:units = "degrees_east" ;
+		longitude:standard_name = "longitude" ;
+}
diff --git a/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/no_aux_cs.cdl b/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/no_aux_cs.cdl
new file mode 100644
index 000000000..9b5aeb2df
--- /dev/null
+++ b/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/no_aux_cs.cdl
@@ -0,0 +1,34 @@
+dimensions:
+	projection_x_coordinate = 4 ;
+	projection_y_coordinate = 3 ;
+variables:
+	int64 air_pressure_anomaly(projection_y_coordinate, projection_x_coordinate) ;
+		air_pressure_anomaly:standard_name = "air_pressure_anomaly" ;
+		air_pressure_anomaly:grid_mapping = "transverse_mercator: projection_x_coordinate projection_y_coordinate" ;
+		air_pressure_anomaly:coordinates = "latitude longitude" ;
+	int transverse_mercator ;
+		transverse_mercator:grid_mapping_name = "transverse_mercator" ;
+		transverse_mercator:longitude_of_prime_meridian = 0. ;
+		transverse_mercator:semi_major_axis = 6377563.396 ;
+		transverse_mercator:semi_minor_axis = 6356256.909 ;
+		transverse_mercator:longitude_of_central_meridian = -2. ;
+		transverse_mercator:latitude_of_projection_origin = 49. ;
+		transverse_mercator:false_easting = -400000. ;
+		transverse_mercator:false_northing = 100000. ;
+		transverse_mercator:scale_factor_at_central_meridian = 0.9996012717 ;
+		transverse_mercator:crs_wkt = "PROJCRS[\"unknown\",BASEGEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6377563.396,299.324961266495,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]]],CONVERSION[\"unknown\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",49,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-2,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996012717,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",-400000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",100000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]]" ;
+	int64 projection_y_coordinate(projection_y_coordinate) ;
+		projection_y_coordinate:axis = "Y" ;
+		projection_y_coordinate:units = "m" ;
+		projection_y_coordinate:standard_name = "projection_y_coordinate" ;
+	int64 projection_x_coordinate(projection_x_coordinate) ;
+		projection_x_coordinate:axis = "X" ;
+		projection_x_coordinate:units = "m" ;
+		projection_x_coordinate:standard_name = "projection_x_coordinate" ;
+	int64 latitude(projection_y_coordinate, projection_x_coordinate) ;
+		latitude:units = "degrees_north" ;
+		latitude:standard_name = "latitude" ;
+	int64 longitude(projection_y_coordinate, projection_x_coordinate) ;
+		longitude:units = "degrees_east" ;
+		longitude:standard_name = "longitude" ;
+}
diff --git a/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/no_cs.cdl b/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/no_cs.cdl
new file mode 100644
index 000000000..ce8825b16
--- /dev/null
+++ b/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/no_cs.cdl
@@ -0,0 +1,22 @@
+dimensions:
+	projection_x_coordinate = 4 ;
+	projection_y_coordinate = 3 ;
+variables:
+	int64 air_pressure_anomaly(projection_y_coordinate, projection_x_coordinate) ;
+		air_pressure_anomaly:standard_name = "air_pressure_anomaly" ;
+		air_pressure_anomaly:coordinates = "latitude longitude" ;
+	int64 projection_y_coordinate(projection_y_coordinate) ;
+		projection_y_coordinate:axis = "Y" ;
+		projection_y_coordinate:units = "m" ;
+		projection_y_coordinate:standard_name = "projection_y_coordinate" ;
+	int64 projection_x_coordinate(projection_x_coordinate) ;
+		projection_x_coordinate:axis = "X" ;
+		projection_x_coordinate:units = "m" ;
+		projection_x_coordinate:standard_name = "projection_x_coordinate" ;
+	int64 latitude(projection_y_coordinate, projection_x_coordinate) ;
+		latitude:units = "degrees_north" ;
+		latitude:standard_name = "latitude" ;
+	int64 longitude(projection_y_coordinate, projection_x_coordinate) ;
+		longitude:units = "degrees_east" ;
+		longitude:standard_name = "longitude" ;
+}
diff --git a/lib/iris/tests/unit/coord_systems/test_ObliqueMercator.py b/lib/iris/tests/unit/coord_systems/test_ObliqueMercator.py
index d80b4311c..95069ba37 100644
--- a/lib/iris/tests/unit/coord_systems/test_ObliqueMercator.py
+++ b/lib/iris/tests/unit/coord_systems/test_ObliqueMercator.py
@@ -5,7 +5,6 @@
 """Unit tests for the :class:`iris.coord_systems.ObliqueMercator` class."""
 
 from typing import List, NamedTuple
-from unittest.mock import Mock
 
 from cartopy import crs as ccrs
 import pytest
@@ -130,18 +129,20 @@ class TestArgs:
     def instance(self):
         return self.make_instance()
 
+    @pytest.fixture()
+    def mock_ccrs(self, mocker):
+        return mocker.patch("cartopy.crs.ObliqueMercator", autospec=True)
+
     def test_instantiate(self):
         _ = self.make_instance()
 
-    def test_cartopy_crs(self, instance):
-        ccrs.ObliqueMercator = Mock()
+    def test_cartopy_crs(self, instance, mock_ccrs):
         instance.as_cartopy_crs()
-        ccrs.ObliqueMercator.assert_called_with(**self.cartopy_kwargs_expected)
+        mock_ccrs.assert_called_with(**self.cartopy_kwargs_expected)
 
-    def test_cartopy_projection(self, instance):
-        ccrs.ObliqueMercator = Mock()
+    def test_cartopy_projection(self, instance, mock_ccrs):
         instance.as_cartopy_projection()
-        ccrs.ObliqueMercator.assert_called_with(**self.cartopy_kwargs_expected)
+        mock_ccrs.assert_called_with(**self.cartopy_kwargs_expected)
 
     @pytest.fixture()
     def label_class(self, instance):
diff --git a/lib/iris/tests/unit/fileformats/nc_load_rules/actions/__init__.py b/lib/iris/tests/unit/fileformats/nc_load_rules/actions/__init__.py
index 7c7974002..296765f85 100644
--- a/lib/iris/tests/unit/fileformats/nc_load_rules/actions/__init__.py
+++ b/lib/iris/tests/unit/fileformats/nc_load_rules/actions/__init__.py
@@ -65,7 +65,7 @@ class Mixin__nc_load_actions:
         # Destroy a temp directory for temp files.
         shutil.rmtree(cls.temp_dirpath)
 
-    def load_cube_from_cdl(self, cdl_string, cdl_path, nc_path):
+    def load_cube_from_cdl(self, cdl_string, cdl_path, nc_path, mocker=None):
         """Load the 'phenom' data variable in a CDL testcase, as a cube.
 
         Using ncgen, CFReader and the _load_cube call.
@@ -85,8 +85,12 @@ class Mixin__nc_load_actions:
 
             # If debug enabled, switch on the activation summary debug output.
             # Use 'patch' so it is restored after the test.
-            self.patch("iris.fileformats.netcdf.loader.DEBUG", self.debug_info)
-
+            if mocker:
+                # If pytest mocker provided, use that to patch:
+                mocker.patch("iris.fileformats.netcdf.loader.DEBUG", self.debug_info)
+            elif hasattr(self, "patch"):
+                # keep old UnitTest patch working
+                self.patch("iris.fileformats.netcdf.loader.DEBUG", self.debug_info)
             with warnings.catch_warnings():
                 warnings.filterwarnings(
                     "ignore",
diff --git a/lib/iris/tests/unit/fileformats/nc_load_rules/actions/test__grid_mappings.py b/lib/iris/tests/unit/fileformats/nc_load_rules/actions/test__grid_mappings.py
index 8731049ea..202fb0fa1 100644
--- a/lib/iris/tests/unit/fileformats/nc_load_rules/actions/test__grid_mappings.py
+++ b/lib/iris/tests/unit/fileformats/nc_load_rules/actions/test__grid_mappings.py
@@ -9,8 +9,12 @@ Here, *specifically* testcases relating to grid-mappings and dim-coords.
 
 """
 
+import iris.coord_systems
+
 import iris.tests as tests  # isort: skip
 
+import pytest
+
 import iris.coord_systems as ics
 import iris.fileformats._nc_load_rules.helpers as hh
 from iris.loading import LOAD_PROBLEMS
@@ -27,6 +31,7 @@ class Mixin__grid_mapping(Mixin__nc_load_actions):
         mapping_missingradius=False,
         mapping_type_name=None,
         mapping_scalefactor=None,
+        phenom_grid_mapping="grid",
         yco_values=None,
         xco_name=None,
         yco_name=None,
@@ -226,7 +231,7 @@ class Mixin__grid_mapping(Mixin__nc_load_actions):
                 double phenom({ydim_name}, {xdim_name}) ;
                     phenom:standard_name = "air_temperature" ;
                     phenom:units = "K" ;
-                    phenom:grid_mapping = "grid" ;
+                    phenom:grid_mapping = "{phenom_grid_mapping}" ;
 {phenom_coords_string}
                 double yco({ydim_name}) ;
                     yco:axis = "Y" ;
@@ -325,6 +330,7 @@ class Mixin__grid_mapping(Mixin__nc_load_actions):
             self.assertIsNone(yco_cs)
             self.assertIsNone(xco_cs)
         else:
+            self.assertIsNotNone(cube_cs)
             if cube_cstype is not None:
                 self.assertIsInstance(cube_cs, cube_cstype)
             if xco_no_cs:
@@ -755,6 +761,45 @@ class Test__grid_mapping(Mixin__grid_mapping, tests.IrisTest):
         )
         self.check_result(result, cube_no_cs=True, xco_stdname=False, yco_stdname=False)
 
+    def test_extended_mapping_basic_latlon(self):
+        # A basic reference example with a lat-long grid, but using extended
+        # grid mapping syntax on phenomenon.
+        #
+        # Rules Triggered:
+        #     001 : fc_default
+        #     002 : fc_provides_grid_mapping_(latitude_longitude)
+        #     003 : fc_provides_coordinate_(latitude)
+        #     004 : fc_provides_coordinate_(longitude)
+        #     005 : fc_build_coordinate_(latitude)
+        #     006 : fc_build_coordinate_(longitude)
+        # Notes:
+        #     * grid-mapping identified : regular latlon
+        #     * dim-coords identified : lat+lon
+        #     * coords built : standard latlon (with latlon coord-system)
+        result = self.run_testcase(phenom_grid_mapping="grid: yco xco")
+        self.check_result(result)
+
+    def test_extended_mapping_basic_latlon_missing_coords(self):
+        # A basic reference example with a lat-long grid, but using extended
+        # grid mapping syntax on phenomenon.
+        #
+        # Rules Triggered:
+        #     001 : fc_default
+        #     002 : fc_provides_grid_mapping_(latitude_longitude)
+        #     003 : fc_provides_coordinate_(latitude)
+        #     004 : fc_provides_coordinate_(longitude)
+        #     005 : fc_build_coordinate_(latitude)
+        #     006 : fc_build_coordinate_(longitude)
+        # Notes:
+        #     * grid-mapping identified : regular latlon
+        #     * dim-coords identified : lat+lon
+        #     * coords built : standard latlon (with latlon coord-system)
+        result = self.run_testcase(
+            phenom_grid_mapping="grid: yco bad_coord",
+            warning_regex="Missing CF-netCDF coordinate variable 'bad_coord'",
+        )
+        self.check_result(result, xco_no_cs=True)
+
 
 class Test__aux_latlons(Mixin__grid_mapping, tests.IrisTest):
     # Testcases for translating auxiliary latitude+longitude variables
@@ -838,6 +883,55 @@ class Test__aux_latlons(Mixin__grid_mapping, tests.IrisTest):
         )
         self.check_result(result, yco_is_aux=True, yco_no_cs=True)
 
+    def test_extended_grid_mapping_aux_lon(self):
+        # Change the name of xdim, and put xco on the coords list.
+        # Uses extended grid mapping syntax.
+        # In this case, the Aux coord WILL have a coordinate system
+        # as extended grid mapping allows for specification of
+        # explicit coordinate_systems for individual coordinate.
+        #
+        # Rules Triggered:
+        #     001 : fc_default
+        #     002 : fc_provides_grid_mapping_(latitude_longitude)
+        #     003 : fc_provides_coordinate_(latitude)
+        #     004 : fc_build_coordinate_(latitude)
+        #     005 : fc_build_auxiliary_coordinate_longitude
+        result = self.run_testcase(
+            xco_is_dim=False, phenom_grid_mapping="grid: yco xco"
+        )
+        self.check_result(result, xco_is_aux=True, xco_no_cs=False)
+
+    def test_extended_grid_mapping_aux_lat(self):
+        # As previous, but with the Y coord.
+        # Uses extended grid mapping syntax
+        #
+        # Rules Triggered:
+        #     001 : fc_default
+        #     002 : fc_provides_grid_mapping_(latitude_longitude)
+        #     003 : fc_provides_coordinate_(longitude)
+        #     004 : fc_build_coordinate_(longitude)
+        #     005 : fc_build_auxiliary_coordinate_latitude
+        result = self.run_testcase(
+            yco_is_dim=False, phenom_grid_mapping="grid: yco xco"
+        )
+        self.check_result(result, yco_is_aux=True, yco_no_cs=False)
+
+    def test_extended_grid_mapping_aux_lat_and_lon(self):
+        # Make *both* X and Y coords into aux-coords.
+        # Uses extended grid mapping syntax; allows coord_system
+        # to be added to an AuxCoord, so cube.coord_system() returns
+        # a valid coordinate_system.
+        #
+        # Rules Triggered:
+        #     001 : fc_default
+        #     002 : fc_provides_grid_mapping_(latitude_longitude)
+        #     003 : fc_build_auxiliary_coordinate_longitude
+        #     004 : fc_build_auxiliary_coordinate_latitude
+        result = self.run_testcase(
+            xco_is_dim=False, yco_is_dim=False, phenom_grid_mapping="grid: yco xco"
+        )
+        self.check_result(result, xco_is_aux=True, yco_is_aux=True, cube_no_cs=False)
+
 
 class Test__nondimcoords(Mixin__grid_mapping, tests.IrisTest):
     @classmethod
@@ -869,5 +963,279 @@ class Test__nondimcoords(Mixin__grid_mapping, tests.IrisTest):
         self.check_result(result, yco_is_aux=True, load_problems_regex=error)
 
 
+@pytest.fixture
+def latlon_cs():
+    return """
+        int crsWGS84 ;
+            crsWGS84:grid_mapping_name = "latitude_longitude" ;
+            crsWGS84:longitude_of_prime_meridian = 0. ;
+            crsWGS84:semi_major_axis = 6378137. ;
+            crsWGS84:inverse_flattening = 298.257223563 ;
+        double lat(x,y) ;
+            lat:standard_name = "latitude" ;
+            lat:units = "degrees_north" ;
+        double lon(x,y) ;
+            lon:standard_name = "longitude" ;
+            lon:units = "degrees_east" ;
+    """
+
+
+@pytest.fixture
+def latlon_cs_missing_coord(latlon_cs):
+    # Trin last 3 lines of latlon_cs to remove longitude coord
+    return "\n".join(latlon_cs.splitlines()[:-4])
+
+
+@pytest.fixture
+def osgb_cs():
+    return """
+        int crsOSGB ;
+            crsOSGB:grid_mapping_name = "transverse_mercator" ;
+            crsOSGB:semi_major_axis = 6377563.396 ;
+            crsOSGB:inverse_flattening = 299.3249646 ;
+            crsOSGB:longitude_of_prime_meridian = 0. ;
+            crsOSGB:latitude_of_projection_origin = 49. ;
+            crsOSGB:longitude_of_central_meridian = -2. ;
+            crsOSGB:scale_factor_at_central_meridian = 0.9996012717 ;
+            crsOSGB:false_easting = 400000. ;
+            crsOSGB:false_northing = -100000. ;
+            crsOSGB:unit = "metre" ;
+        double x(x) ;
+            x:standard_name = "projection_x_coordinate" ;
+            x:long_name = "Easting" ;
+            x:units = "m" ;
+        double y(y) ;
+            y:standard_name = "projection_y_coordinate" ;
+            y:long_name = "Northing" ;
+            y:units = "m" ;
+    """
+
+
+class Test_multi_coordinate_system_grid_mapping(Mixin__nc_load_actions):
+    def test_two_coord_systems(self, osgb_cs, latlon_cs, mocker, tmp_path):
+        """Test load a well described multi coordinate system variable."""
+        cdl = f"""
+        netcdf tmp {{
+        dimensions:
+            x = 4 ;
+            y = 3 ;
+        variables:
+            float phenom(y, x) ;
+                phenom:standard_name = "air_pressure" ;
+                phenom:units = "Pa" ;
+                phenom:coordinates = "lat lon" ;
+                phenom:grid_mapping = "crsOSGB: x y crsWGS84: lat lon" ;
+            {osgb_cs}
+            {latlon_cs}
+        }}
+        """
+        nc_path = str(tmp_path / "tmp.nc")
+        cube = self.load_cube_from_cdl(cdl, None, nc_path, mocker)
+
+        assert len(cube.coord_systems()) == 2
+
+        for coord_name in ["projection_x_coordinate", "projection_x_coordinate"]:
+            assert isinstance(
+                cube.coord(coord_name).coord_system,
+                iris.coord_systems.TransverseMercator,
+            )
+        for coord_name in ["longitude", "latitude"]:
+            assert isinstance(
+                cube.coord(coord_name).coord_system, iris.coord_systems.GeogCS
+            )
+
+        # Loading multiple coord systems or using extended grid mapping implies ordered axes:
+        assert cube.extended_grid_mapping is True
+
+    def test_two_coord_systems_missing_coord(
+        self, osgb_cs, latlon_cs_missing_coord, mocker, tmp_path
+    ):
+        """Test missing coord in grid_mapping raises warning."""
+        cdl = f"""
+        netcdf tmp {{
+        dimensions:
+            x = 4 ;
+            y = 3 ;
+        variables:
+            float phenom(y, x) ;
+                phenom:standard_name = "air_pressure" ;
+                phenom:units = "Pa" ;
+                phenom:coordinates = "lat" ;
+                phenom:grid_mapping = "crsOSGB: x y crsWGS84: lat lon" ;
+            {osgb_cs}
+            {latlon_cs_missing_coord}
+        }}
+        """
+        nc_path = str(tmp_path / "tmp.nc")
+
+        # loader will warn that it can't find the longitude coord
+        with pytest.warns(
+            iris.warnings.IrisCfWarning,
+            match="Missing CF-netCDF coordinate variable 'lon'",
+        ):
+            cube = self.load_cube_from_cdl(cdl, None, nc_path, mocker)
+
+        assert len(cube.coords()) == 3
+        assert len(cube.coord_systems()) == 2
+
+        for coord_name in ["projection_x_coordinate", "projection_x_coordinate"]:
+            assert isinstance(
+                cube.coord(coord_name).coord_system,
+                iris.coord_systems.TransverseMercator,
+            )
+        for coord_name in ["latitude"]:
+            assert isinstance(
+                cube.coord(coord_name).coord_system, iris.coord_systems.GeogCS
+            )
+
+        assert cube.extended_grid_mapping is True
+
+    def test_two_coord_systems_missing_aux_crs(
+        self, osgb_cs, latlon_cs, mocker, tmp_path
+    ):
+        """Test invalid coordinate system mapping for Aux coords."""
+        cdl = f"""
+        netcdf tmp {{
+        dimensions:
+            x = 4 ;
+            y = 3 ;
+        variables:
+            float phenom(y, x) ;
+                phenom:standard_name = "air_pressure" ;
+                phenom:units = "Pa" ;
+                phenom:coordinates = "lat lon" ;
+                phenom:grid_mapping = "crsOSGB: x y non_existent_crs: lat lon" ;
+            {osgb_cs}
+            {latlon_cs}
+        }}
+        """
+        nc_path = str(tmp_path / "tmp.nc")
+
+        # loader will warn that it can't find the longitude coord
+        with pytest.warns(
+            iris.warnings.IrisCfWarning,
+            match="Missing CF-netCDF grid mapping variable 'non_existent_crs'",
+        ):
+            cube = self.load_cube_from_cdl(cdl, None, nc_path, mocker)
+
+        assert len(cube.coords()) == 4
+        assert len(cube.coord_systems()) == 1
+
+        for coord_name in ["projection_x_coordinate", "projection_x_coordinate"]:
+            assert isinstance(
+                cube.coord(coord_name).coord_system,
+                iris.coord_systems.TransverseMercator,
+            )
+        for coord_name in ["latitude", "longitude"]:
+            assert cube.coord(coord_name).coord_system is None
+
+        assert cube.extended_grid_mapping is True
+
+    def test_two_coord_systems_missing_dim_crs(
+        self, osgb_cs, latlon_cs, mocker, tmp_path
+    ):
+        """Test invalid coordinate system mapping for Dim coords."""
+        cdl = f"""
+        netcdf tmp {{
+        dimensions:
+            x = 4 ;
+            y = 3 ;
+        variables:
+            float phenom(y, x) ;
+                phenom:standard_name = "air_pressure" ;
+                phenom:units = "Pa" ;
+                phenom:coordinates = "lat lon" ;
+                phenom:grid_mapping = "non_existent_crs: x y crsWGS84: lat lon" ;
+            {osgb_cs}
+            {latlon_cs}
+        }}
+        """
+        nc_path = str(tmp_path / "tmp.nc")
+
+        # loader will warn that it can't find the longitude coord
+        with pytest.warns(
+            iris.warnings.IrisCfWarning,
+            match="Missing CF-netCDF grid mapping variable 'non_existent_crs'",
+        ):
+            cube = self.load_cube_from_cdl(cdl, None, nc_path, mocker)
+
+        assert len(cube.coords()) == 2  # DimCoords won't be found
+        assert len(cube.coord_systems()) == 1
+
+        for coord_name in ["latitude", "longitude"]:
+            assert isinstance(
+                cube.coord(coord_name).coord_system, iris.coord_systems.GeogCS
+            )
+
+        assert cube.extended_grid_mapping is True
+
+    def test_two_coord_systems_invalid_grid_mapping(
+        self, osgb_cs, latlon_cs, mocker, tmp_path
+    ):
+        """Test invalid grid mapping doesn't load any coord systems and warns."""
+        cdl = f"""
+        netcdf tmp {{
+        dimensions:
+            x = 4 ;
+            y = 3 ;
+        variables:
+            float phenom(y, x) ;
+                phenom:standard_name = "air_pressure" ;
+                phenom:units = "Pa" ;
+                phenom:coordinates = "lat lon" ;
+                phenom:grid_mapping = "crsOSGB: crsWGS84:" ;
+            {osgb_cs}
+            {latlon_cs}
+        }}
+        """
+        nc_path = str(tmp_path / "tmp.nc")
+
+        # loader will warn that grid_mapping is invalid
+        with pytest.warns(
+            iris.warnings.IrisCfWarning, match="Error parsing `grid_mapping` attribute"
+        ):
+            cube = self.load_cube_from_cdl(cdl, None, nc_path, mocker)
+
+        assert len(cube.coords()) == 2
+        assert len(cube.coord_systems()) == 0
+        for coord in cube.coords():
+            assert coord.coord_system is None
+
+        assert cube.extended_grid_mapping is False
+
+    def test_one_coord_system_simple(self, osgb_cs, latlon_cs, mocker, tmp_path):
+        """Make sure the simple coord system syntax still works."""
+        cdl = f"""
+        netcdf tmp {{
+        dimensions:
+            x = 4 ;
+            y = 3 ;
+        variables:
+            float phenom(y, x) ;
+                phenom:standard_name = "air_pressure" ;
+                phenom:units = "Pa" ;
+                phenom:coordinates = "lat lon" ;
+                phenom:grid_mapping = "crsOSGB" ;
+            {osgb_cs}
+            {latlon_cs}
+        }}
+        """
+        nc_path = str(tmp_path / "tmp.nc")
+        cube = self.load_cube_from_cdl(cdl, None, nc_path, mocker)
+
+        assert len(cube.coord_systems()) == 1
+
+        for coord_name in ["projection_x_coordinate", "projection_x_coordinate"]:
+            assert isinstance(
+                cube.coord(coord_name).coord_system,
+                iris.coord_systems.TransverseMercator,
+            )
+        for coord_name in ["longitude", "latitude"]:
+            assert cube.coord(coord_name).coord_system is None
+
+        # Loading multiple coord systems or using extended grid mapping implies ordered axes:
+        assert cube.extended_grid_mapping is False
+
+
 if __name__ == "__main__":
     tests.main()
diff --git a/lib/iris/tests/unit/fileformats/netcdf/saver/test_Saver.py b/lib/iris/tests/unit/fileformats/netcdf/saver/test_Saver.py
index 292a21b18..0905c3d2a 100644
--- a/lib/iris/tests/unit/fileformats/netcdf/saver/test_Saver.py
+++ b/lib/iris/tests/unit/fileformats/netcdf/saver/test_Saver.py
@@ -16,6 +16,7 @@ from unittest import mock
 
 import numpy as np
 from numpy import ma
+import pytest
 
 import iris
 from iris.coord_systems import (
@@ -35,6 +36,7 @@ from iris.coord_systems import (
 from iris.coords import AncillaryVariable, AuxCoord, DimCoord
 from iris.cube import Cube
 from iris.fileformats.netcdf import Saver, _thread_safe_nc
+from iris.tests._shared_utils import assert_CDL
 import iris.tests.stock as stock
 
 
@@ -906,12 +908,122 @@ class Test__cf_coord_identity(tests.IrisTest):
         )
 
 
-class Test__create_cf_grid_mapping(tests.IrisTest):
+@pytest.fixture
+def transverse_mercator_cube_multi_cs():
+    """A transverse mercator cube with an auxiliary GeogGS coordinate system."""
+    data = np.arange(12).reshape(3, 4)
+    cube = Cube(data, "air_pressure_anomaly")
+    cube.extended_grid_mapping = True
+
+    geog_cs = GeogCS(6377563.396, 6356256.909)
+    trans_merc = TransverseMercator(
+        49.0, -2.0, -400000.0, 100000.0, 0.9996012717, geog_cs
+    )
+    coord = DimCoord(
+        np.arange(3),
+        "projection_y_coordinate",
+        units="m",
+        coord_system=trans_merc,
+    )
+    cube.add_dim_coord(coord, 0)
+    coord = DimCoord(
+        np.arange(4),
+        "projection_x_coordinate",
+        units="m",
+        coord_system=trans_merc,
+    )
+    cube.add_dim_coord(coord, 1)
+
+    # Add auxiliary lat/lon coords with a GeogCS coord system
+    coord = AuxCoord(
+        np.arange(3 * 4).reshape((3, 4)),
+        "longitude",
+        units="degrees",
+        coord_system=geog_cs,
+    )
+    cube.add_aux_coord(coord, (0, 1))
+
+    coord = AuxCoord(
+        np.arange(3 * 4).reshape((3, 4)),
+        "latitude",
+        units="degrees",
+        coord_system=geog_cs,
+    )
+    cube.add_aux_coord(coord, (0, 1))
+
+    return cube
+
+
+class Test_write_extended_grid_mapping:
+    def test_multi_cs(self, transverse_mercator_cube_multi_cs, tmp_path, request):
+        """Test writing a cube with multiple coordinate systems.
+        Should generate a grid mapping using extended syntax that references
+        both coordinate systems and the coords.
+        """
+        cube = transverse_mercator_cube_multi_cs
+        nc_path = tmp_path / "tmp.nc"
+        with Saver(nc_path, "NETCDF4") as saver:
+            saver.write(cube)
+        assert_CDL(request, nc_path)
+
+    def test_no_aux_cs(self, transverse_mercator_cube_multi_cs, tmp_path, request):
+        """Test when DimCoords have coord system, but AuxCoords do not.
+        Should write extended grid mapping for just DimCoords.
+        """
+        cube = transverse_mercator_cube_multi_cs
+        cube.coord("latitude").coord_system = None
+        cube.coord("longitude").coord_system = None
+
+        nc_path = tmp_path / "tmp.nc"
+        with Saver(nc_path, "NETCDF4") as saver:
+            saver.write(cube)
+        assert_CDL(request, nc_path)
+
+    def test_multi_cs_missing_coord(
+        self, transverse_mercator_cube_multi_cs, tmp_path, request
+    ):
+        """Test when we have a missing coordinate.
+        Grid mapping will fall back to simple mapping to DimCoord CS (no coords referenced).
+        """
+        cube = transverse_mercator_cube_multi_cs
+        cube.remove_coord("latitude")
+        nc_path = tmp_path / "tmp.nc"
+        with Saver(nc_path, "NETCDF4") as saver:
+            saver.write(cube)
+        assert_CDL(request, nc_path)
+
+    def test_no_cs(self, transverse_mercator_cube_multi_cs, tmp_path, request):
+        """Test when no coordinate systems associated with cube coords.
+        Grid mapping will not be generated at all.
+        """
+        cube = transverse_mercator_cube_multi_cs
+        for coord in cube.coords():
+            coord.coord_system = None
+
+        nc_path = tmp_path / "tmp.nc"
+        with Saver(nc_path, "NETCDF4") as saver:
+            saver.write(cube)
+        assert_CDL(request, nc_path)
+
+
+class Test_create_cf_grid_mapping:
+    """Tests correct generation of CF grid_mapping variable attributes.
+
+    Note: The first 3 tests are run with the "extended grid" mapping
+    both enabled (the default for all these tests) and disabled. This
+    controls the output of the WKT attribute.
+    """
+
+    @pytest.fixture(autouse=True)
+    def _setup(self):
+        self._extended_grid_mapping = True  # forces WKT strings to be written
+
     def _cube_with_cs(self, coord_system):
         """Return a simple 2D cube that uses the given coordinate system."""
         cube = stock.lat_lon_cube()
         x, y = cube.coord("longitude"), cube.coord("latitude")
         x.coord_system = y.coord_system = coord_system
+        cube.extended_grid_mapping = self._extended_grid_mapping
         return cube
 
     def _grid_mapping_variable(self, coord_system):
@@ -930,14 +1042,15 @@ class Test__create_cf_grid_mapping(tests.IrisTest):
         grid_variable = NCMock(name="NetCDFVariable")
         create_var_fn = mock.Mock(side_effect=[grid_variable])
         dataset = mock.Mock(variables=[], createVariable=create_var_fn)
-        saver = mock.Mock(spec=Saver, _coord_systems=[], _dataset=dataset)
         variable = NCMock()
 
-        # This is the method we're actually testing!
-        Saver._create_cf_grid_mapping(saver, cube, variable)
+        saver = Saver(dataset, "NETCDF4", compute=False)
+
+        # The method we want to test:
+        saver._create_cf_grid_mapping(cube, variable)
 
-        self.assertEqual(create_var_fn.call_count, 1)
-        self.assertEqual(variable.grid_mapping, grid_variable.grid_mapping_name)
+        assert create_var_fn.call_count == 1
+        assert variable.grid_mapping, grid_variable.grid_mapping_name
         return grid_variable
 
     def _variable_attributes(self, coord_system):
@@ -957,12 +1070,14 @@ class Test__create_cf_grid_mapping(tests.IrisTest):
         actual = self._variable_attributes(coord_system)
 
         # To see obvious differences, check that they keys are the same.
-        self.assertEqual(sorted(actual.keys()), sorted(expected.keys()))
+        assert sorted(actual.keys()) == sorted(expected.keys())
         # Now check that the values are equivalent.
-        self.assertEqual(actual, expected)
+        assert actual == expected
 
-    def test_rotated_geog_cs(self):
+    def test_rotated_geog_cs(self, extended_grid_mapping):
+        self._extended_grid_mapping = extended_grid_mapping
         coord_system = RotatedGeogCS(37.5, 177.5, ellipsoid=GeogCS(6371229.0))
+
         expected = {
             "grid_mapping_name": b"rotated_latitude_longitude",
             "north_pole_grid_longitude": 0.0,
@@ -971,18 +1086,47 @@ class Test__create_cf_grid_mapping(tests.IrisTest):
             "longitude_of_prime_meridian": 0.0,
             "earth_radius": 6371229.0,
         }
+        if extended_grid_mapping:
+            expected["crs_wkt"] = (
+                'GEOGCRS["unnamed",BASEGEOGCRS["unknown",DATUM["unknown",'
+                'ELLIPSOID["unknown",6371229,0,LENGTHUNIT["metre",1,ID['
+                '"EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",'
+                '0.0174532925199433],ID["EPSG",8901]]],DERIVINGCONVERSION['
+                '"unknown",METHOD["PROJ ob_tran o_proj=latlon"],PARAMETER['
+                '"o_lon_p",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG"'
+                ',9122]]],PARAMETER["o_lat_p",37.5,ANGLEUNIT["degree",'
+                '0.0174532925199433,ID["EPSG",9122]]],PARAMETER["lon_0",357.5'
+                ',ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CS['
+                'ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT['
+                '"degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude"'
+                ',north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID['
+                '"EPSG",9122]]]]'
+            )
+
         self._test(coord_system, expected)
 
-    def test_spherical_geog_cs(self):
+    def test_spherical_geog_cs(self, extended_grid_mapping):
+        self._extended_grid_mapping = extended_grid_mapping
         coord_system = GeogCS(6371229.0)
         expected = {
             "grid_mapping_name": b"latitude_longitude",
             "longitude_of_prime_meridian": 0.0,
             "earth_radius": 6371229.0,
         }
+        if extended_grid_mapping:
+            expected["crs_wkt"] = (
+                'GEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6371229'
+                ',0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich"'
+                ',0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]],CS'
+                '[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT['
+                '"degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude"'
+                ',north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID['
+                '"EPSG",9122]]]]'
+            )
         self._test(coord_system, expected)
 
-    def test_elliptic_geog_cs(self):
+    def test_elliptic_geog_cs(self, extended_grid_mapping):
+        self._extended_grid_mapping = extended_grid_mapping
         coord_system = GeogCS(637, 600)
         expected = {
             "grid_mapping_name": b"latitude_longitude",
@@ -990,6 +1134,17 @@ class Test__create_cf_grid_mapping(tests.IrisTest):
             "semi_minor_axis": 600.0,
             "semi_major_axis": 637.0,
         }
+        if extended_grid_mapping:
+            expected["crs_wkt"] = (
+                'GEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",637,'
+                '17.2162162162162,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],'
+                'PRIMEM["Reference meridian",0,ANGLEUNIT["degree",'
+                '0.0174532925199433,ID["EPSG",9122]]],CS[ellipsoidal,2],AXIS'
+                '["longitude",east,ORDER[1],ANGLEUNIT["degree",'
+                '0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,'
+                'ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",'
+                "9122]]]]"
+            )
         self._test(coord_system, expected)
 
     def test_lambert_conformal(self):
@@ -1002,6 +1157,25 @@ class Test__create_cf_grid_mapping(tests.IrisTest):
             ellipsoid=GeogCS(6371000),
         )
         expected = {
+            "crs_wkt": (
+                'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",'
+                'ELLIPSOID["unknown",6371000,0,LENGTHUNIT["metre",1,ID["EPSG"'
+                ',9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",'
+                '0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",'
+                'METHOD["Lambert Conic Conformal (2SP)",ID["EPSG",9802]],'
+                'PARAMETER["Latitude of false origin",44,ANGLEUNIT["degree",'
+                '0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of '
+                'false origin",2,ANGLEUNIT["degree",0.0174532925199433],ID['
+                '"EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",'
+                '38,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],'
+                'PARAMETER["Latitude of 2nd standard parallel",50,ANGLEUNIT'
+                '["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER['
+                '"Easting at false origin",-2,LENGTHUNIT["metre",1],ID["EPSG",'
+                '8826]],PARAMETER["Northing at false origin",-5,LENGTHUNIT['
+                '"metre",1],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["(E)",east,'
+                'ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",'
+                'north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'
+            ),
             "grid_mapping_name": b"lambert_conformal_conic",
             "latitude_of_projection_origin": 44,
             "longitude_of_central_meridian": 2,
@@ -1022,6 +1196,21 @@ class Test__create_cf_grid_mapping(tests.IrisTest):
             ellipsoid=GeogCS(6377563.396, 6356256.909),
         )
         expected = {
+            "crs_wkt": (
+                'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIP'
+                'SOID["unknown",6377563.396,299.324961266495,LENGTHUNIT["metre'
+                '",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree'
+                '",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",'
+                'METHOD["Lambert Azimuthal Equal Area",ID["EPSG",9820]],PARAME'
+                'TER["Latitude of natural origin",52,ANGLEUNIT["degree",0.0174'
+                '532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natura'
+                'l origin",10,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG"'
+                ',8802]],PARAMETER["False easting",100,LENGTHUNIT["metre",1],I'
+                'D["EPSG",8806]],PARAMETER["False northing",200,LENGTHUNIT["me'
+                'tre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORD'
+                'ER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north'
+                ',ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'
+            ),
             "grid_mapping_name": b"lambert_azimuthal_equal_area",
             "latitude_of_projection_origin": 52,
             "longitude_of_projection_origin": 10,
@@ -1043,6 +1232,25 @@ class Test__create_cf_grid_mapping(tests.IrisTest):
             ellipsoid=GeogCS(6377563.396, 6356256.909),
         )
         expected = {
+            "crs_wkt": (
+                'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIP'
+                'SOID["unknown",6377563.396,299.324961266495,LENGTHUNIT["metre'
+                '",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree'
+                '",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",'
+                'METHOD["Albers Equal Area",ID["EPSG",9822]],PARAMETER["Latitu'
+                'de of false origin",52,ANGLEUNIT["degree",0.0174532925199433]'
+                ',ID["EPSG",8821]],PARAMETER["Longitude of false origin",10,AN'
+                'GLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMET'
+                'ER["Latitude of 1st standard parallel",38,ANGLEUNIT["degree",'
+                '0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2'
+                'nd standard parallel",50,ANGLEUNIT["degree",0.017453292519943'
+                '3],ID["EPSG",8824]],PARAMETER["Easting at false origin",100,L'
+                'ENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at '
+                'false origin",200,LENGTHUNIT["metre",1],ID["EPSG",8827]]],CS['
+                'Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID'
+                '["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",'
+                '1,ID["EPSG",9001]]]]'
+            ),
             "grid_mapping_name": b"albers_conical_equal_area",
             "latitude_of_projection_origin": 52,
             "longitude_of_central_meridian": 10,
@@ -1075,6 +1283,24 @@ class Test__create_cf_grid_mapping(tests.IrisTest):
             ellipsoid=ellipsoid,
         )
         expected = {
+            "crs_wkt": (
+                'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIP'
+                'SOID["unknown",6377563.396,299.324961266495,LENGTHUNIT["metre'
+                '",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree'
+                '",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",'
+                'METHOD["Vertical Perspective",ID["EPSG",9838]],PARAMETER["Lat'
+                'itude of topocentric origin",1,ANGLEUNIT["degree",0.017453292'
+                '5199433],ID["EPSG",8834]],PARAMETER["Longitude of topocentric'
+                ' origin",2,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8'
+                '835]],PARAMETER["Ellipsoidal height of topocentric origin",0,'
+                'LENGTHUNIT["metre",1],ID["EPSG",8836]],PARAMETER["Viewpoint h'
+                'eight",2000000,LENGTHUNIT["metre",1],ID["EPSG",8840]],PARAMET'
+                'ER["False easting",100,LENGTHUNIT["metre",1],ID["EPSG",8806]]'
+                ',PARAMETER["False northing",200,LENGTHUNIT["metre",1],ID["EPS'
+                'G",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNI'
+                'T["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGT'
+                'HUNIT["metre",1,ID["EPSG",9001]]]]'
+            ),
             "grid_mapping_name": b"vertical_perspective",
             "latitude_of_projection_origin": latitude_of_projection_origin,
             "longitude_of_projection_origin": longitude_of_projection_origin,
@@ -1109,6 +1335,23 @@ class Test__create_cf_grid_mapping(tests.IrisTest):
             ellipsoid=ellipsoid,
         )
         expected = {
+            "crs_wkt": (
+                'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIP'
+                'SOID["unknown",6377563.396,299.324961266495,LENGTHUNIT["metre'
+                '",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree'
+                '",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",'
+                'METHOD["Geostationary Satellite (Sweep X)"],PARAMETER["Longit'
+                'ude of natural origin",2,ANGLEUNIT["degree",0.017453292519943'
+                '3],ID["EPSG",8802]],PARAMETER["Satellite Height",2000000,LENG'
+                'THUNIT["metre",1,ID["EPSG",9001]]],PARAMETER["False easting",'
+                '100,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False n'
+                'orthing",200,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Carte'
+                'sian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPS'
+                'G",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID['
+                '"EPSG",9001]]],REMARK["PROJ CRS string: +proj=geos +a=6377563'
+                ".396 +b=6356256.909 +lon_0=2.0 +lat_0=0.0 +h=2000000.0 +x_0=1"
+                '00.0 +y_0=200.0 +units=m +sweep=x +no_defs"]]'
+            ),
             "grid_mapping_name": b"geostationary",
             "latitude_of_projection_origin": latitude_of_projection_origin,
             "longitude_of_projection_origin": longitude_of_projection_origin,
@@ -1126,8 +1369,30 @@ class Test__create_cf_grid_mapping(tests.IrisTest):
         # Some none-default settings to confirm all parameters are being
         #  handled.
 
+        wkt_template = (
+            'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIP'
+            'SOID["unknown",1,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PR'
+            'IMEM["Reference meridian",0,ANGLEUNIT["degree",0.017453292519'
+            '9433,ID["EPSG",9122]]]],CONVERSION["unknown",METHOD["Hotine O'
+            'blique Mercator (variant B)",ID["EPSG",9815]],PARAMETER["Lati'
+            'tude of projection centre",89.9,ANGLEUNIT["degree",0.01745329'
+            '25199433],ID["EPSG",8811]],PARAMETER["Longitude of projection'
+            ' centre",45,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",'
+            '8812]],PARAMETER["Azimuth at projection centre",{angle},ANGLEUNIT['
+            '"degree",0.0174532925199433],ID["EPSG",8813]],PARAMETER["Angl'
+            'e from Rectified to Skew Grid",{angle},ANGLEUNIT["degree",0.017453'
+            '2925199433],ID["EPSG",8814]],PARAMETER["Scale factor at proje'
+            'ction centre",0.939692620786,SCALEUNIT["unity",1],ID["EPSG",8'
+            '815]],PARAMETER["Easting at projection centre",1000000,LENGTH'
+            'UNIT["metre",1],ID["EPSG",8816]],PARAMETER["Northing at proje'
+            'ction centre",-2000000,LENGTHUNIT["metre",1],ID["EPSG",8817]]'
+            '],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre"'
+            ',1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["me'
+            'tre",1,ID["EPSG",9001]]]]'
+        )
+
         kwargs_rotated = dict(
-            latitude_of_projection_origin=90.0,
+            latitude_of_projection_origin=89.9,
             longitude_of_projection_origin=45.0,
             false_easting=1000000.0,
             false_northing=-2000000.0,
@@ -1142,8 +1407,9 @@ class Test__create_cf_grid_mapping(tests.IrisTest):
         expected_rotated = dict(
             # Automatically converted to oblique_mercator in line with CF 1.11 .
             grid_mapping_name=b"oblique_mercator",
-            # Azimuth should be automatically populated.
+            # Azimuth and crs_wkt should be automatically populated.
             azimuth_of_central_line=90.0,
+            crs_wkt=wkt_template.format(angle="89.999"),
             **kwargs_rotated,
         )
         # Convert the ellipsoid
@@ -1156,6 +1422,7 @@ class Test__create_cf_grid_mapping(tests.IrisTest):
 
         # Same as rotated, but different azimuth.
         expected_oblique = dict(expected_rotated, **oblique_azimuth)
+        expected_oblique["crs_wkt"] = wkt_template.format(angle="45")
 
         oblique = ObliqueMercator(**kwargs_oblique)
         rotated = RotatedMercator(**kwargs_rotated)
@@ -1167,5 +1434,16 @@ class Test__create_cf_grid_mapping(tests.IrisTest):
             self._test(coord_system, expected)
 
 
+@pytest.fixture(
+    params=[
+        pytest.param(True, id="extended_grid_mapping"),
+        pytest.param(False, id="no_extended_grid_mapping"),
+    ]
+)
+def extended_grid_mapping(request):
+    """Fixture for enabling/disabling extended grid mapping."""
+    return request.param
+
+
 if __name__ == "__main__":
     tests.main()