Repository: ESMValGroup/ESMValCore

Issue:
Preprocessor shape 'extract_shape' breaks if cell measure is time dependent
**Describe the bug**
Running extract_shape with `decompose: True` on data which has a time dependent cell measure (like `volcello`) attached to it fails with the following error:
```
iris.exceptions.CannotAddError: Unequal lengths. Cube dimension 1 => 3; metadata 'ocean_volume' dimension 0 => 2.
```
This happens because a new dimension is added between `time` and the spatial coordinates `depth lat lon`, which leads to a cell_measure with unmatched dimensions with respect to the cube because it's missing the region coordinate.

**Please attach**
  - The recipe that you are trying to run, you can find a copy in the `run` directory in the output directory
  - The `main_log_debug.txt` file, this can also be found in the `run` directory in the output directory


Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2348?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 94.11%. Comparing base [(`bb7866e`)](https://app.codecov.io/gh/ESMValGroup/ESMValCore/commit/bb7866e6a136f6e011bc6419c4e3ed9eec25721e?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) to head [(`faed177`)](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2348?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup).


<details><summary>Additional details and impacted files</summary>


```diff
@@           Coverage Diff           @@
##             main    #2348   +/-   ##
=======================================
  Coverage   94.11%   94.11%           
=======================================
  Files         241      241           
  Lines       13376    13381    +5     
=======================================
+ Hits        12589    12594    +5     
  Misses        787      787           
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2348?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup).

sloosvel: Thanks for the comments @schlunma , I added a comment and changed a bit the condition to also make it work when only one region is extracted and the result cube has no `shape_id` dimension. Here is a recipe to test both cases:

```yaml
  description: |
    This a recipe to monitor experiments by automatizing the creation of plots
    from raw preprocessor output.
  authors:
    - vegas-regidor_javier
  maintainer:
    - loosveldt-tomas_saskia


datasets:
- dataset: EC-Earth3-CC
  ensemble: r1i1p1f1
  grid: gn
  exp: historical
  mip: Omon
  project: CMIP6

preprocessors:
  extract_two_regions:
    extract_shape:
      shapefile: AR6
      method: contains
      crop: False
      decomposed: True
      ids: {Acronym: [EAO, WAF]}
    volume_statistics:
      operator: mean
  extract_one_region:
    extract_shape:
      shapefile: AR6
      method: contains
      crop: False
      decomposed: True
      ids: {Acronym: [EAO]}
    volume_statistics:
      operator: mean

diagnostics:
  timeseries:
    variables:
      thetao_1:
        short_name: thetao
        timerange: '1950/1952'
        preprocessor: extract_two_regions
      thetao_2:
        short_name: thetao
        timerange: '1950/1952'
        preprocessor: extract_one_region
    scripts: null
```


Patch:
diff --git a/esmvalcore/cmor/tables/custom/CMOR_clhmtisccp.dat b/esmvalcore/cmor/tables/custom/CMOR_clhmtisccp.dat
index 504b69164..3c8b049ec 100644
--- a/esmvalcore/cmor/tables/custom/CMOR_clhmtisccp.dat
+++ b/esmvalcore/cmor/tables/custom/CMOR_clhmtisccp.dat
@@ -17,6 +17,5 @@ comment:           at the top of the atmosphere (to be compared with satellite m
 !----------------------------------
 dimensions:        longitude latitude time
 type:              real
-positive:          up
 !----------------------------------
 !
diff --git a/esmvalcore/cmor/tables/custom/CMOR_clhtkisccp.dat b/esmvalcore/cmor/tables/custom/CMOR_clhtkisccp.dat
index 6ecb43a91..d34d39782 100644
--- a/esmvalcore/cmor/tables/custom/CMOR_clhtkisccp.dat
+++ b/esmvalcore/cmor/tables/custom/CMOR_clhtkisccp.dat
@@ -17,6 +17,5 @@ comment:           at the top of the atmosphere (to be compared with satellite m
 !----------------------------------
 dimensions:        longitude latitude time
 type:              real
-positive:          up
 !----------------------------------
 !
diff --git a/esmvalcore/cmor/tables/custom/CMOR_cllmtisccp.dat b/esmvalcore/cmor/tables/custom/CMOR_cllmtisccp.dat
index a16279a42..86438496b 100644
--- a/esmvalcore/cmor/tables/custom/CMOR_cllmtisccp.dat
+++ b/esmvalcore/cmor/tables/custom/CMOR_cllmtisccp.dat
@@ -17,6 +17,5 @@ comment:           at the top of the atmosphere (to be compared with satellite m
 !----------------------------------
 dimensions:        longitude latitude time
 type:              real
-positive:          up
 !----------------------------------
 !
diff --git a/esmvalcore/cmor/tables/custom/CMOR_clltkisccp.dat b/esmvalcore/cmor/tables/custom/CMOR_clltkisccp.dat
index c15c19d12..a1b989f98 100644
--- a/esmvalcore/cmor/tables/custom/CMOR_clltkisccp.dat
+++ b/esmvalcore/cmor/tables/custom/CMOR_clltkisccp.dat
@@ -17,6 +17,5 @@ comment:           at the top of the atmosphere (to be compared with satellite m
 !----------------------------------
 dimensions:        longitude latitude time
 type:              real
-positive:          up
 !----------------------------------
 !
diff --git a/esmvalcore/cmor/tables/custom/CMOR_clmmtisccp.dat b/esmvalcore/cmor/tables/custom/CMOR_clmmtisccp.dat
index 707709843..a7482c24b 100644
--- a/esmvalcore/cmor/tables/custom/CMOR_clmmtisccp.dat
+++ b/esmvalcore/cmor/tables/custom/CMOR_clmmtisccp.dat
@@ -17,6 +17,5 @@ comment:           at the top of the atmosphere (to be compared with satellite m
 !----------------------------------
 dimensions:        longitude latitude time
 type:              real
-positive:          up
 !----------------------------------
 !
diff --git a/esmvalcore/cmor/tables/custom/CMOR_clmtkisccp.dat b/esmvalcore/cmor/tables/custom/CMOR_clmtkisccp.dat
index 7ec19cf1a..1bdfe4133 100644
--- a/esmvalcore/cmor/tables/custom/CMOR_clmtkisccp.dat
+++ b/esmvalcore/cmor/tables/custom/CMOR_clmtkisccp.dat
@@ -17,6 +17,5 @@ comment:           at the top of the atmosphere (to be compared with satellite m
 !----------------------------------
 dimensions:        longitude latitude time
 type:              real
-positive:          up
 !----------------------------------
 !
diff --git a/esmvalcore/cmor/tables/custom/CMOR_lwcre.dat b/esmvalcore/cmor/tables/custom/CMOR_lwcre.dat
index 67a613e5b..6818af1b9 100644
--- a/esmvalcore/cmor/tables/custom/CMOR_lwcre.dat
+++ b/esmvalcore/cmor/tables/custom/CMOR_lwcre.dat
@@ -17,6 +17,6 @@ comment:           at the top of the atmosphere (to be compared with satellite m
 !----------------------------------
 dimensions:        longitude latitude time
 type:              real
-positive:          up
+positive:          down
 !----------------------------------
 !
diff --git a/esmvalcore/cmor/tables/custom/CMOR_netcre.dat b/esmvalcore/cmor/tables/custom/CMOR_netcre.dat
index 689afdf8a..6bce0c491 100644
--- a/esmvalcore/cmor/tables/custom/CMOR_netcre.dat
+++ b/esmvalcore/cmor/tables/custom/CMOR_netcre.dat
@@ -17,6 +17,6 @@ comment:           at the top of the atmosphere (to be compared with satellite m
 !----------------------------------
 dimensions:        longitude latitude time
 type:              real
-positive:          up
+positive:          down
 !----------------------------------
 !
diff --git a/esmvalcore/cmor/tables/custom/CMOR_rsnt.dat b/esmvalcore/cmor/tables/custom/CMOR_rsnt.dat
index f3c84d1d7..38b43c7f4 100644
--- a/esmvalcore/cmor/tables/custom/CMOR_rsnt.dat
+++ b/esmvalcore/cmor/tables/custom/CMOR_rsnt.dat
@@ -17,6 +17,6 @@ comment:           at the top of the atmosphere (to be compared with satellite m
 !----------------------------------
 dimensions:        longitude latitude time
 type:              real
-positive:          up
+positive:          down
 !----------------------------------
 !
diff --git a/esmvalcore/cmor/tables/custom/CMOR_swcre.dat b/esmvalcore/cmor/tables/custom/CMOR_swcre.dat
index 38caf0a63..9cfe7df0e 100644
--- a/esmvalcore/cmor/tables/custom/CMOR_swcre.dat
+++ b/esmvalcore/cmor/tables/custom/CMOR_swcre.dat
@@ -17,6 +17,6 @@ comment:           at the top of the atmosphere (to be compared with satellite m
 !----------------------------------
 dimensions:        longitude latitude time
 type:              real
-positive:          up
+positive:          down
 !----------------------------------
 !
diff --git a/esmvalcore/preprocessor/_area.py b/esmvalcore/preprocessor/_area.py
index 7d7bbe0c9..61f0402b3 100644
--- a/esmvalcore/preprocessor/_area.py
+++ b/esmvalcore/preprocessor/_area.py
@@ -938,6 +938,23 @@ def _mask_cube(cube: Cube, masks: dict[str, np.ndarray]) -> Cube:
     result = fix_coordinate_ordering(cubelist.merge_cube())
     if cube.cell_measures():
         for measure in cube.cell_measures():
+            # Cell measures that are time-dependent, with 4 dimension and
+            # an original shape of (time, depth, lat, lon), need to be
+            # broadcasted to the cube with 5 dimensions and shape
+            # (time, shape_id, depth, lat, lon)
+            if measure.ndim > 3 and result.ndim > 4:
+                data = measure.core_data()
+                data = da.expand_dims(data, axis=(1,))
+                data = da.broadcast_to(data, result.shape)
+                measure = iris.coords.CellMeasure(
+                    data,
+                    standard_name=measure.standard_name,
+                    long_name=measure.long_name,
+                    units=measure.units,
+                    measure=measure.measure,
+                    var_name=measure.var_name,
+                    attributes=measure.attributes,
+                )
             add_cell_measure(result, measure, measure.measure)
     if cube.ancillary_variables():
         for ancillary_variable in cube.ancillary_variables():
diff --git a/esmvalcore/preprocessor/_derive/lwcre.py b/esmvalcore/preprocessor/_derive/lwcre.py
index ed7b1ada7..718956028 100644
--- a/esmvalcore/preprocessor/_derive/lwcre.py
+++ b/esmvalcore/preprocessor/_derive/lwcre.py
@@ -31,5 +31,6 @@ class DerivedVariable(DerivedVariableBase):
 
         lwcre_cube = rlutcs_cube - rlut_cube
         lwcre_cube.units = rlut_cube.units
+        lwcre_cube.attributes['positive'] = 'down'
 
         return lwcre_cube
diff --git a/esmvalcore/preprocessor/_derive/netcre.py b/esmvalcore/preprocessor/_derive/netcre.py
index 806a79000..9c65e0ec7 100644
--- a/esmvalcore/preprocessor/_derive/netcre.py
+++ b/esmvalcore/preprocessor/_derive/netcre.py
@@ -43,5 +43,6 @@ class DerivedVariable(DerivedVariableBase):
 
         netcre_cube = lwcre_cube + swcre_cube
         netcre_cube.units = lwcre_cube.units
+        netcre_cube.attributes['positive'] = 'down'
 
         return netcre_cube
diff --git a/esmvalcore/preprocessor/_derive/rsnt.py b/esmvalcore/preprocessor/_derive/rsnt.py
index 7838c3a0b..dbd04ca2c 100644
--- a/esmvalcore/preprocessor/_derive/rsnt.py
+++ b/esmvalcore/preprocessor/_derive/rsnt.py
@@ -30,5 +30,7 @@ class DerivedVariable(DerivedVariableBase):
             Constraint(name='toa_outgoing_shortwave_flux'))
 
         rsnt_cube = rsdt_cube - rsut_cube
+        rsnt_cube.units = rsdt_cube.units
+        rsnt_cube.attributes['positive'] = 'down'
 
         return rsnt_cube
diff --git a/esmvalcore/preprocessor/_derive/swcre.py b/esmvalcore/preprocessor/_derive/swcre.py
index 9fcf056f4..c1873f622 100644
--- a/esmvalcore/preprocessor/_derive/swcre.py
+++ b/esmvalcore/preprocessor/_derive/swcre.py
@@ -30,5 +30,7 @@ class DerivedVariable(DerivedVariableBase):
             Constraint(name='toa_outgoing_shortwave_flux_assuming_clear_sky'))
 
         swcre_cube = rsutcs_cube - rsut_cube
+        swcre_cube.units = rsut_cube.units
+        swcre_cube.attributes['positive'] = 'down'
 
         return swcre_cube


Test patch:
diff --git a/tests/unit/preprocessor/_area/test_area.py b/tests/unit/preprocessor/_area/test_area.py
index 8835014ff..f7f507e49 100644
--- a/tests/unit/preprocessor/_area/test_area.py
+++ b/tests/unit/preprocessor/_area/test_area.py
@@ -1399,5 +1399,51 @@ def test_meridional_statistics_invalid_norm_fail(make_testcube):
         meridional_statistics(make_testcube, 'sum', normalize='x')
 
 
+def test_time_dependent_volcello():
+    coord_sys = iris.coord_systems.GeogCS(iris.fileformats.pp.EARTH_RADIUS)
+    data = np.ma.ones((2, 3, 2, 2))
+
+    time = iris.coords.DimCoord([15, 45],
+                                standard_name='time',
+                                bounds=[[1., 30.], [30., 60.]],
+                                units=Unit('days since 1950-01-01',
+                                calendar='gregorian'))
+
+    zcoord = iris.coords.DimCoord([0.5, 5., 50.],
+                                  long_name='zcoord',
+                                  bounds=[[0., 2.5], [2.5, 25.],
+                                          [25., 250.]],
+                                  units='m',
+                                  attributes={'positive': 'down'})
+    lons = iris.coords.DimCoord([1.5, 2.5],
+                                standard_name='longitude',
+                                bounds=[[1., 2.], [2., 3.]],
+                                units='degrees_east',
+                                coord_system=coord_sys)
+    lats = iris.coords.DimCoord([1.5, 2.5],
+                                standard_name='latitude',
+                                bounds=[[1., 2.], [2., 3.]],
+                                units='degrees_north',
+                                coord_system=coord_sys)
+    coords_spec4 = [(time, 0), (zcoord, 1), (lats, 2), (lons, 3)]
+    cube = iris.cube.Cube(data, dim_coords_and_dims=coords_spec4)
+    volcello = iris.coords.CellMeasure(
+        data,
+        standard_name='ocean_volume',
+        units='m3',
+        measure='volume')
+    cube.add_cell_measure(volcello, range(0, volcello.ndim))
+    cube = extract_shape(
+        cube,
+        'AR6',
+        method='contains',
+        crop=False,
+        decomposed=True,
+        ids={'Acronym': ['EAO', 'WAF']},
+    )
+
+    assert cube.shape == cube.cell_measure('ocean_volume').shape
+
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/unit/preprocessor/_derive/test_lwcre.py b/tests/unit/preprocessor/_derive/test_lwcre.py
new file mode 100644
index 000000000..c70807a14
--- /dev/null
+++ b/tests/unit/preprocessor/_derive/test_lwcre.py
@@ -0,0 +1,28 @@
+"""Test derivation of `lwcre`."""
+import numpy as np
+import pytest
+from iris.cube import Cube, CubeList
+
+import esmvalcore.preprocessor._derive.lwcre as lwcre
+
+
+@pytest.fixture
+def cubes():
+    rlut_cube = Cube(
+        3, standard_name='toa_outgoing_longwave_flux', units='W m-2'
+    )
+    rlutcs_cube = Cube(
+        1,
+        standard_name='toa_outgoing_longwave_flux_assuming_clear_sky',
+        units='W m-2',
+    )
+    return CubeList([rlut_cube, rlutcs_cube])
+
+
+def test_lwcre_calculation(cubes):
+    """Test calculation of `lwcre`."""
+    derived_var = lwcre.DerivedVariable()
+    out_cube = derived_var.calculate(cubes)
+    np.testing.assert_equal(out_cube.data, -2)
+    assert out_cube.units == 'W m-2'
+    assert out_cube.attributes['positive'] == 'down'
diff --git a/tests/unit/preprocessor/_derive/test_netcre.py b/tests/unit/preprocessor/_derive/test_netcre.py
new file mode 100644
index 000000000..f520b1bc1
--- /dev/null
+++ b/tests/unit/preprocessor/_derive/test_netcre.py
@@ -0,0 +1,36 @@
+"""Test derivation of `netcre`."""
+import numpy as np
+import pytest
+from iris.cube import Cube, CubeList
+
+import esmvalcore.preprocessor._derive.netcre as netcre
+
+
+@pytest.fixture
+def cubes():
+    rlut_cube = Cube(
+        3, standard_name='toa_outgoing_longwave_flux', units='W m-2'
+    )
+    rlutcs_cube = Cube(
+        1,
+        standard_name='toa_outgoing_longwave_flux_assuming_clear_sky',
+        units='W m-2',
+    )
+    rsut_cube = Cube(
+        3, standard_name='toa_outgoing_shortwave_flux', units='W m-2'
+    )
+    rsutcs_cube = Cube(
+        1,
+        standard_name='toa_outgoing_shortwave_flux_assuming_clear_sky',
+        units='W m-2',
+    )
+    return CubeList([rlut_cube, rlutcs_cube, rsut_cube, rsutcs_cube])
+
+
+def test_netcre_calculation(cubes):
+    """Test calculation of `netcre`."""
+    derived_var = netcre.DerivedVariable()
+    out_cube = derived_var.calculate(cubes)
+    np.testing.assert_equal(out_cube.data, -4)
+    assert out_cube.units == 'W m-2'
+    assert out_cube.attributes['positive'] == 'down'
diff --git a/tests/unit/preprocessor/_derive/test_rsnt.py b/tests/unit/preprocessor/_derive/test_rsnt.py
new file mode 100644
index 000000000..4cc16f170
--- /dev/null
+++ b/tests/unit/preprocessor/_derive/test_rsnt.py
@@ -0,0 +1,26 @@
+"""Test derivation of `rsnt`."""
+import numpy as np
+import pytest
+from iris.cube import Cube, CubeList
+
+import esmvalcore.preprocessor._derive.rsnt as rsnt
+
+
+@pytest.fixture
+def cubes():
+    rsdt_cube = Cube(
+        3, standard_name='toa_incoming_shortwave_flux', units='W m-2'
+    )
+    rsut_cube = Cube(
+        1, standard_name='toa_outgoing_shortwave_flux', units='W m-2'
+    )
+    return CubeList([rsdt_cube, rsut_cube])
+
+
+def test_rsnt_calculation(cubes):
+    """Test calculation of `rsnt`."""
+    derived_var = rsnt.DerivedVariable()
+    out_cube = derived_var.calculate(cubes)
+    np.testing.assert_equal(out_cube.data, 2)
+    assert out_cube.units == 'W m-2'
+    assert out_cube.attributes['positive'] == 'down'
diff --git a/tests/unit/preprocessor/_derive/test_swcre.py b/tests/unit/preprocessor/_derive/test_swcre.py
new file mode 100644
index 000000000..165477061
--- /dev/null
+++ b/tests/unit/preprocessor/_derive/test_swcre.py
@@ -0,0 +1,28 @@
+"""Test derivation of `swcre`."""
+import numpy as np
+import pytest
+from iris.cube import Cube, CubeList
+
+import esmvalcore.preprocessor._derive.swcre as swcre
+
+
+@pytest.fixture
+def cubes():
+    rsut_cube = Cube(
+        3, standard_name='toa_outgoing_shortwave_flux', units='W m-2'
+    )
+    rsutcs_cube = Cube(
+        1,
+        standard_name='toa_outgoing_shortwave_flux_assuming_clear_sky',
+        units='W m-2',
+    )
+    return CubeList([rsut_cube, rsutcs_cube])
+
+
+def test_swcre_calculation(cubes):
+    """Test calculation of `swcre`."""
+    derived_var = swcre.DerivedVariable()
+    out_cube = derived_var.calculate(cubes)
+    np.testing.assert_equal(out_cube.data, -2)
+    assert out_cube.units == 'W m-2'
+    assert out_cube.attributes['positive'] == 'down'