Implement convective condensation level (CCL)

Fix docstring for ccl function Fix inline comment for ccl
Unidata · Jun 29, 2022 · 39b222f · 39b222f
1 parent c7ae7cf
commit 39b222f
Show file tree

Hide file tree

Showing 3 changed files with 223 additions and 1 deletion.
diff --git 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
@@ -442,6 +442,108 @@ 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).
+
+    This function is implemented by simplifying the mixing ratio and
+    the saturation vapor pressure equation to
+
+    .. math:: A(P) = ln(\frac{rP}{6.112(r+\epsilon)}) = \frac{17.67T}{T+243.5}
+
+    where `P` is pressure, and `T` is temperature in degrees C.
+
+    This gives us a direct relationship of how temperature changes with pressure
+    in the atmosphere with a constant mixing ratio.
+
+    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`
+
+    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.
+
+    """
+    pressure, temperature, dewpoint = _remove_nans(pressure, temperature, dewpoint)
+
+    dewpoint, temperature = dewpoint.to(units.degC), temperature.to(units.degC)
+    pressure = pressure.to(units.hPa)
+
+    # 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]
+
+    a_p = np.log(r_start * pressure / (
+        mpconsts.default.sat_pressure_0c * (mpconsts.nounit.epsilon + r_start)))
+    # rt_profile is the temperature-pressure profile with a fixed mixing ratio
+    rt_profile = units.Quantity((243.5 * a_p / (17.67 - a_p)).magnitude, 'degC')
+
+    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('Invalid option for "which". Valid options are "top", "bottom", '
+                         'and "all".')
+
+    x, y = x.to(pressure.units), y.to(temperature.units)
+    return x, y, potential_temperature(x, y).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
@@ -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,125 @@ 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, 38.616596 * 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, 38.616596 * 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, (38.616596 + 273.15) * units.kelvin, 3)
+
+
+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, 31.566319 * 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, 13.158340 * 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([13.158340, 31.566319]) * 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.804624, 19.332071, 23.576331, 32.782670]) * units.degC, 5)
+
+
 def test_lfc_basic():
     """Test LFC calculation."""
     levels = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar