EPPT-2510: Adding calculation of saturation vapour pressure derivative in air

improver/psychrometric_calculations/psychrometric_calculations.py

@@ -16,6 +16,9 @@
 
 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,
 )
@@ -57,6 +60,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 +107,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 +159,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 * temp_C * svp
+    )
+    return svp_derivative * derivative_correction_term.astype(np.float32)
+
+
 def dry_adiabatic_temperature(
     initial_temperature: ndarray, initial_pressure: ndarray, final_pressure: ndarray
 ) -> ndarray:

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([[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 = [[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)
+
+    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 = [[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)
+
+
+if __name__ == "__main__":
+    unittest.main()