el() - Code Metrics - Inspection of "EL calculation" - Unidata/MetPy - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Pull Request — master (#338)

unknown

created 2017-03-06 22:07 UTC

el() B

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	35
Ratio	100 %

Importance

Changes	1
Bugs	0	Features	0

Metric	Value
cc	2
c	1
b	0
f	0
dl	35
loc	35
rs	8.8571

# Copyright (c) 2008-2015 MetPy Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
"""Contains a collection of thermodynamic calculations."""

from __future__ import division

import numpy as np
import scipy.integrate as si

from .tools import find_intersections
from ..constants import Cp_d, epsilon, kappa, Lv, P0, Rd
from ..package_tools import Exporter
from ..units import atleast_1d, concatenate, units

exporter = Exporter(globals())

sat_pressure_0c = 6.112 * units.millibar


@exporter.export
def potential_temperature(pressure, temperature):
    r"""Calculate the potential temperature.

    Uses the Poisson equation to calculation the potential temperature
    given `pressure` and `temperature`.

    Parameters
    ----------
    pressure : `pint.Quantity`
        The total atmospheric pressure
    temperature : `pint.Quantity`
        The temperature

    Returns
    -------
    `pint.Quantity`
        The potential temperature corresponding to the the temperature and
        pressure.

    See Also
    --------
    dry_lapse

    Notes
    -----
    Formula:

    .. math:: \Theta = T (P_0 / P)^\kappa

    Examples
    --------
    >>> from metpy.units import units
    >>> metpy.calc.potential_temperature(800. * units.mbar, 273. * units.kelvin)
    290.9814150577374
    """
    return temperature * (P0 / pressure).to('dimensionless')**kappa


@exporter.export
def dry_lapse(pressure, temperature):
    r"""Calculate the temperature at a level assuming only dry processes.

    This function lifts a parcel starting at `temperature`, conserving
    potential temperature. The starting pressure should be the first item in
    the `pressure` array.

    Parameters
    ----------
    pressure : `pint.Quantity`
        The atmospheric pressure level(s) of interest
    temperature : `pint.Quantity`
        The starting temperature

    Returns
    -------
    `pint.Quantity`
       The resulting parcel temperature at levels given by `pressure`

    See Also
    --------
    moist_lapse : Calculate parcel temperature assuming liquid saturation
                  processes
    parcel_profile : Calculate complete parcel profile
    potential_temperature
    """
    return temperature * (pressure / pressure[0])**kappa


@exporter.export
def moist_lapse(pressure, temperature):
    r"""Calculate the temperature at a level assuming liquid saturation processes.

    This function lifts a parcel starting at `temperature`. The starting
    pressure should be the first item in the `pressure` array. Essentially,
    this function is calculating moist pseudo-adiabats.

    Parameters
    ----------
    pressure : `pint.Quantity`
        The atmospheric pressure level(s) of interest
    temperature : `pint.Quantity`
        The starting temperature

    Returns
    -------
    `pint.Quantity`
       The temperature corresponding to the the starting temperature and
       pressure levels.

    See Also
    --------
    dry_lapse : Calculate parcel temperature assuming dry adiabatic processes
    parcel_profile : Calculate complete parcel profile

    Notes
    -----
    This function is implemented by integrating the following differential
    equation:

    .. math:: \frac{dT}{dP} = \frac{1}{P} \frac{R_d T + L_v r_s}
                                {C_{pd} + \frac{L_v^2 r_s \epsilon}{R_d T^2}}

    This equation comes from [1]_.

    References
    ----------
    .. [1] Bakhshaii, A. and R. Stull, 2013: Saturated Pseudoadiabats--A
           Noniterative Approximation. J. Appl. Meteor. Clim., 52, 5-15.
    """
    def dt(t, p):
        t = units.Quantity(t, temperature.units)
        p = units.Quantity(p, pressure.units)
        rs = saturation_mixing_ratio(p, t)
        frac = ((Rd * t + Lv * rs) /
                (Cp_d + (Lv * Lv * rs * epsilon / (Rd * t * t)))).to('kelvin')
        return frac / p
    return units.Quantity(si.odeint(dt, atleast_1d(temperature).squeeze(),
                                    pressure.squeeze()).T.squeeze(), temperature.units)


@exporter.export
def lcl(pressure, temperature, dewpt, max_iters=50, eps=1e-2):
    r"""Calculate the lifted condensation level (LCL) using from the starting point.

    The starting state for the parcel is defined by `temperature`, `dewpt`,
    and `pressure`.

    Parameters
    ----------
    pressure : `pint.Quantity`
        The starting atmospheric pressure
    temperature : `pint.Quantity`
        The starting temperature
    dewpt : `pint.Quantity`
        The starting dew point

    Returns
    -------
    `(pint.Quantity, pint.Quantity)`
        The LCL pressure and temperature

    Other Parameters
    ----------------
    max_iters : int, optional
        The maximum number of iterations to use in calculation, defaults to 50.
    eps : float, optional
        The desired absolute error in the calculated value, defaults to 1e-2.

    See Also
    --------
    parcel_profile

    Notes
    -----
    This function is implemented using an iterative approach to solve for the
    LCL. The basic algorithm is:

    1. Find the dew point from the LCL pressure and starting mixing ratio
    2. Find the LCL pressure from the starting temperature and dewpoint
    3. Iterate until convergence

    The function is guaranteed to finish by virtue of the `max_iters` counter.
    """
    w = mixing_ratio(saturation_vapor_pressure(dewpt), pressure)
    new_p = p = pressure
    eps = units.Quantity(eps, p.units)
    while max_iters:
        td = dewpoint(vapor_pressure(p, w))
        new_p = pressure * (td / temperature) ** (1. / kappa)
        if np.abs(new_p - p).max() < eps:
            break
        p = new_p
        max_iters -= 1

    else:
        # We have not converged
        raise RuntimeError('LCL calculation has not converged.')

    return new_p, td


@exporter.export

def lfc(pressure, temperature, dewpt):
    r"""Calculate the level of free convection (LFC).

    This works by finding the first intersection of the ideal parcel path and
    the measured parcel temperature.

    Parameters
    ----------
    pressure : `pint.Quantity`
        The atmospheric pressure
    temperature : `pint.Quantity`
        The temperature at the levels given by `pressure`
    dewpt : `pint.Quantity`
        The dew point at the levels given by `pressure`

    Returns
    -------
    `pint.Quantity`
        The LFC

    See Also
    --------
    parcel_profile
    """
    ideal_profile = parcel_profile(pressure, temperature[0], dewpt[0]).to('degC')

    # The parcel profile and data have the same first data point, so we ignore
    # that point to get the real first intersection for the LFC calculation.
    x, y = find_intersections(pressure[1:], ideal_profile[1:], temperature[1:])
    return x[0], y[0]


@exporter.export

def el(pressure, temperature, dewpt):
    r"""Calculate the equilibrium level (EL).

    This works by finding the last intersection of the ideal parcel path and
    the measured parcel temperature. If there is one or fewer intersections, there is
    no equilibrium level.

    Parameters
    ----------
    pressure : `pint.Quantity`
        The atmospheric pressure
    temperature : `pint.Quantity`
        The temperature at the levels given by `pressure`
    dewpt : `pint.Quantity`
        The dew point at the levels given by `pressure`

    Returns
    -------
    `pint.Quantity, pint.Quantity`
        The EL pressure and temperature

    See Also
    --------
    parcel_profile
    """
    ideal_profile = parcel_profile(pressure, temperature[0], dewpt[0]).to('degC')
    x, y = find_intersections(pressure[1:], ideal_profile[1:], temperature[1:])

    # If there is only one intersection, it's the LFC and we return None.
    if len(x) <= 1:
        return None, None

    else:
        return x[-1], y[-1]


@exporter.export
def parcel_profile(pressure, temperature, dewpt):
    r"""Calculate the profile a parcel takes through the atmosphere.

    The parcel starts at `temperature`, and `dewpt`, lifted up
    dry adiabatically to the LCL, and then moist adiabatically from there.
    `pressure` specifies the pressure levels for the profile.

    Parameters
    ----------
    pressure : `pint.Quantity`
        The atmospheric pressure level(s) of interest. The first entry should be the starting
        point pressure.
    temperature : `pint.Quantity`
        The starting temperature
    dewpt : `pint.Quantity`
        The starting dew point

    Returns
    -------
    `pint.Quantity`
        The parcel temperatures at the specified pressure levels.

    See Also
    --------
    lcl, moist_lapse, dry_lapse
    """
    # Find the LCL
    l = lcl(pressure[0], temperature, dewpt)[0].to(pressure.units)

    # Find the dry adiabatic profile, *including* the LCL. We need >= the LCL in case the
    # LCL is included in the levels. It's slightly redundant in that case, but simplifies
    # the logic for removing it later.
    press_lower = concatenate((pressure[pressure >= l], l))
    t1 = dry_lapse(press_lower, temperature)

    # Find moist pseudo-adiabatic profile starting at the LCL
    press_upper = concatenate((l, pressure[pressure < l]))
    t2 = moist_lapse(press_upper, t1[-1]).to(t1.units)

    # Return LCL *without* the LCL point
    return concatenate((t1[:-1], t2[1:]))


@exporter.export
def vapor_pressure(pressure, mixing):
    r"""Calculate water vapor (partial) pressure.

    Given total `pressure` and water vapor `mixing` ratio, calculates the
    partial pressure of water vapor.

    Parameters
    ----------
    pressure : `pint.Quantity`
        total atmospheric pressure
    mixing : `pint.Quantity`
        dimensionless mass mixing ratio

    Returns
    -------
    `pint.Quantity`
        The ambient water vapor (partial) pressure in the same units as
        `pressure`.

    Notes
    -----
    This function is a straightforward implementation of the equation given in many places,
    such as [2]_:

    .. math:: e = p \frac{r}{r + \epsilon}

    References
    ----------
    .. [2] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
           Survey. 71.

    See Also
    --------
    saturation_vapor_pressure, dewpoint
    """
    return pressure * mixing / (epsilon + mixing)


@exporter.export
def saturation_vapor_pressure(temperature):
    r"""Calculate the saturation water vapor (partial) pressure.

    Parameters
    ----------
    temperature : `pint.Quantity`
        The temperature

    Returns
    -------
    `pint.Quantity`
        The saturation water vapor (partial) pressure

    See Also
    --------
    vapor_pressure, dewpoint

    Notes
    -----
    Instead of temperature, dewpoint may be used in order to calculate
    the actual (ambient) water vapor (partial) pressure.

    The formula used is that from Bolton 1980 [3]_ for T in degrees Celsius:

    .. math:: 6.112 e^\frac{17.67T}{T + 243.5}

    References
    ----------
    .. [3] Bolton, D., 1980: The Computation of Equivalent Potential
           Temperature. Mon. Wea. Rev., 108, 1046-1053.
    """
    # Converted from original in terms of C to use kelvin. Using raw absolute values of C in
    # a formula plays havoc with units support.
    return sat_pressure_0c * np.exp(17.67 * (temperature - 273.15 * units.kelvin) /
                                    (temperature - 29.65 * units.kelvin))


@exporter.export
def dewpoint_rh(temperature, rh):
    r"""Calculate the ambient dewpoint given air temperature and relative humidity.

    Parameters
    ----------
    temperature : `pint.Quantity`
        Air temperature
    rh : `pint.Quantity`
        Relative humidity expressed as a ratio in the range [0, 1]

    Returns
    -------
    `pint.Quantity`
        The dew point temperature

    See Also
    --------
    dewpoint, saturation_vapor_pressure
    """
    return dewpoint(rh * saturation_vapor_pressure(temperature))


@exporter.export
def dewpoint(e):
    r"""Calculate the ambient dewpoint given the vapor pressure.

    Parameters
    ----------
    e : `pint.Quantity`
        Water vapor partial pressure

    Returns
    -------
    `pint.Quantity`
        Dew point temperature

    See Also
    --------
    dewpoint_rh, saturation_vapor_pressure, vapor_pressure

    Notes
    -----
    This function inverts the Bolton 1980 [4]_ formula for saturation vapor
    pressure to instead calculate the temperature. This yield the following
    formula for dewpoint in degrees Celsius:

    .. math:: T = \frac{243.5 log(e / 6.112)}{17.67 - log(e / 6.112)}

    References
    ----------
    .. [4] Bolton, D., 1980: The Computation of Equivalent Potential
           Temperature. Mon. Wea. Rev., 108, 1046-1053.
    """
    val = np.log(e / sat_pressure_0c)
    return 0. * units.degC + 243.5 * units.delta_degC * val / (17.67 - val)


@exporter.export
def mixing_ratio(part_press, tot_press, molecular_weight_ratio=epsilon):
    r"""Calculate the mixing ratio of a gas.

    This calculates mixing ratio given its partial pressure and the total pressure of
    the air. There are no required units for the input arrays, other than that
    they have the same units.

    Parameters
    ----------
    part_press : `pint.Quantity`
        Partial pressure of the constituent gas
    tot_press : `pint.Quantity`
        Total air pressure
    molecular_weight_ratio : `pint.Quantity` or float, optional
        The ratio of the molecular weight of the constituent gas to that assumed
        for air. Defaults to the ratio for water vapor to dry air
        (:math:`\epsilon\approx0.622`).

    Returns
    -------
    `pint.Quantity`
        The (mass) mixing ratio, dimensionless (e.g. Kg/Kg or g/g)

    Notes
    -----
    This function is a straightforward implementation of the equation given in many places,
    such as [5]_:

    .. math:: r = \epsilon \frac{e}{p - e}

    References
    ----------
    .. [5] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
           Survey. 73.

    See Also
    --------
    saturation_mixing_ratio, vapor_pressure
    """
    return molecular_weight_ratio * part_press / (tot_press - part_press)


@exporter.export
def saturation_mixing_ratio(tot_press, temperature):
    r"""Calculate the saturation mixing ratio of water vapor.

    This calculation is given total pressure and the temperature. The implementation
    uses the formula outlined in [6]_.

    Parameters
    ----------
    tot_press: `pint.Quantity`
        Total atmospheric pressure
    temperature: `pint.Quantity`
        The temperature

    Returns
    -------
    `pint.Quantity`
        The saturation mixing ratio, dimensionless

    References
    ----------
    .. [6] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
           Survey. 73.
    """
    return mixing_ratio(saturation_vapor_pressure(temperature), tot_press)


@exporter.export
def equivalent_potential_temperature(pressure, temperature):
    r"""Calculate equivalent potential temperature.

    This calculation must be given an air parcel's pressure and temperature.
    The implementation uses the formula outlined in [7]_.

    Parameters
    ----------
    pressure: `pint.Quantity`
        Total atmospheric pressure
    temperature: `pint.Quantity`
        The temperature

    Returns
    -------
    `pint.Quantity`
        The corresponding equivalent potential temperature of the parcel

    Notes
    -----
    .. math:: \Theta_e = \Theta e^\frac{L_v r_s}{C_{pd} T}

    References
    ----------
    .. [7] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
           Survey. 78-79.
    """
    pottemp = potential_temperature(pressure, temperature)
    smixr = saturation_mixing_ratio(pressure, temperature)
    return pottemp * np.exp(Lv * smixr / (Cp_d * temperature))


@exporter.export
def virtual_temperature(temperature, mixing, molecular_weight_ratio=epsilon):
    r"""Calculate virtual temperature.

    This calculation must be given an air parcel's temperature and mixing ratio.
    The implementation uses the formula outlined in [8]_.

    Parameters
    ----------
    temperature: `pint.Quantity`
        The temperature
    mixing : `pint.Quantity`
        dimensionless mass mixing ratio
    molecular_weight_ratio : `pint.Quantity` or float, optional
        The ratio of the molecular weight of the constituent gas to that assumed
        for air. Defaults to the ratio for water vapor to dry air
        (:math:`\epsilon\approx0.622`).

    Returns
    -------
    `pint.Quantity`
        The corresponding virtual temperature of the parcel

    Notes
    -----
    .. math:: T_v = T \frac{\text{w} + \epsilon}{\epsilon\,(1 + \text{w})}

    References
    ----------
    .. [8] Hobbs, Peter V. and Wallace, John M., 2006: Atmospheric Science, an Introductory
           Survey. 2nd ed. 80.
    """
    return temperature * ((mixing + molecular_weight_ratio) /
                          (molecular_weight_ratio * (1 + mixing)))


@exporter.export
def virtual_potential_temperature(pressure, temperature, mixing,
                                  molecular_weight_ratio=epsilon):
    r"""Calculate virtual potential temperature.

    This calculation must be given an air parcel's pressure, temperature, and mixing ratio.
    The implementation uses the formula outlined in [9]_.

    Parameters
    ----------
    pressure: `pint.Quantity`
        Total atmospheric pressure
    temperature: `pint.Quantity`
        The temperature
    mixing : `pint.Quantity`
        dimensionless mass mixing ratio
    molecular_weight_ratio : `pint.Quantity` or float, optional
        The ratio of the molecular weight of the constituent gas to that assumed
        for air. Defaults to the ratio for water vapor to dry air
        (:math:`\epsilon\approx0.622`).

    Returns
    -------
    `pint.Quantity`
        The corresponding virtual potential temperature of the parcel

    Notes
    -----
    .. math:: \Theta_v = \Theta \frac{\text{w} + \epsilon}{\epsilon\,(1 + \text{w})}

    References
    ----------
    .. [9] Markowski, Paul and Richardson, Yvette, 2010: Mesoscale Meteorology in the
           Midlatitudes. 13.
    """
    pottemp = potential_temperature(pressure, temperature)
    return virtual_temperature(pottemp, mixing, molecular_weight_ratio)


@exporter.export
def density(pressure, temperature, mixing, molecular_weight_ratio=epsilon):
    r"""Calculate density.

    This calculation must be given an air parcel's pressure, temperature, and mixing ratio.
    The implementation uses the formula outlined in [10]_.

    Parameters
    ----------
    temperature: `pint.Quantity`
        The temperature
    pressure: `pint.Quantity`
        Total atmospheric pressure
    mixing : `pint.Quantity`
        dimensionless mass mixing ratio
    molecular_weight_ratio : `pint.Quantity` or float, optional
        The ratio of the molecular weight of the constituent gas to that assumed
        for air. Defaults to the ratio for water vapor to dry air
        (:math:`\epsilon\approx0.622`).

    Returns
    -------
    `pint.Quantity`
        The corresponding density of the parcel

    Notes
    -----
    .. math:: \rho = \frac{p}{R_dT_v}

    References
    ----------
    .. [10] Hobbs, Peter V. and Wallace, John M., 2006: Atmospheric Science, an Introductory
           Survey. 2nd ed. 67.
    """
    virttemp = virtual_temperature(temperature, mixing, molecular_weight_ratio)
    return (pressure / (Rd * virttemp)).to(units.kilogram / units.meter ** 3)


1		# Copyright (c) 2008-2015 MetPy Developers.
2		# Distributed under the terms of the BSD 3-Clause License.
3		# SPDX-License-Identifier: BSD-3-Clause
4		"""Contains a collection of thermodynamic calculations."""
5
6		from __future__ import division
7
8		import numpy as np
9		import scipy.integrate as si
10
11		from .tools import find_intersections
12		from ..constants import Cp_d, epsilon, kappa, Lv, P0, Rd
13		from ..package_tools import Exporter
14		from ..units import atleast_1d, concatenate, units
15
16		exporter = Exporter(globals())
17
18		sat_pressure_0c = 6.112 * units.millibar
19
20
21		@exporter.export
22		def potential_temperature(pressure, temperature):
23		r"""Calculate the potential temperature.
24
25		Uses the Poisson equation to calculation the potential temperature
26		given `pressure` and `temperature`.
27
28		Parameters
29		----------
30		pressure : `pint.Quantity`
31		The total atmospheric pressure
32		temperature : `pint.Quantity`
33		The temperature
34
35		Returns
36		-------
37		`pint.Quantity`
38		The potential temperature corresponding to the the temperature and
39		pressure.
40
41		See Also
42		--------
43		dry_lapse
44
45		Notes
46		-----
47		Formula:
48
49		.. math:: \Theta = T (P_0 / P)^\kappa
50
51		Examples
52		--------
53		>>> from metpy.units import units
54		>>> metpy.calc.potential_temperature(800. * units.mbar, 273. * units.kelvin)
55		290.9814150577374
56		"""
57		return temperature * (P0 / pressure).to('dimensionless')**kappa
58
59
60		@exporter.export
61		def dry_lapse(pressure, temperature):
62		r"""Calculate the temperature at a level assuming only dry processes.
63
64		This function lifts a parcel starting at `temperature`, conserving
65		potential temperature. The starting pressure should be the first item in
66		the `pressure` array.
67
68		Parameters
69		----------
70		pressure : `pint.Quantity`
71		The atmospheric pressure level(s) of interest
72		temperature : `pint.Quantity`
73		The starting temperature
74
75		Returns
76		-------
77		`pint.Quantity`
78		The resulting parcel temperature at levels given by `pressure`
79
80		See Also
81		--------
82		moist_lapse : Calculate parcel temperature assuming liquid saturation
83		processes
84		parcel_profile : Calculate complete parcel profile
85		potential_temperature
86		"""
87		return temperature * (pressure / pressure[0])**kappa
88
89
90		@exporter.export
91		def moist_lapse(pressure, temperature):
92		r"""Calculate the temperature at a level assuming liquid saturation processes.
93
94		This function lifts a parcel starting at `temperature`. The starting
95		pressure should be the first item in the `pressure` array. Essentially,
96		this function is calculating moist pseudo-adiabats.
97
98		Parameters
99		----------
100		pressure : `pint.Quantity`
101		The atmospheric pressure level(s) of interest
102		temperature : `pint.Quantity`
103		The starting temperature
104
105		Returns
106		-------
107		`pint.Quantity`
108		The temperature corresponding to the the starting temperature and
109		pressure levels.
110
111		See Also
112		--------
113		dry_lapse : Calculate parcel temperature assuming dry adiabatic processes
114		parcel_profile : Calculate complete parcel profile
115
116		Notes
117		-----
118		This function is implemented by integrating the following differential
119		equation:
120
121		.. math:: \frac{dT}{dP} = \frac{1}{P} \frac{R_d T + L_v r_s}
122		{C_{pd} + \frac{L_v^2 r_s \epsilon}{R_d T^2}}
123
124		This equation comes from [1]_.
125
126		References
127		----------
128		.. [1] Bakhshaii, A. and R. Stull, 2013: Saturated Pseudoadiabats--A
129		Noniterative Approximation. J. Appl. Meteor. Clim., 52, 5-15.
130		"""
131		def dt(t, p):
132		t = units.Quantity(t, temperature.units)
133		p = units.Quantity(p, pressure.units)
134		rs = saturation_mixing_ratio(p, t)
135		frac = ((Rd * t + Lv * rs) /
136		(Cp_d + (Lv * Lv * rs * epsilon / (Rd * t * t)))).to('kelvin')
137		return frac / p
138		return units.Quantity(si.odeint(dt, atleast_1d(temperature).squeeze(),
139		pressure.squeeze()).T.squeeze(), temperature.units)
140
141
142		@exporter.export
143		def lcl(pressure, temperature, dewpt, max_iters=50, eps=1e-2):
144		r"""Calculate the lifted condensation level (LCL) using from the starting point.
145
146		The starting state for the parcel is defined by `temperature`, `dewpt`,
147		and `pressure`.
148
149		Parameters
150		----------
151		pressure : `pint.Quantity`
152		The starting atmospheric pressure
153		temperature : `pint.Quantity`
154		The starting temperature
155		dewpt : `pint.Quantity`
156		The starting dew point
157
158		Returns
159		-------
160		`(pint.Quantity, pint.Quantity)`
161		The LCL pressure and temperature
162
163		Other Parameters
164		----------------
165		max_iters : int, optional
166		The maximum number of iterations to use in calculation, defaults to 50.
167		eps : float, optional
168		The desired absolute error in the calculated value, defaults to 1e-2.
169
170		See Also
171		--------
172		parcel_profile
173
174		Notes
175		-----
176		This function is implemented using an iterative approach to solve for the
177		LCL. The basic algorithm is:
178
179		1. Find the dew point from the LCL pressure and starting mixing ratio
180		2. Find the LCL pressure from the starting temperature and dewpoint
181		3. Iterate until convergence
182
183		The function is guaranteed to finish by virtue of the `max_iters` counter.
184		"""
185		w = mixing_ratio(saturation_vapor_pressure(dewpt), pressure)
186		new_p = p = pressure
187		eps = units.Quantity(eps, p.units)
188		while max_iters:
189		td = dewpoint(vapor_pressure(p, w))
190		new_p = pressure * (td / temperature) ** (1. / kappa)
191		if np.abs(new_p - p).max() < eps:
192		break
193		p = new_p
194		max_iters -= 1
195
196		else:
197		# We have not converged
198		raise RuntimeError('LCL calculation has not converged.')
199
200		return new_p, td
201
202
203	View Code Duplication	@exporter.export
		0 ignored issues – show Duplication introduced 2017-03-06 20:17 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
204		def lfc(pressure, temperature, dewpt):
205		r"""Calculate the level of free convection (LFC).
206
207		This works by finding the first intersection of the ideal parcel path and
208		the measured parcel temperature.
209
210		Parameters
211		----------
212		pressure : `pint.Quantity`
213		The atmospheric pressure
214		temperature : `pint.Quantity`
215		The temperature at the levels given by `pressure`
216		dewpt : `pint.Quantity`
217		The dew point at the levels given by `pressure`
218
219		Returns
220		-------
221		`pint.Quantity`
222		The LFC
223
224		See Also
225		--------
226		parcel_profile
227		"""
228		ideal_profile = parcel_profile(pressure, temperature[0], dewpt[0]).to('degC')
229
230		# The parcel profile and data have the same first data point, so we ignore
231		# that point to get the real first intersection for the LFC calculation.
232		x, y = find_intersections(pressure[1:], ideal_profile[1:], temperature[1:])
233		return x[0], y[0]
234
235
236	View Code Duplication	@exporter.export
		0 ignored issues – show Duplication introduced 2017-03-06 20:17 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
237		def el(pressure, temperature, dewpt):
238		r"""Calculate the equilibrium level (EL).
239
240		This works by finding the last intersection of the ideal parcel path and
241		the measured parcel temperature. If there is one or fewer intersections, there is
242		no equilibrium level.
243
244		Parameters
245		----------
246		pressure : `pint.Quantity`
247		The atmospheric pressure
248		temperature : `pint.Quantity`
249		The temperature at the levels given by `pressure`
250		dewpt : `pint.Quantity`
251		The dew point at the levels given by `pressure`
252
253		Returns
254		-------
255		`pint.Quantity, pint.Quantity`
256		The EL pressure and temperature
257
258		See Also
259		--------
260		parcel_profile
261		"""
262		ideal_profile = parcel_profile(pressure, temperature[0], dewpt[0]).to('degC')
263		x, y = find_intersections(pressure[1:], ideal_profile[1:], temperature[1:])
264
265		# If there is only one intersection, it's the LFC and we return None.
266		if len(x) <= 1:
267		return None, None
268
269		else:
270		return x[-1], y[-1]
271
272
273		@exporter.export
274		def parcel_profile(pressure, temperature, dewpt):
275		r"""Calculate the profile a parcel takes through the atmosphere.
276
277		The parcel starts at `temperature`, and `dewpt`, lifted up
278		dry adiabatically to the LCL, and then moist adiabatically from there.
279		`pressure` specifies the pressure levels for the profile.
280
281		Parameters
282		----------
283		pressure : `pint.Quantity`
284		The atmospheric pressure level(s) of interest. The first entry should be the starting
285		point pressure.
286		temperature : `pint.Quantity`
287		The starting temperature
288		dewpt : `pint.Quantity`
289		The starting dew point
290
291		Returns
292		-------
293		`pint.Quantity`
294		The parcel temperatures at the specified pressure levels.
295
296		See Also
297		--------
298		lcl, moist_lapse, dry_lapse
299		"""
300		# Find the LCL
301		l = lcl(pressure[0], temperature, dewpt)[0].to(pressure.units)
302
303		# Find the dry adiabatic profile, including the LCL. We need >= the LCL in case the
304		# LCL is included in the levels. It's slightly redundant in that case, but simplifies
305		# the logic for removing it later.
306		press_lower = concatenate((pressure[pressure >= l], l))
307		t1 = dry_lapse(press_lower, temperature)
308
309		# Find moist pseudo-adiabatic profile starting at the LCL
310		press_upper = concatenate((l, pressure[pressure < l]))
311		t2 = moist_lapse(press_upper, t1[-1]).to(t1.units)
312
313		# Return LCL without the LCL point
314		return concatenate((t1[:-1], t2[1:]))
315
316
317		@exporter.export
318		def vapor_pressure(pressure, mixing):
319		r"""Calculate water vapor (partial) pressure.
320
321		Given total `pressure` and water vapor `mixing` ratio, calculates the
322		partial pressure of water vapor.
323
324		Parameters
325		----------
326		pressure : `pint.Quantity`
327		total atmospheric pressure
328		mixing : `pint.Quantity`
329		dimensionless mass mixing ratio
330
331		Returns
332		-------
333		`pint.Quantity`
334		The ambient water vapor (partial) pressure in the same units as
335		`pressure`.
336
337		Notes
338		-----
339		This function is a straightforward implementation of the equation given in many places,
340		such as [2]_:
341
342		.. math:: e = p \frac{r}{r + \epsilon}
343
344		References
345		----------
346		.. [2] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
347		Survey. 71.
348
349		See Also
350		--------
351		saturation_vapor_pressure, dewpoint
352		"""
353		return pressure * mixing / (epsilon + mixing)
354
355
356		@exporter.export
357		def saturation_vapor_pressure(temperature):
358		r"""Calculate the saturation water vapor (partial) pressure.
359
360		Parameters
361		----------
362		temperature : `pint.Quantity`
363		The temperature
364
365		Returns
366		-------
367		`pint.Quantity`
368		The saturation water vapor (partial) pressure
369
370		See Also
371		--------
372		vapor_pressure, dewpoint
373
374		Notes
375		-----
376		Instead of temperature, dewpoint may be used in order to calculate
377		the actual (ambient) water vapor (partial) pressure.
378
379		The formula used is that from Bolton 1980 [3]_ for T in degrees Celsius:
380
381		.. math:: 6.112 e^\frac{17.67T}{T + 243.5}
382
383		References
384		----------
385		.. [3] Bolton, D., 1980: The Computation of Equivalent Potential
386		Temperature. Mon. Wea. Rev., 108, 1046-1053.
387		"""
388		# Converted from original in terms of C to use kelvin. Using raw absolute values of C in
389		# a formula plays havoc with units support.
390		return sat_pressure_0c * np.exp(17.67 * (temperature - 273.15 * units.kelvin) /
391		(temperature - 29.65 * units.kelvin))
392
393
394		@exporter.export
395		def dewpoint_rh(temperature, rh):
396		r"""Calculate the ambient dewpoint given air temperature and relative humidity.
397
398		Parameters
399		----------
400		temperature : `pint.Quantity`
401		Air temperature
402		rh : `pint.Quantity`
403		Relative humidity expressed as a ratio in the range [0, 1]
404
405		Returns
406		-------
407		`pint.Quantity`
408		The dew point temperature
409
410		See Also
411		--------
412		dewpoint, saturation_vapor_pressure
413		"""
414		return dewpoint(rh * saturation_vapor_pressure(temperature))
415
416
417		@exporter.export
418		def dewpoint(e):
419		r"""Calculate the ambient dewpoint given the vapor pressure.
420
421		Parameters
422		----------
423		e : `pint.Quantity`
424		Water vapor partial pressure
425
426		Returns
427		-------
428		`pint.Quantity`
429		Dew point temperature
430
431		See Also
432		--------
433		dewpoint_rh, saturation_vapor_pressure, vapor_pressure
434
435		Notes
436		-----
437		This function inverts the Bolton 1980 [4]_ formula for saturation vapor
438		pressure to instead calculate the temperature. This yield the following
439		formula for dewpoint in degrees Celsius:
440
441		.. math:: T = \frac{243.5 log(e / 6.112)}{17.67 - log(e / 6.112)}
442
443		References
444		----------
445		.. [4] Bolton, D., 1980: The Computation of Equivalent Potential
446		Temperature. Mon. Wea. Rev., 108, 1046-1053.
447		"""
448		val = np.log(e / sat_pressure_0c)
449		return 0. * units.degC + 243.5 * units.delta_degC * val / (17.67 - val)
450
451
452		@exporter.export
453		def mixing_ratio(part_press, tot_press, molecular_weight_ratio=epsilon):
454		r"""Calculate the mixing ratio of a gas.
455
456		This calculates mixing ratio given its partial pressure and the total pressure of
457		the air. There are no required units for the input arrays, other than that
458		they have the same units.
459
460		Parameters
461		----------
462		part_press : `pint.Quantity`
463		Partial pressure of the constituent gas
464		tot_press : `pint.Quantity`
465		Total air pressure
466		molecular_weight_ratio : `pint.Quantity` or float, optional
467		The ratio of the molecular weight of the constituent gas to that assumed
468		for air. Defaults to the ratio for water vapor to dry air
469		(:math:`\epsilon\approx0.622`).
470
471		Returns
472		-------
473		`pint.Quantity`
474		The (mass) mixing ratio, dimensionless (e.g. Kg/Kg or g/g)
475
476		Notes
477		-----
478		This function is a straightforward implementation of the equation given in many places,
479		such as [5]_:
480
481		.. math:: r = \epsilon \frac{e}{p - e}
482
483		References
484		----------
485		.. [5] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
486		Survey. 73.
487
488		See Also
489		--------
490		saturation_mixing_ratio, vapor_pressure
491		"""
492		return molecular_weight_ratio * part_press / (tot_press - part_press)
493
494
495		@exporter.export
496		def saturation_mixing_ratio(tot_press, temperature):
497		r"""Calculate the saturation mixing ratio of water vapor.
498
499		This calculation is given total pressure and the temperature. The implementation
500		uses the formula outlined in [6]_.
501
502		Parameters
503		----------
504		tot_press: `pint.Quantity`
505		Total atmospheric pressure
506		temperature: `pint.Quantity`
507		The temperature
508
509		Returns
510		-------
511		`pint.Quantity`
512		The saturation mixing ratio, dimensionless
513
514		References
515		----------
516		.. [6] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
517		Survey. 73.
518		"""
519		return mixing_ratio(saturation_vapor_pressure(temperature), tot_press)
520
521
522		@exporter.export
523		def equivalent_potential_temperature(pressure, temperature):
524		r"""Calculate equivalent potential temperature.
525
526		This calculation must be given an air parcel's pressure and temperature.
527		The implementation uses the formula outlined in [7]_.
528
529		Parameters
530		----------
531		pressure: `pint.Quantity`
532		Total atmospheric pressure
533		temperature: `pint.Quantity`
534		The temperature
535
536		Returns
537		-------
538		`pint.Quantity`
539		The corresponding equivalent potential temperature of the parcel
540
541		Notes
542		-----
543		.. math:: \Theta_e = \Theta e^\frac{L_v r_s}{C_{pd} T}
544
545		References
546		----------
547		.. [7] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
548		Survey. 78-79.
549		"""
550		pottemp = potential_temperature(pressure, temperature)
551		smixr = saturation_mixing_ratio(pressure, temperature)
552		return pottemp * np.exp(Lv * smixr / (Cp_d * temperature))
553
554
555		@exporter.export
556		def virtual_temperature(temperature, mixing, molecular_weight_ratio=epsilon):
557		r"""Calculate virtual temperature.
558
559		This calculation must be given an air parcel's temperature and mixing ratio.
560		The implementation uses the formula outlined in [8]_.
561
562		Parameters
563		----------
564		temperature: `pint.Quantity`
565		The temperature
566		mixing : `pint.Quantity`
567		dimensionless mass mixing ratio
568		molecular_weight_ratio : `pint.Quantity` or float, optional
569		The ratio of the molecular weight of the constituent gas to that assumed
570		for air. Defaults to the ratio for water vapor to dry air
571		(:math:`\epsilon\approx0.622`).
572
573		Returns
574		-------
575		`pint.Quantity`
576		The corresponding virtual temperature of the parcel
577
578		Notes
579		-----
580		.. math:: T_v = T \frac{\text{w} + \epsilon}{\epsilon\,(1 + \text{w})}
581
582		References
583		----------
584		.. [8] Hobbs, Peter V. and Wallace, John M., 2006: Atmospheric Science, an Introductory
585		Survey. 2nd ed. 80.
586		"""
587		return temperature * ((mixing + molecular_weight_ratio) /
588		(molecular_weight_ratio * (1 + mixing)))
589
590
591		@exporter.export
592		def virtual_potential_temperature(pressure, temperature, mixing,
593		molecular_weight_ratio=epsilon):
594		r"""Calculate virtual potential temperature.
595
596		This calculation must be given an air parcel's pressure, temperature, and mixing ratio.
597		The implementation uses the formula outlined in [9]_.
598
599		Parameters
600		----------
601		pressure: `pint.Quantity`
602		Total atmospheric pressure
603		temperature: `pint.Quantity`
604		The temperature
605		mixing : `pint.Quantity`
606		dimensionless mass mixing ratio
607		molecular_weight_ratio : `pint.Quantity` or float, optional
608		The ratio of the molecular weight of the constituent gas to that assumed
609		for air. Defaults to the ratio for water vapor to dry air
610		(:math:`\epsilon\approx0.622`).
611
612		Returns
613		-------
614		`pint.Quantity`
615		The corresponding virtual potential temperature of the parcel
616
617		Notes
618		-----
619		.. math:: \Theta_v = \Theta \frac{\text{w} + \epsilon}{\epsilon\,(1 + \text{w})}
620
621		References
622		----------
623		.. [9] Markowski, Paul and Richardson, Yvette, 2010: Mesoscale Meteorology in the
624		Midlatitudes. 13.
625		"""
626		pottemp = potential_temperature(pressure, temperature)
627		return virtual_temperature(pottemp, mixing, molecular_weight_ratio)
628
629
630		@exporter.export
631		def density(pressure, temperature, mixing, molecular_weight_ratio=epsilon):
632		r"""Calculate density.
633
634		This calculation must be given an air parcel's pressure, temperature, and mixing ratio.
635		The implementation uses the formula outlined in [10]_.
636
637		Parameters
638		----------
639		temperature: `pint.Quantity`
640		The temperature
641		pressure: `pint.Quantity`
642		Total atmospheric pressure
643		mixing : `pint.Quantity`
644		dimensionless mass mixing ratio
645		molecular_weight_ratio : `pint.Quantity` or float, optional
646		The ratio of the molecular weight of the constituent gas to that assumed
647		for air. Defaults to the ratio for water vapor to dry air
648		(:math:`\epsilon\approx0.622`).
649
650		Returns
651		-------
652		`pint.Quantity`
653		The corresponding density of the parcel
654
655		Notes
656		-----
657		.. math:: \rho = \frac{p}{R_dT_v}
658
659		References
660		----------
661		.. [10] Hobbs, Peter V. and Wallace, John M., 2006: Atmospheric Science, an Introductory
662		Survey. 2nd ed. 67.
663		"""
664		virttemp = virtual_temperature(temperature, mixing, molecular_weight_ratio)
665		return (pressure / (Rd * virttemp)).to(units.kilogram / units.meter ** 3)
666

Unidata / MetPy

Pull Request — master (#338)

el() B

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like