diff --git a/pyproject.toml b/pyproject.toml
index fa3124743..b1a6a83f2 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -66,7 +66,7 @@ markers = [
 ]
 minversion = "7"
 pythonpath = ["src"]
-testpaths = ["tests"]
+testpaths = ["tests", "utils"]
 
 [tool.coverage.run]
 branch = true
@@ -83,7 +83,7 @@ skip = "build,*.css,*.ipynb,*.js,*.html,*.svg,*.xml,.git"
 
 [tool.ruff]
 line-length = 88
-src = ["src", "test"]
+src = ["src", "test", "utils"]
 
 [tool.ruff.lint]
 extend-select = ["B", "D"]
diff --git a/src/CSET/operators/read.py b/src/CSET/operators/read.py
index 96791d383..fe64eecbc 100644
--- a/src/CSET/operators/read.py
+++ b/src/CSET/operators/read.py
@@ -413,7 +413,7 @@ def _loading_callback(cube: iris.cube.Cube, field, filename: str) -> iris.cube.C
     _lfric_time_callback(cube)
     _lfric_forecast_period_callback(cube)
     cube = _fix_no_time_coords_callback(cube)
-    _normalise_ML_varname(cube)
+    _normalise_longname(cube)
     return cube
 
 
@@ -1098,8 +1098,8 @@ def _fix_no_time_coords_callback(cube: iris.cube.Cube):
     return cube
 
 
-def _normalise_ML_varname(cube: iris.cube.Cube):
-    """Fix plev variable names to standard names."""
+def _normalise_longname(cube: iris.cube.Cube):
+    """Normalise long_name to the LFRic standard list."""
     if cube.coords("pressure"):
         if cube.name() == "x_wind":
             cube.long_name = "zonal_wind_at_pressure_levels"
@@ -1116,6 +1116,8 @@ def _normalise_ML_varname(cube: iris.cube.Cube):
             cube.long_name = "eastward_wind_at_10m"
         if cube.name() == "y_wind" and cube.var_name == "v_wind_at_10m":
             cube.long_name = "northward_wind_at_10m"
+    if cube.name() == "air_pressure_at_sea_level":
+        cube.long_name = "air_pressure_at_mean_sea_level"
 
 
 def _check_combine_point_observations(cubes: iris.cube.CubeList):
diff --git a/tests/operators/test_read.py b/tests/operators/test_read.py
index ad1b79304..c59a39580 100644
--- a/tests/operators/test_read.py
+++ b/tests/operators/test_read.py
@@ -1166,11 +1166,11 @@ def test_fix_no_time_coords_callback(cube):
     assert cube.coord("time").units == "hours since 0001-01-01 00:00:00"
 
 
-def test_normalise_ML_varname(transect_source_cube):
+def test_normalise_longname(transect_source_cube):
     """Check that pressure varname is changed."""
     cube = transect_source_cube.copy()
     cube.rename = "air_temperature"
-    read._normalise_ML_varname(cube)
+    read._normalise_longname(cube)
     assert cube.long_name == "temperature_at_pressure_levels"
 
 
diff --git a/utils/proc_reanalysis/README.md b/utils/proc_reanalysis/README.md
new file mode 100644
index 000000000..7ea2571a6
--- /dev/null
+++ b/utils/proc_reanalysis/README.md
@@ -0,0 +1,100 @@
+# process_reanalysis
+
+## About
+
+The script `process_reanalysis.py` is a utility for converting atmospheric reanalysis datasets into a forecast-style format that can be directly compared with numerical weather prediction (NWP) model forecasts.
+
+The primary motivation for this tool is to be able to analyse reanalysis alongside model forecasts in CSET, by choosing reanalysis to be the base model for verification and evaluation. CSET expects forecast data to contain forecast metadata such as:
+
+- `forecast_reference_time`
+- `forecast_period`
+- `time`
+
+Reanalysis datasets typically only contain valid time `time`, and is a series of files where the `forecast_reference_time` changes every 6 hours, and the `forecast_period` is zero.
+
+This script resolves this issue by transforming reanalysis data into an effective forecast representation. Rather than treating reanalysis as a special data source in CSET, the transformed output can be ingested directly into the CSET workflow. This allows reanalysis to be treated as another "model" within CSET.
+
+The script currently supports datasets that can be loaded by Iris and has primarily been developed and tested using:
+
+- ERA5 reanalysis
+- Unified Model (UM) analysis data
+
+Other model analyses may work, as only the time dimensions are manipulated, but this has not been tested.
+No scientific changes are made to the meteorological fields themselves.
+
+For each requested forecast cycle, for each variable found in the reanalysis the script:
+1. Extracts the required period of time from reanalysis data.
+2. Treats the start of that extraction as a forecast initialisation.
+3. Generates a forecast period coordinate.
+4. Generates a forecast reference time coordinate.
+5. Preserves the original valid time information.
+6. Saves the result as a forecast-style NetCDF file.
+
+> [!TIP]
+> This script does not download reanalysis data - a user must fetch this first from archive/api.
+
+## Usage
+
+The python script requires the Iris package to be installed and available to python.
+
+Run it with:
+
+```
+python process_reanalysis.py \
+    --files "<input_files>" \
+    --cyclestart YYYYMMDDTHHMMZ" \
+    --cycleend YYYYMMDDTHHMMZ" \
+    --cyclefreq <hours> \
+    --forecastlength <hours> \
+    --outpath "<output_directory>"
+```
+
+Required Arguments:
+
+- `--files`: Path to the input reanalysis data. This can be a single file or wildcard expression understood by Iris. If a wildcard is used, then quote the input to prevent the shell expanding the filelist as arguments to python.
+- `--cyclestart`: First forecast initialisation time that you want the reanalysis to simulate, in format <year><month><day>T<hour><minute>Z.
+- `--cycleend`: Final forecast initialisation time, inclusive, that you want the reanalysis to simulate, in format <year><month><day>T<hour><minute>Z.
+- `--cyclefreq`: Frequency between forecast cycles, in hours, as an integer.
+- `--forecastlength`: Length of the forecast you want the reanalysis to simulate, in hours, as an integer.
+- `--outpath`: Path of where to store the output data. The code will write a file per forecast initialisation, in the format of `reanalysis_%Y%m%dT%H%MZ_.nc`.
+
+## Examples
+
+1. A single forecast that goes out to 48h, initialised on the 1st January 2024 at 00Z.
+
+```
+python process_reanalysis.py \
+    --files "/data/era5/*.nc" \
+    --cyclestart "20240101T0000Z" \
+    --cycleend "20240101T0000Z" \
+    --cyclefreq 6 \
+    --forecastlength 48 \
+    --outpath /my/output/path/
+```
+Producing one file `my/output/path/reanalysis_20240101T0000Z.nc`
+
+2. Produce 6-hourly analysis across one day.
+
+```
+python process_reanalysis.py \
+    --files "/data/era5/*.nc" \
+    --cyclestart "20240101T0000Z" \
+    --cycleend "20240101T1800Z" \
+    --cyclefreq 6 \
+    --forecastlength 48 \
+    --outpath /my/output/path/
+```
+Produces
+
+```
+/my/output/path/reanalysis_20240101T0000Z.nc
+/my/output/path/reanalysis_20240101T0600Z.nc
+/my/output/path/reanalysis_20240101T1200Z.nc
+/my/output/path/reanalysis_20240101T1800Z.nc
+```
+
+## Owners
+
+The following people should be contacted for queries or issues with this utility:
+
+* [@jwarner8](https://github.com/jwarner8)
diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py
new file mode 100755
index 000000000..19d8566b9
--- /dev/null
+++ b/utils/proc_reanalysis/process_reanalysis.py
@@ -0,0 +1,273 @@
+#!/usr/bin/python3
+
+"""
+Code that restructures reanalysis data to give it an effective forecast_period.
+
+As a result, reanalysis can be directly compared to model forecasts in CSET and treated as another model.
+The code base currently supports UM and ERA5, and does not perform any additional metadata
+correction beyond time axis and removing some surplus coords/attributes.
+
+Please see README for further information on how to run the script.
+"""
+
+import iris
+import iris.cube
+
+iris.FUTURE.date_microseconds = True
+iris.FUTURE.save_split_attrs = True
+import argparse
+from datetime import datetime, timedelta
+
+
+def identify_number_of_cycles_required(
+    cyclestart: datetime, cycleend: datetime, cyclefreq: timedelta
+) -> list:
+    """Generate forecast initialisation datetimes between two cycle bounds.
+
+    Parameters
+    ----------
+    cyclestart : datetime
+        First forecast cycle time.
+    cycleend : datetime
+        Last forecast cycle time.
+    cyclefreq : timedelta
+        Frequency between forecast cycles in hours.
+
+    Returns
+    -------
+    forecast_initialisations: list
+        Forecast initialisation datetimes from start_dt to end_dt,
+        inclusive, separated by cyclefreq hours.
+    """
+    # To store initialisation times
+    forecast_initialisations = []
+    current = cyclestart
+
+    # Iterate over all initiations within the bounds, using the cyclefreq to determine interval.
+    while current <= cycleend:
+        forecast_initialisations.append(current)
+        current += cyclefreq
+
+    return forecast_initialisations
+
+
+def create_forecasts(
+    reanalysis: iris.cube.CubeList,
+    forecast_initialisations: list,
+    forecastlength: timedelta,
+    outpath: str,
+) -> None:
+    """Create forecast files from reanalysis data.
+
+    For each forecast initialisation time, extract the corresponding
+    analysis period from each input cube and convert it into a
+    forecast-style representation. This includes generating
+    forecast-period and forecast-reference-time coordinates and
+    writing the resulting cubes to disk.
+
+    Parameters
+    ----------
+    reanalysis: iris.cube.CubeList
+        Collection of reanalysis cubes from which forecast periods
+        will be extracted.
+    forecast_initialisations: list
+        Forecast initialisation times to process.
+    forecastlength: timedelta
+        Forecast length in hours.
+    outpath: str
+        Directory to which the generated forecast files will be saved.
+
+    Returns
+    -------
+    None
+    """
+    # Iterate over all forecast initialisations sequentially.
+    for forecast in forecast_initialisations:
+        print(
+            f"Working on forecast initialisation {forecast} out to {forecastlength.total_seconds() / 3600}H"
+        )
+
+        # Work out start and end time
+        start = forecast
+        end = forecast + forecastlength
+
+        cutouts = iris.cube.CubeList()
+
+        # For each cube (variable) loaded
+        for cube in reanalysis:
+            print(f"{cube.name()}...")
+
+            # Work out minimum, maximum valid times in analysis
+            an_min = cube.coord("time").units.num2date(cube.coord("time").points[0])
+            an_max = cube.coord("time").units.num2date(cube.coord("time").points[-1])
+
+            # Check reanalysis spans what we are looking for time wise, otherwise ignore
+            if start < an_min or end > an_max:
+                print(
+                    f"Warning: Required time {start} {end} outside that found in analysis {an_min} {an_max}"
+                )
+            else:
+                # Generate time constraint object inclusive of time bounds.
+                time_constraint = iris.Constraint(
+                    time=lambda cell, start=start, end=end: start <= cell.point <= end
+                )
+
+                # Extract required timeslice.
+                cube_slice = cube.extract(time_constraint)
+
+                # Remove unnecessary coords and attributes
+                coords_attrs_to_remove = [
+                    "forecast_period",
+                    "forecast_reference_time",
+                    "originating_centre",
+                    "source",
+                    "um_version",
+                ]
+                for item in coords_attrs_to_remove:
+                    if cube_slice.coords(item):
+                        cube_slice.remove_coord(item)
+                    if item in cube_slice.attributes:
+                        del cube_slice.attributes[item]
+
+                # Get a copy of time coord
+                time_coord = cube_slice.coord("time")
+
+                # Work out units of forecast_period and adjust if necessary.
+                units_str = str(time_coord.units)
+
+                # Forecast periods relative to initialisation.
+                fp_points = time_coord.points - time_coord.points[0]
+
+                if units_str.startswith("seconds since"):
+                    fp_points = fp_points / 3600.0
+                    fp_units = "hours"
+                elif units_str.startswith("minutes since"):
+                    fp_points = fp_points / 60.0
+                    fp_units = "hours"
+                elif units_str.startswith("hours since"):
+                    fp_units = "hours"
+                else:
+                    raise ValueError(f"Unhandled time units: {time_coord.units}")
+
+                # Get time points and dimension that time corresponds to.
+                time_coord_points = time_coord.points
+                time_dim = cube_slice.coord_dims("time")[0]
+
+                # Create forecast period dimension
+                fp_coord = iris.coords.DimCoord(
+                    fp_points,
+                    standard_name="forecast_period",
+                    units=fp_units,
+                )
+
+                # Remove time dimension temporarily, as forecast_period will be lead dimension
+                cube_slice.remove_coord("time")
+                cube_slice.add_dim_coord(fp_coord, time_dim)
+
+                # Add auxiliary forecast initialisation dimension
+                cube_slice.add_aux_coord(
+                    iris.coords.AuxCoord(
+                        time_coord.units.date2num(start),
+                        standard_name="forecast_reference_time",
+                        units=time_coord.units,
+                    )
+                )
+
+                # Create auxiliary time dimension of valid time.
+                new_time_coord = iris.coords.AuxCoord(
+                    time_coord_points,
+                    standard_name="time",
+                    units=time_coord.units,
+                )
+
+                # Add this dimension to the cube, tied to the forecast_period dimension.
+                cube_slice.add_aux_coord(
+                    new_time_coord,
+                    data_dims=(cube_slice.coord_dims("forecast_period")[0],),
+                )
+
+                # Append slice to cutout list, ready for saving.
+                cutouts.append(cube_slice)
+                print(f"{cube.name()}...done.")
+
+        # Once all cubes processed, save to disk.
+        if len(cutouts) > 0:
+            filename = f"{outpath}/reanalysis_{start.strftime('%Y%m%dT%H%MZ')}.nc"
+            print(f"Saving {filename}")
+            iris.save(cutouts, filename)
+        else:
+            raise ValueError("No suitable cubes found for saving!")
+
+
+def main() -> None:
+    """Generate forecast-like datasets from reanalysis data.
+
+    Parse command-line arguments, create forecast initialisation times,
+    process the input reanalysis data, and write the resulting forecast
+    files to disk.
+    """
+    parser = argparse.ArgumentParser(description="Process arguments.")
+
+    parser.add_argument("--files", required=True, help="Path to file(s) quoted")
+    parser.add_argument(
+        "--cyclestart",
+        type=datetime.fromisoformat,
+        required=True,
+        help="First forecast initiation/cycle, in format YYYYMMDDTHHMMZ",
+    )
+    parser.add_argument(
+        "--cycleend",
+        type=datetime.fromisoformat,
+        required=True,
+        help="Final forecast initiation/cycle, in format YYYYMMDDTHHMMZ",
+    )
+    parser.add_argument(
+        "--cyclefreq",
+        type=int,
+        required=True,
+        help="Hours between forecast initiations/cycles",
+    )
+    parser.add_argument(
+        "--forecastlength",
+        type=int,
+        required=True,
+        help="Forecast length in hours, i.e. 48",
+    )
+    parser.add_argument(
+        "--outpath", type=str, required=True, help="Where to write output data"
+    )
+
+    args = parser.parse_args()
+
+    # Populate required variables
+    filepath = args.files
+    cyclestart = args.cyclestart
+    cycleend = args.cycleend
+    cyclefreq = timedelta(hours=args.cyclefreq)
+    forecastlength = timedelta(hours=args.forecastlength)
+    outpath = args.outpath
+
+    print()
+    print("Starting process_reanalysis.py...")
+
+    # Get all forecast initiations
+    forecast_initialisations = identify_number_of_cycles_required(
+        cyclestart, cycleend, cyclefreq
+    )
+
+    # Load all reanalysis supplied
+    print(f"Loading reanalysis from {filepath}")
+    reanalysis = iris.load(filepath)
+    print()
+    print("Found the following cubes...")
+    print(reanalysis)
+
+    print()
+    print("Creating postprocessed files...")
+    create_forecasts(reanalysis, forecast_initialisations, forecastlength, outpath)
+
+    print("Done")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/utils/proc_reanalysis/test_process_reanalysis.py b/utils/proc_reanalysis/test_process_reanalysis.py
new file mode 100644
index 000000000..cec5b4be5
--- /dev/null
+++ b/utils/proc_reanalysis/test_process_reanalysis.py
@@ -0,0 +1,229 @@
+"""Unit tests for process_reanalysis.py."""
+
+from datetime import datetime, timedelta
+
+import iris
+import numpy as np
+import pytest
+from iris.coords import DimCoord
+from iris.cube import Cube
+from utils.proc_reanalysis import process_reanalysis as proc_reanalysis
+
+
+def test_single_cycle():
+    """Assert single datetime returned if one initialisation."""
+    result = proc_reanalysis.identify_number_of_cycles_required(
+        datetime.fromisoformat("2024-01-01 00:00:00"),
+        datetime.fromisoformat("2024-01-01 00:00:00"),
+        timedelta(hours=6),
+    )
+
+    assert result == [datetime(2024, 1, 1, 0, 0)]
+
+
+def test_multiple_cycles():
+    """Test handling of multiple cycles identified."""
+    result = proc_reanalysis.identify_number_of_cycles_required(
+        datetime.fromisoformat("2024-01-01 00:00:00"),
+        datetime.fromisoformat("2024-01-01 12:00:00"),
+        timedelta(hours=6),
+    )
+
+    assert result == [
+        datetime(2024, 1, 1, 0),
+        datetime(2024, 1, 1, 6),
+        datetime(2024, 1, 1, 12),
+    ]
+
+
+def test_non_divisible_interval():
+    """Check end point is not exceeded."""
+    result = proc_reanalysis.identify_number_of_cycles_required(
+        datetime.fromisoformat("2024-01-01 00:00:00"),
+        datetime.fromisoformat("2024-01-01 10:00:00"),
+        timedelta(hours=6),
+    )
+
+    assert result == [
+        datetime(2024, 1, 1, 0),
+        datetime(2024, 1, 1, 6),
+    ]
+
+
+def make_cube(
+    name="air_temperature",
+    units="hours since 2024-01-01 00:00:00",
+):
+    """Create a minimal synthetic reanalysis cube."""
+    time = DimCoord(
+        np.arange(5),
+        standard_name="time",
+        units=units,
+    )
+
+    lat = DimCoord([50.0], standard_name="latitude", units="degrees")
+    lon = DimCoord([0.0], standard_name="longitude", units="degrees")
+
+    data = np.arange(5).reshape(5, 1, 1)
+
+    cube = Cube(
+        data,
+        standard_name=name,
+        dim_coords_and_dims=[
+            (time, 0),
+            (lat, 1),
+            (lon, 2),
+        ],
+    )
+
+    return cube
+
+
+def test_forecast_period_created(tmp_path):
+    """Check forecast period constructed correctly."""
+    cube = make_cube()
+
+    proc_reanalysis.create_forecasts(
+        iris.cube.CubeList([cube]),
+        [datetime(2024, 1, 1, 0)],
+        forecastlength=timedelta(hours=4),
+        outpath=str(tmp_path),
+    )
+
+    outfile = tmp_path / "reanalysis_20240101T0000Z.nc"
+
+    cubes = iris.load(str(outfile))
+
+    result = cubes[0]
+
+    fp = result.coord("forecast_period")
+
+    assert fp.points.tolist() == [0, 1, 2, 3, 4]
+    assert str(fp.units) == "hours"
+
+
+def test_seconds_converted_to_hours(tmp_path):
+    """Check conversion from seconds to hours in forecast period."""
+    cube = make_cube(units="seconds since 2024-01-01 00:00:00")
+
+    cube.coord("time").points = [0, 3600, 7200, 10800, 14400]
+
+    proc_reanalysis.create_forecasts(
+        iris.cube.CubeList([cube]),
+        [datetime(2024, 1, 1)],
+        forecastlength=timedelta(hours=2),
+        outpath=str(tmp_path),
+    )
+
+    cubes = iris.load(str(tmp_path / "reanalysis_20240101T0000Z.nc"))
+
+    fp = cubes[0].coord("forecast_period")
+
+    assert fp.points.tolist() == [0, 1, 2]
+
+
+def test_minutes_converted_to_hours(tmp_path):
+    """Check conversion from minutes to hours in forecast period."""
+    cube = make_cube(units="minutes since 2024-01-01 00:00:00")
+
+    cube.coord("time").points = [0, 60, 120, 180, 240]
+
+    proc_reanalysis.create_forecasts(
+        iris.cube.CubeList([cube]),
+        [datetime(2024, 1, 1)],
+        forecastlength=timedelta(hours=2),
+        outpath=str(tmp_path),
+    )
+
+    cubes = iris.load(str(tmp_path / "reanalysis_20240101T0000Z.nc"))
+
+    fp = cubes[0].coord("forecast_period")
+
+    assert fp.points.tolist() == [0, 1, 2]
+
+
+def test_unknown_time_units_raise(tmp_path):
+    """Check that error raised if time units unhandled."""
+    cube = make_cube(units="days since 2024-01-01 00:00:00")
+
+    with pytest.raises(ValueError, match="Unhandled time units"):
+        proc_reanalysis.create_forecasts(
+            iris.cube.CubeList([cube]),
+            [datetime(2024, 1, 1)],
+            forecastlength=timedelta(hours=1),
+            outpath=str(tmp_path),
+        )
+
+
+def test_forecast_reference_time_created(tmp_path):
+    """Check that forecast reference time has been set correctly."""
+    init_time = datetime(2024, 1, 1, 0)
+
+    cube = make_cube()
+
+    proc_reanalysis.create_forecasts(
+        iris.cube.CubeList([cube]),
+        [init_time],
+        forecastlength=timedelta(hours=4),
+        outpath=str(tmp_path),
+    )
+
+    cubes = iris.load(str(tmp_path / "reanalysis_20240101T0000Z.nc"))
+
+    frt = cubes[0].coord("forecast_reference_time")
+
+    recovered = frt.units.num2date(frt.points[0])
+
+    assert str(recovered) == "2024-01-01 00:00:00"
+
+
+def test_forecast_attributes_removed(tmp_path):
+    """Check that common analysis attributes have been removed."""
+    cube = make_cube()
+
+    cube.attributes["source"] = "ERA5"
+    cube.attributes["um_version"] = "13.0"
+
+    proc_reanalysis.create_forecasts(
+        iris.cube.CubeList([cube]),
+        [datetime(2024, 1, 1)],
+        forecastlength=timedelta(hours=4),
+        outpath=str(tmp_path),
+    )
+
+    cubes = iris.load(str(tmp_path / "reanalysis_20240101T0000Z.nc"))
+
+    attrs = cubes[0].attributes
+
+    assert "source" not in attrs
+    assert "um_version" not in attrs
+
+
+def test_cube_skipped_if_insufficient_data(tmp_path):
+    """Check that nothing save if analysis doesn't overlap with target forecast."""
+    cube = make_cube()
+
+    with pytest.raises(ValueError, match="No suitable cubes found for saving!"):
+        proc_reanalysis.create_forecasts(
+            iris.cube.CubeList([cube]),
+            [datetime(2024, 1, 1)],
+            forecastlength=timedelta(hours=10),
+            outpath=str(tmp_path),
+        )
+
+
+def test_multiple_cubes_processed(tmp_path):
+    """Check working with multiple cubes."""
+    cube1 = make_cube("air_temperature")
+    cube2 = make_cube("air_pressure")
+
+    proc_reanalysis.create_forecasts(
+        iris.cube.CubeList([cube1, cube2]),
+        [datetime(2024, 1, 1)],
+        forecastlength=timedelta(hours=4),
+        outpath=str(tmp_path),
+    )
+
+    cubes = iris.load(str(tmp_path / "reanalysis_20240101T0000Z.nc"))
+
+    assert len(cubes) == 2