PyMPDATA_examples.Shipway_and_Hill_2012.formulae 

  1from collections import namedtuple
  2
  3import numpy as np
  4from scipy import constants
  5
  6
  7def convert_to(value, unit):
  8    value /= unit
  9
 10
 11si = namedtuple(
 12    "si",
 13    (
 14        "kg",
 15        "m",
 16        "s",
 17        "metres",
 18        "second",
 19        "um",
 20        "hPa",
 21        "micrometre",
 22        "minutes",
 23        "km",
 24        "dimensionless",
 25        "kelvin",
 26        "mg",
 27    ),
 28)(
 29    kg=1.0,
 30    m=1.0,
 31    s=1.0,
 32    metres=1.0,
 33    second=1.0,
 34    um=1e-6,
 35    hPa=100.0,
 36    micrometre=1e-6,
 37    minutes=60.0,
 38    km=1000.0,
 39    dimensionless=1.0,
 40    kelvin=1.0,
 41    mg=1e-6,
 42)
 43
 44_Mv = 0.018015
 45_Md = 0.028970
 46
 47const = namedtuple(
 48    "const",
 49    (
 50        "eps",
 51        "g",
 52        "p1000",
 53        "Rd",
 54        "Rv",
 55        "c_pd",
 56        "c_pv",
 57        "lv",
 58        "rho_l",
 59        "T0",
 60        "ARM_C1",
 61        "ARM_C2",
 62        "ARM_C3",
 63    ),
 64)(
 65    eps=_Mv / _Md,
 66    g=constants.g,
 67    p1000=1000 * si.hPa,
 68    Rd=287.0027,
 69    Rv=constants.R / _Mv,
 70    c_pd=1005,
 71    c_pv=1850,
 72    lv=2.5e6,
 73    rho_l=1e3 * si.kg / si.m**3,
 74    T0=constants.zero_Celsius,
 75    ARM_C1=6.1094 * si.hPa,
 76    ARM_C2=17.625 * si.dimensionless,
 77    ARM_C3=243.04 * si.kelvin,
 78)
 79
 80
 81def rho_d(p, qv, theta_std):
 82    return (
 83        p
 84        * (1 - 1 / (1 + const.eps / qv))
 85        / (np.power(p / const.p1000, const.Rd / const.c_pd) * const.Rd * theta_std)
 86    )
 87
 88
 89def drho_dz(g, p, T, qv, lv, dql_dz=0):
 90    Rq = const.Rv / (1 / qv + 1) + const.Rd / (1 + qv)
 91    cp = const.c_pv / (1 / qv + 1) + const.c_pd / (1 + qv)
 92    rho = p / Rq / T
 93    return (g / T * rho * (Rq / cp - 1) - p * lv / cp / T**2 * dql_dz) / Rq
 94
 95
 96# A14 in libcloudph++ 1.0 paper
 97def temperature(rhod, thd):
 98    return thd * np.power(
 99        rhod * thd / const.p1000 * const.Rd,
100        const.Rd / const.c_pd / (1 - const.Rd / const.c_pd),
101    )
102
103
104# A15 in libcloudph++ 1.0 paper
105def pressure(rhod, T, qv):
106    return rhod * (1 + qv) * (const.Rv / (1 / qv + 1) + const.Rd / (1 + qv)) * T
107
108
109def th_dry(th_std, qv):
110    return th_std * np.power(1 + qv / const.eps, const.Rd / const.c_pd)
111
112
113def pvs_Celsius(T):
114    return const.ARM_C1 * np.exp((const.ARM_C2 * T) / (T + const.ARM_C3))
115
116
117def pv(p, qv):
118    return p * qv / (qv + const.eps)
def convert_to(value, unit): 
8def convert_to(value, unit):
9    value /= unit
si = si(kg=1.0, m=1.0, s=1.0, metres=1.0, second=1.0, um=1e-06, hPa=100.0, micrometre=1e-06, minutes=60.0, km=1000.0, dimensionless=1.0, kelvin=1.0, mg=1e-06)
const = const(eps=0.6218501898515706, g=9.80665, p1000=100000.0, Rd=287.0027, Rv=461.52998157091315, c_pd=1005, c_pv=1850, lv=2500000.0, rho_l=1000.0, T0=273.15, ARM_C1=610.9399999999999, ARM_C2=17.625, ARM_C3=243.04)
def rho_d(p, qv, theta_std): 
82def rho_d(p, qv, theta_std):
83    return (
84        p
85        * (1 - 1 / (1 + const.eps / qv))
86        / (np.power(p / const.p1000, const.Rd / const.c_pd) * const.Rd * theta_std)
87    )
def drho_dz(g, p, T, qv, lv, dql_dz=0): 
90def drho_dz(g, p, T, qv, lv, dql_dz=0):
91    Rq = const.Rv / (1 / qv + 1) + const.Rd / (1 + qv)
92    cp = const.c_pv / (1 / qv + 1) + const.c_pd / (1 + qv)
93    rho = p / Rq / T
94    return (g / T * rho * (Rq / cp - 1) - p * lv / cp / T**2 * dql_dz) / Rq
def temperature(rhod, thd): 
 98def temperature(rhod, thd):
 99    return thd * np.power(
100        rhod * thd / const.p1000 * const.Rd,
101        const.Rd / const.c_pd / (1 - const.Rd / const.c_pd),
102    )
def pressure(rhod, T, qv): 
106def pressure(rhod, T, qv):
107    return rhod * (1 + qv) * (const.Rv / (1 / qv + 1) + const.Rd / (1 + qv)) * T
def th_dry(th_std, qv): 
110def th_dry(th_std, qv):
111    return th_std * np.power(1 + qv / const.eps, const.Rd / const.c_pd)
def pvs_Celsius(T): 
114def pvs_Celsius(T):
115    return const.ARM_C1 * np.exp((const.ARM_C2 * T) / (T + const.ARM_C3))
def pv(p, qv): 
118def pv(p, qv):
119    return p * qv / (qv + const.eps)