diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst index f6eeeaf3724..063cd061d36 100644 --- a/docs/_templates/overrides/metpy.calc.rst +++ b/docs/_templates/overrides/metpy.calc.rst @@ -70,6 +70,7 @@ Soundings bulk_shear bunkers_storm_motion cape_cin + ccl critical_angle cross_totals el diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index a18c4243d4e..a4c68c358a6 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -442,6 +442,109 @@ def _lcl_iter(p, p0, w, t): return lcl_p, globals()['dewpoint']._nounit(vapor_pressure._nounit(lcl_p, w)) +@exporter.export +@preprocess_and_wrap() +@check_units('[pressure]', '[temperature]', '[temperature]') +def ccl(pressure, temperature, dewpoint, height=None, mixed_layer_depth=None, which='top'): + r"""Calculate the convective condensation level (CCL) and convective temperature. + + This function is implemented directly based on the definition of the CCL, + as in [USAF1990]_, and finding where the ambient temperature profile intersects + the line of constant mixing ratio starting at the surface, using the surface dewpoint + or the average dewpoint of a shallow layer near the surface. + + Parameters + ---------- + pressure : `pint.Quantity` + Atmospheric pressure profile + + temperature : `pint.Quantity` + Temperature at the levels given by `pressure` + + dewpoint : `pint.Quantity` + Dewpoint at the levels given by `pressure` + + height : `pint.Quantity`, optional + Atmospheric heights at the levels given by `pressure`. + Only needed when specifying a mixed layer depth as a height. + + mixed_layer_depth : `pint.Quantity`, optional + The thickness of the mixed layer as a pressure or height above the bottom + of the layer (default None). + + which: str, optional + Pick which CCL value to return; must be one of 'top', 'bottom', or 'all'. + 'top' returns the lowest-pressure CCL (default), + 'bottom' returns the highest-pressure CCL, + 'all' returns every CCL in a `Pint.Quantity` array. + + Returns + ------- + `pint.Quantity` + CCL Pressure + + `pint.Quantity` + CCL Temperature + + `pint.Quantity` + Convective Temperature + + See Also + -------- + lcl, lfc, el + + Notes + ----- + Only functions on 1D profiles (not higher-dimension vertical cross sections or grids). + Since this function returns scalar values when given a profile, this will return Pint + Quantities even when given xarray DataArray profiles. + + Examples + -------- + >>> import metpy.calc as mpcalc + >>> from metpy.units import units + >>> pressure = [993, 957, 925, 886, 850, 813, 798, 732, 716, 700] * units.mbar + >>> temperature = [34.6, 31.1, 27.8, 24.3, 21.4, 19.6, 18.7, 13, 13.5, 13] * units.degC + >>> dewpoint = [19.6, 18.7, 17.8, 16.3, 12.4, -0.4, -3.8, -6, -13.2, -11] * units.degC + >>> ccl_p, ccl_t, t_c = mpcalc.ccl(pressure, temperature, dewpoint) + >>> ccl_p, t_c + (, ) + """ + pressure, temperature, dewpoint = _remove_nans(pressure, temperature, dewpoint) + + # If the mixed layer is not defined, take the starting dewpoint to be the + # first element of the dewpoint array and calculate the corresponding mixing ratio. + if mixed_layer_depth is None: + p_start, dewpoint_start = pressure[0], dewpoint[0] + vapor_pressure_start = saturation_vapor_pressure(dewpoint_start) + r_start = mixing_ratio(vapor_pressure_start, p_start) + + # Else, calculate the mixing ratio of the mixed layer. + else: + vapor_pressure_profile = saturation_vapor_pressure(dewpoint) + r_profile = mixing_ratio(vapor_pressure_profile, pressure) + r_start = mixed_layer(pressure, r_profile, height=height, + depth=mixed_layer_depth)[0] + + # rt_profile is the temperature-pressure profile with a fixed mixing ratio + rt_profile = globals()['dewpoint'](vapor_pressure(pressure, r_start)) + + x, y = find_intersections(pressure, rt_profile, temperature, + direction='increasing', log_x=True) + + # In the case of multiple CCLs, select which to return + if which == 'top': + x, y = x[-1], y[-1] + elif which == 'bottom': + x, y = x[0], y[0] + elif which not in ['top', 'bottom', 'all']: + raise ValueError(f'Invalid option for "which": {which}. Valid options are ' + '"top", "bottom", and "all".') + + x, y = x.to(pressure.units), y.to(temperature.units) + return x, y, dry_lapse(pressure[0], y, x).to(temperature.units) + + @exporter.export @preprocess_and_wrap() @check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]') diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index cc8e5827bfc..99d90ee940c 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -10,7 +10,7 @@ import xarray as xr from metpy.calc import (brunt_vaisala_frequency, brunt_vaisala_frequency_squared, - brunt_vaisala_period, cape_cin, cross_totals, density, dewpoint, + brunt_vaisala_period, cape_cin, ccl, cross_totals, density, dewpoint, dewpoint_from_relative_humidity, dewpoint_from_specific_humidity, dry_lapse, dry_static_energy, el, equivalent_potential_temperature, exner_function, gradient_richardson_number, InvalidSoundingError, @@ -399,6 +399,129 @@ def test_lcl_nans(): np.nan, 18.82281982535794]) * units.degC) +def test_ccl_basic(): + """First test of CCL calculation. Data: ILX, June 17 2022 00Z.""" + pressure = np.array([993.0, 984.0, 957.0, 948.0, 925.0, 917.0, 886.0, 868.0, 850.0, + 841.0, 813.0, 806.0, 798.0, 738.0, 732.0, 723.0, 716.0, 711.0, + 700.0, 623.0, 621.0, 582.0, 541.0, 500.0, 468.0]) * units.mbar + temperature = np.array([34.6, 33.7, 31.1, 30.1, 27.8, 27.1, 24.3, 22.6, 21.4, + 20.8, 19.6, 19.4, 18.7, 13.0, 13.0, 13.4, 13.5, 13.6, + 13.0, 5.2, 5.0, 1.5, -2.4, -6.7, -10.7]) * units.degC + dewpoint = np.array([19.6, 19.4, 18.7, 18.4, 17.8, 17.5, 16.3, 15.6, 12.4, 10.8, + -0.4, -3.6, -3.8, -5.0, -6.0, -15.6, -13.2, -11.4, -11.0, + -5.8, -6.2, -14.8, -24.3, -34.7, -38.1]) * units.degC + ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint) + assert_almost_equal(ccl_p, 763.006048 * units.mbar, 5) + assert_almost_equal(ccl_t, 15.429946 * units.degC, 5) + assert_almost_equal(t_c, 37.991498 * units.degC, 5) + + +def test_ccl_nans(): + """Tests CCL handles nans.""" + pressure = np.array([993.0, 984.0, 957.0, np.nan, 925.0, 917.0, np.nan, 868.0, 850.0, + 841.0, 813.0, 806.0, 798.0, 738.0, 732.0, 723.0, 716.0, 711.0, + 700.0, 623.0, 621.0, 582.0, 541.0, 500.0, 468.0]) * units.mbar + temperature = np.array([34.6, np.nan, 31.1, np.nan, 27.8, 27.1, 24.3, 22.6, 21.4, + 20.8, 19.6, 19.4, 18.7, 13.0, 13.0, 13.4, 13.5, 13.6, + 13.0, 5.2, 5.0, 1.5, -2.4, -6.7, -10.7]) * units.degC + dewpoint = np.array([19.6, 19.4, 18.7, np.nan, 17.8, 17.5, 16.3, 15.6, 12.4, 10.8, + -0.4, -3.6, -3.8, -5.0, -6.0, -15.6, -13.2, -11.4, -11.0, + -5.8, -6.2, -14.8, -24.3, -34.7, -38.1]) * units.degC + ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint) + assert_almost_equal(ccl_p, 763.006048 * units.mbar, 5) + assert_almost_equal(ccl_t, 15.429946 * units.degC, 5) + assert_almost_equal(t_c, 37.991498 * units.degC, 5) + + +def test_ccl_unit(): + """Tests CCL pressure and temperature is returned in the correct unit.""" + pressure = (np.array([993.0, 984.0, 957.0, 948.0, 925.0, 917.0, 886.0, 868.0, 850.0, + 841.0, 813.0, 806.0, 798.0, 738.0, 732.0, 723.0, 716.0, 711.0, + 700.0, 623.0, 621.0, 582.0, 541.0, 500.0, 468.0]) * 100) * units.Pa + temperature = (np.array([34.6, 33.7, 31.1, 30.1, 27.8, 27.1, 24.3, 22.6, 21.4, + 20.8, 19.6, 19.4, 18.7, 13.0, 13.0, 13.4, 13.5, 13.6, + 13.0, 5.2, 5.0, 1.5, -2.4, -6.7, -10.7]) + 273.15) * units.kelvin + dewpoint = (np.array([19.6, 19.4, 18.7, 18.4, 17.8, 17.5, 16.3, 15.6, 12.4, 10.8, + -0.4, -3.6, -3.8, -5.0, -6.0, -15.6, -13.2, -11.4, -11.0, + -5.8, -6.2, -14.8, -24.3, -34.7, -38.1]) + 273.15) * units.kelvin + + ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint) + assert_almost_equal(ccl_p, (763.006048 * 100) * units.Pa, 3) + assert_almost_equal(ccl_t, (15.429946 + 273.15) * units.kelvin, 3) + assert_almost_equal(t_c, (37.991498 + 273.15) * units.kelvin, 3) + + assert ccl_p.units == pressure.units + assert ccl_t.units == temperature.units + assert t_c.units == temperature.units + + +def test_multiple_ccl(): + """Tests the case where there are multiple CCLs. Data: BUF, May 18 2022 12Z.""" + pressure = np.array([992.0, 990.0, 983.0, 967.0, 950.0, 944.0, 928.0, 925.0, 922.0, + 883.0, 877.7, 858.0, 853.0, 850.0, 835.0, 830.0, 827.0, 826.0, + 813.6, 808.0, 799.0, 784.0, 783.3, 769.0, 760.0, 758.0, 754.0, + 753.0, 738.0, 725.7, 711.0, 704.0, 700.0, 685.0, 672.0, 646.6, + 598.6, 596.0, 587.0, 582.0, 567.0, 560.0, 555.0, 553.3, 537.0, + 526.0, 521.0, 519.0, 515.0, 500.0]) * units.mbar + temperature = np.array([6.8, 6.2, 7.8, 7.6, 7.2, 7.6, 6.6, 6.4, 6.2, 3.2, 2.8, 1.2, + 1.0, 0.8, -0.3, -0.1, 0.4, 0.6, 0.9, 1.0, 0.6, -0.3, -0.3, + -0.7, -1.5, -1.3, 0.2, 0.2, -1.1, -2.1, -3.3, -2.3, -1.7, 0.2, + -0.9, -3.0, -7.3, -7.5, -8.1, -8.3, -9.5, -10.1, -10.7, + -10.8, -12.1, -12.5, -12.7, -12.9, -13.5, -15.5]) * units.degC + dewpoint = np.array([5.1, 5.0, 4.2, 2.7, 2.2, 0.6, -2.4, -2.6, -2.8, -3.8, -3.6, + -3.1, -5.0, -4.2, -1.8, -4.3, -7.6, -6.4, -8.2, -9.0, -10.4, + -9.3, -9.6, -14.7, -11.5, -12.3, -25.8, -25.8, -19.1, -19.6, + -20.3, -42.3, -39.7, -46.8, -46.8, -46.7, -46.5, -46.5, + -52.1, -36.3, -47.5, -30.1, -29.7, -30.4, -37.1, -49.5, + -36.7, -28.9, -28.5, -22.5]) * units.degC + + ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint) + assert_almost_equal(ccl_p, 680.191653 * units.mbar, 5) + assert_almost_equal(ccl_t, -0.204408 * units.degC, 5) + assert_almost_equal(t_c, 30.8678258 * units.degC, 5) + + ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint, which='bottom') + assert_almost_equal(ccl_p, 886.835325 * units.mbar, 5) + assert_almost_equal(ccl_t, 3.500840 * units.degC, 5) + assert_almost_equal(t_c, 12.5020423 * units.degC, 5) + + ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint, which='all') + assert_array_almost_equal(ccl_p, np.array([886.835325, 680.191653]) * units.mbar, 5) + assert_array_almost_equal(ccl_t, np.array([3.500840, -0.204408]) * units.degC, 5) + assert_array_almost_equal(t_c, np.array([12.5020423, 30.8678258]) * units.degC, 5) + + +def test_ccl_with_ml(): + """Test CCL calculation with a specified mixed-layer depth.""" + pressure = np.array([992.0, 990.0, 983.0, 967.0, 950.0, 944.0, 928.0, 925.0, 922.0, + 883.0, 877.7, 858.0, 853.0, 850.0, 835.0, 830.0, 827.0, 826.0, + 813.6, 808.0, 799.0, 784.0, 783.3, 769.0, 760.0, 758.0, 754.0, + 753.0, 738.0, 725.7, 711.0, 704.0, 700.0, 685.0, 672.0, 646.6, + 598.6, 596.0, 587.0, 582.0, 567.0, 560.0, 555.0, 553.3, 537.0, + 526.0, 521.0, 519.0, 515.0, 500.0]) * units.mbar + temperature = np.array([6.8, 6.2, 7.8, 7.6, 7.2, 7.6, 6.6, 6.4, 6.2, 3.2, 2.8, 1.2, + 1.0, 0.8, -0.3, -0.1, 0.4, 0.6, 0.9, 1.0, 0.6, -0.3, -0.3, + -0.7, -1.5, -1.3, 0.2, 0.2, -1.1, -2.1, -3.3, -2.3, -1.7, 0.2, + -0.9, -3.0, -7.3, -7.5, -8.1, -8.3, -9.5, -10.1, -10.7, + -10.8, -12.1, -12.5, -12.7, -12.9, -13.5, -15.5]) * units.degC + dewpoint = np.array([5.1, 5.0, 4.2, 2.7, 2.2, 0.6, -2.4, -2.6, -2.8, -3.8, -3.6, + -3.1, -5.0, -4.2, -1.8, -4.3, -7.6, -6.4, -8.2, -9.0, -10.4, + -9.3, -9.6, -14.7, -11.5, -12.3, -25.8, -25.8, -19.1, -19.6, + -20.3, -42.3, -39.7, -46.8, -46.8, -46.7, -46.5, -46.5, + -52.1, -36.3, -47.5, -30.1, -29.7, -30.4, -37.1, -49.5, + -36.7, -28.9, -28.5, -22.5]) * units.degC + + ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint, + mixed_layer_depth=500 * units.m, which='all') + + assert_array_almost_equal(ccl_p, np.array( + [850.600930, 784.325312, 737.767377, 648.076147]) * units.mbar, 5) + assert_array_almost_equal(ccl_t, np.array( + [0.840118, -0.280299, -1.118757, -2.875716]) * units.degC, 5) + assert_array_almost_equal(t_c, np.array( + [13.146845, 18.661621, 22.896152, 32.081388]) * units.degC, 5) + + def test_lfc_basic(): """Test LFC calculation.""" levels = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar