Saturation Vapor Pressure Comparisons

svp.py

import numpy as np
import matplotlib.pyplot as plt
from metpy.units import units
from metpy.calc import saturation_vapor_pressure
import timeit

def goffgratch46(Tin):
    """Computes saturation vapor pressure (mb) from temperature using the
    Goff and Gratch (1946) formulation as quoted by Alduchov and Eskridge
    (1996)"""

    T = Tin.to('degK').m

    Ts = 373.16
    ews = 1013.246

    loge = -7.90298*(Ts/T-1) + 5.02808*np.log10(Ts/T) -  \
           1.3816e-7*(10**(11.344*(1-T/Ts)) - 1) +  \
           8.1328e-3*(10**(-3.49149*(Ts/T-1)) - 1) + np.log10(ews)
    return 10**loge * units.millibar

def goff65(Tin):
    """Computes saturation vapor pressure (mb) from temperature using the
    Goff (1965) formulation as quoted by Gibbins (1990)"""

    T = Tin.to('degK').m

    e1 = 1013.25
    T01 = 273.16

    T01T = T01/T
    TT01 = T/T01

    ew = 10.79586 * (1-T01T) - 5.02808*np.log10(TT01) +  \
         1.50474e-4*(1 - 10**(-8.29692*(TT01-1))) +  \
         0.42873e-3*(10**(4.76955*(1-T01T)) - 1) - 2.2195983
    return e1 * 10**ew * units.millibar

def wexler76(Tin):
    """Computes saturation vapor pressure (mb) from temperature (K) using
    Wexler's 1976 formula as quoted by Flatau et al. (1992)"""

    g = (-0.29912729e4, -0.60170128e4, 0.1887643854e2, -0.28354721e-1,
         0.17838301e-4, -0.84150417e-9, 0.44412543e-12, 0.2858487e1)

    T = Tin.to('degK').m

    lne = (g[0] + (g[1] + (g[2] + g[7]*np.log(T)
                + (g[3] + (g[4] + (g[5] + g[6]*T)*T)*T)*T)*T)*T) / T**2
    return .01 * np.exp(lne) * units.millibar

def eval_poly(a, x):
    """Evaluate polynomial at x with coefficients a via Horner's method"""

    p = a[-1]
    for i in range(len(a)-2, -1, -1):
        p = p * x + a[i]
    return p

def flatau6(T):
    """Compute saturation vapor pressure (mb) from temperature (K) using
    the 6th-order polynomial fit (Table 3) from Flatau et al. (1992)"""

    a = (6.11176750, 0.443986062, 0.143053301e-1,
         0.265027242e-3, 0.302246994e-5, 0.203886313e-7, 0.638780966e-10)
    
    x = T.to('degC').m

    # Surprisingly, eval_poly is faster than np.polyval
    return eval_poly(a, x) * units.millibar


def gueymard7(T):
    """Compute saturation vapor pressure (mb) from temperature (K) using        
    the 6th-order polynomial fit (Table 3) from Flatau et al. (1992)"""

    a = (6.1104546, 0.4442351, 1.4302099e-2, 2.6454708e-4,
         3.0357098e-6, 2.0972268e-8, 6.0487594e-11, -1.469687e-13)
    
    x = T.to('degC').m

    return eval_poly(a, x) * units.millibar


def murphykoop05(Tin):
    """Compute saturation vapor pressure (mb) from temperature (K) using
    the Murphy and Koop (2005) formulation"""

    T = Tin.to('degK').m
    res = np.exp(54.842763 - 6763.22/T - 4.21*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))
    return .01 * res * units.millibar

def wagnerpruss02(Tin):
    """Compute saturation vapor pressure (mb) from temperature (K) using
    the Wagner and Pruss (2002) formulation"""

    T = Tin.to('degK').m
    Tc = 647.096
    theta = 1 - T/Tc
    pc = 22.064e6
    a = (None, -7.85951783, 1.84408259, -11.7866497, 22.6807411, -15.9618719,
         1.80122502)
    res = pc * np.exp(Tc * (a[1]*theta + a[2]*theta**1.5 + a[3]*theta**3
                + a[4]*theta**3.5 + a[5]*theta**4 + a[6]*theta**7.5) / T)
    return .01 * res * units.millibar


def koutsoyiannis(Tin):
    """Compute saturation vapor pressure (mb) from temperature (K) using
    the Koutsoyiannis (2012) formulation"""

    T = Tin.to('degK').m
    T0 = 273.16
    p0 = 6.11657

    res = p0*np.exp(24.921*(1-T0/T)) * (T0/T)**5.06
    return res * units.millibar


def romps17(Tin):
    """Compute saturation vapor pressure (mb) from temperature (K) using
    the Romps (2017) formulation"""
    
    T = Tin.to('degK').m

    cvv = 1418
    ptrip = 611.65
    Ttrip = 273.16
    E0v = 2.374e6
    E0s = 0.3337e6
    Rv = 461
    cvl = 4119
    cvs = 1861

    cpv = cvv + Rv
    c1 = (cpv - cvl) / Rv
    c2 = (E0v - (cvv-cvl)*Ttrip) / Rv

    res = ptrip * (T/Ttrip)**c1 * np.exp(c2 * (1/Ttrip - 1/T))
    return .01 * res * units.millibar


def alduchov96(Tin):
    """Compute saturation vapor pressure (mb) from temperature using the
    Alduchov and Eskridge (1996) formulation"""

    T = Tin.to('degC').m

    return 6.1094 * np.exp(17.625*T/(243.04+T)) * units.millibar


def huang18(Tin):
    """Compute saturation vapor pressure (mb) from temperature using the
    Huang (2018) formulation"""

    T = Tin.to('degC').m

    res  = np.exp(34.494 - 4924.99 / (T + 237.1)) / (T + 105)**1.57
    return .01 * res * units.millibar



def compare(x, curves, names):
    fig, ax = plt.subplots()
    
    for i, y in enumerate(curves):
        ax.plot(x, y, label=names[i])
    ax.legend()
    ax.set_title('Saturation Vapor Pressure [mb]')
    plt.show()

def compare_abs(f, names, x, svp):
    fig, ax = plt.subplots()
    
    for n in names:
        d = f[n]['svp'] - svp
        ax.plot(x, d, label=n)
    ax.axhline(color='k', ls='--')
    ax.legend()
    ax.set_title('Difference of SVP from Reference Mean')
    plt.show()

def compare_rel(f, names, x, svp):
    fig, ax = plt.subplots()

    for n in names:
        err = 100*np.abs((f[n]['svp'] - svp) / svp)
        ax.semilogy(x, err, label=n)
    ax.legend()
    ax.set_title('Relative Error [%] using Reference Mean as Truth')
    plt.show()

def calc_efficiency(f):
    ref_fun = 'Flatau'
    ref_call = f[ref_fun]['fun'].__name__ + '(Ts)'
    t0 = timeit.timeit(ref_call, number=10, globals=globals())
    print(f'Time for {ref_fun}: {t0:.4f}')
    for name, g in f.items():
        ref_call = g['fun'].__name__ + '(Ts)'
        t = timeit.timeit(ref_call, number=10, globals=globals())
        print(f"Time for {name} normalized by {ref_fun}: {t/t0:.2f}")


def compute_svps(f, T):
    for g in f:
        svp = f[g]['fun'](T)
        f[g]['svp'] = svp


def compute_mean(f):
    ref_vals = []
    for g in f:
        if f[g]['ref']:
            ref_vals.append(f[g]['svp'].m)
    ref_vals = np.asarray(ref_vals)
    res = np.mean(ref_vals, axis=0)

    return res * units.millibar
        

def table_diagnostics(f):
    t = np.arange(-40,60,10) * units.celsius

    print('                   Saturation Vapor Pressure Table [mb]')
    print('                                Temperature [C]')
    print('   Method    |  -40  |  -30  |  -20  |  -10  |   0   |   10  |'
          '   20  |   30  |   40  |   50')
    print(94*'-')

    svp_list = []
    for name, fun in f.items():
        if fun['ref']:
            svp = fun['fun'](t).m
            svp_list.append(svp)
            print(f'{name:^13}', end='')
            for s in svp:
                print(f'|{s:7.4f}', end='')
            print()
    svp_mean = np.mean(np.asarray(svp_list), axis=0)
    print('    Mean     ', end='')
    for s in svp_mean:
        print(f'|{s:7.4f}', end='')
    print()

    for name, fun in f.items():
        if not fun['ref']:
            svp = fun['fun'](t).m
            svp_list.append(svp)
            print(f'{name:^13}', end='')
            for s in svp:
                print(f'|{s:7.4f}', end='')
            print()

    print()
    print('       Errors [mb]')
    print('   Method    |  -40  |  -30  |  -20  |  -10  |   0   |   10  |'
          '   20  |   30  |   40  |   50')
    print(94*'-')

    for name, fun in f.items():
        svp = fun['fun'](t).m
        print(f'{name:^13}', end='')
        for i, s in enumerate(svp):
            print(f'|{(s-svp_mean[i]):7.4f}', end='')
        print()

# Store info on the various functions
functions = {'Goff': {'fun': goff65, 'ref': True},
             'Wexler': {'fun': wexler76, 'ref': True},
             'Murphy Koop': {'fun': murphykoop05, 'ref': True},
             'Wagner Pruss': {'fun': wagnerpruss02, 'ref': True},
             'Flatau': {'fun': flatau6, 'ref': False},
             'MetPy': {'fun': saturation_vapor_pressure, 'ref': False},
             'Koutsoyiannis': {'fun': koutsoyiannis, 'ref': False},
             'Romps': {'fun': romps17, 'ref': False},
             'Alduchov': {'fun': alduchov96, 'ref': False},
             'Gueymard': {'fun': gueymard7, 'ref': False},
             'Huang': {'fun': huang18, 'ref': False}}

Ts = np.linspace(-55,55,4000000) * units.celsius

calc_efficiency(functions)

compute_svps(functions, Ts)
#for f in functions:
#    print(functions[f]['svp'][999999])

ref_svp = compute_mean(functions)
#print(ref_svp[999999])

ref_funs = [f for f in functions if functions[f]['ref']]
compare_abs(functions, ref_funs, Ts, ref_svp)
compare_rel(functions, ref_funs, Ts, ref_svp)

apx_funs = [f for f in functions if not functions[f]['ref']]
compare_abs(functions, apx_funs, Ts, ref_svp)
compare_rel(functions, apx_funs, Ts, ref_svp)

table_diagnostics(functions)