From 82be35fc91e84b8d7bde069e0f436e5759ff8014 Mon Sep 17 00:00:00 2001
From: Henry Wright <henrywright@sky.com>
Date: Fri, 28 Feb 2025 17:49:54 +0000
Subject: [PATCH 01/11] reuploading changes for bounds of derived variables

---
 lib/iris/__init__.py       |   5 ++
 lib/iris/fileformats/cf.py | 132 ++++++++++++++++++++++++++++++-------
 2 files changed, 113 insertions(+), 24 deletions(-)

diff --git a/lib/iris/__init__.py b/lib/iris/__init__.py
index d622bd18b0..0dbdc40510 100644
--- a/lib/iris/__init__.py
+++ b/lib/iris/__init__.py
@@ -169,6 +169,7 @@ def __init__(
         pandas_ndim=False,
         save_split_attrs=False,
         date_microseconds=False,
+        derived_bounds=False,
     ):
         """Container for run-time options controls.
 
@@ -202,6 +203,8 @@ def __init__(
             behaviour, such as when using :class:`~iris.Constraint`, and you
             may need to defend against floating point precision issues where
             you didn't need to before.
+        derived_bounds : bool, default=False
+            When True, uses the new form for deriving bounds with the load.
 
         """
         # The flag 'example_future_flag' is provided as a reference for the
@@ -215,6 +218,7 @@ def __init__(
         self.__dict__["pandas_ndim"] = pandas_ndim
         self.__dict__["save_split_attrs"] = save_split_attrs
         self.__dict__["date_microseconds"] = date_microseconds
+        self.__dict__["derived_bounds"] = derived_bounds
 
         # TODO: next major release: set IrisDeprecation to subclass
         #  DeprecationWarning instead of UserWarning.
@@ -228,6 +232,7 @@ def __repr__(self):
             self.pandas_ndim,
             self.save_split_attrs,
             self.date_microseconds,
+            self.derived_bounds,
         )
 
     # deprecated_options = {'example_future_flag': 'warning',}
diff --git a/lib/iris/fileformats/cf.py b/lib/iris/fileformats/cf.py
index 82be010c6e..69072f274a 100644
--- a/lib/iris/fileformats/cf.py
+++ b/lib/iris/fileformats/cf.py
@@ -16,6 +16,7 @@
 
 from abc import ABCMeta, abstractmethod
 from collections.abc import Iterable, MutableMapping
+import contextlib
 import os
 import re
 from typing import ClassVar, Optional
@@ -94,6 +95,8 @@ def __init__(self, name, data):
         #: CF-netCDF formula terms that his variable participates in.
         self.cf_terms_by_root = {}
 
+        self._to_be_promoted = False
+
         self.cf_attrs_reset()
 
     @staticmethod
@@ -1416,16 +1419,66 @@ def _translate(self, variables):
         # Identify and register all CF formula terms.
         formula_terms = _CFFormulaTermsVariable.identify(variables)
 
-        for cf_var in formula_terms.values():
-            for cf_root, cf_term in cf_var.cf_terms_by_root.items():
-                # Ignore formula terms owned by a bounds variable.
-                if cf_root not in self.cf_group.bounds:
-                    cf_name = cf_var.cf_name
-                    if cf_var.cf_name not in self.cf_group:
-                        self.cf_group[cf_name] = CFAuxiliaryCoordinateVariable(
-                            cf_name, cf_var.cf_data
-                        )
-                    self.cf_group[cf_name].add_formula_term(cf_root, cf_term)
+        if iris.FUTURE.derived_bounds:
+            # cf_var = CFFormulaTermsVariable (loops through everything that appears in formula terms)
+            for cf_var in formula_terms.values():
+                # eg. eta:'a' | cf_root = eta and cf_term = a. cf_var.cf_terms_by_root = {'eta': 'a'} (looking at all appearances in formula terms)
+                for cf_root, cf_term in cf_var.cf_terms_by_root.items():
+                    # gets set to the bounds of the coord from cf_root_coord
+                    bounds_name = None
+                    # cf_root_coord = CFCoordinateVariable of the coordinate relating to the root
+                    cf_root_coord = self.cf_group.coordinates.get(cf_root)
+                    if cf_root_coord is None:
+                        cf_root_coord = self.cf_group.auxiliary_coordinates.get(cf_root)
+                    with contextlib.suppress(AttributeError):
+                        # Copes with cf_root_coord not existing, OR not having
+                        #  `bounds` attribute.
+                        bounds_name = cf_root_coord.bounds
+
+                    if bounds_name is not None:
+                        try:
+                            # This will error if more or less than 1 variable is found.
+                            # TODO: try a try/except here or logical alternative
+                            (bounds_var,) = [
+                                # loop through all formula terms and add them if they have a cf_term_by_root
+                                # where (bounds of cf_root): cf_term (same as before)
+                                f
+                                for f in formula_terms.values()
+                                if f.cf_terms_by_root.get(bounds_name) == cf_term
+                            ]
+                            if bounds_var != cf_var:
+                                cf_var.bounds = bounds_var.cf_name
+                                new_var = CFBoundaryVariable(
+                                    bounds_var.cf_name, bounds_var.cf_data
+                                )
+                                new_var.add_formula_term(bounds_name, cf_term)
+                                self.cf_group[bounds_var.cf_name] = new_var
+                        except ValueError:
+                            # Modify the boundary_variable set _to_be_promoted to True
+                            self.cf_group.get(bounds_name)._to_be_promoted = True
+
+                    if cf_root not in self.cf_group.bounds:
+                        cf_name = cf_var.cf_name
+                        if cf_var.cf_name not in self.cf_group:
+                            new_var = CFAuxiliaryCoordinateVariable(
+                                cf_name, cf_var.cf_data
+                            )
+                            if hasattr(cf_var, "bounds"):
+                                new_var.bounds = cf_var.bounds
+                            new_var.add_formula_term(cf_root, cf_term)
+                            self.cf_group[cf_name] = new_var
+
+        else:
+            for cf_var in formula_terms.values():
+                for cf_root, cf_term in cf_var.cf_terms_by_root.items():
+                    # Ignore formula terms owned by a bounds variable.
+                    if cf_root not in self.cf_group.bounds:
+                        cf_name = cf_var.cf_name
+                        if cf_var.cf_name not in self.cf_group:
+                            self.cf_group[cf_name] = CFAuxiliaryCoordinateVariable(
+                                cf_name, cf_var.cf_data
+                            )
+                        self.cf_group[cf_name].add_formula_term(cf_root, cf_term)
 
         # Determine the CF data variables.
         data_variable_names = (
@@ -1454,7 +1507,7 @@ def _span_check(
                 """Sanity check dimensionality."""
                 var = self.cf_group[var_name]
                 # No span check is necessary if variable is attached to a mesh.
-                if is_mesh_var or var.spans(cf_variable):
+                if (is_mesh_var or var.spans(cf_variable)) and not var._to_be_promoted:
                     cf_group[var_name] = var
                 else:
                     # Register the ignored variable.
@@ -1498,6 +1551,14 @@ def _span_check(
                 for cf_name in match:
                     _span_check(cf_name)
 
+            if iris.FUTURE.derived_bounds:
+                if hasattr(cf_variable, "bounds"):
+                    if cf_variable.bounds not in cf_group:
+                        bounds_var = self.cf_group[cf_variable.bounds]
+                        # TODO: warning if span fails
+                        if bounds_var.spans(cf_variable):
+                            cf_group[cf_variable.bounds] = bounds_var
+
             # Build CF data variable relationships.
             if isinstance(cf_variable, CFDataVariable):
                 # Add global netCDF attributes.
@@ -1539,19 +1600,42 @@ def _span_check(
         # Determine whether there are any formula terms that
         # may be promoted to a CFDataVariable and restrict promotion to only
         # those formula terms that are reference surface/phenomenon.
-        for cf_var in self.cf_group.formula_terms.values():
-            for cf_root, cf_term in cf_var.cf_terms_by_root.items():
-                cf_root_var = self.cf_group[cf_root]
-                name = cf_root_var.standard_name or cf_root_var.long_name
-                terms = reference_terms.get(name, [])
-                if isinstance(terms, str) or not isinstance(terms, Iterable):
-                    terms = [terms]
-                cf_var_name = cf_var.cf_name
-                if cf_term in terms and cf_var_name not in self.cf_group.promoted:
-                    data_var = CFDataVariable(cf_var_name, cf_var.cf_data)
-                    self.cf_group.promoted[cf_var_name] = data_var
-                    _build(data_var)
-                    break
+        if iris.FUTURE.derived_bounds:
+            for cf_var in self.cf_group.formula_terms.values():
+                if self.cf_group[cf_var.cf_name] is CFBoundaryVariable:
+                    continue
+                else:
+                    for cf_root, cf_term in cf_var.cf_terms_by_root.items():
+                        cf_root_var = self.cf_group[cf_root]
+                        if not hasattr(cf_root_var, "standard_name"):
+                            continue
+                        name = cf_root_var.standard_name or cf_root_var.long_name
+                        terms = reference_terms.get(name, [])
+                        if isinstance(terms, str) or not isinstance(terms, Iterable):
+                            terms = [terms]
+                        cf_var_name = cf_var.cf_name
+                        if (
+                            cf_term in terms
+                            and cf_var_name not in self.cf_group.promoted
+                        ):
+                            data_var = CFDataVariable(cf_var_name, cf_var.cf_data)
+                            self.cf_group.promoted[cf_var_name] = data_var
+                            _build(data_var)
+                            break
+        else:
+            for cf_var in self.cf_group.formula_terms.values():
+                for cf_root, cf_term in cf_var.cf_terms_by_root.items():
+                    cf_root_var = self.cf_group[cf_root]
+                    name = cf_root_var.standard_name or cf_root_var.long_name
+                    terms = reference_terms.get(name, [])
+                    if isinstance(terms, str) or not isinstance(terms, Iterable):
+                        terms = [terms]
+                    cf_var_name = cf_var.cf_name
+                    if cf_term in terms and cf_var_name not in self.cf_group.promoted:
+                        data_var = CFDataVariable(cf_var_name, cf_var.cf_data)
+                        self.cf_group.promoted[cf_var_name] = data_var
+                        _build(data_var)
+                        break
         # Promote any ignored variables.
         promoted = set()
         not_promoted = ignored.difference(promoted)

From 4d424393f2b1c17717b6dca0ead898279b48db63 Mon Sep 17 00:00:00 2001
From: Patrick Peglar <patrick.peglar@metoffice.gov.uk>
Date: Wed, 14 May 2025 13:59:23 +0100
Subject: [PATCH 02/11] New testing facilities (WIP temporary).

---
 .../netcdf/derived_bounds/__init__.py         |   6 +
 .../temp_nc_sources/a_new_file.nc             | Bin 0 -> 17045 bytes
 .../temp_nc_sources/a_new_file.ncdump.txt     |  47 ++++++
 .../temp_nc_sources/hybrid_height.ncdump.txt  |  72 +++++++++
 .../derived_bounds/test_bounds_files.py       | 138 ++++++++++++++++++
 5 files changed, 263 insertions(+)
 create mode 100644 lib/iris/tests/integration/netcdf/derived_bounds/__init__.py
 create mode 100644 lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/a_new_file.nc
 create mode 100644 lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/a_new_file.ncdump.txt
 create mode 100644 lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/hybrid_height.ncdump.txt
 create mode 100644 lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py

diff --git a/lib/iris/tests/integration/netcdf/derived_bounds/__init__.py b/lib/iris/tests/integration/netcdf/derived_bounds/__init__.py
new file mode 100644
index 0000000000..979cb7bf6a
--- /dev/null
+++ b/lib/iris/tests/integration/netcdf/derived_bounds/__init__.py
@@ -0,0 +1,6 @@
+# Copyright Iris contributors
+#
+# This file is part of Iris and is released under the BSD license.
+# See LICENSE in the root of the repository for full licensing details.
+"""Integration tests for loading and saving netcdf files with derived coordinate bounds."""
+
diff --git a/lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/a_new_file.nc b/lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/a_new_file.nc
new file mode 100644
index 0000000000000000000000000000000000000000..148d8aa8238e2dd144a6322af0de3e3094cf7e39
GIT binary patch
literal 17045
zcmeHPeQZ?65ubN=b`Ek71I8whG!L4lR-vv<2qq3R_St9b;uy<Efj=4@&%QIh$iBPY
zoheBZ2vOBEZPQjFB~?`NDV0Q3R98)+T9MLLO{%mgQa@BF>W6Ce18Ab85Va~OYSWqB
z{oR-E3^)c+@)o#xZ+3Qfc6N4l-^?x#ba!?wn7@2JN%H&2M>XmZKku3J&7j8M9r?TP
zPyuf3?i}z0raY@uzBgdVdOoTjes7@Fpx@_v$dad}*DxnUUW?FRDGYo{=qy@De*P2b
z>*!A>Q?Yb59?LL(3ICHORb>qK#L}5~GO=mxnnWx+6dm5QuCb-@#x*0+;hQ!!H#Rpl
z-b5Rij7qT;Xz4B?mgLdEx8h8#>gs~D$66*h%hGLw13+)BFty0OmARJGM(eJmb?<qJ
zs%j+1<MCKJM{*zC><b9C5aVAgIbLgnks0zJvxPKqwNMZRmuom;ApII#XuSJ-c}O6-
z^Tf4U`?g!>q<t%|Q`3vcrg}fs?`JVo`JsIxlNz=HUOLjQ<@2@%xmj6Nc~&+FsYTUl
zR(ABZhr|3Xvr)_nw6C-0woqTVXKP<1+|k|}qQEZ6v#3gy>TAC(#Kfvp8$*=ZHWJU+
z{MWI04tI2R+0pn|EFsg{9vsiwiDWjPwZ}=vox#yqqm9%|fps!n#bN}wok*WD1nWBE
zWAP}rQp3|oOb)S3B%xIef4Vn#%Y~5enN6xO`XEbKo6uNb-*|Zm8xsptOU8nplMegn
zMP=8C{;LI@)8Xp>Lmgf_OC9FLh0_a1dwp8M&s&ySJ1mdZZs<BBUHU$EYgU_Cv#`m)
zn&quL8vVk%MYTL@m*r+{M>3g?#uH9ft@dSRuUUQ3BS5K3q!H-R7Dw+MhLh>B@lhv|
zjitx1dswMj6?!74)ozarCZZX8u+<Lgr&OlZ?ho6kCI(Gnz{Ip%#b$09<;`1fPk0-V
ztj5yedj2PiHMRy|3Q)v8F^_*Z39j&6m_D*GwS05O)?PI{MuO^@^Q2x$TMPnqmg2Q4
z<mn%**5$P6%i~c`mCl}4zbQ5ie?<vB>IWZysNZyVT0MG0I-Q>@E7XsTCF=jV>F9zI
zdf*q9;8&ZDUR8nzpHNuT5B#DMdepD>Xq;C{FH4~$;6Y!j(K$T$F=(m*!$7#um4GV&
zR|2jCTnV@ma3$bMz?Fb20apU91TJd{bPaU(tB*vFaw(X;A~@1&zW%=khd(~+0S!IA
zy3ZymE|x2)t8QSgMA$~TAKkm_C&6vqoy-wC9OPR#+Aw9>o=0dU82s%GifW!v$rm|O
zU8C#BQ~6g--rp7Gx(L01Zvzmm`lz4wQUky0oVB*jSz!nOgOk6!Hn^>~i@AGLvRbw8
zr}t^QhUx^-BtJwiA3wB_d{S4K{oo8<-TR3F+98NriCx)5CqplX>5w3}RdI2wJ@(X{
z)Wr{Uz2+kFl>=X<I|QM;%tW!kClAs~Qp9hmKG*&8trzKap;)7(df&=L|9*4Hzu<9@
z95>Ci1l0J?D;n;iaVr|LA~h)*t0MI%QjH?DC>rp^k1YEBTzAK<?>r^2wWE7W`M_$R
ze>Of_nb497A$+#t)7$;jz5x@JQgnGwpLv}IrUjlV1djh<T&rmL*>}GB#>MQNECg?S
zch#w`a5zL*U@W%PKU2Qu$;C@PM6~O*C-!%BcLsAIXw}{4T8=$?zA?da(LiAV-H&A*
zem=igo-kw4d|nUyG@4B8jwGD1nBq`J8^vi`P`~YD?cwp!QF}DDCpMa~ctMkrsbnUe
zjqgE06{jrZESmgW0_PGkStk*7(owV|?%&k~tvK1SWF|EdOUEK3cMhiG(MTq~d(4TX
z(y>ftoUuc>KnJQU7p#!T0vMNIg^E+yMIj1{*Dc|GFgebF4BcN*SQKM<s=$T>F|r)Q
z7zlNR210!up@<G*)bS{b6Fll<olGptPctV3J90r<LoRL3lqyXJ1RMY{9Z*-87vjS7
z=Z_S=5p6JwJ+Ql*@ucDNywtmCE+n91j2z{f<vAERdq8y9j9<vGRQViCKkvyTL9@;w
zruZ45pUD#NCN4KK9D67k8y=1i#bXIJ2+6d~r-`!hA!n2`5QZ$V6;l7Q#?xn%vCOQM
zSFZ36yj_&}_VBWn`xNFZKVsgsG+w$+lVaWlRtkTxspE-wmZ!9`L{dNIbn}TnJ@QO$
zIFy*)f5=xe^vZh8&=)kI%EPnM438RBrpEAiE^^3s>CtYr`!)>EHBcXBvP75tFyoO|
z396vtyMj?`r_PgKgfa7&@=v_p<Vxz%wt3ViZwB?MAHV)>^~qHwbP_a;HoRigFpz7d
zW648rDynLsk`EE464mFFr@G;AqEX|f{hl;NbIQQe8n%7CO_8?k*D&?4hC~LUoL|4U
zJ!a^RFg(;jZgM5yO2CzXD*;ymt^_Vg3Fuwrhu`NWe~MiOlCZ6<yo(&Uz?=NRS#0vp
zaT8v5u;<1Lvk7l8XIsy{2~S-joA7&nFmY5n`LQF$w+f(E2=S1Mli%i7of9$NgNSSZ
ze&?~b|E`?8RZnbT&)EZ9pq4y`=QYr2eq2i5+Sl7a?fh6rUwnM>WAqn={q)Dz-^`hU
z{o6aOn`r?zQcr&!pBkbajIAf{j)S}D7yNiNy*Ba32bI&Njza%g{UjY$*oiGCPEdj$
zSJ3UjyWXW2__3C5ANbyH_`yvVRRU`K;TV!ZmvPl<!;j)`tIK~bX<fiMSTo%GJ*AT#
z(JH)Zb;?z(haP8K&6xW^tknfkFQaMmk*n1a=U0u!g+5>RGzeTcTM5kI2H&OqBy)E#
zTXaEG3r>GghT%J=PR-22+#U{>5NaV0ayEGTfrpMWxJYT{-il{x?j$ZVsIa#(g?@Ex
zVxRDUyKo4AV=e~R)l;RoxD<WaFg}3rWlJw5!FP>!_^R<<(|F*3in8Dur4+hiY<Q=Y
zvtwTN=i`jJrDEGc;~jOfIB{<I;6d_oEClDxEgzWj$%O|T-?x0I*^fX$j%gpPrPU56
zyyb(*tkxGNay|0Fcv?Yt<SidYty<r$$aU8TTk7c#Pu?OOU`_DBqI&w9-r5Kx?5Zyt
z8%uo<eEE7u*-+so4>Y9|R=&xj-MT$3TA*j&QQq62DKU0%>t3vzJkm2!;J_-e4&^-I
zl8gTLJExv|JiiN9Oe3$dOv(q2D}SVXse(*qxW&xe3zz!7WT79wn0x9Ly17cnuYdC0
za@G@bJn*XQzc1BDph!j%V9qjVG-385|DNrgr(}J|_3npq4(L+t9kdd#w0F<+G~KV~
z{xQqiAYHcjmb`-MH?Q-kp%R4Nm$?h9?wrd_bS<-!ce9D+YVd|t_ZvpvAz9b}U**uw
zHqGhmZu4Pj0z|K6qyucZ#17~iedXfr65oI<oOqgV`;GB`7~o_%%r0N@vPGLd4W$4I
yZ3(2HIp=fD58MJub=W=P!~x#rd_ZA-Ex)(i{k!#--)!ZMsf5`IA$Nba+W$XBgmt_C

literal 0
HcmV?d00001

diff --git a/lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/a_new_file.ncdump.txt b/lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/a_new_file.ncdump.txt
new file mode 100644
index 0000000000..ff9f9f2d4f
--- /dev/null
+++ b/lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/a_new_file.ncdump.txt
@@ -0,0 +1,47 @@
+netcdf a_new_file {
+dimensions:
+	eta = 1 ;
+	lat = 1 ;
+	lon = 1 ;
+	bnds = 2 ;
+variables:
+	double eta(eta) ;
+		eta:long_name = "eta at full levels" ;
+		eta:positive = "down" ;
+		eta:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
+		eta:formula_terms = "a: A b: B ps: PS p0: P0" ;
+		eta:bounds = "eta_bnds" ;
+	double eta_bnds(eta, bnds) ;
+		eta_bnds:formula_terms = "a: A_bnds b: B_bnds ps: PS p0: P0" ;
+	double A(eta) ;
+		A:long_name = "a coefficient for vertical coordinate at full levels" ;
+		A:units = "1" ;
+	double A_bnds(eta, bnds) ;
+	double B(eta) ;
+		B:long_name = "b coefficient for vertical coordinate at full levels" ;
+		B:units = "1" ;
+	double B_bnds(eta, bnds) ;
+	double PS(lat, lon) ;
+		PS:units = "Pa" ;
+	double P0 ;
+		P0:units = "Pa" ;
+	float temp(eta, lat, lon) ;
+		temp:standard_name = "air_temperature" ;
+		temp:units = "K" ;
+
+data:
+ eta = 1 ;
+ eta_bnds =
+  0.5, 1.5 ;
+ A = 1000 ;
+ A_bnds =
+  500, 1500 ;
+ B = 2 ;
+ B_bnds =
+  1, 3 ;
+ PS =
+  2000 ;
+ P0 = 3000 ;
+ temp =
+  300 ;
+}
diff --git a/lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/hybrid_height.ncdump.txt b/lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/hybrid_height.ncdump.txt
new file mode 100644
index 0000000000..e54b16d2c0
--- /dev/null
+++ b/lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/hybrid_height.ncdump.txt
@@ -0,0 +1,72 @@
+netcdf hybrid_height {
+dimensions:
+	model_level_number = UNLIMITED ; // (15 currently)
+	grid_latitude = 100 ;
+	grid_longitude = 100 ;
+	bnds = 2 ;
+variables:
+	float air_potential_temperature(model_level_number, grid_latitude, grid_longitude) ;
+		air_potential_temperature:standard_name = "air_potential_temperature" ;
+		air_potential_temperature:units = "K" ;
+		air_potential_temperature:ukmo__um_stash_source = "m01s00i004" ;
+		air_potential_temperature:source = "Data from Met Office Unified Model 7.04" ;
+		air_potential_temperature:grid_mapping = "rotated_latitude_longitude" ;
+		air_potential_temperature:coordinates = "forecast_period forecast_reference_time level_height sigma surface_altitude time" ;
+	int rotated_latitude_longitude ;
+		rotated_latitude_longitude:grid_mapping_name = "rotated_latitude_longitude" ;
+		rotated_latitude_longitude:longitude_of_prime_meridian = 0. ;
+		rotated_latitude_longitude:semi_major_axis = 6371229. ;
+		rotated_latitude_longitude:semi_minor_axis = 6371229. ;
+		rotated_latitude_longitude:grid_north_pole_latitude = 37.5 ;
+		rotated_latitude_longitude:grid_north_pole_longitude = 177.5 ;
+		rotated_latitude_longitude:north_pole_grid_longitude = 0. ;
+	int model_level_number(model_level_number) ;
+		model_level_number:axis = "Z" ;
+		model_level_number:units = "1" ;
+		model_level_number:standard_name = "model_level_number" ;
+		model_level_number:positive = "up" ;
+	float grid_latitude(grid_latitude) ;
+		grid_latitude:axis = "Y" ;
+		grid_latitude:bounds = "grid_latitude_bnds" ;
+		grid_latitude:units = "degrees" ;
+		grid_latitude:standard_name = "grid_latitude" ;
+	float grid_latitude_bnds(grid_latitude, bnds) ;
+	float grid_longitude(grid_longitude) ;
+		grid_longitude:axis = "X" ;
+		grid_longitude:bounds = "grid_longitude_bnds" ;
+		grid_longitude:units = "degrees" ;
+		grid_longitude:standard_name = "grid_longitude" ;
+	float grid_longitude_bnds(grid_longitude, bnds) ;
+	double forecast_period ;
+		forecast_period:units = "hours" ;
+		forecast_period:standard_name = "forecast_period" ;
+	double forecast_reference_time ;
+		forecast_reference_time:units = "hours since 1970-01-01 00:00:00" ;
+		forecast_reference_time:standard_name = "forecast_reference_time" ;
+		forecast_reference_time:calendar = "gregorian" ;
+	float level_height(model_level_number) ;
+		level_height:bounds = "level_height_bnds" ;
+		level_height:axis = "Z" ;
+		level_height:formula_terms = "a: level_height b: sigma orog: surface_altitude" ;
+		level_height:units = "m" ;
+		level_height:standard_name = "atmosphere_hybrid_height_coordinate" ;
+		level_height:long_name = "level_height" ;
+		level_height:positive = "up" ;
+	float level_height_bnds(model_level_number, bnds) ;
+	float sigma(model_level_number) ;
+		sigma:bounds = "sigma_bnds" ;
+		sigma:units = "1" ;
+		sigma:long_name = "sigma" ;
+	float sigma_bnds(model_level_number, bnds) ;
+	float surface_altitude(grid_latitude, grid_longitude) ;
+		surface_altitude:units = "m" ;
+		surface_altitude:standard_name = "surface_altitude" ;
+		surface_altitude:ukmo__um_stash_source = "m01s00i033" ;
+		surface_altitude:source = "Data from Met Office Unified Model 7.04" ;
+	double time ;
+		time:units = "hours since 1970-01-01 00:00:00" ;
+		time:standard_name = "time" ;
+		time:calendar = "gregorian" ;
+// global attributes:
+		:Conventions = "CF-1.5" ;
+}
diff --git a/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py b/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
new file mode 100644
index 0000000000..fab71cc3c2
--- /dev/null
+++ b/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
@@ -0,0 +1,138 @@
+# Copyright Iris contributors
+#
+# This file is part of Iris and is released under the BSD license.
+# See LICENSE in the root of the repository for full licensing details.
+"""
+Integration tests for loading of data with bounds of derived coordinates.
+
+This includes both out "legacy" sample data (which was wrongly recorded), and more
+modern testdata (which does not load bounds fully prior to Iris 3.13).
+
+Both "legacy" and "newstyle" loading behaviour needs to be tested, which depends on the
+"FUTURE.derived_bounds" flag setting.
+"""
+import numpy as np
+import pytest
+
+import iris
+from iris import sample_data_path, load, FUTURE
+from iris.aux_factory import HybridHeightFactory
+
+from pathlib import Path
+db_testfile_path = Path(__file__).parent / "temp_nc_sources" / "a_new_file.nc"
+legacy_filepath = sample_data_path("hybrid_height.nc")
+
+
+@pytest.fixture(params=[True, False], ids=["with_db", "without_db"])
+def derived_bounds(request):
+    db = request.param
+    with FUTURE.context(derived_bounds=db):
+        yield db
+
+
+
+def test_load_legacy_hh(derived_bounds):
+    cubes = iris.load(legacy_filepath)
+
+    # DEBUG, for now
+    print("")
+    print(cubes)
+    for i_cube, cube in enumerate(cubes):
+        print(f"\n{i_cube:02d} : {cube.summary(shorten=True)}:")
+        print(cube)
+
+    return
+
+    main_cube = cubes.extract_cube("air_potential_temperature")
+    assert len(cubes) == 2
+    (other_cube,) = [cube for cube in cubes if cube != main_cube]
+
+    if not derived_bounds:
+        assert other_cube.name() == "surface_altitude"
+        altitude_coord = main_cube.coord("altitude")
+        assert np.all(other_cube.data == main_cube.coord("surface_altitude").points)
+    else:
+        assert not main_cube.coords("altitude")
+
+        #
+        # OK, confirm for now
+        #   but ***THIS BIT*** is surely wrong ????
+        #
+        assert other_cube.name() == "level_height_bnds"
+
+        # FOR NOW: fix our "problem" by adding a factory "manually"
+        sigma = main_cube.coord("sigma")
+        delta = main_cube.coord("atmosphere_hybrid_height_coordinate")
+        orog = main_cube.coord("surface_altitude")
+        fact = HybridHeightFactory(sigma=sigma, delta=delta, orography=orog)
+        main_cube.add_aux_factory(fact)
+        assert main_cube.coords("altitude")
+        altitude_coord = main_cube.coord("altitude")
+
+    assert altitude_coord.has_bounds()
+    assert altitude_coord.has_lazy_bounds()
+    assert altitude_coord.shape == main_cube.shape
+
+    # Also confirm what we expect from the other factory components (dependencies)
+    (factory,) = main_cube.aux_factories
+    for coord in factory.dependencies.values():
+        assert coord in main_cube.coords()
+        if coord.name() in ("atmosphere_hybrid_height_coordinate", "sigma"):
+            assert coord.has_bounds()
+        else:
+            assert coord.name() == "surface_altitude"
+            assert not coord.has_bounds()
+
+
+def test_load_primary_cf_style(derived_bounds):
+    cubes = iris.load(db_testfile_path)
+    print("")
+    print(cubes)
+    main_cube = cubes.extract_cube("air_temperature")
+    print(main_cube)
+    pressure_coord = main_cube.coord("air_pressure")
+    if not derived_bounds:
+        # We don't expect this case to work "properly"
+        assert not pressure_coord.has_bounds()
+        other_cubes = [cube for cube in cubes if cube != main_cube]
+        assert len(other_cubes) == 3
+        assert set(cube.name() for cube in other_cubes) == {"PS", "A_bnds", "B_bnds"}
+        # TODO: what happened to other 'bits', i.e. a/b (we have only the bounds) ??
+
+    else:
+        assert pressure_coord.has_bounds()
+        # We expect to get an extra "promoted" surface-pressure cube
+        assert len(cubes) == 2
+        (other_cube,) = [cube for cube in cubes if cube != main_cube]
+        assert other_cube.name() == "PS"
+        assert other_cube.units == "Pa"
+        # It should match the reference coord in the main cube
+        assert np.all(main_cube.coord("PS").points == other_cube.data)
+
+        # Check all the dependency coords
+        (factory,) = main_cube.aux_factories
+        assert isinstance(factory, iris.aux_factory.HybridPressureFactory)
+
+        # #
+        # # Fix this : ***something here is wrong*** ???
+        # #
+        # a_coord = main_cube.coord("vertical pressure")
+        # a_coord.rename("a coefficient for vertical coordinate at full levels")
+
+        dep_ids = {
+            dep_name: coord.name()
+            for dep_name, coord in factory.dependencies.items()
+        }
+        assert dep_ids == {
+            #
+            # TODO: ***FIX THIS???*** something seems wrong
+            #
+            # "delta": "a coefficient for vertical coordinate at full levels",
+            "delta": "vertical pressure",
+            "sigma": "b coefficient for vertical coordinate at full levels",
+            "surface_air_pressure": "PS"
+        }
+        for dep_name, coord in factory.dependencies.items():
+            assert coord in main_cube.coords()
+            assert coord.has_bounds() == (dep_name != "surface_air_pressure")
+

From 905493f48053b00256e70454b2cbac911935e226 Mon Sep 17 00:00:00 2001
From: Patrick Peglar <patrick.peglar@metoffice.gov.uk>
Date: Fri, 16 May 2025 14:16:18 +0100
Subject: [PATCH 03/11] Fixes?

---
 lib/iris/fileformats/cf.py                    |  5 +--
 .../derived_bounds/test_bounds_files.py       | 32 ++++++++++---------
 2 files changed, 20 insertions(+), 17 deletions(-)

diff --git a/lib/iris/fileformats/cf.py b/lib/iris/fileformats/cf.py
index 69072f274a..24756deb4e 100644
--- a/lib/iris/fileformats/cf.py
+++ b/lib/iris/fileformats/cf.py
@@ -1459,15 +1459,16 @@ def _translate(self, variables):
 
                     if cf_root not in self.cf_group.bounds:
                         cf_name = cf_var.cf_name
-                        if cf_var.cf_name not in self.cf_group:
+                        if cf_name not in self.cf_group:
                             new_var = CFAuxiliaryCoordinateVariable(
                                 cf_name, cf_var.cf_data
                             )
                             if hasattr(cf_var, "bounds"):
                                 new_var.bounds = cf_var.bounds
-                            new_var.add_formula_term(cf_root, cf_term)
                             self.cf_group[cf_name] = new_var
 
+                        self.cf_group[cf_name].add_formula_term(cf_root, cf_term)
+
         else:
             for cf_var in formula_terms.values():
                 for cf_root, cf_term in cf_var.cf_terms_by_root.items():
diff --git a/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py b/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
index fab71cc3c2..84274e68f6 100644
--- a/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
+++ b/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
@@ -41,9 +41,11 @@ def test_load_legacy_hh(derived_bounds):
         print(f"\n{i_cube:02d} : {cube.summary(shorten=True)}:")
         print(cube)
 
-    return
-
     main_cube = cubes.extract_cube("air_potential_temperature")
+    print("---\nmain cube coords...")
+    for coord in main_cube.coords():
+        print(f'  {coord.var_name.rjust(20)!s} : {coord.summary(shorten=True)}')
+
     assert len(cubes) == 2
     (other_cube,) = [cube for cube in cubes if cube != main_cube]
 
@@ -54,20 +56,20 @@ def test_load_legacy_hh(derived_bounds):
     else:
         assert not main_cube.coords("altitude")
 
+        # #
+        # # OK, confirm for now
+        # #   but ***THIS BIT*** is surely wrong ????
+        # #
+        # assert other_cube.name() == "level_height_bnds"
         #
-        # OK, confirm for now
-        #   but ***THIS BIT*** is surely wrong ????
-        #
-        assert other_cube.name() == "level_height_bnds"
-
-        # FOR NOW: fix our "problem" by adding a factory "manually"
-        sigma = main_cube.coord("sigma")
-        delta = main_cube.coord("atmosphere_hybrid_height_coordinate")
-        orog = main_cube.coord("surface_altitude")
-        fact = HybridHeightFactory(sigma=sigma, delta=delta, orography=orog)
-        main_cube.add_aux_factory(fact)
-        assert main_cube.coords("altitude")
-        altitude_coord = main_cube.coord("altitude")
+        # # FOR NOW: fix our "problem" by adding a factory "manually"
+        # sigma = main_cube.coord("sigma")
+        # delta = main_cube.coord("atmosphere_hybrid_height_coordinate")
+        # orog = main_cube.coord("surface_altitude")
+        # fact = HybridHeightFactory(sigma=sigma, delta=delta, orography=orog)
+        # main_cube.add_aux_factory(fact)
+        # assert main_cube.coords("altitude")
+        # altitude_coord = main_cube.coord("altitude")
 
     assert altitude_coord.has_bounds()
     assert altitude_coord.has_lazy_bounds()

From cc6183ea778b1ad906cb00b4fb3f10b7d270ed4e Mon Sep 17 00:00:00 2001
From: Patrick Peglar <patrick.peglar@metoffice.gov.uk>
Date: Wed, 28 May 2025 17:12:24 +0100
Subject: [PATCH 04/11] Remove links for invalid derived bounds.

---
 lib/iris/fileformats/cf.py | 21 +++++++++++++++++----
 1 file changed, 17 insertions(+), 4 deletions(-)

diff --git a/lib/iris/fileformats/cf.py b/lib/iris/fileformats/cf.py
index 24756deb4e..36b57afba0 100644
--- a/lib/iris/fileformats/cf.py
+++ b/lib/iris/fileformats/cf.py
@@ -1421,9 +1421,11 @@ def _translate(self, variables):
 
         if iris.FUTURE.derived_bounds:
             # cf_var = CFFormulaTermsVariable (loops through everything that appears in formula terms)
+            all_roots = set()
             for cf_var in formula_terms.values():
                 # eg. eta:'a' | cf_root = eta and cf_term = a. cf_var.cf_terms_by_root = {'eta': 'a'} (looking at all appearances in formula terms)
                 for cf_root, cf_term in cf_var.cf_terms_by_root.items():
+                    all_roots.add(cf_root)
                     # gets set to the bounds of the coord from cf_root_coord
                     bounds_name = None
                     # cf_root_coord = CFCoordinateVariable of the coordinate relating to the root
@@ -1469,6 +1471,16 @@ def _translate(self, variables):
 
                         self.cf_group[cf_name].add_formula_term(cf_root, cf_term)
 
+
+            for cf_root in all_roots:
+                # Invalidate "broken" bounds connections
+                root_var = self.cf_group[cf_root]
+                if all(key in root_var.ncattrs() for key in ("bounds", "formula_terms")):
+                    bounds_var = self.cf_group.get(root_var.bounds)
+                    if bounds_var is not None and "formula_terms" not in bounds_var.ncattrs():
+                        # This means it is *not* a valid bounds var
+                        root_var.bounds = None
+
         else:
             for cf_var in formula_terms.values():
                 for cf_root, cf_term in cf_var.cf_terms_by_root.items():
@@ -1555,10 +1567,11 @@ def _span_check(
             if iris.FUTURE.derived_bounds:
                 if hasattr(cf_variable, "bounds"):
                     if cf_variable.bounds not in cf_group:
-                        bounds_var = self.cf_group[cf_variable.bounds]
-                        # TODO: warning if span fails
-                        if bounds_var.spans(cf_variable):
-                            cf_group[cf_variable.bounds] = bounds_var
+                        bounds_var = self.cf_group.get(cf_variable.bounds)
+                        if bounds_var:
+                            # TODO: warning if span fails
+                            if bounds_var.spans(cf_variable):
+                                cf_group[cf_variable.bounds] = bounds_var
 
             # Build CF data variable relationships.
             if isinstance(cf_variable, CFDataVariable):

From c35bdd18d254c24b141613679a84899a2fa939a1 Mon Sep 17 00:00:00 2001
From: Patrick Peglar <patrick.peglar@metoffice.gov.uk>
Date: Wed, 28 May 2025 17:13:08 +0100
Subject: [PATCH 05/11] Updated tests for derived bounds : WIP working but more
 to do.

---
 .../derived_bounds/test_bounds_files.py       | 58 ++++++++-----------
 1 file changed, 25 insertions(+), 33 deletions(-)

diff --git a/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py b/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
index 84274e68f6..cc599bd7c0 100644
--- a/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
+++ b/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
@@ -44,46 +44,33 @@ def test_load_legacy_hh(derived_bounds):
     main_cube = cubes.extract_cube("air_potential_temperature")
     print("---\nmain cube coords...")
     for coord in main_cube.coords():
-        print(f'  {coord.var_name.rjust(20)!s} : {coord.summary(shorten=True)}')
+        var_name = coord.var_name or "-"
+        print(f'  {var_name.rjust(20)!s} : {coord.summary(shorten=True)}')
 
-    assert len(cubes) == 2
-    (other_cube,) = [cube for cube in cubes if cube != main_cube]
-
-    if not derived_bounds:
-        assert other_cube.name() == "surface_altitude"
-        altitude_coord = main_cube.coord("altitude")
-        assert np.all(other_cube.data == main_cube.coord("surface_altitude").points)
+    cube_names = sorted([cube.name() for cube in cubes])
+    if derived_bounds:
+        # get an extra promoted cube for the lost 'level-height bounds"
+        assert cube_names == ["air_potential_temperature", "level_height_bnds", "surface_altitude"]
     else:
-        assert not main_cube.coords("altitude")
-
-        # #
-        # # OK, confirm for now
-        # #   but ***THIS BIT*** is surely wrong ????
-        # #
-        # assert other_cube.name() == "level_height_bnds"
-        #
-        # # FOR NOW: fix our "problem" by adding a factory "manually"
-        # sigma = main_cube.coord("sigma")
-        # delta = main_cube.coord("atmosphere_hybrid_height_coordinate")
-        # orog = main_cube.coord("surface_altitude")
-        # fact = HybridHeightFactory(sigma=sigma, delta=delta, orography=orog)
-        # main_cube.add_aux_factory(fact)
-        # assert main_cube.coords("altitude")
-        # altitude_coord = main_cube.coord("altitude")
+        assert cube_names == ["air_potential_temperature", "surface_altitude"]
 
+    altitude_coord = main_cube.coord("altitude")
     assert altitude_coord.has_bounds()
     assert altitude_coord.has_lazy_bounds()
     assert altitude_coord.shape == main_cube.shape
 
-    # Also confirm what we expect from the other factory components (dependencies)
-    (factory,) = main_cube.aux_factories
-    for coord in factory.dependencies.values():
-        assert coord in main_cube.coords()
-        if coord.name() in ("atmosphere_hybrid_height_coordinate", "sigma"):
-            assert coord.has_bounds()
-        else:
-            assert coord.name() == "surface_altitude"
-            assert not coord.has_bounds()
+    level_height_coord = main_cube.coord("atmosphere_hybrid_height_coordinate")
+    sigma_coord = main_cube.coord("sigma")
+    surface_altitude_coord = main_cube.coord("surface_altitude")
+    assert sigma_coord.has_bounds()
+    assert not surface_altitude_coord.has_bounds()
+
+    if not derived_bounds:
+        other_cube = cubes.extract_cube("surface_altitude")
+        assert np.all(other_cube.data == surface_altitude_coord.points)
+        assert level_height_coord.has_bounds()
+    else:
+        assert not level_height_coord.has_bounds()
 
 
 def test_load_primary_cf_style(derived_bounds):
@@ -92,6 +79,11 @@ def test_load_primary_cf_style(derived_bounds):
     print(cubes)
     main_cube = cubes.extract_cube("air_temperature")
     print(main_cube)
+    print("---\nmain cube coords...")
+    for coord in main_cube.coords():
+        var_name = coord.var_name or "-"
+        print(f'  {var_name.rjust(20)!s} : {coord.summary(shorten=True)}')
+
     pressure_coord = main_cube.coord("air_pressure")
     if not derived_bounds:
         # We don't expect this case to work "properly"

From 917818b793111bcd6947a9285e98b5324bfffc9c Mon Sep 17 00:00:00 2001
From: Patrick Peglar <patrick.peglar@metoffice.gov.uk>
Date: Wed, 28 May 2025 19:31:01 +0100
Subject: [PATCH 06/11] Completed checks, tidied up and removed debug.

---
 .../derived_bounds/test_bounds_files.py       | 131 +++++++++---------
 1 file changed, 63 insertions(+), 68 deletions(-)

diff --git a/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py b/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
index cc599bd7c0..672d993438 100644
--- a/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
+++ b/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
@@ -34,26 +34,15 @@ def derived_bounds(request):
 def test_load_legacy_hh(derived_bounds):
     cubes = iris.load(legacy_filepath)
 
-    # DEBUG, for now
-    print("")
-    print(cubes)
-    for i_cube, cube in enumerate(cubes):
-        print(f"\n{i_cube:02d} : {cube.summary(shorten=True)}:")
-        print(cube)
-
-    main_cube = cubes.extract_cube("air_potential_temperature")
-    print("---\nmain cube coords...")
-    for coord in main_cube.coords():
-        var_name = coord.var_name or "-"
-        print(f'  {var_name.rjust(20)!s} : {coord.summary(shorten=True)}')
-
     cube_names = sorted([cube.name() for cube in cubes])
     if derived_bounds:
         # get an extra promoted cube for the lost 'level-height bounds"
-        assert cube_names == ["air_potential_temperature", "level_height_bnds", "surface_altitude"]
+        expected_cube_names = ["air_potential_temperature", "level_height_bnds", "surface_altitude"]
     else:
-        assert cube_names == ["air_potential_temperature", "surface_altitude"]
+        expected_cube_names = ["air_potential_temperature", "surface_altitude"]
+    assert cube_names == expected_cube_names
 
+    main_cube = cubes.extract_cube("air_potential_temperature")
     altitude_coord = main_cube.coord("altitude")
     assert altitude_coord.has_bounds()
     assert altitude_coord.has_lazy_bounds()
@@ -64,10 +53,10 @@ def test_load_legacy_hh(derived_bounds):
     surface_altitude_coord = main_cube.coord("surface_altitude")
     assert sigma_coord.has_bounds()
     assert not surface_altitude_coord.has_bounds()
+    surface_altitude_cube = cubes.extract_cube("surface_altitude")
+    assert np.all(surface_altitude_cube.data == surface_altitude_coord.points)
 
     if not derived_bounds:
-        other_cube = cubes.extract_cube("surface_altitude")
-        assert np.all(other_cube.data == surface_altitude_coord.points)
         assert level_height_coord.has_bounds()
     else:
         assert not level_height_coord.has_bounds()
@@ -75,58 +64,64 @@ def test_load_legacy_hh(derived_bounds):
 
 def test_load_primary_cf_style(derived_bounds):
     cubes = iris.load(db_testfile_path)
-    print("")
-    print(cubes)
-    main_cube = cubes.extract_cube("air_temperature")
-    print(main_cube)
-    print("---\nmain cube coords...")
-    for coord in main_cube.coords():
-        var_name = coord.var_name or "-"
-        print(f'  {var_name.rjust(20)!s} : {coord.summary(shorten=True)}')
 
+    cube_names = sorted([cube.name() for cube in cubes])
+    if derived_bounds:
+        expected_cube_names = ["PS", "air_temperature"]
+    else:
+        # In this case, the bounds are not properly connected
+        expected_cube_names = ["A_bnds", "B_bnds", "PS", "air_temperature"]
+    assert cube_names == expected_cube_names
+
+    # Check all the coords on the cube, including whether they have bounds
+    main_cube = cubes.extract_cube("air_temperature")
+    assert set(co.name() for co in main_cube.coords()) == set([
+        "a coefficient for vertical coordinate at full levels",
+        "b coefficient for vertical coordinate at full levels",
+        "atmosphere_hybrid_sigma_pressure_coordinate",
+        "vertical pressure",
+        "PS",
+        "air_pressure",
+        "P0",
+    ])
+
+    # First, the main hybrid coord
     pressure_coord = main_cube.coord("air_pressure")
-    if not derived_bounds:
-        # We don't expect this case to work "properly"
-        assert not pressure_coord.has_bounds()
-        other_cubes = [cube for cube in cubes if cube != main_cube]
-        assert len(other_cubes) == 3
-        assert set(cube.name() for cube in other_cubes) == {"PS", "A_bnds", "B_bnds"}
-        # TODO: what happened to other 'bits', i.e. a/b (we have only the bounds) ??
+    assert pressure_coord.has_bounds() == derived_bounds
+    assert pressure_coord.var_name is None
+    assert main_cube.coord_dims(pressure_coord) == (0, 1, 2)
+
+    co_a = main_cube.coord("a coefficient for vertical coordinate at full levels")
+    assert co_a.var_name == "A"
+    assert co_a.has_bounds() == derived_bounds
+    assert main_cube.coord_dims(co_a) == (0,)
+
+    co_b = main_cube.coord("b coefficient for vertical coordinate at full levels")
+    assert co_b.var_name == "B"
+    assert co_b.has_bounds() == derived_bounds
+    assert main_cube.coord_dims(co_b) == (0,)
+
+    co_eta = main_cube.coord("atmosphere_hybrid_sigma_pressure_coordinate")
+    assert co_eta.var_name == "eta"
+    assert co_eta.has_bounds()
+    assert main_cube.coord_dims(co_eta) == (0,)
+
+    # N.B. this coord is 'made up' by the factory, and does *not* come from a variable
+    co_VP = main_cube.coord("vertical pressure")
+    assert co_VP.var_name == "ap"
+    assert co_VP.has_bounds() == derived_bounds
+    assert main_cube.coord_dims(co_VP) == (0,)
+
+    # This is the surface pressure
+    co_PS = main_cube.coord("PS")
+    assert co_PS.var_name == "PS"
+    assert not co_PS.has_bounds()
+    assert main_cube.coord_dims(co_PS) == (1, 2)
+
+    # The scalar reference
+    co_P0 = main_cube.coord("P0")
+    assert co_P0.var_name == "P0"
+    assert not co_P0.has_bounds()
+    assert main_cube.coord_dims(co_P0) == ()
 
-    else:
-        assert pressure_coord.has_bounds()
-        # We expect to get an extra "promoted" surface-pressure cube
-        assert len(cubes) == 2
-        (other_cube,) = [cube for cube in cubes if cube != main_cube]
-        assert other_cube.name() == "PS"
-        assert other_cube.units == "Pa"
-        # It should match the reference coord in the main cube
-        assert np.all(main_cube.coord("PS").points == other_cube.data)
-
-        # Check all the dependency coords
-        (factory,) = main_cube.aux_factories
-        assert isinstance(factory, iris.aux_factory.HybridPressureFactory)
-
-        # #
-        # # Fix this : ***something here is wrong*** ???
-        # #
-        # a_coord = main_cube.coord("vertical pressure")
-        # a_coord.rename("a coefficient for vertical coordinate at full levels")
-
-        dep_ids = {
-            dep_name: coord.name()
-            for dep_name, coord in factory.dependencies.items()
-        }
-        assert dep_ids == {
-            #
-            # TODO: ***FIX THIS???*** something seems wrong
-            #
-            # "delta": "a coefficient for vertical coordinate at full levels",
-            "delta": "vertical pressure",
-            "sigma": "b coefficient for vertical coordinate at full levels",
-            "surface_air_pressure": "PS"
-        }
-        for dep_name, coord in factory.dependencies.items():
-            assert coord in main_cube.coords()
-            assert coord.has_bounds() == (dep_name != "surface_air_pressure")
 

From ee744991a64880bddfc23302cfed301ba13a23f1 Mon Sep 17 00:00:00 2001
From: Patrick Peglar <patrick.peglar@metoffice.gov.uk>
Date: Thu, 29 May 2025 15:15:10 +0100
Subject: [PATCH 07/11] Replace test datafile with CDL in test.

---
 .../temp_nc_sources/a_new_file.nc             | Bin 17045 -> 0 bytes
 .../temp_nc_sources/a_new_file.ncdump.txt     |  47 ------------
 .../temp_nc_sources/hybrid_height.ncdump.txt  |  72 ------------------
 .../derived_bounds/test_bounds_files.py       |  55 ++++++++++++-
 4 files changed, 53 insertions(+), 121 deletions(-)
 delete mode 100644 lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/a_new_file.nc
 delete mode 100644 lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/a_new_file.ncdump.txt
 delete mode 100644 lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/hybrid_height.ncdump.txt

diff --git a/lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/a_new_file.nc b/lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/a_new_file.nc
deleted file mode 100644
index 148d8aa8238e2dd144a6322af0de3e3094cf7e39..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001

literal 17045
zcmeHPeQZ?65ubN=b`Ek71I8whG!L4lR-vv<2qq3R_St9b;uy<Efj=4@&%QIh$iBPY
zoheBZ2vOBEZPQjFB~?`NDV0Q3R98)+T9MLLO{%mgQa@BF>W6Ce18Ab85Va~OYSWqB
z{oR-E3^)c+@)o#xZ+3Qfc6N4l-^?x#ba!?wn7@2JN%H&2M>XmZKku3J&7j8M9r?TP
zPyuf3?i}z0raY@uzBgdVdOoTjes7@Fpx@_v$dad}*DxnUUW?FRDGYo{=qy@De*P2b
z>*!A>Q?Yb59?LL(3ICHORb>qK#L}5~GO=mxnnWx+6dm5QuCb-@#x*0+;hQ!!H#Rpl
z-b5Rij7qT;Xz4B?mgLdEx8h8#>gs~D$66*h%hGLw13+)BFty0OmARJGM(eJmb?<qJ
zs%j+1<MCKJM{*zC><b9C5aVAgIbLgnks0zJvxPKqwNMZRmuom;ApII#XuSJ-c}O6-
z^Tf4U`?g!>q<t%|Q`3vcrg}fs?`JVo`JsIxlNz=HUOLjQ<@2@%xmj6Nc~&+FsYTUl
zR(ABZhr|3Xvr)_nw6C-0woqTVXKP<1+|k|}qQEZ6v#3gy>TAC(#Kfvp8$*=ZHWJU+
z{MWI04tI2R+0pn|EFsg{9vsiwiDWjPwZ}=vox#yqqm9%|fps!n#bN}wok*WD1nWBE
zWAP}rQp3|oOb)S3B%xIef4Vn#%Y~5enN6xO`XEbKo6uNb-*|Zm8xsptOU8nplMegn
zMP=8C{;LI@)8Xp>Lmgf_OC9FLh0_a1dwp8M&s&ySJ1mdZZs<BBUHU$EYgU_Cv#`m)
zn&quL8vVk%MYTL@m*r+{M>3g?#uH9ft@dSRuUUQ3BS5K3q!H-R7Dw+MhLh>B@lhv|
zjitx1dswMj6?!74)ozarCZZX8u+<Lgr&OlZ?ho6kCI(Gnz{Ip%#b$09<;`1fPk0-V
ztj5yedj2PiHMRy|3Q)v8F^_*Z39j&6m_D*GwS05O)?PI{MuO^@^Q2x$TMPnqmg2Q4
z<mn%**5$P6%i~c`mCl}4zbQ5ie?<vB>IWZysNZyVT0MG0I-Q>@E7XsTCF=jV>F9zI
zdf*q9;8&ZDUR8nzpHNuT5B#DMdepD>Xq;C{FH4~$;6Y!j(K$T$F=(m*!$7#um4GV&
zR|2jCTnV@ma3$bMz?Fb20apU91TJd{bPaU(tB*vFaw(X;A~@1&zW%=khd(~+0S!IA
zy3ZymE|x2)t8QSgMA$~TAKkm_C&6vqoy-wC9OPR#+Aw9>o=0dU82s%GifW!v$rm|O
zU8C#BQ~6g--rp7Gx(L01Zvzmm`lz4wQUky0oVB*jSz!nOgOk6!Hn^>~i@AGLvRbw8
zr}t^QhUx^-BtJwiA3wB_d{S4K{oo8<-TR3F+98NriCx)5CqplX>5w3}RdI2wJ@(X{
z)Wr{Uz2+kFl>=X<I|QM;%tW!kClAs~Qp9hmKG*&8trzKap;)7(df&=L|9*4Hzu<9@
z95>Ci1l0J?D;n;iaVr|LA~h)*t0MI%QjH?DC>rp^k1YEBTzAK<?>r^2wWE7W`M_$R
ze>Of_nb497A$+#t)7$;jz5x@JQgnGwpLv}IrUjlV1djh<T&rmL*>}GB#>MQNECg?S
zch#w`a5zL*U@W%PKU2Qu$;C@PM6~O*C-!%BcLsAIXw}{4T8=$?zA?da(LiAV-H&A*
zem=igo-kw4d|nUyG@4B8jwGD1nBq`J8^vi`P`~YD?cwp!QF}DDCpMa~ctMkrsbnUe
zjqgE06{jrZESmgW0_PGkStk*7(owV|?%&k~tvK1SWF|EdOUEK3cMhiG(MTq~d(4TX
z(y>ftoUuc>KnJQU7p#!T0vMNIg^E+yMIj1{*Dc|GFgebF4BcN*SQKM<s=$T>F|r)Q
z7zlNR210!up@<G*)bS{b6Fll<olGptPctV3J90r<LoRL3lqyXJ1RMY{9Z*-87vjS7
z=Z_S=5p6JwJ+Ql*@ucDNywtmCE+n91j2z{f<vAERdq8y9j9<vGRQViCKkvyTL9@;w
zruZ45pUD#NCN4KK9D67k8y=1i#bXIJ2+6d~r-`!hA!n2`5QZ$V6;l7Q#?xn%vCOQM
zSFZ36yj_&}_VBWn`xNFZKVsgsG+w$+lVaWlRtkTxspE-wmZ!9`L{dNIbn}TnJ@QO$
zIFy*)f5=xe^vZh8&=)kI%EPnM438RBrpEAiE^^3s>CtYr`!)>EHBcXBvP75tFyoO|
z396vtyMj?`r_PgKgfa7&@=v_p<Vxz%wt3ViZwB?MAHV)>^~qHwbP_a;HoRigFpz7d
zW648rDynLsk`EE464mFFr@G;AqEX|f{hl;NbIQQe8n%7CO_8?k*D&?4hC~LUoL|4U
zJ!a^RFg(;jZgM5yO2CzXD*;ymt^_Vg3Fuwrhu`NWe~MiOlCZ6<yo(&Uz?=NRS#0vp
zaT8v5u;<1Lvk7l8XIsy{2~S-joA7&nFmY5n`LQF$w+f(E2=S1Mli%i7of9$NgNSSZ
ze&?~b|E`?8RZnbT&)EZ9pq4y`=QYr2eq2i5+Sl7a?fh6rUwnM>WAqn={q)Dz-^`hU
z{o6aOn`r?zQcr&!pBkbajIAf{j)S}D7yNiNy*Ba32bI&Njza%g{UjY$*oiGCPEdj$
zSJ3UjyWXW2__3C5ANbyH_`yvVRRU`K;TV!ZmvPl<!;j)`tIK~bX<fiMSTo%GJ*AT#
z(JH)Zb;?z(haP8K&6xW^tknfkFQaMmk*n1a=U0u!g+5>RGzeTcTM5kI2H&OqBy)E#
zTXaEG3r>GghT%J=PR-22+#U{>5NaV0ayEGTfrpMWxJYT{-il{x?j$ZVsIa#(g?@Ex
zVxRDUyKo4AV=e~R)l;RoxD<WaFg}3rWlJw5!FP>!_^R<<(|F*3in8Dur4+hiY<Q=Y
zvtwTN=i`jJrDEGc;~jOfIB{<I;6d_oEClDxEgzWj$%O|T-?x0I*^fX$j%gpPrPU56
zyyb(*tkxGNay|0Fcv?Yt<SidYty<r$$aU8TTk7c#Pu?OOU`_DBqI&w9-r5Kx?5Zyt
z8%uo<eEE7u*-+so4>Y9|R=&xj-MT$3TA*j&QQq62DKU0%>t3vzJkm2!;J_-e4&^-I
zl8gTLJExv|JiiN9Oe3$dOv(q2D}SVXse(*qxW&xe3zz!7WT79wn0x9Ly17cnuYdC0
za@G@bJn*XQzc1BDph!j%V9qjVG-385|DNrgr(}J|_3npq4(L+t9kdd#w0F<+G~KV~
z{xQqiAYHcjmb`-MH?Q-kp%R4Nm$?h9?wrd_bS<-!ce9D+YVd|t_ZvpvAz9b}U**uw
zHqGhmZu4Pj0z|K6qyucZ#17~iedXfr65oI<oOqgV`;GB`7~o_%%r0N@vPGLd4W$4I
yZ3(2HIp=fD58MJub=W=P!~x#rd_ZA-Ex)(i{k!#--)!ZMsf5`IA$Nba+W$XBgmt_C

diff --git a/lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/a_new_file.ncdump.txt b/lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/a_new_file.ncdump.txt
deleted file mode 100644
index ff9f9f2d4f..0000000000
--- a/lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/a_new_file.ncdump.txt
+++ /dev/null
@@ -1,47 +0,0 @@
-netcdf a_new_file {
-dimensions:
-	eta = 1 ;
-	lat = 1 ;
-	lon = 1 ;
-	bnds = 2 ;
-variables:
-	double eta(eta) ;
-		eta:long_name = "eta at full levels" ;
-		eta:positive = "down" ;
-		eta:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
-		eta:formula_terms = "a: A b: B ps: PS p0: P0" ;
-		eta:bounds = "eta_bnds" ;
-	double eta_bnds(eta, bnds) ;
-		eta_bnds:formula_terms = "a: A_bnds b: B_bnds ps: PS p0: P0" ;
-	double A(eta) ;
-		A:long_name = "a coefficient for vertical coordinate at full levels" ;
-		A:units = "1" ;
-	double A_bnds(eta, bnds) ;
-	double B(eta) ;
-		B:long_name = "b coefficient for vertical coordinate at full levels" ;
-		B:units = "1" ;
-	double B_bnds(eta, bnds) ;
-	double PS(lat, lon) ;
-		PS:units = "Pa" ;
-	double P0 ;
-		P0:units = "Pa" ;
-	float temp(eta, lat, lon) ;
-		temp:standard_name = "air_temperature" ;
-		temp:units = "K" ;
-
-data:
- eta = 1 ;
- eta_bnds =
-  0.5, 1.5 ;
- A = 1000 ;
- A_bnds =
-  500, 1500 ;
- B = 2 ;
- B_bnds =
-  1, 3 ;
- PS =
-  2000 ;
- P0 = 3000 ;
- temp =
-  300 ;
-}
diff --git a/lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/hybrid_height.ncdump.txt b/lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/hybrid_height.ncdump.txt
deleted file mode 100644
index e54b16d2c0..0000000000
--- a/lib/iris/tests/integration/netcdf/derived_bounds/temp_nc_sources/hybrid_height.ncdump.txt
+++ /dev/null
@@ -1,72 +0,0 @@
-netcdf hybrid_height {
-dimensions:
-	model_level_number = UNLIMITED ; // (15 currently)
-	grid_latitude = 100 ;
-	grid_longitude = 100 ;
-	bnds = 2 ;
-variables:
-	float air_potential_temperature(model_level_number, grid_latitude, grid_longitude) ;
-		air_potential_temperature:standard_name = "air_potential_temperature" ;
-		air_potential_temperature:units = "K" ;
-		air_potential_temperature:ukmo__um_stash_source = "m01s00i004" ;
-		air_potential_temperature:source = "Data from Met Office Unified Model 7.04" ;
-		air_potential_temperature:grid_mapping = "rotated_latitude_longitude" ;
-		air_potential_temperature:coordinates = "forecast_period forecast_reference_time level_height sigma surface_altitude time" ;
-	int rotated_latitude_longitude ;
-		rotated_latitude_longitude:grid_mapping_name = "rotated_latitude_longitude" ;
-		rotated_latitude_longitude:longitude_of_prime_meridian = 0. ;
-		rotated_latitude_longitude:semi_major_axis = 6371229. ;
-		rotated_latitude_longitude:semi_minor_axis = 6371229. ;
-		rotated_latitude_longitude:grid_north_pole_latitude = 37.5 ;
-		rotated_latitude_longitude:grid_north_pole_longitude = 177.5 ;
-		rotated_latitude_longitude:north_pole_grid_longitude = 0. ;
-	int model_level_number(model_level_number) ;
-		model_level_number:axis = "Z" ;
-		model_level_number:units = "1" ;
-		model_level_number:standard_name = "model_level_number" ;
-		model_level_number:positive = "up" ;
-	float grid_latitude(grid_latitude) ;
-		grid_latitude:axis = "Y" ;
-		grid_latitude:bounds = "grid_latitude_bnds" ;
-		grid_latitude:units = "degrees" ;
-		grid_latitude:standard_name = "grid_latitude" ;
-	float grid_latitude_bnds(grid_latitude, bnds) ;
-	float grid_longitude(grid_longitude) ;
-		grid_longitude:axis = "X" ;
-		grid_longitude:bounds = "grid_longitude_bnds" ;
-		grid_longitude:units = "degrees" ;
-		grid_longitude:standard_name = "grid_longitude" ;
-	float grid_longitude_bnds(grid_longitude, bnds) ;
-	double forecast_period ;
-		forecast_period:units = "hours" ;
-		forecast_period:standard_name = "forecast_period" ;
-	double forecast_reference_time ;
-		forecast_reference_time:units = "hours since 1970-01-01 00:00:00" ;
-		forecast_reference_time:standard_name = "forecast_reference_time" ;
-		forecast_reference_time:calendar = "gregorian" ;
-	float level_height(model_level_number) ;
-		level_height:bounds = "level_height_bnds" ;
-		level_height:axis = "Z" ;
-		level_height:formula_terms = "a: level_height b: sigma orog: surface_altitude" ;
-		level_height:units = "m" ;
-		level_height:standard_name = "atmosphere_hybrid_height_coordinate" ;
-		level_height:long_name = "level_height" ;
-		level_height:positive = "up" ;
-	float level_height_bnds(model_level_number, bnds) ;
-	float sigma(model_level_number) ;
-		sigma:bounds = "sigma_bnds" ;
-		sigma:units = "1" ;
-		sigma:long_name = "sigma" ;
-	float sigma_bnds(model_level_number, bnds) ;
-	float surface_altitude(grid_latitude, grid_longitude) ;
-		surface_altitude:units = "m" ;
-		surface_altitude:standard_name = "surface_altitude" ;
-		surface_altitude:ukmo__um_stash_source = "m01s00i033" ;
-		surface_altitude:source = "Data from Met Office Unified Model 7.04" ;
-	double time ;
-		time:units = "hours since 1970-01-01 00:00:00" ;
-		time:standard_name = "time" ;
-		time:calendar = "gregorian" ;
-// global attributes:
-		:Conventions = "CF-1.5" ;
-}
diff --git a/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py b/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
index 672d993438..c512e3c98d 100644
--- a/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
+++ b/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
@@ -17,6 +17,7 @@
 import iris
 from iris import sample_data_path, load, FUTURE
 from iris.aux_factory import HybridHeightFactory
+from iris.tests.stock.netcdf import ncgen_from_cdl
 
 from pathlib import Path
 db_testfile_path = Path(__file__).parent / "temp_nc_sources" / "a_new_file.nc"
@@ -29,6 +30,56 @@ def derived_bounds(request):
     with FUTURE.context(derived_bounds=db):
         yield db
 
+@pytest.fixture()
+def cf_primary_sample_path(tmp_path_factory):
+    cdl = """
+        netcdf a_new_file {
+        dimensions:
+            eta = 1 ;
+            lat = 1 ;
+            lon = 1 ;
+            bnds = 2 ;
+        variables:
+            double eta(eta) ;
+                eta:long_name = "eta at full levels" ;
+                eta:positive = "down" ;
+                eta:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
+                eta:formula_terms = "a: A b: B ps: PS p0: P0" ;
+                eta:bounds = "eta_bnds" ;
+            double eta_bnds(eta, bnds) ;
+                eta_bnds:formula_terms = "a: A_bnds b: B_bnds ps: PS p0: P0" ;
+            double A(eta) ;
+                A:long_name = "a coefficient for vertical coordinate at full levels" ;
+                A:units = "1" ;
+            double A_bnds(eta, bnds) ;
+            double B(eta) ;
+                B:long_name = "b coefficient for vertical coordinate at full levels" ;
+                B:units = "1" ;
+            double B_bnds(eta, bnds) ;
+            double PS(lat, lon) ;
+                PS:units = "Pa" ;
+            double P0 ;
+                P0:units = "Pa" ;
+            float temp(eta, lat, lon) ;
+                temp:standard_name = "air_temperature" ;
+                temp:units = "K" ;
+        
+        data:
+         eta = 1 ;
+         eta_bnds = 0.5, 1.5 ;
+         A = 1000 ;
+         A_bnds = 500, 1500 ;
+         B = 2 ;
+         B_bnds = 1, 3 ;
+         PS = 2000 ;
+         P0 = 3000 ;
+         temp = 300 ;
+        }
+    """
+    dirpath = tmp_path_factory.mktemp("tmp")
+    filepath = dirpath / "tmp.nc"
+    ncgen_from_cdl(cdl_str=cdl, cdl_path=None, nc_path=filepath)
+    return filepath
 
 
 def test_load_legacy_hh(derived_bounds):
@@ -62,8 +113,8 @@ def test_load_legacy_hh(derived_bounds):
         assert not level_height_coord.has_bounds()
 
 
-def test_load_primary_cf_style(derived_bounds):
-    cubes = iris.load(db_testfile_path)
+def test_load_primary_cf_style(derived_bounds, cf_primary_sample_path):
+    cubes = iris.load(cf_primary_sample_path)
 
     cube_names = sorted([cube.name() for cube in cubes])
     if derived_bounds:

From ad29c2bac2e6f8be1417806e46f6d28b0bf0fb8b Mon Sep 17 00:00:00 2001
From: Patrick Peglar <patrick.peglar@metoffice.gov.uk>
Date: Thu, 29 May 2025 15:50:55 +0100
Subject: [PATCH 08/11] Remove redundant imports.

---
 .../integration/netcdf/derived_bounds/test_bounds_files.py     | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py b/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
index c512e3c98d..9defc00ba9 100644
--- a/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
+++ b/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
@@ -15,8 +15,7 @@
 import pytest
 
 import iris
-from iris import sample_data_path, load, FUTURE
-from iris.aux_factory import HybridHeightFactory
+from iris import sample_data_path, FUTURE
 from iris.tests.stock.netcdf import ncgen_from_cdl
 
 from pathlib import Path

From b7dd7bce017d2c2f579df5a5f60f17dc95ac23e6 Mon Sep 17 00:00:00 2001
From: Patrick Peglar <patrick.peglar@metoffice.gov.uk>
Date: Thu, 29 May 2025 23:58:33 +0100
Subject: [PATCH 09/11] Remerge split code sections for
 derived-bounds-handling/not.

---
 lib/iris/fileformats/cf.py | 119 ++++++++++++++++---------------------
 1 file changed, 51 insertions(+), 68 deletions(-)

diff --git a/lib/iris/fileformats/cf.py b/lib/iris/fileformats/cf.py
index 36b57afba0..5f50711b4d 100644
--- a/lib/iris/fileformats/cf.py
+++ b/lib/iris/fileformats/cf.py
@@ -1420,27 +1420,35 @@ def _translate(self, variables):
         formula_terms = _CFFormulaTermsVariable.identify(variables)
 
         if iris.FUTURE.derived_bounds:
-            # cf_var = CFFormulaTermsVariable (loops through everything that appears in formula terms)
+            # Keep track of all the root vars so we can unpick invalid bounds
             all_roots = set()
-            for cf_var in formula_terms.values():
-                # eg. eta:'a' | cf_root = eta and cf_term = a. cf_var.cf_terms_by_root = {'eta': 'a'} (looking at all appearances in formula terms)
-                for cf_root, cf_term in cf_var.cf_terms_by_root.items():
+
+        # cf_var = CFFormulaTermsVariable (loops through everything that appears in formula terms)
+        for cf_var in formula_terms.values():
+            # eg. eta:'a' | cf_root = eta and cf_term = a. cf_var.cf_terms_by_root = {'eta': 'a'} (looking at all appearances in formula terms)
+            for cf_root, cf_term in cf_var.cf_terms_by_root.items():
+
+                if iris.FUTURE.derived_bounds:
                     all_roots.add(cf_root)
-                    # gets set to the bounds of the coord from cf_root_coord
-                    bounds_name = None
+
+                    # For the "newstyle" derived-bounds implementation, find vars which appear in derived bounds terms
+                    #  and turn them into bounds vars (though they don't appear in a "bounds" attribute)
+
                     # cf_root_coord = CFCoordinateVariable of the coordinate relating to the root
                     cf_root_coord = self.cf_group.coordinates.get(cf_root)
                     if cf_root_coord is None:
                         cf_root_coord = self.cf_group.auxiliary_coordinates.get(cf_root)
-                    with contextlib.suppress(AttributeError):
-                        # Copes with cf_root_coord not existing, OR not having
-                        #  `bounds` attribute.
-                        bounds_name = cf_root_coord.bounds
 
-                    if bounds_name is not None:
+                    bounds_name = getattr(cf_root_coord, "bounds", None)
+                    # TODO: rename to "root_bounds_name" ?
+                    # N.B. cf_root_coord may here be None, if the root var was not a
+                    #  coord - that is ok, it will not have a 'bounds', we will skip it.
+                    if bounds_name in self.cf_group:
+                        # Found a valid *root* bounds variable : search for a corresponding *term* bounds variable
                         try:
                             # This will error if more or less than 1 variable is found.
                             # TODO: try a try/except here or logical alternative
+                            # TODO: rename as "root_bounds_var" ?
                             (bounds_var,) = [
                                 # loop through all formula terms and add them if they have a cf_term_by_root
                                 # where (bounds of cf_root): cf_term (same as before)
@@ -1459,19 +1467,21 @@ def _translate(self, variables):
                             # Modify the boundary_variable set _to_be_promoted to True
                             self.cf_group.get(bounds_name)._to_be_promoted = True
 
-                    if cf_root not in self.cf_group.bounds:
-                        cf_name = cf_var.cf_name
-                        if cf_name not in self.cf_group:
-                            new_var = CFAuxiliaryCoordinateVariable(
-                                cf_name, cf_var.cf_data
-                            )
-                            if hasattr(cf_var, "bounds"):
-                                new_var.bounds = cf_var.bounds
-                            self.cf_group[cf_name] = new_var
-
-                        self.cf_group[cf_name].add_formula_term(cf_root, cf_term)
+                if cf_root not in self.cf_group.bounds:
+                    # TODO: explain this section ?
+                    cf_name = cf_var.cf_name
+                    if cf_name not in self.cf_group:
+                        new_var = CFAuxiliaryCoordinateVariable(
+                            cf_name, cf_var.cf_data
+                        )
+                        if iris.FUTURE.derived_bounds and hasattr(cf_var, "bounds"):
+                            # Implement "new-style" derived bounds link
+                            new_var.bounds = cf_var.bounds
+                        self.cf_group[cf_name] = new_var
 
+                    self.cf_group[cf_name].add_formula_term(cf_root, cf_term)
 
+        if iris.FUTURE.derived_bounds:
             for cf_root in all_roots:
                 # Invalidate "broken" bounds connections
                 root_var = self.cf_group[cf_root]
@@ -1481,18 +1491,6 @@ def _translate(self, variables):
                         # This means it is *not* a valid bounds var
                         root_var.bounds = None
 
-        else:
-            for cf_var in formula_terms.values():
-                for cf_root, cf_term in cf_var.cf_terms_by_root.items():
-                    # Ignore formula terms owned by a bounds variable.
-                    if cf_root not in self.cf_group.bounds:
-                        cf_name = cf_var.cf_name
-                        if cf_var.cf_name not in self.cf_group:
-                            self.cf_group[cf_name] = CFAuxiliaryCoordinateVariable(
-                                cf_name, cf_var.cf_data
-                            )
-                        self.cf_group[cf_name].add_formula_term(cf_root, cf_term)
-
         # Determine the CF data variables.
         data_variable_names = (
             set(netcdf_variable_names) - self.cf_group.non_data_variable_names
@@ -1565,6 +1563,7 @@ def _span_check(
                     _span_check(cf_name)
 
             if iris.FUTURE.derived_bounds:
+                # TODO: explain this section
                 if hasattr(cf_variable, "bounds"):
                     if cf_variable.bounds not in cf_group:
                         bounds_var = self.cf_group.get(cf_variable.bounds)
@@ -1614,42 +1613,26 @@ def _span_check(
         # Determine whether there are any formula terms that
         # may be promoted to a CFDataVariable and restrict promotion to only
         # those formula terms that are reference surface/phenomenon.
-        if iris.FUTURE.derived_bounds:
-            for cf_var in self.cf_group.formula_terms.values():
+        for cf_var in self.cf_group.formula_terms.values():
+            if iris.FUTURE.derived_bounds:
                 if self.cf_group[cf_var.cf_name] is CFBoundaryVariable:
                     continue
-                else:
-                    for cf_root, cf_term in cf_var.cf_terms_by_root.items():
-                        cf_root_var = self.cf_group[cf_root]
-                        if not hasattr(cf_root_var, "standard_name"):
-                            continue
-                        name = cf_root_var.standard_name or cf_root_var.long_name
-                        terms = reference_terms.get(name, [])
-                        if isinstance(terms, str) or not isinstance(terms, Iterable):
-                            terms = [terms]
-                        cf_var_name = cf_var.cf_name
-                        if (
-                            cf_term in terms
-                            and cf_var_name not in self.cf_group.promoted
-                        ):
-                            data_var = CFDataVariable(cf_var_name, cf_var.cf_data)
-                            self.cf_group.promoted[cf_var_name] = data_var
-                            _build(data_var)
-                            break
-        else:
-            for cf_var in self.cf_group.formula_terms.values():
-                for cf_root, cf_term in cf_var.cf_terms_by_root.items():
-                    cf_root_var = self.cf_group[cf_root]
-                    name = cf_root_var.standard_name or cf_root_var.long_name
-                    terms = reference_terms.get(name, [])
-                    if isinstance(terms, str) or not isinstance(terms, Iterable):
-                        terms = [terms]
-                    cf_var_name = cf_var.cf_name
-                    if cf_term in terms and cf_var_name not in self.cf_group.promoted:
-                        data_var = CFDataVariable(cf_var_name, cf_var.cf_data)
-                        self.cf_group.promoted[cf_var_name] = data_var
-                        _build(data_var)
-                        break
+            for cf_root, cf_term in cf_var.cf_terms_by_root.items():
+                cf_root_var = self.cf_group[cf_root]
+                if iris.FUTURE.derived_bounds:
+                    if not hasattr(cf_root_var, "standard_name"):
+                        continue
+                name = cf_root_var.standard_name or cf_root_var.long_name
+                terms = reference_terms.get(name, [])
+                if isinstance(terms, str) or not isinstance(terms, Iterable):
+                    terms = [terms]
+                cf_var_name = cf_var.cf_name
+                if cf_term in terms and cf_var_name not in self.cf_group.promoted:
+                    data_var = CFDataVariable(cf_var_name, cf_var.cf_data)
+                    self.cf_group.promoted[cf_var_name] = data_var
+                    _build(data_var)
+                    break
+
         # Promote any ignored variables.
         promoted = set()
         not_promoted = ignored.difference(promoted)

From 306151464863c356e638778749f49e1b7f36c708 Mon Sep 17 00:00:00 2001
From: Patrick Peglar <patrick.peglar@metoffice.gov.uk>
Date: Fri, 30 May 2025 11:09:23 +0100
Subject: [PATCH 10/11] Replace try-except; adjust some varnames and comments.

---
 lib/iris/fileformats/cf.py | 53 +++++++++++++++++++-------------------
 1 file changed, 26 insertions(+), 27 deletions(-)

diff --git a/lib/iris/fileformats/cf.py b/lib/iris/fileformats/cf.py
index 5f50711b4d..46b4bdad8a 100644
--- a/lib/iris/fileformats/cf.py
+++ b/lib/iris/fileformats/cf.py
@@ -1420,7 +1420,7 @@ def _translate(self, variables):
         formula_terms = _CFFormulaTermsVariable.identify(variables)
 
         if iris.FUTURE.derived_bounds:
-            # Keep track of all the root vars so we can unpick invalid bounds
+            # Keep track of all the root vars so we can unpick invalid bounds vars
             all_roots = set()
 
         # cf_var = CFFormulaTermsVariable (loops through everything that appears in formula terms)
@@ -1429,43 +1429,42 @@ def _translate(self, variables):
             for cf_root, cf_term in cf_var.cf_terms_by_root.items():
 
                 if iris.FUTURE.derived_bounds:
-                    all_roots.add(cf_root)
-
                     # For the "newstyle" derived-bounds implementation, find vars which appear in derived bounds terms
                     #  and turn them into bounds vars (though they don't appear in a "bounds" attribute)
+                    all_roots.add(cf_root)
 
                     # cf_root_coord = CFCoordinateVariable of the coordinate relating to the root
                     cf_root_coord = self.cf_group.coordinates.get(cf_root)
                     if cf_root_coord is None:
                         cf_root_coord = self.cf_group.auxiliary_coordinates.get(cf_root)
 
-                    bounds_name = getattr(cf_root_coord, "bounds", None)
-                    # TODO: rename to "root_bounds_name" ?
+                    root_bounds_name = getattr(cf_root_coord, "bounds", None)
                     # N.B. cf_root_coord may here be None, if the root var was not a
                     #  coord - that is ok, it will not have a 'bounds', we will skip it.
-                    if bounds_name in self.cf_group:
+                    if root_bounds_name in self.cf_group:
                         # Found a valid *root* bounds variable : search for a corresponding *term* bounds variable
-                        try:
-                            # This will error if more or less than 1 variable is found.
-                            # TODO: try a try/except here or logical alternative
-                            # TODO: rename as "root_bounds_var" ?
-                            (bounds_var,) = [
-                                # loop through all formula terms and add them if they have a cf_term_by_root
-                                # where (bounds of cf_root): cf_term (same as before)
-                                f
-                                for f in formula_terms.values()
-                                if f.cf_terms_by_root.get(bounds_name) == cf_term
-                            ]
-                            if bounds_var != cf_var:
-                                cf_var.bounds = bounds_var.cf_name
+                        term_bounds_vars = [
+                            # loop through all formula terms and add them if they have a cf_term_by_root
+                            # where (bounds of cf_root): cf_term (same as before)
+                            f
+                            for f in formula_terms.values()
+                            if f.cf_terms_by_root.get(root_bounds_name) == cf_term
+                        ]
+                        if len(term_bounds_vars) == 1:
+                            (term_bounds_var,) = term_bounds_vars
+                            if term_bounds_var != cf_var:
+                                # N.B. bounds==main-var is valid CF for *no* bounds
+                                cf_var.bounds = term_bounds_var.cf_name
                                 new_var = CFBoundaryVariable(
-                                    bounds_var.cf_name, bounds_var.cf_data
+                                    term_bounds_var.cf_name, term_bounds_var.cf_data
                                 )
-                                new_var.add_formula_term(bounds_name, cf_term)
-                                self.cf_group[bounds_var.cf_name] = new_var
-                        except ValueError:
+                                new_var.add_formula_term(root_bounds_name, cf_term)
+                                # "Reclassify" this var as a bounds variable
+                                self.cf_group[term_bounds_var.cf_name] = new_var
+                        else:
+                            # Found 0 term bounds, or >1 : discard the term
                             # Modify the boundary_variable set _to_be_promoted to True
-                            self.cf_group.get(bounds_name)._to_be_promoted = True
+                            self.cf_group.get(root_bounds_name)._to_be_promoted = True
 
                 if cf_root not in self.cf_group.bounds:
                     # TODO: explain this section ?
@@ -1486,9 +1485,9 @@ def _translate(self, variables):
                 # Invalidate "broken" bounds connections
                 root_var = self.cf_group[cf_root]
                 if all(key in root_var.ncattrs() for key in ("bounds", "formula_terms")):
-                    bounds_var = self.cf_group.get(root_var.bounds)
-                    if bounds_var is not None and "formula_terms" not in bounds_var.ncattrs():
-                        # This means it is *not* a valid bounds var
+                    root_bounds_var = self.cf_group.get(root_var.bounds)
+                    if root_bounds_var is not None and "formula_terms" not in root_bounds_var.ncattrs():
+                        # This means it is *not* a valid bounds var, according to CF
                         root_var.bounds = None
 
         # Determine the CF data variables.

From 8b6bc0c58bc56fec23ae6bff759c3c0a3aa13c3b Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Fri, 30 May 2025 10:11:00 +0000
Subject: [PATCH 11/11] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 lib/iris/fileformats/cf.py                    | 14 ++++---
 .../netcdf/derived_bounds/__init__.py         |  1 -
 .../derived_bounds/test_bounds_files.py       | 38 +++++++++++--------
 3 files changed, 30 insertions(+), 23 deletions(-)

diff --git a/lib/iris/fileformats/cf.py b/lib/iris/fileformats/cf.py
index 46b4bdad8a..a47c126e9b 100644
--- a/lib/iris/fileformats/cf.py
+++ b/lib/iris/fileformats/cf.py
@@ -1427,7 +1427,6 @@ def _translate(self, variables):
         for cf_var in formula_terms.values():
             # eg. eta:'a' | cf_root = eta and cf_term = a. cf_var.cf_terms_by_root = {'eta': 'a'} (looking at all appearances in formula terms)
             for cf_root, cf_term in cf_var.cf_terms_by_root.items():
-
                 if iris.FUTURE.derived_bounds:
                     # For the "newstyle" derived-bounds implementation, find vars which appear in derived bounds terms
                     #  and turn them into bounds vars (though they don't appear in a "bounds" attribute)
@@ -1470,9 +1469,7 @@ def _translate(self, variables):
                     # TODO: explain this section ?
                     cf_name = cf_var.cf_name
                     if cf_name not in self.cf_group:
-                        new_var = CFAuxiliaryCoordinateVariable(
-                            cf_name, cf_var.cf_data
-                        )
+                        new_var = CFAuxiliaryCoordinateVariable(cf_name, cf_var.cf_data)
                         if iris.FUTURE.derived_bounds and hasattr(cf_var, "bounds"):
                             # Implement "new-style" derived bounds link
                             new_var.bounds = cf_var.bounds
@@ -1484,9 +1481,14 @@ def _translate(self, variables):
             for cf_root in all_roots:
                 # Invalidate "broken" bounds connections
                 root_var = self.cf_group[cf_root]
-                if all(key in root_var.ncattrs() for key in ("bounds", "formula_terms")):
+                if all(
+                    key in root_var.ncattrs() for key in ("bounds", "formula_terms")
+                ):
                     root_bounds_var = self.cf_group.get(root_var.bounds)
-                    if root_bounds_var is not None and "formula_terms" not in root_bounds_var.ncattrs():
+                    if (
+                        root_bounds_var is not None
+                        and "formula_terms" not in root_bounds_var.ncattrs()
+                    ):
                         # This means it is *not* a valid bounds var, according to CF
                         root_var.bounds = None
 
diff --git a/lib/iris/tests/integration/netcdf/derived_bounds/__init__.py b/lib/iris/tests/integration/netcdf/derived_bounds/__init__.py
index 979cb7bf6a..0e88f6563d 100644
--- a/lib/iris/tests/integration/netcdf/derived_bounds/__init__.py
+++ b/lib/iris/tests/integration/netcdf/derived_bounds/__init__.py
@@ -3,4 +3,3 @@
 # This file is part of Iris and is released under the BSD license.
 # See LICENSE in the root of the repository for full licensing details.
 """Integration tests for loading and saving netcdf files with derived coordinate bounds."""
-
diff --git a/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py b/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
index 9defc00ba9..968e74b9b3 100644
--- a/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
+++ b/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py
@@ -2,8 +2,7 @@
 #
 # This file is part of Iris and is released under the BSD license.
 # See LICENSE in the root of the repository for full licensing details.
-"""
-Integration tests for loading of data with bounds of derived coordinates.
+"""Integration tests for loading of data with bounds of derived coordinates.
 
 This includes both out "legacy" sample data (which was wrongly recorded), and more
 modern testdata (which does not load bounds fully prior to Iris 3.13).
@@ -11,14 +10,16 @@
 Both "legacy" and "newstyle" loading behaviour needs to be tested, which depends on the
 "FUTURE.derived_bounds" flag setting.
 """
+
+from pathlib import Path
+
 import numpy as np
 import pytest
 
 import iris
-from iris import sample_data_path, FUTURE
+from iris import FUTURE, sample_data_path
 from iris.tests.stock.netcdf import ncgen_from_cdl
 
-from pathlib import Path
 db_testfile_path = Path(__file__).parent / "temp_nc_sources" / "a_new_file.nc"
 legacy_filepath = sample_data_path("hybrid_height.nc")
 
@@ -29,6 +30,7 @@ def derived_bounds(request):
     with FUTURE.context(derived_bounds=db):
         yield db
 
+
 @pytest.fixture()
 def cf_primary_sample_path(tmp_path_factory):
     cdl = """
@@ -87,7 +89,11 @@ def test_load_legacy_hh(derived_bounds):
     cube_names = sorted([cube.name() for cube in cubes])
     if derived_bounds:
         # get an extra promoted cube for the lost 'level-height bounds"
-        expected_cube_names = ["air_potential_temperature", "level_height_bnds", "surface_altitude"]
+        expected_cube_names = [
+            "air_potential_temperature",
+            "level_height_bnds",
+            "surface_altitude",
+        ]
     else:
         expected_cube_names = ["air_potential_temperature", "surface_altitude"]
     assert cube_names == expected_cube_names
@@ -125,15 +131,17 @@ def test_load_primary_cf_style(derived_bounds, cf_primary_sample_path):
 
     # Check all the coords on the cube, including whether they have bounds
     main_cube = cubes.extract_cube("air_temperature")
-    assert set(co.name() for co in main_cube.coords()) == set([
-        "a coefficient for vertical coordinate at full levels",
-        "b coefficient for vertical coordinate at full levels",
-        "atmosphere_hybrid_sigma_pressure_coordinate",
-        "vertical pressure",
-        "PS",
-        "air_pressure",
-        "P0",
-    ])
+    assert set(co.name() for co in main_cube.coords()) == set(
+        [
+            "a coefficient for vertical coordinate at full levels",
+            "b coefficient for vertical coordinate at full levels",
+            "atmosphere_hybrid_sigma_pressure_coordinate",
+            "vertical pressure",
+            "PS",
+            "air_pressure",
+            "P0",
+        ]
+    )
 
     # First, the main hybrid coord
     pressure_coord = main_cube.coord("air_pressure")
@@ -173,5 +181,3 @@ def test_load_primary_cf_style(derived_bounds, cf_primary_sample_path):
     assert co_P0.var_name == "P0"
     assert not co_P0.has_bounds()
     assert main_cube.coord_dims(co_P0) == ()
-
-