VaporPressure is called for the calculation of specific humidity of air and sea surface

VaporPressure.py

+import numpy as np
+import logging
+"""
+    Calculate the saturation vapor pressure
+
+    Example: Calculate the vapor pressure over liquid water using the WMO
+    formula
+    temp = FINDGEN(100) - 100.
+    P = VaporPressure(temp, 'liquid', 'WMO')
+
+    For temperatures above 0 deg C the vapor pressure over liquid water
+    is calculated.
+
+    The optional parameter 'liquid' changes the calculation to vapor pressure
+    over liquid water over the entire temperature range.
+    The formula to be used can be selected with the appropriate keyword
+
+    The current default fomulas are Hyland and Wexler for liquid and
+    Goff Gratch for ice.
+
+    Parameters
+    ----------
+    temp : float
+        Temperature in [deg C]
+    phase : str
+        'liquid' : Calculate vapor pressure over liqiud water or
+        'ice' : Calculate vapor pressure over ice
+    meth : str
+        formula to be used
+        MartiMauersberger   : vaporpressure formula from Marti Mauersberger
+        MagnusTetens        : vaporpressure formula from Magnus Tetens
+        GoffGratch          : vaporpressure formula from Goff Gratch
+        Buck_original       : vaporpressure formula from Buck (original)
+        Buck_manual         : vaporpressure formula from the Buck manual
+        CIMO                : vaporpressure formula recommended by CIMO
+        Vaisala             : vaporpressure formula from Vaisala
+        WMO_Goff            : vaporpressure formula from WMO, which should have been Goff
+        WMO2000             : vaporpressure formula from WMO (2000) containing a typo
+        Wexler              : vaporpressure formula from Wexler (1976)
+        Sonntag             : vaporpressure formula from Sonntag (1994)
+        Bolton              : vaporpressure formula from Bolton (1980)
+        HylandWexler        : vaporpressure formula from Hyland and Wexler (1983)
+        IAPWS               : vaporpressure formula from IAPWS (2002)
+        Preining            : vaporpressure formula from Preining (2002)
+        MurphyKoop          : vaporpressure formula from Murphy and Koop (2005)
+
+    Returns
+    -------
+    P : float
+        Saturation vapor pressure in [hPa]
+
+
+    Author
+    ------
+    Holger Voemel
+    NCAR / EOL
+    PO BOX 3000
+    Boulder, CO 80303
+    USA
+    Ported to Python and modified by S. Biri
+
+    COPYRIGHT
+    ---------
+    Copyright (c) 2015, Holger Voemel
+"""
+def VaporPressure(temp, P, phase, meth):
+    logging.basicConfig(filename='VaporPressure.log',
+                        format='%(asctime)s %(message)s',level=logging.info)
+
+    Psat = np.zeros(temp.size)*np.nan
+    if (np.nanmin(temp) > 200):  # if Ta in Kelvin convert to Celsius
+        temp = temp-273.16
+    T = np.copy(temp)+273.16    #  Most formulas use T in [K]
+                                #  Formulas using [C] use the variable temp
+    #  Default uses Hyland and Wexler over liquid. While this may not be
+    #  the best formula, it is consistent with what Vaisala uses in their
+    #  system
+    #  Calculate saturation pressure over liquid water
+    if (phase == 'liquid'):
+        if (meth == 'HylandWexler' or meth == ''):
+            logging.info("Source 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.")
+            Psat = np.exp(-0.58002206e4/T+0.13914993e1-0.48640239e-1*T +
+                          0.41764768e-4*np.power(T, 2) -
+                          0.14452093e-7*np.power(T, 3) +
+                          0.65459673e1*np.log(T))/100
+        if (meth == 'Hardy'):
+            logging.info("Source Hardy, B., 1998, ITS-90 Formulations \
+                         for Vapor Pressure, Frostpoint Temperature, \
+                         Dewpoint Temperature, and Enhancement Factors \
+                         in the Range -100 to +100°C, The Proceedings of \
+                         the Third International Symposium on Humidity & \
+                         Moisture, London, England")
+            Psat = np.exp(-2.8365744e3/np.power(T, 2)-6.028076559e3/T +
+                          1.954263612e1-2.737830188e-2*T +
+                          1.6261698e-5*np.power(T, 2) +
+                          7.0229056e-10*np.power(T, 3) -
+                          1.8680009e-13*np.power(T, 4) +
+                          2.7150305*np.log(T))/100
+        if (meth == 'Preining'):
+            logging.info("Source : Vehkamaeki, H., M. Kulmala, I. Napari, \
+                         K. E. J. Lehtinen, C.Timmreck, M. Noppel, and \
+                         A. Laaksonen (2002), J. Geophys. Res., 107, \
+                         doi:10.1029/2002JD002184.")
+            Psat = np.exp(-7235.424651/T+77.34491296+5.7113e-3*T -
+                          8.2*np.log(T))/100
+        if (meth == 'Wexler'):
+            logging.info("Wexler, A., Vapor Pressure Formulation for Water \
+                         in Range 0 to 100 C. A Revision, Journal of Research \
+                         of the National Bureau of Standards - A. Physics and \
+                         Chemistry, September - December 1976, Vol. 80A, \
+                         Nos.5 and 6, 775-785")
+            Psat = np.exp(-0.29912729e4*np.power(T, -2) -
+                          0.60170128e4*np.power(T, -1) +
+                          0.1887643854e2-0.28354721e-1*T +
+                          0.17838301e-4*np.power(T, 2) -
+                          0.84150417e-9*np.power(T, 3) +
+                          0.44412543e-12*np.power(T, 4) +  # This line was corrected from '-' to '+' following the original citation. (HV 20140819). The change makes only negligible difference
+                          2.858487*np.log(T))/100
+        if ((meth == 'GoffGratch') or (meth == 'MartiMauersberger')):
+            # logging.info("Marti and Mauersberger don't have a vapor pressure \
+            #              curve over liquid. Using Goff Gratch instead; \
+            #              Goff Gratch formulation Source : Smithsonian \
+            #              Meteorological Tables, 5th edition, p. 350, 1984 \
+            #              From original source: Goff & Gratch (1946), p. 107) \
+            #              Goff Gratch formulation; Source : Smithsonian \
+            #              Meteorological Tables, 5th edition, p. 350, 1984 \
+            #              From 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
+            Psat = np.power(10, -7.90298*(Ts/T-1)+5.02808*np.log10(Ts/T) -
+                            1.3816e-7*(np.power(10, 11.344*(1-T/Ts))-1) +
+                            8.1328e-3*(np.power(10, -3.49149*(Ts/T-1))-1) +
+                            np.log10(ews))
+        if (meth == 'CIMO'):
+            logging.info("Source: Annex 4B, Guide to Meteorological \
+                         Instruments and Methods of Observation, \
+                         WMO Publication No 8, 7th edition, Geneva, 2008. \
+                         (CIMO Guide)")
+            Psat = (6.112*np.exp(17.62*temp/(243.12+temp)) *
+                    (1.0016+3.15e-6*P-0.074/P))
+        if (meth == 'MagnusTetens'):
+            logging.info("Source: Murray, F. W., On the computation of \
+                         saturation vapor pressure, J. Appl. Meteorol., \
+                         6, 203-204, 1967.")
+            Psat = np.power(10, 7.5*(temp)/(temp+237.5)+0.7858)
+            # Murray quotes this as the original formula and
+            Psat = 6.1078*np.exp(17.269388*temp/(temp+237.3))
+            # this as the mathematical aquivalent in the form of base e.
+        if (meth == 'Buck'):
+            logging.info("Bucks vapor pressure formulation based on \
+                          Tetens formula. Source: Buck, A. L., \
+                          New equations for computing vapor pressure and \
+                          enhancement factor, J. Appl. Meteorol., 20, \
+                          1527-1532, 1981.")
+            Psat = (6.1121*np.exp(17.502*temp/(240.97+temp)) *
+                    (1.0007+(3.46e-6*P)))
+        if (meth == 'Buck2'):
+            logging.info("Bucks vapor pressure formulation based on \
+                         Tetens formula. Source: Buck Research, \
+                         Model CR-1A Hygrometer Operating Manual, Sep 2001")
+            Psat = (6.1121*np.exp((18.678-(temp)/234.5)*(temp)/(257.14+temp)) *
+                    (1+1e-4*(7.2+P*(0.0320)+5.9e-6*np.power(T, 2))))
+        if (meth == 'WMO'):
+            logging.info("Intended WMO formulation, originally published by \
+                         Goff (1957) incorrectly referenced by WMO technical \
+                         regulations, WMO-NO 49, Vol I, General \
+                         Meteorological Standards and Recommended Practices, \
+                         App. A, Corrigendum Aug 2000. and incorrectly \
+                         referenced by WMO technical regulations, \
+                         WMO-NO 49, Vol I, General Meteorological Standards \
+                         and Recommended Practices, App. A, 1988.")
+            Ts = 273.16  # triple point temperature in K
+            Psat = np.power(10, 10.79574*(1-Ts/T)-5.028*np.log10(T/Ts) +
+                            1.50475e-4*(1-np.power(10, -8.2969*(T/Ts-1))) +
+                            0.42873e-3*(np.power(10, -4.76955*(1-Ts/T))-1) +
+                            0.78614)
+        if (meth == 'WMO2000'):
+            logging.info("WMO formulation, which is very similar to Goff \
+                         Gratch. Source : WMO technical regulations, \
+                         WMO-NO 49, Vol I, General Meteorological Standards \
+                         and Recommended Practices, App. A, \
+                         Corrigendum Aug 2000.")
+            Ts = 273.16  # triple point temperature in K
+            Psat = np.power(10, 10.79574*(1-Ts/T)-5.028*np.log10(T/Ts) +
+                            1.50475e-4*(1-np.power(10, -8.2969*(T/Ts-1))) +
+                            0.42873e-3*(np.power(10, -4.76955*(1.-Ts/T))-1) +
+                            0.78614)
+        if (meth == 'Sonntag'):
+            logging.info("Source: Sonntag, D., Advancements in the field of \
+                         hygrometry, Meteorol. Z., N. F., 3, 51-66, 1994.")
+            Psat = np.exp(-6096.9385*np.power(T, -1)+16.635794 -
+                          2.711193e-2*T+1.673952e-5*np.power(T, 2) +
+                          2.433502*np.log(T))#*(1.0016+P*3.15e-6-0.074/P)
+        if (meth == 'Bolton'):
+            logging.info("Source: Bolton, D., The computation of equivalent \
+                         potential temperature, Monthly Weather Report, 108, \
+                         1046-1053, 1980. equation (10)")
+            Psat = 6.112*np.exp(17.67*temp/(temp+243.5))
+        if (meth == 'IAPWS'):
+            logging.info("Source: 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 : Vapor pressure at the critical point
+            nu = (1-T/Tc)
+            a1, a2, a3 = -7.85951783, 1.84408259, -11.7866497
+            a4, a5, a6 = 22.6807411, -15.9618719, 1.80122502
+            Psat = (Pc*np.exp(Tc/T*(a1*nu+a2*np.power(nu, 1.5) +
+                    a3*np.power(nu, 3)+a4*np.power(nu, 3.5) +
+                    a5*np.power(nu, 4)+ a6*np.power(nu, 7.5))))
+        if (meth == 'MurphyKoop'):
+            logging.info("Source : 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.")
+            Psat = np.exp(54.842763-6763.22/T-4.210*np.log(T)+0.000367*T +
+                          np.tanh(0.0415*(T-218.8))*(53.878-1331.22/T -
+                          9.44523*np.log(T)+0.014025*T))/100
+    # Calculate saturation pressure over ice ----------------------------------
+    elif (phase == 'ice'):
+        logging.info("Default uses Goff Gratch over ice. There is little \
+                     ambiguity in the ice saturation curve. Goff Gratch \
+                     is widely used.")
+        if (meth == 'MartiMauersberger'):
+            logging.info("Source : Marti, J. and K Mauersberger, A survey and \
+                         new measurements of ice vapor pressure at \
+                         temperatures between 170 and 250 K, GRL 20, \
+                         363-366, 1993.")
+            Psat = np.power(10, -2663.5/T+12.537)/100
+        if (meth == 'HylandWexler'):
+            logging.info("Source 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.")
+            Psat = np.exp(-0.56745359e4/T+0.63925247e1-0.96778430e-2*T +
+                          0.62215701e-6*np.power(T, 2) +
+                          0.20747825e-8*np.power(T, 3) -
+                          0.9484024e-12*np.power(T, 4) +
+                          0.41635019e1*np.log(T))/100
+        if (meth == 'Wexler'):
+            logging.info("Wexler, A., Vapor pressure formulation for ice, \
+                         Journal of Research of the National Bureau of \
+                         Standards-A. 81A, 5-20, 1977.")
+            Psat = np.exp(-0.58653696e4*np.power(T, -1)+0.22241033e2 +
+                          0.13749042e-1*T-0.34031775e-4*np.power(T, 2) +
+                          0.26967687e-7*np.power(T, 3) +
+                          0.6918651*np.log(T))/100
+        if (meth == 'Hardy'):
+            logging.info("Source Hardy, B., 1998, ITS-90 Formulations for \
+                         Vapor Pressure, Frostpoint Temperature, Dewpoint \
+                         Temperature, and Enhancement Factors in the Range \
+                         -100 to +100°C, The Proceedings of the Third \
+                         International Symposium on Humidity & Moisture, \
+                         London, England. These coefficients are updated to \
+                         ITS90 based on the work by Bob Hardy at Thunder \
+                         Scientific: http://www.thunderscientific.com/\
+                         tech_info/reflibrary/its90formulas.pdf \
+                         The difference to the older ITS68 coefficients used \
+                         by Wexler is academic.")
+            Psat = np.exp(-0.58666426e4*np.power(T, -1)+0.2232870244e2 +
+                          0.139387003e-1*T-0.34262402e-4*np.power(T, 2) +
+                          0.27040955e-7*np.power(T, 3) +
+                          0.67063522e-1*np.log(T))/100
+        if (meth == 'GoffGratch' or meth == '' or meth == 'IAPWS'):
+            logging.info("IAPWS does not provide a vapor pressure formulation \
+                         over ice use Goff Gratch instead.\
+                         Source : Smithsonian Meteorological Tables, \
+                         5th edition, p. 350, 1984")
+            ei0 = 6.1071  # mbar
+            T0 = 273.16   # triple point in K
+            Psat = np.power(10, -9.09718*(T0/T-1)-3.56654*np.log10(T0/T) +
+                            0.876793*(1-T/T0)+np.log10(ei0))
+        if (meth == 'MagnusTetens'):
+            logging.info("Source: Murray, F. W., On the computation of \
+                         saturation vapor pressure, J. Appl. Meteorol., 6, \
+                         203-204, 1967.")
+            Psat = np.power(10, 9.5*temp/(265.5+temp)+0.7858)
+            #  Murray quotes this as the original formula and
+            Psat = 6.1078*np.exp(21.8745584*(T-273.16)/(T-7.66))
+            # this as the mathematical aquivalent in the form of base e.
+        if (meth == 'Buck'):
+            logging.info("Bucks vapor pressure formulation based on Tetens \
+                         formula. Source: Buck, A. L., New equations for \
+                         computing vapor pressure and enhancement factor, \
+                         J. Appl. Meteorol., 20, 1527-1532, 1981.")
+            Psat = (6.1115*np.exp(22.452*temp/(272.55+temp)) *
+                    (1.0003+(4.18e-6*P)))
+        if (meth == 'Buck2'):
+            logging.info("Bucks vapor pressure formulation based on Tetens \
+                         formula. Source: Buck Research, Model CR-1A \
+                         Hygrometer Operating Manual, Sep 2001")
+            Psat = (6.1115*np.exp((23.036-temp/333.7)*temp/(279.82+temp)) *
+                    (1+1e-4*(2.2+P*(0.0383+6.4e-6*np.power(T, 2)))))
+        if (meth == 'CIMO'):
+            logging.info("Source: Annex 4B, Guide to Meteorological \
+                         Instruments and Methods of Observation, \
+                         WMO Publication No 8, 7th edition, Geneva, 2008. \
+                         (CIMO Guide)")
+            Psat = (6.112*np.exp(22.46*temp/(272.62+temp)) *
+                    (1.0016+3.15e-6*P-0.074/P))
+        if (meth == 'WMO' or meth == 'WMO2000'):
+            logging.info("There is no typo issue in the WMO formulations for \
+                         ice. WMO formulation, which is very similar to Goff \
+                         Gratch. Source : WMO technical regulations, \
+                         WMO-NO 49, Vol I, General Meteorological Standards \
+                         and Recommended Practices, Aug 2000, App. A.")
+            T0 = 273.16  # triple point temperature in K
+            Psat = np.power(10, -9.09685*(T0/T-1)-3.56654*np.log10(T0/T) +
+                            0.87682*(1-T/T0)+0.78614)
+        if (meth == 'Sonntag'):
+            logging.info("Source: Sonntag, D., Advancements in the field of \
+                         hygrometry, Meteorol. Z., N. F., 3, 51-66, 1994.")
+            Psat = np.exp(-6024.5282*np.power(T, -1)+24.721994+1.0613868e-2*T -
+                          1.3198825e-5*np.power(T, 2)-0.49382577*np.log(T))
+        if (meth == 'MurphyKoop'):
+            logging.info("Source : 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.")
+            Psat = np.exp(9.550426-5723.265/T+3.53068*np.log(T) -
+                          0.00728332*T)/100
+        s = np.where(temp > 0)
+        if (s.size[0] >= 1):
+            logging.info("Independent of the formula used for ice, use \
+                         Hyland Wexler (water) for temperatures above freezing\
+                         (see above). Source 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.")
+            Psat_w = np.exp(-0.58002206e4/T+0.13914993e1-0.48640239e-1*T +
+                            0.41764768e-4*np.power(T, 2) -
+                            0.14452093e-7*np.power(T, 3) +
+                            0.65459673e1*np.log(T))/100
+            Psat[s] = Psat_w[s]
+
+    return Psat