Compare revisions

bjorn-stevens · bjorn-stevens · bjorn-stevens · bjorn-stevens · bjorn-stevens · bjorn-stevens
--- a/examples/examples.ipynb
+++ b/examples/examples.ipynb
--- a/examples/saturation-water-vapor.ipynb
+++ b/examples/saturation-water-vapor.ipynb
 %% Cell type:markdown id: tags:

 # Notes on calculations of vapor pressures and the specific heats #

 These notes do not use the moist_thermodynamic libararies, rather they were the basis for the choice for the particular formulations for the specific heats and saturation vapor pressures within that library.

 There are a large number of expressions for the saturation and sublimation vapor pressure in the literature, and many of these, even recent ones, seem to reference previous studies in a haphazard way.  So how much do these differ, is there a standard, and by what criteria should one judge them by.  Here I try to develop an intuition for the answers.

 The first thing to note is that there is a community that concerns itself with this question.  They call themselves the international association for the physical properties of water and steam, and mostly concern themselves with the behavior of water at high temperature.  The approach of the IAPWS is to develop an empirical equation of state for water, in the form of a specification of its Helmholtz free energy, or potential, from which all other properties can be derived.  The standard reference for the IAPWS equation of state is the publication by Wagner and Pru{\ss} (Thermodynamic Properties of Ordinary Water) published in 2002 and which describes the IAPWS-95 approved formulation.  Minor corrections have since been made to this, which as best I can tell are relevant at high temperatures.  The most substantial change has been the TEOS-10 work by Rainer Feistel of IOW, which extends these approaches to composite systems, thereby allowing for representations of sea-water and moist air.  By working with an equation of state, all properties of water, from the specific heats to the gas constants to the phase-change enthalpies can be derived consistently.  The disadvantage of this approach is that the equation is derived by positing an analytic form that is then fit to a very wide and diverse abundance of existing data.  The resultant equation is described in an ideal part, which involves a summation of nine terms and thirteen coefficients, and a residual part, with more than 50 terms and over 200 constants.

 For the case of the saturation vapor pressure over water Wagner and Pru{\ss} suggest a much simpler equation that is described in terms of only six coefficients. First, below I compare the relative error to the IAPWS standard as has been formlated and distributed in the iapws python package, version (1.4).  There has been some discussion on the web of its implementation, but the similarity with the Wagner and Pru{\ss} formulation gives me confidence.  Next we look at the TEOS-10 Sea-Ice-Air formulations of Feistel et al (Ocean Sci., 6, 91–141, 2010) which are distributed as FORTRAN90 code which I downloaded, ran, and tabulated to assess some empirical fits later used as a reference.

 %% Cell type:code id: tags:

 ``` python
 import numpy as np
 import matplotlib.pyplot as plt
 import seaborn as sns
 import os
 from scipy import interpolate, optimize

 os.makedirs("plots", exist_ok=True)

 !%matplotlib inline
 ```

 %% Output

    zsh:fg:1: no job control in this shell.

 %% Cell type:code id: tags:

 ``` python
 gravity = 9.8076

 cpd     = 1006.
 Rd      = 287.05

 Rv      = 461.53    # IAPWS97 at 273.15
 cpv     = 1865.01   # ''
 lv0     = 2500.93e3 # IAPWS97 at 273.15
 lf0     =  333.42e3 #''

 cl      = 4179.57   # IAPWS97 at 305 and P=0.1 MPa (chosen to give a good fit for es over ice)
 ci      = 1905.43   # IAPWS97 at 247.065 and P=0.1 MPa (chosen to give a good fit for es over ice)

 eps1     = Rd/Rv
 eps2     = Rv/Rd -1.

 P0      = 100000.  # Standard Pressure
 T0      = 273.15   # Standard Temperature
 PvC     = 22.064e6 # Critical pressure of water vapor
 TvC     = 647.096  # Critical temperature of water vapor
 TvT     = 273.16   # Triple point temperature of water
 PvT     = 611.655
 lvT     = lv0 + (cpv-cl)*(TvT-T0)
 lfT     = lf0 + (cpv-ci)*(TvT-T0)
 lsT     = lvT + lfT

 es_default = 'sonntag'

 def thermo_input(x, xtype='none'):

    import numpy as np

    x = np.asarray(x).flatten()
    scalar_input = False
    if x.ndim == 0:
        x = x[None]  # Makes x 1D
        scalar_input = True

    if (xtype == 'Kelvin' and x.max() < 100 ): x = x+273.15
    if (xtype == 'Celcius'and x.max() > 100 ): x = x-273.15
    if (xtype == 'Pascal' and x.max() < 1200): x = x*100.
    if (xtype == 'kg/kg'  and x.max() > 1.0) : x = x/1000.
    if (xtype == 'meter'  and x.max() < 10.0): print('Warning: input should be in meters, max value less than 10, not corrected')

    return x, scalar_input

 def eslf(T, formula=es_default):
    """ Returns the saturation vapour pressure [Pa] over liquid given
 	the temperature.  Temperatures can be in Celcius or Kelvin.
 	Formulas supported are
 	  - Goff-Gratch (1994 Smithsonian Tables)
 	  - Sonntag (1994)
 	  - Flatau
 	  - Magnus Tetens (MT)
      - Romps (2017)
      - Murphy-Koop
      - Bolton
      - Wagner and Pruss (WP, 2002) is the default
 	>>> eslf(273.16)
 	611.657
    """
    import numpy as np

    x,  scalar_input = thermo_input(T, 'Kelvin')

    if formula == "flatau":
        if (np.min(x) > 100): x = x-273.16
        np.maximum(x,-80.)
        c_es= np.asarray([0.6105851e+03, 0.4440316e+02, 0.1430341e+01, 0.2641412e-01,
                           0.2995057e-03,0.2031998e-05,0.6936113e-08,0.2564861e-11,-0.3704404e-13])
        es = np.polyval(c_es[::-1],x)
    elif formula == "bolton":
        if (np.min(x) > 100): x = x-273.15
        es = 611.2*np.exp((17.67*x)/(243.5+x))
    elif formula == "sonntag":
        xx = -6096.9385/x + 16.635794 - 2.711193e-2*x + 1.673952e-5*x*x + 2.433502 * np.log(x)
        es = 100.*np.exp(xx)
    elif formula =='goff-gratch':
        x1 = 273.16/x
        x2 = 373.16/x
        xl = np.log10(1013.246 ) - 7.90298*(x2 - 1) + 5.02808*np.log10(x2) - 1.3816e-7*(10**(11.344*(1.-1./x2)) - 1.0) + 8.1328e-3 * (10**(-3.49149*(x2-1)) - 1.0)
        es =10**(xl+2) # plus 2 converts from hPa to Pa
    elif formula == 'wagner-pruss':
        vt = 1.-x/TvC
        es = PvC * np.exp(TvC/x * (-7.85951783*vt + 1.84408259*vt**1.5 - 11.7866497*vt**3 + 22.6807411*vt**3.5 - 15.9618719*vt**4 + 1.80122502*vt**7.5))
    elif formula == 'hardy98':
        y  = -2.8365744e+3/(x*x) - 6.028076559e+3/x + 19.54263612 - 2.737830188e-2*x + 1.6261698e-5*x**2 + 7.0229056e-10*x**3 - 1.8680009e-13*x**4 + 2.7150305 * np.log(x)
        es = np.exp(y)
    elif formula == 'romps':
        Rr    = 461.
        cvl_r = 4119
        cvv_r = 1418
        cpv_r = cvv_r + Rr
        es = 611.65 * (x/TvT) **((cpv_r-cvl_r)/Rr) * np.exp((2.37403e6 - (cvv_r-cvl_r)*TvT)*(1/TvT - 1/x)/Rr)
    elif formula == "murphy-koop":
        es = np.exp(54.842763 - 6763.22/x - 4.210*np.log(x) + 0.000367*x + np.tanh(0.0415*(x - 218.8)) * (53.878 - 1331.22/x - 9.44523 * np.log(x) + 0.014025*x))
    elif formula == "standard-analytic":
        c1 = (cpv-cl)/Rv
        c2 = lvT/(Rv*TvT) - c1
        es = PvT * np.exp(c2*(1.-TvT/x)) * (x/TvT)**c1
    else:
        exit("formula not supported")

    es = np.maximum(es,0)
    if scalar_input:
        return np.squeeze(es)
    return es

 def esif(T, formula=es_default):
    """ Returns the saturation vapour pressure [Pa] over ice given
 	the temperature.  Temperatures can be in Celcius or Kelvin.
 	uses the Goff-Gratch (1994 Smithsonian Tables) formula
 	>>> esli(273.15)
 	6.112
 m    """
    import numpy as np

    x,  scalar_input = thermo_input(T, 'Kelvin')

    if formula == "sonntag":
        es = 100 * np.exp(24.7219 - 6024.5282/x + 0.010613868*x - 0.000013198825*x**2 - 0.49382577*np.log(x))
    elif formula == "goff-gratch":
        x1 = 273.16/x
        xi = np.log10(   6.1071) - 9.09718*(x1 - 1) - 3.56654*np.log10(x1) + 0.876793*(1 - 1./x1)
        es = 10**(xi+2)
    elif formula == "wagner-pruss": #(actually wagner et al, 2011)
        a1 = -0.212144006e+2
        a2 =  0.273203819e+2
        a3 = -0.610598130e+1
        b1 =  0.333333333e-2
        b2 =  0.120666667e+1
        b3 =  0.170333333e+1
        theta = T/TvT
        es = PvT * np.exp((a1*theta**b1 + a2 * theta**b2 + a3 * theta**b3)/theta)
    elif formula == "murphy-koop":
        es = np.exp(9.550426 - 5723.265/x + 3.53068 * np.log(x) - 0.00728332*x)
    elif formula == "romps":
        Rr    = 461.
        cvv_r = 1418.
        cvs_r = 1861.
        cpv_r = cvv_r + Rr
        es = 611.65 * (x/TvT) **((cpv_r-cvs_r)/Rr) * np.exp((2.37403e6 + 0.33373e6 - (cvv_r-cvs_r)*TvT)*(1/TvT - 1/x)/Rr)
    elif formula == "standard-analytic":
        c1 = (cpv-ci)/Rv
        c2 = lsT/(Rv*TvT) - c1
        es = PvT * np.exp(c2*(1.-TvT/x)) * (x/TvT)**c1
    else:
        exit("formula not supported")

    es = np.maximum(es,0)
    if scalar_input:
        return np.squeeze(es)
    return es

 def esilf(T,formula=es_default):
    import numpy as np
    return np.minimum(esif(T,formula),eslf(T,formula))

 def es(T,formula=es_default,state='liq'):

    import numpy as np

    x,  scalar_input = thermo_input(T, 'Kelvin')

    if (state == 'liq'):
        return eslf(x,formula)
    if (state == 'ice'):
        return esif(x,formula)
    if (state == 'mxd'):
        return esilf(x,formula)

 def des(T,formula=es_default,state='liq'):

    import numpy as np

    x,  scalar_input = thermo_input(T, 'Kelvin')
    dx = 0.01; xp = x+dx/2; xm = x-dx/2
    return (es(xp,formula,state)-es(xm,formula,state))/dx

 def dlnesdlnT(T,formula=es_default,state='liq'):

    import numpy as np

    x,  scalar_input = thermo_input(T, 'Kelvin')
    dx = 0.01; xp = x+dx/2; xm = x-dx/2
    return ((es(xp,formula,state)-es(xm,formula,state))/es(x,formula,state) * (x/dx))

 def phase_change_enthalpy(Tx,fusion=False):
    """ Returns the enthlapy [J/g] of vaporization (default) of water vapor or
    (if fusion=True) the fusion anthalpy.  Input temperature can be in degC or Kelvin
    >>> phase_change_enthalpy(273.15)
    2500.8e3
    """
    import numpy as np

    TC, scalar_input = thermo_input(Tx, 'Celcius')
    TK, scalar_input = thermo_input(Tx, 'Kelvin')

    if (fusion):
        el = lfT + (cl-ci)*(TK-TvT)
    else:
        el = lvT + (cpv-cl)*(TK-TvT)

    if scalar_input:
        return np.squeeze(el)
    return el

 ```

 %% Cell type:markdown id: tags:

 ## Selected properties of IAWPS water ##

 These routines only provide information on $c_{p,\mathrm{liq}}$ to a temperature of 253 K or -20$^\circ$C, but already demonstrate its divergent behavior (exponential increase) with increased super-cooling.  They also demonstrate the near linearity of $c_{p,\mathrm{ice}}.$

 %% Cell type:code id: tags:

 ``` python
 import iapws

 print ('Using IAPWS Version %s\n'%(iapws.__version__,))
 T = np.arange(183.15,313.15)
 ci_iapws = np.full(len(T),np.nan)
 cl_iapws = np.full(len(T),np.nan)
 for i,Tx in enumerate(T):
    if (Tx < 283): ci_iapws[i] =  iapws._iapws._Ice(Tx, 0.1)['cp']*1000 / ci
    if (Tx > 253.15): cl_iapws[i] =  iapws._iapws._Liquid(Tx, 0.1)['cp']*1000 / cl

 fig = plt.figure(figsize=(4,3))

 ax1 = plt.subplot(1,1,1)
 ax1.set_xlabel('$T$ / K')
 ax1.set_ylabel('$c_\mathrm{i}$ / %5.2f, $c_\mathrm{l}$ / %5.2f'%(ci,cl))
 ax1.set_xticks([185,247.07,273.15,305.00])
 plt.scatter([247.065],[1.])
 plt.scatter([305.000],[1.])
 plt.plot(T,ci_iapws)
 plt.plot(T,cl_iapws)

 sns.set_context("paper", font_scale=1.2)
 sns.despine(offset=10)

 plt.tight_layout()
 fig.savefig(plot_dir+'cp-Tdependance.pdf')

 TK = np.arange(273.15,315.15,0.01)
 es_iapws = np.zeros(len(TK))
 for i, x in enumerate(TK):
    es_iapws[i] = iapws.iapws97._PSat_T(x) *1.e6 #Temperature, [K]; Returns:Pressure, [MPa]
 ```

 %% Output

    Using IAPWS Version 1.5.2
    

    /Users/m219063/opt/miniforge3/lib/python3.9/site-packages/iapws/_iapws.py:124: UserWarning: Metastable ice
      warnings.warn("Metastable ice")

    ---------------------------------------------------------------------------
    NameError                                 Traceback (most recent call last)
 Input     In [3], in <cell line: 26>()
         23 sns.despine(offset=10)
         25 plt.tight_layout()
    ---> 26 fig.savefig(plot_dir+'cp-Tdependance.pdf')
         28 TK = np.arange(273.15,315.15,0.01)
         29 es_iapws = np.zeros(len(TK))
    NameError: name 'plot_dir' is not defined



 %% Cell type:markdown id: tags:

 ## Comparison with Sea-Air-Ice library of Feistel et al (2010) v4.0.1

 Here we compare different thermodynamic constants or empirical formula to the IAPWS-10 and TEOS-10 standards taken from the Sea-Air-Ice library.  These are calculated based on fits to potential functions as described above.  The libraries are run off-line and the output is tabulated for comparison. The Feistel et al formulation extends to the IAPWS formulation shown above to allow for representations of liquid to the homogeneous freezing point of ice.

 %% Cell type:code id: tags:

 ``` python
 import pandas as pd
 sia_dir = './data/'
 d0    = pd.read_csv(sia_dir + 'psat.dat',sep=' ',names=['T','Psat_liq','Psat_ice'])
 x0    = d0.to_xarray()
 d1    = pd.read_csv(sia_dir + 'cp-liq-vap.dat',sep=' ',names=['T','liq_density','vap_density','cp_liq','cp_vap','lv'])
 x1    = d1.to_xarray()
 d2    = pd.read_csv(sia_dir + 'cp-ice-vap.dat',sep=' ',names=['T','liq_density','vap_density','cp_ice','cp_vap','ls'])
 x2    = d2.to_xarray()

 fig = plt.figure(figsize=(7,14/1.610834))

 ax1 = plt.subplot(2,1,2)
 ax1.set_ylabel('difference / %')
 ax1.set_xlabel('T / K')

 formula = 'wagner-pruss'
 es_r = es(x0.T,formula=formula,state='ice')
 diff = es_r/x0.Psat_ice
 plt.plot(x0.T,100*(diff.where(diff>1.e-4)-1),c='dodgerblue',label=formula)
 es_r = es(x0.T,formula=formula,state='liq')
 diff = es_r/x0.Psat_liq
 plt.plot(x0.T,100*(diff.where(diff>1.e-4)-1),c='orange',label=formula)

 formula = 'romps'
 es_r = es(x0.T,formula=formula,state='ice')
 diff = es_r/x0.Psat_ice
 plt.plot(x0.T,100*(diff.where(diff>1.e-4)-1),c='dodgerblue',ls='dashed',label=formula)
 es_r = es(x0.T,formula=formula,state='liq')
 diff = es_r/x0.Psat_liq
 plt.plot(x0.T,100*(diff.where(diff>1.e-4)-1),c='orange',ls='dashed',label=formula)

 formula = 'murphy-koop'
 es_r = es(x0.T,formula=formula,state='ice')
 diff = es_r/x0.Psat_ice
 plt.plot(x0.T,100*(diff.where(diff>1.e-4)-1),c='dodgerblue',ls='dotted',label=formula)
 es_r = es(x0.T,formula=formula,state='liq')
 diff = es_r/x0.Psat_liq
 plt.plot(x0.T,100*(diff.where(diff>1.e-4)-1),c='orange',ls='dotted',label=formula)
 plt.legend(loc="best",ncol=3)

 ###
 ax2 = plt.subplot(2,1,1)
 ax2.set_ylabel('difference / %')
 ax2.set_xlabel('T / K')
 ax2.set_ylim(0.8,1.2)

 plt.plot(x1.T,x1.cp_vap/cpv,c='green',ls='solid',label='$c_{p,\mathrm{vap}}$ (vap-liq)')
 plt.plot(x2.T,x2.cp_vap/cpv,c='green',ls='dashed',label='$c_{p,\mathrm{vap}}$ (vap-ice)')

 plt.plot(x1.T,x1.cp_liq/cl,c='orange',ls='solid',label='$c_{p,\mathrm{liq}}$')
 plt.plot(x2.T,x2.cp_ice/ci,c='dodgerblue',ls='solid',label='$c_{p,\mathrm{ice}}$')

 lvx = phase_change_enthalpy(x1.T,fusion=False)
 lsx = phase_change_enthalpy(x2.T,fusion=True) + phase_change_enthalpy(x2.T,fusion=False)

 plt.plot(x1.T,x1.lv/lvx,c='gray',ls='solid',label='$\\ell_v$')
 plt.plot(x2.T,x2.ls/lsx,c='gray',ls='dashed',label='$\\ell_s$')

 plt.legend(loc="lower right",ncol=3)
 sns.set_context("paper")
 sns.despine(offset=10)
 plt.tight_layout()
 ```

+%% Output
+
+
+
 %% Cell type:markdown id: tags:

 ## Comparing formulations of saturation vapor pressure ##

 This comparison of relative error suggests that the Wagner-Pru{\ss}, Murphy and Koop, Hardy, and Sonntag formulations lie closest to the IAPWS-97 reference.  Romps (2017) and Bolton (1980) are similarly accurate and may have advantages.  Hardy is interesting as it appears in a technical document and is rarely mentioned in the subsequent literature, but used by Vaisala in the calibration of their sondes

 %% Cell type:markdown id: tags:

 ### Temperatures above the triple point ###

 %% Cell type:code id: tags:

 ``` python
 state = 'liq'
 fig = plt.figure(figsize=(10,5))
 ax1 = plt.subplot(1,1,1)
 ax1.set_xlabel('$T$ / K')
 ax1.set_ylabel('$e_{\mathrm{s,x}}/e_{\mathrm{s,ref}} - 1$')
 ax1.set_yscale('log')

 es_ref = es_iapws
 es_w = es(TK,formula="wagner-pruss",state=state)
 es_r = es(TK,formula='romps',state=state)
 es_g = es(TK,formula='goff-gratch',state=state)
 es_m = es(TK,formula='murphy-koop',state=state)
 es_s = es(TK,formula='sonntag',state=state)
 es_b = es(TK,formula='bolton',state=state)
 es_f = es(TK,formula='flatau',state=state)
 es_h = es(TK,formula='hardy98',state=state)
 es_a = es(TK,formula='standard-analytic',state=state)

 plt.plot(TK,np.abs(es_h/es_ref-1),c='tab:blue',ls='solid',label='Hardy (1998)')
 plt.plot(TK,np.abs(es_f/es_ref-1),c='tab:orange',label='Flatau (1992)')
 plt.plot(TK,np.abs(es_g/es_ref-1),c='tab:green',label='Goff-Gratch (1957)')
 plt.plot(TK,np.abs(es_b/es_ref-1),c='tab:red',ls='dotted',label='Bolton (1980)')

 plt.plot(TK,np.abs(es_r/es_ref-1),c='tab:purple',label='Romps (2017)')
 plt.plot(TK,np.abs(es_s/es_ref-1),c='tab:grey',label='Sonntag (1990)')
 plt.plot(TK,np.abs(es_m/es_ref-1),c='tab:pink',label='Murphy-Koop (2005)')
 plt.plot(TK,np.abs(es_w/es_ref-1),c='tab:brown',label='Wagner-Pruss (2002)')
 plt.plot(TK,np.abs(es_a/es_ref-1),c='tab:purple',ls='dotted',label='Analytic')

 plt.legend(loc="lower right",ncol=3)

 sns.set_context("paper", font_scale=1.2)
 sns.despine(offset=10)

 fig.savefig(plot_dir+'es_l-error.pdf')
 ```

+%% Output
+
+    ---------------------------------------------------------------------------
+    NameError                                 Traceback (most recent call last)
+Input     In [5], in <cell line: 8>()
+          5 ax1.set_ylabel('$e_{\mathrm{s,x}}/e_{\mathrm{s,ref}} - 1$')
+          6 ax1.set_yscale('log')
+    ----> 8 es_ref = es_iapws
+          9 es_w = es(TK,formula="wagner-pruss",state=state)
+         10 es_r = es(TK,formula='romps',state=state)
+    NameError: name 'es_iapws' is not defined
+
+
+
 %% Cell type:markdown id: tags:

 ### Extension to the freezing regime ###

 To extend over the entire temperature range a different reference is required, for this any of the Hardy, Sonntag, Murphy-Koop and Wagner-Pru{\ss} formulations could suffice.  We choose Wagner-Pru{\ss} because Wagner's group is responsible for the standard, and has also developed the IAPWS standard for saturation vapor pressure over ice.  Below the results are plooted with respect to this standard over a much larger temperature range.

 It is not clear how accurate Wagner and Pru{\ss} wis hen extended well beyond the IAPWS range, based on which it might be that the grouping of errors of similar magnitude from the Bolton, Flatau and Goff-Gratch formulations are indicative of a low temperature bias in the Wagner-Pru{\ss} formualtion.  I doubt that this is the case, as the poor performance of all these formulations in the higher temperature range, and the simplicity of their formulation make it unlikely.   The agreement of the Murphy-Koop formulation with these simpler formulations at low temperature may be indicative of Murphy and Koops focus on saturation over ice rather than liquid.

 %% Cell type:code id: tags:

 ``` python
 state = 'liq'
 fig = plt.figure(figsize=(10,5))
 ax1 = plt.subplot(1,1,1)
 ax1.set_xlabel('$T$ / K')
 ax1.set_ylabel('$e_{\mathrm{s,x}}/e_{\mathrm{s,ref}} - 1$')
 ax1.set_yscale('log')

 TK = np.arange(180,320,0.5)

 es_w = es(TK,formula="wagner-pruss",state=state)
 es_r = es(TK,formula='romps',state=state)
 es_g = es(TK,formula='goff-gratch',state=state)
 es_m = es(TK,formula='murphy-koop',state=state)
 es_s = es(TK,formula='sonntag',state=state)
 es_b = es(TK,formula='bolton',state=state)
 es_f = es(TK,formula='flatau',state=state)
 es_h = es(TK,formula='hardy98',state=state)
 es_a = es(TK,formula='standard-analytic',state=state)

 es_ref = es_w

 plt.plot(TK,np.abs(es_h/es_ref-1),c='tab:blue',ls='solid',label='Hardy (1998)')
 plt.plot(TK,np.abs(es_f/es_ref-1),c='tab:orange',label='Flatau (1992)')
 plt.plot(TK,np.abs(es_g/es_ref-1),c='tab:green',label='Goff-Gratch (1957)')
 plt.plot(TK,np.abs(es_b/es_ref-1),c='tab:red',ls='dotted',label='Bolton (1980)')

 plt.plot(TK,np.abs(es_r/es_ref-1),c='tab:purple',label='Romps (2017)')
 plt.plot(TK,np.abs(es_s/es_ref-1),c='tab:grey',label='Sonntag (1990)')
 plt.plot(TK,np.abs(es_m/es_ref-1),c='tab:pink',label='Murphy-Koop (2005)')
 plt.plot(TK,np.abs(es_a/es_ref-1),c='tab:purple',ls='dotted',label='Analytic')

 #plt.plot(TK,np.abs(es_w/es_ref-1),c='tab:olive',label='Wagner-Pruss (2002)')

 plt.legend(loc="lower left",ncol=2)

 sns.set_context("paper", font_scale=1.2)
 sns.despine(offset=10)

 fig.savefig(plot_dir+'es_lsc-error.pdf')
 ```

 %% Cell type:markdown id: tags:

 ## Sublimation vapor pressure ##

 A subset of the formulations also postulate forms for the saturation vapor pressure over ice.  For the reference in this quantity we use Wagner et al., (2011) as this has been adopted as the IAPWS standard.   Here is seems that Murphy and Koop's (2005) formulation behaves very well in comparision to Wagner et al., but Sonntag is also quite adequate, particularly at lower ($T<273.15$ K) temperatures where it is likely to be applied.

 %% Cell type:code id: tags:

 ``` python
 state = 'ice'
 fig = plt.figure(figsize=(10,5))
 ax1 = plt.subplot(1,1,1)
 ax1.set_xlabel('$T$ / K')
 ax1.set_ylabel('$e_{\mathrm{s,x}}/e_{\mathrm{s,ref}} - 1$')
 ax1.set_yscale('log')

 TK = np.arange(180,320,0.5)

 es_w = es(TK,formula="wagner-pruss",state=state)
 es_r = es(TK,formula='romps',state=state)
 es_g = es(TK,formula='goff-gratch',state=state)
 es_m = es(TK,formula='murphy-koop',state=state)
 es_s = es(TK,formula='sonntag',state=state)
 es_a = es(TK,formula='standard-analytic',state=state)
 es_ref = es_w

 plt.plot(TK,np.abs(es_g/es_ref-1),c='tab:green',label='Goff-Gratch (1957)')
 plt.plot(TK,np.abs(es_r/es_ref-1),c='tab:purple',label='Romps (2017)')
 plt.plot(TK,np.abs(es_s/es_ref-1),c='tab:grey',label='Sonntag (1990)')
 plt.plot(TK,np.abs(es_m/es_ref-1),c='tab:pink',label='Murphy-Koop (2005)')
 plt.plot(TK,np.abs(es_a/es_ref-1),c='tab:purple',ls='dotted',label='Analytic')

 #plt.plot(TK,np.abs(es_w/es_ref-1),c='tab:olive',label='Wagner-Pruss (2002)')

 plt.legend(loc="lower left",ncol=2)

 sns.set_context("paper", font_scale=1.2)
 sns.despine(offset=10)

 fig.savefig(plot_dir+'es_i-error.pdf')
 ```

 %% Cell type:markdown id: tags:

 ## Clausius Clapeyron ##

 Often over looked is that many conceptual models are built on the application of the Clausius-Clapeyron equation,
 \begin{equation}
 \frac{\mathrm{d} \ln e_\mathrm{s}}{\mathrm{d \ln T}} \left(\frac{\ell_\mathrm{v}}{R_\mathrm{v} T}\right)^{-1} = 1
 \end{equation}
 with the assumption that the vaporization enthalpy, $\ell_\mathrm{v}$ that appears in this equation, is linear in temperature following Kirchoff's relation.  This is similar to assuming that the specific heats are independent of temeprature, an idealization which is, unfortunately, just that, and idealization.

 But because of this it is interesting to compare this expression as given by the above formulation of the saturation vapor pressure (through their numerical derivative) and independent expressions of $\ell_\mathrm{v}$ based on the assumption of constant specific heats.

 This is shown below for ice and liquid saturation.  The analytic expression, which has larger errors for es is constructued to satisfy this relationship and is exact to the precision of the numerical calculations.  The various formulations using more accurate expressions for $e_s$ which implicityl don't assume constancy in specific heats are similarly accurate, with the exception of Goff-Gratch, and Romps for Ice.  Hardy is only shown for water.  For ice Sonntag does not behave well for $T> 290$ K, but it is not likely to be used at these temperatures.  Note that Romps would be perfect had we adopted his modified specific heats.

 Based on the above my recommendation is to use the formulations by Wagner's group, unless one is interested in very low temperatures ($T<180$K) in which case the formulation of Koop and Murphy may be desirable.  For just liquid processes Hardy might be a good choice, it is less well known but used by Vaisala for its sondes.   There may be advantages to using Sonntag if there is interest in liquid and ice as it might allow more efficient implementations, but for my tests all formulations were within 30% of one another.

 Another alternative, would be to use the analytic approach, either using Romps' formulae if getthing the staturation vapor pressure as close to measurements as possible is preferred, or using the analytic formula with the correct (at the standard temperature and pressure) specific heats and gast constants.

 %% Cell type:code id: tags:

 ``` python
 state = 'liq'

 fig = plt.figure(figsize=(10,10))

 ax1 = plt.subplot(2,1,1)

 ax1.set_ylabel('$|\mathrm{CC}_\mathrm{liq} - 1|$')
 ax1.set_yscale('log')
 ax1.set_xticklabels([])

 TK = np.arange(180,320,0.5)

 lv   = phase_change_enthalpy(TK)
 if (state == 'ice'): lv += phase_change_enthalpy(TK,fusion=True)

 y    = lv/(Rv * TK)
 cc_w = dlnesdlnT(TK,formula="wagner-pruss",state=state) / y
 cc_r = dlnesdlnT(TK,formula='romps',state=state) /y
 cc_g = dlnesdlnT(TK,formula='goff-gratch',state=state) /y
 cc_m = dlnesdlnT(TK,formula='murphy-koop',state=state) /y
 cc_s = dlnesdlnT(TK,formula='sonntag',state=state) /y
 cc_h = dlnesdlnT(TK,formula='hardy98',state=state) /y
 cc_a = dlnesdlnT(TK,formula='standard-analytic',state=state) /y

 plt.plot(TK,np.abs(cc_h/1 -1.),c='tab:blue',label='Hardy (1998)')
 plt.plot(TK,np.abs(cc_g/1 -1.),c='tab:green',label='Goff-Gratch (1957)')
 plt.plot(TK,np.abs(cc_r/1 -1.),c='tab:purple',label='Romps (2017)')
 plt.plot(TK,np.abs(cc_s/1 -1.),c='tab:grey',label='Sonntag (1990)')
 plt.plot(TK,np.abs(cc_m/1 -1.),c='tab:pink',label='Murphy-Koop (2005)')
 plt.plot(TK,np.abs(cc_w/1 -1.),c='tab:olive',label='Wagner-Pruss (2002)')
 plt.plot(TK,np.abs(cc_a/1 -1.),c='tab:purple',ls='dotted',label='Analytic')

 plt.legend(loc="lower left",ncol=1)

 state = 'ice'
 TK = np.arange(180,320,0.5)

 lv   = phase_change_enthalpy(TK)
 if (state == 'ice'): lv   = phase_change_enthalpy(TK,fusion=True) + phase_change_enthalpy(TK)

 y    = lv/(Rv * TK)
 cc_w = dlnesdlnT(TK,formula="wagner-pruss",state=state) / y
 cc_r = dlnesdlnT(TK,formula='romps',state=state) /y
 cc_g = dlnesdlnT(TK,formula='goff-gratch',state=state) /y
 cc_m = dlnesdlnT(TK,formula='murphy-koop',state=state) /y
 cc_s = dlnesdlnT(TK,formula='sonntag',state=state) /y
 cc_a = dlnesdlnT(TK,formula='standard-analytic',state=state) /y

 ax2 = plt.subplot(2,1,2)
 ax2.set_xlabel('$T$ / K')
 ax2.set_ylabel('$|\mathrm{CC}_\mathrm{ice} - 1|$')
 ax2.set_yscale('log')

 plt.plot(TK,np.abs(cc_g/1 -1.),c='tab:green',label='Goff-Gratch (1957)')
 plt.plot(TK,np.abs(cc_r/1 -1.),c='tab:purple',label='Romps (2017)')
 plt.plot(TK,np.abs(cc_s/1 -1.),c='tab:grey',label='Sonntag (1990)')
 plt.plot(TK,np.abs(cc_m/1 -1.),c='tab:pink',label='Murphy-Koop (2005)')
 plt.plot(TK,np.abs(cc_w/1 -1.),c='tab:olive',label='Wagner-Pruss (2002)')
 plt.plot(TK,np.abs(cc_a/1. -1.),c='tab:purple',ls='dotted',label='Analytic')

 sns.set_context("paper", font_scale=1.2)
 sns.despine(offset=10)

 fig.savefig(plot_dir+'cc-error.pdf')
 ```

 %% Cell type:markdown id: tags:

 ## Optimizing analytic fits for saturation and sublimation vapor pressure ##

 Romps suggests modifying the specific heats of liquid, ice and the gas constant of vapor to arrive at an optimal fit for the saturation vapor pressure using the analytic form.  One can do almost as good by just modifying the specific heat of the condensate phases.  Here we show how the maximum error in the fit depends on the specific heat of the condensate phases as compared to the reference, and how we arrive at our optimal fit by only manipulating the condensate phase specific heats to values that they anyway adopt within the range of temperatures spanned by the atmosphere.  This justifys the default choice for saturation vapor pressure and the specific heats used in aes_thermo.py

 %% Cell type:code id: tags:

 ``` python
 fig = plt.figure(figsize=(10,5))

 cl_1 = (iapws._iapws._Liquid(265, 0.1)['cp'])*1000.
 cl_2 = (iapws._iapws._Liquid(305, 0.1)['cp'])*1000
 ci_1 = (iapws._iapws._Ice(193, 0.01)['cp'])*1000.
 ci_2 = (iapws._iapws._Ice(273, 0.10)['cp'])*1000

 cls = np.arange(cl_2,cl_1)
 err = np.zeros(len(cls))

 ax1 = plt.subplot(1,2,1)
 ax1.set_xlabel('$c_\mathrm{liq}$ / Jkg$^{-1}$K$^{-1}$')
 ax1.set_ylabel('$(e_{\mathrm{s,x}}/e_{\mathrm{s,ref}} - 1)_\mathrm{max}$ / %')
 ax1.set_yscale('log')

 state = 'liq'
 TK = np.arange(260,300,0.5)
 es_ref = es(TK,formula="wagner-pruss",state=state)
 for i,cx in enumerate(cls):
    c1 = (cpv-cx)/Rv
    c2 = lvT/(Rv*TvT) - c1
    es_a = PvT * np.exp(c2*(1.-TvT/TK)) * (TK/TvT)**c1
    err[i] = np.max(np.abs(es_a/es_ref -1.))*100.
 ax1.plot(cls,err,c='tab:purple',ls='dotted',label='Analytic $c_\mathrm{liq}$ for $T\in$ (260K,305K)')
 ax1.legend(loc="upper left",ncol=2)

 cis = np.arange(ci_1,ci_2)
 err = np.zeros(len(cis))

 ax2 = plt.subplot(1,2,2)
 ax2.set_xlabel('$c_\mathrm{ice}$ / Jkg$^{-1}$K$^{-1}$')
 ax2.set_ylabel('$(e_{\mathrm{s,x}}/e_{\mathrm{s,ref}} - 1)_\mathrm{max}$ / %')
 ax2.set_yscale('log')

 state = 'ice'
 TK = np.arange(180,273,0.5)
 es_ref = es(TK,formula="wagner-pruss",state=state)
 for i,cx in enumerate(cis):
    c1 = (cpv-cx)/Rv
    c2 = lsT/(Rv*TvT) - c1
    es_a = PvT * np.exp(c2*(1.-TvT/TK)) * (TK/TvT)**c1
    err[i] = np.max(np.abs(es_a/es_ref -1.))*100.
 ax2.plot(cis,err,c='tab:purple',ls='dotted',label='Analytic $c_\mathrm{ice}$ for $T\in$ (193K,273K)')
 ax2.legend(loc="upper right",ncol=2)

 sns.set_context("paper", font_scale=1.2)
 sns.despine(offset=10)

 fig.savefig(plot_dir+'es-analytic-fits.pdf')
 Tfit = 305
 print ('Taking fit for $c_\mathrm{liq}=$ %3.2f J/(kg K) at $T=$ %3.2f K'%(iapws._iapws._Liquid(Tfit, 0.1)['cp']*1000.,Tfit))
 Tfit = 247.065
 print ('Taking fit for $c_\mathrm{ice}=$ %3.2f J/(kg K) at $T=$ %3.2f K'%(iapws._iapws._Ice(Tfit, 0.1)['cp']*1000.,Tfit))
 ```

 %% Cell type:markdown id: tags:

 ## RCEMIP comparision ##

 During RCEMIP (Wing et al.) different models output different RH, differing in ways of calculating it and also whether or not it was calculated relative to liquid or ice.  In this analysis we create a python implementation of the intial RCEMIP sounding and then for the given state estimate the RH using different formulat and different assumptions regarding the reference condensate (liquid/ice).  We also show the difference associated with 1 K of temperature.

 %% Cell type:code id: tags:

 ``` python
 def rcemip_on_z(z,SST):
    # function [T,q,p] = rcemip_on_z(z,SST)
    #
    # Inputs:
    # z: array of heights (low to high, m)
    # SST: sea surface temperature (K)
    #
    # Outputs:
    T = np.zeros(len(z)) # temperature (K)
    q = np.zeros(len(z)) # specific humidity (g/g)
    p = np.zeros(len(z)) # pressure (Pa)

    ## Constants
    g = 9.79764 #m/s^2
    Rd = 287.04 #J/kgK

    ## Parameters
    p0 = 101480   #Pa surface pressure
    qt = 10**(-11) #g/g specific humidity at tropopause
    zq1 = 4000 #m
    zq2 = 7500 #m
    zt = 15000 #m tropopause height
    gamma = 0.0067 #K/m lapse rate

    ## Scratch
    Tv = np.zeros(len(z)) # temperature (K)

    if SST == 295:
        q0 = 0.01200; #g/g specific humidity at surface (adjusted from 300K value so RH near surface approx 80%)
    elif SST == 300:
        q0 = 0.01865; #g/g specific humidity at surface
    elif SST == 305:
        q0 = 0.02400 #g/g specific humidity at surface (adjusted from 300K value so RH near surface approx 80%)

    T0 = SST - 0 #surface air temperature adjusted to be 0K less than SST

    ## Virtual Temperature at surface and tropopause
    Tv0 = T0*(1 + 0.608*q0) #virtual temperature at surface
    Tvt = Tv0 - gamma*zt #virtual temperature at tropopause z=zt

    ## Pressure
    pt = p0*(Tvt/Tv0)**(g/(Rd*gamma)); #pressure at tropopause z=zt
    p  = p0*((Tv0-gamma*z)/Tv0)**(g/(Rd*gamma)) #0 <= z <= zt
    p[z>zt] = pt*np.exp(-g*(z[z>zt]-zt)/(Rd*Tvt)) #z > zt

    ## Specific humidity
    q = q0*np.exp(-z/zq1)*np.exp(-(z/zq2)**2)
    q[z>zt] = qt #z > zt

    ## Temperature
    #Virtual Temperature
    Tv = Tv0 - gamma*z #0 <= z <= zt
    Tv[z>zt] = Tvt #z > zt

    #Absolute Temperature at all heights
    T = Tv/(1 + 0.608*q)

    return T, q, p

 z = np.arange(0,17000,100)
 T, q , p = rcemip_on_z(z,300)
 ```

 %% Cell type:code id: tags:

 ``` python
 def get_rh (T,q,p,formula='wagner-pruss',state='liq'):
    es_w = es(T,formula=formula,state=state)
    x = es_w * eps1/(p-es_w)
    return 100.*q*(1+x)/x

 fig = plt.figure(figsize=(4,5))

 ax1 = plt.subplot(1,1,1)
 ax1.set_ylabel('$z$ / km')
 ax1.set_xlabel('RH / %')
 ax1.set_ylim(0,14.5)
 ax1.set_yticks([0,4,8,12])

 plt.plot(get_rh(T,q,p,state='mxd'),z/1000.,label = 'Wagner Pruss (ice/liq)')
 plt.plot(get_rh(T+1,q,p,state='mxd'),z/1000.,label = 'Wagner Pruss (ice/liq) + 1 K')
 plt.plot(get_rh(T,q,p,state='ice'),z/1000.,label = 'Wagner Pruss (ice)')
 plt.plot(get_rh(T,q,p,formula='romps',state='mxd'),z/1000.,label = 'Romps (ice/liq)')
 plt.plot(get_rh(T,q,p),z/1000.,label = 'Wagner Pruss (liq)')
 plt.plot(get_rh(T,q,p,formula='flatau'),z/1000.,label = 'Flatau (liq)')

 plt.legend(loc="lower left",ncol=1)

 sns.set_context("paper")
 sns.despine(offset=10)
 plt.tight_layout()

 fig.savefig(plot_dir+'RCEMIP-RHerror.pdf')
 ```

 %% Cell type:markdown id: tags:

 ## lifting-condensation level approximations

 Here we compare the LCL base predictions to those proposed by Romps and Bolton as well as the differences between density potential temperatures.

 For the estimation of the LCL we modify the Romps expressions (using his code) to output pressure at the LCL, as this eliminates an assumption as to how pressure is distributed in the atmosphere, and thus only depends on the parcel state.  What we find is that the much simpler Bolton expression is as good as the more complex expression by Romps, and differences between the two are commensurate with those arising from slight differences in how the saturation vapor pressure is calculated.

 %% Cell type:code id: tags:

 ``` python
 # Version 1.0 released by David Romps on September 12, 2017.
 #
 # When using this code, please cite:
 #
 # @article{16lcl,
 #   Title   = {Exact expression for the lifting condensation level},
 #   Author  = {David M. Romps},
 #   Journal = {Journal of the Atmospheric Sciences},
 #   Year    = {2017},
 #   Volume  = {in press},
 # }
 #
 # This lcl function returns the height of the lifting condensation level
 # (LCL) in meters.  The inputs are:
 # - p in Pascals
 # - T in Kelvins
 # - Exactly one of rh, rhl, and rhs (dimensionless, from 0 to 1):
 #    * The value of rh is interpreted to be the relative humidity with
 #      respect to liquid water if T >= 273.15 K and with respect to ice if
 #      T < 273.15 K.
 #    * The value of rhl is interpreted to be the relative humidity with
 #      respect to liquid water
 #    * The value of rhs is interpreted to be the relative humidity with
 #      respect to ice
 # - ldl is an optional logical flag.  If true, the lifting deposition
 #   level (LDL) is returned instead of the LCL.
 # - min_lcl_ldl is an optional logical flag.  If true, the minimum of the
 #   LCL and LDL is returned.

 def lcl(p,T,rh=None,rhl=None,rhs=None,return_ldl=False,return_min_lcl_ldl=False):

    import math
    import scipy.special

    # Parameters
    Ttrip = 273.16     # K
    ptrip = 611.65     # Pa
    E0v   = 2.3740e6   # J/kg
    E0s   = 0.3337e6   # J/kg
    ggr   = 9.81       # m/s^2
    rgasa = 287.04     # J/kg/K
    rgasv = 461        # J/kg/K
    cva   = 719        # J/kg/K
    cvv   = 1418       # J/kg/K
    cvl   = 4119       # J/kg/K
    cvs   = 1861       # J/kg/K
    cpa   = cva + rgasa
    cpv   = cvv + rgasv

    # The saturation vapor pressure over liquid water
    def pvstarl(T):
        return ptrip * (T/Ttrip)**((cpv-cvl)/rgasv) * math.exp( (E0v - (cvv-cvl)*Ttrip) / rgasv * (1/Ttrip - 1/T) )
    # The saturation vapor pressure over solid ice
    def pvstars(T):
        return ptrip * (T/Ttrip)**((cpv-cvs)/rgasv) *  math.exp( (E0v + E0s - (cvv-cvs)*Ttrip) / rgasv * (1/Ttrip - 1/T))

    # Calculate pv from rh, rhl, or rhs
    rh_counter = 0
    if rh  is not None:
        rh_counter = rh_counter + 1
    if rhl is not None:
        rh_counter = rh_counter + 1
    if rhs is not None:
        rh_counter = rh_counter + 1
    if rh_counter != 1:
        print(rh_counter)
        exit('Error in lcl: Exactly one of rh, rhl, and rhs must be specified')
    if rh is not None:
    # The variable rh is assumed to be
    # with respect to liquid if T > Ttrip and
    # with respect to solid if T < Ttrip
        if T > Ttrip:
            pv = rh * pvstarl(T)
        else:
            pv = rh * pvstars(T)
            rhl = pv / pvstarl(T)
            rhs = pv / pvstars(T)
    elif rhl is not None:
        pv = rhl * pvstarl(T)
        rhs = pv / pvstars(T)
        if T > Ttrip:
            rh = rhl
        else:
            rh = rhs
    elif rhs is not None:
        pv = rhs * pvstars(T)
        rhl = pv / pvstarl(T)
        if T > Ttrip:
            rh = rhl
        else:
            rh = rhs
    if pv > p:
        return N

 # Calculate lcl_liquid and lcl_solid
    qv = rgasa*pv / (rgasv*p + (rgasa-rgasv)*pv)
    rgasm = (1-qv)*rgasa + qv*rgasv
    cpm = (1-qv)*cpa + qv*cpv
    if rh == 0:
        return cpm*T/ggr
    aL = -(cpv-cvl)/rgasv + cpm/rgasm
    bL = -(E0v-(cvv-cvl)*Ttrip)/(rgasv*T)
    cL = pv/pvstarl(T)*math.exp(-(E0v-(cvv-cvl)*Ttrip)/(rgasv*T))
    aS = -(cpv-cvs)/rgasv + cpm/rgasm
    bS = -(E0v+E0s-(cvv-cvs)*Ttrip)/(rgasv*T)
    cS = pv/pvstars(T)*math.exp(-(E0v+E0s-(cvv-cvs)*Ttrip)/(rgasv*T))
    X  = bL/(aL*scipy.special.lambertw(bL/aL*cL**(1/aL),-1).real)
    Y  = bS/(aS*scipy.special.lambertw(bS/aS*cS**(1/aS),-1).real)

    lcl = cpm*T/ggr*( 1 - X)
    ldl = cpm*T/ggr*( 1 - Y)

    # Modifications of the code to output Plcl or Pldl
    Plcl = PPa * X**(cpm/rgasm)
    Pldl = PPa * X**(cpm/rgasm)
    # Return either lcl or ldl
    if return_ldl and return_min_lcl_ldl:
        exit('return_ldl and return_min_lcl_ldl cannot both be true')
    elif return_ldl:
        return Pldl
    elif return_min_lcl_ldl:
        return min(Plcl,Pldl)
    else:
        return Plcl
 ```

 %% Cell type:code id: tags:

 ``` python
 import aes_thermo as mt
 PPa  = 101325.

 qt = np.arange(2.5e-3,8e-3,0.2e-3)
 TK = 285.
 Plcl_X = mt.plcl(TK,PPa,qt)
 Plcl_B = mt.plcl_boloton(TK,PPa,qt)
 Plcl_R = np.zeros(len(Plcl_X))

 for i,x in enumerate(qt):
    if (x>0.1): x = x/1000.
    RH        = mt.mr2pp(x/(1.-x),PPa)/mt.es(TK)
    Plcl_R[i] = lcl(PPa,TK,RH)

 del1 = (Plcl_B-Plcl_X)/100.
 del2 = (Plcl_R-Plcl_X)/100.

 fig = plt.figure(figsize=(10,5))
 ax1 = plt.subplot(1,2,1)
 ax1.set_ylabel('$P$ / hPa')
 ax1.set_xlabel('$q_\mathrm{t}$ / g/kg')
 #ax1.set_ylim(-1.2,1.2)

 plt.plot(qt*1.e3,del1,label='$\\delta_\mathrm{B}$, $T$=285K')
 plt.plot(qt*1.e3,del2,label='$\\delta_\mathrm{R}$, $T$=285K')
 #plt.gca().invert_yaxis()
 plt.legend(loc="best")

 qt = np.arange(0.5e-3,28e-3,0.2e-3)
 TK = 310.
 Plcl_X = mt.get_Plcl(TK,PPa,qt,iterate=True)
 Plcl_B = mt.get_Plcl(TK,PPa,qt)
 Plcl_R = np.zeros(len(Plcl_X))

 for i,x in enumerate(qt):
    if (x>0.1): x = x/1000.
    RH        = mt.mr2pp(x/(1.-x),PPa)/mt.es(TK)
    Plcl_R[i] = lcl(PPa,TK,RH)

 del1 = (Plcl_B-Plcl_X)/100.
 del2 = (Plcl_R-Plcl_X)/100.

 ax2 = plt.subplot(1,2,2)
 ax2.set_xlabel('$q_\mathrm{t}$ / g/kg')
 #ax2.set_ylim(-1.2,1.2)
 ax2.set_yticklabels([])

 plt.plot(qt*1.e3,del1,label='$\\delta_\mathrm{B}$, $T$=310K')
 plt.plot(qt*1.e3,del2,label='$\\delta_\mathrm{R}$, $T$=310K')
 #plt.gca().invert_yaxis()
 plt.legend(loc="best")

 sns.set_context("talk", font_scale=1.2)
 sns.despine(offset=10)
 fig.savefig(plot_dir+'Plcl.pdf')
 ```

 %% Cell type:markdown id: tags:

 ## Acknowledgments ##

 Jiawei Bao, Geet George, and Hauke Schulz are thanked for comments on these notes, and the identification of some errors in earlier versions. Axel Seifert is thanked for his comments and insights, and for pointing out the TEOS routines of Feisel et al. (2010) which extend the IAPWS libraries.

 %% Cell type:markdown id: tags:

 # Notes on calculations of vapor pressures and the specific heats #

 These notes do not use the moist_thermodynamic libararies, rather they were the basis for the choice for the particular formulations for the specific heats and saturation vapor pressures within that library.

 There are a large number of expressions for the saturation and sublimation vapor pressure in the literature, and many of these, even recent ones, seem to reference previous studies in a haphazard way.  So how much do these differ, is there a standard, and by what criteria should one judge them by.  Here I try to develop an intuition for the answers.

 The first thing to note is that there is a community that concerns itself with this question.  They call themselves the international association for the physical properties of water and steam, and mostly concern themselves with the behavior of water at high temperature.  The approach of the IAPWS is to develop an empirical equation of state for water, in the form of a specification of its Helmholtz free energy, or potential, from which all other properties can be derived.  The standard reference for the IAPWS equation of state is the publication by Wagner and Pru{\ss} (Thermodynamic Properties of Ordinary Water) published in 2002 and which describes the IAPWS-95 approved formulation.  Minor corrections have since been made to this, which as best I can tell are relevant at high temperatures.  The most substantial change has been the TEOS-10 work by Rainer Feistel of IOW, which extends these approaches to composite systems, thereby allowing for representations of sea-water and moist air.  By working with an equation of state, all properties of water, from the specific heats to the gas constants to the phase-change enthalpies can be derived consistently.  The disadvantage of this approach is that the equation is derived by positing an analytic form that is then fit to a very wide and diverse abundance of existing data.  The resultant equation is described in an ideal part, which involves a summation of nine terms and thirteen coefficients, and a residual part, with more than 50 terms and over 200 constants.

 For the case of the saturation vapor pressure over water Wagner and Pru{\ss} suggest a much simpler equation that is described in terms of only six coefficients. First, below I compare the relative error to the IAPWS standard as has been formlated and distributed in the iapws python package, version (1.4).  There has been some discussion on the web of its implementation, but the similarity with the Wagner and Pru{\ss} formulation gives me confidence.  Next we look at the TEOS-10 Sea-Ice-Air formulations of Feistel et al (Ocean Sci., 6, 91–141, 2010) which are distributed as FORTRAN90 code which I downloaded, ran, and tabulated to assess some empirical fits later used as a reference.

 %% Cell type:code id: tags:

 ``` python
 import numpy as np
 import matplotlib.pyplot as plt
 import seaborn as sns
 import os
 from scipy import interpolate, optimize

 os.makedirs("plots", exist_ok=True)

 !%matplotlib inline
 ```

 %% Output

    zsh:fg:1: no job control in this shell.

 %% Cell type:code id: tags:

 ``` python
 gravity = 9.8076

 cpd     = 1006.
 Rd      = 287.05

 Rv      = 461.53    # IAPWS97 at 273.15
 cpv     = 1865.01   # ''
 lv0     = 2500.93e3 # IAPWS97 at 273.15
 lf0     =  333.42e3 #''

 cl      = 4179.57   # IAPWS97 at 305 and P=0.1 MPa (chosen to give a good fit for es over ice)
 ci      = 1905.43   # IAPWS97 at 247.065 and P=0.1 MPa (chosen to give a good fit for es over ice)

 eps1     = Rd/Rv
 eps2     = Rv/Rd -1.

 P0      = 100000.  # Standard Pressure
 T0      = 273.15   # Standard Temperature
 PvC     = 22.064e6 # Critical pressure of water vapor
 TvC     = 647.096  # Critical temperature of water vapor
 TvT     = 273.16   # Triple point temperature of water
 PvT     = 611.655
 lvT     = lv0 + (cpv-cl)*(TvT-T0)
 lfT     = lf0 + (cpv-ci)*(TvT-T0)
 lsT     = lvT + lfT

 es_default = 'sonntag'

 def thermo_input(x, xtype='none'):

    import numpy as np

    x = np.asarray(x).flatten()
    scalar_input = False
    if x.ndim == 0:
        x = x[None]  # Makes x 1D
        scalar_input = True

    if (xtype == 'Kelvin' and x.max() < 100 ): x = x+273.15
    if (xtype == 'Celcius'and x.max() > 100 ): x = x-273.15
    if (xtype == 'Pascal' and x.max() < 1200): x = x*100.
    if (xtype == 'kg/kg'  and x.max() > 1.0) : x = x/1000.
    if (xtype == 'meter'  and x.max() < 10.0): print('Warning: input should be in meters, max value less than 10, not corrected')

    return x, scalar_input

 def eslf(T, formula=es_default):
    """ Returns the saturation vapour pressure [Pa] over liquid given
 	the temperature.  Temperatures can be in Celcius or Kelvin.
 	Formulas supported are
 	  - Goff-Gratch (1994 Smithsonian Tables)
 	  - Sonntag (1994)
 	  - Flatau
 	  - Magnus Tetens (MT)
      - Romps (2017)
      - Murphy-Koop
      - Bolton
      - Wagner and Pruss (WP, 2002) is the default
 	>>> eslf(273.16)
 	611.657
    """
    import numpy as np

    x,  scalar_input = thermo_input(T, 'Kelvin')

    if formula == "flatau":
        if (np.min(x) > 100): x = x-273.16
        np.maximum(x,-80.)
        c_es= np.asarray([0.6105851e+03, 0.4440316e+02, 0.1430341e+01, 0.2641412e-01,
                           0.2995057e-03,0.2031998e-05,0.6936113e-08,0.2564861e-11,-0.3704404e-13])
        es = np.polyval(c_es[::-1],x)
    elif formula == "bolton":
        if (np.min(x) > 100): x = x-273.15
        es = 611.2*np.exp((17.67*x)/(243.5+x))
    elif formula == "sonntag":
        xx = -6096.9385/x + 16.635794 - 2.711193e-2*x + 1.673952e-5*x*x + 2.433502 * np.log(x)
        es = 100.*np.exp(xx)
    elif formula =='goff-gratch':
        x1 = 273.16/x
        x2 = 373.16/x
        xl = np.log10(1013.246 ) - 7.90298*(x2 - 1) + 5.02808*np.log10(x2) - 1.3816e-7*(10**(11.344*(1.-1./x2)) - 1.0) + 8.1328e-3 * (10**(-3.49149*(x2-1)) - 1.0)
        es =10**(xl+2) # plus 2 converts from hPa to Pa
    elif formula == 'wagner-pruss':
        vt = 1.-x/TvC
        es = PvC * np.exp(TvC/x * (-7.85951783*vt + 1.84408259*vt**1.5 - 11.7866497*vt**3 + 22.6807411*vt**3.5 - 15.9618719*vt**4 + 1.80122502*vt**7.5))
    elif formula == 'hardy98':
        y  = -2.8365744e+3/(x*x) - 6.028076559e+3/x + 19.54263612 - 2.737830188e-2*x + 1.6261698e-5*x**2 + 7.0229056e-10*x**3 - 1.8680009e-13*x**4 + 2.7150305 * np.log(x)
        es = np.exp(y)
    elif formula == 'romps':
        Rr    = 461.
        cvl_r = 4119
        cvv_r = 1418
        cpv_r = cvv_r + Rr
        es = 611.65 * (x/TvT) **((cpv_r-cvl_r)/Rr) * np.exp((2.37403e6 - (cvv_r-cvl_r)*TvT)*(1/TvT - 1/x)/Rr)
    elif formula == "murphy-koop":
        es = np.exp(54.842763 - 6763.22/x - 4.210*np.log(x) + 0.000367*x + np.tanh(0.0415*(x - 218.8)) * (53.878 - 1331.22/x - 9.44523 * np.log(x) + 0.014025*x))
    elif formula == "standard-analytic":
        c1 = (cpv-cl)/Rv
        c2 = lvT/(Rv*TvT) - c1
        es = PvT * np.exp(c2*(1.-TvT/x)) * (x/TvT)**c1
    else:
        exit("formula not supported")

    es = np.maximum(es,0)
    if scalar_input:
        return np.squeeze(es)
    return es

 def esif(T, formula=es_default):
    """ Returns the saturation vapour pressure [Pa] over ice given
 	the temperature.  Temperatures can be in Celcius or Kelvin.
 	uses the Goff-Gratch (1994 Smithsonian Tables) formula
 	>>> esli(273.15)
 	6.112
 m    """
    import numpy as np

    x,  scalar_input = thermo_input(T, 'Kelvin')

    if formula == "sonntag":
        es = 100 * np.exp(24.7219 - 6024.5282/x + 0.010613868*x - 0.000013198825*x**2 - 0.49382577*np.log(x))
    elif formula == "goff-gratch":
        x1 = 273.16/x
        xi = np.log10(   6.1071) - 9.09718*(x1 - 1) - 3.56654*np.log10(x1) + 0.876793*(1 - 1./x1)
        es = 10**(xi+2)
    elif formula == "wagner-pruss": #(actually wagner et al, 2011)
        a1 = -0.212144006e+2
        a2 =  0.273203819e+2
        a3 = -0.610598130e+1
        b1 =  0.333333333e-2
        b2 =  0.120666667e+1
        b3 =  0.170333333e+1
        theta = T/TvT
        es = PvT * np.exp((a1*theta**b1 + a2 * theta**b2 + a3 * theta**b3)/theta)
    elif formula == "murphy-koop":
        es = np.exp(9.550426 - 5723.265/x + 3.53068 * np.log(x) - 0.00728332*x)
    elif formula == "romps":
        Rr    = 461.
        cvv_r = 1418.
        cvs_r = 1861.
        cpv_r = cvv_r + Rr
        es = 611.65 * (x/TvT) **((cpv_r-cvs_r)/Rr) * np.exp((2.37403e6 + 0.33373e6 - (cvv_r-cvs_r)*TvT)*(1/TvT - 1/x)/Rr)
    elif formula == "standard-analytic":
        c1 = (cpv-ci)/Rv
        c2 = lsT/(Rv*TvT) - c1
        es = PvT * np.exp(c2*(1.-TvT/x)) * (x/TvT)**c1
    else:
        exit("formula not supported")

    es = np.maximum(es,0)
    if scalar_input:
        return np.squeeze(es)
    return es

 def esilf(T,formula=es_default):
    import numpy as np
    return np.minimum(esif(T,formula),eslf(T,formula))

 def es(T,formula=es_default,state='liq'):

    import numpy as np

    x,  scalar_input = thermo_input(T, 'Kelvin')

    if (state == 'liq'):
        return eslf(x,formula)
    if (state == 'ice'):
        return esif(x,formula)
    if (state == 'mxd'):
        return esilf(x,formula)

 def des(T,formula=es_default,state='liq'):

    import numpy as np

    x,  scalar_input = thermo_input(T, 'Kelvin')
    dx = 0.01; xp = x+dx/2; xm = x-dx/2
    return (es(xp,formula,state)-es(xm,formula,state))/dx

 def dlnesdlnT(T,formula=es_default,state='liq'):

    import numpy as np

    x,  scalar_input = thermo_input(T, 'Kelvin')
    dx = 0.01; xp = x+dx/2; xm = x-dx/2
    return ((es(xp,formula,state)-es(xm,formula,state))/es(x,formula,state) * (x/dx))

 def phase_change_enthalpy(Tx,fusion=False):
    """ Returns the enthlapy [J/g] of vaporization (default) of water vapor or
    (if fusion=True) the fusion anthalpy.  Input temperature can be in degC or Kelvin
    >>> phase_change_enthalpy(273.15)
    2500.8e3
    """
    import numpy as np

    TC, scalar_input = thermo_input(Tx, 'Celcius')
    TK, scalar_input = thermo_input(Tx, 'Kelvin')

    if (fusion):
        el = lfT + (cl-ci)*(TK-TvT)
    else:
        el = lvT + (cpv-cl)*(TK-TvT)

    if scalar_input:
        return np.squeeze(el)
    return el

 ```

 %% Cell type:markdown id: tags:

 ## Selected properties of IAWPS water ##

 These routines only provide information on $c_{p,\mathrm{liq}}$ to a temperature of 253 K or -20$^\circ$C, but already demonstrate its divergent behavior (exponential increase) with increased super-cooling.  They also demonstrate the near linearity of $c_{p,\mathrm{ice}}.$

 %% Cell type:code id: tags:

 ``` python
 import iapws

 print ('Using IAPWS Version %s\n'%(iapws.__version__,))
 T = np.arange(183.15,313.15)
 ci_iapws = np.full(len(T),np.nan)
 cl_iapws = np.full(len(T),np.nan)
 for i,Tx in enumerate(T):
    if (Tx < 283): ci_iapws[i] =  iapws._iapws._Ice(Tx, 0.1)['cp']*1000 / ci
    if (Tx > 253.15): cl_iapws[i] =  iapws._iapws._Liquid(Tx, 0.1)['cp']*1000 / cl

 fig = plt.figure(figsize=(4,3))

 ax1 = plt.subplot(1,1,1)
 ax1.set_xlabel('$T$ / K')
 ax1.set_ylabel('$c_\mathrm{i}$ / %5.2f, $c_\mathrm{l}$ / %5.2f'%(ci,cl))
 ax1.set_xticks([185,247.07,273.15,305.00])
 plt.scatter([247.065],[1.])
 plt.scatter([305.000],[1.])
 plt.plot(T,ci_iapws)
 plt.plot(T,cl_iapws)

 sns.set_context("paper", font_scale=1.2)
 sns.despine(offset=10)

 plt.tight_layout()
 fig.savefig(plot_dir+'cp-Tdependance.pdf')

 TK = np.arange(273.15,315.15,0.01)
 es_iapws = np.zeros(len(TK))
 for i, x in enumerate(TK):
    es_iapws[i] = iapws.iapws97._PSat_T(x) *1.e6 #Temperature, [K]; Returns:Pressure, [MPa]
 ```

 %% Output

    Using IAPWS Version 1.5.2
    

    /Users/m219063/opt/miniforge3/lib/python3.9/site-packages/iapws/_iapws.py:124: UserWarning: Metastable ice
      warnings.warn("Metastable ice")

    ---------------------------------------------------------------------------
    NameError                                 Traceback (most recent call last)
 Input     In [3], in <cell line: 26>()
         23 sns.despine(offset=10)
         25 plt.tight_layout()
    ---> 26 fig.savefig(plot_dir+'cp-Tdependance.pdf')
         28 TK = np.arange(273.15,315.15,0.01)
         29 es_iapws = np.zeros(len(TK))
    NameError: name 'plot_dir' is not defined



 %% Cell type:markdown id: tags:

 ## Comparison with Sea-Air-Ice library of Feistel et al (2010) v4.0.1

 Here we compare different thermodynamic constants or empirical formula to the IAPWS-10 and TEOS-10 standards taken from the Sea-Air-Ice library.  These are calculated based on fits to potential functions as described above.  The libraries are run off-line and the output is tabulated for comparison. The Feistel et al formulation extends to the IAPWS formulation shown above to allow for representations of liquid to the homogeneous freezing point of ice.

 %% Cell type:code id: tags:

 ``` python
 import pandas as pd
 sia_dir = './data/'
 d0    = pd.read_csv(sia_dir + 'psat.dat',sep=' ',names=['T','Psat_liq','Psat_ice'])
 x0    = d0.to_xarray()
 d1    = pd.read_csv(sia_dir + 'cp-liq-vap.dat',sep=' ',names=['T','liq_density','vap_density','cp_liq','cp_vap','lv'])
 x1    = d1.to_xarray()
 d2    = pd.read_csv(sia_dir + 'cp-ice-vap.dat',sep=' ',names=['T','liq_density','vap_density','cp_ice','cp_vap','ls'])
 x2    = d2.to_xarray()

 fig = plt.figure(figsize=(7,14/1.610834))

 ax1 = plt.subplot(2,1,2)
 ax1.set_ylabel('difference / %')
 ax1.set_xlabel('T / K')

 formula = 'wagner-pruss'
 es_r = es(x0.T,formula=formula,state='ice')
 diff = es_r/x0.Psat_ice
 plt.plot(x0.T,100*(diff.where(diff>1.e-4)-1),c='dodgerblue',label=formula)
 es_r = es(x0.T,formula=formula,state='liq')
 diff = es_r/x0.Psat_liq
 plt.plot(x0.T,100*(diff.where(diff>1.e-4)-1),c='orange',label=formula)

 formula = 'romps'
 es_r = es(x0.T,formula=formula,state='ice')
 diff = es_r/x0.Psat_ice
 plt.plot(x0.T,100*(diff.where(diff>1.e-4)-1),c='dodgerblue',ls='dashed',label=formula)
 es_r = es(x0.T,formula=formula,state='liq')
 diff = es_r/x0.Psat_liq
 plt.plot(x0.T,100*(diff.where(diff>1.e-4)-1),c='orange',ls='dashed',label=formula)

 formula = 'murphy-koop'
 es_r = es(x0.T,formula=formula,state='ice')
 diff = es_r/x0.Psat_ice
 plt.plot(x0.T,100*(diff.where(diff>1.e-4)-1),c='dodgerblue',ls='dotted',label=formula)
 es_r = es(x0.T,formula=formula,state='liq')
 diff = es_r/x0.Psat_liq
 plt.plot(x0.T,100*(diff.where(diff>1.e-4)-1),c='orange',ls='dotted',label=formula)
 plt.legend(loc="best",ncol=3)

 ###
 ax2 = plt.subplot(2,1,1)
 ax2.set_ylabel('difference / %')
 ax2.set_xlabel('T / K')
 ax2.set_ylim(0.8,1.2)

 plt.plot(x1.T,x1.cp_vap/cpv,c='green',ls='solid',label='$c_{p,\mathrm{vap}}$ (vap-liq)')
 plt.plot(x2.T,x2.cp_vap/cpv,c='green',ls='dashed',label='$c_{p,\mathrm{vap}}$ (vap-ice)')

 plt.plot(x1.T,x1.cp_liq/cl,c='orange',ls='solid',label='$c_{p,\mathrm{liq}}$')
 plt.plot(x2.T,x2.cp_ice/ci,c='dodgerblue',ls='solid',label='$c_{p,\mathrm{ice}}$')

 lvx = phase_change_enthalpy(x1.T,fusion=False)
 lsx = phase_change_enthalpy(x2.T,fusion=True) + phase_change_enthalpy(x2.T,fusion=False)

 plt.plot(x1.T,x1.lv/lvx,c='gray',ls='solid',label='$\\ell_v$')
 plt.plot(x2.T,x2.ls/lsx,c='gray',ls='dashed',label='$\\ell_s$')

 plt.legend(loc="lower right",ncol=3)
 sns.set_context("paper")
 sns.despine(offset=10)
 plt.tight_layout()
 ```

+%% Output
+
+
+
 %% Cell type:markdown id: tags:

 ## Comparing formulations of saturation vapor pressure ##

 This comparison of relative error suggests that the Wagner-Pru{\ss}, Murphy and Koop, Hardy, and Sonntag formulations lie closest to the IAPWS-97 reference.  Romps (2017) and Bolton (1980) are similarly accurate and may have advantages.  Hardy is interesting as it appears in a technical document and is rarely mentioned in the subsequent literature, but used by Vaisala in the calibration of their sondes

 %% Cell type:markdown id: tags:

 ### Temperatures above the triple point ###

 %% Cell type:code id: tags:

 ``` python
 state = 'liq'
 fig = plt.figure(figsize=(10,5))
 ax1 = plt.subplot(1,1,1)
 ax1.set_xlabel('$T$ / K')
 ax1.set_ylabel('$e_{\mathrm{s,x}}/e_{\mathrm{s,ref}} - 1$')
 ax1.set_yscale('log')

 es_ref = es_iapws
 es_w = es(TK,formula="wagner-pruss",state=state)
 es_r = es(TK,formula='romps',state=state)
 es_g = es(TK,formula='goff-gratch',state=state)
 es_m = es(TK,formula='murphy-koop',state=state)
 es_s = es(TK,formula='sonntag',state=state)
 es_b = es(TK,formula='bolton',state=state)
 es_f = es(TK,formula='flatau',state=state)
 es_h = es(TK,formula='hardy98',state=state)
 es_a = es(TK,formula='standard-analytic',state=state)

 plt.plot(TK,np.abs(es_h/es_ref-1),c='tab:blue',ls='solid',label='Hardy (1998)')
 plt.plot(TK,np.abs(es_f/es_ref-1),c='tab:orange',label='Flatau (1992)')
 plt.plot(TK,np.abs(es_g/es_ref-1),c='tab:green',label='Goff-Gratch (1957)')
 plt.plot(TK,np.abs(es_b/es_ref-1),c='tab:red',ls='dotted',label='Bolton (1980)')

 plt.plot(TK,np.abs(es_r/es_ref-1),c='tab:purple',label='Romps (2017)')
 plt.plot(TK,np.abs(es_s/es_ref-1),c='tab:grey',label='Sonntag (1990)')
 plt.plot(TK,np.abs(es_m/es_ref-1),c='tab:pink',label='Murphy-Koop (2005)')
 plt.plot(TK,np.abs(es_w/es_ref-1),c='tab:brown',label='Wagner-Pruss (2002)')
 plt.plot(TK,np.abs(es_a/es_ref-1),c='tab:purple',ls='dotted',label='Analytic')

 plt.legend(loc="lower right",ncol=3)

 sns.set_context("paper", font_scale=1.2)
 sns.despine(offset=10)

 fig.savefig(plot_dir+'es_l-error.pdf')
 ```

+%% Output
+
+    ---------------------------------------------------------------------------
+    NameError                                 Traceback (most recent call last)
+Input     In [5], in <cell line: 8>()
+          5 ax1.set_ylabel('$e_{\mathrm{s,x}}/e_{\mathrm{s,ref}} - 1$')
+          6 ax1.set_yscale('log')
+    ----> 8 es_ref = es_iapws
+          9 es_w = es(TK,formula="wagner-pruss",state=state)
+         10 es_r = es(TK,formula='romps',state=state)
+    NameError: name 'es_iapws' is not defined
+
+
+
 %% Cell type:markdown id: tags:

 ### Extension to the freezing regime ###

 To extend over the entire temperature range a different reference is required, for this any of the Hardy, Sonntag, Murphy-Koop and Wagner-Pru{\ss} formulations could suffice.  We choose Wagner-Pru{\ss} because Wagner's group is responsible for the standard, and has also developed the IAPWS standard for saturation vapor pressure over ice.  Below the results are plooted with respect to this standard over a much larger temperature range.

 It is not clear how accurate Wagner and Pru{\ss} wis hen extended well beyond the IAPWS range, based on which it might be that the grouping of errors of similar magnitude from the Bolton, Flatau and Goff-Gratch formulations are indicative of a low temperature bias in the Wagner-Pru{\ss} formualtion.  I doubt that this is the case, as the poor performance of all these formulations in the higher temperature range, and the simplicity of their formulation make it unlikely.   The agreement of the Murphy-Koop formulation with these simpler formulations at low temperature may be indicative of Murphy and Koops focus on saturation over ice rather than liquid.

 %% Cell type:code id: tags:

 ``` python
 state = 'liq'
 fig = plt.figure(figsize=(10,5))
 ax1 = plt.subplot(1,1,1)
 ax1.set_xlabel('$T$ / K')
 ax1.set_ylabel('$e_{\mathrm{s,x}}/e_{\mathrm{s,ref}} - 1$')
 ax1.set_yscale('log')

 TK = np.arange(180,320,0.5)

 es_w = es(TK,formula="wagner-pruss",state=state)
 es_r = es(TK,formula='romps',state=state)
 es_g = es(TK,formula='goff-gratch',state=state)
 es_m = es(TK,formula='murphy-koop',state=state)
 es_s = es(TK,formula='sonntag',state=state)
 es_b = es(TK,formula='bolton',state=state)
 es_f = es(TK,formula='flatau',state=state)
 es_h = es(TK,formula='hardy98',state=state)
 es_a = es(TK,formula='standard-analytic',state=state)

 es_ref = es_w

 plt.plot(TK,np.abs(es_h/es_ref-1),c='tab:blue',ls='solid',label='Hardy (1998)')
 plt.plot(TK,np.abs(es_f/es_ref-1),c='tab:orange',label='Flatau (1992)')
 plt.plot(TK,np.abs(es_g/es_ref-1),c='tab:green',label='Goff-Gratch (1957)')
 plt.plot(TK,np.abs(es_b/es_ref-1),c='tab:red',ls='dotted',label='Bolton (1980)')

 plt.plot(TK,np.abs(es_r/es_ref-1),c='tab:purple',label='Romps (2017)')
 plt.plot(TK,np.abs(es_s/es_ref-1),c='tab:grey',label='Sonntag (1990)')
 plt.plot(TK,np.abs(es_m/es_ref-1),c='tab:pink',label='Murphy-Koop (2005)')
 plt.plot(TK,np.abs(es_a/es_ref-1),c='tab:purple',ls='dotted',label='Analytic')

 #plt.plot(TK,np.abs(es_w/es_ref-1),c='tab:olive',label='Wagner-Pruss (2002)')

 plt.legend(loc="lower left",ncol=2)

 sns.set_context("paper", font_scale=1.2)
 sns.despine(offset=10)

 fig.savefig(plot_dir+'es_lsc-error.pdf')
 ```

 %% Cell type:markdown id: tags:

 ## Sublimation vapor pressure ##

 A subset of the formulations also postulate forms for the saturation vapor pressure over ice.  For the reference in this quantity we use Wagner et al., (2011) as this has been adopted as the IAPWS standard.   Here is seems that Murphy and Koop's (2005) formulation behaves very well in comparision to Wagner et al., but Sonntag is also quite adequate, particularly at lower ($T<273.15$ K) temperatures where it is likely to be applied.

 %% Cell type:code id: tags:

 ``` python
 state = 'ice'
 fig = plt.figure(figsize=(10,5))
 ax1 = plt.subplot(1,1,1)
 ax1.set_xlabel('$T$ / K')
 ax1.set_ylabel('$e_{\mathrm{s,x}}/e_{\mathrm{s,ref}} - 1$')
 ax1.set_yscale('log')

 TK = np.arange(180,320,0.5)

 es_w = es(TK,formula="wagner-pruss",state=state)
 es_r = es(TK,formula='romps',state=state)
 es_g = es(TK,formula='goff-gratch',state=state)
 es_m = es(TK,formula='murphy-koop',state=state)
 es_s = es(TK,formula='sonntag',state=state)
 es_a = es(TK,formula='standard-analytic',state=state)
 es_ref = es_w

 plt.plot(TK,np.abs(es_g/es_ref-1),c='tab:green',label='Goff-Gratch (1957)')
 plt.plot(TK,np.abs(es_r/es_ref-1),c='tab:purple',label='Romps (2017)')
 plt.plot(TK,np.abs(es_s/es_ref-1),c='tab:grey',label='Sonntag (1990)')
 plt.plot(TK,np.abs(es_m/es_ref-1),c='tab:pink',label='Murphy-Koop (2005)')
 plt.plot(TK,np.abs(es_a/es_ref-1),c='tab:purple',ls='dotted',label='Analytic')

 #plt.plot(TK,np.abs(es_w/es_ref-1),c='tab:olive',label='Wagner-Pruss (2002)')

 plt.legend(loc="lower left",ncol=2)

 sns.set_context("paper", font_scale=1.2)
 sns.despine(offset=10)

 fig.savefig(plot_dir+'es_i-error.pdf')
 ```

 %% Cell type:markdown id: tags:

 ## Clausius Clapeyron ##

 Often over looked is that many conceptual models are built on the application of the Clausius-Clapeyron equation,
 \begin{equation}
 \frac{\mathrm{d} \ln e_\mathrm{s}}{\mathrm{d \ln T}} \left(\frac{\ell_\mathrm{v}}{R_\mathrm{v} T}\right)^{-1} = 1
 \end{equation}
 with the assumption that the vaporization enthalpy, $\ell_\mathrm{v}$ that appears in this equation, is linear in temperature following Kirchoff's relation.  This is similar to assuming that the specific heats are independent of temeprature, an idealization which is, unfortunately, just that, and idealization.

 But because of this it is interesting to compare this expression as given by the above formulation of the saturation vapor pressure (through their numerical derivative) and independent expressions of $\ell_\mathrm{v}$ based on the assumption of constant specific heats.

 This is shown below for ice and liquid saturation.  The analytic expression, which has larger errors for es is constructued to satisfy this relationship and is exact to the precision of the numerical calculations.  The various formulations using more accurate expressions for $e_s$ which implicityl don't assume constancy in specific heats are similarly accurate, with the exception of Goff-Gratch, and Romps for Ice.  Hardy is only shown for water.  For ice Sonntag does not behave well for $T> 290$ K, but it is not likely to be used at these temperatures.  Note that Romps would be perfect had we adopted his modified specific heats.

 Based on the above my recommendation is to use the formulations by Wagner's group, unless one is interested in very low temperatures ($T<180$K) in which case the formulation of Koop and Murphy may be desirable.  For just liquid processes Hardy might be a good choice, it is less well known but used by Vaisala for its sondes.   There may be advantages to using Sonntag if there is interest in liquid and ice as it might allow more efficient implementations, but for my tests all formulations were within 30% of one another.

 Another alternative, would be to use the analytic approach, either using Romps' formulae if getthing the staturation vapor pressure as close to measurements as possible is preferred, or using the analytic formula with the correct (at the standard temperature and pressure) specific heats and gast constants.

 %% Cell type:code id: tags:

 ``` python
 state = 'liq'

 fig = plt.figure(figsize=(10,10))

 ax1 = plt.subplot(2,1,1)

 ax1.set_ylabel('$|\mathrm{CC}_\mathrm{liq} - 1|$')
 ax1.set_yscale('log')
 ax1.set_xticklabels([])

 TK = np.arange(180,320,0.5)

 lv   = phase_change_enthalpy(TK)
 if (state == 'ice'): lv += phase_change_enthalpy(TK,fusion=True)

 y    = lv/(Rv * TK)
 cc_w = dlnesdlnT(TK,formula="wagner-pruss",state=state) / y
 cc_r = dlnesdlnT(TK,formula='romps',state=state) /y
 cc_g = dlnesdlnT(TK,formula='goff-gratch',state=state) /y
 cc_m = dlnesdlnT(TK,formula='murphy-koop',state=state) /y
 cc_s = dlnesdlnT(TK,formula='sonntag',state=state) /y
 cc_h = dlnesdlnT(TK,formula='hardy98',state=state) /y
 cc_a = dlnesdlnT(TK,formula='standard-analytic',state=state) /y

 plt.plot(TK,np.abs(cc_h/1 -1.),c='tab:blue',label='Hardy (1998)')
 plt.plot(TK,np.abs(cc_g/1 -1.),c='tab:green',label='Goff-Gratch (1957)')
 plt.plot(TK,np.abs(cc_r/1 -1.),c='tab:purple',label='Romps (2017)')
 plt.plot(TK,np.abs(cc_s/1 -1.),c='tab:grey',label='Sonntag (1990)')
 plt.plot(TK,np.abs(cc_m/1 -1.),c='tab:pink',label='Murphy-Koop (2005)')
 plt.plot(TK,np.abs(cc_w/1 -1.),c='tab:olive',label='Wagner-Pruss (2002)')
 plt.plot(TK,np.abs(cc_a/1 -1.),c='tab:purple',ls='dotted',label='Analytic')

 plt.legend(loc="lower left",ncol=1)

 state = 'ice'
 TK = np.arange(180,320,0.5)

 lv   = phase_change_enthalpy(TK)
 if (state == 'ice'): lv   = phase_change_enthalpy(TK,fusion=True) + phase_change_enthalpy(TK)

 y    = lv/(Rv * TK)
 cc_w = dlnesdlnT(TK,formula="wagner-pruss",state=state) / y
 cc_r = dlnesdlnT(TK,formula='romps',state=state) /y
 cc_g = dlnesdlnT(TK,formula='goff-gratch',state=state) /y
 cc_m = dlnesdlnT(TK,formula='murphy-koop',state=state) /y
 cc_s = dlnesdlnT(TK,formula='sonntag',state=state) /y
 cc_a = dlnesdlnT(TK,formula='standard-analytic',state=state) /y

 ax2 = plt.subplot(2,1,2)
 ax2.set_xlabel('$T$ / K')
 ax2.set_ylabel('$|\mathrm{CC}_\mathrm{ice} - 1|$')
 ax2.set_yscale('log')

 plt.plot(TK,np.abs(cc_g/1 -1.),c='tab:green',label='Goff-Gratch (1957)')
 plt.plot(TK,np.abs(cc_r/1 -1.),c='tab:purple',label='Romps (2017)')
 plt.plot(TK,np.abs(cc_s/1 -1.),c='tab:grey',label='Sonntag (1990)')
 plt.plot(TK,np.abs(cc_m/1 -1.),c='tab:pink',label='Murphy-Koop (2005)')
 plt.plot(TK,np.abs(cc_w/1 -1.),c='tab:olive',label='Wagner-Pruss (2002)')
 plt.plot(TK,np.abs(cc_a/1. -1.),c='tab:purple',ls='dotted',label='Analytic')

 sns.set_context("paper", font_scale=1.2)
 sns.despine(offset=10)

 fig.savefig(plot_dir+'cc-error.pdf')
 ```

 %% Cell type:markdown id: tags:

 ## Optimizing analytic fits for saturation and sublimation vapor pressure ##

 Romps suggests modifying the specific heats of liquid, ice and the gas constant of vapor to arrive at an optimal fit for the saturation vapor pressure using the analytic form.  One can do almost as good by just modifying the specific heat of the condensate phases.  Here we show how the maximum error in the fit depends on the specific heat of the condensate phases as compared to the reference, and how we arrive at our optimal fit by only manipulating the condensate phase specific heats to values that they anyway adopt within the range of temperatures spanned by the atmosphere.  This justifys the default choice for saturation vapor pressure and the specific heats used in aes_thermo.py

 %% Cell type:code id: tags:

 ``` python
 fig = plt.figure(figsize=(10,5))

 cl_1 = (iapws._iapws._Liquid(265, 0.1)['cp'])*1000.
 cl_2 = (iapws._iapws._Liquid(305, 0.1)['cp'])*1000
 ci_1 = (iapws._iapws._Ice(193, 0.01)['cp'])*1000.
 ci_2 = (iapws._iapws._Ice(273, 0.10)['cp'])*1000

 cls = np.arange(cl_2,cl_1)
 err = np.zeros(len(cls))

 ax1 = plt.subplot(1,2,1)
 ax1.set_xlabel('$c_\mathrm{liq}$ / Jkg$^{-1}$K$^{-1}$')
 ax1.set_ylabel('$(e_{\mathrm{s,x}}/e_{\mathrm{s,ref}} - 1)_\mathrm{max}$ / %')
 ax1.set_yscale('log')

 state = 'liq'
 TK = np.arange(260,300,0.5)
 es_ref = es(TK,formula="wagner-pruss",state=state)
 for i,cx in enumerate(cls):
    c1 = (cpv-cx)/Rv
    c2 = lvT/(Rv*TvT) - c1
    es_a = PvT * np.exp(c2*(1.-TvT/TK)) * (TK/TvT)**c1
    err[i] = np.max(np.abs(es_a/es_ref -1.))*100.
 ax1.plot(cls,err,c='tab:purple',ls='dotted',label='Analytic $c_\mathrm{liq}$ for $T\in$ (260K,305K)')
 ax1.legend(loc="upper left",ncol=2)

 cis = np.arange(ci_1,ci_2)
 err = np.zeros(len(cis))

 ax2 = plt.subplot(1,2,2)
 ax2.set_xlabel('$c_\mathrm{ice}$ / Jkg$^{-1}$K$^{-1}$')
 ax2.set_ylabel('$(e_{\mathrm{s,x}}/e_{\mathrm{s,ref}} - 1)_\mathrm{max}$ / %')
 ax2.set_yscale('log')

 state = 'ice'
 TK = np.arange(180,273,0.5)
 es_ref = es(TK,formula="wagner-pruss",state=state)
 for i,cx in enumerate(cis):
    c1 = (cpv-cx)/Rv
    c2 = lsT/(Rv*TvT) - c1
    es_a = PvT * np.exp(c2*(1.-TvT/TK)) * (TK/TvT)**c1
    err[i] = np.max(np.abs(es_a/es_ref -1.))*100.
 ax2.plot(cis,err,c='tab:purple',ls='dotted',label='Analytic $c_\mathrm{ice}$ for $T\in$ (193K,273K)')
 ax2.legend(loc="upper right",ncol=2)

 sns.set_context("paper", font_scale=1.2)
 sns.despine(offset=10)

 fig.savefig(plot_dir+'es-analytic-fits.pdf')
 Tfit = 305
 print ('Taking fit for $c_\mathrm{liq}=$ %3.2f J/(kg K) at $T=$ %3.2f K'%(iapws._iapws._Liquid(Tfit, 0.1)['cp']*1000.,Tfit))
 Tfit = 247.065
 print ('Taking fit for $c_\mathrm{ice}=$ %3.2f J/(kg K) at $T=$ %3.2f K'%(iapws._iapws._Ice(Tfit, 0.1)['cp']*1000.,Tfit))
 ```

 %% Cell type:markdown id: tags:

 ## RCEMIP comparision ##

 During RCEMIP (Wing et al.) different models output different RH, differing in ways of calculating it and also whether or not it was calculated relative to liquid or ice.  In this analysis we create a python implementation of the intial RCEMIP sounding and then for the given state estimate the RH using different formulat and different assumptions regarding the reference condensate (liquid/ice).  We also show the difference associated with 1 K of temperature.

 %% Cell type:code id: tags:

 ``` python
 def rcemip_on_z(z,SST):
    # function [T,q,p] = rcemip_on_z(z,SST)
    #
    # Inputs:
    # z: array of heights (low to high, m)
    # SST: sea surface temperature (K)
    #
    # Outputs:
    T = np.zeros(len(z)) # temperature (K)
    q = np.zeros(len(z)) # specific humidity (g/g)
    p = np.zeros(len(z)) # pressure (Pa)

    ## Constants
    g = 9.79764 #m/s^2
    Rd = 287.04 #J/kgK

    ## Parameters
    p0 = 101480   #Pa surface pressure
    qt = 10**(-11) #g/g specific humidity at tropopause
    zq1 = 4000 #m
    zq2 = 7500 #m
    zt = 15000 #m tropopause height
    gamma = 0.0067 #K/m lapse rate

    ## Scratch
    Tv = np.zeros(len(z)) # temperature (K)

    if SST == 295:
        q0 = 0.01200; #g/g specific humidity at surface (adjusted from 300K value so RH near surface approx 80%)
    elif SST == 300:
        q0 = 0.01865; #g/g specific humidity at surface
    elif SST == 305:
        q0 = 0.02400 #g/g specific humidity at surface (adjusted from 300K value so RH near surface approx 80%)

    T0 = SST - 0 #surface air temperature adjusted to be 0K less than SST

    ## Virtual Temperature at surface and tropopause
    Tv0 = T0*(1 + 0.608*q0) #virtual temperature at surface
    Tvt = Tv0 - gamma*zt #virtual temperature at tropopause z=zt

    ## Pressure
    pt = p0*(Tvt/Tv0)**(g/(Rd*gamma)); #pressure at tropopause z=zt
    p  = p0*((Tv0-gamma*z)/Tv0)**(g/(Rd*gamma)) #0 <= z <= zt
    p[z>zt] = pt*np.exp(-g*(z[z>zt]-zt)/(Rd*Tvt)) #z > zt

    ## Specific humidity
    q = q0*np.exp(-z/zq1)*np.exp(-(z/zq2)**2)
    q[z>zt] = qt #z > zt

    ## Temperature
    #Virtual Temperature
    Tv = Tv0 - gamma*z #0 <= z <= zt
    Tv[z>zt] = Tvt #z > zt

    #Absolute Temperature at all heights
    T = Tv/(1 + 0.608*q)

    return T, q, p

 z = np.arange(0,17000,100)
 T, q , p = rcemip_on_z(z,300)
 ```

 %% Cell type:code id: tags:

 ``` python
 def get_rh (T,q,p,formula='wagner-pruss',state='liq'):
    es_w = es(T,formula=formula,state=state)
    x = es_w * eps1/(p-es_w)
    return 100.*q*(1+x)/x

 fig = plt.figure(figsize=(4,5))

 ax1 = plt.subplot(1,1,1)
 ax1.set_ylabel('$z$ / km')
 ax1.set_xlabel('RH / %')
 ax1.set_ylim(0,14.5)
 ax1.set_yticks([0,4,8,12])

 plt.plot(get_rh(T,q,p,state='mxd'),z/1000.,label = 'Wagner Pruss (ice/liq)')
 plt.plot(get_rh(T+1,q,p,state='mxd'),z/1000.,label = 'Wagner Pruss (ice/liq) + 1 K')
 plt.plot(get_rh(T,q,p,state='ice'),z/1000.,label = 'Wagner Pruss (ice)')
 plt.plot(get_rh(T,q,p,formula='romps',state='mxd'),z/1000.,label = 'Romps (ice/liq)')
 plt.plot(get_rh(T,q,p),z/1000.,label = 'Wagner Pruss (liq)')
 plt.plot(get_rh(T,q,p,formula='flatau'),z/1000.,label = 'Flatau (liq)')

 plt.legend(loc="lower left",ncol=1)

 sns.set_context("paper")
 sns.despine(offset=10)
 plt.tight_layout()

 fig.savefig(plot_dir+'RCEMIP-RHerror.pdf')
 ```

 %% Cell type:markdown id: tags:

 ## lifting-condensation level approximations

 Here we compare the LCL base predictions to those proposed by Romps and Bolton as well as the differences between density potential temperatures.

 For the estimation of the LCL we modify the Romps expressions (using his code) to output pressure at the LCL, as this eliminates an assumption as to how pressure is distributed in the atmosphere, and thus only depends on the parcel state.  What we find is that the much simpler Bolton expression is as good as the more complex expression by Romps, and differences between the two are commensurate with those arising from slight differences in how the saturation vapor pressure is calculated.

 %% Cell type:code id: tags:

 ``` python
 # Version 1.0 released by David Romps on September 12, 2017.
 #
 # When using this code, please cite:
 #
 # @article{16lcl,
 #   Title   = {Exact expression for the lifting condensation level},
 #   Author  = {David M. Romps},
 #   Journal = {Journal of the Atmospheric Sciences},
 #   Year    = {2017},
 #   Volume  = {in press},
 # }
 #
 # This lcl function returns the height of the lifting condensation level
 # (LCL) in meters.  The inputs are:
 # - p in Pascals
 # - T in Kelvins
 # - Exactly one of rh, rhl, and rhs (dimensionless, from 0 to 1):
 #    * The value of rh is interpreted to be the relative humidity with
 #      respect to liquid water if T >= 273.15 K and with respect to ice if
 #      T < 273.15 K.
 #    * The value of rhl is interpreted to be the relative humidity with
 #      respect to liquid water
 #    * The value of rhs is interpreted to be the relative humidity with
 #      respect to ice
 # - ldl is an optional logical flag.  If true, the lifting deposition
 #   level (LDL) is returned instead of the LCL.
 # - min_lcl_ldl is an optional logical flag.  If true, the minimum of the
 #   LCL and LDL is returned.

 def lcl(p,T,rh=None,rhl=None,rhs=None,return_ldl=False,return_min_lcl_ldl=False):

    import math
    import scipy.special

    # Parameters
    Ttrip = 273.16     # K
    ptrip = 611.65     # Pa
    E0v   = 2.3740e6   # J/kg
    E0s   = 0.3337e6   # J/kg
    ggr   = 9.81       # m/s^2
    rgasa = 287.04     # J/kg/K
    rgasv = 461        # J/kg/K
    cva   = 719        # J/kg/K
    cvv   = 1418       # J/kg/K
    cvl   = 4119       # J/kg/K
    cvs   = 1861       # J/kg/K
    cpa   = cva + rgasa
    cpv   = cvv + rgasv

    # The saturation vapor pressure over liquid water
    def pvstarl(T):
        return ptrip * (T/Ttrip)**((cpv-cvl)/rgasv) * math.exp( (E0v - (cvv-cvl)*Ttrip) / rgasv * (1/Ttrip - 1/T) )
    # The saturation vapor pressure over solid ice
    def pvstars(T):
        return ptrip * (T/Ttrip)**((cpv-cvs)/rgasv) *  math.exp( (E0v + E0s - (cvv-cvs)*Ttrip) / rgasv * (1/Ttrip - 1/T))

    # Calculate pv from rh, rhl, or rhs
    rh_counter = 0
    if rh  is not None:
        rh_counter = rh_counter + 1
    if rhl is not None:
        rh_counter = rh_counter + 1
    if rhs is not None:
        rh_counter = rh_counter + 1
    if rh_counter != 1:
        print(rh_counter)
        exit('Error in lcl: Exactly one of rh, rhl, and rhs must be specified')
    if rh is not None:
    # The variable rh is assumed to be
    # with respect to liquid if T > Ttrip and
    # with respect to solid if T < Ttrip
        if T > Ttrip:
            pv = rh * pvstarl(T)
        else:
            pv = rh * pvstars(T)
            rhl = pv / pvstarl(T)
            rhs = pv / pvstars(T)
    elif rhl is not None:
        pv = rhl * pvstarl(T)
        rhs = pv / pvstars(T)
        if T > Ttrip:
            rh = rhl
        else:
            rh = rhs
    elif rhs is not None:
        pv = rhs * pvstars(T)
        rhl = pv / pvstarl(T)
        if T > Ttrip:
            rh = rhl
        else:
            rh = rhs
    if pv > p:
        return N

 # Calculate lcl_liquid and lcl_solid
    qv = rgasa*pv / (rgasv*p + (rgasa-rgasv)*pv)
    rgasm = (1-qv)*rgasa + qv*rgasv
    cpm = (1-qv)*cpa + qv*cpv
    if rh == 0:
        return cpm*T/ggr
    aL = -(cpv-cvl)/rgasv + cpm/rgasm
    bL = -(E0v-(cvv-cvl)*Ttrip)/(rgasv*T)
    cL = pv/pvstarl(T)*math.exp(-(E0v-(cvv-cvl)*Ttrip)/(rgasv*T))
    aS = -(cpv-cvs)/rgasv + cpm/rgasm
    bS = -(E0v+E0s-(cvv-cvs)*Ttrip)/(rgasv*T)
    cS = pv/pvstars(T)*math.exp(-(E0v+E0s-(cvv-cvs)*Ttrip)/(rgasv*T))
    X  = bL/(aL*scipy.special.lambertw(bL/aL*cL**(1/aL),-1).real)
    Y  = bS/(aS*scipy.special.lambertw(bS/aS*cS**(1/aS),-1).real)

    lcl = cpm*T/ggr*( 1 - X)
    ldl = cpm*T/ggr*( 1 - Y)

    # Modifications of the code to output Plcl or Pldl
    Plcl = PPa * X**(cpm/rgasm)
    Pldl = PPa * X**(cpm/rgasm)
    # Return either lcl or ldl
    if return_ldl and return_min_lcl_ldl:
        exit('return_ldl and return_min_lcl_ldl cannot both be true')
    elif return_ldl:
        return Pldl
    elif return_min_lcl_ldl:
        return min(Plcl,Pldl)
    else:
        return Plcl
 ```

 %% Cell type:code id: tags:

 ``` python
 import aes_thermo as mt
 PPa  = 101325.

 qt = np.arange(2.5e-3,8e-3,0.2e-3)
 TK = 285.
 Plcl_X = mt.plcl(TK,PPa,qt)
 Plcl_B = mt.plcl_boloton(TK,PPa,qt)
 Plcl_R = np.zeros(len(Plcl_X))

 for i,x in enumerate(qt):
    if (x>0.1): x = x/1000.
    RH        = mt.mr2pp(x/(1.-x),PPa)/mt.es(TK)
    Plcl_R[i] = lcl(PPa,TK,RH)

 del1 = (Plcl_B-Plcl_X)/100.
 del2 = (Plcl_R-Plcl_X)/100.

 fig = plt.figure(figsize=(10,5))
 ax1 = plt.subplot(1,2,1)
 ax1.set_ylabel('$P$ / hPa')
 ax1.set_xlabel('$q_\mathrm{t}$ / g/kg')
 #ax1.set_ylim(-1.2,1.2)

 plt.plot(qt*1.e3,del1,label='$\\delta_\mathrm{B}$, $T$=285K')
 plt.plot(qt*1.e3,del2,label='$\\delta_\mathrm{R}$, $T$=285K')
 #plt.gca().invert_yaxis()
 plt.legend(loc="best")

 qt = np.arange(0.5e-3,28e-3,0.2e-3)
 TK = 310.
 Plcl_X = mt.get_Plcl(TK,PPa,qt,iterate=True)
 Plcl_B = mt.get_Plcl(TK,PPa,qt)
 Plcl_R = np.zeros(len(Plcl_X))

 for i,x in enumerate(qt):
    if (x>0.1): x = x/1000.
    RH        = mt.mr2pp(x/(1.-x),PPa)/mt.es(TK)
    Plcl_R[i] = lcl(PPa,TK,RH)

 del1 = (Plcl_B-Plcl_X)/100.
 del2 = (Plcl_R-Plcl_X)/100.

 ax2 = plt.subplot(1,2,2)
 ax2.set_xlabel('$q_\mathrm{t}$ / g/kg')
 #ax2.set_ylim(-1.2,1.2)
 ax2.set_yticklabels([])

 plt.plot(qt*1.e3,del1,label='$\\delta_\mathrm{B}$, $T$=310K')
 plt.plot(qt*1.e3,del2,label='$\\delta_\mathrm{R}$, $T$=310K')
 #plt.gca().invert_yaxis()
 plt.legend(loc="best")

 sns.set_context("talk", font_scale=1.2)
 sns.despine(offset=10)
 fig.savefig(plot_dir+'Plcl.pdf')
 ```

 %% Cell type:markdown id: tags:

 ## Acknowledgments ##

 Jiawei Bao, Geet George, and Hauke Schulz are thanked for comments on these notes, and the identification of some errors in earlier versions. Axel Seifert is thanked for his comments and insights, and for pointing out the TEOS routines of Feisel et al. (2010) which extend the IAPWS libraries.

--- a/moist_thermodynamics/functions.py
+++ b/moist_thermodynamics/functions.py
@@ -8,101 +8,17 @@ copygright, bjorn stevens Max Planck Institute for Meteorology, Hamburg
 License: BSD-3C
 """
 #
-from . import constants
 import numpy as np
 from scipy import interpolate, optimize

+from . import constants
+from . import saturation_vapor_pressures

-def planck(T, nu):
-    """Planck source function (J/m2 per steradian per Hz)
-
-    Args:
-        T: temperature in kelvin
-        nu: frequency in Hz
-
-    Returns:
-        Returns the radiance in the differential frequency interval per unit steradian. Usually we
-        multiply by $\pi$ to convert to irradiances
-
-    >>> planck(300,1000*constants.c)
-    8.086837160291128e-15
-    """
-    c = constants.speed_of_light
-    h = constants.planck_constant
-    kB = constants.boltzmann_constant
-    return (2 * h * nu**3 / c**2) / (np.exp(h * nu / (kB * T)) - 1)
-
-
-def es_liq(T):
-    """Returns saturation vapor pressure (Pa) over planer liquid water
-
-    Encodes the empirical fits of Wagner and Pruss (2002), Eq 2.5a (page 399). Their formulation
-    is compared to other fits in the example scripts used in this package, and deemed to be the
-    best reference.
-
-    Args:
-        T: temperature in kelvin
-
-    Reference:
-        W. Wagner and A. Pruß , "The IAPWS Formulation 1995 for the Thermodynamic Properties
-    of Ordinary Water Substance for General and Scientific Use", Journal of Physical and Chemical
-    Reference Data 31, 387-535 (2002) https://doi.org/10.1063/1.1461829
-
-    >>> es_liq(np.asarray([273.16,305.]))
-    array([ 611.65706974, 4719.32683147])
-    """
-    TvC = constants.temperature_water_vapor_critical_point
-    PvC = constants.pressure_water_vapor_critical_point
-
-    vt = 1.0 - T / TvC
-    es = PvC * np.exp(
-        TvC
-        / T
-        * (
-            -7.85951783 * vt
-            + 1.84408259 * vt**1.5
-            - 11.7866497 * vt**3
-            + 22.6807411 * vt**3.5
-            - 15.9618719 * vt**4
-            + 1.80122502 * vt**7.5
-        )
-    )
-    return es
-
-
-def es_ice(T):
-    """Returns sublimation vapor pressure (Pa) over simple (Ih) ice
-
-    Encodes the emperical fits of Wagner et al., (2011) which also define the IAPWS standard for
-    sublimation vapor pressure over ice-Ih
-
-    Args:
-        T: temperature in kelvin
-
-    Reference:
-        Wagner, W., Riethmann, T., Feistel, R. & Harvey, A. H. New Equations for the Sublimation
-        Pressure and Melting Pressure of H 2 O Ice Ih. Journal of Physical and Chemical Reference
-        Data 40, 043103 (2011).
-
-
-    >>> es_ice(np.asarray([273.16,260.]))
-    array([611.655     , 195.80103377])
-    """
-    TvT = constants.temperature_water_vapor_triple_point
-    PvT = constants.pressure_water_vapor_triple_point
-
-    a1 = -0.212144006e2
-    a2 = 0.273203819e2
-    a3 = -0.610598130e1
-    b1 = 0.333333333e-2
-    b2 = 0.120666667e1
-    b3 = 0.170333333e1
-    theta = T / TvT
-    es = PvT * np.exp((a1 * theta**b1 + a2 * theta**b2 + a3 * theta**b3) / theta)
-    return es
+es_liq_default = saturation_vapor_pressures.liq_wagner_pruss
+es_ice_default = saturation_vapor_pressures.ice_wagner_etal


-def es_mxd(T):
+def es_mxd(T, es_liq=es_liq_default, es_ice=es_ice_default):
    """Returns the minimum of the sublimation and saturation vapor pressure

    Calculates both the sublimation vapor pressure over ice Ih using es_ice and that over planar
@@ -120,152 +36,27 @@ def es_mxd(T):
    return np.minimum(es_liq(T), es_ice(T))


-def es_liq_murphykoop(T):
-    """Returns saturation vapor pressure (Pa) over liquid water
-
-    Encodes the empirical fit (Eq. 10) of Murphy and Koop (2011) which improves on the Wagner and
-    Pruß fits for supercooled conditions.
-
-    Args:
-        T: temperature in kelvin
-
-    Reference:
-        Murphy, D. M. & Koop, T. Review of the vapour pressures of ice and supercooled water for
-        atmospheric applications. Q. J. R. Meteorol. Soc. 131, 1539–1565 (2005).
-
-    >>> es_liq_murphykoop(np.asarray([273.16,140.]))
-    array([6.11657044e+02, 9.39696372e-07])
-    """
-
-    X = np.tanh(0.0415 * (T - 218.8)) * (
-        53.878 - 1331.22 / T - 9.44523 * np.log(T) + 0.014025 * T
-    )
-    return np.exp(54.842763 - 6763.22 / T - 4.210 * np.log(T) + 0.000367 * T + X)
-
-
-def es_liq_hardy(T):
-    """Returns satruation vapor pressure (Pa) over liquid water
-
-    Encodes the empirical fit (Eq. 10) of Hardy (1998) which is often used in the postprocessing
-    of radiosondes
-
-    Args:
-        T: temperature in kelvin
-
-    Reference:
-        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
-
-    >>> es_liq_hardy(np.asarray([273.16,260.]))
-    array([611.65715494, 222.65143353])
-    """
-    X = (
-        -2.8365744e3 / (T * T)
-        - 6.028076559e3 / T
-        + 19.54263612
-        - 2.737830188e-2 * T
-        + 1.6261698e-5 * T**2
-        + 7.0229056e-10 * T**3
-        - 1.8680009e-13 * T**4
-        + 2.7150305 * np.log(T)
-    )
-    return np.exp(X)
-
-
-def es_liq_analytic(T, delta_cl=constants.delta_cl):
-    """Analytic approximation for saturation vapor pressure over iquid
-
-    Uses the rankine (constant specific heat, negligible condensate volume) approximations to
-    calculate the saturation vapor pressure over liquid.  The procedure is described in Eq(4) of
-    Romps (2017) and best approximates the actual value for specific heats that differ slightly
-    from the best estimates of these quantities which are provided as default quantities.
-    Romps recommends cl = 4119 J/kg/K, and cpv = 1861 J/kg/K.
-
-    Args:
-        T: temperature in kelvin
-        delta_cl: differnce between isobaric specific heat capacity of vapor and that of liquid.
-
-    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.
-
-    >>> es_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 = delta_cl / Rv
-    c2 = lvT / (Rv * TvT) - c1
-    es = PvT * np.exp(c2 * (1.0 - TvT / T)) * (T / TvT) ** c1
-    return es
-
-
-def es_ice_analytic(T, delta_ci=constants.delta_ci):
-    """Analytic approximation for saturation vapor pressure over ice
-
-    Uses the rankine (constant specific heat, negligible condensate volume) approximations to
-    calculate the saturation vapor pressure over ice.  The procedure is described in Eq(4) of
-    Romps (2017) and best approximates the actual value for specific heats that differ slightly
-    from the best estimates of these quantities which are provided as default quantities.
-    Romps recommends ci = 1861 J/kg/K, and cpv = 1879 J/kg/K.
-
-    Args:
-        T: temperature in kelvin
-        delta_cl: differnce between isobaric specific heat capacity of vapor and that of liquid.
-
-    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.
-
-
-    >>> es_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
-
-    c1 = delta_ci / Rv
-    c2 = lsT / (Rv * TvT) - c1
-    es = PvT * np.exp(c2 * (1.0 - TvT / T)) * (T / TvT) ** c1
-    return es
-
-
-def es_mxd_analytic(T, delta_cl=constants.delta_cl, delta_ci=constants.delta_ci):
-    """Returns the minimum of the analytic sublimation and saturation vapor pressure
-
-    Calculates both the sublimation vapor pressure over ice Ih using es_ice_analytic and
-    that over planar water using es_liq_analytic, and returns the minimum of the two
-    quantities.
+def planck(T, nu):
+    """Planck source function (J/m2 per steradian per Hz)

    Args:
        T: temperature in kelvin
+        nu: frequency in Hz

    Returns:
-        value of es_ice_analytic(T) for T < 273.15 and es_liq_analytic(T) otherwise
+        Returns the radiance in the differential frequency interval per unit steradian. Usually we
+        multiply by $\pi$ to convert to irradiances

-    >>> es_ice_analytic(np.asarray([273.16,260.]))
-    array([611.655     , 195.99959431])
+    >>> planck(300,1000*constants.c)
+    8.086837160291128e-15
    """
-    return np.minimum(es_liq_analytic(T, delta_cl), es_ice_analytic(T, delta_ci))
+    c = constants.speed_of_light
+    h = constants.planck_constant
+    kB = constants.boltzmann_constant
+    return (2 * h * nu**3 / c**2) / (np.exp(h * nu / (kB * T)) - 1)


-def vaporization_enthalpy(TK, delta_cl=constants.delta_cl):
+def vaporization_enthalpy(T, delta_cl=constants.delta_cl):
    """Returns the vaporization enthlapy of water (J/kg)

    The vaporization enthalpy is calculated from a linear depdence on temperature about a
@@ -281,10 +72,10 @@ def vaporization_enthalpy(TK, delta_cl=constants.delta_cl):
    """
    T0 = constants.standard_temperature
    lv0 = constants.vaporization_enthalpy_stp
-    return lv0 + delta_cl * (TK - T0)
+    return lv0 + delta_cl * (T - T0)


-def sublimation_enthalpy(TK, delta_ci=constants.delta_ci):
+def sublimation_enthalpy(T, delta_ci=constants.delta_ci):
    """Returns the sublimation enthlapy of water (J/kg)

    The sublimation enthalpy is calculated from a linear depdence on temperature about a
@@ -301,13 +92,13 @@ def sublimation_enthalpy(TK, delta_ci=constants.delta_ci):
    """
    T0 = constants.standard_temperature
    ls0 = constants.sublimation_enthalpy_stp
-    return ls0 + delta_ci * (TK - T0)
+    return ls0 + delta_ci * (T - T0)


 def partial_pressure_to_mixing_ratio(pp, p):
    """Returns the mass mixing ratio given the partial pressure and pressure

-    >>> partial_pressure_to_mixing_ratio(es_liq(300.),60000.)
+    >>> partial_pressure_to_mixing_ratio(es_liq_default(300.),60000.)
    0.0389569254590098
    """
    eps1 = constants.rd_over_rv
@@ -338,7 +129,7 @@ def partial_pressure_to_specific_humidity(pp, p):
    situations where condensate is present one should instead calculate
    $q = r*(1-qt)$ which would require an additional argument

-    >>> partial_pressure_to_specific_humidity(es_liq(300.),60000.)
+    >>> partial_pressure_to_specific_humidity(es_liq_default(300.),60000.)
    0.037496189210922945
    """
    r = partial_pressure_to_mixing_ratio(pp, p)
@@ -358,31 +149,60 @@ def saturation_partition(P, ps, qt):
    return np.minimum(qt, qs)


-def static_energy(TK, qv=0.0, ql=0.0, qi=0.0):
-    """Returns the static energy, defaulting to that for a dry atmosphere
+def static_energy(T, Z, qv=0, ql=0, qi=0, hv0=constants.cpv * constants.T0):
+    """Returns the static energy

-    The moist static energy is calculated so that it includes the effects of composition
-    on the specific heat if specific humidities are included, but defaults to the dry static energy
+    The static energy is calculated so that it includes the effects of composition on the
+    specific heat if specific humidities are included.  Different common forms of the static
+    energy arise from different choices of the reference state and condensate loading:
+        - hv0 = cpv*T0      -> frozen, liquid moist static energy
+        - hv0 = ls0 + ci*T0 -> frozen moist static energy
+        - hv0 = cpv*T0      -> liquid water static energy if qi= 0 (default if qv /= 0)
+        - hv0 = lv0 + cl*T0 -> moist static energy if qi= 0.
+        - qv=ql=q0=0        -> dry static energy (default)
+
+    Because the composition weights the reference enthalpies, different choices do not differ by
+    a constant, but rather by a constant weighted by the specific masses of the different water
+    phases.

    Args:
-        TK: temperature in kelvin
-        PPa: pressure in pascal
+        T: temperature in kelvin
+        Z: altitude (above mean sea-level) in meters
        qv: specific vapor mass
        ql: specific liquid mass
        qi: specific ice mass
+        hv0: reference vapor enthalpy
+
+        >>> static_energy(300.,600.,15.e-3,hv0=constants.lv0 + constants.cl * constants.T0)
+        358162.78621841426
+
    """
    cpd = constants.isobaric_dry_air_specific_heat
    cpv = constants.isobaric_water_vapor_specific_heat
    cl = constants.liquid_water_specific_heat
    ci = constants.frozen_water_specific_heat
+    lv0 = constants.lv0
+    ls0 = constants.ls0
+    T0 = constants.T0
    g = constants.gravity_earth

-    qt = qv + ql + qi
-    cp = cpd + qt * (cl - cpd)
-    return TK * cp + qv * lv + gz
+    qd = 1.0 - qv - ql - qi
+    cp = qd * cpd + qv * cpv + ql * cl + qi * ci
+
+    h = (
+        qd * cpd * T
+        + qv * cpv * T
+        + ql * cl * T
+        + qi * ci * T
+        + qv * (hv0 - cpv * T0)
+        + ql * (hv0 - lv0 - cl * T0)
+        + qi * (hv0 - ls0 - ci * T0)
+        + g * Z
+    )
+    return h


-def theta(TK, PPa, qv=0.0, ql=0.0, qi=0.0):
+def theta(T, P, qv=0.0, ql=0.0, qi=0.0):
    """Returns the potential temperature for an unsaturated moist fluid

    This expressed the potential temperature in away that makes it possible to account
@@ -390,8 +210,8 @@ def theta(TK, PPa, qv=0.0, ql=0.0, qi=0.0):
    adiabatic factor R/cp.  The default is the usualy dry potential temperature.

    Args:
-        TK: temperature in kelvin
-        PPa: pressure in pascal
+        T: temperature in kelvin
+        P: pressure in pascal
        qv: specific vapor mass
        ql: specific liquid mass
        qi: specific ice mass
@@ -407,18 +227,18 @@ def theta(TK, PPa, qv=0.0, ql=0.0, qi=0.0):

    qd = 1.0 - qv - ql - qi
    kappa = (qd * Rd + qv * Rv) / (qd * cpd + qv * cpv + ql * cl + qi * ci)
-    return TK * (P0 / PPa) ** kappa
+    return T * (P0 / P) ** kappa


-def theta_e_bolton(TK, PPa, qt, es=es_liq):
+def theta_e_bolton(T, P, qt, es=es_liq_default):
    """Returns the pseudo equivalent potential temperature.

    Following Eq. 43 in Bolton (1980) the (pseudo) equivalent potential temperature
    is calculated and returned by this function

    Args:
-        TK: temperature in kelvin
-        PPa: pressure in pascal
+        T: temperature in kelvin
+        P: pressure in pascal
        qt: specific total water mass
        es: form of the saturation vapor pressure to use

@@ -431,19 +251,19 @@ def theta_e_bolton(TK, PPa, qt, es=es_liq):
    r2p = mixing_ratio_to_partial_pressure

    rv = np.minimum(
-        qt / (1.0 - qt), p2r(es(TK), PPa)
+        qt / (1.0 - qt), p2r(es(T), P)
    )  # mixing ratio of vapor (not gas Rv)
-    pv = r2p(rv, PPa)
+    pv = r2p(rv, P)

-    TL = 55.0 + 2840.0 / (3.5 * np.log(TK) - np.log(pv / 100.0) - 4.805)
+    TL = 55.0 + 2840.0 / (3.5 * np.log(T) - np.log(pv / 100.0) - 4.805)
    return (
-        TK
-        * (P0 / PPa) ** (0.2854 * (1.0 - 0.28 * rv))
+        T
+        * (P0 / P) ** (0.2854 * (1.0 - 0.28 * rv))
        * np.exp((3376.0 / TL - 2.54) * rv * (1 + 0.81 * rv))
    )


-def theta_e(TK, PPa, qt, es=es_liq):
+def theta_e(T, P, qt, es=es_liq_default):
    """Returns the equivalent potential temperature

    Follows Eq. 11 in Marquet and Stevens (2022). The closed form solutionis derived for a
@@ -452,8 +272,8 @@ def theta_e(TK, PPa, qt, es=es_liq):
    accurate, but more consistent, formulations are on the order of millikelvin

    Args:
-        TK: temperature in kelvin
-        PPa: pressure in pascal
+        T: temperature in kelvin
+        P: pressure in pascal
        qt: total water specific humidity (unitless)
        es: form of the saturation vapor pressure

@@ -469,20 +289,20 @@ def theta_e(TK, PPa, qt, es=es_liq):
    cl = constants.liquid_water_specific_heat
    lv = vaporization_enthalpy

-    ps = es(TK)
-    qv = saturation_partition(PPa, ps, qt)
+    ps = es(T)
+    qv = saturation_partition(P, ps, qt)

    Re = (1.0 - qt) * Rd
    R = Re + qv * Rv
-    pv = qv * (Rv / R) * PPa
+    pv = qv * (Rv / R) * P
    RH = pv / ps
    cpe = cpd + qt * (cl - cpd)
    omega_e = RH ** (-qv * Rv / cpe) * (R / Re) ** (Re / cpe)
-    theta_e = TK * (P0 / PPa) ** (Re / cpe) * omega_e * np.exp(qv * lv(TK) / (cpe * TK))
+    theta_e = T * (P0 / P) ** (Re / cpe) * omega_e * np.exp(qv * lv(T) / (cpe * T))
    return theta_e


-def theta_l(TK, PPa, qt, es=es_liq):
+def theta_l(T, P, qt, es=es_liq_default):
    """Returns the liquid-water potential temperature

    Follows Eq. 16 in Marquet and Stevens (2022). The closed form solutionis derived for a
@@ -491,8 +311,8 @@ def theta_l(TK, PPa, qt, es=es_liq):
    accurate, but more consistent, formulations are on the order of millikelvin

    Args:
-        TK: temperature in kelvin
-        PPa: pressure in pascal
+        T: temperature in kelvin
+        P: pressure in pascal
        qt: total water specific humidity (unitless)
        es: form of the saturation vapor pressure

@@ -508,8 +328,8 @@ def theta_l(TK, PPa, qt, es=es_liq):
    cpv = constants.isobaric_water_vapor_specific_heat
    lv = vaporization_enthalpy

-    ps = es(TK)
-    qv = saturation_partition(PPa, ps, qt)
+    ps = es(T)
+    qv = saturation_partition(P, ps, qt)
    ql = qt - qv

    R = Rd * (1 - qt) + qv * Rv
@@ -517,13 +337,11 @@ def theta_l(TK, PPa, qt, es=es_liq):
    cpl = cpd + qt * (cpv - cpd)

    omega_l = (R / Rl) ** (Rl / cpl) * (qt / (qv + 1.0e-15)) ** (qt * Rv / cpl)
-    theta_l = (
-        (TK * (P0 / PPa) ** (Rl / cpl)) * omega_l * np.exp(-ql * lv(TK) / (cpl * TK))
-    )
+    theta_l = (T * (P0 / P) ** (Rl / cpl)) * omega_l * np.exp(-ql * lv(T) / (cpl * T))
    return theta_l


-def theta_s(TK, PPa, qt, es=es_liq):
+def theta_s(T, P, qt, es=es_liq_default):
    """Returns the entropy potential temperature

    Follows Eq. 18 in Marquet and Stevens (2022). The closed form solutionis derived for a
@@ -532,8 +350,8 @@ def theta_s(TK, PPa, qt, es=es_liq):
    accurate, but more consistent, formulations are on the order of millikelvin

    Args:
-        TK: temperature in kelvin
-        PPa: pressure in pascal
+        T: temperature in kelvin
+        P: pressure in pascal
        qt: total water specific humidity (unitless)
        es: form of the saturation vapor pressure

@@ -566,18 +384,18 @@ def theta_s(TK, PPa, qt, es=es_liq):
    gamma = kappa / eps1
    r0 = e0 / (P0 - e0) / eta

-    ps = es(TK)
-    qv = saturation_partition(PPa, ps, qt)
+    ps = es(T)
+    qv = saturation_partition(P, ps, qt)
    ql = qt - qv

    R = Rd + qv * (Rv - Rd)
-    pv = qv * (Rv / R) * PPa
+    pv = qv * (Rv / R) * P
    RH = pv / ps
    rv = qv / (1 - qv)

    x1 = (
-        (TK / T0) ** (lmbd * qt)
-        * (P0 / PPa) ** (kappa * delta * qt)
+        (T / T0) ** (lmbd * qt)
+        * (P0 / P) ** (kappa * delta * qt)
        * (rv / r0) ** (-gamma * qt)
        * RH ** (gamma * ql)
    )
@@ -585,8 +403,8 @@ def theta_s(TK, PPa, qt, es=es_liq):
        -kappa * delta * qt
    )
    theta_s = (
-        (TK * (P0 / PPa) ** (kappa))
-        * np.exp(-ql * lv(TK) / (cpd * TK))
+        (T * (P0 / P) ** (kappa))
+        * np.exp(-ql * lv(T) / (cpd * T))
        * np.exp(qt * Lmbd)
        * x1
        * x2
@@ -594,15 +412,15 @@ def theta_s(TK, PPa, qt, es=es_liq):
    return theta_s


-def theta_es(TK, PPa, es=es_liq):
+def theta_es(T, P, es=es_liq_default):
    """Returns the saturated equivalent potential temperature

    Adapted from Eq. 11 in Marquet and Stevens (2022) with the assumption that the gas quanta is
    everywhere just saturated.

    Args:
-        TK: temperature in kelvin
-        PPa: pressure in pascal
+        T: temperature in kelvin
+        P: pressure in pascal
        qt: total water specific humidity (unitless)
        es: form of the saturation vapor pressure

@@ -618,20 +436,18 @@ def theta_es(TK, PPa, es=es_liq):
    p2q = partial_pressure_to_specific_humidity
    lv = vaporization_enthalpy

-    ps = es(TK)
-    qs = p2q(ps, PPa)
+    ps = es(T)
+    qs = p2q(ps, P)

    Re = (1.0 - qs) * Rd
    R = Re + qs * Rv
    cpe = cpd + qs * (cl - cpd)
    omega_e = (R / Re) ** (Re / cpe)
-    theta_es = (
-        TK * (P0 / PPa) ** (Re / cpe) * omega_e * np.exp(qs * lv(TK) / (cpe * TK))
-    )
+    theta_es = T * (P0 / P) ** (Re / cpe) * omega_e * np.exp(qs * lv(T) / (cpe * T))
    return theta_es


-def theta_rho(TK, PPa, qt, es=es_liq):
+def theta_rho(T, P, qt, es=es_liq_default):
    """Returns the density liquid-water potential temperature

    calculates $\theta_\mathrm{l} R/R_\mathrm{d}$ where $R$ is the gas constant of a
@@ -639,21 +455,21 @@ def theta_rho(TK, PPa, qt, es=es_liq):
    temperature baswed on the two component fluid thermodynamic constants.

    Args:
-        TK: temperature in kelvin
-        PPa: pressure in pascal
+        T: temperature in kelvin
+        P: pressure in pascal
        qt: total water specific humidity (unitless)
        es: form of the saturation vapor pressure
    """
    Rd = constants.dry_air_gas_constant
    Rv = constants.water_vapor_gas_constant

-    ps = es(TK)
-    qv = saturation_partition(PPa, ps, qt)
-    theta_rho = theta_l(TK, PPa, qt, es) * (1.0 - qt + qv * Rv / Rd)
+    ps = es(T)
+    qv = saturation_partition(P, ps, qt)
+    theta_rho = theta_l(T, P, qt, es) * (1.0 - qt + qv * Rv / Rd)
    return theta_rho


-def invert_for_temperature(f, f_val, P, qt, es=es_liq):
+def invert_for_temperature(f, f_val, P, qt, es=es_liq_default):
    """Returns temperature for an atmosphere whose state is given by f, P and qt

        Infers the temperature from a state description (f,P,qt), where
@@ -678,7 +494,7 @@ def invert_for_temperature(f, f_val, P, qt, es=es_liq):
    return optimize.newton(zero, 280.0, args=(f_val,))


-def invert_for_pressure(f, f_val, T, qt, es=es_liq):
+def invert_for_pressure(f, f_val, T, qt, es=es_liq_default):
    """Returns pressure for an atmosphere whose state is given by f, T and qt

        Infers the pressure from a state description (f,T,qt), where
@@ -703,7 +519,7 @@ def invert_for_pressure(f, f_val, T, qt, es=es_liq):
    return optimize.newton(zero, 80000.0, args=(f_val,))


-def plcl(TK, PPa, qt, es=es_liq):
+def plcl(T, P, qt, es=es_liq_default):
    """Returns the pressure at the lifting condensation level

    Calculates the lifting condensation level pressure using an interative solution under the
@@ -711,8 +527,8 @@ def plcl(TK, PPa, qt, es=es_liq):
    which depends on the expression for the saturation vapor pressure

    Args:
-        TK: temperature in kelvin
-        PPa: pressure in pascal
+        T: temperature in kelvin
+        P: pressure in pascal
        qt: specific total water mass

        >>> plcl(300.,102000.,17e-3)
@@ -721,23 +537,23 @@ def plcl(TK, PPa, qt, es=es_liq):

    def zero(P, Tl):
        p2r = partial_pressure_to_mixing_ratio
-        TK = invert_for_temperature(theta_l, Tl, P, qt, es=es)
-        qs = p2r(es(TK), P) * (1.0 - qt)
+        T = invert_for_temperature(theta_l, Tl, P, qt, es=es)
+        qs = p2r(es(T), P) * (1.0 - qt)
        return np.abs(qs / qt - 1.0)

-    Tl = theta_l(TK, PPa, qt, es=es)
+    Tl = theta_l(T, P, qt, es=es)
    return optimize.fsolve(zero, 80000.0, args=(Tl,))


-def plcl_bolton(TK, PPa, qt):
+def plcl_bolton(T, P, qt):
    """Returns the pressure at the lifting condensation level

    Following Bolton (1980) the lifting condensation level pressure is derived from the state
    of an air parcel.  Usually accurate to within about 10 Pa, or about 1 m

    Args:
-        TK: temperature in kelvin
-        PPa: pressure in pascal
+        T: temperature in kelvin
+        P: pressure in pascal
        qt: specific total water mass

    Reference:
@@ -755,9 +571,9 @@ def plcl_bolton(TK, PPa, qt):

    cp = cpd + qt * (cpv - cpd)
    R = Rd + qt * (Rv - Rd)
-    pv = r2p(qt / (1.0 - qt), PPa)
-    Tl = 55 + 2840.0 / (3.5 * np.log(TK) - np.log(pv / 100.0) - 4.805)
-    return PPa * (Tl / TK) ** (cp / R)
+    pv = r2p(qt / (1.0 - qt), P)
+    Tl = 55 + 2840.0 / (3.5 * np.log(T) - np.log(pv / 100.0) - 4.805)
+    return P * (Tl / T) ** (cp / R)


 def zlcl(Plcl, T, P, qt, z):
@@ -793,7 +609,14 @@ from scipy.integrate import ode


 def moist_adiabat(
-    Tbeg, Pbeg, Pend, dP, qt, cc=constants.cl, l=vaporization_enthalpy, es=es_liq
+    Tbeg,
+    Pbeg,
+    Pend,
+    dP,
+    qt,
+    cc=constants.cl,
+    l=vaporization_enthalpy,
+    es=es_liq_default,
 ):
    """Returns the temperature and pressure by integrating along a moist adiabat


--- a/moist_thermodynamics/saturation_vapor_pressures.py
+++ b/moist_thermodynamics/saturation_vapor_pressures.py
+# -*- coding: utf-8 -*-
+"""
+Provides a collection of fits for saturation and sublimation vapor pressure
+
+Author: Bjorn Stevens (bjorn.stevens@mpimet.mpg.de)
+copygright, bjorn stevens Max Planck Institute for Meteorology, Hamburg
+
+License: BSD-3C
+"""
+#
+from . import constants
+import numpy as np
+
+
+def liq_wagner_pruss(T):
+    """Returns saturation vapor pressure (Pa) over planer liquid water
+
+    Encodes the empirical fits of Wagner and Pruss (2002), Eq 2.5a (page 399). Their formulation
+    is compared to other fits in the example scripts used in this package, and deemed to be the
+    best reference.
+
+    The fit has been verified for TvT <= T < = TvC.  For super cooled water (T<TvT) it deviates
+    from the results of Murphy and Koop where were developed for super-cooled water.  It is about
+    10% larger at 200 K, 25 % larter at 150 K, and then decreases again so it is 12% smaller at
+    the limit (123K) of the Murphy and Koop fit.  For accurate fits for super-cooled water the
+    function of Murphy and Koop should be used.
+
+    Args:
+        T: temperature in kelvin
+
+    Reference:
+        W. Wagner and A. Pruß , "The IAPWS Formulation 1995 for the Thermodynamic Properties
+    of Ordinary Water Substance for General and Scientific Use", Journal of Physical and Chemical
+    Reference Data 31, 387-535 (2002) https://doi.org/10.1063/1.1461829
+
+    >>> liq_wagner_pruss(np.asarray([273.16,305.]))
+    array([ 611.65706974, 4719.32683147])
+    """
+    TvC = constants.temperature_water_vapor_critical_point
+    PvC = constants.pressure_water_vapor_critical_point
+
+    vt = 1.0 - T / TvC
+    es = PvC * np.exp(
+        TvC
+        / T
+        * (
+            -7.85951783 * vt
+            + 1.84408259 * vt**1.5
+            - 11.7866497 * vt**3
+            + 22.6807411 * vt**3.5
+            - 15.9618719 * vt**4
+            + 1.80122502 * vt**7.5
+        )
+    )
+    return es
+
+
+def ice_wagner_etal(T):
+    """Returns sublimation vapor pressure (Pa) over simple (Ih) ice
+
+    Encodes the emperical fits of Wagner et al., (2011) which also define the IAPWS standard for
+    sublimation vapor pressure over ice-Ih
+
+    Args:
+        T: temperature in kelvin
+
+    Reference:
+        Wagner, W., Riethmann, T., Feistel, R. & Harvey, A. H. New Equations for the Sublimation
+        Pressure and Melting Pressure of H 2 O Ice Ih. Journal of Physical and Chemical Reference
+        Data 40, 043103 (2011).
+
+
+    >>> ice_wagner_etal(np.asarray([273.16,260.]))
+    array([611.655     , 195.80103377])
+    """
+    TvT = constants.temperature_water_vapor_triple_point
+    PvT = constants.pressure_water_vapor_triple_point
+
+    a1 = -0.212144006e2
+    a2 = 0.273203819e2
+    a3 = -0.610598130e1
+    b1 = 0.333333333e-2
+    b2 = 0.120666667e1
+    b3 = 0.170333333e1
+    theta = T / TvT
+    es = PvT * np.exp((a1 * theta**b1 + a2 * theta**b2 + a3 * theta**b3) / theta)
+    return es
+
+
+def liq_murphy_koop(T):
+    """Returns saturation vapor pressure (Pa) over liquid water
+
+    Encodes the empirical fit (Eq. 10) of Murphy and Koop (2011) which improves on the Wagner and
+    Pruß fits for supercooled conditions.
+
+    The fit has been verified for 123K <= T < = 332 K
+
+    Args:
+        T: temperature in kelvin
+
+    Reference:
+        Murphy, D. M. & Koop, T. Review of the vapour pressures of ice and supercooled water for
+        atmospheric applications. Q. J. R. Meteorol. Soc. 131, 1539–1565 (2005).
+
+    >>> liq_murphy_koop(np.asarray([273.16,140.]))
+    array([6.11657044e+02, 9.39696372e-07])
+    """
+
+    X = np.tanh(0.0415 * (T - 218.8)) * (
+        53.878 - 1331.22 / T - 9.44523 * np.log(T) + 0.014025 * T
+    )
+    return np.exp(54.842763 - 6763.22 / T - 4.210 * np.log(T) + 0.000367 * T + X)
+
+
+def liq_hardy(T):
+    """Returns satruation vapor pressure (Pa) over liquid water
+
+    Encodes the empirical fit (Eq. 10) of Hardy (1998) which is often used in the postprocessing
+    of radiosondes
+
+    Args:
+        T: temperature in kelvin
+
+    Reference:
+        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
+
+    >>> liq_hardy(np.asarray([273.16,260.]))
+    array([611.65715494, 222.65143353])
+    """
+    X = (
+        -2.8365744e3 / (T * T)
+        - 6.028076559e3 / T
+        + 19.54263612
+        - 2.737830188e-2 * T
+        + 1.6261698e-5 * T**2
+        + 7.0229056e-10 * T**3
+        - 1.8680009e-13 * T**4
+        + 2.7150305 * np.log(T)
+    )
+    return np.exp(X)
+
+
+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
+    calculate the saturation vapor pressure over liquid.  The procedure is described in Eq(4) of
+    Romps (2017) and best approximates the actual value for specific heats that differ slightly
+    from the best estimates of these quantities which are provided as default quantities.
+    Romps recommends cl = 4119 J/kg/K, and cpv = 1861 J/kg/K.
+
+    Args:
+        T: temperature in kelvin
+        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
+
+    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.
+
+    >>> liq_analytic(np.asarray([273.16,305.]))
+    array([ 611.655     , 4711.13161169])
+    """
+    return analytic(T, lx, cx)
+
+
+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
+    calculate the saturation vapor pressure over ice.  The procedure is described in Eq(4) of
+    Romps (2017) and best approximates the actual value for specific heats that differ slightly
+    from the best estimates of these quantities which are provided as default quantities.
+    Romps recommends ci = 1861 J/kg/K, and cpv = 1879 J/kg/K.
+
+    Args:
+        T: temperature in kelvin
+        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 ice 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.
+
+
+    >>> ice_analytic(np.asarray([273.16,260.]))
+    array([611.655     , 195.99959431])
+    """
+    return analytic(T, lx, cx)
+
+
+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,
+    in which case PvT and TvT need to be substituted by Pv0 and T0.
+
+    Args:
+        T: temperature in kelvin
+
+    >>> tetens(285.,17.269,35.86)
+    1389.7114123472836
+    """
+
+    es = constants.PvT * np.exp(a * (T - constants.TvT) / (T - b))
+    return es
+
+
+def liq_tetens(T):
+    """Returns saturation vapor pressure over liquid using the Magnus-Teten's formula
+
+    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
+    revised constants (a,b) chosen to better agree with the fits of Wagner and Pruss
+
+    Args:
+        T: temperature in kelvin
+
+    Reference:
+        Murray, F. W. On the Computation of Saturation Vapor Pressure. Journal of Applied Meteorology
+        and Climatology 6, 203–204 (1967).
+
+    >>> liq_tetens(np.asarray([273.16,305.]))
+    array([ 611.655     , 4719.73680592])
+    """
+    a = 17.41463775
+    b = 33.6393413
+
+    return tetens(T, a, b)
+
+
+def ice_tetens(T):
+    """Returns saturation vapor pressure over liquid using the Magnus-Teten's formula
+
+    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
+    revised constants (a,b) chosen to better agree with the fits of Wagner and Pruss
+
+    Args:
+        T: temperature in kelvin
+
+    Reference:
+        Murray, F. W. On the Computation of Saturation Vapor Pressure. Journal of Applied Meteorology
+        and Climatology 6, 203–204 (1967).
+
+    >>> ice_tetens(np.asarray([273.16,260.]))
+    array([611.655     , 196.10072658])
+    """
+    a = 22.0419977
+    b = 5.0
+
+    return tetens(T, a, b)
--- 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.5",
    description="Constants and functions for the treatment of moist atmospheric thermodynamics",
    packages=find_packages(),
    install_requires=[
No results found