diff --git a/esmvalcore/cmor/_fixes/icon/icon_xpp.py b/esmvalcore/cmor/_fixes/icon/icon_xpp.py index d26e3a15af..baeb9bee8d 100644 --- a/esmvalcore/cmor/_fixes/icon/icon_xpp.py +++ b/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 diff --git a/esmvalcore/cmor/tables/cmip5-custom/CMOR_phcint.dat b/esmvalcore/cmor/tables/cmip5-custom/CMOR_phcint.dat new file mode 100644 index 0000000000..b28c63cb6a --- /dev/null +++ b/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 +!---------------------------------- +! diff --git a/esmvalcore/cmor/tables/cmip6-custom/CMIP6_custom.json b/esmvalcore/cmor/tables/cmip6-custom/CMIP6_custom.json index d399838a75..fa0b8bdd8c 100644 --- a/esmvalcore/cmor/tables/cmip6-custom/CMIP6_custom.json +++ b/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", diff --git a/esmvalcore/config/configurations/defaults/extra_facets_icon.yml b/esmvalcore/config/configurations/defaults/extra_facets_icon.yml index ab8418d5bb..bc346b096f 100644 --- a/esmvalcore/config/configurations/defaults/extra_facets_icon.yml +++ b/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} diff --git a/esmvalcore/preprocessor/_derive/amoc.py b/esmvalcore/preprocessor/_derive/amoc.py index 67b179f0dd..f6ccef14ef 100644 --- a/esmvalcore/preprocessor/_derive/amoc.py +++ b/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) diff --git a/esmvalcore/preprocessor/_derive/phcint.py b/esmvalcore/preprocessor/_derive/phcint.py new file mode 100644 index 0000000000..f19e7f3416 --- /dev/null +++ b/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 diff --git a/tests/integration/cmor/_fixes/icon/test_icon_xpp.py b/tests/integration/cmor/_fixes/icon/test_icon_xpp.py index aa50be5ebb..369bd161d7 100644 --- a/tests/integration/cmor/_fixes/icon/test_icon_xpp.py +++ b/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) diff --git a/tests/unit/preprocessor/_derive/test_amoc.py b/tests/unit/preprocessor/_derive/test_amoc.py index 3b2150d830..a6fb43b9a4 100644 --- a/tests/unit/preprocessor/_derive/test_amoc.py +++ b/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"): diff --git a/tests/unit/preprocessor/_derive/test_phcint.py b/tests/unit/preprocessor/_derive/test_phcint.py new file mode 100644 index 0000000000..75e118144a --- /dev/null +++ b/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)