From d7fc75c15b896a5ffe9e22a03d7e31a22c84fc03 Mon Sep 17 00:00:00 2001 From: Israel Silber Date: Fri, 13 May 2022 10:38:38 -0400 Subject: [PATCH 1/8] More accurate e_s calculation following MK2005 --- src/metpy/calc/thermo.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index a85c4d6df4d..fc402836e97 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1035,14 +1035,19 @@ def saturation_vapor_pressure(temperature): Instead of temperature, dewpoint may be used in order to calculate the actual (ambient) water vapor (partial) pressure. - The formula used is that from [Bolton1980]_ for T in degrees Celsius: + The formula used is from eq. 10 in [MurphyKoop2005]_ for T in degrees Kelvin: - .. math:: 6.112 e^\frac{17.67T}{T + 243.5} + .. math:: e^[54.842763 - \frac{6763.22}{T} - 4.210\:ln(T) + 0.000367T + + tanh(0.0415(T - 218.8))(53.878 - \frac{1331.22}{T} + - 9.44523\:ln(T) + 0.014025T)] """ - # Converted from original in terms of C to use kelvin. - return mpconsts.nounit.sat_pressure_0c * np.exp( - 17.67 * (temperature - 273.15) / (temperature - 29.65) + # valid for 123 < T < 332 + return np.exp( + 54.842763 - 6763.22 / temperature - 4.210 * np.log(temperature) + + 0.000367 * temperature + np.tanh(0.0415 * (temperature - 218.8)) + * (53.878 - 1331.22 / temperature - 9.44523 * np.log(temperature) + + 0.014025 * temperature) ) From 36abca7e11914d343fd6cdc6a69417b0078dcd13 Mon Sep 17 00:00:00 2001 From: Israel Silber Date: Fri, 13 May 2022 10:47:04 -0400 Subject: [PATCH 2/8] ADD: reference to Murphy and Koop (2005) --- docs/api/references.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/api/references.rst b/docs/api/references.rst index 5c64c0603cb..9bfd417ade7 100644 --- a/docs/api/references.rst +++ b/docs/api/references.rst @@ -141,6 +141,10 @@ References *Journal of Geodesy* **74**, 128-133, doi:`10.1007/s001900050278 `_. +.. [MurphyKoop2005] Murphy, D.M. and T. Koop, 2005: Review of the vapour pressures of ice and + supercooled water for atmospheric applications. *Q.J.R. Meteorol. Soc.*, **131**, + 1539-1565, doi:`10.1256/qj.04.94 `_. + .. [NOAA1976] National Oceanic and Atmospheric Administration, National Aeronautics and Space Administration, and U. S. Air Force, 1976: `U. S. Standard Atmosphere 1976 `_, From d12c224ce585cb2a6add9728c780eabb92557440 Mon Sep 17 00:00:00 2001 From: Israel Silber Date: Fri, 13 May 2022 11:25:36 -0400 Subject: [PATCH 3/8] ADD: method to calculate ice equilibrium pressure of vapor --- src/metpy/calc/thermo.py | 34 +++++++++++++++++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index fc402836e97..fda580eee5a 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1014,7 +1014,7 @@ def vapor_pressure(pressure, mixing_ratio): @preprocess_and_wrap(wrap_like='temperature') @process_units({'temperature': '[temperature]'}, '[pressure]') def saturation_vapor_pressure(temperature): - r"""Calculate the saturation water vapor (partial) pressure. + r"""Calculate the saturation (equilibrium) water vapor (partial) pressure. Parameters ---------- @@ -1051,6 +1051,38 @@ def saturation_vapor_pressure(temperature): ) +@exporter.export +@preprocess_and_wrap(wrap_like='temperature') +@process_units({'temperature': '[temperature]'}, '[pressure]') +def ice_saturation_vapor_pressure(temperature): + r"""Calculate the ice saturation (equilibrium) water vapor (partial) pressure. + + Parameters + ---------- + temperature : `pint.Quantity` + Air temperature + + Returns + ------- + `pint.Quantity` + Ice saturation water vapor (partial) pressure + + See Also + -------- + vapor_pressure + + The formula used is from eq. 7 in [MurphyKoop2005]_ for T in degrees Kelvin: + + .. math:: e^[9.550426 - \frac{5723.265}{T} + 3.53068\:ln(T) - 0.00728332T] + + """ + # valid for T > 110 K + return np.exp( + 9.550426 - 5723.265 / temperature + 3.53068 * np.log(temperature) + - 0.00728332 * temperature + ) + + @exporter.export @preprocess_and_wrap(wrap_like='temperature', broadcast=('temperature', 'relative_humidity')) @check_units('[temperature]', '[dimensionless]') From 5b8ffdc236a975940dca7eee397b343d3da7cfac Mon Sep 17 00:00:00 2001 From: Israel Silber Date: Fri, 13 May 2022 11:51:31 -0400 Subject: [PATCH 4/8] pep8 for thermo --- src/metpy/calc/thermo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index fda580eee5a..909cff99b8f 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1047,7 +1047,7 @@ def saturation_vapor_pressure(temperature): 54.842763 - 6763.22 / temperature - 4.210 * np.log(temperature) + 0.000367 * temperature + np.tanh(0.0415 * (temperature - 218.8)) * (53.878 - 1331.22 / temperature - 9.44523 * np.log(temperature) - + 0.014025 * temperature) + + 0.014025 * temperature) ) From 1e65c9a44a94cd3ec314eb3450c727ee008e0ff8 Mon Sep 17 00:00:00 2001 From: Israel Silber Date: Fri, 13 May 2022 12:07:43 -0400 Subject: [PATCH 5/8] ADD ice_saturation_vapor_pressure to overrides --- docs/_templates/overrides/metpy.calc.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst index f6eeeaf3724..99fd4c28eaa 100644 --- a/docs/_templates/overrides/metpy.calc.rst +++ b/docs/_templates/overrides/metpy.calc.rst @@ -35,6 +35,7 @@ Moist Thermodynamics dewpoint_from_relative_humidity dewpoint_from_specific_humidity equivalent_potential_temperature + ice_saturation_vapor_pressure mixing_ratio mixing_ratio_from_relative_humidity mixing_ratio_from_specific_humidity From 151e0fe43296dec351c40f2ed3cbd62e8a689ec0 Mon Sep 17 00:00:00 2001 From: Israel Silber Date: Fri, 13 May 2022 17:56:18 -0400 Subject: [PATCH 6/8] ADD: unit test for ice_saturation_vapor_pressure --- tests/calc/test_thermo.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index cc8e5827bfc..6a853aa2691 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -25,7 +25,7 @@ relative_humidity_from_specific_humidity, relative_humidity_wet_psychrometric, saturation_equivalent_potential_temperature, saturation_mixing_ratio, - saturation_vapor_pressure, showalter_index, + saturation_vapor_pressure, ice_saturation_vapor_pressure, showalter_index, specific_humidity_from_dewpoint, specific_humidity_from_mixing_ratio, static_stability, surface_based_cape_cin, temperature_from_potential_temperature, thickness_hydrostatic, @@ -295,6 +295,16 @@ def test_sat_vapor_pressure(): assert_array_almost_equal(saturation_vapor_pressure(temp), real_es, 2) +def test_ice_sat_vapor_pressure(): + """Test ice_saturation_vapor_pressure calculation.""" + temp = (np.arange(1, 38) * units.degC).to(units.degK) + assert np.all(calcs.saturation_vapor_pressure(temp) + < calcs.ice_saturation_vapor_pressure(temp)) + temp = (temp.to(units.degC) * (-1)).to(units.degK) + assert np.all(calcs.saturation_vapor_pressure(temp) + > calcs.ice_saturation_vapor_pressure(temp)) + + def test_sat_vapor_pressure_scalar(): """Test saturation_vapor_pressure handles scalar values.""" es = saturation_vapor_pressure(0 * units.degC) From 240269406e80f9b3c987654eb1396ed2d128c4e8 Mon Sep 17 00:00:00 2001 From: Israel Silber Date: Fri, 13 May 2022 18:00:40 -0400 Subject: [PATCH 7/8] FIX: typo in test_ice_saturation_vapor_pressure --- tests/calc/test_thermo.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index 6a853aa2691..2b86c487841 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -298,11 +298,11 @@ def test_sat_vapor_pressure(): def test_ice_sat_vapor_pressure(): """Test ice_saturation_vapor_pressure calculation.""" temp = (np.arange(1, 38) * units.degC).to(units.degK) - assert np.all(calcs.saturation_vapor_pressure(temp) - < calcs.ice_saturation_vapor_pressure(temp)) + assert np.all(saturation_vapor_pressure(temp) + < ice_saturation_vapor_pressure(temp)) temp = (temp.to(units.degC) * (-1)).to(units.degK) - assert np.all(calcs.saturation_vapor_pressure(temp) - > calcs.ice_saturation_vapor_pressure(temp)) + assert np.all(saturation_vapor_pressure(temp) + > ice_saturation_vapor_pressure(temp)) def test_sat_vapor_pressure_scalar(): From 9c9e6e7b8a5a38ed40f0fff6ed417c236b91d819 Mon Sep 17 00:00:00 2001 From: Israel Silber Date: Fri, 13 May 2022 18:13:47 -0400 Subject: [PATCH 8/8] pep8 --- tests/calc/test_thermo.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index 2b86c487841..c71f7df059d 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -25,7 +25,8 @@ relative_humidity_from_specific_humidity, relative_humidity_wet_psychrometric, saturation_equivalent_potential_temperature, saturation_mixing_ratio, - saturation_vapor_pressure, ice_saturation_vapor_pressure, showalter_index, + saturation_vapor_pressure, ice_saturation_vapor_pressure, + showalter_index, specific_humidity_from_dewpoint, specific_humidity_from_mixing_ratio, static_stability, surface_based_cape_cin, temperature_from_potential_temperature, thickness_hydrostatic,