# Source code for hesseflux.esat

```#!/usr/bin/env python
"""
esat : Saturation vapour pressure of water and ice.

This module was written by Matthias Cuntz while at Department of
Computational Hydrosystems, Helmholtz Centre for Environmental
Research - UFZ, Leipzig, Germany, and continued while at Institut
National de Recherche pour l'Agriculture, l'Alimentation et
l'Environnement (INRAE), Nancy, France.

Copyright (c) 2012-2020 Matthias Cuntz - mc (at) macu (dot) de

* Written Jan 2012 by Matthias Cuntz (mc (at) macu (dot) de)
* Ported to Python 3, Feb 2013, Matthias Cuntz
* Changed handling of masked arrays, Oct 2013, Matthias Cuntz
* Assert T>0, Apr 2014, Matthias Cuntz
* Using numpy docstring format, May 2020, Matthias Cuntz

.. moduleauthor:: Matthias Cuntz

The following functions are provided

.. autosummary::
esat
"""
from __future__ import division, absolute_import, print_function
import numpy as np

__all__ = ['esat']

[docs]def esat(T, liquid=False, formula='GoffGratch'):
"""
Calculates the saturation vapour pressure of water and/or ice.

For temperatures above (and equal) 0 degree C (273.15 K),
the vapour pressure over liquid water is calculated.
For temperatures below 0 degree C, the vapour pressure over ice is calculated.

The optional parameter liquid=True changes the calculation to vapour pressure
over liquid water over the entire temperature range.

Parameters
----------
T : float or array_like
Temperature [K]
liquid : bool, optional
If True, use liquid formula for all temperatures.
formula : str, optional
Name of reference to use for calculations, case-insensitive (default: GoffGratch).

Note that several formulations do not provide a vapour pressure formulation over ice
and Marti and Mauersberger do not provide a formula over liquid: GoffGratch is used in theses cases.

GoffGratch:        Smithsonian Tables, 1984; after Goff and Gratch, 1946 (default)

MartiMauersberger: Marti and Mauersberger, 1993

MagnusTeten:       Murray, 1967

Buck_original:     Buck, 1981

Buck:              Buck Research Manual, 1996

WMO:               Goff, 1957; WMO 1988, 2000

Wexler:            Wexler, 1977

Sonntag:           Sonntag, 1994

Bolton:            Bolton, 1980

Fukuta:            Fukuta, N. and C. M. Gramada, 2003

HylandWexler:      Hyland and Wexler, 1983

IAPWS:             Wagner and Pruss, 2002

MurphyKoop:        Murphy and Koop, 2005

Returns
-------
float or array_like
Saturation water pressure at temperature T in Pascal [Pa].

Notes
-----

Referred literature cited in code.

Examples
--------
>>> print('{:.3f}'.format(esat(293.15)))
2335.847
>>> print('{:.3f}'.format(esat(253.15)))
103.074

>>> print('{:.3f} {:.3f}'.format(*esat([293.15,253.15])))
2335.847 103.074
>>> print('{:.3f} {:.3f}'.format(*esat([293.15,253.15],formula='GoffGratch')))
2335.847 103.074
>>> print('{:.3f} {:.3f}'.format(*esat([293.15,253.15],formula='MartiMauersberger')))
2335.847 103.650
>>> print('{:.3f} {:.3f}'.format(*esat([293.15,253.15],formula='MagnusTeten')))
2335.201 102.771
>>> print('{:.3f} {:.3f}'.format(*esat([293.15,253.15],formula='buck')))
2338.340 103.286
>>> print('{:.3f} {:.3f}'.format(*esat([293.15,253.15],formula='Buck_original')))
2337.282 103.267
>>> print('{:.3f} {:.3f}'.format(*esat([293.15,253.15],formula='wmo')))
2337.080 103.153
>>> print('{:.3f} {:.3f}'.format(*esat([293.15,253.15],formula='WEXLER')))
2323.254 103.074
>>> print('{:.3f} {:.3f}'.format(*esat([293.15,253.15],formula='Sonntag')))
2339.249 103.249
>>> print('{:.3f} {:.3f}'.format(*esat([293.15,253.15],formula='Bolton')))
2336.947 103.074
>>> print('{:.3f} {:.3f}'.format(*esat([293.15,253.15],formula='Fukuta')))
2335.847 103.074
>>> print('{:.3f} {:.3f}'.format(*esat([293.15,253.15],formula='HylandWexler')))
2338.804 103.260
>>> print('{:.3f} {:.3f}'.format(*esat([293.15,253.15],formula='IAPWS')))
2339.194 103.074
>>> print('{:.3f} {:.3f}'.format(*esat([293.15,253.15],formula='MurphyKoop')))
2339.399 103.252

>>> print('{:.3f} {:.3f}'.format(*esat(np.array([293.15,253.15]), liquid=True)))
2335.847 125.292
>>> print('{:.3f} {:.3f}'.format(*esat([293.15,253.15],formula='Fukuta', liquid=True)))
2335.847 125.079

>>> print('{:.3f} {:.3f}'.format(*esat(np.array([293.15,393.15]))))
esat.py: UserWarning: T>373.15 K - something might be wrong with T.
2335.847 198473.378
>>> print('{:.3f} {:.3f}'.format(*esat(np.array([293.15,93.15]))))
esat.py: UserWarning: T<100 - T probably given in Celsius instead of Kelvin.
2335.847 0.000

>>> print('{:.3f} {:.3f}'.format(*out.filled(-9999.)))
103.074 -9999.000

History
-------
Written,  Matthias Cuntz, Jan 2012
Modified, Matthias Cuntz, Feb 2013 - ported to Python 3
Matthias Cuntz, Oct 2013 - changed masked array handling
Matthias Cuntz, Apr 2014 - assert
Matthias Cuntz, May 2020 - numpy docstring format
"""
#
# Constants
T0 = 273.15 # Celcius <-> Kelvin [K]
knownforms = (['Buck', 'Buck_original', 'Bolton', 'Fukuta', 'GoffGratch', 'HylandWexler',
'IAPWS', 'MagnusTeten', 'MartiMauersberger', 'MurphyKoop',
'Sonntag', 'Vaisala', 'Wexler', 'WMO'])
lknown = [i.lower() for i in knownforms]
#
# Check input
T = np.ma.array(T)
assert not np.ma.any(T <= 0.), 'T<0 - T probably given in Celsius instead of Kelvin.'
if np.ma.any(T < 100.):
print("esat.py: UserWarning: T<100 - T probably given in Celsius instead of Kelvin.")
if np.ma.any(T > (T0+100.)):
print("esat.py: UserWarning: T>373.15 K - something might be wrong with T.")
form = formula.lower()
if form not in lknown:
raise ValueError('Formula not know. Known formulas are {:s}'.format(knownforms))
#
# Split input into masked arrays
if liquid == True:
Tlim = 1e-3
else:
Tlim = T0

if T.size > 1:
isone = False
ii = np.ma.where(T>=Tlim)[0]
jj = np.ma.where(T<Tlim)[0]
if np.size(ii) > 0:
T_liq  = T[ii]
if np.size(jj) > 0:
T_ice  = T[jj]
else:
isone = True
if T>=Tlim:
ii = [0]
jj = []
T_liq = T
else:
ii = []
jj = [0]
T_ice = T
esat_out = T.copy() # to conserve mask
#
# Calc
#
# Liquid
if np.size(ii) > 0:
TC_liq = T_liq - T0
if form == 'buck':
'''Bucks vapour pressure formulation based on Tetens formula
Buck Research, Model CR-1A Hygrometer Operating Manual, Sep 2001'''
esat_liq = 6.1121 * np.ma.exp((18.678 - (TC_liq) / 234.5) * (TC_liq) / (257.14+TC_liq))
elif form == 'buck_original':
'''Bucks vapour pressure formulation based on Tetens formula
Buck, A. L., New equations for computing vapour pressure and enhancement factor,
J. Appl. Meteorol., 20, 1527-1532, 1981.'''
esat_liq = 6.1121 * np.ma.exp(17.502 * TC_liq / (240.97+TC_liq))
elif form == 'bolton':
'''Bolton, D., The computation of equivalent potential temperature
Monthly Weather Report, 108, 1046-1053, 1980. equation (10)'''
esat_liq = 6.112 * np.ma.exp(17.67 * TC_liq / (TC_liq+243.5))
elif form == 'fukuta':
'''Fukuta, N. and C. M. Gramada, Vapour pressure measurement of supercooled water
J. Atmos. Sci., 60, 1871-1875, 2003.
This paper does not give a vapour pressure formulation, but rather a correction over the Smithsonian Tables.
Thus calculate the table value first, then use the correction to get to the measured value.
This is done only for -39<TC<0.'''
Ts    = 373.16    # steam point temperature in K
ews   = 1013.246  # saturation pressure at steam point temperature, normal atmosphere
esat_liq = (10.**(-7.90298*(Ts/T_liq-1.) + 5.02808 * np.ma.log10(Ts/T_liq)
- 1.3816e-7 * (10.**(11.344*(1.-T_liq/Ts))-1.)
+ 8.1328e-3 * (10.**(-3.49149*(Ts/T_liq-1))-1.) + np.ma.log10(ews)))
mm = (TC_liq < 0.) & (TC_liq > -39.)
if np.ma.any(mm):
x = TC_liq + 19.
esat_liq = (np.where(mm, esat_liq * (0.9992 + 7.113e-4*x - 1.847e-4*x**2 + 1.189e-5*x**3
+ 1.130e-7*x**4 - 1.743e-8*x**5), esat_liq))
elif ((form == 'goffgratch') | (form == 'martimauersberger')):
'''Goff Gratch formulation, Smithsonian Meteorological Tables, 5th edition, p. 350, 1984
Original source: Goff and Gratch (1946), p. 107.'''
Ts    = 373.16    # steam point temperature in K
ews   = 1013.246  # saturation pressure at steam point temperature, normal atmosphere
esat_liq = (10.**(-7.90298*(Ts/T_liq-1.) + 5.02808 * np.ma.log10(Ts/T_liq)
- 1.3816e-7 * (10.**(11.344*(1.-T_liq/Ts))-1.)
+ 8.1328e-3 * (10.**(-3.49149*(Ts/T_liq-1.))-1.) + np.ma.log10(ews)))
elif form == 'hylandwexler':
'''Hyland, R. W. and A. Wexler, Formulations for the Thermodynamic Properties of the
saturated Phases of H2O from 173.15K to 473.15K, ASHRAE Trans, 89(2A), 500-519, 1983.'''
esat_liq = (np.ma.exp(- 0.58002206e4/T_liq + 0.13914993e1 - 0.48640239e-1*T_liq
+ 0.41764768e-4*T_liq**2 - 0.14452093e-7*T_liq**3
+ 0.65459673e1 * np.ma.log(T_liq)) / 100.)
elif form == 'iapws':
'''Wagner W. and A. Pruss (2002)
The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance
for general and scientific use, J. Phys. Chem. Ref. Data, 31(2), 387-535.
This is the 'official' formulation from the International Association for the Properties of Water and Steam
The valid range of this formulation is 273.16 <= T <= 647.096 K and is based on the ITS90 temperature scale.'''
Tc = 647.096      # K   : Temperature at the critical point
Pc = 22.064e4     # hPa : Vapour pressure at the critical point
nu = (1-T_liq/Tc)
a1 = -7.85951783
a2 = 1.84408259
a3 = -11.7866497
a4 = 22.6807411
a5 = -15.9618719
a6 = 1.80122502
esat_liq = Pc * np.ma.exp(Tc/T_liq * (a1*nu + a2*nu**1.5 + a3*nu**3 + a4*nu**3.5 + a5*nu**4 + a6*nu**7.5))
elif form == 'magnusteten':
'''Murray, F. W., On the computation of saturation vapour pressure, J. Appl. Meteorol., 6, 203-204, 1967.'''
esat_liq = 10.**(7.5*(TC_liq)/(TC_liq+237.5) + 0.7858)
elif form == 'murphykoop':
'''Murphy and Koop, Review of the vapour pressure of ice and supercooled water for atmospheric applications
Q. J. R. Meteorol. Soc (2005), 131, pp. 1539-1565.'''
esat_liq = (np.exp(54.842763 - 6763.22 / T_liq - 4.210 * np.ma.log(T_liq) + 0.000367 * T_liq
+ np.tanh(0.0415 * (T_liq - 218.8)) * (53.878 - 1331.22 / T_liq -
9.44523 * np.ma.log(T_liq) + 0.014025 * T_liq)) / 100.)
elif form == 'sonntag':
'''Sonntag, D., Advancements in the field of hygrometry, Meteorol. Z., N. F., 3, 51-66, 1994.'''
esat_liq = (np.ma.exp( -6096.9385 * 1./T_liq + 16.635794 - 2.711193e-2 * T_liq
+ 1.673952e-5 * T_liq**2 + 2.433502 * np.ma.log(T_liq)))
elif form == 'wexler':
'''Wexler, A., Vapour pressure formulation for ice
Journal of Research of the National Bureau of Standards-A. 81A, 5-20, 1977.'''
esat_liq = (np.ma.exp(- 2.9912729e3 * 1./T_liq**2 - 6.0170128e3 * 1./T_liq
+ 1.887643854e1 - 2.8354721e-2 * T_liq**1
+ 1.7838301e-5 * T_liq**2 - 8.4150417e-10 * T_liq**3
- 4.4412543e-13 * T_liq**4 + 2.858487 * np.ma.log(T_liq)) / 100.)
elif form == 'wmo':
incorrectly referenced by WMO technical regulations, WMO-NO 49, Vol I
General Meteorological Standards and Recommended Practices, App. A, 1988, Corrigendum Aug 2000.'''
Ts   = 273.16    # steam point temperature in K
esat_liq = (10.**(10.79574*(1.-Ts/T_liq) - 5.02800 * np.ma.log10(T_liq/Ts)
+ 1.50475e-4 * (1.-10.**(-8.2969*(T_liq/Ts-1.)))
+ 0.42873e-3 * (10.**(+4.76955*(1.-Ts/T_liq))-1.)
+ 0.78614))
else:
raise ValueError("formulation not known: {:s}".format(formula))
esat_liq *= 100.
if isone:
esat_out = esat_liq
else:
esat_out[ii] = esat_liq
#
# Ice
if np.size(jj) > 0:
TC_ice = T_ice - T0
if form == 'buck':
'''Bucks vapour pressure formulation based on Tetens formula
Buck Research, Model CR-1A Hygrometer Operating Manual, Sep 2001'''
esat_ice = 6.1115 * np.exp((23.036 - TC_ice / 333.7) * TC_ice / (279.82+TC_ice))
elif form == 'buck_original':
'''Bucks vapour pressure formulation based on Tetens formula
Buck, A. L., New equations for computing vapour pressure and enhancement factor,
J. Appl. Meteorol., 20, 1527-1532, 1981.'''
esat_ice = 6.1115 * np.exp(22.452 * TC_ice / (272.55+TC_ice))
elif ((form == 'goffgratch') | (form == 'bolton') | (form == 'fukuta') | (form == 'iapws') | (form == 'wexler')):
'''Smithsonian Meteorological Tables, 5th edition, p. 350, 1984'''
ei0    = 6.1071   # mbar
Ts     = 273.16   # freezing point in K
esat_ice = np.ma.exp(np.log(10.)*(-9.09718 * (Ts/T_ice-1.) - 3.56654 * np.ma.log10(Ts/T_ice)
+ 0.876793 * (1.-T_ice/Ts) + np.log10(ei0)))
elif (form == 'hylandwexler'):
'''Hyland, R. W. and A. Wexler, Formulations for the Thermodynamic Properties of
the saturated Phases of H2O from 173.15K to 473.15K, ASHRAE Trans, 89(2A), 500-519, 1983.'''
esat_ice = (np.exp(- 0.56745359E4 / T_ice + 0.63925247E1 - 0.96778430E-2 * T_ice
+ 0.62215701E-6 * T_ice**2. + 0.20747825E-8 * T_ice**3.
- 0.94840240E-12 * T_ice**4. + 0.41635019E1 * np.log(T_ice)) / 100.)
elif form == 'magnusteten':
'''Murray, F. W., On the computation of saturation vapo rpressure, J. Appl. Meteorol., 6, 203-204, 1967.'''
esat_ice = 10.**(9.5 * TC_ice/(265.5+TC_ice) + 0.7858)
elif form == 'martimauersberger':
'''Marti, J. and K Mauersberger, A survey and new measurements of ice vapour pressure
at temperatures between 170 and 250 K, GRL 20, 363-366, 1993.'''
esat_ice = 10.**(-2663.5/T_ice + 12.537) / 100.
elif form == 'murphykoop':
'''Murphy and Koop, Review of the vapour pressure of ice and supercooled water
for atmospheric applications, Q. J. R. Meteorol. Soc (2005), 131, pp. 1539-1565.'''
esat_ice = np.exp(9.550426 - 5723.265/T_ice + 3.53068 * np.log(T_ice) - 0.00728332 * T_ice) / 100.
elif form == 'sonntag':
'''Sonntag, D., Advancements in the field of hygrometry, Meteorol. Z., N. F., 3, 51-66, 1994.'''
esat_ice = (np.exp(- 6024.5282 * 1./T_ice + 24.721994 + 1.0613868E-2 * T_ice
- 1.3198825E-5 * T_ice**2 - 0.49382577 * np.log(T_ice)))
elif form == 'wmo':
'''WMO formulation, which is very similar to Goff Gratch
WMO technical regulations, WMO-NO 49, Vol I, General Meteorological Standards
and Recommended Practices, Aug 2000, App. A.'''
Ts    = 273.16    # steam point temperature in K
esat_ice = (10.**(-9.09685 * (Ts/T_ice-1.) - 3.56654 * np.log10(Ts/T_ice)
+ 0.87682 * (1.-T_ice/Ts) + 0.78614))
else:
raise ValueError("formulation not known: {:s}".format(formula))
esat_ice *= 100.
if isone:
esat_out = esat_ice
else:
esat_out[jj] = esat_ice
#
# Finish
return esat_out

if __name__ == '__main__':
import doctest
doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE)
```