Simple data analysis techniques

The major technical difference in the EC-Earth 4 output, compared to EC-Earth 3, is in the output files from the atmosphere. In contrast to EC-Earth 3, OpenIFS writes output via XIOS, the XML-I/O server. The implication of this is that the output is no longer configured via ppt files and that the output files will be in NetCDF format and not GRIB.

The resulting NetCDF files with atmosphere data will technically contain an unstructured grid, representing the reduced Gaussian grid used by OpenIFS. Depending on the grid type (Tco or Tl) and horizontal resolution, the grids will be of different dimensions and cells will map differently to the earth surface.

Processing atmosphere data with CDO

CDO can remap grids from reduced to regular Gaussian grids with the setgridtype,regular command. However, the netCDF files produced by XIOS cannot be processed directly by CDO because they are supposedly unstructured. The trick is to add a proper grid description before applying the setgridtype command. There are 2 ways to achieve this:

(.ECE4) > cdo -L setgridtype,regular -setgrid,<OIFS_INIGG_FILE> output.nc regular_output.nc

where <OIFS_INIGG_FILE> is the ICMGG*INIT file from your simulation. Don’t worry that the initial file is in GRIB, CDO will handle that.

If you have to process many files it could be helpful to prepare a grid description file instead of reading the same GRIB file every time you process an output file. This grid description file is created with

(.ECE4) > cdo griddes <OIFS_INIGG_FILE> > griddes.txt

Then we add this grid description on the fly to process the output.nc file:

(.ECE4) > cdo -L setgridtype,regular -setgrid,griddes.txt output.nc regular_output.nc

The grid description file can be re-used every time you process output on the same grid.

The -L flag added in either method is optional but helps avoiding I/O errors that frequently occur with netCDF4/HDF5 files.

The above method works with linear reduced as well as with cubic orthogonal grids, all the grid information is in the ICMGG*INIT file. For spectral output (spherical harmonics) we need to distinguish between the Tl and Tco case. To process specral fields on the Tl grid one would use

(.ECE4) > cdo -f nc4c -z zip -L sp2gp,linear output regular_output.nc

and for fields on the Tco grid

(.ECE4) > cdo -f nc4c -z zip -L sp2gp,cubic output regular_output.nc

Note that older versions of CDO had sp2gpl which is just a short version of sp2gp,linear, the short form is obsolete and will disappear in the future.

Processing atmosphere data with Iris

The Python Iris package can be used to process OpenIFS data from EC-Earth 4. A simple example is presented here how to load OpenIFS data, make a simple mean calculation, and plot.

The following Python packages are needed for this analysis (the pathlib package is used for convenience):

import iris
import matplotlib.pyplot as plt
from pathlib import Path

A function is defined for loading one or more variables from a given output file, possibly removing coordinates (see below) and applying an Iris constraint.

def load_oifs_data(path, varname, remove_coords=None, new_bounds=None, extract=None):
    cube = iris.load_cube(str(path), varname)
    if remove_coords:
        for c in remove_coords:
            cube.remove_coord(cube.coord(var_name=c))
    if extract:
        cube = cube.extract(extract)
    if new_bounds:
        for c in new_bounds:
            cube.coord(var_name=c).bounds = None
            cube.coord(var_name=c).guess_bounds()
    return cube

Another function helps to plot one or more Iris cubes:

def plot_cubes(*cubes, figsize=(14.5, 4.8)):  # default figsize good for 1x2 plots
    _, axes = plt.subplots(1, len(cubes), figsize=figsize)
    for c, a in zip(cubes, axes if len(cubes)>1 else [axes]):
        lons = c.coord('longitude').points
        lats = c.coord('latitude').points
        s = a.scatter(lons, lats, c=c.data, s=0.7)
        a.legend(*s.legend_elements(num=5), title=c.name())
        a.grid(True)

Some variables for the EC-Earth 4 experiment id and the model output path for OpenIFS (first leg):

expid = 'TUT1'  # ECE4 experiment id
output_dir = Path(f'../run/{expid}/output/001/oifs')

Finally, the data is loaded and plotted. In this example, we load load mean sea level pressure (msl) and (2t) from the monthly output:

msl_monthly = load_oifs_data(
    output_dir / f'{expid}_monthly.nc',
    'msl',
    remove_coords=['time_centered'],
    extract=iris.Constraint(time=lambda cell: cell.point.month==2),  # February
)
tas_monthly = load_oifs_data(
    output_dir / f'{expid}_monthly.nc',
    '2t',
    remove_coords=['time_centered'],
    extract=iris.Constraint(time=lambda cell: cell.point.month==2),  # February
)
plot_cubes(msl_monthly, tas_monthly)

As a simple analysis, the mean sea level pressure is loaded from the daily file as well and the monthly mean is computed and compared to the earlier result:

msl_1d_avg = load_oifs_data(
    output_dir / f'{expid}_daily.nc',
    'msl',
    remove_coords=['time_centered'],
    new_bounds=['time_counter'],
    extract=iris.Constraint(time=lambda cell: cell.point.month==2),
).collapsed('time', iris.analysis.MEAN)
diff = msl_1d_avg - msl_monthly
plot_cubes(msl_1d_avg, diff)
print(
    f'Difference: min: {diff.data.min()}, mean: {diff.data.mean()}, max: {diff.data.max()}'
)