From a5a67e9b8b71db6b81fe8a36a8328332f93b00b8 Mon Sep 17 00:00:00 2001
From: bjorn-stevens <64255981+bjorn-stevens@users.noreply.github.com>
Date: Tue, 16 Aug 2022 15:14:58 +0200
Subject: [PATCH] v0.4: includes reformatting of svp analytic
Released as version 0.4, includes reformatting of analytic formulation
of saturation vapor pressure to better emphasize the commont form of the
calculation for ice and liquid, similar to what is done for teten's
formulas. Now both it and Teten's serve as 'master' formlae, with liq
and ice being particular implementations.
---
examples/examples.ipynb | 165 +++++++++++++-----
.../saturation_vapor_pressures.py | 85 +++++----
setup.py | 2 +-
3 files changed, 171 insertions(+), 81 deletions(-)
diff --git a/examples/examples.ipynb b/examples/examples.ipynb
index 109f024..d019579 100644
--- a/examples/examples.ipynb
+++ b/examples/examples.ipynb
@@ -574,30 +574,37 @@
"source": [
"from scipy.optimize import curve_fit\n",
"\n",
- "def liq_error(T,a,b):\n",
- " return np.abs(svp.tetens(T,a,b)/svp.liq_wagner_pruss(T) -1.)\n",
"\n",
- "def ice_error(T,a,b):\n",
- " return np.abs(svp.tetens(T,a,b)/svp.ice_wagner_etal(T) -1.)\n",
+ "def liq_error(T, a, b):\n",
+ " return np.abs(svp.tetens(T, a, b) / svp.liq_wagner_pruss(T) - 1.0)\n",
"\n",
- "T = np.arange(270.,310.,0.1)\n",
+ "\n",
+ "def ice_error(T, a, b):\n",
+ " return np.abs(svp.tetens(T, a, b) / svp.ice_wagner_etal(T) - 1.0)\n",
+ "\n",
+ "\n",
+ "T = np.arange(270.0, 310.0, 0.1)\n",
"\n",
"rng = np.random.default_rng()\n",
"y_noise = 0.001 * rng.normal(size=T.size)\n",
- "ydata = y_noise\n",
- "popt, pcov = curve_fit(liq_error, T, ydata, bounds = ((16.,33.), (19.,36.)), method='dogbox')\n",
+ "ydata = y_noise\n",
+ "popt, pcov = curve_fit(\n",
+ " liq_error, T, ydata, bounds=((16.0, 33.0), (19.0, 36.0)), method=\"dogbox\"\n",
+ ")\n",
"a_liq = popt[0]\n",
"b_liq = popt[1]\n",
- "print (f'Best fit parameters for liquid: a={a_liq:.4f}, b={b_liq:.4f}')\n",
+ "print(f\"Best fit parameters for liquid: a={a_liq:.4f}, b={b_liq:.4f}\")\n",
"\n",
- "T = np.arange(230.,260.,0.01)\n",
+ "T = np.arange(230.0, 260.0, 0.01)\n",
"rng = np.random.default_rng()\n",
"y_noise = 0.001 * rng.normal(size=T.size)\n",
"ydata = y_noise\n",
- "popt, pcov = curve_fit(ice_error, T, ydata, bounds = ((20.,5.), (23.,8.)), method='dogbox')\n",
+ "popt, pcov = curve_fit(\n",
+ " ice_error, T, ydata, bounds=((20.0, 5.0), (23.0, 8.0)), method=\"dogbox\"\n",
+ ")\n",
"a_ice = popt[0]\n",
"b_ice = popt[1]\n",
- "print (f'Best fit parameters for ice: a={a_ice:.4f}, b={b_ice:.4f}')"
+ "print(f\"Best fit parameters for ice: a={a_ice:.4f}, b={b_ice:.4f}\")"
]
},
{
@@ -621,54 +628,116 @@
],
"source": [
"sns.set_context(\"talk\")\n",
- "fig, ax = plt.subplots(1,3,figsize=(12, 3), constrained_layout=True)\n",
- "\n",
- "T = np.arange(200.,273.16)\n",
- "es_c1 = (svp.ice_analytic(T)/svp.ice_wagner_etal(T) - 1.)*100.\n",
- "es_c2 = (svp.ice_analytic(T,ci=1861.)/svp.ice_wagner_etal(T) - 1.)*100.\n",
- "es_c3 = (svp.tetens(T,21.875,7.66)/svp.ice_wagner_etal(T) - 1.)*100.\n",
- "es_c4 = (svp.tetens(T,a_ice,b_ice)/svp.ice_wagner_etal(T) - 1.)*100.\n",
- "\n",
- "ax[0].plot(T,np.abs(es_c1),c='k',label='Romps')\n",
- "ax[0].plot(T,np.abs(es_c2),c='k',ls='dotted',label='Romps best fit')\n",
- "ax[0].plot(T,np.abs(es_c3),c='violet',label='Teten-Murray')\n",
- "ax[0].plot(T,np.abs(es_c4),c='violet',ls='dotted',label='Teten best fit')\n",
- "T = np.arange(235.,273.16)\n",
- "es_sc1 = (svp.liq_analytic(T)/svp.liq_murphy_koop(T) - 1.)*100.\n",
- "es_sc2 = (svp.liq_analytic(T,cl=4119.)/svp.liq_murphy_koop(T) - 1.)*100.\n",
- "es_sc3 = (svp.tetens(T,17.269,35.86)/svp.liq_murphy_koop(T) - 1.)*100.\n",
- "es_sc4 = (svp.tetens(T,a_liq,b_liq)/svp.liq_murphy_koop(T) - 1.)*100.\n",
- "\n",
- "ax[1].plot(T,np.abs(es_sc1),c='k')\n",
- "ax[1].plot(T,np.abs(es_sc2),c='k',ls='dotted')\n",
- "ax[1].plot(T,np.abs(es_sc3),c='violet')\n",
- "ax[1].plot(T,np.abs(es_sc4),c='violet',ls='dotted')\n",
- "\n",
- "T = np.arange(273.,330.)\n",
- "es_w1 = (svp.liq_analytic(T)/svp.liq_wagner_pruss(T) - 1.)*100.\n",
- "es_w2 = (svp.liq_analytic(T,cl=4119.)/svp.liq_wagner_pruss(T) - 1.)*100.\n",
- "es_w3 = (svp.tetens(T,17.269,35.86)/svp.liq_wagner_pruss(T) - 1.)*100.\n",
- "es_w4 = (svp.tetens(T,a_liq,b_liq)/svp.liq_wagner_pruss(T) - 1.)*100.\n",
- "\n",
- "ax[2].plot(T,np.abs(es_w1),c='k',label='Romps')\n",
- "ax[2].plot(T,np.abs(es_w2),c='k',ls='dotted',label='Romps best fit')\n",
- "ax[2].plot(T,np.abs(es_w3),c='violet',label='Teten-Murray')\n",
- "ax[2].plot(T,np.abs(es_w4),c='violet',ls='dotted',label='Teten best fit')\n",
+ "fig, ax = plt.subplots(1, 3, figsize=(12, 3), constrained_layout=True)\n",
+ "\n",
+ "T = np.arange(200.0, 273.16)\n",
+ "es_c1 = (svp.ice_analytic(T) / svp.ice_wagner_etal(T) - 1.0) * 100.0\n",
+ "es_c2 = (svp.ice_analytic(T, ci=1861.0) / svp.ice_wagner_etal(T) - 1.0) * 100.0\n",
+ "es_c3 = (svp.tetens(T, 21.875, 7.66) / svp.ice_wagner_etal(T) - 1.0) * 100.0\n",
+ "es_c4 = (svp.tetens(T, a_ice, b_ice) / svp.ice_wagner_etal(T) - 1.0) * 100.0\n",
+ "\n",
+ "ax[0].plot(T, np.abs(es_c1), c=\"k\", label=\"Romps\")\n",
+ "ax[0].plot(T, np.abs(es_c2), c=\"k\", ls=\"dotted\", label=\"Romps best fit\")\n",
+ "ax[0].plot(T, np.abs(es_c3), c=\"violet\", label=\"Teten-Murray\")\n",
+ "ax[0].plot(T, np.abs(es_c4), c=\"violet\", ls=\"dotted\", label=\"Teten best fit\")\n",
+ "T = np.arange(235.0, 273.16)\n",
+ "es_sc1 = (svp.liq_analytic(T) / svp.liq_murphy_koop(T) - 1.0) * 100.0\n",
+ "es_sc2 = (svp.liq_analytic(T, cl=4119.0) / svp.liq_murphy_koop(T) - 1.0) * 100.0\n",
+ "es_sc3 = (svp.tetens(T, 17.269, 35.86) / svp.liq_murphy_koop(T) - 1.0) * 100.0\n",
+ "es_sc4 = (svp.tetens(T, a_liq, b_liq) / svp.liq_murphy_koop(T) - 1.0) * 100.0\n",
+ "\n",
+ "ax[1].plot(T, np.abs(es_sc1), c=\"k\")\n",
+ "ax[1].plot(T, np.abs(es_sc2), c=\"k\", ls=\"dotted\")\n",
+ "ax[1].plot(T, np.abs(es_sc3), c=\"violet\")\n",
+ "ax[1].plot(T, np.abs(es_sc4), c=\"violet\", ls=\"dotted\")\n",
+ "\n",
+ "T = np.arange(273.0, 330.0)\n",
+ "es_w1 = (svp.liq_analytic(T) / svp.liq_wagner_pruss(T) - 1.0) * 100.0\n",
+ "es_w2 = (svp.liq_analytic(T, cl=4119.0) / svp.liq_wagner_pruss(T) - 1.0) * 100.0\n",
+ "es_w3 = (svp.tetens(T, 17.269, 35.86) / svp.liq_wagner_pruss(T) - 1.0) * 100.0\n",
+ "es_w4 = (svp.tetens(T, a_liq, b_liq) / svp.liq_wagner_pruss(T) - 1.0) * 100.0\n",
+ "\n",
+ "ax[2].plot(T, np.abs(es_w1), c=\"k\", label=\"Romps\")\n",
+ "ax[2].plot(T, np.abs(es_w2), c=\"k\", ls=\"dotted\", label=\"Romps best fit\")\n",
+ "ax[2].plot(T, np.abs(es_w3), c=\"violet\", label=\"Teten-Murray\")\n",
+ "ax[2].plot(T, np.abs(es_w4), c=\"violet\", ls=\"dotted\", label=\"Teten best fit\")\n",
"\n",
"ax[0].legend()\n",
"for a in ax:\n",
- " a.set_ylabel('abs error / %')\n",
- " a.set_xlabel('T / K')\n",
+ " a.set_ylabel(\"abs error / %\")\n",
+ " a.set_xlabel(\"T / K\")\n",
"\n",
- "sns.despine (offset=10)"
+ "sns.despine(offset=10)"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 36,
"id": "fef2f42d-8685-4cd8-9a4c-e0014f721f8c",
"metadata": {},
"outputs": [],
+ "source": [
+ "def analytic(T, lx, cx):\n",
+ " \"\"\"returns saturation vapor pressure over a given phase\n",
+ "\n",
+ " Uses the rankine (constant specific heat, negligible condensate volume) approximations to\n",
+ " calculate the saturation vapor pressure over a phase with the specific heat cx, and phase\n",
+ " change enthalpy (from vapor) lx, at temperature T.\n",
+ "\n",
+ " Args:\n",
+ " T: temperature in kelvin\n",
+ " lx: phase change enthalpy between vapor and given phase (liquid, ice)\n",
+ " cx: specific heat capacity of given phase (liquid, ice)\n",
+ "\n",
+ " Returns:\n",
+ " value of saturation vapor pressure over liquid water in Pa\n",
+ "\n",
+ " Reference:\n",
+ " Romps, D. M. Exact Expression for the Lifting Condensation Level. Journal of the Atmospheric\n",
+ " Sciences 74, 3891–3900 (2017).\n",
+ " Romps, D. M. Accurate expressions for the dew point and frost point derived from the Rankine-\n",
+ " Kirchhoff approximations. Journal of the Atmospheric Sciences (2021) doi:10.1175/JAS-D-20-0301.1.\n",
+ "\n",
+ " >>> analytic(np.asarray([273.16,305.]))\n",
+ " array([ 611.655 , 4711.13161169])\n",
+ " \"\"\"\n",
+ " TvT = constants.temperature_water_vapor_triple_point\n",
+ " PvT = constants.pressure_water_vapor_triple_point\n",
+ " Rv = constants.water_vapor_gas_constant\n",
+ "\n",
+ " c1 = (constants.cpv - cx) / Rv\n",
+ " c2 = lx / (Rv * TvT) - c1\n",
+ " es = PvT * np.exp(c2 * (1.0 - TvT / T)) * (T / TvT) ** c1\n",
+ " return es"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 37,
+ "id": "93db947c-bb33-4ea1-8d8a-b8edf9f6aec1",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "4711.131611687174"
+ ]
+ },
+ "execution_count": 37,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "analytic(305.0, constants.lvT, constants.cl)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "86a706f5-ec56-44ca-97ce-5190cd27843c",
+ "metadata": {},
+ "outputs": [],
"source": []
}
],
diff --git a/moist_thermodynamics/saturation_vapor_pressures.py b/moist_thermodynamics/saturation_vapor_pressures.py
index 9893ddc..c89710a 100644
--- a/moist_thermodynamics/saturation_vapor_pressures.py
+++ b/moist_thermodynamics/saturation_vapor_pressures.py
@@ -142,7 +142,41 @@ def liq_hardy(T):
return np.exp(X)
-def liq_analytic(T, cl=constants.cl):
+def analytic(T, lx, cx):
+ """returns saturation vapor pressure over a given phase
+
+ Uses the rankine (constant specific heat, negligible condensate volume) approximations to
+ calculate the saturation vapor pressure over a phase with the specific heat cx, and phase
+ change enthalpy (from vapor) lx, at temperature T.
+
+ Args:
+ T: temperature in kelvin
+ lx: phase change enthalpy between vapor and given phase (liquid, ice)
+ cx: specific heat capacity of given phase (liquid, ice)
+
+ Returns:
+ value of saturation vapor pressure over liquid water in Pa
+
+ Reference:
+ Romps, D. M. Exact Expression for the Lifting Condensation Level. Journal of the Atmospheric
+ Sciences 74, 3891–3900 (2017).
+ Romps, D. M. Accurate expressions for the dew point and frost point derived from the Rankine-
+ Kirchhoff approximations. Journal of the Atmospheric Sciences (2021) doi:10.1175/JAS-D-20-0301.1.
+
+ >>> analytic(305.,constants.lvT,constants.cl)
+ 4711.131611687174
+ """
+ TvT = constants.temperature_water_vapor_triple_point
+ PvT = constants.pressure_water_vapor_triple_point
+ Rv = constants.water_vapor_gas_constant
+
+ c1 = (constants.cpv - cx) / Rv
+ c2 = lx / (Rv * TvT) - c1
+ es = PvT * np.exp(c2 * (1.0 - TvT / T)) * (T / TvT) ** c1
+ return es
+
+
+def liq_analytic(T, lx=constants.lvT, cx=constants.cl):
"""Analytic approximation for saturation vapor pressure over iquid
Uses the rankine (constant specific heat, negligible condensate volume) approximations to
@@ -153,7 +187,8 @@ def liq_analytic(T, cl=constants.cl):
Args:
T: temperature in kelvin
- delta_cl: differnce between isobaric specific heat capacity of vapor and that of liquid.
+ lx: enthalpy of vaporization, at triple point, default constants.lvT
+ cl: specific heat capacity of liquid at triple point
Returns:
value of saturation vapor pressure over liquid water in Pa
@@ -167,18 +202,10 @@ def liq_analytic(T, cl=constants.cl):
>>> liq_analytic(np.asarray([273.16,305.]))
array([ 611.655 , 4711.13161169])
"""
- TvT = constants.temperature_water_vapor_triple_point
- PvT = constants.pressure_water_vapor_triple_point
- lvT = constants.vaporization_enthalpy_triple_point
- Rv = constants.water_vapor_gas_constant
-
- c1 = (constants.cpv - cl) / Rv
- c2 = lvT / (Rv * TvT) - c1
- es = PvT * np.exp(c2 * (1.0 - TvT / T)) * (T / TvT) ** c1
- return es
+ return analytic(T, lx, cx)
-def ice_analytic(T, ci=constants.ci):
+def ice_analytic(T, lx=constants.lsT, cx=constants.ci):
"""Analytic approximation for saturation vapor pressure over ice
Uses the rankine (constant specific heat, negligible condensate volume) approximations to
@@ -189,10 +216,11 @@ def ice_analytic(T, ci=constants.ci):
Args:
T: temperature in kelvin
- delta_cl: differnce between isobaric specific heat capacity of vapor and that of liquid.
+ lx: enthalpy of sublimation, at triple point, default constants.lsT
+ ci: specific heat capacity of ice Ih, at triple point
Returns:
- value of saturation vapor pressure over liquid water in Pa
+ value of saturation vapor pressure over ice in Pa
Reference:
Romps, D. M. Exact Expression for the Lifting Condensation Level. Journal of the Atmospheric
@@ -204,23 +232,16 @@ def ice_analytic(T, ci=constants.ci):
>>> ice_analytic(np.asarray([273.16,260.]))
array([611.655 , 195.99959431])
"""
- TvT = constants.temperature_water_vapor_triple_point
- PvT = constants.pressure_water_vapor_triple_point
- lsT = constants.sublimation_enthalpy_triple_point
- Rv = constants.water_vapor_gas_constant
+ return analytic(T, lx, cx)
- c1 = (constants.cpv-ci) / Rv
- c2 = lsT / (Rv * TvT) - c1
- es = PvT * np.exp(c2 * (1.0 - TvT / T)) * (T / TvT) ** c1
- return es
-def tetens(T,a,b):
+def tetens(T, a, b):
"""Returns saturation vapor pressure over liquid using the Magnus-Teten's formula
- This equation is written in a general form, with the constants a and b determining the fit. As
- such it can be specified for either ice or water, or adapted as originally impelemented in ICON,
+ This equation is written in a general form, with the constants a and b determining the fit. As
+ such it can be specified for either ice or water, or adapted as originally impelemented in ICON,
in which case PvT and TvT need to be substituted by Pv0 and T0.
-
+
Args:
T: temperature in kelvin
@@ -237,7 +258,7 @@ def liq_tetens(T):
This equation is what is used in the ICON code, hence its inclusion in this library. The original
ICON implementation followed Murray's choice of constants (T0=273.15, Pv0=610.78, a=17.269, b=35.86).
- This implementation is referenced to the triple point values of temperature and vapor and with
+ This implementation is referenced to the triple point values of temperature and vapor and with
revised constants (a,b) chosen to better agree with the fits of Wagner and Pruss
Args:
@@ -250,10 +271,10 @@ def liq_tetens(T):
>>> liq_tetens(np.asarray([273.16,305.]))
array([ 611.655 , 4719.73680592])
"""
- a = 17.41463775
+ a = 17.41463775
b = 33.6393413
- return tetens(T,a,b)
+ return tetens(T, a, b)
def ice_tetens(T):
@@ -261,7 +282,7 @@ def ice_tetens(T):
This equation is what is used in the ICON code, hence its inclusion in this library. The original
ICON implementation followed Murray's choice of constants (T0=273.15, Pv0=610.78, a=21.875, b=7.66).
- This implementation is referenced to the triple point values of temperature and vapor and with
+ This implementation is referenced to the triple point values of temperature and vapor and with
revised constants (a,b) chosen to better agree with the fits of Wagner and Pruss
Args:
@@ -275,6 +296,6 @@ def ice_tetens(T):
array([611.655 , 196.10072658])
"""
a = 22.0419977
- b = 5.
+ b = 5.0
- return tetens(T,a,b)
+ return tetens(T, a, b)
diff --git a/setup.py b/setup.py
index 10519fc..a558231 100644
--- a/setup.py
+++ b/setup.py
@@ -3,7 +3,7 @@ from setuptools import setup, find_packages
setup(
name="moist_thermodynamics",
- version="0.3",
+ version="0.4",
description="Constants and functions for the treatment of moist atmospheric thermodynamics",
packages=find_packages(),
install_requires=[
--
GitLab