diff --git a/localbuild/meta.yaml b/localbuild/meta.yaml index 00ed5bb56..698519e1d 100644 --- a/localbuild/meta.yaml +++ b/localbuild/meta.yaml @@ -80,6 +80,7 @@ requirements: - sqlalchemy <1.4.0 - sqlite <3.35.1 - gpxpy >=1.4.2 + - metpy test: imports: diff --git a/mslib/thermolib.py b/mslib/thermolib.py index dd4636913..d9a094212 100644 --- a/mslib/thermolib.py +++ b/mslib/thermolib.py @@ -26,12 +26,11 @@ limitations under the License. """ -# The function sat_vapour_pressure() has been ported from the IDL function -# 'VaporPressure' by Holger Voemel, available at http://cires.colorado.edu/~voemel/vp.html. - import numpy import scipy.integrate import logging +import metpy.calc as mpcalc +from metpy.units import units class VapourPressureError(Exception): @@ -43,353 +42,21 @@ def __init__(self, error_string): logging.debug("%s", error_string) -def sat_vapour_pressure(t, liquid='HylandWexler', ice='GoffGratch', - force_phase='None'): - """ - Compute the saturation vapour pressure over liquid water and over ice - with a variety of formulations. - - This function is a direct port from the IDL function 'VaporPressure' by - Holger Voemel, available at http://cires.colorado.edu/~voemel/vp.html. - - By default, for temperatures > 0 degC, the saturation pressure over - liquid water is computed; from temperatures <= 0 degC a formulation - over ice is used. - - The current default fomulas are Hyland and Wexler for liquid and - Goff Gratch for ice. (hv20040521) +def sat_vapour_pressure(t): + """Compute saturation vapour presure in Pa from temperature. Arguments: - t -- Temperature in [K]. Can be a scaler or an n-dimensional NumPy array. - liquid -- Optional; specify the formulation for computing the saturation - pressure over liquid water. Can be one of: - - HylandWexler, GoffGratch, Wexler, MagnusTeten, Buck_original, - Buck_manual, WMO_Goff, WMO2000, Sonntag, Bolton, [Fukuta (N/A)], - IAPWS, MurphyKoop. - - ice -- Optional; specify the formulation for computing the saturation - pressure over ice. Can be one of: + t -- temperature in [K] - MartiMauersberger, HylandWexler, GoffGratch, MagnusTeten, - Buck_original, Buck_manual, WMO_Goff, Sonntag, MurphyKoop. + Returns: Saturation Vapour Pressure in [Pa], in the same dimensions as the input. + """ + v_pr = mpcalc.saturation_vapor_pressure(t * units.kelvin) - Please have a look at the source code for further information - about the formulations. + # Convert return value units from mbar to Pa. + return v_pr.to('Pa').magnitude - force_phase -- Optional; force liquid or ice phase to avoid automatic - switching of formulations at 0 degC. Can be 'liquid' - or 'ice'. - Returns: - Saturation vapor pressure [Pa], in the same dimensions as the input. - """ - - # Make sure the input is a NumPy array. - if numpy.isscalar(t): - t = numpy.array([t]) - input_scalar = True - else: - t = numpy.array(t) - input_scalar = False - - # Get indexes of input temperatures above and below freezing, to select - # the appropriate method for each temperature. - if force_phase == "ice": - idx_ice = () # numpy.where(t is not None) - idx_liq = None - elif force_phase == "liquid": - idx_liq = () # numpy.where(t is not None) - idx_ice = None - elif force_phase == "None": - idx_ice = numpy.where(t <= 273.15) - idx_liq = numpy.where(t > 273.15) - else: - raise VapourPressureError("Cannot recognize the force_phase " - f"keyword: '{force_phase}' (valid are ice, liquid, None)") - - # Initialise output field. - e_sat = numpy.zeros(numpy.shape(t)) - - # ============================================================================= - # Calculate saturation pressure over liquid water ---------------------------- - if not force_phase == 'ice': - - if liquid == 'MartiMauersberger': - raise VapourPressureError("Marti and Mauersberger don't " - "have a vapour pressure curve over liquid.") - - elif liquid == 'HylandWexler': - # Source: Hyland, R. W. and A. Wexler, Formulations for the - # Thermodynamic Properties of the saturated Phases of H2O - # from 173.15K to 473.15K, ASHRAE Trans, 89(2A), 500-519, 1983. - e_sat[idx_liq] = (numpy.exp((-0.58002206E4 / t[idx_liq]) + - 0.13914993E1 - - 0.48640239E-1 * t[idx_liq] + - 0.41764768E-4 * t[idx_liq] ** 2. - - 0.14452093E-7 * t[idx_liq] ** 3. + - 0.65459673E1 * numpy.log(t[idx_liq])) / 100.) - - elif liquid == 'Wexler': - # Wexler, A., Vapor pressure formulation for ice, Journal of - # Research of the National Bureau of Standards-A. 81A, 5-20, 1977. - e_sat[idx_liq] = (numpy.exp(-2.9912729E3 * t[idx_liq] ** (-2.) - - 6.0170128E3 * t[idx_liq] ** (-1.) + - 1.887643854E1 * t[idx_liq] ** 0. - - 2.8354721E-2 * t[idx_liq] ** 1. + - 1.7838301E-5 * t[idx_liq] ** 2. - - 8.4150417E-10 * t[idx_liq] ** 3. - - 4.4412543E-13 * t[idx_liq] ** 4. + - 2.858487 * numpy.log(t[idx_liq])) / 100.) - - elif liquid == 'GoffGratch': - # Goff Gratch formulation. - # Source: Smithsonian Meteorological Tables, 5th edition, - # p. 350, 1984 - # From original source: Goff and Gratch (1946), p. 107. - ts = 373.16 # steam point temperature in K - ews = 1013.246 # saturation pressure at steam point - # temperature, normal atmosphere - e_sat[idx_liq] = 10. ** (-7.90298 * ((ts / t[idx_liq]) - 1.) + - 5.02808 * numpy.log10((ts / t[idx_liq])) - - 1.3816E-7 * (10. ** (11.344 * (1. - (t[idx_liq] / ts))) - 1.) + - 8.1328E-3 * (10. ** (-3.49149 * ((ts / t[idx_liq]) - 1)) - 1.) + - numpy.log10(ews)) - - elif liquid == 'MagnusTeten': - # Source: Murray, F. W., On the computation of saturation - # vapor pressure, J. Appl. Meteorol., 6, 203-204, 1967. - tc = t - 273.15 - e_sat[idx_liq] = 10. ** (7.5 * (tc[idx_liq]) / (tc[idx_liq] + 237.5) + 0.7858) - - elif liquid == 'Buck_original': - # Bucks vapor pressure formulation based on Tetens formula - # Source: Buck, A. L., New equations for computing vapor - # pressure and enhancement factor, J. Appl. Meteorol., 20, - # 1527-1532, 1981. - tc = t - 273.15 - e_sat[idx_liq] = 6.1121 * numpy.exp(17.502 * tc[idx_liq] / (240.97 + tc[idx_liq])) - - elif liquid == 'Buck_manual': - # Bucks vapor pressure formulation based on Tetens formula - # Source: Buck Research, Model CR-1A Hygrometer Operating - # Manual, Sep 2001 - tc = t - 273.15 - e_sat[idx_liq] = 6.1121 * numpy.exp((18.678 - (tc[idx_liq] / 234.5)) * - (tc[idx_liq]) / (257.14 + tc[idx_liq])) - - elif liquid == 'WMO_Goff': - # Intended WMO formulation, originally published by Goff (1957) - # incorrectly referenced by WMO technical regulations, WMO-NO 49, - # Vol I, General Meteorological Standards and Recommended - # Practices, App. A, Corrigendum Aug 2000. - # and incorrectly referenced by WMO technical regulations, - # WMO-NO 49, Vol I, General Meteorological Standards and - # Recommended Practices, App. A, 1988. - ts = 273.16 # steam point temperature in K - e_sat[idx_liq] = 10. ** (10.79574 * (1. - (ts / t[idx_liq])) - - 5.02800 * numpy.log10((t[idx_liq] / ts)) + - 1.50475E-4 * (1. - 10. ** (-8.2969 * ((t[idx_liq] / ts) - 1.))) + - 0.42873E-3 * (10. ** (+4.76955 * (1. - (ts / t[idx_liq]))) - 1.) + - 0.78614) - - elif liquid == 'WMO2000': - # WMO formulation, which is very similar to Goff Gratch - # Source: WMO technical regulations, WMO-NO 49, Vol I, - # General Meteorological Standards and Recommended Practices, - # App. A, Corrigendum Aug 2000. - ts = 273.16 # steam point temperature in K - e_sat[idx_liq] = 10. ** (10.79574 * (1. - (ts / t[idx_liq])) - - 5.02800 * numpy.log10((t[idx_liq] / ts)) + - 1.50475E-4 * (1. - 10. ** (-8.2969 * ((t[idx_liq] / ts) - 1.))) + - 0.42873E-3 * (10. ** (-4.76955 * (1. - (ts / t[idx_liq]))) - 1.) + - 0.78614) - - elif liquid == 'Sonntag': - # Source: Sonntag, D., Advancements in the field of hygrometry, - # Meteorol. Z., N. F., 3, 51-66, 1994. - e_sat[idx_liq] = numpy.exp(-6096.9385 * t[idx_liq] ** (-1.) + - 16.635794 - - 2.711193E-2 * t[idx_liq] ** 1. + - 1.673952E-5 * t[idx_liq] ** 2. + - 2.433502 * numpy.log(t[idx_liq])) - - elif liquid == 'Bolton': - # Source: Bolton, D., The computation of equivalent potential - # temperature, Monthly Weather Report, 108, 1046-1053, 1980. - # equation (10) - tc = t - 273.15 - e_sat[idx_liq] = 6.112 * numpy.exp(17.67 * tc[idx_liq] / (tc[idx_liq] + 243.5)) - - # THIS CURVE LOOKS WRONG! - # elif liquid == 'Fukuta': - # # Source: Fukuta, N. and C. M. Gramada, Vapor pressure - # # measurement of supercooled water, J. Atmos. Sci., 60, - # # 1871-1875, 2003. - # # This paper does not give a vapor pressure formulation, - # # but rather a correction over the Smithsonian Tables. - # # Thus calculate the table value first, then use the - # # correction to get to the measured value. - # ts = 373.16 # steam point temperature in K - # ews = 1013.246 # saturation pressure at steam point - # # temperature, normal atmosphere - - # e_sat[idx_liq] = 10.**(-7.90298*(ts/t[idx_liq]-1.) - # + 5.02808 * numpy.log10(ts/t[idx_liq]) - # - 1.3816E-7 * (10.**(11.344*(1.-t[idx_liq]/ts))-1.) - # + 8.1328E-3*(10.**(-3.49149*(ts/t[idx_liq]-1)) -1.) - # + numpy.log10(ews)) - - # tc = t - 273.15 - # x = tc[idx_liq] + 19 - # e_sat[idx_liq] = e_sat[idx_liq] * (0.9992 + 7.113E-4*x - # - 1.847E-4*x**2. - # + 1.189E-5*x**3. - # + 1.130E-7*x**4. - # - 1.743E-8*x**5.) - - # e_sat[numpy.where(tc < -39.)] = None - - elif liquid == 'IAPWS': - # Source: Wagner W. and A. Pruss (2002), The IAPWS - # formulation 1995 for the thermodynamic properties - # of ordinary water substance for general and scientific - # use, J. Phys. Chem. Ref. Data, 31(2), 387-535. - # This is the 'official' formulation from the International - # Association for the Properties of Water and Steam - # The valid range of this formulation is 273.16 <= T <= - # 647.096 K and is based on the ITS90 temperature scale. - Tc = 647.096 # K : Temperature at the critical point - Pc = 22.064 * 10 ** 4 # hPa : Vapor pressure at the critical point - nu = (1. - (t[idx_liq] / Tc)) - a1 = -7.85951783 - a2 = 1.84408259 - a3 = -11.7866497 - a4 = 22.6807411 - a5 = -15.9618719 - a6 = 1.80122502 - e_sat[idx_liq] = Pc * numpy.exp(Tc / t[idx_liq] * - (a1 * nu + a2 * nu ** 1.5 + a3 * nu ** 3. + - a4 * nu ** 3.5 + a5 * nu ** 4. + a6 * nu ** 7.5)) - - elif liquid == 'MurphyKoop': - # Source : Murphy and Koop, Review of the vapour pressure - # of ice and supercooled water for atmospheric applications, - # Q. J. R. Meteorol. Soc (2005), 131, pp. 1539-1565. - e_sat[idx_liq] = (numpy.exp(54.842763 - (6763.22 / t[idx_liq]) - - 4.210 * numpy.log(t[idx_liq]) + - 0.000367 * t[idx_liq] + - numpy.tanh(0.0415 * (t[idx_liq] - 218.8)) * - (53.878 - (1331.22 / t[idx_liq]) - - 9.44523 * numpy.log(t[idx_liq]) + - 0.014025 * t[idx_liq])) / 100.) - - else: - raise VapourPressureError("Unkown method for computing " - f"the vapour pressure curve over liquid: {liquid}") - - # ============================================================================= - # Calculate saturation pressure over ice ------------------------------------- - if not force_phase == 'liquid': - - if ice == 'WMO2000': - ice = 'WMO_Goff' - - if ice == 'IAWPS': - raise VapourPressureError("IAPWS does not provide a vapour " - "pressure formulation over ice") - - elif ice == 'MartiMauersberger': - # Source: Marti, J. and K Mauersberger, A survey and new - # measurements of ice vapor pressure at temperatures between - # 170 and 250 K, GRL 20, 363-366, 1993. - e_sat[idx_ice] = (10. ** ((-2663.5 / t[idx_ice]) + 12.537) / 100.) - - elif ice == 'HylandWexler': - # Source Hyland, R. W. and A. Wexler, Formulations for the - # Thermodynamic Properties of the saturated Phases of H2O - # from 173.15K to 473.15K, ASHRAE Trans, 89(2A), 500-519, 1983. - e_sat[idx_ice] = (numpy.exp((-0.56745359E4 / t[idx_ice]) + - 0.63925247E1 - - 0.96778430E-2 * t[idx_ice] + - 0.62215701E-6 * t[idx_ice] ** 2. + - 0.20747825E-8 * t[idx_ice] ** 3. - - 0.94840240E-12 * t[idx_ice] ** 4. + - 0.41635019E1 * numpy.log(t[idx_ice])) / 100.) - - elif ice == 'GoffGratch': - # Source: Smithsonian Meteorological Tables, 5th edition, - # p. 350, 1984 - - ei0 = 6.1071 # mbar - T0 = 273.16 # freezing point in K - - e_sat[idx_ice] = 10. ** (-9.09718 * ((T0 / t[idx_ice]) - 1.) - - 3.56654 * numpy.log10((T0 / t[idx_ice])) + - 0.876793 * (1. - (t[idx_ice] / T0)) + - numpy.log10(ei0)) - - elif ice == 'MagnusTeten': - # Source: Murray, F. W., On the computation of saturation - # vapour pressure, J. Appl. Meteorol., 6, 203-204, 1967. - tc = t - 273.15 - e_sat[idx_ice] = 10. ** (9.5 * tc[idx_ice] / (265.5 + tc[idx_ice]) + 0.7858) - - elif ice == 'Buck_original': - # Bucks vapor pressure formulation based on Tetens formula - # Source: Buck, A. L., New equations for computing vapor - # pressure and enhancement factor, J. Appl. Meteorol., 20, - # 1527-1532, 1981. - tc = t - 273.15 - e_sat[idx_ice] = 6.1115 * numpy.exp(22.452 * tc[idx_ice] / (272.55 + tc[idx_ice])) - - elif ice == 'Buck_manual': - # Bucks vapor pressure formulation based on Tetens formula - # Source: Buck Research, Model CR-1A Hygrometer Operating - # Manual, Sep 2001 - tc = t - 273.15 - e_sat[idx_ice] = 6.1115 * numpy.exp((23.036 - (tc[idx_ice] / 333.7)) * - tc[idx_ice] / (279.82 + tc[idx_ice])) - - elif ice == 'WMO_Goff': - # WMO formulation, which is very similar to Goff Gratch - # Source: WMO technical regulations, WMO-NO 49, Vol I, - # General Meteorological Standards and Recommended Practices, - # Aug 2000, App. A. - - T0 = 273.16 # steam point temperature in K - - e_sat[idx_ice] = 10. ** (-9.09685 * ((T0 / t[idx_ice]) - 1.) - - 3.56654 * numpy.log10((T0 / t[idx_ice])) + - 0.87682 * (1. - (t[idx_ice] / T0)) + 0.78614) - - elif ice == 'Sonntag': - # Source: Sonntag, D., Advancements in the field of hygrometry, - # Meteorol. Z., N. F., 3, 51-66, 1994. - e_sat[idx_ice] = numpy.exp(-6024.5282 * t[idx_ice] ** (-1.) + - 24.721994 + - 1.0613868E-2 * t[idx_ice] ** 1. - - 1.3198825E-5 * t[idx_ice] ** 2. - - 0.49382577 * numpy.log(t[idx_ice])) - - elif ice == 'MurphyKoop': - # Source: Murphy and Koop, Review of the vapour pressure of ice - # and supercooled water for atmospheric applications, Q. J. R. - # Meteorol. Soc (2005), 131, pp. 1539-1565. - e_sat[idx_ice] = (numpy.exp(9.550426 - (5723.265 / t[idx_ice]) + - 3.53068 * numpy.log(t[idx_ice]) - - 0.00728332 * t[idx_ice]) / 100.) - - else: - raise VapourPressureError("Unkown method for computing " - f"the vapour pressure curve over ice: {ice}") - - # Convert return value units from hPa to Pa. - return e_sat * 100. if not input_scalar else e_sat[0] * 100. - - -def rel_hum(p, t, q, liquid='HylandWexler', ice='GoffGratch', - force_phase='None'): +def rel_hum(p, t, q): """Compute relative humidity in [%] from pressure, temperature, and specific humidity. @@ -398,39 +65,17 @@ def rel_hum(p, t, q, liquid='HylandWexler', ice='GoffGratch', t -- temperature in [K] q -- specific humidity in [kg/kg] - p, t and q can be scalars of NumPy arrays. They just have to either all - scalars, or all arrays. - - liquid, ice, force_phase -- optional keywords to control the calculation - of the saturation vapour pressure; see - help of function 'sat_vapour_pressure()' for - details. - Returns: Relative humidity in [%]. Same dimension as input fields. """ - if not (numpy.isscalar(p) or numpy.isscalar(t) or numpy.isscalar(q)): - if not isinstance(p, numpy.ndarray): - p = numpy.array(p) - if not isinstance(t, numpy.ndarray): - t = numpy.array(t) - if not isinstance(q, numpy.ndarray): - q = numpy.array(q) - - # Compute mixing ratio w from specific humidiy q. - w = q / (1. - q) - - # Compute saturation vapour pressure from temperature t. - e_sat = sat_vapour_pressure(t, liquid=liquid, ice=ice, - force_phase=force_phase) + p = units.Quantity(p, "Pa") + t = units.Quantity(t, "K") + rel_humidity = mpcalc.relative_humidity_from_specific_humidity(p, t, q) - # Compute saturation mixing ratio from e_sat and pressure p. - w_sat = 0.622 * e_sat / (p - e_sat) + # Return specific humidity in [%]. + return rel_humidity * 100 - # Return the relative humidity, computed from w and w_sat. - return 100. * w / w_sat - -def virt_temp(t, q, method='exact'): +def virt_temp(t, q): """ Compute virtual temperature in [K] from temperature and specific humidity. @@ -442,33 +87,12 @@ def virt_temp(t, q, method='exact'): t and q can be scalars of NumPy arrays. They just have to either all scalars, or all arrays. - method -- optional keyword to specify the equation used. Default is - 'exact', which uses - Tv = T * (q + 0.622(1-q)) / 0.622, - 'approx' uses - Tv = T * (1 + 0.61w), - with w = q/(1-q) being the water vapour mixing ratio. - - Reference: Wallace&Hobbs 2nd ed., eq. 3.16, 3.59, and 3.60 - (substitute w=q/(1-q) in 3.16 and 3.59 to obtain the exact - formula). - Returns: Virtual temperature in [K]. Same dimension as input fields. """ - if not (numpy.isscalar(t) or numpy.isscalar(q)): - if not isinstance(t, numpy.ndarray): - t = numpy.array(t) - if not isinstance(q, numpy.ndarray): - q = numpy.array(q) - - if method == 'exact': - return t * (q + 0.622 * (1. - q)) / 0.622 - elif method == 'approx': - # Compute mixing ratio w from specific humidiy q. - w = q / (1. - q) - return t * (1. + 0.61 * w) - else: - raise TypeError('virtual temperature method not understood') + t = units.Quantity(t, "K") + mix_rat = mpcalc.mixing_ratio_from_specific_humidity(q) + v_temp = mpcalc.virtual_temperature(t, mix_rat) + return v_temp def geop_difference(p, t, method='trapz', axis=-1): @@ -521,129 +145,6 @@ def geop_difference(p, t, method='trapz', axis=-1): raise TypeError('integration method for geopotential not understood') -def geop_thickness(p, t, q=None, cumulative=False, axis=-1): - """ - Compute the geopotential thickness in [m] between the pressure levels - given by the first and last element in p (= pressure). - - Implements the hypsometric equation (1.18) from Holton, 3rd edition (or - alternatively (3.24) in Wallace and Hobbs, 2nd ed.). - - Arguments: - p -- pressure in [Pa] - t -- temperature in [K] - q -- [optional] specific humidity in [kg/kg]. If q is given, T will - be converted to virtual temperature to account for the effects - of moisture in the air. - - All inputs need to be NumPy arrays with at least 2 elements. - - cumulative -- optional keyword to specify whether the single geopotential - thickness between p[0] and p[-1] is returned (False, default), - or whether an array containing the thicknesses between - p[0] and all other elements in p is returned (True). The - latter option is useful for computing the geopotential height - of all model levels. - - axis -- see geop_difference(). - - Uses geop_difference() for the integral in the above equations. - - Returns: Geopotential thickness between p[0] and p[-1] in [m]. - If 'cumtrapz' is specified, an array of dimension dim(p)-1 - will be returned, in which value n represents the geopotential - thickness between p[0] and p[n+1]. - """ - - # Check whether humidity effects should be considered. If q is specified, - # simply evaluate the hypsometric equation with virtual temperature instead - # of absolute temperature (see Wallace and Hobbs, 2nd ed., section 3.2.1). - if q is None: - tv = t - else: - tv = virt_temp(t, q) - - # Evaluate equation 3.24 in Wallace and Hobbs: - # delta Z = -Rd/g0 * int( Tv, d ln(p), p1, p2 ), - # where Z denotes the geopotential height, Z = phi/g0. - return -1. / 9.80665 * geop_difference(p, tv, method='cumtrapz' if cumulative else 'trapz', axis=axis) - - -def spec_hum_from_pTd(p, td, liquid='HylandWexler'): - """ - Computes specific humidity in [kg/kg] from pressure and dew point - temperature. - - Arguments: - p -- pressure in [Pa] - td -- dew point temperature in [K] - - p and td can be scalars or NumPy arrays. They just have to either both - scalars, or both arrays. - - liquid -- optional keyword to specify the method used for computing the - saturation water wapour. See sat_vapour_pressure() for - further details. - - Returns: specific humidity in [kg/kg]. Same dimensions as the inputs. - - Method: - Specific humidity q = w / (1+w), with w = mixing ratio. (Wallace & Hobbs, - 2nd ed., (3.57)). W&H write: 'The dew point [Td] is the temperature at - which the saturation mixing ratio ws with respect to liquid water becomes - equal to the actual mixing ratio w.'. Hence we need ws(Td). - From W&H 3.62, we get ws = 0.622 * es / (p-es). Plugging this into the - above equation for q and simplifying, we get - q = 0.622 * es / (p + es * [0.622-1.]) - """ - # Compute saturation vapour pressure from dew point temperature td. - e_sat = sat_vapour_pressure(td, liquid=liquid) - - return 0.622 * e_sat / (p + e_sat * (0.622 - 1.)) - - -def dewpoint_approx(p, q, method='Bolton'): - """ - Computes dew point in [K] from pressure and specific humidity. - - Arguments: - p -- pressure in [Pa] - q -- specific humidity in [kg/kg] - - p and q can be scalars or NumPy arrays. They just have to either both - scalars, or both arrays. - - method -- optional keyword to specify the method used to approximate - the dew point temperature. Valid values are: - - 'Bolton' (default): Use the inversion of Bolton (1980), eq. - 10, to compute dewpoint. According to Bolton, this is accurate - to 0.03 K in the range 238..308 K. See also Emanuel (1994, - 'Atmospheric Convection', eq. 4.6.2). - - Returns: dew point temperature in [K]. - """ - if not (numpy.isscalar(p) or numpy.isscalar(q)): - if not isinstance(p, numpy.ndarray): - p = numpy.array(p) - if not isinstance(q, numpy.ndarray): - q = numpy.array(q) - - # Compute mixing ratio w from specific humidiy q. - w = q / (1. - q) - - # Compute vapour pressure from pressure and mixing ratio - # (Wallace and Hobbs 2nd ed. eq. 3.59). - e_q = w / (w + 0.622) * p - - if method == 'Bolton': - td = (243.5 / ((17.67 / numpy.log(e_q / 100. / 6.112)) - 1)) + 273.15 - else: - raise ValueError(f"invalid dew point method '{method}'") - - return td - - def pot_temp(p, t): """ Computes potential temperature in [K] from pressure and temperature. @@ -656,16 +157,14 @@ def pot_temp(p, t): scalars, or both arrays. Returns: potential temperature in [K]. Same dimensions as the inputs. - - Method: - theta = T * (p0/p)^(R/cp) - with p0 = 100000. Pa, R = 287.058 JK-1kg-1, cp = 1004 JK-1kg-1. """ - return t * (100000. / p) ** (287.058 / 1004.) + p = units.Quantity(p, "Pa") + t = units.Quantity(t, "K") + potential_temp = mpcalc.potential_temperature(p, t) + return potential_temp -def eqpt_approx(p, t, q, liquid='HylandWexler', ice='GoffGratch', - force_phase='None'): +def eqpt_approx(p, t, q): """ Computes equivalent potential temperature in [K] from pressure, temperature and specific humidity. @@ -679,31 +178,12 @@ def eqpt_approx(p, t, q, liquid='HylandWexler', ice='GoffGratch', Returns: equivalent potential temperature in [K]. Same dimensions as the inputs. - - Method: - theta_e = theta * exp((Lv*w_sat)/(cp*T)) - with theta = potential temperature (see pot_temp()), Lv = 2.25e6 Jkg-1, - cp = 1004 JK-1kg-1. - - Reference: Wallace & Hobbs, 2nd ed., eq. 3.71 """ - # Compute potential temperature from p and t. - theta = pot_temp(p, t) - - # Compute saturation vapour pressure from temperature t. - e_sat = sat_vapour_pressure(t, liquid=liquid, ice=ice, - force_phase=force_phase) - - # Compute saturation mixing ratio from e_sat and pressure p. - w_sat = 0.622 * e_sat / (p - e_sat) - - # Latent heat of evaporation. - Lv = 2.25 * 1.e6 - cp = 1004. - - # Equation 3.71 from Wallace & Hobbs, 2nd ed. - theta_e = theta * numpy.exp((Lv * w_sat) / (cp * t)) - return theta_e + p = units.Quantity(p, "Pa") + t = units.Quantity(t, "K") + dew_temp = mpcalc.dewpoint_from_specific_humidity(p, t, q) + eqpt_temp = mpcalc.equivalent_potential_temperature(p, t, dew_temp) + return eqpt_temp.to('degC').magnitude def omega_to_w(omega, p, t): @@ -718,18 +198,12 @@ def omega_to_w(omega, p, t): All inputs can be scalars or NumPy arrays. Returns the vertical velocity in geometric coordinates, [m/s]. - - For all grid points, the pressure vertical velocity in Pa/s is converted - to m/s via - w[m/s] =(approx) omega[Pa/s] / (-g*rho) - rho = p / R*T - with R = 287.058 JK-1kg-1, g = 9.80665 m2s-2. - (see p.13 of 'Introduction to circulating atmospheres' by Ian N. James). - - NOTE: Please check the resulting values, especially in the upper atmosphere! """ - rho = p / (287.058 * t) - return (omega / (-9.80665 * rho)) + omega = units.Quantity(omega, "Pa/s") + p = units.Quantity(p, "Pa") + t = units.Quantity(t, "K") + om_w = mpcalc.vertical_velocity(omega, p, t) + return om_w def flightlevel2pressure(flightlevel):