Incorrect hybrid height levels when reading UM fields that also have pseudolevels

At cf-python v3.13.0, When pseudo levels and atmosphere hybrid height coordinates are both present in UM/PP file fields, then the vertical coordinates are not mapped correctly to cf-python constructs:

>>> import cf
>>> x = cf.read('atmosa.pa19810901_00', select="stash_code=2530")[0]
>>> print(x)
Field: id%UM_m01s02i530_vn1108 (ncvar%UM_m01s02i530_vn1108)
-----------------------------------------------------------
Data            : id%UM_m01s02i530_vn1108(id%UM_pseudolevel(6), atmosphere_hybrid_height_coordinate(38), time(2), latitude(72), longitude(96))
Cell methods    : time(2): point
Dimension coords: time(2) = [1981-09-01 00:30:00, 1981-09-01 02:00:00] 360_day
                : atmosphere_hybrid_height_coordinate(38) = [0.4000157022846824, ..., 0.9639196605247963]
                : latitude(72) = [-88.75, ..., 88.75] degrees_north
                : longitude(96) = [1.875, ..., 358.125] degrees_east
                : id%UM_pseudolevel(6) = [1, ..., 6] 1
Auxiliary coords: model_level_number(atmosphere_hybrid_height_coordinate(38)) = [1, ..., 29] 1
Coord references: standard_name:atmosphere_hybrid_height_coordinate
Domain ancils   : id%UM_atmosphere_hybrid_height_coordinate_a(atmosphere_hybrid_height_coordinate(38)) = [20.000337706971997, ..., 16875.311437270288] m
                : id%UM_atmosphere_hybrid_height_coordinate_b(atmosphere_hybrid_height_coordinate(38)) = [0.9977164621117602, ..., 0.0013017908966459443] 1
>>> x.coordinate('model_level_number').array
array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 37, 11, 12, 38, 13, 14, 15,
       16, 36, 17, 18, 19, 20, 35, 21, 22, 23, 34, 24, 25, 33, 26, 32, 27,
       31, 28, 30, 29])
# WRONG! Expected
# array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
#        18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
#        35, 36, 37, 38])
>>> x.coordinate('atmosphere_hybrid_height_coordinate').array
array([0.4000157 , 0.61539391, 0.72000126, 0.78048944, 0.81967245,
       0.84705866, 0.8672563 , 0.88275791, 0.895027  , 0.90497666,
       0.91205889, 0.91320829, 0.92012832, 0.92521266, 0.9260267 ,
       0.9311166 , 0.93555108, 0.9394496 , 0.94061245, 0.94290374,
       0.94598535, 0.94875156, 0.95124885, 0.95277226, 0.95351337,
       0.95557733, 0.95746588, 0.95836068, 0.95920047, 0.96079918,
       0.96108679, 0.96225022, 0.9624934 , 0.96320978, 0.96327341,
       0.96374978, 0.96375322, 0.96391966])
# WRONG! Expected:
# array([4.71395848e-04, 1.88558339e-03, 4.24247011e-03, 7.54224105e-03,
#        1.17847112e-02, 1.69699730e-02, 2.30980265e-02, 3.01687791e-02,
#        3.81824160e-02, 4.71387521e-02, 5.70379724e-02, 6.78798919e-02,
#        7.96645105e-02, 9.23920133e-02, 1.06062308e-01, 1.20675302e-01,
#        1.36231087e-01, 1.52729664e-01, 1.70171033e-01, 1.88555193e-01,
#        2.07882053e-01, 2.28151705e-01, 2.49364148e-01, 2.71519382e-01,
#        2.94617409e-01, 3.18658227e-01, 3.43660619e-01, 3.69913158e-01,
#        3.97740871e-01, 4.27516435e-01, 4.59674144e-01, 4.94726010e-01,
#        5.33283045e-01, 5.76466976e-01, 6.26560302e-01, 6.88676026e-01,
#        7.75637967e-01, 9.25212656e-01])

Why

This is because the underlying C library that does a pre-aggregation (prior to the cf.aggregate process) can only reliably do so in two dimensions, not three.

Fix

Whilst updating the C library would be ideal, it might be a major job, and this fix is needed soon. A work-around is to better control the creation of coordinate and coordinate reference constructs in umread.py, so that cf.aggregate doesn't get confused.

Thanks to @theabro for persevering with tracking this down.