From dea48047c40ef07b3238eba2621a65221a3f900a Mon Sep 17 00:00:00 2001 From: Robert Neal Date: Wed, 23 Jul 2025 11:33:49 +0100 Subject: [PATCH 1/4] Code changes to add the SVP derivative corrected for in air --- .../psychrometric_calculations.py | 77 +++++++++++++++++++ .../test_calculate_svp_derivative_in_air.py | 50 ++++++++++++ 2 files changed, 127 insertions(+) create mode 100644 improver_tests/psychrometric_calculations/test_calculate_svp_derivative_in_air.py diff --git a/improver/psychrometric_calculations/psychrometric_calculations.py b/improver/psychrometric_calculations/psychrometric_calculations.py index f4177cef63..fdca79cfad 100644 --- a/improver/psychrometric_calculations/psychrometric_calculations.py +++ b/improver/psychrometric_calculations/psychrometric_calculations.py @@ -18,6 +18,7 @@ from improver import BasePlugin from improver.generate_ancillaries.generate_svp_table import ( SaturatedVapourPressureTable, + SaturatedVapourPressureTableDerivative, ) from improver.metadata.utilities import ( create_new_diagnostic_cube, @@ -57,6 +58,26 @@ def _svp_table() -> ndarray: return svp_data.data +@functools.lru_cache() +def _svp_derivative_table() -> ndarray: + """ + Calculate a saturated vapour pressure (SVP) derivative lookup table. + The lru_cache decorator caches this table on first call to this function, + so that the table does not need to be re-calculated if used multiple times. + + A value of SVP derivative for any temperature between T_MIN and T_MAX (inclusive) can be + obtained by interpolating through the table, as is done in the _svp_derivative_from_lookup + function. + + Returns: + Array of first derivative saturated vapour pressures (Pa). + """ + svp_derivative_data = SaturatedVapourPressureTableDerivative( + t_min=SVP_T_MIN, t_max=SVP_T_MAX, t_increment=SVP_T_INCREMENT + ).process() + return svp_derivative_data.data + + def _svp_from_lookup(temperature: ndarray) -> ndarray: """ Gets value for saturation vapour pressure in a pure water vapour system @@ -84,6 +105,33 @@ def _svp_from_lookup(temperature: ndarray) -> ndarray: ] + interpolation_factor * svp_table_data[table_index + 1] +def _svp_derivative_from_lookup(temperature: ndarray) -> ndarray: + """ + Gets value for saturation vapour pressure derivative in a pure water vapour system + from a pre-calculated lookup table. Interpolates linearly between points in + the table to the temperatures required. + + Args: + temperature: + Array of air temperatures (K). + + Returns: + Array of first derivative saturated vapour pressures (Pa). + """ + # where temperatures are outside the SVP derivative table range, clip data to + # within the available range + t_clipped = np.clip(temperature, SVP_T_MIN, SVP_T_MAX - SVP_T_INCREMENT) + + # interpolate between bracketing values + table_position = (t_clipped - SVP_T_MIN) / SVP_T_INCREMENT + table_index = table_position.astype(int) + interpolation_factor = table_position - table_index + svp_derivative_table_data = _svp_derivative_table() + return (1.0 - interpolation_factor) * svp_derivative_table_data[ + table_index + ] + interpolation_factor * svp_derivative_table_data[table_index + 1] + + def calculate_svp_in_air(temperature: ndarray, pressure: ndarray) -> ndarray: """ Calculates the saturation vapour pressure in air. Looks up the saturation @@ -109,6 +157,35 @@ def calculate_svp_in_air(temperature: ndarray, pressure: ndarray) -> ndarray: return svp * correction.astype(np.float32) +def calculate_svp_derivative_in_air(temperature: ndarray, pressure: ndarray) -> ndarray: + """ + Calculates the saturation vapour pressure derivative in air. Looks up the saturation + vapour pressure derivative in a pure water vapour system, and pressure-corrects the + result to obtain the saturation vapour pressure derivative in air. + + Args: + temperature: + Array of air temperatures (K). + pressure: + Array of pressure (Pa). + + Returns: + Saturation vapour pressure derivative in air (Pa). + + References: + Atmosphere-Ocean Dynamics, Adrian E. Gill, International Geophysics + Series, Vol. 30; Equation A4.7. + """ + svp = _svp_from_lookup(temperature) + svp_derivative = _svp_derivative_from_lookup(temperature) + temp_C = temperature + consts.ABSOLUTE_ZERO + correction = 1.0 + 1.0e-8 * pressure * (4.5 + 6.0e-4 * temp_C * temp_C) + derivative_correction_term = correction * svp_derivative + ( + 2 * 1.0e-8 * 6.0e-4 * pressure * svp + ) + return svp_derivative * derivative_correction_term.astype(np.float32) + + def dry_adiabatic_temperature( initial_temperature: ndarray, initial_pressure: ndarray, final_pressure: ndarray ) -> ndarray: diff --git a/improver_tests/psychrometric_calculations/test_calculate_svp_derivative_in_air.py b/improver_tests/psychrometric_calculations/test_calculate_svp_derivative_in_air.py new file mode 100644 index 0000000000..5a7bf29ea6 --- /dev/null +++ b/improver_tests/psychrometric_calculations/test_calculate_svp_derivative_in_air.py @@ -0,0 +1,50 @@ +# (C) Crown Copyright, Met Office. All rights reserved. +# +# This file is part of 'IMPROVER' and is released under the BSD 3-Clause license. +# See LICENSE in the root of the repository for full licensing details. +"""Unit tests for psychrometric_calculations calculate_svp_derivative_in_air""" + +import unittest + +import numpy as np +from iris.tests import IrisTest + +from improver.psychrometric_calculations.psychrometric_calculations import ( + _svp_derivative_from_lookup, + calculate_svp_derivative_in_air, +) + + +class Test_calculate_svp_derivative_in_air(IrisTest): + """Test the calculate_svp_derivative_in_air function""" + + def setUp(self): + """Set up test data""" + self.temperature = np.array([[185.0, 260.65, 338.15]], dtype=np.float32) + self.pressure = np.array([[1.0e5, 9.9e4, 9.8e4]], dtype=np.float32) + + def test_calculate_svp_derivative_in_air(self): + """Test pressure-corrected SVP derivative values""" + expected = np.array([[0.01362905, 208.47170252, 25187.76423485]]) + result = calculate_svp_derivative_in_air(self.temperature, self.pressure) + np.testing.assert_allclose(result, expected, rtol=1e-5, atol=1e-5) + + def test_values(self): + """Basic extraction of SVP derivative values from lookup table""" + self.temperature[0, 1] = 260.56833 + expected = [[1.350531e-02, 2.06000274e02, 2.501530e04]] + result = _svp_derivative_from_lookup(self.temperature) + np.testing.assert_allclose(result, expected, rtol=1e-5, atol=1e-5) + + def test_beyond_table_bounds(self): + """Extracting SVP derivative values from the table with temperatures beyond + its valid range. Should return the nearest end of the table.""" + self.temperature[0, 0] = 150.0 + self.temperature[0, 2] = 400.0 + expected = [[9.664590e-03, 2.075279e02, 2.501530e04]] + result = _svp_derivative_from_lookup(self.temperature) + np.testing.assert_allclose(result, expected, rtol=1e-5, atol=1e-5) + + +if __name__ == "__main__": + unittest.main() From 583c96c654dc2802a9e26f988df853831f1120c8 Mon Sep 17 00:00:00 2001 From: Robert Neal Date: Fri, 25 Jul 2025 10:14:17 +0100 Subject: [PATCH 2/4] Adding a missing temperature term to the correction calculation --- .../psychrometric_calculations/psychrometric_calculations.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/improver/psychrometric_calculations/psychrometric_calculations.py b/improver/psychrometric_calculations/psychrometric_calculations.py index fdca79cfad..c617eb654d 100644 --- a/improver/psychrometric_calculations/psychrometric_calculations.py +++ b/improver/psychrometric_calculations/psychrometric_calculations.py @@ -181,7 +181,7 @@ def calculate_svp_derivative_in_air(temperature: ndarray, pressure: ndarray) -> temp_C = temperature + consts.ABSOLUTE_ZERO correction = 1.0 + 1.0e-8 * pressure * (4.5 + 6.0e-4 * temp_C * temp_C) derivative_correction_term = correction * svp_derivative + ( - 2 * 1.0e-8 * 6.0e-4 * pressure * svp + 2 * 1.0e-8 * 6.0e-4 * pressure * temp_C * svp ) return svp_derivative * derivative_correction_term.astype(np.float32) From c6ddfaf241fce677c4d5b57190d65b11a7615495 Mon Sep 17 00:00:00 2001 From: Robert Neal Date: Fri, 25 Jul 2025 12:57:08 +0100 Subject: [PATCH 3/4] Updating expected values in unit tests and changing one of the imports so it correctly finds the function --- .../psychrometric_calculations.py | 4 +++- .../test_calculate_svp_derivative_in_air.py | 6 +++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/improver/psychrometric_calculations/psychrometric_calculations.py b/improver/psychrometric_calculations/psychrometric_calculations.py index c617eb654d..2fe96d6063 100644 --- a/improver/psychrometric_calculations/psychrometric_calculations.py +++ b/improver/psychrometric_calculations/psychrometric_calculations.py @@ -16,9 +16,11 @@ import improver.constants as consts from improver import BasePlugin +from improver.generate_ancillaries.generate_svp_derivative_table import ( + SaturatedVapourPressureTableDerivative, +) from improver.generate_ancillaries.generate_svp_table import ( SaturatedVapourPressureTable, - SaturatedVapourPressureTableDerivative, ) from improver.metadata.utilities import ( create_new_diagnostic_cube, diff --git a/improver_tests/psychrometric_calculations/test_calculate_svp_derivative_in_air.py b/improver_tests/psychrometric_calculations/test_calculate_svp_derivative_in_air.py index 5a7bf29ea6..f891af7856 100644 --- a/improver_tests/psychrometric_calculations/test_calculate_svp_derivative_in_air.py +++ b/improver_tests/psychrometric_calculations/test_calculate_svp_derivative_in_air.py @@ -25,14 +25,14 @@ def setUp(self): def test_calculate_svp_derivative_in_air(self): """Test pressure-corrected SVP derivative values""" - expected = np.array([[0.01362905, 208.47170252, 25187.76423485]]) + expected = np.array([[5.89845229e-06, 3.54367486e02, 1.26270031e06]]) result = calculate_svp_derivative_in_air(self.temperature, self.pressure) np.testing.assert_allclose(result, expected, rtol=1e-5, atol=1e-5) def test_values(self): """Basic extraction of SVP derivative values from lookup table""" self.temperature[0, 1] = 260.56833 - expected = [[1.350531e-02, 2.06000274e02, 2.501530e04]] + expected = [[2.41833058e-03, 1.86569996e01, 1.11889656e03]] result = _svp_derivative_from_lookup(self.temperature) np.testing.assert_allclose(result, expected, rtol=1e-5, atol=1e-5) @@ -41,7 +41,7 @@ def test_beyond_table_bounds(self): its valid range. Should return the nearest end of the table.""" self.temperature[0, 0] = 150.0 self.temperature[0, 2] = 400.0 - expected = [[9.664590e-03, 2.075279e02, 2.501530e04]] + expected = [[1.76528667e-03, 1.86569996e01, 1.11889656e03]] result = _svp_derivative_from_lookup(self.temperature) np.testing.assert_allclose(result, expected, rtol=1e-5, atol=1e-5) From 47d90205234d4bf9ecb86caa7100c767aa1b4a00 Mon Sep 17 00:00:00 2001 From: Robert Neal Date: Fri, 25 Jul 2025 13:50:16 +0100 Subject: [PATCH 4/4] Updating one of the ecpected values within the unit tests --- .../test_calculate_svp_derivative_in_air.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/improver_tests/psychrometric_calculations/test_calculate_svp_derivative_in_air.py b/improver_tests/psychrometric_calculations/test_calculate_svp_derivative_in_air.py index f891af7856..b200846232 100644 --- a/improver_tests/psychrometric_calculations/test_calculate_svp_derivative_in_air.py +++ b/improver_tests/psychrometric_calculations/test_calculate_svp_derivative_in_air.py @@ -41,7 +41,7 @@ def test_beyond_table_bounds(self): its valid range. Should return the nearest end of the table.""" self.temperature[0, 0] = 150.0 self.temperature[0, 2] = 400.0 - expected = [[1.76528667e-03, 1.86569996e01, 1.11889656e03]] + expected = [[1.76528667e-03, 1.87835245e01, 1.11889656e03]] result = _svp_derivative_from_lookup(self.temperature) np.testing.assert_allclose(result, expected, rtol=1e-5, atol=1e-5)