Add support for variables phcint and amoc to ICON-XPP CMORizer (#3025)

Co-authored-by: Felicity Chun <32269066+flicj191@users.noreply.github.com>

Author: schlunma

esmvalcore/cmor/_fixes/icon/icon_xpp.py

@@ -2,6 +2,7 @@
 
 import logging
 
+from iris.coords import AuxCoord
 from iris.cube import CubeList
 from scipy import constants
 
@@ -56,6 +57,51 @@ def fix_metadata(self, cubes: CubeList) -> CubeList:
 Hfss = NegateData
 
 
+class Msftmz(IconFix):
+    """Fixes for ``msftmz``."""
+
+    def fix_metadata(self, cubes: CubeList) -> CubeList:
+        """Fix metadata."""
+        preprocessed_cubes = CubeList([])
+        basin_coord = AuxCoord(
+            "placeholder",
+            standard_name="region",
+            long_name="ocean basin",
+            var_name="basin",
+        )
+        var_names = {
+            "atlantic_moc": "atlantic_arctic_ocean",
+            "pacific_moc": "indian_pacific_ocean",
+            "global_moc": "global_ocean",
+        }
+        for var_name, basin in var_names.items():
+            cube = self.get_cube(cubes, var_name=var_name)
+            cube.var_name = "msftmz"
+            cube.long_name = None
+            cube.attributes.locals = {}
+
+            # Remove longitude coordinate (with length 1)
+            cube = cube[..., 0]
+            cube.remove_coord("longitude")
+
+            # Add scalar basin coordinate
+            cube.add_aux_coord(basin_coord.copy(basin), ())
+            preprocessed_cubes.append(cube)
+
+        msftmz_cube = preprocessed_cubes.merge_cube()
+
+        # Swap time and basin coordinates
+        msftmz_cube.transpose([1, 0, 2, 3])
+
+        # By default, merge_cube() sorts the coordinate alphabetically (i.e.,
+        # atlantic_arctic_ocean -> global_ocean -> indian_pacific_ocean). Thus,
+        # we need to restore the desired order (atlantic_arctic_ocean ->
+        # indian_pacific_ocean -> global_ocean).
+        msftmz_cube = msftmz_cube[:, [0, 2, 1], ...]
+
+        return CubeList([msftmz_cube])
+
+
 Rlut = NegateData
 
 

esmvalcore/cmor/tables/cmip5-custom/CMOR_phcint.dat

@@ -0,0 +1,22 @@
+SOURCE: CMIP7
+generic_levels:    olevel
+!============
+variable_entry:    phcint
+!============
+modeling_realm:    ocean
+!----------------------------------
+! Variable attributes:
+!----------------------------------
+standard_name:     integral_wrt_depth_of_sea_water_potential_temperature_expressed_as_heat_content
+units:             J m-2
+cell_methods:      area: time: mean where sea
+cell_measures:     area: areacello
+long_name:         Integrated Ocean Heat Content from Potential Temperature
+comment:           This is the vertically-integrated heat content derived from potential temperature (thetao).
+!----------------------------------
+! Additional variable information:
+!----------------------------------
+dimensions:        longitude latitude time olevel
+type:              real
+!----------------------------------
+!

esmvalcore/cmor/tables/cmip6-custom/CMIP6_custom.json

@@ -767,6 +767,18 @@
             "type": "real",
             "units": "kg s-1"
         },
+        "phcint": {
+            "cell_measures": "area: areacello",
+            "cell_methods": "area: time: mean where sea",
+            "comment": "This is the vertically-integrated heat content derived from potential temperature (thetao).",
+            "dimensions": "longitude latitude time olevel",
+            "long_name": "Integrated Ocean Heat Content from Potential Temperature",
+            "modeling_realm": "ocean",
+            "out_name": "phcint",
+            "standard_name": "integral_wrt_depth_of_sea_water_potential_temperature_expressed_as_heat_content",
+            "type": "real",
+            "units": "J m-2"
+        },
         "ptype": {
             "cell_measures": "area: areacella",
             "cell_methods": "time: mean",

esmvalcore/config/configurations/defaults/extra_facets_icon.yml

@@ -132,6 +132,9 @@ projects:
           gpp: {raw_name: assimi_gross_assimilation_box, var_type: jsb_2d_ml}
           lai: {raw_name: pheno_lai_box, var_type: jsb_2d_ml, raw_units: '1'}
 
+          # 1D ocean variables
+          amoc: {var_type: oce_moc}
+
           # 2D ocean variables
           hfds: {raw_name: HeatFlux_Total, var_type: oce_dbg}
           mlotst: {raw_name: mld, var_type: oce_dbg}
@@ -143,6 +146,8 @@ projects:
           zos: {raw_name: zos, var_type: oce_dbg}
 
           # 3D ocean variables
+          msftmz: {var_type: oce_moc, lat_var: lat}
+          phcint: {var_type: oce_def}
           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}

esmvalcore/preprocessor/_derive/amoc.py

@@ -21,6 +21,8 @@ def required(project):
                 {"short_name": "msftmz", "optional": True},
                 {"short_name": "msftyz", "optional": True},
             ]
+        elif project == "ICON":
+            required = [{"short_name": "msftmz", "mip": "Omon"}]
         else:
             msg = f"Project {project} can not be used for Amoc derivation."
             raise ValueError(msg)

esmvalcore/preprocessor/_derive/phcint.py

@@ -0,0 +1,66 @@
+"""Derivation of variable ``phcint``."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from cf_units import Unit
+from iris import NameConstraint
+
+from esmvalcore.preprocessor._shared import get_coord_weights
+
+from ._baseclass import DerivedVariableBase
+
+if TYPE_CHECKING:
+    from iris.cube import Cube, CubeList
+
+    from esmvalcore.typing import Facets
+
+RHO_CP = 4.09169e6
+RHO_CP_UNIT = Unit("kg m-3 J kg-1 K-1")
+
+
+class DerivedVariable(DerivedVariableBase):
+    """Derivation of variable `ohc`."""
+
+    @staticmethod
+    def required(project: str) -> list[Facets]:  # noqa: ARG004
+        """Declare the variables needed for derivation."""
+        return [{"short_name": "thetao"}]
+
+    @staticmethod
+    def calculate(cubes: CubeList) -> Cube:
+        """Compute vertically-integrated heat content.
+
+        Use c_p * rho_0 = 4.09169e+6 J m-3 K-1 (Kuhlbrodt et al., 2015, Clim.
+        Dyn.)
+
+        Arguments
+        ---------
+        cubes:
+           Input cubes.
+
+        Returns
+        -------
+        :
+            Output cube.
+
+        """
+        cube = cubes.extract_cube(NameConstraint(var_name="thetao"))
+        cube.convert_units("K")
+
+        # In the following, we modify the cube's data and units instead of the
+        # cube directly to avoid dropping cell measures and ancillary variables
+        # (https://scitools-iris.readthedocs.io/en/stable/further_topics/lenient_maths.html#finer-detail)
+
+        # Multiply by c_p * rho_0 -> J m-3
+        cube.data = cube.core_data() * RHO_CP
+        cube.units *= RHO_CP_UNIT
+
+        # Multiply by layer depth -> J m-2
+        z_coord = cube.coord(axis="z")
+        layer_depth = get_coord_weights(cube, z_coord, broadcast=True)
+        cube.data = cube.core_data() * layer_depth
+        cube.units *= z_coord.units
+
+        return cube

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

@@ -5,8 +5,9 @@
 import pytest
 from cf_units import Unit
 from iris import NameConstraint
-from iris.coords import DimCoord
+from iris.coords import AuxCoord, DimCoord
 from iris.cube import Cube, CubeList
+from iris.util import new_axis
 
 import esmvalcore.cmor._fixes.icon.icon_xpp
 from esmvalcore.cmor._fixes.fix import GenericFix
@@ -18,6 +19,7 @@
     Gpp,
     Hfls,
     Hfss,
+    Msftmz,
     Rlut,
     Rlutcs,
     Rsutcs,
@@ -704,6 +706,64 @@ def test_hfss_fix(cubes_regular_grid):
     np.testing.assert_allclose(fixed_cube.data, [[[0.0, -1.0], [-2.0, -3.0]]])
 
 
+# Test msftmz (for extra fix)
+
+
+def test_get_msftmz_fix():
+    """Test getting of fix."""
+    fix = Fix.get_fixes("ICON", "ICON-XPP", "Omon", "msftmz")
+    assert fix == [Msftmz(None), AllVars(None), GenericFix(None)]
+
+
+def test_msftmz_fix(cubes_regular_grid):
+    """Test fix."""
+    depth_coord = AuxCoord(
+        10.0,
+        standard_name="depth",
+        long_name="depth below sea",
+        units="m",
+        attributes={"positive": "down"},
+    )
+    cube = cubes_regular_grid[0][..., [0]]
+    cube.coord("latitude").var_name = "lat"
+    cube.add_aux_coord(depth_coord, ())
+    cube = new_axis(cube, "depth")
+    cube.transpose([1, 0, 2, 3])
+    cubes = CubeList([cube.copy() * 0.0, cube.copy() * 1.0, cube.copy() * 2.0])
+    cubes[0].var_name = "atlantic_moc"
+    cubes[0].units = "kg s-1"
+    cubes[1].var_name = "pacific_moc"
+    cubes[1].units = "kg s-1"
+    cubes[2].var_name = "global_moc"
+    cubes[2].units = "kg s-1"
+
+    fixed_cubes = fix_metadata(cubes, "Omon", "msftmz")
+
+    assert len(fixed_cubes) == 1
+    cube = fixed_cubes[0]
+    assert cube.var_name == "msftmz"
+    assert (
+        cube.standard_name
+        == "ocean_meridional_overturning_mass_streamfunction"
+    )
+    assert cube.long_name == "Ocean Meridional Overturning Mass Streamfunction"
+
+    assert cube.units == "kg s-1"
+    assert "positive" not in cube.attributes
+    assert "invalid_units" not in cube.attributes
+
+    np.testing.assert_equal(
+        cube.coord("region").points,
+        ["atlantic_arctic_ocean", "indian_pacific_ocean", "global_ocean"],
+    )
+
+    assert cube.shape == (1, 3, 1, 2)
+    np.testing.assert_allclose(
+        cube.data,
+        [[[[0.0, 0.0]], [[0.0, 2.0]], [[0.0, 4.0]]]],
+    )
+
+
 # Test rlut (for extra fix)
 
 

tests/unit/preprocessor/_derive/test_amoc.py

@@ -48,11 +48,16 @@ def cubes():
 def test_amoc_preamble(cubes):
     derived_var = amoc.DerivedVariable()
 
-    cmip5_required = derived_var.required("CMIP5")
-    assert cmip5_required[0]["short_name"] == "msftmyz"
-    cmip6_required = derived_var.required("CMIP6")
-    assert cmip6_required[0]["short_name"] == "msftmz"
-    assert cmip6_required[1]["short_name"] == "msftyz"
+    assert derived_var.required("CMIP5") == [
+        {"short_name": "msftmyz", "mip": "Omon"},
+    ]
+    assert derived_var.required("CMIP6") == [
+        {"short_name": "msftmz", "optional": True},
+        {"short_name": "msftyz", "optional": True},
+    ]
+    assert derived_var.required("ICON") == [
+        {"short_name": "msftmz", "mip": "Omon"},
+    ]
 
     # if project s neither CMIP5 nor CMIP6
     with pytest.raises(ValueError, match="Project CMIPX can not be used"):

tests/unit/preprocessor/_derive/test_phcint.py

@@ -0,0 +1,73 @@
+"""Test derivation of ``phcint``."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+import numpy as np
+import pytest
+from iris.coords import DimCoord
+from iris.cube import CubeList
+
+from esmvalcore.preprocessor._derive import derive, get_required
+
+if TYPE_CHECKING:
+    from iris.cube import Cube
+
+
+@pytest.fixture
+def cubes(realistic_4d_cube: Cube) -> CubeList:
+    depth_coord = DimCoord(
+        [500.0],
+        bounds=[[0.0, 1000.0]],
+        standard_name="depth",
+        units="m",
+        attributes={"positive": "down"},
+    )
+    realistic_4d_cube.remove_coord("air_pressure")
+    realistic_4d_cube.add_dim_coord(depth_coord, 1)
+    realistic_4d_cube.var_name = "thetao"
+    return CubeList([realistic_4d_cube])
+
+
+@pytest.mark.parametrize("project", ["CMIP3", "CMIP5", "CMIP6", "CMIP7"])
+def test_get_required(project: str) -> None:
+    assert get_required("phcint", project) == [{"short_name": "thetao"}]
+
+
+def test_derive(cubes: CubeList) -> None:
+    short_name = "phcint"
+    long_name = "Integrated Ocean Heat Content from Potential Temperature"
+    units = "J m-2"
+    standard_name = "integral_wrt_depth_of_sea_water_potential_temperature_expressed_as_heat_content"
+
+    derived_cube = derive(
+        cubes,
+        short_name=short_name,
+        long_name=long_name,
+        units=units,
+        standard_name=standard_name,
+    )
+
+    assert derived_cube.standard_name == standard_name
+    assert derived_cube.long_name == long_name
+    assert derived_cube.var_name == short_name
+    assert derived_cube.units == units
+    assert derived_cube.shape == (2, 1, 2, 3)
+    expected_data = np.ma.masked_invalid(
+        [
+            [
+                [
+                    [0.0, np.nan, np.nan],
+                    [np.nan, 16366760000.0, 20458450000.0],
+                ],
+            ],
+            [
+                [
+                    [24550140000.0, 28641830000.0, 32733520000.0],
+                    [36825210000.0, 40916900000.0, 45008590000.0],
+                ],
+            ],
+        ],
+    )
+    np.testing.assert_allclose(derived_cube.data, expected_data)