Status update 2025-04-14

I have recently investigated this, and found that the problem relates to multiple orography fields,
It is "cured" by adopting iris.LOAD_POLICY.set("legacy") (as you might expect -- i.e. this gives results equivalent to Iris 3.10).

Unfortunately I haven't (yet) managed to produce a standalone testcase, but am loading a single 500Mb file
(for reference, prodm_op_gl-mn_20250327_06_012.pp)
to reproduce the problem.
This (slightly) reduced testcase has 3 timesteps, and contains an orography field for each timestep, i.e. each with a different timepoint, though in fact the payload content is (unsurprisingly) identical.

In these cases, the "old" Iris simply ignoring the varying orography, producing 3 cubes like this :

>>> print(cubes)
0: mass_fraction_of_cloud_ice_in_air / (kg kg-1) (time: 3; model_level_number: 71; latitude: 1920; longitude: 2560)
1: mass_fraction_of_cloud_liquid_water_in_air / (kg kg-1) (time: 3; model_level_number: 71; latitude: 1920; longitude: 2560)
2: surface_altitude / (m)              (time: 3; latitude: 1920; longitude: 2560)

In these legacy loads, the main data cubes mya have a time dimension, but they never have a time-dependent orography.
But the separate orography (i.e. "surface_altitude") cubes are time-dependent, if that appears in the input.

(
Aside :
I'm surprised that it is not issuing a "multiple references" warning for this in all cases.
... /iris/fileformats/rules.py:44: IrisUserWarning: Multiple reference cubes for orography
I'm seeing that it does do that if I load the whole original set of files, but not for my reduced single-file case.
There's a distinction here that I don't understand.
)

Summary

• We do expect the "newer" Iris (version >=3.11) to load this treating the orography itself as time-dependent.
  • What is unexpected is not that it loads differently, but that it fails.
• This looks like a hitherto unknown concatenate bug. Which we need to investigate
• UsingLOAD_POLICY.set("legacy") cures the problem in these cases,
  • this may in fact be correct + appropriate for these input cases where, in effect, there is duplication in the orography input fields which should not be respected
  • there may be a case here for altering the default loading behaviour, where we knowingly adopted a backwards-incompatible approach as "more correct".

@pp-mo has explained to me an opportunity for improving chunk awareness in hybrid height coordinates. Regardless of the problems with prodm_op_gl-mn_20250327_06_012.pp, this may become necessary when we start seeing more time-dependent hybrid heights due to over-large chunk sizes.

My loose understanding

• The way 'aux factory' coordinates are constructed means they are NOT chunked.
• This is usually workable since their maximum size is 72 * 1920 * 2560 (height * lat * lon).
• Adding a time dimension now risks the chunk being larger than available memory.
• When we create the final array for the hybrid height coordinate, we need chunking to be calculated by _optimum_chunksize:
  

  iris/lib/iris/_lazy_data.py

  Line 216 in 1285efb

  def _optimum_chunksize(

• Final array creation happens somewhere around here:
  

  iris/lib/iris/aux_factory.py

  Lines 486 to 492 in 1285efb

  	# Build the points array
  	nd_points_by_key = self._remap(dependency_dims, derived_dims)
  	points = self._derive(
  	nd_points_by_key["pressure_at_top"],
  	nd_points_by_key["sigma"],
  	nd_points_by_key["surface_air_pressure"],
  	)

• Unfortunately array creation is repeated for every type of aux factory, and must also be implemented again if the user defines their own aux factory, so fixing the chunking will need some refactoring so that we can use common code.

WIP Status update 2025-06-04

Fix prototype is currently up on branch "concatmem_prob" : see changes here

• test/exercise in /home/users/patrick.peglar/worklog/concat_bug/sample_concat_fail.py
• key change is in _aux_factory.py
  • rework _derive to call a new _inner_calc
  • LATEST made switchable
  • TODO make this a generic solution shared across all factory coord types
• currently mixed with various debug-related changes
  • _concatenate.py : array hashing debug (including more human-friendly array hash naming)
  • _data_manager.py : no-touch realisation control in
  • _cube.py : make derived coords tuple from list instead of iterator
• also a separate fix proposed in rules.py/_deferefence_args() --> use core_data so reference targets can be lazy

NOTE : although this problem was seen as potentially due to the change in loading policy introducing a "concatenate" step, in practice I don't now think this is the key difference : It is just that the chunks got bigger since the derived coordinate now varies on the time dimension as well - i.e. it is the time-varying nature that matters.
Meantime, I think that the one-chunk nature of the derived coord has always been there, but I think it only just got bad enough to notice.

ASIDE : there is a lesson here about dask chunking in general.
E.G. we have y * x where they broadcast together as (1, N) and (M, 1) and are chunked as (1, N) and (M, 1) -- both quite small individually (e.g. dim coords, single chunks ).
But the result can then end up with chunks of (M, N), i.e. a size of M*N, which may be rather huge.
So we need to rechunk to get, say, (Q, N) for manageable size : and this will requires rechunking x from (M, 1) to (Q, 1).