Repository: ESMValGroup/ESMValCore

Issue:
concatenate experiments handles time ordering incorrectly
I am analysing the 1pctCO2 and esm-1pct-brch-1000PgC experiments, where the latter branches from the former after about 70 years. The simulation time for 1pctCO2 is longer than 70 years, so there will be several years that are present in both experiments. When concatenating both experiments I need to take the first 70 years from 1pctCO2 and the remainder from esm-1pct-brch-1000PgC, essentially discarding the years in 1pctCO2 after year 70. 

Thats the default behaviour of ESMValTool and iris, except the above logic operates on a per-file basis rather than on per-experiment. So if the experiment datasets consist of multiple files (which they frequently do), we get nonsensical output that depends on how the files were chunked and saved.

Below is an MRE:

```python
import esmvalcore.dataset
import iris
import iris.quickplot as qplt

ds = esmvalcore.dataset.Dataset( **{
    "dataset": "CESM2",
    "project": "CMIP6", 
    "exp": ["1pctCO2", "esm-1pct-brch-1000PgC"], 
    "ensemble": "r1i1p1f1", 
    "grid": "gn", 
    "timerange": "0001/0217",
    "short_name":"tas",
    "mip":"Amon",
})

qplt.plot(ds.load().collapsed(['latitude', 'longitude'], iris.analysis.MEAN))
```

which gives the following plot:
<img width="681" alt="image" src="https://github.com/ESMValGroup/ESMValTool/assets/3708393/140d57c3-a244-49f9-8313-b9e757bce5b0">

One can see the discontinuity in the temperature value that is caused by taking some timeranges from the wrong experiment.

There is a single unambiguous way to concatenate 2 experiments with overlapping timestamps in the CMIP context and that is at the branching time. I have a local bugfix that first concatenates the experiments-wise and then concatenates the result of that, which I have pushed to https://github.com/ESMValGroup/ESMValCore/tree/experiments_first (will open a PR now)

There are several issues pertaining to concatenation, e.g. this one #1423 which considerable activity over the years. However the underlying issue persists.

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2343?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) Report
Attention: Patch coverage is `75.00000%` with `1 lines` in your changes are missing coverage. Please review.
> 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 [(`0fe5c76`)](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2343?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup).
> Report is 2 commits behind head on main.

| [Files](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2343?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) | Patch % | Lines |
|---|---|---|
| [esmvalcore/preprocessor/\_io.py](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2343?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup#diff-ZXNtdmFsY29yZS9wcmVwcm9jZXNzb3IvX2lvLnB5) | 75.00% | [1 Missing :warning: ](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2343?src=pr&el=tree&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    #2343      +/-   ##
==========================================
- Coverage   94.11%   94.11%   -0.01%     
==========================================
  Files         241      241              
  Lines       13376    13380       +4     
==========================================
+ Hits        12589    12592       +3     
- Misses        787      788       +1     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2343?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).

dhohn: I need help with adding tests.
bouweandela: @dhohn I simplified the code a bit and added a test for you. Could you please check if it still does what you expect?
bouweandela: It's a bit more tricky than I thought, I think the changes I made are not doing the job. Let me have another look.
bouweandela: @dhohn It should work now, could you please check?
dhohn: Dear @bouweandela 
Thanks for all these changes! Your code looks much better than mine. I tested this version it works as expected.

Patch:
diff --git a/esmvalcore/preprocessor/_io.py b/esmvalcore/preprocessor/_io.py
index 900d02694..649bb577a 100644
--- a/esmvalcore/preprocessor/_io.py
+++ b/esmvalcore/preprocessor/_io.py
@@ -20,6 +20,7 @@ from iris.cube import CubeList
 
 from esmvalcore.cmor.check import CheckLevels
 from esmvalcore.iris_helpers import merge_cube_attributes
+from esmvalcore.esgf.facets import FACETS
 
 from .._task import write_ncl_settings
 
@@ -306,6 +307,36 @@ def _sort_cubes_by_time(cubes):
     return cubes
 
 
+def _concatenate_cubes_by_experiment(
+    cubes: list[iris.cube.Cube],
+) -> list[iris.cube.Cube]:
+    """Concatenate cubes by experiment.
+
+    This ensures overlapping (branching) experiments are handled correctly.
+    """
+    # get the possible facet names in CMIP3, 5, 6 for exp
+    # currently these are 'experiment', 'experiment_id'
+    exp_facet_names = {
+        project["exp"] for project in FACETS.values() if "exp" in project
+    }
+
+    def get_exp(cube: iris.cube.Cube) -> str:
+        for key in exp_facet_names:
+            if key in cube.attributes:
+                return cube.attributes[key]
+        return ""
+
+    experiments = {get_exp(cube) for cube in cubes}
+    if len(experiments) > 1:
+        # first do experiment-wise concatenation, then time-based
+        cubes = [
+            concatenate([cube for cube in cubes if get_exp(cube) == exp])
+            for exp in experiments
+        ]
+
+    return cubes
+
+
 def concatenate(cubes, check_level=CheckLevels.DEFAULT):
     """Concatenate all cubes after fixing metadata.
 
@@ -331,6 +362,8 @@ def concatenate(cubes, check_level=CheckLevels.DEFAULT):
     if len(cubes) == 1:
         return cubes[0]
 
+    cubes = _concatenate_cubes_by_experiment(cubes)
+
     merge_cube_attributes(cubes)
     cubes = _sort_cubes_by_time(cubes)
     _fix_calendars(cubes)


Test patch:
diff --git a/tests/integration/preprocessor/_io/test_concatenate.py b/tests/integration/preprocessor/_io/test_concatenate.py
index d8d5d680c..7076f3a60 100644
--- a/tests/integration/preprocessor/_io/test_concatenate.py
+++ b/tests/integration/preprocessor/_io/test_concatenate.py
@@ -11,6 +11,7 @@ from iris.cube import Cube, CubeList
 
 from esmvalcore.cmor.check import CheckLevels
 from esmvalcore.preprocessor import _io
+from tests import assert_array_equal
 
 
 def get_hybrid_pressure_cube():
@@ -253,6 +254,33 @@ class TestConcatenate(unittest.TestCase):
             concatenated.coord('time').points, np.array([1., 2., 5., 7.,
                                                          100.]))
 
+    def test_concatenate_by_experiment_first(self):
+        """Test that data from experiments does not get mixed."""
+        historical_1 = Cube(
+            np.zeros(2),
+            dim_coords_and_dims=([
+                DimCoord(np.arange(2),
+                         var_name='time',
+                         standard_name='time',
+                         units='days since 1950-01-01'), 0
+            ], ),
+            attributes={'experiment_id': 'historical'},
+        )
+        historical_2 = historical_1.copy()
+        historical_2.coord('time').points = np.arange(2, 4)
+        historical_3 = historical_1.copy()
+        historical_3.coord('time').points = np.arange(4, 6)
+        ssp585_1 = historical_1.copy(np.ones(2))
+        ssp585_1.coord('time').points = np.arange(3, 5)
+        ssp585_1.attributes['experiment_id'] = 'ssp585'
+        ssp585_2 = ssp585_1.copy()
+        ssp585_2.coord('time').points = np.arange(5, 7)
+        result = _io.concatenate(
+            [historical_1, historical_2, historical_3, ssp585_1, ssp585_2]
+        )
+        assert_array_equal(result.coord('time').points, np.arange(7))
+        assert_array_equal(result.data, np.array([0, 0, 0, 1, 1, 1, 1]))
+
     def test_concatenate_differing_attributes(self):
         """Test concatenation of cubes with different attributes."""
         cubes = CubeList(self.raw_cubes)