virtual_temperature() - Code Metrics - Inspection of "Added thermo calculations" - Unidata/MetPy - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Pull Request — master (#337)

unknown

created 2017-03-04 04:55 UTC

virtual_temperature() B

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes	2
Bugs	0	Features	0

Metric	Value
cc	1
c	2
b	0
f	0
dl	0
loc	34
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 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{\textsl{w} + \epsilon}{\epsilon\,(1 + \textsl{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{\textsl{w} + \epsilon}{\epsilon\,(1 + \textsl{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			@exporter.export
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			@exporter.export
237			def parcel_profile(pressure, temperature, dewpt):
238			r"""Calculate the profile a parcel takes through the atmosphere.
239
240			The parcel starts at `temperature`, and `dewpt`, lifted up
241			dry adiabatically to the LCL, and then moist adiabatically from there.
242			`pressure` specifies the pressure levels for the profile.
243
244			Parameters
245			----------
246			pressure : `pint.Quantity`
247			The atmospheric pressure level(s) of interest. The first entry should be the starting
248			point pressure.
249			temperature : `pint.Quantity`
250			The starting temperature
251			dewpt : `pint.Quantity`
252			The starting dew point
253
254			Returns
255			-------
256			`pint.Quantity`
257			The parcel temperatures at the specified pressure levels.
258
259			See Also
260			--------
261			lcl, moist_lapse, dry_lapse
262			"""
263			# Find the LCL
264			l = lcl(pressure[0], temperature, dewpt)[0].to(pressure.units)
265
266			# Find the dry adiabatic profile, including the LCL. We need >= the LCL in case the
267			# LCL is included in the levels. It's slightly redundant in that case, but simplifies
268			# the logic for removing it later.
269			press_lower = concatenate((pressure[pressure >= l], l))
270			t1 = dry_lapse(press_lower, temperature)
271
272			# Find moist pseudo-adiabatic profile starting at the LCL
273			press_upper = concatenate((l, pressure[pressure < l]))
274			t2 = moist_lapse(press_upper, t1[-1]).to(t1.units)
275
276			# Return LCL without the LCL point
277			return concatenate((t1[:-1], t2[1:]))
278
279
280			@exporter.export
281			def vapor_pressure(pressure, mixing):
282			r"""Calculate water vapor (partial) pressure.
283
284			Given total `pressure` and water vapor `mixing` ratio, calculates the
285			partial pressure of water vapor.
286
287			Parameters
288			----------
289			pressure : `pint.Quantity`
290			total atmospheric pressure
291			mixing : `pint.Quantity`
292			dimensionless mass mixing ratio
293
294			Returns
295			-------
296			`pint.Quantity`
297			The ambient water vapor (partial) pressure in the same units as
298			`pressure`.
299
300			Notes
301			-----
302			This function is a straightforward implementation of the equation given in many places,
303			such as [2]_:
304
305			.. math:: e = p \frac{r}{r + \epsilon}
306
307			References
308			----------
309			.. [2] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
310			Survey. 71.
311
312			See Also
313			--------
314			saturation_vapor_pressure, dewpoint
315			"""
316			return pressure * mixing / (epsilon + mixing)
317
318
319			@exporter.export
320			def saturation_vapor_pressure(temperature):
321			r"""Calculate the saturation water vapor (partial) pressure.
322
323			Parameters
324			----------
325			temperature : `pint.Quantity`
326			The temperature
327
328			Returns
329			-------
330			`pint.Quantity`
331			The saturation water vapor (partial) pressure
332
333			See Also
334			--------
335			vapor_pressure, dewpoint
336
337			Notes
338			-----
339			Instead of temperature, dewpoint may be used in order to calculate
340			the actual (ambient) water vapor (partial) pressure.
341
342			The formula used is that from Bolton 1980 [3]_ for T in degrees Celsius:
343
344			.. math:: 6.112 e^\frac{17.67T}{T + 243.5}
345
346			References
347			----------
348			.. [3] Bolton, D., 1980: The Computation of Equivalent Potential
349			Temperature. Mon. Wea. Rev., 108, 1046-1053.
350			"""
351			# Converted from original in terms of C to use kelvin. Using raw absolute values of C in
352			# a formula plays havoc with units support.
353			return sat_pressure_0c * np.exp(17.67 * (temperature - 273.15 * units.kelvin) /
354			(temperature - 29.65 * units.kelvin))
355
356
357			@exporter.export
358			def dewpoint_rh(temperature, rh):
359			r"""Calculate the ambient dewpoint given air temperature and relative humidity.
360
361			Parameters
362			----------
363			temperature : `pint.Quantity`
364			Air temperature
365			rh : `pint.Quantity`
366			Relative humidity expressed as a ratio in the range [0, 1]
367
368			Returns
369			-------
370			`pint.Quantity`
371			The dew point temperature
372
373			See Also
374			--------
375			dewpoint, saturation_vapor_pressure
376			"""
377			return dewpoint(rh * saturation_vapor_pressure(temperature))
378
379
380			@exporter.export
381			def dewpoint(e):
382			r"""Calculate the ambient dewpoint given the vapor pressure.
383
384			Parameters
385			----------
386			e : `pint.Quantity`
387			Water vapor partial pressure
388
389			Returns
390			-------
391			`pint.Quantity`
392			Dew point temperature
393
394			See Also
395			--------
396			dewpoint_rh, saturation_vapor_pressure, vapor_pressure
397
398			Notes
399			-----
400			This function inverts the Bolton 1980 [4]_ formula for saturation vapor
401			pressure to instead calculate the temperature. This yield the following
402			formula for dewpoint in degrees Celsius:
403
404			.. math:: T = \frac{243.5 log(e / 6.112)}{17.67 - log(e / 6.112)}
405
406			References
407			----------
408			.. [4] Bolton, D., 1980: The Computation of Equivalent Potential
409			Temperature. Mon. Wea. Rev., 108, 1046-1053.
410			"""
411			val = np.log(e / sat_pressure_0c)
412			return 0. * units.degC + 243.5 * units.delta_degC * val / (17.67 - val)
413
414
415			@exporter.export
416			def mixing_ratio(part_press, tot_press, molecular_weight_ratio=epsilon):
417			r"""Calculate the mixing ratio of a gas.
418
419			This calculates mixing ratio given its partial pressure and the total pressure of
420			the air. There are no required units for the input arrays, other than that
421			they have the same units.
422
423			Parameters
424			----------
425			part_press : `pint.Quantity`
426			Partial pressure of the constituent gas
427			tot_press : `pint.Quantity`
428			Total air pressure
429			molecular_weight_ratio : `pint.Quantity` or float, optional
430			The ratio of the molecular weight of the constituent gas to that assumed
431			for air. Defaults to the ratio for water vapor to dry air
432			(:math:`\epsilon\approx0.622`).
433
434			Returns
435			-------
436			`pint.Quantity`
437			The (mass) mixing ratio, dimensionless (e.g. Kg/Kg or g/g)
438
439			Notes
440			-----
441			This function is a straightforward implementation of the equation given in many places,
442			such as [5]_:
443
444			.. math:: r = \epsilon \frac{e}{p - e}
445
446			References
447			----------
448			.. [5] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
449			Survey. 73.
450
451			See Also
452			--------
453			saturation_mixing_ratio, vapor_pressure
454			"""
455			return molecular_weight_ratio * part_press / (tot_press - part_press)
456
457
458			@exporter.export
459			def saturation_mixing_ratio(tot_press, temperature):
460			r"""Calculate the saturation mixing ratio of water vapor.
461
462			This calculation is given total pressure and the temperature. The implementation
463			uses the formula outlined in [6]_.
464
465			Parameters
466			----------
467			tot_press: `pint.Quantity`
468			Total atmospheric pressure
469			temperature: `pint.Quantity`
470			The temperature
471
472			Returns
473			-------
474			`pint.Quantity`
475			The saturation mixing ratio, dimensionless
476
477			References
478			----------
479			.. [6] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
480			Survey. 73.
481			"""
482			return mixing_ratio(saturation_vapor_pressure(temperature), tot_press)
483
484
485			@exporter.export
486			def equivalent_potential_temperature(pressure, temperature):
487			r"""Calculate equivalent potential temperature.
488
489			This calculation must be given an air parcel's pressure and temperature.
490			The implementation uses the formula outlined in [7]_.
491
492			Parameters
493			----------
494			pressure: `pint.Quantity`
495			Total atmospheric pressure
496			temperature: `pint.Quantity`
497			The temperature
498
499			Returns
500			-------
501			`pint.Quantity`
502			The corresponding equivalent potential temperature of the parcel
503
504			Notes
505			-----
506			.. math:: \Theta_e = \Theta e^\frac{L_v r_s}{C_{pd} T}
507
508			References
509			----------
510			.. [7] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
511			Survey. 78-79.
512			"""
513			pottemp = potential_temperature(pressure, temperature)
514			smixr = saturation_mixing_ratio(pressure, temperature)
515			return pottemp * np.exp(Lv * smixr / (Cp_d * temperature))
516
517
518			@exporter.export
519			def virtual_temperature(temperature, mixing, molecular_weight_ratio=epsilon):
520			r"""Calculate virtual temperature.
521
522			This calculation must be given an air parcel's temperature and mixing ratio.
523			The implementation uses the formula outlined in [8]_.
524
525			Parameters
526			----------
527			temperature: `pint.Quantity`
528			The temperature
529			mixing : `pint.Quantity`
530			dimensionless mass mixing ratio
531			molecular_weight_ratio : `pint.Quantity` or float, optional
532			The ratio of the molecular weight of the constituent gas to that assumed
533			for air. Defaults to the ratio for water vapor to dry air
534			(:math:`\epsilon\approx0.622`).
535
536			Returns
537			-------
538			`pint.Quantity`
539			The corresponding virtual temperature of the parcel
540
541			Notes
542			-----
543			.. math:: T_v = T \frac{\textsl{w} + \epsilon}{\epsilon\,(1 + \textsl{w})}
544
545			References
546			----------
547			.. [8] Hobbs, Peter V. and Wallace, John M., 2006: Atmospheric Science, an Introductory
548			Survey. 2nd ed. 80.
549			"""
550			return temperature * ((mixing + molecular_weight_ratio) /
551			(molecular_weight_ratio * (1 + mixing)))
552
553
554			@exporter.export
555			def virtual_potential_temperature(pressure, temperature, mixing,
556			molecular_weight_ratio=epsilon):
557			r"""Calculate virtual potential temperature.
558
559			This calculation must be given an air parcel's pressure, temperature, and mixing ratio.
560			The implementation uses the formula outlined in [9]_.
561
562			Parameters
563			----------
564			pressure: `pint.Quantity`
565			Total atmospheric pressure
566			temperature: `pint.Quantity`
567			The temperature
568			mixing : `pint.Quantity`
569			dimensionless mass mixing ratio
570			molecular_weight_ratio : `pint.Quantity` or float, optional
571			The ratio of the molecular weight of the constituent gas to that assumed
572			for air. Defaults to the ratio for water vapor to dry air
573			(:math:`\epsilon\approx0.622`).
574
575			Returns
576			-------
577			`pint.Quantity`
578			The corresponding virtual potential temperature of the parcel
579
580			Notes
581			-----
582			.. math:: \Theta_v = \Theta \frac{\textsl{w} + \epsilon}{\epsilon\,(1 + \textsl{w})}
583
584			References
585			----------
586			.. [9] Markowski, Paul and Richardson, Yvette, 2010: Mesoscale Meteorology in the
587			Midlatitudes. 13.
588			"""
589			pottemp = potential_temperature(pressure, temperature)
590			return virtual_temperature(pottemp, mixing, molecular_weight_ratio)
591
592
593			@exporter.export
594			def density(pressure, temperature, mixing, molecular_weight_ratio=epsilon):
595			r"""Calculate density.
596
597			This calculation must be given an air parcel's pressure, temperature, and mixing ratio.
598			The implementation uses the formula outlined in [10]_.
599
600			Parameters
601			----------
602			temperature: `pint.Quantity`
603			The temperature
604			pressure: `pint.Quantity`
605			Total atmospheric pressure
606			mixing : `pint.Quantity`
607			dimensionless mass mixing ratio
608			molecular_weight_ratio : `pint.Quantity` or float, optional
609			The ratio of the molecular weight of the constituent gas to that assumed
610			for air. Defaults to the ratio for water vapor to dry air
611			(:math:`\epsilon\approx0.622`).
612
613			Returns
614			-------
615			`pint.Quantity`
616			The corresponding density of the parcel
617
618			Notes
619			-----
620			.. math:: \rho = \frac{p}{R_dT_v}
621
622			References
623			----------
624			.. [10] Hobbs, Peter V. and Wallace, John M., 2006: Atmospheric Science, an Introductory
625			Survey. 2nd ed. 67.
626			"""
627			virttemp = virtual_temperature(temperature, mixing, molecular_weight_ratio)
628			return (pressure / (Rd * virttemp)).to(units.kilogram / units.meter ** 3)
629

Unidata / MetPy

Pull Request — master (#337)

virtual_temperature() B

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like