add q and theta to L2

[Geet-George]

No description provided.

[Geet-George]

For the interpolation needed in L3, it is better to interpolate specific humidity and potential temperature (because they are conserved) than the relative humidity and temperature (which are the values in L2). Therefore, we first add the variables q and theta. We could use the following functions from JOANNE (with the necessary modifications):

import numpy as np
import metpy.calc as mpcalc
from metpy.units import units

def calc_saturation_pressure(temperature_K, method='hardy1998'):
    """
    Calculate saturation water vapor pressure
    
    Input
    -----
    temperature_K : array
        array of temperature in Kevlin or dew point temperature for actual vapor pressure
    method : str
        Formula used for calculating the saturation pressure
            'hardy1998' : ITS-90 Formulations for Vapor Pressure, Frostpoint Temperature,
                Dewpoint Temperature, and Enhancement Factors in the Range –100 to +100 C,
                Bob Hardy, Proceedings of the Third International Symposium on Humidity and Moisture,
                1998 (same as used in Aspen software after May 2018)
    
    Return
    ------
    e_sw : array
        saturation pressure (Pa)
    
    Examples
    --------
    >>> calc_saturation_pressure([273.15])
    array([ 611.2129107])
    
    >>> calc_saturation_pressure([273.15, 293.15, 253.15])
    array([  611.2129107 ,  2339.26239586,   125.58350529])
    """

    if method == 'hardy1998':
        g = np.empty(8)
        g[0] = -2.8365744*10**3
        g[1] = -6.028076559*10**3
        g[2] = 1.954263612*10**1
        g[3] = -2.737830188*10**(-2)
        g[4] = 1.6261698*10**(-5)
        g[5] = 7.0229056*10**(-10)
        g[6] = -1.8680009*10**(-13)
        g[7] = 2.7150305

        e_sw = np.zeros_like(temperature_K)

        for t, temp in enumerate(temperature_K):
            ln_e_sw = np.sum([g[i]*temp**(i-2) for i in range(0,7)]) + g[7]*np.log(temp)
            e_sw[t] = np.exp(ln_e_sw)
        return e_sw

def calc_q_from_rh(ds):
    """
    Input :

        ds : Dataset 

    Output :

        q : Specific humidity values
    
    Function to estimate specific humidity from the relative humidity,
    temperature and pressure in the given dataset. 
    """ 
    e_s = calc_saturation_pressure(ds.ta.values)
    w_s = mpcalc.mixing_ratio(e_s * units.Pa, ds.p.values * units.Pa).magnitude
    w = ds.rh.values * w_s
    q = w / (1 + w)

    return q


def calc_theta_from_T(dataset):
    """
    Input :

        dataset : Dataset

    Output :

        theta : Potential temperature values 

    Function to estimate potential temperature from the
    temperature and pressure in the given dataset. This function uses MetPy's
    functions to get theta:

    (i) mpcalc.potential_temperature()
    
    """
    theta = mpcalc.potential_temperature(
        dataset.p.values * units.Pa, dataset.ta.values * units.kelvin
    ).magnitude

    return theta
    

[Geet-George]

Fixed with #119