xmet/lib/xmet/thermo.py

import math

from typing import Callable

from xmet.series import Series

LAPSE_RATE_DRY   = 9.8 / 1000 # degrees C per km
LAPSE_RATE_MOIST = 4.0 / 1000

PRESSURE_MAX  = 1050 # millibar
PRESSURE_MSL  = 1013.25
PRESSURE_MIN  =  100
PRESSURE_STEP =   50

def kelvin(c: float) -> float:
    return 273.15 + c

def celsius(k: float) -> float:
    return k - 273.15

def vapor_pressure(dewpoint: float) -> float:
    """
    Return the pressure, in millibar, of vapor in a parcel of a given
    dewpoint.
    """
    return 6.11 * 10 ** ((7.5 * dewpoint) / (237.3 + dewpoint))

def mixing_ratio(dewpoint: float, pressure: float) -> float:
    """
    Return the amount, in kilograms, of water vapor versus dry air in a parcel
    of a given dewpoint and pressure.
    """
    e = vapor_pressure(dewpoint)

    return (0.62197 * e) / (pressure - e)

def saturated_mixing_ratio(temp: float, pressure: float) -> float:
    """
    Return the maximum amount, in kilograms, of water vapor a parcel of a
    given temperature and pressure can hold.
    """
    es = vapor_pressure(temp)

    return (0.62197 * es) / (pressure - es)

def mixing_ratio_temp(ratio: float, pressure: float) -> float:
    """
    Return the temperature, in degrees celsius, of a given mixing ratio of
    kilograms of moisture to kilograms of air, at a specified pressure, in
    millibar.
    """
    c1 =  0.0498646455
    c2 =  2.4082965
    c3 =  7.07475
    c4 = 38.9114
    c5 =  0.0915
    c6 =  1.2035

    w = ratio * 1000.0
    x = math.log10(w * pressure / (622.0 + w))

    return celsius(math.pow(10.0, ((c1 * x) + c2)) - c3 + \
             (c4 * math.pow((math.pow(10, (c5 * x)) - c6), 2)))

def lcl(temp: float, dewpoint: float) -> float:
    """
    Return the height, in meters, at which a parcel of the given temperature
    is cooled to the given dewpoint.
    """
    return (temp - dewpoint) / 0.008

def pressure_height(pressure: float, surface: float=PRESSURE_MSL) -> float:
    """
    Return the approximate altitude, in meters, for a given pressure in
    millibar.
    """
    return (1 - (pressure / surface) ** 0.190284) * 145366.45 * 0.3048

def lapse(temp: float, delta: float, rate=LAPSE_RATE_DRY) -> float:
    """
    Return the temperature of a parcel cooled at the dry lapse rate for a
    given increase in height (in meters).
    """
    return temp - (rate * delta)

def moist_lapse_rate(temp: float, pressure: float) -> float:
    g   = 9.8076
    Hv  = 2501000
    Rsd = 287
    Rsw = 461.5
    Cpd = 1003.5
    r   = saturated_mixing_ratio(temp, pressure)
    T   = kelvin(temp)

    return g * (1 + (Hv * r) / (Rsd * T)) \
         / (Cpd + (((Hv**2) * r) / (Rsw * (T**2))))

def loft_parcel(start_temp: float,
                start_pressure: float,
                lapse_rate: Callable,
                step: float=1.0):
    """
    Loft a parcel of air from a given pressure, at a given temperature,
    yielding a Tuple containing the temperature and pressure of that parcel
    at each increment of elevation.  A Callable which, given a temperature,
    dewpoint, pressure value, will be dispatched to obtain the lapse rate
    for each height.
    """
    last_height  = None

    temp     = start_temp
    pressure = start_pressure

    yield temp, pressure

    while pressure >= PRESSURE_MIN:
        height = pressure_height(pressure)

        if last_height is not None:
            try:
                rate = lapse_rate(temp, pressure-step/2)
            except OverflowError:
                break

            temp = lapse(temp, height - last_height, rate)

            yield temp, pressure

        if pressure == PRESSURE_MIN:
            break
        elif pressure - step < PRESSURE_MIN:
            pressure = PRESSURE_MIN
        else:
            pressure -= step

        last_height = height

def follow_dry_adiabat(temp: float, pressure: float) -> Series:
    """
    Follow a dry adiabat starting at a given temp and pressure level, returning
    a Series object depicting the data points in descending pressure order.
    """
    series = Series()

    for level in loft_parcel(temp, pressure, lambda t, p: LAPSE_RATE_DRY):
        t2, p2 = level

        series[p2] = t2

    return series

def follow_moist_adiabat(temp: float, pressure: float) -> Series:
    """
    Follow a moist adiabat starting at a given temp and pressure level,
    returning a Series object depicting the data points in descending pressure
    order.
    """
    series = Series()

    for level in loft_parcel(temp, pressure, moist_lapse_rate):
        t2, p2 = level

        series[p2] = t2

    return series

def follow_saturated_mixing_ratio(temp: float, pressure: float, step: float=1.0) -> Series:
    """
    Follow a line of constant saturated mixing ratio calculated from a given
    temp and pressure level, returning a Series object depicting the data
    points in descending pressure order.
    """
    series = dict()

    ratio = saturated_mixing_ratio(temp, pressure)

    p2 = pressure

    while p2 >= PRESSURE_MIN:
        series[p2] = mixing_ratio_temp(ratio, p2)

        p2 -= step

    return series