Combining UGRID fields erroneously creates an extra axis and broadcasts over it

At v3.16.2, combining UGRID fields (e.g. x = f + g) creates an erroneous axis and broadcasts over it:

>>> import cf
>>> f = cf.example_field(8)
>>> f
<CF Field: air_temperature(time(2), ncdim%nMesh2_face(3)) K>
>>> g = f + f
>>> g  # WRONG
<CF Field: air_temperature(time(2), ncdim%nMesh2_face(3), key%domainaxis3(3)) K>

What it should do is:

>>> g = f + f
>>> g  # CORRECT
<CF Field: air_temperature(time(2), ncdim%nMesh2_face(3)) K>

The reason for this is that the code which works out if combining fields is OK and how do to it, both in a metadata-aware fashion, was never designed to cope with the case of physical nature of a discrete axis is defined by auxiliary coordinates (as opposed to dimension coordinates).

>>> cf.environment(paths=False)
Platform: Linux-5.15.0-112-generic-x86_64-with-glibc2.35
HDF5 library: 1.12.2
netcdf library: 4.9.3-development
udunits2 library: /home/david/miniconda3/lib/libudunits2.so.0
esmpy/ESMF: 8.6.1
Python: 3.12.2
dask: 2024.5.0
netCDF4: 1.6.5
psutil: 5.9.8
packaging: 23.2
numpy: 1.26.4
scipy: 1.13.0
matplotlib: 3.8.4
cftime: 1.6.3
cfunits: 3.3.7
cfplot: 3.3.0
cfdm: 1.11.1.0
cf: 3.16.2