Coverage for src/CSET/operators/fluxes.py: 75%
86 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-06 12:50 +0000
1# © Crown copyright, Met Office (2022-2025) and CSET contributors.
2#
3# Licensed under the Apache License, Version 2.0 (the "License");
4# you may not use this file except in compliance with the License.
5# You may obtain a copy of the License at
6#
7# http://www.apache.org/licenses/LICENSE-2.0
8#
9# Unless required by applicable law or agreed to in writing, software
10# distributed under the License is distributed on an "AS IS" BASIS,
11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12# See the License for the specific language governing permissions and
13# limitations under the License.
15"""Operators to calculate kinematic heat fluxes from covariances."""
17import iris
18import iris.cube
19from cf_units import Unit
20from iris.cube import Cube, CubeList
22from CSET._common import iter_maybe
23from CSET.operators._atmospheric_constants import CPD, LV, RD
26def _exactly_one(matches, role):
27 if len(matches) == 0:
28 raise ValueError(f"sensible_heat_units could not identify a unique {role} cube")
29 if len(matches) > 1: 29 ↛ 30line 29 didn't jump to line 30 because the condition on line 29 was never true
30 names = [getattr(c, "var_name", None) or c.name() for c in matches]
31 raise ValueError(
32 f"sensible_heat_units found multiple possible {role} cubes: {names}"
33 )
34 return matches[0]
37def _is_p_cube(cube):
38 return (
39 cube.units is not None
40 and not cube.units.is_unknown()
41 and cube.units.is_convertible(Unit("Pa"))
42 )
45def _is_T_cube(cube):
46 if cube.units is None or cube.units.is_unknown(): 46 ↛ 47line 46 didn't jump to line 47 because the condition on line 46 was never true
47 return False
48 return cube.units.is_convertible(Unit("K")) or cube.units.is_convertible(
49 Unit("degC")
50 )
53def _is_wt_covar_cube(cube):
54 if cube.units is None or cube.units.is_unknown(): 54 ↛ 55line 54 didn't jump to line 55 because the condition on line 54 was never true
55 return False
56 # turbulence covariance may be recorded either as K m s-1
57 # or degC m s-1; these are equivalent
58 return cube.units.is_convertible(Unit("K m s-1")) or cube.units.is_convertible(
59 Unit("degC m s-1")
60 )
63def sensible_heat_flux_from_covariance(cubes, **kwargs):
64 """
65 Convert turbulent temperature covariance into sensible heat flux.
67 This operator computes surface upward sensible heat flux (SHF) from
68 temperature covariance using:
70 SHF = ρ * CPD * (w'T')
72 where air density is calculated from pressure and temperature via the
73 ideal gas law.
75 The required input cubes are identified primarily from their physical
76 units:
78 - temperature covariance (e.g. K m s-1 or degC m s-1)
79 - air temperature (convertible to K or degC)
80 - air pressure (convertible to Pa)
82 If multiple physically plausible candidates are found, CF metadata
83 (e.g. standard names) are used as a secondary disambiguation step.
84 A ValueError is raised if the required cubes cannot be uniquely
85 identified.
87 Parameters
88 ----------
89 cubes : Cube or CubeList
90 Input cube(s) containing exactly one identifiable covariance,
91 temperature and pressure cube. Additional cubes are passed through
92 unchanged.
94 **kwargs : dict, optional
95 Additional keyword arguments.
97 Returns
98 -------
99 Cube or CubeList
100 Input cubes with the pressure, temperature and covariance cubes
101 removed and a new
102 ``surface_upward_sensible_heat_flux`` cube added. Unrelated cubes
103 are passed through unchanged.
105 Notes
106 -----
107 - Pressure is converted internally to Pa and temperature to K.
108 - Covariance units of ``degC m s-1`` are treated as numerically
109 equivalent to ``K m s-1`` because temperature offsets cancel when
110 forming fluctuations.
111 - Input cubes are assumed to be physically compatible; no regridding or
112 coordinate alignment is performed.
113 - Identification is unit-based, with metadata used only to resolve
114 ambiguities.
116 Raises
117 ------
118 ValueError
119 If suitable pressure, temperature or covariance cubes cannot be
120 uniquely identified.
121 """
122 from cf_units import Unit
124 cubes = (
125 iris.cube.CubeList(cubes)
126 if not isinstance(cubes, iris.cube.CubeList)
127 else cubes
128 )
130 # Pressure cube
131 p_cand = [c for c in cubes if _is_p_cube(c)]
132 if len(p_cand) > 1:
133 preferred = [c for c in p_cand if c.name() == "barometric_pressure"]
135 if len(preferred) == 1: 135 ↛ 138line 135 didn't jump to line 138 because the condition on line 135 was always true
136 p_cand = preferred
138 print("=== pressure candidates ===")
139 for c in p_cand:
140 print(
141 "name=",
142 c.name(),
143 "var_name=",
144 c.var_name,
145 "standard_name=",
146 c.standard_name,
147 "long_name=",
148 c.long_name,
149 )
150 pressure = _exactly_one(p_cand, "pressure")
152 # Temperature cube
153 T_cand = [c for c in cubes if _is_T_cube(c)]
154 if len(T_cand) > 1: 154 ↛ 155line 154 didn't jump to line 155 because the condition on line 154 was never true
155 preferred = [c for c in T_cand if c.standard_name == "air_temperature"]
157 if len(preferred) == 1:
158 T_cand = preferred
160 temp = _exactly_one(T_cand, "temperature")
162 # Covariance cube
163 covar_cand = [c for c in cubes if _is_wt_covar_cube(c)]
164 if len(covar_cand) > 1: 164 ↛ 165line 164 didn't jump to line 165 because the condition on line 164 was never true
165 preferred = []
166 for cube in covar_cand:
167 text = " ".join(
168 str(x).lower()
169 for x in (
170 cube.standard_name,
171 cube.var_name,
172 cube.long_name,
173 cube.name(),
174 )
175 if x
176 )
177 if "wt" in text or "w't" in text:
178 preferred.append(cube)
180 if len(preferred) == 1:
181 covar_cand = preferred
183 wT = _exactly_one(covar_cand, "w'T' covariance")
185 #
186 # Unit conversions
187 #
188 temp_K = temp.copy()
189 if temp_K.units.is_convertible(Unit("degC")): 189 ↛ 192line 189 didn't jump to line 192 because the condition on line 189 was always true
190 temp_K.convert_units("K")
192 pres_Pa = pressure.copy()
193 pres_Pa.convert_units("Pa")
195 # Treat degC covariance numerically as K covariance
196 wT_cov = wT.copy()
197 if str(wT_cov.units) == "degC m s-1": 197 ↛ 198line 197 didn't jump to line 198 because the condition on line 197 was never true
198 wT_cov.units = Unit("K m s-1")
200 rho_air = pres_Pa.data / (RD * temp_K.data)
202 shf = wT_cov.copy()
203 shf.data = CPD * rho_air * wT_cov.data
204 shf.units = Unit("W m-2")
205 shf.rename("surface_upward_sensible_heat_flux")
206 shf.var_name = "surface_upward_sensible_heat_flux"
208 used_ids = {id(wT), id(temp), id(pressure)}
210 out = iris.cube.CubeList(c for c in cubes if id(c) not in used_ids)
212 out.append(shf)
214 return out[0] if len(out) == 1 else out
217def latent_heat_units(
218 cubes: Cube | CubeList,
219 **kwargs,
220) -> Cube | CubeList:
221 """
222 Convert covariance into latent heat flux units.
224 This operator converts any cube with units convertible to kg m-2 s-1
225 (i.e. water mass flux) into latent heat flux (W m-2) by multiplying
226 by a constant latent heat of vaporisation.
228 No attempt is made to distinguish between turbulent fluxes (e.g. w'q')
229 and other water mass fluxes. This generalisation seems reasonable
230 given that interpreting rainfall or dewfall, for example, as an
231 equivalent heat flux is physically meaningful.
233 This function operates on one or more Iris cubes. Any cube with
234 units convertible to mass flux (kg m-2 s-1) is multiplied by a
235 constant latent heat of vaporisation to produce a latent heat flux.
236 Cubes with incompatible, missing, or unknown units are passed through
237 unchanged.
239 Parameters
240 ----------
241 cubes : Cube or CubeList
242 Input cube(s), typically containing w'q' covariance or other flux-like
243 quantities.
245 **kwargs : dict
246 Unused; accepted for interface consistency with other operators.
248 Returns
249 -------
250 Cube or CubeList
251 Output cube(s) where:
252 - Cubes with units convertible to kg m-2 s-1 are converted to W m-2.
253 - All other cubes are returned unchanged.
254 - The return type matches the input type (single Cube or CubeList).
256 Notes
257 -----
258 - The conversion uses a fixed latent heat of vaporisation:
259 LV = 2.5 × 10^6 J kg-1
260 - In reality, Lc varies with temperature (~5% variation between -20 °C
261 and +40 °C). This dependency is currently neglected but could be
262 included in future improvements.
263 - This function does not attempt to identify specific variables; it relies
264 solely on unit convertibility to determine applicability.
265 """
266 REQUIRED_UNITS = Unit("kg m-2 s-1")
267 OUTPUT_UNITS = Unit("W m-2")
269 out = iris.cube.CubeList()
270 for cube in iter_maybe(cubes):
271 # ACT ON MASS FLUXES
272 if cube.units is None or cube.units.is_unknown():
273 out.append(cube)
274 continue
275 if not cube.units.is_convertible(REQUIRED_UNITS):
276 # e.g. if UM LE or some other diagnostic — leave untouched
277 out.append(cube)
278 continue
280 cube_a = cube.copy()
281 cube_a = cube_a * LV
282 cube_a.units = cube.units * Unit("J kg-1")
283 cube_a.convert_units(OUTPUT_UNITS)
284 out.append(cube_a)
286 return out[0] if len(out) == 1 else out