Refactor `fix_forecast_reference_time` and `convert_proleptic`

[atteggiani]

Overview

The fix_forecast_reference_time function currently:

• Ensures the time points are correct and produces time bounds only if the time reference is < 1600 and calendar is Proleptic Gregorian. Time bounds should be produced only for output files for any calendar type.
• Removes forecast_reference_time and forecast_period coordinates if present
• Fails to always produce files with correct time coordinates (example here)

Solution

Refactor the fix_forecast_reference_time and convert_proleptic functions to:

• Convert time units to days
• Remove forecast_reference_time and forecast_period coordinates if present
• Produce time bounds only for output files (no ancil and restarts) for any calendar.
• Always produce files with the correct time coordinate
• Change proleptic gregorian start date to 15/10/1582
• For time having a proleptic gregorian calendar and forecast reference time before the proleptic gregorian start date, convert time points to be consistent with the proleptic gregorian calendar (iris loads them as gregorian times)
• Time units origin will not change (stay whatever iris sets)

[atteggiani]

Overview of the problem

Based on this comment and on the body of this issue, I am trying to always output files having a time coordinate with the correct dimension (number of values and data).

However, it turns out iris is doing something very strange that makes this more complicated than expected.

The example will use the following ancil file, referenced by @blimlim in the comment linked above:

/g/data/vk83/configurations/inputs/access-esm1p5/modern/unused/atmosphere/climatology/global.N96/2020.05.19/qrclim.slt

The referenced ancil is a monthly climatology (between 1979 and 1996) of a soil temperature, and it has one STASH variable with 4 different soil levels.

• iris fails to load this variable in a single cube with a "level" Dimension coordinate, and reads it instead as 4 different cubes
• iris also fails to set the time dimension correctly, loading a cube with a 13-elements pseudo-time dimension instead of 12.

More details

iris is loading the cubes with a 13-elements forecast_reference_time Dimension coordinate (the pseudo-time dimension referenced above) and a 12-elements time Auxiliary coordinate.
As a consequence, the loaded cube (and output file) has a shape of (13 , <nlat> , <nlon>), instead of (12 , <nlat> , <nlon>), as also documented by @blimlim.

The data on the 12-elements time dimension is actually the same as the one on the first 12 elements of the 13-elements forecast_reference_time dimension, with the exception of the very last field.
The data of the very last field is instead the data of the 13th forecast_reference_time point (the last one).

Loaded cube structure

In our case for iris, the time coordinate below is in fact the forecast_reference_time coordinate, with a length of 13.

Cube 1: time lat lon (This corresponds to level 1)
Cube 2: time lat lon (This corresponds to level 2)
Cube 3: time lat lon (This corresponds to level 3)
Cube 4: time lat lon (This corresponds to level 4)

Mule fields structure

For mule, the time coordinate is the "real" time coordinate with a length of 12.
Each UM field is a 2D (lat/lon) slice and the other dimensions are stored with the following structure (in this case there is only one STASH so the structure is simpler):

fld1: lat lon time1 level1
fld2: lat lon time1 level2
fld3: lat lon time1 level3
fld4: lat lon time1 level4

fld5: lat lon time2 level1
fld6: lat lon time2 level2
fld7: lat lon time2 level3
fld8: lat lon time2 level4
...

Since the time coordinate has 12 elements, the total number of fields is 48.

Correspondence between mule field and iris cube

In our case with 4 levels, the correspondence between field and cube can be found in Python with the following:

mule_fields[N] == Cubes[N%4][N//4, ...]

where % is the mod (remainder) operation and // is the integer division.

For example, the 26th field (N = 25) will correspond to Cubes[1][6,...] (25%4 = 1, 25//4 = 6), that is the second Cube, seventh time.

Comparison

I compared the first 12 elements of the loaded iris Cubes data (13-elements "time-dim-coord") with the data of the fields loaded with mule (12-elements "time-dim-coord"):

np.where(
    [
        ~np.all(
            ( field.get_data() - cubes[ifld%4][ifld//4,...].data.data ) == 0
        ) 
        for ifld, field in enumerate(ff.fields)
    ]
)
# (array([47]),)

So the data in the first 12 element of the 13-elements-time cubes is the same as the data of the 12-elements-time coordinate for all fields apart from the very last one.

If you then compare the last field with data from the last cube at the last "13-elements-time" point, they are the same:

np.all(ff.fields[-1].get_data() - cubes[-1][-1,...].data.data) == 0
# True

The bottom line

In this case, iris seems to keep the same exact data for the forecast_reference_time coordinate as the one for the "original" time, with the exception of the very last time point. The very last time point is equal to the last (13th) element of the forecast_reference_time coordinate.

Plan

I don't know what is exactly going on as I didn't check iris source code yet.
My plan is now to dig deeper within iris source code, to check where the issue happens and if we have ways to patch it within um2nc and avoid having to fix this problem in some more messy ways.

@blimlim @MartinDix

[atteggiani]

Based on recent discussion with @MartinDix and @blimlim, the following tasks will be carried out to try and solve the issue in the message above:

• @MartinDix to provide a list of ancillary files to be checked
• @atteggiani to verify if any of the ancillary files in the list produce the same errors of the one referenced in the message above.
• @atteggiani to dig deeper in iris source code and find out why the cube merging (Cubelist.merge method) fails to properly join times correctly for the file referenced in the message above
• @atteggiani to find out if it is possible to use a callback at the time of loading of the iris cube, to fix the time-related errors referenced in the message above

[MartinDix]

Ancillary file collections

/g/data/access/TIDS/UM/ancil/atmos/
/g/data/access/TIDS/CMIP6_ANCIL/data/ancils/
/g/data/access/projects/access/data/ancil/access_cm2_n96e/O1
/g/data/tm70/mrd599/ancil/n96e_momO1_20201102

[MartinDix]

You can remove the forecast_reference_time with something like

def remove_reftime(cube, field, filename):
    cube.remove_coord('forecast_period')
    cube.remove_coord('forecast_reference_time')
cubes = iris.load('/g/data/access/projects/access/umdir/ancil/atmos/GC5/n320e/soil_temp/amip/v1/qrclim.slt', callback=remove_reftime)

Doesn't actually help because iris still thinks there are 13 times for some reason and still won't create a proper time dimension.

With this file xconv thinks there are 3 levels and 24 times so is even more wrong!