Add functions to convert covariances into heat fluxes

src/CSET/operators/__init__.py

@@ -35,6 +35,7 @@
     ensembles,
     feature,
     filters,
+    fluxes,
     humidity,
     imageprocessing,
     mesoscale,
@@ -64,6 +65,7 @@
     "execute_recipe",
     "feature",
     "filters",
+    "fluxes",
     "humidity",
     "get_operator",
     "imageprocessing",

src/CSET/operators/fluxes.py

@@ -0,0 +1,159 @@
+# © Crown copyright, Met Office (2022-2025) and CSET contributors.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Operators to calculate kinematic heat fluxes from covariances."""
+
+import iris
+import iris.cube
+from cf_units import Unit
+from iris.cube import Cube, CubeList
+
+from CSET._common import iter_maybe
+from CSET.operators._atmospheric_constants import CPD, LV, RD
+
+
+def sensible_heat_flux_from_covariance(
+    wt_flux: iris.cube.Cube | iris.cube.CubeList,
+    air_temperature: iris.cube.Cube | iris.cube.CubeList,
+    pressure: iris.cube.Cube | iris.cube.CubeList,
+) -> iris.cube.CubeList:
+    """
+    Convert turbulent temperature covariance into sensible heat flux.
+
+    Computes:
+        SHF = rho * CPD * w'T'
+
+    where:
+        rho = pressure / (RD * temperature)
+
+    Parameters
+    ----------
+    wt_flux : iris.cube.Cube or iris.cube.CubeList
+        Turbulent temperature covariance (w'T').
+
+    air_temperature : iris.cube.Cube or iris.cube.CubeList
+        Air temperature.
+
+    pressure : iris.cube.Cube or iris.cube.CubeList
+        Air pressure.
+
+    Returns
+    -------
+    iris.cube.CubeList
+        Surface upward sensible heat flux cube(s).
+    """
+    from cf_units import Unit
+
+    out = iris.cube.CubeList()
+    for wt_cube, temp_cube, pres_cube in zip(
+        iter_maybe(wt_flux),
+        iter_maybe(air_temperature),
+        iter_maybe(pressure),
+        strict=True,
+    ):
+        # Unit conversions
+        temp_K = temp_cube.copy()
+        temp_K.convert_units("K")
+        pres_Pa = pres_cube.copy()
+        pres_Pa.convert_units("Pa")
+        wt_cov = wt_cube.copy()
+
+        # Treat degC covariance numerically as K covariance
+        if str(wt_cov.units) == "degC m s-1":
+            wt_cov.units = Unit("K m s-1")
+
+        # Density and SHF
+        rho_air = pres_Pa / (RD * temp_K)
+        shf = CPD * rho_air * wt_cov
+
+        shf.units = Unit("W m-2")
+        shf.rename("surface_upward_sensible_heat_flux")
+        # Keep compatibility with existing tests/output
+        shf.var_name = "surface_upward_sensible_heat_flux"
+
+        out.append(shf)
+
+    return out[0] if len(out) == 1 else out
+
+
+def latent_heat_units(
+    cubes: Cube | CubeList,
+    **kwargs,
+) -> Cube | CubeList:
+    """
+    Convert covariance into latent heat flux units.
+
+    This operator converts any cube with units convertible to kg m-2 s-1
+    (i.e. water mass flux) into latent heat flux (W m-2) by multiplying
+    by a constant latent heat of vaporisation.
+
+    No attempt is made to distinguish between turbulent fluxes (e.g. w'q')
+    and other water mass fluxes. This generalisation seems reasonable
+    given that interpreting rainfall or dewfall, for example, as an
+    equivalent heat flux is physically meaningful.
+
+    This function operates on one or more Iris cubes. Any cube with
+    units convertible to mass flux (kg m-2 s-1) is multiplied by a
+    constant latent heat of vaporisation to produce a latent heat flux.
+    Cubes with incompatible, missing, or unknown units are passed through
+    unchanged.
+
+    Parameters
+    ----------
+    cubes : Cube or CubeList
+        Input cube(s), typically containing w'q' covariance or other flux-like
+        quantities.
+
+    **kwargs : dict
+        Unused; accepted for interface consistency with other operators.
+
+    Returns
+    -------
+    Cube or CubeList
+        Output cube(s) where:
+        - Cubes with units convertible to kg m-2 s-1 are converted to W m-2.
+        - All other cubes are returned unchanged.
+        - The return type matches the input type (single Cube or CubeList).
+
+    Notes
+    -----
+    - The conversion uses a fixed latent heat of vaporisation:
+          LV = 2.5 × 10^6 J kg-1
+    - In reality, Lc varies with temperature (~5% variation between -20 °C
+      and +40 °C). This dependency is currently neglected but could be
+      included in future improvements.
+    - This function does not attempt to identify specific variables; it relies
+      solely on unit convertibility to determine applicability.
+    """
+    REQUIRED_UNITS = Unit("kg m-2 s-1")
+    OUTPUT_UNITS = Unit("W m-2")
+
+    out = iris.cube.CubeList()
+    for cube in iter_maybe(cubes):
+        # ACT ON MASS FLUXES
+        if cube.units is None or cube.units.is_unknown():
+            out.append(cube)
+            continue
+        if not cube.units.is_convertible(REQUIRED_UNITS):
+            # e.g. if UM LE or some other diagnostic — leave untouched
+            out.append(cube)
+            continue
+
+        cube_a = cube.copy()
+        cube_a = cube_a * LV
+        cube_a.units = cube.units * Unit("J kg-1")
+        cube_a.convert_units(OUTPUT_UNITS)
+        out.append(cube_a)
+
+    return out[0] if len(out) == 1 else out

tests/operators/test_fluxes.py

@@ -0,0 +1,180 @@
+# © Crown copyright, Met Office (2022-2025) and CSET contributors.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Test turbulent flux operators."""
+
+import iris
+import iris.coords
+import iris.cube
+import numpy as np
+from cf_units import Unit
+
+from CSET.operators import fluxes
+from CSET.operators._atmospheric_constants import CPD, LV, RD
+
+
+def _make_scalar_cube(
+    value,
+    var_name,
+    units=None,
+    standard_name=None,
+    long_name=None,
+):
+    """Make tiny 1x1 cube with metadata and units."""
+    data = np.array([[value]], dtype=np.float64)
+    lat = iris.coords.DimCoord([0.0], standard_name="latitude", units="degrees")
+    lon = iris.coords.DimCoord([0.0], standard_name="longitude", units="degrees")
+
+    cube = iris.cube.Cube(
+        data,
+        dim_coords_and_dims=[(lat, 0), (lon, 1)],
+        units=units if units is not None else Unit("unknown"),
+        standard_name=standard_name,
+        long_name=long_name,
+    )
+    cube.var_name = var_name
+    return cube
+
+
+def test_sensible_heat_flux_core_calculation():
+    """Compute sensible heat flux from covariance, temperature and pressure."""
+    wT = _make_scalar_cube(
+        0.1,
+        "wt_covariance",
+        units=Unit("K m s-1"),
+    )
+    temp = _make_scalar_cube(
+        20.0,
+        "air_temperature",
+        units=Unit("degC"),
+        standard_name="air_temperature",
+    )
+    pressure = _make_scalar_cube(
+        1000.0,
+        "surface_air_pressure",
+        units=Unit("hPa"),
+    )
+    out = fluxes.sensible_heat_flux_from_covariance(
+        wt_flux=wT,
+        air_temperature=temp,
+        pressure=pressure,
+    )
+    arr = out.data
+    if hasattr(arr, "compute"):
+        arr = arr.compute()
+
+    T = 20.0 + 273.15
+    pPa = 1000.0 * 100.0
+    rho = pPa / (RD * T)
+
+    expected = CPD * rho * 0.1
+    assert np.isclose(arr[0, 0], expected)
+    assert out.name() == "surface_upward_sensible_heat_flux"
+    assert out.units == Unit("W m-2")
+
+
+def test_sensible_heat_flux_accepts_degc_covariance():
+    """DegC m s-1 covariance should be treated as K m s-1."""
+    wT = _make_scalar_cube(
+        0.1,
+        "wt_covariance_2m",
+        units=Unit("degC m s-1"),
+    )
+    temp = _make_scalar_cube(
+        20.0,
+        "air_temperature_rtd_1p2m",
+        units=Unit("degC"),
+    )
+    pressure = _make_scalar_cube(
+        1000.0,
+        "surface_pressure",
+        units=Unit("hPa"),
+    )
+    out = fluxes.sensible_heat_flux_from_covariance(
+        wt_flux=wT,
+        air_temperature=temp,
+        pressure=pressure,
+    )
+    assert out.units == Unit("W m-2")
+
+
+def _make_scalar_latent_cube(
+    value,
+    var_name,
+    units=None,
+):
+    """Make tiny 1x1 cube with var_name and units."""
+    data = np.array([[value]], dtype=np.float64)
+    lat = iris.coords.DimCoord([0.0], standard_name="latitude", units="degrees")
+    lon = iris.coords.DimCoord([0.0], standard_name="longitude", units="degrees")
+
+    cube = iris.cube.Cube(
+        data,
+        dim_coords_and_dims=[(lat, 0), (lon, 1)],
+        units=units if units is not None else Unit("unknown"),
+    )
+    cube.var_name = var_name
+    return cube
+
+
+def test_latent_heat_units_conversion():
+    """Test unit conversion of latent heat units."""
+    wq = _make_scalar_latent_cube(0.001, "wq_covariance", units=Unit("kg m-2 s-1"))
+
+    out = fluxes.latent_heat_units(wq)
+    arr = out.data
+    if hasattr(arr, "compute"):
+        arr = arr.compute()
+
+    expected = LV * 0.001
+    assert np.isclose(arr[0, 0], expected)
+    assert out.units == Unit("W m-2")
+
+
+def test_latent_heat_units_passthrough_non_convertible():
+    """Test operator returns unchanged cube if not convertible."""
+    cube = _make_scalar_latent_cube(5.0, "not_flux", units=Unit("K"))
+    out = fluxes.latent_heat_units(cube)
+    # should be unchanged
+    assert out == cube
+    assert out.units == cube.units
+
+
+def test_latent_heat_units_cubelist_mixed():
+    """Test operator with mixed cubelist."""
+    wq = _make_scalar_latent_cube(0.001, "wq_covariance", units=Unit("kg m-2 s-1"))
+    other = _make_scalar_latent_cube(10.0, "temperature", units=Unit("K"))
+    out = fluxes.latent_heat_units([wq, other])
+
+    assert isinstance(out, iris.cube.CubeList)
+    assert len(out) == 2
+    # check conversion happened only for first cube
+    converted = [c for c in out if c.units.is_convertible("W m-2")]
+    assert len(converted) == 1
+    converted = converted[0]
+    arr = converted.data
+    if hasattr(arr, "compute"):
+        arr = arr.compute()
+
+    assert np.isclose(arr[0, 0], LV * 0.001)
+    assert converted.units == "W m-2"
+    # check passthrough
+    assert any(c is other for c in out)
+
+
+def test_latent_heat_units_unknown_units_passthrough():
+    """Test operator if units unknown."""
+    cube = _make_scalar_latent_cube(1.0, "unknown", units=None)
+    out = fluxes.latent_heat_units(cube)
+    assert out is cube