#723: Add soil thickness processor (#758)

* #723: add a soil thickness calculation function.

* #723: Create functional test for slthick processor

Author: mo-laurenboon

cdds/cdds/common/constants.py

@@ -1,4 +1,4 @@
-# (C) British Crown Copyright 2019-2025, Met Office.
+# (C) British Crown Copyright 2019-2026, Met Office.
 # Please see LICENSE.md for license details.
 """The :mod:`constants` module contains constants (values that should never be changed by a user and exist for
 readability and maintainability purposes) for all CDDS components.
@@ -8,6 +8,7 @@
 
 # Note that orography (m01s00i033) must not go in this list
 ANCIL_VARIABLES = [
+    # soil moisture content
     # root depth
     'm01s00i009',
     'm01s00i216',

mip_convert/mip_convert/plugins/base/data/processors.py

@@ -1,4 +1,4 @@
-# (C) British Crown Copyright 2016-2025, Met Office.
+# (C) British Crown Copyright 2016-2026, Met Office.
 # Please see LICENSE.md for license details.
 
 """Module containing processor functions.  These processors can be referred
@@ -560,7 +560,7 @@ def fix_packing_division(numerator, denominator):
 def _z_axis(cube):
     z_coords = [coord for coord in cube.coords() if guess_coord_axis(coord) == 'Z']
     if len(z_coords) != 1:
-        raise ValueError("cube should have exactly one 'Z' axis")
+        raise ValueError(f"cube should have exactly one 'Z' axis but found {len(z_coords)}")
     return z_coords[0]
 
 
@@ -2321,19 +2321,77 @@ def tos_ORCA12(tos_con, tossq_con):
     return tos_con
 
 
-def calc_rootd(soil_cube, frac_cube, ice_class=None):
-    """
-    Returns a cube of the maximum rooting depth calculated from a cube on a model's standard latitude-longitude grid and
+def calc_slthick(soil_cube, frac_cube, ice_class=None):
+    """Returns a cube of soil layer thickness calculated from a cube on a model's standard latitude-longitude grid and
     on soil levels.
 
     Parameters
     ----------
     soil_cube: :class:`iris.cube.Cube`
         A cube containing data on soil levels for that model.
-
     frac_cube: :class:`iris.cube.Cube`
         A cube containing JULES tile fractions for that model.
+    ice_class: str
+        The tile ID string for land ice. If present, the soil layer thickness is set to zero on land ice points.
+
+    Returns
+    -------
+    :class:`iris.cube.Cube`
+        A cube containing the thickness of each soil level of that model.
 
+    Raises
+    ------
+    RuntimeError
+        If the given soil cube does not contain cell bounds on the z axis or the soil cube is not masked.
+    ValueError
+        If the soil is identified as having a negative thickness.
+    """
+    # Calculate the thickness of a soil layer from the cell bounds of the vertical coord of the input cube.
+    depth_coord = _z_axis(soil_cube)
+    slthick_data = soil_cube.data.copy()
+    for cell, data in zip(depth_coord, slthick_data):
+        if not cell.has_bounds():
+            raise RuntimeError("The provided cube does not contian cell bounds on the identified Z axis.")
+        data[:] = cell.bounds[0][1] - cell.bounds[0][0]
+    if np.any(data < 0):
+        raise ValueError("Soil cannot have a negative thickness.")
+
+    # Set max root depth to 0.0 m on land ice points.
+    if ice_class:
+        ice_frac = frac_cube.extract(_pseudo_constraint(ice_class))
+        ice_cube = _collapse_pseudo(ice_frac, SUM)
+        ice_mask = ice_cube.data.data > 1e-6
+        slthick_data[:, ice_mask] = 0.0
+
+    # Copy the land/sea mask from the source cube to the new array.
+    if not is_masked(soil_cube.data):
+        raise RuntimeError("The soil cube is not masked, any land/sea masks will not be copied to the new array.")
+
+    slthick_data.mask = soil_cube.data.mask
+
+    # Create a Cube of the soil level thickness data.
+    dim_coords_and_dims = [(coord, k) for (k, coord) in enumerate(soil_cube.dim_coords)]
+
+    slthick_cube = iris.cube.Cube(slthick_data,
+                                  standard_name="cell_thickness",
+                                  long_name="Thickness of Soil Layers",
+                                  units=depth_coord.units,
+                                  dim_coords_and_dims=dim_coords_and_dims,
+                                  )
+
+    return slthick_cube
+
+
+def calc_rootd(soil_cube, frac_cube, ice_class=None):
+    """Returns a cube of the maximum rooting depth calculated from a cube on a model's standard latitude-longitude grid
+    and on soil levels.
+
+    Parameters
+    ----------
+    soil_cube: :class:`iris.cube.Cube`
+        A cube containing data on soil levels for that model.
+    frac_cube: :class:`iris.cube.Cube`
+        A cube containing JULES tile fractions for that model.
     ice_class: str
         Tile ID string for land ice.  If present, max root depth is set to zero on land ice points.
 

mip_convert/mip_convert/plugins/hadgem3_gc5/data/HadGEM3GC5_mappings.cfg

@@ -1,4 +1,4 @@
-# (C) British Crown Copyright 2018-2025, Met Office.
+# (C) British Crown Copyright 2018-2026, Met Office.
 # Please see LICENSE.md for license details.
 #
 # This 'model to MIP mappings' configuration file contains sections
@@ -42,8 +42,14 @@ expression = m01s03i210[lbproc=128]
 mip_table_id = atmos
 units = m s-1
 
+[slthick_ti-sl-hxy-lnd]
+dimension = longitude latitude sdepth
+expression = calc_slthick(m01s00i009, m01s00i216, ice_class='ice')
+mip_table_id = land
+units = m
+
 [rootd_ti-u-hxy-lnd]
 dimension = longitude latitude
 expression = calc_rootd(m01s00i009, m01s00i216, ice_class='ice')
 mip_table_id = land
-units = m
+units = m

mip_convert/mip_convert/tests/test_functional/test_functional_cmip7/test_cmip7_fx_slthick_ti_sl_hxy_lnd.py

@@ -0,0 +1,54 @@
+# (C) British Crown Copyright 2022-2026, Met Office.
+# Please see LICENSE.md for license details.
+import os
+
+import pytest
+
+from mip_convert.tests.test_functional.test_command_line import AbstractFunctionalTests
+from mip_convert.tests.test_functional.utils.configurations import Cmip7TestData, SpecificInfo
+from mip_convert.tests.test_functional.utils.directories import (get_cmor_log, get_output_dir,
+                                                                 MODEL_OUTPUT_DIR,
+                                                                 ROOT_OUTPUT_CASES_DIR)
+
+
+class TestCmip7_fx_slthick_ti_sl_hxy_lnd(AbstractFunctionalTests):
+
+    def get_test_data(self):
+        test_location = os.path.join(ROOT_OUTPUT_CASES_DIR, 'test_CMIP7_land_slthick_ti-sl-hxy-lnd')
+        return Cmip7TestData(
+            mip_table='land',
+            variables=['slthick_ti-sl-hxy-lnd'],
+            specific_info=SpecificInfo(
+                common={
+                    'test_location': test_location
+                },
+                cmor_setup={
+                    'cmor_log_file': get_cmor_log(test_location)
+                },
+                cmor_dataset={
+                    'output_dir': get_output_dir(test_location)
+                },
+                request={
+                    'ancil_files': ' '.join([
+                        os.path.join(MODEL_OUTPUT_DIR, 'u-dk469', 'ancil', 'test_input_land_fx_u-dk469.pp'),
+                    ]),
+                    'model_output_dir': MODEL_OUTPUT_DIR,
+                    'run_bounds': '2257-01-01T00:00:00 2257-02-01T00:00:00',
+                    'suite_id': 'u-dk469',
+                    'mip_convert_plugin': 'HadGEM3GC5'
+                },
+                streams={
+                    'ancil': {'CMIP7_land@fx': 'slthick_ti-sl-hxy-lnd'}
+                },
+                other={
+                    'reference_version': 'v1',
+                    'filenames': ['slthick_ti-sl-hxy-lnd_fx_glb_gn_PCMDI-test-1-0_1pctCO2_r1i1dp1f1.nc'],
+                    'ignore_history': True,
+                }
+            )
+        )
+
+    @pytest.mark.slow
+    def test_cmip7_efx_slthick(self):
+        self.maxDiff = True
+        self.check_convert(plugin_id="CMIP7")

mip_convert/mip_convert/tests/test_functional/utils/directories.py

@@ -1,4 +1,4 @@
-# (C) British Crown Copyright 2022-2025, Met Office.
+# (C) British Crown Copyright 2022-2026, Met Office.
 # Please see LICENSE.md for license details.
 import os
 
@@ -17,7 +17,7 @@
 CORDEX_MIP_TABLE_DIR = '{}/CORDEX/cordex-cmip6-cmor-tables/Tables'.format(ROOT_MIP_TABLES_DIR)
 ARISE_MIP_TABLE_DIR = '{}/ARISE/for_functional_tests'.format(ROOT_MIP_TABLES_DIR)
 CMIP6_MIP_TABLE_DIR = '{}/CMIP6/for_functional_tests'.format(ROOT_MIP_TABLES_DIR)
-CMIP7_MIP_TABLE_DIR = '{}/CMIP7/DR-1.2.2.2-v1.0.1'.format(ROOT_MIP_TABLES_DIR)
+CMIP7_MIP_TABLE_DIR = '{}/CMIP7/for_functional_tests'.format(ROOT_MIP_TABLES_DIR)
 SEASONAL_MIP_TABLE_DIR = '{}/SEASONAL/for_functional_tests'.format(ROOT_MIP_TABLES_DIR)
 NAHOSMIP_MIP_TABLE_DIR = '{}/GCModelDev/for_functional_tests'.format(ROOT_MIP_TABLES_DIR)
 

mip_convert/mip_convert/tests/test_process/test_processors.py

@@ -1,4 +1,4 @@
-# (C) British Crown Copyright 2016-2025, Met Office.
+# (C) British Crown Copyright 2016-2026, Met Office.
 # Please see LICENSE.md for license details.
 # pylint: disable = missing-docstring, invalid-name, too-many-public-methods
 # pylint: disable = no-member, no-value-for-parameter
@@ -16,12 +16,9 @@
 import numpy as np
 
 from mip_convert.plugins.base.data.processors import (
-    area_mean, areacella, calc_rho_mean, calc_zostoga,
-    combine_cubes_to_basin_coord, eos_insitu, fix_clmisr_height,
-    land_class_area, land_class_mean, level_sum, mask_copy,
-    mask_zeros, mask_polar_column_zonal_means,
-    mon_mean_from_day,
-    ocean_quasi_barotropic_streamfunc, tile_ids_for_class, volcello, vortmean,
+    area_mean, areacella, calc_rho_mean, calc_zostoga, combine_cubes_to_basin_coord, eos_insitu, fix_clmisr_height,
+    land_class_area, land_class_mean, level_sum, mask_copy, mask_zeros, mask_polar_column_zonal_means,
+    mon_mean_from_day, ocean_quasi_barotropic_streamfunc, tile_ids_for_class, volcello, vortmean,
     annual_from_monthly_2d, annual_from_monthly_3d, calculate_thkcello_weights, check_data_is_monthly)
 from mip_convert.tests.common import dummy_cube
 from functools import reduce