Repository: ESMValGroup/ESMValCore

Issue:
ICON output is shifted in time
After doing some evaluation of monthly ICON data I noticed that the ICON data is shifted by half a month. This becomes very obvious for the top-of-the-atmosphere incoming solar radiation, e.g.:

![grafik](https://github.com/ESMValGroup/ESMValCore/assets/32543114/8576a7a6-f6bb-495f-aed1-53d78323d485)

Note that other models and other OBS coincide exactly with ERA5.

Together with a colleague we found that ICON reports aggregated data at the end of the interval, e.g., the 2020 January mean is reported as 1 February 2020 ( `20200201`). Since the raw output does not contain time bounds, this is impossible to tell from the data alone.

This also affects all other time frequencies (e.g., daily data is shifted by half a day, 6hr data by 3 hours, etc.). I will add a fix that subtracts half of the corresponding time interval to the ICON CMORizer.

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2038?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) Report
> Merging [#2038](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2038?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) (fcd712b) into [main](https://app.codecov.io/gh/ESMValGroup/ESMValCore/commit/f192e400f114b482152f83941332bdd4743b020b?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) (f192e40) will **increase** coverage by `0.03%`.
> The diff coverage is `100.00%`.

```diff
@@            Coverage Diff             @@
##             main    #2038      +/-   ##
==========================================
+ Coverage   92.78%   92.81%   +0.03%     
==========================================
  Files         236      236              
  Lines       12464    12517      +53     
==========================================
+ Hits        11565    11618      +53     
  Misses        899      899              
```


| [Impacted Files](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2038?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) | Coverage Δ | |
|---|---|---|
| [esmvalcore/cmor/\_fixes/icon/icon.py](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2038?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup#diff-ZXNtdmFsY29yZS9jbW9yL19maXhlcy9pY29uL2ljb24ucHk=) | `100.00% <100.00%> (ø)` | |



Patch:
diff --git a/doc/quickstart/find_data.rst b/doc/quickstart/find_data.rst
index 948cb7f62..5a810081b 100644
--- a/doc/quickstart/find_data.rst
+++ b/doc/quickstart/find_data.rst
@@ -342,8 +342,33 @@ Please note the duplication of the name ``ICON`` in ``project`` and
 ``dataset``, which is necessary to comply with ESMValTool's data finding and
 CMORizing functionalities.
 A variable-specific default for the facet ``var_type`` is given in the extra
-facets (see next paragraph) for many variables, but this can be overwritten in
-the recipe.
+facets (see below) for many variables, but this can be overwritten in the
+recipe.
+This is necessary if your ICON output is structured in one variable per file.
+For example, if your output is stored in files called
+``<exp>_<variable_name>_atm_2d_ml_YYYYMMDDThhmmss.nc``, use ``var_type:
+<variable_name>_atm_2d_ml`` in the recipe for this variable.
+
+Usually, ICON reports aggregated values at the end of the corresponding time
+output intervals.
+For example, for monthly output, ICON reports the month February as "1 March".
+Thus, by default, ESMValCore shifts all time points back by 1/2 of the output
+time interval so that the new time point corresponds to the center of the
+interval.
+This can be disabled by using ``shift_time: false`` in the recipe or the extra
+facets (see below).
+For point measurements (identified by ``cell_methods = "time: point"``), this
+is always disabled.
+
+.. warning::
+
+   To get all desired time points, do **not** use ``start_year`` and
+   ``end_year`` in the recipe, but rather ``timerange`` with at least 8 digits.
+   For example, to get data for the years 2000 and 2001, use ``timerange:
+   20000101/20020101`` instead of ``timerange: 2000/2001`` or ``start_year:
+   2000``, ``end_year: 2001``.
+   See :ref:`timerange_examples` for more information on the ``timerange``
+   option.
 
 ESMValCore can automatically make native ICON data `UGRID
 <https://ugrid-conventions.github.io/ugrid-conventions/>`__-compliant when
@@ -429,6 +454,9 @@ Key                 Description                      Default value if not specif
 ``raw_name``        Variable name of the             CMOR variable name of the
                     variable in the raw input        corresponding variable
                     file
+``shift_time``      Shift time points back by 1/2 of ``True``
+                    the corresponding output time
+                    interval
 ``raw_units``       Units of the variable in the     If specified, the value given by
                     raw input file                   the ``units`` attribute in the
                                                      raw input file; otherwise
diff --git a/doc/recipe/overview.rst b/doc/recipe/overview.rst
index 9f3654738..e0d63dc06 100644
--- a/doc/recipe/overview.rst
+++ b/doc/recipe/overview.rst
@@ -329,7 +329,7 @@ When using the ``timerange`` tag to specify the start and end points, possible v
 +---------------------------------------+---------------------------------------------+
 | timerange                             | effect                                      |
 +=======================================+=============================================+
-| ``'1980/1982'``                       | Spans from 01/01/1980 to 31/12/1980         |
+| ``'1980/1982'``                       | Spans from 01/01/1980 to 31/12/1982         |
 +---------------------------------------+---------------------------------------------+
 | ``'198002/198205'``                   | Spans from 01/02/1980 to 31/05/1982         |
 +---------------------------------------+---------------------------------------------+
@@ -547,7 +547,7 @@ data (has ancestors ``diagnostic_1/script_a`` and ``diagnostic_2/precip``).
 Task priority
 -------------
 Tasks are assigned a priority, with tasks appearing earlier on in the recipe
-getting higher priority. The tasks will be executed sequentially or in parellel,
+getting higher priority. The tasks will be executed sequentially or in parallel,
 depending on the setting of ``max_parallel_tasks`` in the :ref:`user configuration file`.
 When there are fewer than ``max_parallel_tasks`` running, tasks will be started
 according to their priority. For obvious reasons, only tasks that are not waiting for
diff --git a/esmvalcore/cmor/_fixes/icon/icon.py b/esmvalcore/cmor/_fixes/icon/icon.py
index 6f4af1982..214cb6e99 100644
--- a/esmvalcore/cmor/_fixes/icon/icon.py
+++ b/esmvalcore/cmor/_fixes/icon/icon.py
@@ -1,13 +1,14 @@
 """On-the-fly CMORizer for ICON."""
 
 import logging
+from datetime import datetime, timedelta
 
-import cf_units
 import dask.array as da
 import iris
 import iris.util
 import numpy as np
 import pandas as pd
+from cf_units import Unit
 from iris import NameConstraint
 from iris.coords import AuxCoord, DimCoord
 from iris.cube import CubeList
@@ -272,55 +273,139 @@ class AllVars(IconFix):
 
         # Fix metadata
         time_coord = self.fix_time_metadata(cube)
-        if 'invalid_units' not in time_coord.attributes:
-            self.guess_coord_bounds(cube, time_coord)
-            return cube
 
         # If necessary, convert invalid time units of the form "day as
         # %Y%m%d.%f" to CF format (e.g., "days since 1850-01-01")
-        # ICON data has no time bounds, let's make sure we remove the bounds
-        # here (they will be added after converting the time points to the
-        # correct units)
-        time_coord.bounds = None
-        time_format = 'day as %Y%m%d.%f'
-        t_unit = time_coord.attributes.pop('invalid_units')
-        if t_unit != time_format:
-            raise ValueError(
-                f"Expected time units '{time_format}' in input file, got "
-                f"'{t_unit}'")
-        new_t_unit = cf_units.Unit('days since 1850-01-01',
-                                   calendar='proleptic_gregorian')
-
-        # New routine to convert time of daily and hourly data. The string %f
-        # (fraction of day) is not a valid format string for datetime.strptime,
-        # so we have to convert it ourselves.
-        time_str = pd.Series(time_coord.points, dtype=str)
+        if 'invalid_units' in time_coord.attributes:
+            self._fix_invalid_time_units(time_coord)
+
+        # ICON usually reports aggregated values at the end of the time period,
+        # e.g., for monthly output, ICON reports the month February as 1 March.
+        # Thus, if not disabled, shift all time points back by 1/2 of the given
+        # time period.
+        if self.extra_facets.get('shift_time', True):
+            self._shift_time_coord(cube, time_coord)
+
+        # If not already present, try to add bounds here. Usually bounds are
+        # set in _shift_time_coord.
+        self.guess_coord_bounds(cube, time_coord)
 
-        # First, extract date (year, month, day) from string and convert it to
-        # datetime object
-        year_month_day_str = time_str.str.extract(r'(\d*)\.?\d*', expand=False)
-        year_month_day = pd.to_datetime(year_month_day_str, format='%Y%m%d')
+        return cube
 
-        # Second, extract day fraction and convert it to timedelta object
-        day_float_str = time_str.str.extract(
-            r'\d*(\.\d*)', expand=False
-        ).fillna('0.0')
-        day_float = pd.to_timedelta(day_float_str.astype(float), unit='D')
+    def _shift_time_coord(self, cube, time_coord):
+        """Shift time points back by 1/2 of given time period (in-place)."""
+        # Do not modify time coordinate for point measurements
+        for cell_method in cube.cell_methods:
+            is_point_measurement = ('time' in cell_method.coord_names and
+                                    'point' in cell_method.method)
+            if is_point_measurement:
+                logger.debug(
+                    "ICON data describes point measurements: time coordinate "
+                    "will not be shifted back by 1/2 of output interval (%s)",
+                    self.extra_facets['frequency'],
+                )
+                return
+
+        # Remove bounds; they will be re-added later after shifting
+        time_coord.bounds = None
 
-        # Finally, add date and day fraction to get final datetime and convert
-        # it to correct units. Note: we also round to next second, otherwise
-        # this results in times that are off by 1s (e.g., 13:59:59 instead of
-        # 14:00:00).
-        new_datetimes = (year_month_day + day_float).round(
-            'S'
-        ).dt.to_pydatetime()
-        new_dt_points = date2num(np.array(new_datetimes), new_t_unit)
+        # For decadal, yearly and monthly data, round datetimes to closest day
+        freq = self.extra_facets['frequency']
+        if 'dec' in freq or 'yr' in freq or 'mon' in freq:
+            time_units = time_coord.units
+            time_coord.convert_units(
+                Unit('days since 1850-01-01', calendar=time_units.calendar)
+            )
+            try:
+                time_coord.points = np.around(time_coord.points)
+            except ValueError as exc:
+                error_msg = (
+                    "Cannot shift time coordinate: Rounding to closest day "
+                    "failed. Most likely you specified the wrong frequency in "
+                    "the recipe (use `frequency: <your_frequency>` to fix "
+                    "this). Alternatively, use `shift_time=false` in the "
+                    "recipe to disable this feature"
+                )
+                raise ValueError(error_msg) from exc
+            time_coord.convert_units(time_units)
+            logger.debug(
+                "Rounded ICON time coordinate to closest day for decadal, "
+                "yearly and monthly data"
+            )
 
-        time_coord.points = new_dt_points
-        time_coord.units = new_t_unit
-        self.guess_coord_bounds(cube, time_coord)
+        # Use original time points to calculate bounds (for a given point,
+        # start of bounds is previous point, end of bounds is point)
+        first_datetime = time_coord.units.num2date(time_coord.points[0])
+        previous_time_point = time_coord.units.date2num(
+            self._get_previous_timestep(first_datetime)
+        )
+        extended_time_points = np.concatenate(
+            ([previous_time_point], time_coord.points)
+        )
+        time_coord.points = (
+            np.convolve(extended_time_points, np.ones(2), 'valid') / 2.0
+        )  # running mean with window length 2
+        time_coord.bounds = np.stack(
+            (extended_time_points[:-1], extended_time_points[1:]), axis=-1
+        )
+        logger.debug(
+            "Shifted ICON time coordinate back by 1/2 of output interval (%s)",
+            self.extra_facets['frequency'],
+        )
 
-        return cube
+    def _get_previous_timestep(self, datetime_point):
+        """Get previous time step."""
+        freq = self.extra_facets['frequency']
+        year = datetime_point.year
+        month = datetime_point.month
+
+        # Invalid input
+        invalid_freq_error_msg = (
+            f"Cannot shift time coordinate: failed to determine previous time "
+            f"step for frequency '{freq}'. Use `shift_time=false` in the "
+            f"recipe to disable this feature"
+        )
+        if 'fx' in freq or 'subhr' in freq:
+            raise ValueError(invalid_freq_error_msg)
+
+        # For decadal, yearly and monthly data, the points needs to be the
+        # first of the month 00:00:00
+        if 'dec' in freq or 'yr' in freq or 'mon' in freq:
+            if datetime_point != datetime(year, month, 1):
+                raise ValueError(
+                    f"Cannot shift time coordinate: expected first of the "
+                    f"month at 00:00:00 for decadal, yearly and monthly data, "
+                    f"got {datetime_point}. Use `shift_time=false` in the "
+                    f"recipe to disable this feature"
+                )
+
+        # Decadal data
+        if 'dec' in freq:
+            return datetime_point.replace(year=year - 10)
+
+        # Yearly data
+        if 'yr' in freq:
+            return datetime_point.replace(year=year - 1)
+
+        # Monthly data
+        if 'mon' in freq:
+            new_month = (month - 2) % 12 + 1
+            new_year = year + (month - 2) // 12
+            return datetime_point.replace(year=new_year, month=new_month)
+
+        # Daily data
+        if 'day' in freq:
+            return datetime_point - timedelta(days=1)
+
+        # Hourly data
+        if 'hr' in freq:
+            (n_hours, _, _) = freq.partition('hr')
+            if not n_hours:
+                n_hours = 1
+            return datetime_point - timedelta(hours=int(n_hours))
+
+        # Unknown input
+        raise ValueError(invalid_freq_error_msg)
 
     def _fix_mesh(self, cube, mesh_idx):
         """Fix mesh."""
@@ -370,6 +455,52 @@ class AllVars(IconFix):
 
         return True
 
+    @staticmethod
+    def _fix_invalid_time_units(time_coord):
+        """Fix invalid time units (in-place)."""
+        # ICON data usually has no time bounds. To be 100% sure, we remove the
+        # bounds here (they will be added at a later stage).
+        time_coord.bounds = None
+        time_format = 'day as %Y%m%d.%f'
+        t_unit = time_coord.attributes.pop('invalid_units')
+        if t_unit != time_format:
+            raise ValueError(
+                f"Expected time units '{time_format}' in input file, got "
+                f"'{t_unit}'"
+            )
+        new_t_units = Unit(
+            'days since 1850-01-01', calendar='proleptic_gregorian'
+        )
+
+        # New routine to convert time of daily and hourly data. The string %f
+        # (fraction of day) is not a valid format string for datetime.strptime,
+        # so we have to convert it ourselves.
+        time_str = pd.Series(time_coord.points, dtype=str)
+
+        # First, extract date (year, month, day) from string and convert it to
+        # datetime object
+        year_month_day_str = time_str.str.extract(r'(\d*)\.?\d*', expand=False)
+        year_month_day = pd.to_datetime(year_month_day_str, format='%Y%m%d')
+
+        # Second, extract day fraction and convert it to timedelta object
+        day_float_str = time_str.str.extract(
+            r'\d*(\.\d*)', expand=False
+        ).fillna('0.0')
+        day_float = pd.to_timedelta(day_float_str.astype(float), unit='D')
+
+        # Finally, add date and day fraction to get final datetime and convert
+        # it to correct units. Note: we also round to next second, otherwise
+        # this results in times that are off by 1s (e.g., 13:59:59 instead of
+        # 14:00:00).
+        new_datetimes = (year_month_day + day_float).round(
+            'S'
+        ).dt.to_pydatetime()
+        new_dt_points = date2num(np.array(new_datetimes), new_t_units)
+
+        # Modify time coordinate in place
+        time_coord.points = new_dt_points
+        time_coord.units = new_t_units
+
 
 class Clwvi(IconFix):
     """Fixes for ``clwvi``."""


Test patch:
diff --git a/tests/integration/cmor/_fixes/icon/test_icon.py b/tests/integration/cmor/_fixes/icon/test_icon.py
index d0aaf0293..61a382bf6 100644
--- a/tests/integration/cmor/_fixes/icon/test_icon.py
+++ b/tests/integration/cmor/_fixes/icon/test_icon.py
@@ -8,7 +8,7 @@ import numpy as np
 import pytest
 from cf_units import Unit
 from iris import NameConstraint
-from iris.coords import AuxCoord, DimCoord
+from iris.coords import AuxCoord, CellMethod, DimCoord
 from iris.cube import Cube, CubeList
 
 import esmvalcore.cmor._fixes.icon.icon
@@ -119,6 +119,7 @@ def _get_fix(mip, short_name, fix_name, session=None):
         short_name=short_name,
     )
     extra_facets = get_extra_facets(dataset, ())
+    extra_facets['frequency'] = 'mon'
     vardef = get_var_info(project='ICON', mip=mip, short_name=short_name)
     cls = getattr(esmvalcore.cmor._fixes.icon.icon, fix_name)
     fix = cls(vardef, extra_facets=extra_facets, session=session)
@@ -190,8 +191,8 @@ def check_time(cube):
     assert time.long_name == 'time'
     assert time.units == Unit('days since 1850-01-01',
                               calendar='proleptic_gregorian')
-    np.testing.assert_allclose(time.points, [54786.0])
-    assert time.bounds is None
+    np.testing.assert_allclose(time.points, [54770.5])
+    np.testing.assert_allclose(time.bounds, [[54755.0, 54786.0]])
     assert time.attributes == {}
 
 
@@ -725,7 +726,7 @@ def test_ta_fix_no_plev_bounds(cubes_3d):
     check_lat_lon(cube)
 
 
-# Test tas (for height2m coordinate)
+# Test tas (for height2m coordinate, no mesh, no shift time)
 
 
 def test_get_tas_fix():
@@ -775,10 +776,10 @@ def test_tas_scalar_height2m_already_present(cubes_2d):
     check_heightxm(cube, 2.0)
 
 
-def test_tas_dim_height2m_already_present(cubes_2d, monkeypatch):
+def test_tas_dim_height2m_already_present(cubes_2d):
     """Test fix."""
     fix = get_allvars_fix('Amon', 'tas')
-    monkeypatch.setitem(fix.extra_facets, 'ugrid', False)
+    fix.extra_facets['ugrid'] = False
     fixed_cubes = fix.fix_metadata(cubes_2d)
 
     cube = check_tas_metadata(fixed_cubes)
@@ -825,6 +826,28 @@ def test_tas_no_mesh(cubes_2d):
     check_heightxm(cube, 2.0)
 
 
+def test_tas_no_shift_time(cubes_2d):
+    """Test fix."""
+    fix = get_allvars_fix('Amon', 'tas')
+    fix.extra_facets['shift_time'] = False
+    fixed_cubes = fix.fix_metadata(cubes_2d)
+
+    cube = check_tas_metadata(fixed_cubes)
+    check_lat_lon(cube)
+    check_heightxm(cube, 2.0)
+
+    assert cube.coords('time', dim_coords=True)
+    time = cube.coord('time', dim_coords=True)
+    assert time.var_name == 'time'
+    assert time.standard_name == 'time'
+    assert time.long_name == 'time'
+    assert time.units == Unit('days since 1850-01-01',
+                              calendar='proleptic_gregorian')
+    np.testing.assert_allclose(time.points, [54786.0])
+    assert time.bounds is None
+    assert time.attributes == {}
+
+
 # Test uas (for height10m coordinate)
 
 
@@ -938,7 +961,7 @@ def test_ch4clim_fix(cubes_regular_grid):
     cube.units = 'mol mol-1'
     cube.coord('time').units = 'no_unit'
     cube.coord('time').attributes['invalid_units'] = 'day as %Y%m%d.%f'
-    cube.coord('time').points = [18500104.0]
+    cube.coord('time').points = [18500201.0]
     cube.coord('time').long_name = 'wrong_time_name'
 
     fix = get_allvars_fix('Amon', 'ch4Clim')
@@ -959,8 +982,8 @@ def test_ch4clim_fix(cubes_regular_grid):
     assert time_coord.units == Unit(
         'days since 1850-01-01', calendar='proleptic_gregorian'
     )
-    np.testing.assert_allclose(time_coord.points, [3.0])
-    assert time_coord.bounds is None
+    np.testing.assert_allclose(time_coord.points, [15.5])
+    np.testing.assert_allclose(time_coord.bounds, [[0.0, 31.0]])
 
 
 # Test fix with empty standard_name
@@ -1227,7 +1250,7 @@ def test_only_time(monkeypatch):
     monkeypatch.setattr(fix.vardef, 'coordinates', {'time': coord_info})
 
     # Create cube with only a single dimension
-    time_coord = DimCoord([0.0, 1.0],
+    time_coord = DimCoord([0.0, 31.0],
                           var_name='time',
                           standard_name='time',
                           long_name='time',
@@ -1254,9 +1277,9 @@ def test_only_time(monkeypatch):
     assert new_time_coord.units == 'days since 1850-01-01'
 
     # Check time data
-    np.testing.assert_allclose(new_time_coord.points, [0.0, 1.0])
+    np.testing.assert_allclose(new_time_coord.points, [-15.5, 15.5])
     np.testing.assert_allclose(new_time_coord.bounds,
-                               [[-0.5, 0.5], [0.5, 1.5]])
+                               [[-31.0, 0.0], [0.0, 31.0]])
 
     # Check that no mesh has been created
     assert cube.mesh is None
@@ -1425,12 +1448,13 @@ def test_invalid_time_units(cubes_2d):
         fix.fix_metadata(cubes_2d)
 
 
-# Test fix with hourly data
+# Test fix with (sub-)hourly data
 
 
 def test_hourly_data(cubes_2d):
     """Test fix."""
     fix = get_allvars_fix('Amon', 'tas')
+    fix.extra_facets['frequency'] = '1hr'
     for cube in cubes_2d:
         cube.coord('time').points = [20041104.5833333]
 
@@ -1438,8 +1462,11 @@ def test_hourly_data(cubes_2d):
 
     cube = check_tas_metadata(fixed_cubes)
     date = cube.coord('time').units.num2date(cube.coord('time').points)
-    np.testing.assert_array_equal(date, [datetime(2004, 11, 4, 14, 0)])
-    assert cube.coord('time').bounds is None
+    date_bnds = cube.coord('time').units.num2date(cube.coord('time').bounds)
+    np.testing.assert_array_equal(date, [datetime(2004, 11, 4, 13, 30)])
+    np.testing.assert_array_equal(
+        date_bnds, [[datetime(2004, 11, 4, 13), datetime(2004, 11, 4, 14)]]
+    )
 
 
 @pytest.mark.parametrize(
@@ -1447,28 +1474,28 @@ def test_hourly_data(cubes_2d):
     [
         None,
         [
-            [20211231.9, 20220101.1],
-            [20220101.1, 20220101.6],
-            [20220101.6, 20220102.4],
+            [20211231.875, 20220101.125],
+            [20220101.125, 20220101.375],
         ],
     ],
 )
-def test_hourly_data_multiple_points(bounds, cubes_2d):
+def test_6hourly_data_multiple_points(bounds):
     """Test fix."""
     time_coord = DimCoord(
-        [20220101, 20220101.2, 20220102],
+        [20220101, 20220101.25],
         bounds=bounds,
         standard_name='time',
         attributes={'invalid_units': 'day as %Y%m%d.%f'},
     )
     cube = Cube(
-        [1, 2, 3],
+        [1, 2],
         var_name='tas',
         units='K',
         dim_coords_and_dims=[(time_coord, 0)],
     )
     cubes = CubeList([cube])
     fix = get_allvars_fix('Amon', 'tas')
+    fix.extra_facets['frequency'] = '6hr'
 
     fixed_cube = fix._fix_time(cube, cubes)
 
@@ -1476,22 +1503,298 @@ def test_hourly_data_multiple_points(bounds, cubes_2d):
     bounds = fixed_cube.coord('time').units.num2date(cube.coord('time').bounds)
     np.testing.assert_array_equal(
         points,
+        [datetime(2021, 12, 31, 21), datetime(2022, 1, 1, 3)],
+    )
+    np.testing.assert_array_equal(
+        bounds,
         [
-            datetime(2022, 1, 1, 0, 0),
-            datetime(2022, 1, 1, 4, 48),
-            datetime(2022, 1, 2, 0, 0),
+            [datetime(2021, 12, 31, 18), datetime(2022, 1, 1)],
+            [datetime(2022, 1, 1), datetime(2022, 1, 1, 6)],
         ],
     )
+
+
+def test_subhourly_data_no_shift():
+    """Test fix."""
+    time_coord = DimCoord(
+        [0.5, 1.0],
+        standard_name='time',
+        units=Unit('hours since 2022-01-01', calendar='proleptic_gregorian'),
+    )
+    cube = Cube(
+        [1, 2],
+        var_name='tas',
+        units='K',
+        dim_coords_and_dims=[(time_coord, 0)],
+    )
+    cubes = CubeList([cube])
+    fix = get_allvars_fix('Amon', 'tas')
+    fix.extra_facets['frequency'] = 'subhr'
+    fix.extra_facets['shift_time'] = False
+
+    fixed_cube = fix._fix_time(cube, cubes)
+
+    points = fixed_cube.coord('time').units.num2date(cube.coord('time').points)
+    bounds = fixed_cube.coord('time').units.num2date(cube.coord('time').bounds)
+    np.testing.assert_array_equal(
+        points,
+        [datetime(2022, 1, 1, 0, 30), datetime(2022, 1, 1, 1)],
+    )
     np.testing.assert_array_equal(
         bounds,
         [
-            [datetime(2021, 12, 31, 21, 36), datetime(2022, 1, 1, 2, 24)],
-            [datetime(2022, 1, 1, 2, 24), datetime(2022, 1, 1, 14, 24)],
-            [datetime(2022, 1, 1, 14, 24), datetime(2022, 1, 2, 9, 36)],
+            [datetime(2022, 1, 1, 0, 15), datetime(2022, 1, 1, 0, 45)],
+            [datetime(2022, 1, 1, 0, 45), datetime(2022, 1, 1, 1, 15)],
         ],
     )
 
 
+# Test _shift_time_coord
+
+
+@pytest.mark.parametrize(
+    'frequency,dt_in,dt_out,bounds',
+    [
+        (
+            'dec',
+            [(2000, 1, 1)],
+            [(1995, 1, 1)],
+            [[(1990, 1, 1), (2000, 1, 1)]],
+        ),
+        (
+            'yr',
+            [(2000, 1, 1), (2001, 1, 1)],
+            [(1999, 7, 2, 12), (2000, 7, 2)],
+            [[(1999, 1, 1), (2000, 1, 1)], [(2000, 1, 1), (2001, 1, 1)]],
+        ),
+        (
+            'mon',
+            [(2000, 1, 1)],
+            [(1999, 12, 16, 12)],
+            [[(1999, 12, 1), (2000, 1, 1)]],
+        ),
+        (
+            'mon',
+            [(2000, 11, 30, 23, 45), (2000, 12, 31, 23)],
+            [(2000, 11, 16), (2000, 12, 16, 12)],
+            [[(2000, 11, 1), (2000, 12, 1)], [(2000, 12, 1), (2001, 1, 1)]],
+        ),
+        (
+            'day',
+            [(2000, 1, 1, 12)],
+            [(2000, 1, 1)],
+            [[(1999, 12, 31, 12), (2000, 1, 1, 12)]],
+        ),
+        (
+            '6hr',
+            [(2000, 1, 5, 14), (2000, 1, 5, 20)],
+            [(2000, 1, 5, 11), (2000, 1, 5, 17)],
+            [
+                [(2000, 1, 5, 8), (2000, 1, 5, 14)],
+                [(2000, 1, 5, 14), (2000, 1, 5, 20)],
+            ],
+        ),
+        (
+            '3hr',
+            [(2000, 1, 1)],
+            [(1999, 12, 31, 22, 30)],
+            [[(1999, 12, 31, 21), (2000, 1, 1)]],
+        ),
+        (
+            '1hr',
+            [(2000, 1, 5, 14), (2000, 1, 5, 15)],
+            [(2000, 1, 5, 13, 30), (2000, 1, 5, 14, 30)],
+            [
+                [(2000, 1, 5, 13), (2000, 1, 5, 14)],
+                [(2000, 1, 5, 14), (2000, 1, 5, 15)],
+            ],
+        ),
+    ],
+)
+def test_shift_time_coord(frequency, dt_in, dt_out, bounds):
+    """Test ``_shift_time_coord``."""
+    cube = Cube(0, cell_methods=[CellMethod('mean', 'time')])
+    datetimes = [datetime(*dt) for dt in dt_in]
+    time_units = Unit('days since 1950-01-01', calendar='proleptic_gregorian')
+    time_coord = DimCoord(
+        time_units.date2num(datetimes),
+        standard_name='time',
+        var_name='time',
+        long_name='time',
+        units=time_units,
+    )
+
+    fix = get_allvars_fix('Amon', 'tas')
+    fix.extra_facets['frequency'] = frequency
+
+    fix._shift_time_coord(cube, time_coord)
+
+    dt_out = [datetime(*dt) for dt in dt_out]
+    bounds = [[datetime(*dt1), datetime(*dt2)] for (dt1, dt2) in bounds]
+    np.testing.assert_allclose(
+        time_coord.points, time_coord.units.date2num(dt_out)
+    )
+    np.testing.assert_allclose(
+        time_coord.bounds, time_coord.units.date2num(bounds)
+    )
+
+
+@pytest.mark.parametrize(
+    'frequency,dt_in',
+    [
+        ('dec', [(2000, 1, 15)]),
+        ('yr', [(2000, 1, 1), (2001, 1, 1)]),
+        ('mon', [(2000, 6, 15)]),
+        ('day', [(2000, 1, 1), (2001, 1, 2)]),
+        ('6hr', [(2000, 6, 15, 12)]),
+        ('3hr', [(2000, 1, 1, 4), (2000, 1, 1, 7)]),
+        ('1hr', [(2000, 1, 1, 4), (2000, 1, 1, 5)]),
+    ],
+)
+def test_shift_time_point_measurement(frequency, dt_in):
+    """Test ``_shift_time_coord``."""
+    cube = Cube(0, cell_methods=[CellMethod('point', 'time')])
+    datetimes = [datetime(*dt) for dt in dt_in]
+    time_units = Unit('days since 1950-01-01', calendar='proleptic_gregorian')
+    time_coord = DimCoord(
+        time_units.date2num(datetimes),
+        standard_name='time',
+        var_name='time',
+        long_name='time',
+        units=time_units,
+    )
+
+    fix = get_allvars_fix('Amon', 'tas')
+    fix.extra_facets['frequency'] = frequency
+
+    fix._shift_time_coord(cube, time_coord)
+
+    np.testing.assert_allclose(
+        time_coord.points, time_coord.units.date2num(datetimes)
+    )
+    assert time_coord.bounds is None
+
+
+@pytest.mark.parametrize(
+    'frequency', ['dec', 'yr', 'yrPt', 'mon', 'monC', 'monPt']
+)
+def test_shift_time_coord_hourly_data_low_freq_fail(frequency):
+    """Test ``_shift_time_coord``."""
+    cube = Cube(0, cell_methods=[CellMethod('mean', 'time')])
+    time_units = Unit('hours since 1950-01-01', calendar='proleptic_gregorian')
+    time_coord = DimCoord(
+        [1, 2, 3],
+        standard_name='time',
+        var_name='time',
+        long_name='time',
+        units=time_units,
+    )
+
+    fix = get_allvars_fix('Amon', 'tas')
+    fix.extra_facets['frequency'] = frequency
+
+    msg = (
+        "Cannot shift time coordinate: Rounding to closest day failed."
+    )
+    with pytest.raises(ValueError, match=msg):
+        fix._shift_time_coord(cube, time_coord)
+
+
+@pytest.mark.parametrize(
+    'frequency', ['dec', 'yr', 'yrPt', 'mon', 'monC', 'monPt']
+)
+def test_shift_time_coord_not_first_of_month(frequency):
+    """Test ``_get_previous_timestep``."""
+    cube = Cube(0, cell_methods=[CellMethod('mean', 'time')])
+    time_units = Unit('days since 1950-01-01', calendar='proleptic_gregorian')
+    time_coord = DimCoord(
+        [1.5],
+        standard_name='time',
+        var_name='time',
+        long_name='time',
+        units=time_units,
+    )
+    fix = get_allvars_fix('Amon', 'tas')
+    fix.extra_facets['frequency'] = frequency
+
+    msg = (
+        "Cannot shift time coordinate: expected first of the month at 00:00:00"
+    )
+    with pytest.raises(ValueError, match=msg):
+        fix._shift_time_coord(cube, time_coord)
+
+
+@pytest.mark.parametrize('frequency', ['fx', 'subhrPt', 'invalid_freq'])
+def test_shift_time_coord_invalid_freq(frequency):
+    """Test ``_get_previous_timestep``."""
+    cube = Cube(0, cell_methods=[CellMethod('mean', 'time')])
+    time_units = Unit('days since 1950-01-01', calendar='proleptic_gregorian')
+    time_coord = DimCoord(
+        [1.5, 2.5],
+        standard_name='time',
+        var_name='time',
+        long_name='time',
+        units=time_units,
+    )
+    fix = get_allvars_fix('Amon', 'tas')
+    fix.extra_facets['frequency'] = frequency
+
+    msg = (
+        "Cannot shift time coordinate: failed to determine previous time step"
+    )
+    with pytest.raises(ValueError, match=msg):
+        fix._shift_time_coord(cube, time_coord)
+
+
+# Test _get_previous_timestep
+
+
+@pytest.mark.parametrize(
+    'frequency,datetime_in,datetime_out',
+    [
+        ('dec', (2000, 1, 1), (1990, 1, 1)),
+        ('yr', (2000, 1, 1), (1999, 1, 1)),
+        ('yrPt', (2001, 6, 1), (2000, 6, 1)),
+        ('mon', (2001, 1, 1), (2000, 12, 1)),
+        ('mon', (2001, 2, 1), (2001, 1, 1)),
+        ('mon', (2001, 3, 1), (2001, 2, 1)),
+        ('mon', (2001, 4, 1), (2001, 3, 1)),
+        ('monC', (2000, 5, 1), (2000, 4, 1)),
+        ('monC', (2000, 6, 1), (2000, 5, 1)),
+        ('monC', (2000, 7, 1), (2000, 6, 1)),
+        ('monC', (2000, 8, 1), (2000, 7, 1)),
+        ('monPt', (2002, 9, 1), (2002, 8, 1)),
+        ('monPt', (2002, 10, 1), (2002, 9, 1)),
+        ('monPt', (2002, 11, 1), (2002, 10, 1)),
+        ('monPt', (2002, 12, 1), (2002, 11, 1)),
+        ('day', (2000, 1, 1), (1999, 12, 31)),
+        ('day', (2000, 3, 1), (2000, 2, 29)),
+        ('day', (2187, 3, 14), (2187, 3, 13)),
+        ('6hr', (2000, 3, 14, 15), (2000, 3, 14, 9)),
+        ('6hrPt', (2000, 1, 1), (1999, 12, 31, 18)),
+        ('6hrCM', (2000, 1, 1, 1), (1999, 12, 31, 19)),
+        ('3hr', (2000, 3, 14, 15), (2000, 3, 14, 12)),
+        ('3hrPt', (2000, 1, 1), (1999, 12, 31, 21)),
+        ('3hrCM', (2000, 1, 1, 1), (1999, 12, 31, 22)),
+        ('1hr', (2000, 3, 14, 15), (2000, 3, 14, 14)),
+        ('1hrPt', (2000, 1, 1), (1999, 12, 31, 23)),
+        ('1hrCM', (2000, 1, 1, 1), (2000, 1, 1)),
+        ('hr', (2000, 3, 14), (2000, 3, 13, 23)),
+    ],
+)
+def test_get_previous_timestep(frequency, datetime_in, datetime_out):
+    """Test ``_get_previous_timestep``."""
+    datetime_in = datetime(*datetime_in)
+    datetime_out = datetime(*datetime_out)
+
+    fix = get_allvars_fix('Amon', 'tas')
+    fix.extra_facets['frequency'] = frequency
+
+    new_datetime = fix._get_previous_timestep(datetime_in)
+
+    assert new_datetime == datetime_out
+
+
 # Test mesh creation raises warning because bounds do not match vertices