diff --git a/metpy/calc/tests/test_thermo.py b/metpy/calc/tests/test_thermo.py index 675174cea41..96f99a49c5d 100644 --- a/metpy/calc/tests/test_thermo.py +++ b/metpy/calc/tests/test_thermo.py @@ -19,7 +19,7 @@ saturation_mixing_ratio, saturation_vapor_pressure, specific_humidity_from_mixing_ratio, - surface_based_cape_cin, vapor_pressure, + surface_based_cape_cin, thickness_hydrostatic, vapor_pressure, virtual_potential_temperature, virtual_temperature) from metpy.calc.thermo import _find_append_zero_crossings from metpy.testing import assert_almost_equal, assert_array_almost_equal, assert_nan @@ -744,3 +744,39 @@ def test_moist_static_energy(): """Tests the moist static energy calculation.""" mse = moist_static_energy(1000 * units.m, 25 * units.degC, 0.012 * units.dimensionless) assert_almost_equal(mse, 339.4594 * units('kJ/kg'), 6) + + +def test_thickness_hydrostatic(): + """Tests the thickness calculation for a moist layer.""" + pressure = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.hPa + temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.degC + mixing = np.array([0.01458, 0.00209, 0.00224, 0.00240, 0.00256, 0.00010]) + thickness = thickness_hydrostatic(pressure, temperature, mixing=mixing) + assert_almost_equal(thickness, 9892.07 * units.m, 2) + + +def test_thickness_hydrostatic_subset(): + """Tests the thickness calculation with a subset of the moist layer.""" + pressure = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.hPa + temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.degC + mixing = np.array([0.01458, 0.00209, 0.00224, 0.00240, 0.00256, 0.00010]) + thickness = thickness_hydrostatic(pressure, temperature, mixing=mixing, + bottom=850 * units.hPa, depth=150 * units.hPa) + assert_almost_equal(thickness, 1630.81 * units.m, 2) + + +def test_thickness_hydrostatic_isothermal(): + """Tests the thickness calculation for a dry isothermal layer at 0 degC.""" + pressure = np.arange(1000, 500 - 1e-10, -10) * units.hPa + temperature = np.zeros_like(pressure) * units.degC + thickness = thickness_hydrostatic(pressure, temperature) + assert_almost_equal(thickness, 5542.12 * units.m, 2) + + +def test_thickness_hydrostatic_isothermal_subset(): + """Tests the thickness calculation for a dry isothermal layer subset at 0 degC.""" + pressure = np.arange(1000, 500 - 1e-10, -10) * units.hPa + temperature = np.zeros_like(pressure) * units.degC + thickness = thickness_hydrostatic(pressure, temperature, bottom=850 * units.hPa, + depth=350 * units.hPa) + assert_almost_equal(thickness, 4242.68 * units.m, 2) diff --git a/metpy/calc/thermo.py b/metpy/calc/thermo.py index 87c5662f5a2..1214648137e 100644 --- a/metpy/calc/thermo.py +++ b/metpy/calc/thermo.py @@ -1508,3 +1508,75 @@ def moist_static_energy(heights, temperature, specific_humidity): """ return (dry_static_energy(heights, temperature) + Lv * specific_humidity).to('kJ/kg') + + +@exporter.export +@check_units('[pressure]', '[temperature]') +def thickness_hydrostatic(pressure, temperature, **kwargs): + r"""Calculate the thickness of a layer via the hypsometric equation. + + This thickness calculation uses the pressure and temperature profiles (and optionally + mixing ratio) via the hypsometric equation with virtual temperature adjustment + + .. math:: Z_2 - Z_1 = -\frac{R_d}{g} \int_{p_1}^{p_2} T_v d\ln p, + + which is based off of Equation 3.24 in [Hobbs2006]_. + + This assumes a hydrostatic atmosphere. + + Layer bottom and depth specified in pressure. + + Parameters + ---------- + pressure : `pint.Quantity` + Atmospheric pressure profile + temperature : `pint.Quantity` + Atmospheric temperature profile + mixing : `pint.Quantity`, optional + Profile of dimensionless mass mixing ratio. If none is given, virtual temperature + is simply set to be the given temperature. + 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`). + bottom : `pint.Quantity`, optional + The bottom of the layer in pressure. Defaults to the first observation. + depth : `pint.Quantity`, optional + The depth of the layer in hPa. Defaults to the full profile if bottom is not given, + and 100 hPa if bottom is given. + + Returns + ------- + `pint.Quantity` + The thickness of the layer in meters. + + See Also + -------- + pressure_to_height_std, virtual_temperature + + """ + mixing = kwargs.pop('mixing', None) + molecular_weight_ratio = kwargs.pop('molecular_weight_ratio', epsilon) + bottom = kwargs.pop('bottom', None) + depth = kwargs.pop('depth', None) + + # Get the data for the layer, conditional upon bottom/depth being specified and mixing + # ratio being given + if bottom is None and depth is None: + if mixing is None: + layer_p, layer_virttemp = pressure, temperature + else: + layer_p = pressure + layer_virttemp = virtual_temperature(temperature, mixing, molecular_weight_ratio) + else: + if mixing is None: + layer_p, layer_virttemp = get_layer(pressure, temperature, bottom=bottom, + depth=depth) + else: + layer_p, layer_temp, layer_w = get_layer(pressure, temperature, mixing, + bottom=bottom, depth=depth) + layer_virttemp = virtual_temperature(layer_temp, layer_w, molecular_weight_ratio) + + # Take the integral (with unit handling) and return the result in meters + return (- Rd / g * np.trapz(layer_virttemp.to('K'), x=np.log(layer_p / units.hPa)) * + units.K).to('m')