Coverage for src/CSET/operators/fluxes.py: 97%
60 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 12:32 +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 sensible_heat_flux_from_covariance(
27 wt_flux: iris.cube.Cube | iris.cube.CubeList,
28 air_temperature: iris.cube.Cube | iris.cube.CubeList,
29 pressure: iris.cube.Cube | iris.cube.CubeList,
30) -> iris.cube.CubeList:
31 """
32 Convert turbulent temperature covariance into sensible heat flux.
34 Computes:
35 SHF = rho * CPD * w'T'
37 where:
38 rho = pressure / (RD * temperature)
40 Parameters
41 ----------
42 wt_flux : iris.cube.Cube or iris.cube.CubeList
43 Turbulent temperature covariance (w'T').
45 air_temperature : iris.cube.Cube or iris.cube.CubeList
46 Air temperature.
48 pressure : iris.cube.Cube or iris.cube.CubeList
49 Air pressure.
51 Returns
52 -------
53 iris.cube.CubeList
54 Surface upward sensible heat flux cube(s).
55 """
56 from cf_units import Unit
58 out = iris.cube.CubeList()
59 for wt_cube, temp_cube, pres_cube in zip(
60 iter_maybe(wt_flux),
61 iter_maybe(air_temperature),
62 iter_maybe(pressure),
63 strict=True,
64 ):
65 # Unit conversions
66 temp_K = temp_cube.copy()
67 temp_K.convert_units("K")
68 pres_Pa = pres_cube.copy()
69 pres_Pa.convert_units("Pa")
70 wt_cov = wt_cube.copy()
72 # Treat degC covariance numerically as K covariance
73 if str(wt_cov.units) == "degC m s-1":
74 wt_cov.units = Unit("K m s-1")
76 # Density and SHF
77 rho_air = pres_Pa / (RD * temp_K)
78 shf = CPD * rho_air * wt_cov
80 shf.units = Unit("W m-2")
81 shf.rename("surface_upward_sensible_heat_flux")
82 # Keep compatibility with existing tests/output
83 shf.var_name = "surface_upward_sensible_heat_flux"
85 out.append(shf)
87 return out[0] if len(out) == 1 else out
90def latent_heat_units(
91 cubes: Cube | CubeList,
92 **kwargs,
93) -> Cube | CubeList:
94 """
95 Convert covariance into latent heat flux units.
97 This operator converts any cube with units convertible to kg m-2 s-1
98 (i.e. water mass flux) into latent heat flux (W m-2) by multiplying
99 by a constant latent heat of vaporisation.
101 No attempt is made to distinguish between turbulent fluxes (e.g. w'q')
102 and other water mass fluxes. This generalisation seems reasonable
103 given that interpreting rainfall or dewfall, for example, as an
104 equivalent heat flux is physically meaningful.
106 This function operates on one or more Iris cubes. Any cube with
107 units convertible to mass flux (kg m-2 s-1) is multiplied by a
108 constant latent heat of vaporisation to produce a latent heat flux.
109 Cubes with incompatible, missing, or unknown units are passed through
110 unchanged.
112 Parameters
113 ----------
114 cubes : Cube or CubeList
115 Input cube(s), typically containing w'q' covariance or other flux-like
116 quantities.
118 **kwargs : dict
119 Unused; accepted for interface consistency with other operators.
121 Returns
122 -------
123 Cube or CubeList
124 Output cube(s) where:
125 - Cubes with units convertible to kg m-2 s-1 are converted to W m-2.
126 - All other cubes are returned unchanged.
127 - The return type matches the input type (single Cube or CubeList).
129 Notes
130 -----
131 - The conversion uses a fixed latent heat of vaporisation:
132 LV = 2.5 × 10^6 J kg-1
133 - In reality, Lc varies with temperature (~5% variation between -20 °C
134 and +40 °C). This dependency is currently neglected but could be
135 included in future improvements.
136 - This function does not attempt to identify specific variables; it relies
137 solely on unit convertibility to determine applicability.
138 """
139 REQUIRED_UNITS = Unit("kg m-2 s-1")
140 OUTPUT_UNITS = Unit("W m-2")
142 out = iris.cube.CubeList()
143 for cube in iter_maybe(cubes):
144 # ACT ON MASS FLUXES
145 if cube.units is None or cube.units.is_unknown():
146 out.append(cube)
147 continue
148 if not cube.units.is_convertible(REQUIRED_UNITS):
149 # e.g. if UM LE or some other diagnostic — leave untouched
150 out.append(cube)
151 continue
153 cube_a = cube.copy()
154 cube_a = cube_a * LV
155 cube_a.units = cube.units * Unit("J kg-1")
156 cube_a.convert_units(OUTPUT_UNITS)
157 out.append(cube_a)
159 return out[0] if len(out) == 1 else out
162def net_radiative_flux(
163 downwelling_flux: iris.cube.Cube | iris.cube.CubeList,
164 upwelling_flux: iris.cube.Cube | iris.cube.CubeList,
165) -> iris.cube.Cube | iris.cube.CubeList:
166 """
167 Calculate net downward radiative flux.
169 Computes:
171 F_net = F_down - F_up
173 where positive values indicate a net downward flux.
174 Will work for a pair of downwelling and upwelling fluxes at any atmospheric level.
176 Parameters
177 ----------
178 downwelling_flux : iris.cube.Cube or iris.cube.CubeList
179 Downwelling radiative flux.
181 upwelling_flux : iris.cube.Cube or iris.cube.CubeList
182 Upwelling radiative flux.
184 Returns
185 -------
186 iris.cube.Cube or iris.cube.CubeList
187 Net downward radiative flux cube(s).
188 Returned cube name is set to either shortwave or longwave
189 (or neither) depending on input names (if specified).
190 """
191 out = iris.cube.CubeList()
193 for down_cube, up_cube in zip(
194 iter_maybe(downwelling_flux),
195 iter_maybe(upwelling_flux),
196 strict=True,
197 ):
198 down_name = down_cube.name().lower()
199 up_name = up_cube.name().lower()
200 if ("shortwave" in down_name and "longwave" in up_name) or (
201 "longwave" in down_name and "shortwave" in up_name
202 ):
203 raise ValueError("Cannot combine shortwave and longwave fluxes")
205 net_flux = down_cube - up_cube
206 name = down_cube.name().lower()
207 if "shortwave" in name:
208 output_name = "net_downward_shortwave_flux"
209 elif "longwave" in name: 209 ↛ 212line 209 didn't jump to line 212 because the condition on line 209 was always true
210 output_name = "net_downward_longwave_flux"
211 else:
212 output_name = "net_downward_radiative_flux"
214 net_flux.rename(output_name)
215 net_flux.units = down_cube.units
216 net_flux.var_name = output_name
218 out.append(net_flux)
220 return out[0] if len(out) == 1 else out