Expand ICON-XPP fix for new model version (#3014)

Author: schlunma

doc/changelog.rst

@@ -201,6 +201,7 @@ Fixes for datasets
 -  Fix ERA5 native6 fix to handle single monthly-averaged NetCDF files. (:pull:`2512`) by :user:`rbeucher`
 -  Expand ICON extra facets (:pull:`2965`) by :user:`schlunma`
 -  Copy fixes for obs4MIPs dataset SSMI RSSv07r00 to RSS-v7 (:pull:`2968`) by :user:`bouweandela`
+-  Expand ICON-XPP fix for new model version (:pull:`3014`) by :user:`schlunma`
 
 Installation
 ~~~~~~~~~~~~

esmvalcore/cmor/_fixes/icon/_base_fixes.py

@@ -584,13 +584,8 @@ def fix_metadata(self, cubes: CubeList) -> CubeList:
         if cube.coords(long_name="depth_below_sea"):
             self._fix_depth(cube, cubes)
 
-        # Remove undesired lev coordinate with length 1
-        lev_coord = DimCoord(0.0, var_name="lev")
-        if cube.coords(lev_coord, dim_coords=True):
-            slices = [slice(None)] * cube.ndim
-            slices[cube.coord_dims(lev_coord)[0]] = 0
-            cube = cube[tuple(slices)]
-            cube.remove_coord(lev_coord)
+        # Remove undesired coordinates of length 1
+        cube = self._remove_undesired_coords(cube)
 
         # Fix latitude
         if self.vardef.has_coord_with_standard_name("latitude"):
@@ -795,19 +790,33 @@ def _fix_height(self, cube: Cube, cubes: CubeList) -> Cube:
 
     def _fix_depth(self, cube: Cube, cubes: CubeList) -> None:
         """Fix ocean depth coordinate."""
-        depth_coord = self.fix_depth_coord_metadata(cube)
+        depth_coord = cube.coord(long_name="depth_below_sea")
+        depth_var_name = depth_coord.var_name
+        depth_coord = self.fix_depth_coord_metadata(cube, coord=depth_coord)
         if depth_coord.has_bounds():
             return
 
-        # Try to get bounds of depth coordinate from depth_2 coordinate that
-        # might be present in other variables loaded from the same file
+        # ICON ocean output usually consists of two different depth
+        # coordinates, "depth" and "depth_2". One of them is the depth at the
+        # cell centers, the other at the cell boundaries. The nomenclature is
+        # not consistent on which coordinate is which. Thus, if possible, we
+        # try to add the cell boundaries from the other depth coordinate here
+        # (potentially available from other variables loaded from same file).
+        if depth_var_name == "depth":
+            other_depth_var_name = "depth_2"
+        else:
+            other_depth_var_name = "depth"
         for other_cube in cubes:
-            if not other_cube.coords(var_name="depth_2"):
+            if not other_cube.coords(var_name=other_depth_var_name):
                 continue
-            depth_2_coord = other_cube.coord(var_name="depth_2")
+            depth_2_coord = other_cube.coord(var_name=other_depth_var_name)
             depth_2_coord.convert_units(depth_coord.units)
             bounds = depth_2_coord.core_points()
-            depth_coord.bounds = np.stack((bounds[:-1], bounds[1:]), axis=-1)
+            if bounds.shape == (depth_coord.shape[0] + 1,):
+                depth_coord.bounds = np.stack(
+                    (bounds[:-1], bounds[1:]),
+                    axis=-1,
+                )
 
     def _fix_lat(self, cube: Cube) -> tuple[int, ...]:
         """Fix latitude coordinate of cube."""
@@ -1121,6 +1130,20 @@ def _fix_invalid_time_units(time_coord: Coord) -> None:
         time_coord.points = new_dt_points
         time_coord.units = new_t_units
 
+    def _remove_undesired_coords(self, cube: Cube) -> Cube:
+        """Remove undesired coordinates of length 1 from cube."""
+        coords_to_remove = [
+            DimCoord(0.0, var_name="lev"),
+            DimCoord(0.0, var_name="ice_class"),
+        ]
+        for coord in coords_to_remove:
+            if cube.coords(coord, dim_coords=True):
+                slices: list[int | slice] = [slice(None)] * cube.ndim
+                slices[cube.coord_dims(coord)[0]] = 0
+                cube = cube[tuple(slices)]
+                cube.remove_coord(coord)
+        return cube
+
 
 class NegateData(IconFix):
     """Base fix to negate data."""

esmvalcore/config/configurations/defaults/extra_facets_icon.yml

@@ -138,11 +138,12 @@ projects:
           siconc: {raw_name: conc, var_type: oce_ice, raw_units: '1'}
           siconca: {raw_name: fr_seaice, var_type: atm_2d_ml}
           sithick: {raw_name: hi, var_type: oce_ice}
+          sos: {var_type: oce_def, raw_units: "0.001"}
           tos: {raw_name: t_seasfc, var_type: atm_2d_ml}
           zos: {raw_name: zos, var_type: oce_dbg}
 
           # 3D ocean variables
-          so: {raw_name: so, var_type: oce_def, raw_units: "0.001"}
+          so: {var_type: oce_def, raw_units: "0.001"}
           thetao: {raw_name: to, var_type: oce_def, raw_units: degC}
           uo: {raw_name: u, var_type: oce_def}
           vo: {raw_name: v, var_type: oce_def}

tests/integration/cmor/_fixes/icon/test_icon_xpp.py

@@ -929,7 +929,8 @@ def test_rtmt_fix(cubes_regular_grid):
     np.testing.assert_allclose(cube.data, [[[0.0, 2.0], [4.0, 6.0]]])
 
 
-# Test siconc (for extra_facets, removal of lev coord and  typesi coordinate)
+# Test siconc (for extra_facets, removal of lev/ice_class coord and typesi
+# coordinate)
 
 
 def test_get_siconc_fix():
@@ -939,18 +940,19 @@ def test_get_siconc_fix():
 
 
 @pytest.mark.online
-def test_siconc_fix(cubes_ocean_3d, session):
+@pytest.mark.parametrize("z_coord_name", ["lev", "ice_class"])
+def test_siconc_fix(z_coord_name, cubes_ocean_3d, session):
     """Test fix."""
     cubes = CubeList(
         [cubes_ocean_3d.extract_cube(NameConstraint(var_name="to")).copy()],
     )
     cubes[0].var_name = "conc"
     cubes[0].units = None
 
-    # Add lev coord to test removal of it
+    # Add Z-coord to test removal of it
     cubes[0] = cubes[0][:, [0], :]
     cubes[0].remove_coord("depth")
-    cubes[0].add_dim_coord(DimCoord(0.0, var_name="lev"), 1)
+    cubes[0].add_dim_coord(DimCoord(0.0, var_name=z_coord_name), 1)
 
     fix = get_allvars_fix("SImon", "siconc", session=session)
     fixed_cubes = fix.fix_metadata(cubes)
@@ -965,7 +967,7 @@ def test_siconc_fix(cubes_ocean_3d, session):
     check_typesi(cube)
 
     assert cube.shape == (1, 8)
-    assert not cube.coords(var_name="lev")
+    assert not cube.coords(var_name=z_coord_name)
 
     assert cube.dtype == np.float32
     np.testing.assert_allclose(
@@ -1109,6 +1111,34 @@ def test_thetao_fix(cubes_ocean_3d, session):
     assert cube.shape == (1, 47, 8)
 
 
+@pytest.mark.online
+def test_thetao_fix_switched_depth_coords(cubes_ocean_3d, session):
+    """Test fix."""
+    to_cube = cubes_ocean_3d.extract_cube(NameConstraint(var_name="to"))
+    w_cube = cubes_ocean_3d.extract_cube(NameConstraint(var_name="w"))
+    to_cube.coord("depth").var_name = "depth_2"
+    w_cube.coord("depth").var_name = "depth"
+    cubes = CubeList([to_cube, w_cube])
+
+    fix = get_allvars_fix("Omon", "thetao", session=session)
+
+    fixed_cubes = fix.fix_metadata(cubes)
+
+    assert len(fixed_cubes) == 1
+    cube = fixed_cubes[0]
+    assert cube.var_name == "thetao"
+    assert cube.standard_name == "sea_water_potential_temperature"
+    assert cube.long_name == "Sea Water Potential Temperature"
+    assert cube.units == "degC"
+    assert "positive" not in cube.attributes
+
+    depth_coord = cube.coord("depth")
+    assert depth_coord.has_bounds()
+
+    assert cube.dtype == np.float32
+    assert cube.shape == (1, 47, 8)
+
+
 @pytest.mark.online
 def test_thetao_fix_already_bounds(cubes_ocean_3d, session):
     """Test fix."""
@@ -1164,6 +1194,36 @@ def test_thetao_fix_no_bounds(cubes_ocean_3d, session):
     assert cube.shape == (1, 47, 8)
 
 
+@pytest.mark.online
+def test_thetao_fix_no_bounds_invalid_depth_2(cubes_ocean_3d, session):
+    """Test fix."""
+    to_cube = cubes_ocean_3d.extract_cube(NameConstraint(var_name="to"))
+    w_cube = cubes_ocean_3d.extract_cube(NameConstraint(var_name="w"))[
+        :,
+        :40,
+        :,
+    ]
+    cubes = CubeList([to_cube, w_cube])
+
+    fix = get_allvars_fix("Omon", "thetao", session=session)
+
+    fixed_cubes = fix.fix_metadata(cubes)
+
+    assert len(fixed_cubes) == 1
+    cube = fixed_cubes[0]
+    assert cube.var_name == "thetao"
+    assert cube.standard_name == "sea_water_potential_temperature"
+    assert cube.long_name == "Sea Water Potential Temperature"
+    assert cube.units == "degC"
+    assert "positive" not in cube.attributes
+
+    depth_coord = cube.coord("depth")
+    assert not depth_coord.has_bounds()
+
+    assert cube.dtype == np.float32
+    assert cube.shape == (1, 47, 8)
+
+
 # Test zg (for extra fix)
 
 

tests/unit/test_dataset.py

@@ -2242,15 +2242,15 @@ def test_get_extra_facets_icon_xpp():
     dataset = Dataset(
         project="ICON",
         mip="Omon",
-        short_name="so",
+        short_name="thetao",
         dataset="ICON-XPP",
     )
 
     extra_facets = dataset._get_extra_facets()
 
     assert extra_facets == {
-        "raw_name": "so",
-        "raw_units": "0.001",
+        "raw_name": "to",
+        "raw_units": "degC",
         "var_type": "oce_def",
     }