Repository: ESMValGroup/ESMValCore

Issue:
Make `bias` preprocessor work with cubes
In it's current form, the [`bias`](https://github.com/ESMValGroup/ESMValCore/blob/291ce38cd6e77176b232de2a443044cb37f3e3fd/esmvalcore/preprocessor/_bias.py#L12) preprocessor can only be used on `PreprocessorFile` instances, which makes it hard to use outside of a recipe. It would be great to make this preprocessor work with cubes, too (e.g., like [`mask_multimodel`](https://github.com/ESMValGroup/ESMValCore/blob/bc5a4ac728c920f6d6c3e3676f8bc7b70a2bee5d/esmvalcore/preprocessor/_mask.py#L525)).

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

```diff
@@            Coverage Diff             @@
##             main    #2183      +/-   ##
==========================================
+ Coverage   93.10%   93.11%   +0.01%     
==========================================
  Files         237      237              
  Lines       12826    12845      +19     
==========================================
+ Hits        11942    11961      +19     
  Misses        884      884              
```


| [Files Changed](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2183?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) | Coverage Δ | |
|---|---|---|
| [esmvalcore/\_recipe/check.py](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2183?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup#diff-ZXNtdmFsY29yZS9fcmVjaXBlL2NoZWNrLnB5) | `92.24% <100.00%> (+0.21%)` | :arrow_up: |
| [esmvalcore/\_recipe/recipe.py](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2183?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup#diff-ZXNtdmFsY29yZS9fcmVjaXBlL3JlY2lwZS5weQ==) | `99.01% <100.00%> (+<0.01%)` | :arrow_up: |
| [esmvalcore/preprocessor/\_bias.py](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2183?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup#diff-ZXNtdmFsY29yZS9wcmVwcm9jZXNzb3IvX2JpYXMucHk=) | `100.00% <100.00%> (ø)` | |


valeriupredoi: I'll have a review done tomorrow, bud
schlunma: Thanks for the review V! Merging this now since I need it to continue with https://github.com/ESMValGroup/ESMValCore/pull/2299.

Patch:
diff --git a/doc/recipe/preprocessor.rst b/doc/recipe/preprocessor.rst
index 8412ceba3..271641c85 100644
--- a/doc/recipe/preprocessor.rst
+++ b/doc/recipe/preprocessor.rst
@@ -2386,9 +2386,9 @@ The bias module contains the following preprocessor functions:
 ``bias``
 --------
 
-This function calculates biases with respect to a given reference dataset. For
-this, exactly one input dataset needs to be declared as ``reference_for_bias:
-true`` in the recipe, e.g.,
+This function calculates biases with respect to a given reference dataset.
+For this, exactly one input dataset needs to be declared as
+``reference_for_bias: true`` in the recipe, e.g.,
 
 .. code-block:: yaml
 
@@ -2400,34 +2400,37 @@ true`` in the recipe, e.g.,
        reference_for_bias: true}
 
 In the example above, ERA-Interim is used as reference dataset for the bias
-calculation. For this preprocessor, all input datasets need to have identical
-dimensional coordinates. This can for example be ensured with the preprocessors
-:func:`esmvalcore.preprocessor.regrid` and/or
-:func:`esmvalcore.preprocessor.regrid_time`.
-
-The ``bias`` preprocessor supports 4 optional arguments:
-
-   * ``bias_type`` (:obj:`str`, default: ``'absolute'``): Bias type that is
-     calculated. Can be ``'absolute'`` (i.e., calculate bias for dataset
-     :math:`X` and reference :math:`R` as :math:`X - R`) or ``relative`` (i.e,
-     calculate bias as :math:`\frac{X - R}{R}`).
-   * ``denominator_mask_threshold`` (:obj:`float`, default: ``1e-3``):
-     Threshold to mask values close to zero in the denominator (i.e., the
-     reference dataset) during the calculation of relative biases. All values
-     in the reference dataset with absolute value less than the given threshold
-     are masked out. This setting is ignored when ``bias_type`` is set to
-     ``'absolute'``. Please note that for some variables with very small
-     absolute values (e.g., carbon cycle fluxes, which are usually :math:`<
-     10^{-6}` kg m :math:`^{-2}` s :math:`^{-1}`) it is absolutely essential to
-     change the default value in order to get reasonable results.
-   * ``keep_reference_dataset`` (:obj:`bool`, default: ``False``): If
-     ``True``, keep the reference dataset in the output. If ``False``, drop the
-     reference dataset.
-   * ``exclude`` (:obj:`list` of :obj:`str`): Exclude specific datasets from
-     this preprocessor. Note that this option is only available in the recipe,
-     not when using :func:`esmvalcore.preprocessor.bias` directly (e.g., in
-     another python script). If the reference dataset has been excluded, an
-     error is raised.
+calculation.
+The reference dataset needs to be broadcastable to all other datasets.
+This supports `iris' rich broadcasting abilities
+<https://scitools-iris.readthedocs.io/en/stable/userguide/cube_maths.
+html#calculating-a-cube-anomaly>`__.
+To ensure this, the preprocessors :func:`esmvalcore.preprocessor.regrid` and/or
+:func:`esmvalcore.preprocessor.regrid_time` might be helpful.
+
+The ``bias`` preprocessor supports 4 optional arguments in the recipe:
+
+* ``bias_type`` (:obj:`str`, default: ``'absolute'``): Bias type that is
+  calculated. Can be ``'absolute'`` (i.e., calculate bias for dataset
+  :math:`X` and reference :math:`R` as :math:`X - R`) or ``relative`` (i.e.,
+  calculate bias as :math:`\frac{X - R}{R}`).
+* ``denominator_mask_threshold`` (:obj:`float`, default: ``1e-3``):
+  Threshold to mask values close to zero in the denominator (i.e., the
+  reference dataset) during the calculation of relative biases. All values
+  in the reference dataset with absolute value less than the given threshold
+  are masked out. This setting is ignored when ``bias_type`` is set to
+  ``'absolute'``. Please note that for some variables with very small
+  absolute values (e.g., carbon cycle fluxes, which are usually :math:`<
+  10^{-6}` kg m :math:`^{-2}` s :math:`^{-1}`) it is absolutely essential to
+  change the default value in order to get reasonable results.
+* ``keep_reference_dataset`` (:obj:`bool`, default: ``False``): If
+  ``True``, keep the reference dataset in the output. If ``False``, drop the
+  reference dataset.
+* ``exclude`` (:obj:`list` of :obj:`str`): Exclude specific datasets from
+  this preprocessor. Note that this option is only available in the recipe,
+  not when using :func:`esmvalcore.preprocessor.bias` directly (e.g., in
+  another python script). If the reference dataset has been excluded, an
+  error is raised.
 
 Example:
 
diff --git a/esmvalcore/_recipe/check.py b/esmvalcore/_recipe/check.py
index 4b1b36cb6..f80487ece 100644
--- a/esmvalcore/_recipe/check.py
+++ b/esmvalcore/_recipe/check.py
@@ -385,6 +385,19 @@ def differing_timeranges(timeranges, required_vars):
             "Set `timerange` to a common value.")
 
 
+def bias_type(settings: dict) -> None:
+    """Check that bias_type for bias preprocessor is valid."""
+    if 'bias' not in settings:
+        return
+    valid_options = ('absolute', 'relative')
+    user_bias_type = settings['bias'].get('bias_type', 'absolute')
+    if user_bias_type not in valid_options:
+        raise RecipeError(
+            f"Expected one of {valid_options} for `bias_type`, got "
+            f"'{user_bias_type}'"
+        )
+
+
 def reference_for_bias_preproc(products):
     """Check that exactly one reference dataset for bias preproc is given."""
     step = 'bias'
diff --git a/esmvalcore/_recipe/recipe.py b/esmvalcore/_recipe/recipe.py
index d245880fa..5ea319981 100644
--- a/esmvalcore/_recipe/recipe.py
+++ b/esmvalcore/_recipe/recipe.py
@@ -656,6 +656,7 @@ def _update_preproc_functions(settings, dataset, datasets, missing_vars):
     if dataset.facets.get('frequency') == 'fx':
         check.check_for_temporal_preprocs(settings)
     check.statistics_preprocessors(settings)
+    check.bias_type(settings)
 
 
 def _get_preprocessor_task(datasets, profiles, task_name):
diff --git a/esmvalcore/preprocessor/_bias.py b/esmvalcore/preprocessor/_bias.py
index 816baa863..a8e0f4f2c 100644
--- a/esmvalcore/preprocessor/_bias.py
+++ b/esmvalcore/preprocessor/_bias.py
@@ -1,36 +1,62 @@
 """Preprocessor functions to calculate biases from data."""
+from __future__ import annotations
+
 import logging
+from collections.abc import Iterable
+from typing import TYPE_CHECKING, Literal, Optional
 
 import dask.array as da
-import iris.cube
+from iris.cube import Cube, CubeList
 
 from ._io import concatenate
 
+if TYPE_CHECKING:
+    from esmvalcore.preprocessor import PreprocessorFile
+
 logger = logging.getLogger(__name__)
 
 
-def bias(products, bias_type='absolute', denominator_mask_threshold=1e-3,
-         keep_reference_dataset=False):
-    """Calculate biases.
+BiasType = Literal['absolute', 'relative']
+
+
+def bias(
+    products: set[PreprocessorFile] | Iterable[Cube],
+    ref_cube: Optional[Cube] = None,
+    bias_type: BiasType = 'absolute',
+    denominator_mask_threshold: float = 1e-3,
+    keep_reference_dataset: bool = False,
+) -> set[PreprocessorFile] | CubeList:
+    """Calculate biases relative to a reference dataset.
+
+    The reference dataset needs to be broadcastable to all input `products`.
+    This supports `iris' rich broadcasting abilities
+    <https://scitools-iris.readthedocs.io/en/stable/userguide/cube_maths.
+    html#calculating-a-cube-anomaly>`__. To ensure this, the preprocessors
+    :func:`esmvalcore.preprocessor.regrid` and/or
+    :func:`esmvalcore.preprocessor.regrid_time` might be helpful.
 
     Notes
     -----
-    This preprocessor requires a reference dataset. For this, exactly one input
-    dataset needs to have the facet ``reference_for_bias: true`` defined in the
-    recipe. In addition, all input datasets need to have identical dimensional
-    coordinates. This can for example be ensured with the preprocessors
-    :func:`esmvalcore.preprocessor.regrid` and/or
-    :func:`esmvalcore.preprocessor.regrid_time`.
+    The reference dataset can be specified with the `ref_cube` argument. If
+    `ref_cube` is ``None``, exactly one input dataset in the `products` set
+    needs to have the facet ``reference_for_bias: true`` defined in the recipe.
+    Please do **not** specify the option `ref_cube` when using this
+    preprocessor function in a recipe.
 
     Parameters
     ----------
-    products: set of esmvalcore.preprocessor.PreprocessorFile
-        Input datasets. Exactly one datasets needs the facet
-        ``reference_for_bias: true``.
-    bias_type: str, optional (default: 'absolute')
+    products:
+        Input datasets/cubes for which the bias is calculated relative to a
+        reference dataset/cube.
+    ref_cube:
+        Cube which is used as reference for the bias calculation. If ``None``,
+        `products` needs to be a :obj:`set` of
+        `~esmvalcore.preprocessor.PreprocessorFile` objects and exactly one
+        dataset in `products` needs the facet ``reference_for_bias: true``.
+    bias_type:
         Bias type that is calculated. Must be one of ``'absolute'`` (dataset -
         ref) or ``'relative'`` ((dataset - ref) / ref).
-    denominator_mask_threshold: float, optional (default: 1e-3)
+    denominator_mask_threshold:
         Threshold to mask values close to zero in the denominator (i.e., the
         reference dataset) during the calculation of relative biases. All
         values in the reference dataset with absolute value less than the given
@@ -40,79 +66,123 @@ def bias(products, bias_type='absolute', denominator_mask_threshold=1e-3,
         :math:`< 10^{-6}` kg m :math:`^{-2}` s :math:`^{-1}`) it is absolutely
         essential to change the default value in order to get reasonable
         results.
-    keep_reference_dataset: bool, optional (default: False)
+    keep_reference_dataset:
         If ``True``, keep the reference dataset in the output. If ``False``,
-        drop the reference dataset.
+        drop the reference dataset. Ignored if `ref_cube` is given.
 
     Returns
     -------
-    set of esmvalcore.preprocessor.PreprocessorFile
-        Output datasets.
+    set[PreprocessorFile] | CubeList
+        Output datasets/cubes. Will be a :obj:`set` of
+        :class:`~esmvalcore.preprocessor.PreprocessorFile` objects if
+        `products` is also one, a :class:`~iris.cube.CubeList` otherwise.
 
     Raises
     ------
     ValueError
         Not exactly one input datasets contains the facet
-        ``reference_for_bias: true``; ``bias_type`` is not one of
+        ``reference_for_bias: true`` if ``ref_cube=None``; ``ref_cube=None``
+        and the input products are given as iterable of
+        :class:`~iris.cube.Cube` objects; ``bias_type`` is not one of
         ``'absolute'`` or ``'relative'``.
 
     """
-    # Get reference product
-    reference_product = []
-    for product in products:
-        if product.attributes.get('reference_for_bias', False):
-            reference_product.append(product)
-    if len(reference_product) != 1:
-        raise ValueError(
-            f"Expected exactly 1 dataset with 'reference_for_bias: true', "
-            f"found {len(reference_product):d}")
-    reference_product = reference_product[0]
+    ref_product = None
+    all_cubes_given = all(isinstance(p, Cube) for p in products)
 
-    # Extract reference cube
-    # Note: For technical reasons, product objects contain the member
-    # ``cubes``, which is a list of cubes. However, this is expected to be a
-    # list with exactly one element due to the call of concatenate earlier in
-    # the preprocessing chain of ESMValTool. To make sure that this
-    # preprocessor can also be used outside the ESMValTool preprocessing chain,
-    # an additional concatenate call is added here.
-    ref_cube = concatenate(reference_product.cubes)
+    # Get reference cube if not explicitly given
+    if ref_cube is None:
+        if all_cubes_given:
+            raise ValueError(
+                "A list of Cubes is given to this preprocessor; please "
+                "specify a `ref_cube`"
+            )
+        (ref_cube, ref_product) = _get_ref(products, 'reference_for_bias')
+    else:
+        ref_product = None
+
+    # Mask reference cube appropriately for relative biases
     if bias_type == 'relative':
         ref_cube = ref_cube.copy()
-        ref_cube.data = da.ma.masked_inside(ref_cube.core_data(),
-                                            -denominator_mask_threshold,
-                                            denominator_mask_threshold)
-
-    # Iterate over all input datasets and calculate bias
+        ref_cube.data = da.ma.masked_inside(
+            ref_cube.core_data(),
+            -denominator_mask_threshold,
+            denominator_mask_threshold,
+        )
+
+    # If input is an Iterable of Cube objects, calculate bias for each element
+    if all_cubes_given:
+        cubes = [_calculate_bias(c, ref_cube, bias_type) for c in products]
+        return CubeList(cubes)
+
+    # Otherwise, iterate over all input products, calculate bias and adapt
+    # metadata and provenance information accordingly
     output_products = set()
     for product in products:
-        if product == reference_product:
+        if product == ref_product:
             continue
         cube = concatenate(product.cubes)
-        cube_metadata = cube.metadata
 
         # Calculate bias
-        if bias_type == 'absolute':
-            cube = cube - ref_cube
-            new_units = str(cube.units)
-        elif bias_type == 'relative':
-            cube = (cube - ref_cube) / ref_cube
-            new_units = '1'
-        else:
-            raise ValueError(
-                f"Expected one of ['absolute', 'relative'] for bias_type, got "
-                f"'{bias_type}'")
+        cube = _calculate_bias(cube, ref_cube, bias_type)
 
-        # Adapt cube metadata and provenance information
-        cube.metadata = cube_metadata
-        cube.units = new_units
-        product.attributes['units'] = new_units
-        product.wasderivedfrom(reference_product)
+        # Adapt metadata and provenance information
+        product.attributes['units'] = str(cube.units)
+        if ref_product is not None:
+            product.wasderivedfrom(ref_product)
 
-        product.cubes = iris.cube.CubeList([cube])
+        product.cubes = CubeList([cube])
         output_products.add(product)
 
     # Add reference dataset to output if desired
-    if keep_reference_dataset:
-        output_products.add(reference_product)
+    if keep_reference_dataset and ref_product is not None:
+        output_products.add(ref_product)
 
     return output_products
+
+
+def _get_ref(products, ref_tag: str) -> tuple[Cube, PreprocessorFile]:
+    """Get reference cube and product."""
+    ref_products = []
+    for product in products:
+        if product.attributes.get(ref_tag, False):
+            ref_products.append(product)
+    if len(ref_products) != 1:
+        raise ValueError(
+            f"Expected exactly 1 dataset with '{ref_tag}: true', found "
+            f"{len(ref_products):d}"
+        )
+    ref_product = ref_products[0]
+
+    # Extract reference cube
+    # Note: For technical reasons, product objects contain the member
+    # ``cubes``, which is a list of cubes. However, this is expected to be a
+    # list with exactly one element due to the call of concatenate earlier in
+    # the preprocessing chain of ESMValTool. To make sure that this
+    # preprocessor can also be used outside the ESMValTool preprocessing chain,
+    # an additional concatenate call is added here.
+    ref_cube = concatenate(ref_product.cubes)
+
+    return (ref_cube, ref_product)
+
+
+def _calculate_bias(cube: Cube, ref_cube: Cube, bias_type: BiasType) -> Cube:
+    """Calculate bias for a single cube relative to a reference cube."""
+    cube_metadata = cube.metadata
+
+    if bias_type == 'absolute':
+        cube = cube - ref_cube
+        new_units = cube.units
+    elif bias_type == 'relative':
+        cube = (cube - ref_cube) / ref_cube
+        new_units = '1'
+    else:
+        raise ValueError(
+            f"Expected one of ['absolute', 'relative'] for bias_type, got "
+            f"'{bias_type}'"
+        )
+
+    cube.metadata = cube_metadata
+    cube.units = new_units
+
+    return cube


Test patch:
diff --git a/tests/integration/recipe/test_recipe.py b/tests/integration/recipe/test_recipe.py
index 8b7bfedb8..eb1c26041 100644
--- a/tests/integration/recipe/test_recipe.py
+++ b/tests/integration/recipe/test_recipe.py
@@ -2767,3 +2767,110 @@ def test_invalid_stat_preproc(tmp_path, patched_datafinder, session):
     msg = "Unknown preprocessor function"
     with pytest.raises(ValueError, match=msg):
         get_recipe(tmp_path, content, session)
+
+
+def test_bias_no_ref(tmp_path, patched_datafinder, session):
+    content = dedent("""
+        preprocessors:
+          test_bias:
+            bias:
+              bias_type: relative
+              denominator_mask_threshold: 5
+
+        diagnostics:
+          diagnostic_name:
+            variables:
+              ta:
+                preprocessor: test_bias
+                project: CMIP6
+                mip: Amon
+                exp: historical
+                timerange: '20000101/20001231'
+                ensemble: r1i1p1f1
+                grid: gn
+                additional_datasets:
+                  - {dataset: CanESM5}
+                  - {dataset: CESM2}
+
+            scripts: null
+        """)
+    msg = (
+        "Expected exactly 1 dataset with 'reference_for_bias: true' in "
+        "products"
+    )
+    with pytest.raises(RecipeError) as exc:
+        get_recipe(tmp_path, content, session)
+    assert str(exc.value) == INITIALIZATION_ERROR_MSG
+    assert msg in exc.value.failed_tasks[0].message
+    assert "found 0" in exc.value.failed_tasks[0].message
+
+
+def test_bias_two_refs(tmp_path, patched_datafinder, session):
+    content = dedent("""
+        preprocessors:
+          test_bias:
+            bias:
+              bias_type: relative
+              denominator_mask_threshold: 5
+
+        diagnostics:
+          diagnostic_name:
+            variables:
+              ta:
+                preprocessor: test_bias
+                project: CMIP6
+                mip: Amon
+                exp: historical
+                timerange: '20000101/20001231'
+                ensemble: r1i1p1f1
+                grid: gn
+                additional_datasets:
+                  - {dataset: CanESM5, reference_for_bias: true}
+                  - {dataset: CESM2, reference_for_bias: true}
+
+            scripts: null
+        """)
+    msg = (
+        "Expected exactly 1 dataset with 'reference_for_bias: true' in "
+        "products"
+    )
+    with pytest.raises(RecipeError) as exc:
+        get_recipe(tmp_path, content, session)
+    assert str(exc.value) == INITIALIZATION_ERROR_MSG
+    assert msg in exc.value.failed_tasks[0].message
+    assert "found 2" in exc.value.failed_tasks[0].message
+
+
+def test_inlvaid_bias_type(tmp_path, patched_datafinder, session):
+    content = dedent("""
+        preprocessors:
+          test_bias:
+            bias:
+              bias_type: INVALID
+              denominator_mask_threshold: 5
+
+        diagnostics:
+          diagnostic_name:
+            variables:
+              ta:
+                preprocessor: test_bias
+                project: CMIP6
+                mip: Amon
+                exp: historical
+                timerange: '20000101/20001231'
+                ensemble: r1i1p1f1
+                grid: gn
+                additional_datasets:
+                  - {dataset: CanESM5}
+                  - {dataset: CESM2, reference_for_bias: true}
+
+            scripts: null
+        """)
+    msg = (
+        "Expected one of ('absolute', 'relative') for `bias_type`, got "
+        "'INVALID'"
+    )
+    with pytest.raises(RecipeError) as exc:
+        get_recipe(tmp_path, content, session)
+    assert str(exc.value) == INITIALIZATION_ERROR_MSG
+    assert exc.value.failed_tasks[0].message == msg
diff --git a/tests/unit/preprocessor/_bias/test_bias.py b/tests/unit/preprocessor/_bias/test_bias.py
index 42ab9f95a..9adab3872 100644
--- a/tests/unit/preprocessor/_bias/test_bias.py
+++ b/tests/unit/preprocessor/_bias/test_bias.py
@@ -1,10 +1,10 @@
 """Unit tests for :mod:`esmvalcore.preprocessor._bias`."""
 
 import iris
-import iris.cube
 import numpy as np
 import pytest
 from cf_units import Unit
+from iris.cube import Cube, CubeList
 
 from esmvalcore.preprocessor._bias import bias
 from tests import PreprocessorFile
@@ -42,8 +42,8 @@ def get_3d_cube(data, **cube_kwargs):
                                 var_name='lon', long_name='longitude',
                                 units='degrees_east')
     coord_specs = [(times, 0), (lats, 1), (lons, 2)]
-    cube = iris.cube.Cube(data.astype('float32'),
-                          dim_coords_and_dims=coord_specs, **cube_kwargs)
+    cube = Cube(data.astype('float32'),
+                dim_coords_and_dims=coord_specs, **cube_kwargs)
     return cube
 
 
@@ -53,7 +53,7 @@ def regular_cubes():
     cube_data = np.arange(8.0).reshape(2, 2, 2)
     cube = get_3d_cube(cube_data, standard_name='air_temperature',
                        var_name='tas', units='K')
-    return iris.cube.CubeList([cube])
+    return CubeList([cube])
 
 
 @pytest.fixture
@@ -63,48 +63,18 @@ def ref_cubes():
     cube_data[1, 1, 1] = 4.0
     cube = get_3d_cube(cube_data, standard_name='air_temperature',
                        var_name='tas', units='K')
-    return iris.cube.CubeList([cube])
+    return CubeList([cube])
 
 
-def test_no_reference_for_bias(regular_cubes, ref_cubes):
-    """Test fail when no reference_for_bias is given."""
-    products = {
-        PreprocessorFile(regular_cubes, 'A', {}),
-        PreprocessorFile(regular_cubes, 'B', {}),
-        PreprocessorFile(ref_cubes, 'REF', {}),
-    }
-    msg = "Expected exactly 1 dataset with 'reference_for_bias: true', found 0"
-    with pytest.raises(ValueError, match=msg):
-        bias(products)
+TEST_BIAS = [
+    ('absolute', [[[-2.0, -1.0], [0.0, 1.0]], [[2.0, 3.0], [4.0, 3.0]]], 'K'),
+    ('relative', [[[-1.0, -0.5], [0.0, 0.5]], [[1.0, 1.5], [2.0, 0.75]]], '1'),
+]
 
 
-def test_two_reference_for_bias(regular_cubes, ref_cubes):
-    """Test fail when two reference_for_bias is given."""
-    products = {
-        PreprocessorFile(regular_cubes, 'A', {'reference_for_bias': False}),
-        PreprocessorFile(ref_cubes, 'REF1', {'reference_for_bias': True}),
-        PreprocessorFile(ref_cubes, 'REF2', {'reference_for_bias': True}),
-    }
-    msg = "Expected exactly 1 dataset with 'reference_for_bias: true', found 2"
-    with pytest.raises(ValueError, match=msg):
-        bias(products)
-
-
-def test_invalid_bias_type(regular_cubes, ref_cubes):
-    """Test fail when invalid bias_type is given."""
-    products = {
-        PreprocessorFile(regular_cubes, 'A', {}),
-        PreprocessorFile(regular_cubes, 'B', {}),
-        PreprocessorFile(ref_cubes, 'REF', {'reference_for_bias': True}),
-    }
-    msg = (r"Expected one of \['absolute', 'relative'\] for bias_type, got "
-           r"'invalid_bias_type'")
-    with pytest.raises(ValueError, match=msg):
-        bias(products, 'invalid_bias_type')
-
-
-def test_absolute_bias(regular_cubes, ref_cubes):
-    """Test calculation of absolute bias."""
+@pytest.mark.parametrize('bias_type,data,units', TEST_BIAS)
+def test_bias_products(regular_cubes, ref_cubes, bias_type, data, units):
+    """Test calculation of bias with products."""
     ref_product = PreprocessorFile(ref_cubes, 'REF',
                                    {'reference_for_bias': True})
     products = {
@@ -112,23 +82,21 @@ def test_absolute_bias(regular_cubes, ref_cubes):
         PreprocessorFile(regular_cubes, 'B', {'dataset': 'b'}),
         ref_product,
     }
-    out_products = bias(products)
+    out_products = bias(products, bias_type=bias_type)
+
+    assert isinstance(out_products, set)
     out_dict = products_set_to_dict(out_products)
     assert len(out_dict) == 2
 
     product_a = out_dict['A']
     assert product_a.filename == 'A'
-    assert product_a.attributes == {'units': 'K', 'dataset': 'a'}
+    assert product_a.attributes == {'units': units, 'dataset': 'a'}
     assert len(product_a.cubes) == 1
     out_cube = product_a.cubes[0]
-    expected_data = [[[-2.0, -1.0],
-                      [0.0, 1.0]],
-                     [[2.0, 3.0],
-                      [4.0, 3.0]]]
-    assert_array_equal(out_cube.data, expected_data)
+    assert_array_equal(out_cube.data, data)
     assert out_cube.var_name == 'tas'
     assert out_cube.standard_name == 'air_temperature'
-    assert out_cube.units == 'K'
+    assert out_cube.units == units
     assert out_cube.dim_coords == regular_cubes[0].dim_coords
     assert out_cube.aux_coords == regular_cubes[0].aux_coords
     product_a.wasderivedfrom.assert_called_once()
@@ -136,82 +104,76 @@ def test_absolute_bias(regular_cubes, ref_cubes):
 
     product_b = out_dict['B']
     assert product_b.filename == 'B'
-    assert product_b.attributes == {'units': 'K', 'dataset': 'b'}
+    assert product_b.attributes == {'units': units, 'dataset': 'b'}
     assert len(product_b.cubes) == 1
     out_cube = product_b.cubes[0]
-    expected_data = [[[-2.0, -1.0],
-                      [0.0, 1.0]],
-                     [[2.0, 3.0],
-                      [4.0, 3.0]]]
-    assert_array_equal(out_cube.data, expected_data)
+    assert_array_equal(out_cube.data, data)
     assert out_cube.var_name == 'tas'
     assert out_cube.standard_name == 'air_temperature'
-    assert out_cube.units == 'K'
+    assert out_cube.units == units
     assert out_cube.dim_coords == regular_cubes[0].dim_coords
     assert out_cube.aux_coords == regular_cubes[0].aux_coords
     product_b.wasderivedfrom.assert_called_once()
     assert product_b.mock_ancestors == {ref_product}
 
 
-def test_relative_bias(regular_cubes, ref_cubes):
-    """Test calculation of relative bias."""
-    ref_product = PreprocessorFile(ref_cubes, 'REF',
-                                   {'reference_for_bias': True})
-    products = {
-        PreprocessorFile(regular_cubes, 'A', {'dataset': 'a'}),
-        PreprocessorFile(regular_cubes, 'B', {'dataset': 'b'}),
-        ref_product,
-    }
-    out_products = bias(products, 'relative')
-    out_dict = products_set_to_dict(out_products)
-    assert len(out_dict) == 2
+@pytest.mark.parametrize('bias_type,data,units', TEST_BIAS)
+def test_bias_cubes(regular_cubes, ref_cubes, bias_type, data, units):
+    """Test calculation of bias with cubes."""
+    ref_cube = ref_cubes[0]
+    out_cubes = bias(regular_cubes, ref_cube, bias_type=bias_type)
 
-    product_a = out_dict['A']
-    assert product_a.filename == 'A'
-    assert product_a.attributes == {'units': '1', 'dataset': 'a'}
-    assert len(product_a.cubes) == 1
-    out_cube = product_a.cubes[0]
-    expected_data = [[[-1.0, -0.5],
-                      [0.0, 0.5]],
-                     [[1.0, 1.5],
-                      [2.0, 0.75]]]
-    assert_array_equal(out_cube.data, expected_data)
+    assert isinstance(out_cubes, CubeList)
+    assert len(out_cubes) == 1
+    out_cube = out_cubes[0]
+
+    assert_array_equal(out_cube.data, data)
     assert out_cube.var_name == 'tas'
     assert out_cube.standard_name == 'air_temperature'
-    assert out_cube.units == '1'
+    assert out_cube.units == units
     assert out_cube.dim_coords == regular_cubes[0].dim_coords
     assert out_cube.aux_coords == regular_cubes[0].aux_coords
-    product_a.wasderivedfrom.assert_called_once()
-    assert product_a.mock_ancestors == {ref_product}
 
-    product_b = out_dict['B']
-    assert product_b.filename == 'B'
-    assert product_b.attributes == {'units': '1', 'dataset': 'b'}
-    assert len(product_b.cubes) == 1
-    out_cube = product_b.cubes[0]
-    expected_data = [[[-1.0, -0.5],
-                      [0.0, 0.5]],
-                     [[1.0, 1.5],
-                      [2.0, 0.75]]]
-    assert_array_equal(out_cube.data, expected_data)
+
+TEST_BIAS_BROADCASTABLE = [
+    ('absolute', [[[-2.0, -1.0], [0.0, 1.0]], [[2.0, 3.0], [4.0, 5.0]]], 'K'),
+    ('relative', [[[-1.0, -0.5], [0.0, 0.5]], [[1.0, 1.5], [2.0, 2.5]]], '1'),
+]
+
+
+@pytest.mark.parametrize('bias_type,data,units', TEST_BIAS_BROADCASTABLE)
+def test_bias_cubes_broadcastable(
+    regular_cubes, ref_cubes, bias_type, data, units
+):
+    """Test calculation of bias with cubes."""
+    ref_cube = ref_cubes[0][0]  # only select one time step
+    out_cubes = bias(regular_cubes, ref_cube, bias_type=bias_type)
+
+    assert isinstance(out_cubes, CubeList)
+    assert len(out_cubes) == 1
+    out_cube = out_cubes[0]
+
+    assert_array_equal(out_cube.data, data)
     assert out_cube.var_name == 'tas'
     assert out_cube.standard_name == 'air_temperature'
-    assert out_cube.units == '1'
+    assert out_cube.units == units
     assert out_cube.dim_coords == regular_cubes[0].dim_coords
     assert out_cube.aux_coords == regular_cubes[0].aux_coords
-    product_b.wasderivedfrom.assert_called_once()
-    assert product_b.mock_ancestors == {ref_product}
 
 
-def test_denominator_mask_threshold(regular_cubes, ref_cubes):
-    """Test denominator_mask_threshold argument."""
+def test_denominator_mask_threshold_products(regular_cubes, ref_cubes):
+    """Test denominator_mask_threshold argument with products."""
     ref_product = PreprocessorFile(ref_cubes, 'REF',
                                    {'reference_for_bias': True})
     products = {
         PreprocessorFile(regular_cubes, 'A', {'dataset': 'a'}),
         ref_product,
     }
-    out_products = bias(products, 'relative', denominator_mask_threshold=3.0)
+    out_products = bias(
+        products, bias_type='relative', denominator_mask_threshold=3.0
+    )
+
+    assert isinstance(out_products, set)
     out_dict = products_set_to_dict(out_products)
     assert len(out_dict) == 1
 
@@ -234,13 +196,44 @@ def test_denominator_mask_threshold(regular_cubes, ref_cubes):
     assert product_a.mock_ancestors == {ref_product}
 
 
-def test_keep_reference_dataset(regular_cubes, ref_cubes):
+def test_denominator_mask_threshold_cubes(regular_cubes, ref_cubes):
+    """Test denominator_mask_threshold argument with cubes."""
+    ref_cube = ref_cubes[0]
+    out_cubes = bias(
+        regular_cubes,
+        ref_cube,
+        bias_type='relative',
+        denominator_mask_threshold=3.0,
+    )
+
+    assert isinstance(out_cubes, CubeList)
+    assert len(out_cubes) == 1
+    out_cube = out_cubes[0]
+
+    expected_data = np.ma.masked_equal([[[42.0, 42.0],
+                                         [42.0, 42.0]],
+                                        [[42.0, 42.0],
+                                         [42.0, 0.75]]], 42.0)
+    assert_array_equal(out_cube.data, expected_data)
+    assert out_cube.var_name == 'tas'
+    assert out_cube.standard_name == 'air_temperature'
+    assert out_cube.units == '1'
+    assert out_cube.dim_coords == regular_cubes[0].dim_coords
+    assert out_cube.aux_coords == regular_cubes[0].aux_coords
+
+
+@pytest.mark.parametrize('bias_type', ['absolute', 'relative'])
+def test_keep_reference_dataset(regular_cubes, ref_cubes, bias_type):
     """Test denominator_mask_threshold argument."""
     products = {
         PreprocessorFile(regular_cubes, 'A', {'dataset': 'a'}),
         PreprocessorFile(ref_cubes, 'REF', {'reference_for_bias': True})
     }
-    out_products = bias(products, keep_reference_dataset=True)
+    out_products = bias(
+        products, bias_type=bias_type, keep_reference_dataset=True
+    )
+
+    assert isinstance(out_products, set)
     out_dict = products_set_to_dict(out_products)
     assert len(out_dict) == 2
 
@@ -249,13 +242,92 @@ def test_keep_reference_dataset(regular_cubes, ref_cubes):
     assert product_ref.attributes == {'reference_for_bias': True}
     assert len(product_ref.cubes) == 1
     out_cube = product_ref.cubes[0]
-    expected_data = [[[2.0, 2.0],
-                      [2.0, 2.0]],
-                     [[2.0, 2.0],
-                      [2.0, 4.0]]]
+    expected_data = [[[2.0, 2.0], [2.0, 2.0]], [[2.0, 2.0], [2.0, 4.0]]]
     assert_array_equal(out_cube.data, expected_data)
     assert out_cube.var_name == 'tas'
     assert out_cube.standard_name == 'air_temperature'
     assert out_cube.units == 'K'
     assert out_cube.dim_coords == ref_cubes[0].dim_coords
     assert out_cube.aux_coords == ref_cubes[0].aux_coords
+
+
+@pytest.mark.parametrize('bias_type,data,units', TEST_BIAS)
+@pytest.mark.parametrize('keep_ref', [True, False])
+def test_bias_products_and_ref_cube(
+    regular_cubes, ref_cubes, keep_ref, bias_type, data, units
+):
+    """Test calculation of bias with products and ref_cube given."""
+    ref_cube = ref_cubes[0]
+    products = set([PreprocessorFile(regular_cubes, 'A', {'dataset': 'a'})])
+
+    out_products = bias(
+        products,
+        ref_cube=ref_cube,
+        bias_type=bias_type,
+        keep_reference_dataset=keep_ref,
+    )
+
+    assert isinstance(out_products, set)
+    out_dict = products_set_to_dict(out_products)
+    assert len(out_dict) == 1
+
+    product_a = out_dict['A']
+    assert product_a.filename == 'A'
+    assert product_a.attributes == {'units': units, 'dataset': 'a'}
+    assert len(product_a.cubes) == 1
+    out_cube = product_a.cubes[0]
+    assert_array_equal(out_cube.data, data)
+    assert out_cube.var_name == 'tas'
+    assert out_cube.standard_name == 'air_temperature'
+    assert out_cube.units == units
+    assert out_cube.dim_coords == regular_cubes[0].dim_coords
+    assert out_cube.aux_coords == regular_cubes[0].aux_coords
+    product_a.wasderivedfrom.assert_not_called()
+    assert product_a.mock_ancestors == set()
+
+
+def test_no_reference_for_bias(regular_cubes, ref_cubes):
+    """Test fail when no reference_for_bias is given."""
+    products = {
+        PreprocessorFile(regular_cubes, 'A', {}),
+        PreprocessorFile(regular_cubes, 'B', {}),
+        PreprocessorFile(ref_cubes, 'REF', {}),
+    }
+    msg = "Expected exactly 1 dataset with 'reference_for_bias: true', found 0"
+    with pytest.raises(ValueError, match=msg):
+        bias(products)
+
+
+def test_two_references_for_bias(regular_cubes, ref_cubes):
+    """Test fail when two reference_for_bias products is given."""
+    products = {
+        PreprocessorFile(regular_cubes, 'A', {'reference_for_bias': False}),
+        PreprocessorFile(ref_cubes, 'REF1', {'reference_for_bias': True}),
+        PreprocessorFile(ref_cubes, 'REF2', {'reference_for_bias': True}),
+    }
+    msg = "Expected exactly 1 dataset with 'reference_for_bias: true', found 2"
+    with pytest.raises(ValueError, match=msg):
+        bias(products)
+
+
+def test_invalid_bias_type(regular_cubes, ref_cubes):
+    """Test fail when invalid bias_type is given."""
+    products = {
+        PreprocessorFile(regular_cubes, 'A', {}),
+        PreprocessorFile(regular_cubes, 'B', {}),
+        PreprocessorFile(ref_cubes, 'REF', {'reference_for_bias': True}),
+    }
+    msg = (r"Expected one of \['absolute', 'relative'\] for bias_type, got "
+           r"'invalid_bias_type'")
+    with pytest.raises(ValueError, match=msg):
+        bias(products, bias_type='invalid_bias_type')
+
+
+def test_ref_cube_non_cubes(regular_cubes):
+    """Test ref_cube=None with with cubes."""
+    msg = (
+        "A list of Cubes is given to this preprocessor; please specify a "
+        "`ref_cube`"
+    )
+    with pytest.raises(ValueError, match=msg):
+        bias(regular_cubes)