Anomalies preprocessor with period 'full' removes important metadata from cube
**Describe the bug**
When the `anomalies` preprocessor function is applied with `period=full`, it removes important metadata such as the variable name, all attributes, etc. This doesn't happen for other periods like `month` or `day`. Consequently, subsequent preprocessors or diagnostics fail.

I think it's happening [here](https://github.com/ESMValGroup/ESMValCore/blob/8306a6561a011f61dec949dc56889568b1af815e/esmvalcore/preprocessor/_time.py#L463):
```
cube = cube - reference
```
 and this seems to be the intended behaviour in Iris, as it is documented [here](https://scitools.org.uk/iris/docs/latest/userguide/cube_maths.html#calculating-a-cube-anomaly). 

Perhaps EMSValTool should only subtract the data instead of the full cube?


**Example**
`print(cube)`:

```
air_temperature / (K)               (time: 1668)
     Dimension coordinates:
          time                           x
     Auxiliary coordinates:
          day_of_year                    x
          month_number                   x
     Scalar coordinates:
          height: 1.5 m
          latitude: 0.0 degrees, bound=(-90.0, 90.0) degrees
          longitude: 179.0625 degrees, bound=(-0.9375, 359.0625) degrees
     Attributes:
          Conventions: CF-1.7
          associated_files: baseURL: http://cmip-pcmdi.llnl.gov/CMIP5/dataLocation gridspecFile: gridspec_atmos_fx_ACCESS1-0_rcp45_r0i0p0.nc...
          branch_time: 732311.0;109207.0;732311.0
          cmor_version: 2.8.0
          contact: The ACCESS wiki: http://wiki.csiro.au/confluence/display/ACCESS/Home. Contact...
          experiment: RCP4.5;historical;RCP4.5
          experiment_id: rcp45;historical;rcp45
          forcing: GHG, Oz, SA, Sl, Vl, BC, OC, (GHG = CO2, N2O, CH4, CFC11, CFC12, CFC113,...
          frequency: mon
          initialization_method: 1
          institute_id: CSIRO-BOM
          institution: CSIRO (Commonwealth Scientific and Industrial Research Organisation, Australia),...
          model_id: ACCESS1-0
          modeling_realm: atmos
          parent_experiment: historical;pre-industrial control;historical
          parent_experiment_id: historical;piControl;historical
          parent_experiment_rip: r1i1p1
          physics_version: 1
          product: output
          project_id: CMIP5
          realization: 1
          references: See http://wiki.csiro.au/confluence/display/ACCESS/ACCESS+Publications
          source: ACCESS1-0 2011. Atmosphere: AGCM v1.0 (N96 grid-point, 1.875 degrees EW...
          table_id: Table Amon (27 April 2011) 9c851218e3842df9a62ef38b1e2575bb
          title: ACCESS1-0 model output prepared for CMIP5 RCP4.5;ACCESS1-0 model output...
          version_number: v20120115
     Cell methods:
          mean: time
          mean: longitude, latitude
```

```
anoms = anomalies(cube, period='full', reference=dict(start_year=1980, start_month=1, start_day=1, end_year=2009, end_month=12, end_day=31))
print(anoms)
```
gives
```
unknown / (K)                       (time: 1668)
     Dimension coordinates:
          time                           x
     Auxiliary coordinates:
          month_number                   x
     Scalar coordinates:
          height: 1.5 m
          latitude: 0.0 degrees, bound=(-90.0, 90.0) degrees
          longitude: 179.0625 degrees, bound=(-0.9375, 359.0625) degrees
```
diff --git a/esmvalcore/preprocessor/_time.py b/esmvalcore/preprocessor/_time.py
index 3338f0b51..217fdd443 100644
--- a/esmvalcore/preprocessor/_time.py
+++ b/esmvalcore/preprocessor/_time.py
@@ -3,18 +3,19 @@
 Allows for selecting data subsets using certain time bounds;
 constructing seasonal and area averages.
 """
+import copy
 import datetime
 import logging
 from warnings import filterwarnings
 
 import dask.array as da
 import iris
-import iris.cube
 import iris.coord_categorisation
+import iris.cube
 import iris.exceptions
 import iris.util
-from iris.time import PartialDateTime
 import numpy as np
+from iris.time import PartialDateTime
 
 from ._shared import get_iris_analysis_operation, operator_accept_weights
 
@@ -460,7 +461,9 @@ def anomalies(cube, period, reference=None, standardize=False):
         reference_cube = extract_time(cube, **reference)
     reference = climate_statistics(reference_cube, period=period)
     if period in ['full']:
+        metadata = copy.deepcopy(cube.metadata)
         cube = cube - reference
+        cube.metadata = metadata
         if standardize:
             cube_stddev = climate_statistics(
                 cube, operator='std_dev', period=period)
diff --git a/tests/unit/preprocessor/_time/test_time.py b/tests/unit/preprocessor/_time/test_time.py
index 159eed146..f9fd445bc 100644
--- a/tests/unit/preprocessor/_time/test_time.py
+++ b/tests/unit/preprocessor/_time/test_time.py
@@ -1,25 +1,26 @@
 """Unit tests for the :func:`esmvalcore.preprocessor._time` module."""
 
+import copy
 import unittest
-import pytest
-import tests
-
-import numpy as np
-from numpy.testing import assert_array_equal, assert_array_almost_equal
 
-from cf_units import Unit
 import iris
 import iris.coord_categorisation
 import iris.coords
+import numpy as np
+import pytest
+from cf_units import Unit
 from iris.cube import Cube
+from numpy.testing import assert_array_almost_equal, assert_array_equal
 
-from esmvalcore.preprocessor._time import (
-    extract_month, extract_season, extract_time,
-    regrid_time,
-    decadal_statistics, annual_statistics, seasonal_statistics,
-    monthly_statistics, daily_statistics, timeseries_filter,
-    climate_statistics, anomalies
-)
+import tests
+from esmvalcore.preprocessor._time import (annual_statistics, anomalies,
+                                           climate_statistics,
+                                           daily_statistics,
+                                           decadal_statistics, extract_month,
+                                           extract_season, extract_time,
+                                           monthly_statistics, regrid_time,
+                                           seasonal_statistics,
+                                           timeseries_filter)
 
 
 def _create_sample_cube():
@@ -1022,6 +1023,20 @@ def test_standardized_anomalies(period, standardize=True):
             )
 
 
+@pytest.mark.parametrize('period, reference', PARAMETERS)
+def test_anomalies_preserve_metadata(period, reference, standardize=False):
+    cube = make_map_data(number_years=2)
+    cube.var_name = "si"
+    cube.units = "m"
+    metadata = copy.deepcopy(cube.metadata)
+    result = anomalies(cube, period, reference, standardize=standardize)
+    assert result.metadata == metadata
+    for coord_cube, coord_res in zip(cube.coords(), result.coords()):
+        if coord_cube.has_bounds() and coord_res.has_bounds():
+            assert_array_equal(coord_cube.bounds, coord_res.bounds)
+        assert coord_cube == coord_res
+
+
 @pytest.mark.parametrize('period, reference', PARAMETERS)
 def test_anomalies(period, reference, standardize=False):
     cube = make_map_data(number_years=2)
['esmvalcore/preprocessor/_time.py', 'tests/unit/preprocessor/_time/test_time.py']