Skip to content

Commit

Permalink
Merge pull request #2550 from Z-Richard/ccl
Browse files Browse the repository at this point in the history 
Implement convective condensation level (CCL)
  • Loading branch information
@dopplershift
dopplershift authored Sep 30, 2022
2 parents c8cecd4 + 79805c5 commit 0f079df
Show file tree
Hide file tree
Showing 3 changed files with 228 additions and 1 deletion.
1 change: 1 addition & 0 deletions docs/_templates/overrides/metpy.calc.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,7 @@ Soundings
bulk_shear
bunkers_storm_motion
cape_cin
ccl
critical_angle
cross_totals
el
Expand Down
103 changes: 103 additions & 0 deletions src/metpy/calc/thermo.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -473,6 +473,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
(<Quantity(758.348093, 'millibar')>, <Quantity(38.4336274, 'degree_Celsius')>)
"""
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]')
Expand Down
125 changes: 124 additions & 1 deletion tests/calc/test_thermo.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 0f079df

Please sign in to comment.