Repository: Unidata/MetPy

Issue:
More examples in the documentation
Give examples of each calculation and examples on the pages

Discussion:

Patch:
diff --git a/ci/requirements.txt b/ci/requirements.txt
index dc8e8ea4e9..da0eb2e5d0 100644
--- a/ci/requirements.txt
+++ b/ci/requirements.txt
@@ -7,4 +7,4 @@ pint==0.19.2
 pyproj==3.4.0
 scipy==1.9.1
 traitlets==5.4.0
-xarray==2022.6.0
+xarray==2022.9.0
diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst
index 88a9d5cdad..3262e11ba7 100644
--- a/docs/_templates/overrides/metpy.calc.rst
+++ b/docs/_templates/overrides/metpy.calc.rst
@@ -23,6 +23,7 @@ Dry Thermodynamics
       static_stability
       temperature_from_potential_temperature
       thickness_hydrostatic
+      weighted_continuous_average
 
 
 Moist Thermodynamics
@@ -70,6 +71,7 @@ Soundings
       bulk_shear
       bunkers_storm_motion
       cape_cin
+      ccl
       critical_angle
       cross_totals
       el
diff --git a/examples/calculations/Bulk_Shear.py b/examples/calculations/Bulk_Shear.py
new file mode 100644
index 0000000000..9b0b4bdde0
--- /dev/null
+++ b/examples/calculations/Bulk_Shear.py
@@ -0,0 +1,49 @@
+# Copyright (c) 2022 MetPy Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+"""
+==========
+Bulk Shear
+==========
+
+Use functions from `metpy.calc` as well as pint's unit support to perform calculations.
+
+The code below uses example data from our test suite to calculate the bulk shear over the
+lowest three kilometers for the provided sounding data.
+"""
+import pandas as pd
+
+from metpy.calc import bulk_shear, wind_components
+from metpy.cbook import get_test_data
+from metpy.units import units
+
+###########################################
+# Upper air data can be obtained using the siphon package, but for this example we will use
+# some of MetPy's sample data.
+
+# Set column names
+col_names = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed']
+
+# Read in test data using col_names
+df = pd.read_fwf(get_test_data('jan20_sounding.txt', as_file_obj=False),
+                 skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
+
+###########################################
+# Drop any rows with all NaN values for T, Td, winds
+df = df.dropna(subset=('temperature', 'dewpoint', 'direction', 'speed'),
+               how='all').reset_index(drop=True)
+
+###########################################
+# Isolate pressure, wind direction, wid speed, and height and add units
+p = df['pressure'].values * units.hPa
+wdir = df['direction'].values * units.degree
+sped = df['speed'].values * units.knot
+height = df['height'].values * units.meter
+
+###########################################
+# Calculate the u and v-components of the wind
+u, v = wind_components(sped, wdir)
+
+###########################################
+# Compute the bulk shear for the lowest three km
+print(bulk_shear(p, u, v, height, depth=3 * units.km, bottom=height[0]))
diff --git a/examples/calculations/Equivalent_Potential_Temperature.py b/examples/calculations/Equivalent_Potential_Temperature.py
new file mode 100644
index 0000000000..403cb2e8cf
--- /dev/null
+++ b/examples/calculations/Equivalent_Potential_Temperature.py
@@ -0,0 +1,57 @@
+# Copyright (c) 2022 MetPy Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+"""
+================================
+Equivalent Potential Temperature
+================================
+
+Use functions from `metpy.calc` as well as pint's unit support to perform calculations.
+
+The code below uses example data from our test suite to calculate the equivalent potential
+temperature over the provided sounding data and plots the values up to 300-hPa.
+"""
+import matplotlib.pyplot as plt
+import pandas as pd
+
+from metpy.calc import equivalent_potential_temperature
+from metpy.cbook import get_test_data
+from metpy.units import units
+
+###########################################
+# Upper air data can be obtained using the siphon package, but for this example we will use
+# some of MetPy's sample data.
+
+# Set column names
+col_names = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed']
+
+# Read in test data using col_names
+df = pd.read_fwf(get_test_data('jan20_sounding.txt', as_file_obj=False),
+                 skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
+
+###########################################
+# Drop any rows with all NaN values for T, Td, winds
+df = df.dropna(subset=('temperature', 'dewpoint', 'direction', 'speed'
+                       ), how='all').reset_index(drop=True)
+
+###########################################
+# Isolate pressure, temperature, dewpoint, and height and add units
+p = df['pressure'].values * units.hPa
+T = df['temperature'].values * units.degC
+Td = df['dewpoint'].values * units.degC
+h = df['height'].values * units.meter
+
+###########################################
+# Calculate the equivalent potential temperature for the whole sounding
+theta_e = equivalent_potential_temperature(p, T, Td)
+
+###########################################
+# Define a layer with pressure greater than 300 hPa
+layer = p > 300 * units.hPa
+
+###########################################
+# Use the layer defined above and plot the equivalent potential temperature
+plt.figure(figsize=(8, 8))
+plt.plot(theta_e[layer], h[layer])
+plt.title('Equivalent Potential Temperature (Kelvin)')
+plt.show()
diff --git a/examples/calculations/Mean_Pressure_Weighted.py b/examples/calculations/Mean_Pressure_Weighted.py
new file mode 100644
index 0000000000..195695a3e8
--- /dev/null
+++ b/examples/calculations/Mean_Pressure_Weighted.py
@@ -0,0 +1,44 @@
+# Copyright (c) 2022 MetPy Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+"""
+======================
+Mean Pressure Weighted
+======================
+
+Use `metpy.calc.mean_pressure_weighted` as well as pint's unit support to perform calculations.
+
+The code below uses example data from our test suite to calculate the pressure-weighted mean
+temperature over a depth of 500 hPa.
+"""
+import pandas as pd
+
+from metpy.calc import mean_pressure_weighted
+from metpy.cbook import get_test_data
+from metpy.units import units
+
+###########################################
+# Upper air data can be obtained using the siphon package, but for this example we will use
+# some of MetPy's sample data.
+
+# Set column names
+col_names = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed']
+
+# Read in test data using col_names
+df = pd.read_fwf(get_test_data('jan20_sounding.txt', as_file_obj=False),
+                 skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
+
+###########################################
+# Drop any rows with all NaN values for T, Td, winds
+df = df.dropna(subset=('temperature', 'dewpoint', 'direction', 'speed'),
+               how='all').reset_index(drop=True)
+
+###########################################
+# Isolate pressure, temperature, and height and add units
+p = df['pressure'].values * units.hPa
+T = df['temperature'].values * units.degC
+h = df['height'].values * units.meters
+
+###########################################
+# Calculate the mean pressure weighted temperature over a depth of 500 hPa
+print(mean_pressure_weighted(p, T, height=h, depth=500 * units.hPa))
diff --git a/examples/calculations/Mountain_Problem.py b/examples/calculations/Mountain_Problem.py
new file mode 100644
index 0000000000..396ec2d512
--- /dev/null
+++ b/examples/calculations/Mountain_Problem.py
@@ -0,0 +1,65 @@
+# Copyright (c) 2022 MetPy Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+"""
+================
+Mountain Problem
+================
+
+Use functions from `metpy.calc` to perform calculations for a rising and sinking parcel.
+
+The code below explores a common problem in meteorology where a parcel can be defined and
+lifted initiatlly dry adibatically until saturation is reached. It then ascends moist
+adiabatically to a desired level before descending back to the original level from which the
+parcel started.
+"""
+import numpy as np
+
+from metpy.calc import dry_lapse, lcl, moist_lapse
+from metpy.units import units
+
+###########################################
+# To set up the classic mountain problem, let's define that the parcel will start at 1000-hPa
+# and ascend to 700-hPa with an initial temperature of 25 Celsius and a dewpoint of 10 Celsius
+
+p = np.linspace(1000, 700, 301) * units.hPa
+T = 25 * units.degC
+Td = 10 * units.degC
+
+###########################################
+# We first need to determine the maximum level of dry ascent. For this we can use the LCL
+# function and retain the pressure level and temperature of the parcel at the LCL
+lclp, lclt = lcl(p[0], T, Td)
+print('Initial dry ascent yields:')
+print(f'  LCL Pressure: {lclp:.2f}')
+print(f'  LCL Temperature: {lclt:.2f}')
+print()
+
+###########################################
+# Knowing that the calculation of the dry ascent is accomplished by the LCL calculation, we
+# know how to begin our moist ascent. Begin by subsetting the pressure to begin at levels
+# less than or equal to the LCL pressure and use the moist_lapse to find the temperature
+# at the top of our ascent (700-hPa in our case)
+moist_ascent_p = p[p <= lclp]
+moist_ascent_t = moist_lapse(moist_ascent_p, lclt)
+print('After moist ascent:')
+print(f'  Temperature at top of ascent: {moist_ascent_t[-1]:.2f}')
+print()
+
+###########################################
+# Now to come "down the mountain" our parcel will warm dry adiabatically, so we can use the
+# dry_laspe function to descend from the lowest pressure (using [::-1] to reverse the order)
+# and convert our solution to Celsius
+dry_descent = dry_lapse(p[::-1], moist_ascent_t[-1]).to('degC')
+
+###########################################
+# Pulling it all together
+print(f'Starting Temperature: {T:.2f}')
+print(f'Starting Dewpoint: {Td:.2f}', end='\n\n')
+print(f'Final Temperature: {dry_descent[-1]:.2f}')
+print(f'Final Dewpoint: {moist_ascent_t[-1]:.2f}')
+
+###########################################
+# So as we expect the parcel has warmed and dried out through the combined dry and moist ascent
+# due to the release of latent heat and subsequent precipitating out of moisture from the
+# parcel
diff --git a/examples/calculations/Sounding_Calculations.py b/examples/calculations/Sounding_Calculations.py
new file mode 100644
index 0000000000..fad5281800
--- /dev/null
+++ b/examples/calculations/Sounding_Calculations.py
@@ -0,0 +1,210 @@
+# Copyright (c) 2022 MetPy Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+"""
+=============================
+Sounding Calculation Examples
+=============================
+
+Use functions from `metpy.calc` to perform a number of calculations using sounding data.
+
+The code below uses example data to perform many sounding calculations for a severe weather
+event on May 22, 2011 from the Norman, OK sounding.
+"""
+import numpy as np
+import pandas as pd
+
+import metpy.calc as mpcalc
+from metpy.cbook import get_test_data
+from metpy.units import units
+
+###########################################
+# Effective Shear Algorithm for use in Supercell Composite Calculation
+
+
+def effective_layer(p, t, td, h, height_layer=False):
+    """A function that determines the effective inflow layer for a convective sounding.
+
+    Uses the default values of Thompason et al. (2004) for CAPE (100 J/kg) and CIN (-250 J/kg).
+
+    Input:
+      - p: sounding pressure with units
+      - T: sounding temperature with units
+      - Td: sounding dewpoint temperature with units
+      - h: sounding heights with units
+
+    Returns:
+      - pbot/hbot, ptop/htop: pressure/height of the bottom level,
+                              pressure/height of the top level
+    """
+    from metpy.calc import cape_cin, parcel_profile
+    from metpy.units import units
+
+    pbot = None
+
+    for i in range(p.shape[0]):
+        prof = parcel_profile(p[i:], t[i], td[i])
+        sbcape, sbcin = cape_cin(p[i:], t[i:], td[i:], prof)
+        if sbcape >= 100 * units('J/kg') and sbcin > -250 * units('J/kg'):
+            pbot = p[i]
+            hbot = h[i]
+            bot_idx = i
+            break
+    if not pbot:
+        return None, None
+
+    for i in range(bot_idx + 1, p.shape[0]):
+        prof = parcel_profile(p[i:], t[i], td[i])
+        sbcape, sbcin = cape_cin(p[i:], t[i:], td[i:], prof)
+        if sbcape < 100 * units('J/kg') or sbcin < -250 * units('J/kg'):
+            ptop = p[i]
+            htop = h[i]
+            break
+
+    if height_layer:
+        return hbot, htop
+    else:
+        return pbot, ptop
+
+
+###########################################
+# Upper air data can be obtained using the siphon package, but for this example we will use
+# some of MetPy's sample data.
+col_names = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed']
+
+df = pd.read_fwf(get_test_data('20110522_OUN_12Z.txt', as_file_obj=False),
+                 skiprows=7, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
+
+# Drop any rows with all NaN values for T, Td, winds
+df = df.dropna(subset=('temperature', 'dewpoint', 'direction', 'speed'
+                       ), how='all').reset_index(drop=True)
+
+###########################################
+# Isolate needed variables from our data file and attach units
+p = df['pressure'].values * units.hPa
+T = df['temperature'].values * units.degC
+Td = df['dewpoint'].values * units.degC
+wdir = df['direction'].values * units.degree
+sped = df['speed'].values * units.knot
+height = df['height'].values * units.meter
+
+###########################################
+# Compute the wind components
+u, v = mpcalc.wind_components(sped, wdir)
+
+###########################################
+# Compute common sounding index parameters
+ctotals = mpcalc.cross_totals(p, T, Td)
+kindex = mpcalc.k_index(p, T, Td)
+showalter = mpcalc.showalter_index(p, T, Td)
+total_totals = mpcalc.total_totals_index(p, T, Td)
+vert_totals = mpcalc.vertical_totals(p, T)
+
+###########################################
+# Compture the parcel profile for a surface-based parcel
+prof = mpcalc.parcel_profile(p, T[0], Td[0])
+
+###########################################
+# Compute the corresponding LI, CAPE, CIN values for a surface parcel
+lift_index = mpcalc.lifted_index(p, T, prof)
+cape, cin = mpcalc.cape_cin(p, T, Td, prof)
+
+###########################################
+# Determine the LCL, LFC, and EL for our surface parcel
+lclp, lclt = mpcalc.lcl(p[0], T[0], Td[0])
+lfcp, _ = mpcalc.lfc(p, T, Td)
+el_pressure, _ = mpcalc.el(p, T, Td, prof)
+
+###########################################
+# Compute the characteristics of a mean layer parcel (50-hPa depth)
+ml_t, ml_td = mpcalc.mixed_layer(p, T, Td, depth=50 * units.hPa)
+ml_p, _, _ = mpcalc.mixed_parcel(p, T, Td, depth=50 * units.hPa)
+mlcape, mlcin = mpcalc.mixed_layer_cape_cin(p, T, prof, depth=50 * units.hPa)
+
+###########################################
+# Compute the characteristics of the most unstable parcel (50-hPa depth)
+mu_p, mu_t, mu_td, _ = mpcalc.most_unstable_parcel(p, T, Td, depth=50 * units.hPa)
+mucape, mucin = mpcalc.most_unstable_cape_cin(p, T, Td, depth=50 * units.hPa)
+
+###########################################
+# Compute the Bunkers Storm Motion vector and use to calculate the critical angle
+(u_storm, v_storm), *_ = mpcalc.bunkers_storm_motion(p, u, v, height)
+critical_angle = mpcalc.critical_angle(p, u, v, height, u_storm, v_storm)
+
+###########################################
+# Work on the calculations needed to compute the significant tornado parameter
+
+# Estimate height of LCL in meters from hydrostatic thickness
+new_p = np.append(p[p > lclp], lclp)
+new_t = np.append(T[p > lclp], lclt)
+lcl_height = mpcalc.thickness_hydrostatic(new_p, new_t)
+
+# Compute Surface-based CAPE
+sbcape, _ = mpcalc.surface_based_cape_cin(p, T, Td)
+
+# Compute SRH, given a motion vector toward the NE at 9.9 m/s
+*_, total_helicity = mpcalc.storm_relative_helicity(height, u, v, depth=1 * units.km,
+                                                    storm_u=u_storm, storm_v=v_storm)
+
+# Copmute Bulk Shear components and then magnitude
+ubshr, vbshr = mpcalc.bulk_shear(p, u, v, height=height, depth=6 * units.km)
+bshear = mpcalc.wind_speed(ubshr, vbshr)
+
+# Use all computed pieces to calculate the Significant Tornado parameter
+sig_tor = mpcalc.significant_tornado(sbcape, lcl_height,
+                                     total_helicity, bshear).to_base_units()
+
+###########################################
+# Compute the supercell composite parameter, if possible
+
+# Determine the top and bottom of the effective layer using our own function
+hbot, htop = effective_layer(p, T, Td, height, height_layer=True)
+
+# Perform the calculation of supercell composite if an effective layer exists
+if hbot:
+    esrh = mpcalc.storm_relative_helicity(height, u, v, depth=htop - hbot, bottom=hbot)
+    eubshr, evbshr = mpcalc.bulk_shear(p, u, v, height=height, depth=htop - hbot, bottom=hbot)
+    ebshear = mpcalc.wind_speed(eubshr, evbshr)
+
+    super_comp = mpcalc.supercell_composite(mucape, esrh[0], ebshear)
+else:
+    super_comp = np.nan
+
+###########################################
+# Print Important Sounding Parameters
+print('Important Sounding Parameters for KOUN on 22 Mary 2011 12 UTC')
+print()
+print(f'        CAPE: {cape:.2f}')
+print(f'         CIN: {cin:.2f}')
+print(f'LCL Pressure: {lclp:.2f}')
+print(f'LFC Pressure: {lfcp:.2f}')
+print(f' EL Pressure: {el_pressure:.2f}')
+print()
+print(f'   Lifted Index: {lift_index:.2f}')
+print(f'        K-Index: {kindex:.2f}')
+print(f'Showalter Index: {showalter:.2f}')
+print(f'   Cross Totals: {ctotals:.2f}')
+print(f'   Total Totals: {total_totals:.2f}')
+print(f'Vertical Totals: {vert_totals:.2f}')
+print()
+print('Mixed Layer - Lowest 50-hPa')
+print(f'     ML Temp: {ml_t:.2f}')
+print(f'     ML Dewp: {ml_td:.2f}')
+print(f'     ML CAPE: {mlcape:.2f}')
+print(f'      ML CIN: {mlcin:.2f}')
+print()
+print('Most Unstable - Lowest 50-hPa')
+print(f'     MU Temp: {mu_t:.2f}')
+print(f'     MU Dewp: {mu_td:.2f}')
+print(f' MU Pressure: {mu_p:.2f}')
+print(f'     MU CAPE: {mucape:.2f}')
+print(f'      MU CIN: {mucin:.2f}')
+print()
+print('Bunkers Storm Motion Vector')
+print(f'  u_storm: {u_storm:.2f}')
+print(f'  v_storm: {v_storm:.2f}')
+print(f'Critical Angle: {critical_angle:.2f}')
+print()
+print(f'Storm Relative Helicity: {total_helicity:.2f}')
+print(f'Significant Tornado Parameter: {sig_tor:.2f}')
+print(f'Supercell Composite Parameter: {super_comp:.2f}')
diff --git a/examples/calculations/Static_Stability.py b/examples/calculations/Static_Stability.py
new file mode 100644
index 0000000000..aa300c0d16
--- /dev/null
+++ b/examples/calculations/Static_Stability.py
@@ -0,0 +1,44 @@
+# Copyright (c) 2022 MetPy Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+"""
+================
+Static Stability
+================
+
+Use `metpy.calc.static_stability` as well as pint's unit support to perform calculations.
+
+The code below uses example data from our test suite to calculate the static stability for each
+layer of the provided sounding data.
+"""
+import pandas as pd
+
+from metpy.calc import static_stability
+from metpy.cbook import get_test_data
+from metpy.constants import Rd
+from metpy.units import units
+
+###########################################
+# Upper air data can be obtained using the siphon package, but for this example we will use
+# some of MetPy's sample data.
+
+# Set column names
+col_names = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed']
+
+# Read in test data using col_names
+df = pd.read_fwf(get_test_data('jan20_sounding.txt', as_file_obj=False),
+                 skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
+
+###########################################
+# Drop any rows with all NaN values for T, Td, winds
+df = df.dropna(subset=('temperature', 'dewpoint', 'direction', 'speed'
+                       ), how='all').reset_index(drop=True)
+
+###########################################
+# Isolate pressure and temperature and add units
+p = df['pressure'].values * units.hPa
+T = df['temperature'].values * units.degC
+
+###########################################
+# Calculate the static stability and convert to K/Pa units
+print((static_stability(p, T) * (p / Rd)).to('K/Pa'))
diff --git a/examples/calculations/Thickness_Hydrostatic.py b/examples/calculations/Thickness_Hydrostatic.py
new file mode 100644
index 0000000000..fdd1061535
--- /dev/null
+++ b/examples/calculations/Thickness_Hydrostatic.py
@@ -0,0 +1,51 @@
+# Copyright (c) 2022 MetPy Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+"""
+=====================
+Hydrostatic Thickness
+=====================
+
+Use functions from `metpy.calc` as well as pint's unit support to perform calculations.
+
+The code below uses example data from our test suite to calculate the hydrostatic thickness
+between the surface and 500-hPa level for the provided sounding data.
+"""
+import pandas as pd
+
+from metpy.calc import (mixing_ratio_from_relative_humidity, relative_humidity_from_dewpoint,
+                        thickness_hydrostatic)
+from metpy.cbook import get_test_data
+from metpy.units import units
+
+###########################################
+# Upper air data can be obtained using the siphon package, but for this example we will use
+# some of MetPy's sample data.
+
+# Set column names
+col_names = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed']
+
+# Read in test data using col_names
+df = pd.read_fwf(get_test_data('jan20_sounding.txt', as_file_obj=False),
+                 skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
+
+###########################################
+# Drop any rows with all NaN values for T, Td, winds
+df = df.dropna(subset=('temperature', 'dewpoint', 'direction', 'speed'
+                       ), how='all').reset_index(drop=True)
+
+###########################################
+# Isolate pressure, temperature, and dewpoint and add units
+p = df['pressure'].values * units.hPa
+T = df['temperature'].values * units.degC
+Td = df['dewpoint'].values * units.degC
+
+###########################################
+# Calculate the relative humidity to compute the mixing ratio
+rh = relative_humidity_from_dewpoint(T, Td)
+mixrat = mixing_ratio_from_relative_humidity(p, T, rh)
+
+###########################################
+# Calculate the thickness from the pressure, temperature, and mixing ratio
+# for the layer from the surface pressure to 500-hPa
+print(thickness_hydrostatic(p, T, mixing_ratio=mixrat, depth=p[0] - 500 * units.hPa))
diff --git a/examples/plots/Sounding_LCL_Dataset.py b/examples/plots/Sounding_LCL_Dataset.py
new file mode 100644
index 0000000000..ffc475b303
--- /dev/null
+++ b/examples/plots/Sounding_LCL_Dataset.py
@@ -0,0 +1,85 @@
+# Copyright (c) 2015-2022 MetPy Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+"""
+===========================
+Sounding as Dataset Example
+===========================
+
+Use MetPy to make a Skew-T LogP plot from an xarray Dataset after computing LCL parcel
+profile.
+
+This example makes a skewT diagram with specified special lines while storing the sounding
+data as an xarray dataset.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+import metpy.calc as mpcalc
+from metpy.cbook import get_test_data
+from metpy.plots import add_metpy_logo, SkewT
+from metpy.units import units
+
+###########################################
+# Upper air data can be obtained using the siphon package, but for this example we will use
+# some of MetPy's sample data.
+
+col_names = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed']
+
+df = pd.read_fwf(get_test_data('20110522_OUN_12Z.txt', as_file_obj=False),
+                 skiprows=7, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
+
+# Drop any rows with all NaN values for T, Td, winds
+df = df.dropna(subset=('temperature', 'dewpoint', 'direction', 'speed'),
+               how='all').reset_index(drop=True)
+
+###########################################
+# We will pull the data out of the example dataset into individual variables and
+# assign units.
+
+p = df['pressure'].values * units.hPa
+T = df['temperature'].values * units.degC
+Td = df['dewpoint'].values * units.degC
+wind_speed = df['speed'].values * units.knots
+wind_dir = df['direction'].values * units.degrees
+u, v = mpcalc.wind_components(wind_speed, wind_dir)
+
+ds = mpcalc.parcel_profile_with_lcl_as_dataset(p, T, Td)
+print(list(ds.variables))
+
+###########################################
+fig = plt.figure(figsize=(10, 8))
+skew = SkewT(fig, rotation=45)
+
+# Plot the data using the data from the xarray Dataset including the parcel temperature with
+# the LCL level included
+skew.plot(ds.isobaric, ds.ambient_temperature, 'r')
+skew.plot(ds.isobaric, ds.ambient_dew_point, 'g')
+skew.plot(ds.isobaric, ds.parcel_temperature.metpy.convert_units('degC'), 'black')
+
+# Plot the wind barbs from the original data
+skew.plot_barbs(p[::2], u[::2], v[::2])
+
+# Add the relevant special lines
+pressure = np.arange(1000, 499, -50) * units('hPa')
+mixing_ratio = np.array([0.1, 0.2, 0.4, 0.6, 1, 1.5, 2, 3, 4,
+                        6, 8, 10, 13, 16, 20, 25, 30, 36, 42]).reshape(-1, 1) * units('g/kg')
+
+skew.plot_dry_adiabats(t0=np.arange(233, 533, 10) * units.K, alpha=0.25,
+                       colors='orangered', linewidths=1)
+skew.plot_moist_adiabats(t0=np.arange(233, 400, 5) * units.K, alpha=0.25,
+                         colors='tab:green', linewidths=1)
+skew.plot_mixing_lines(pressure=pressure, mixing_ratio=mixing_ratio, linestyles='dotted',
+                       colors='dodgerblue', linewidths=1)
+skew.ax.set_ylim(1000, 100)
+
+# Add the MetPy logo!
+add_metpy_logo(fig, 350, 200)
+
+# Add some titles
+plt.title('KOUN', loc='left')
+plt.title('2011-05-22 12 UTC', loc='right')
+
+plt.show()
diff --git a/src/metpy/calc/indices.py b/src/metpy/calc/indices.py
index 102b0f895c..c2b8d2caa9 100644
--- a/src/metpy/calc/indices.py
+++ b/src/metpy/calc/indices.py
@@ -104,10 +104,10 @@ def mean_pressure_weighted(pressure, *args, height=None, bottom=None, depth=None
     Parameters
     ----------
     pressure : `pint.Quantity`
-        Atmospheric pressure profile
+        Atmospheric pressure profile.
 
     args : `pint.Quantity`
-        Parameters for which the pressure-weighted mean is to be calculated
+        Parameters for which the weighted-continuous mean is to be calculated.
 
     height : `pint.Quantity`, optional
         Heights from sounding. Standard atmosphere heights assumed (if needed)
@@ -120,12 +120,16 @@ def mean_pressure_weighted(pressure, *args, height=None, bottom=None, depth=None
         assumed to be the surface.
 
     depth: `pint.Quantity`, optional
-        Depth of the layer in meters or hPa
+        Depth of the layer in meters or hPa.
 
     Returns
     -------
     list of `pint.Quantity`
-        list of layer mean value for each profile in args
+        list of layer mean value for each profile in args.
+
+    See Also
+    --------
+    weighted_continuous_average
 
     Examples
     --------
@@ -157,6 +161,65 @@ def mean_pressure_weighted(pressure, *args, height=None, bottom=None, depth=None
     return [np.trapz(var_prof * pres_prof, x=pres_prof) / pres_int for var_prof in others]
 
 
+@exporter.export
+@preprocess_and_wrap()
+@check_units('[pressure]')
+def weighted_continuous_average(pressure, *args, height=None, bottom=None, depth=None):
+    r"""Calculate weighted-continuous mean of an arbitrary variable through a layer.
+
+    Layer top and bottom specified in height or pressure.
+
+    Formula based on that from [Holton2004]_ pg. 76 and the NCL function _wgt_vertical_n_
+
+    .. math::  WCA = \frac{\int_{p_s}^{p} A dp}{\int_{p_s}^{p} dp},
+
+    * :math:`WCA` is the weighted continuous average of a variable.
+    * :math:`p` is the bottom pressure level.
+    * :math:`p_s` is the top pressure level.
+    * :math:`A` is the variable whose weighted continuous average is being calculated.
+
+    Parameters
+    ----------
+    pressure : `pint.Quantity`
+        Atmospheric pressure profile.
+
+    args : `pint.Quantity`
+        Parameters for which the weighted-continuous mean is to be calculated.
+
+    height : `pint.Quantity`, optional
+        Heights from sounding. Standard atmosphere heights assumed (if needed)
+        if no heights are given.
+
+    bottom: `pint.Quantity`, optional
+        The bottom of the layer in either the provided height coordinate
+        or in pressure. Don't provide in meters AGL unless the provided
+        height coordinate is meters AGL. Default is the first observation,
+        assumed to be the surface.
+
+    depth: `pint.Quantity`, optional
+        Depth of the layer in meters or hPa.
+
+    Returns
+    -------
+    list of `pint.Quantity`
+        list of layer mean value for each profile in args.
+
+    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.
+
+    """
+    # Split pressure profile from other variables to average
+    pres_prof, *others = get_layer(
+        pressure, *args, height=height, bottom=bottom, depth=depth
+    )
+
+    return [np.trapz(var_prof, x=pres_prof) / (pres_prof[-1] - pres_prof[0])
+            for var_prof in others]
+
+
 @exporter.export
 @preprocess_and_wrap()
 @check_units('[pressure]', '[speed]', '[speed]', '[length]')
diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py
index 44bb3aa5d8..d53e5ef538 100644
--- a/src/metpy/calc/thermo.py
+++ b/src/metpy/calc/thermo.py
@@ -473,6 +473,109 @@ def lcl(pressure, temperature, dewpoint, max_iters=50, eps=1e-5):
     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]')
diff --git a/src/metpy/static-data-manifest.txt b/src/metpy/static-data-manifest.txt
index b5e7c06bc2..d494e89d0c 100644
--- a/src/metpy/static-data-manifest.txt
+++ b/src/metpy/static-data-manifest.txt
@@ -1,3 +1,4 @@
+20110522_OUN_12Z.txt 9a831910c077173deede5432b147438c033a460b06a5305697a5a481f9073e21
 2020010600_sao.wmo e8f327f32d1765ae02c2c34c9393287a853a62afbf1987bf9a4f32894c7f24bd
 AK-REGIONAL_8km_3.9_20160408_1445.gini 4919b80abe95a765a42a6e3904ef1e38b8be804dc52276ef5db6b587f85abd55
 CAM_test.nc ab62e6dca907c34eb44b842c67eb3eb817127306feb3524fa6d83670b628727c
diff --git a/staticdata/20110522_OUN_12Z.txt b/staticdata/20110522_OUN_12Z.txt
new file mode 100644
index 0000000000..bde30fa95f
--- /dev/null
+++ b/staticdata/20110522_OUN_12Z.txt
@@ -0,0 +1,77 @@
+72357 OUN Norman Observations at 12Z 22 May 2011
+
+-----------------------------------------------------------------------------
+   PRES   HGHT   TEMP   DWPT   RELH   MIXR   DRCT   SKNT   THTA   THTE   THTV
+    hPa     m      C      C      %    g/kg    deg   knot     K      K      K 
+-----------------------------------------------------------------------------
+ 1000.0     36                                                               
+  966.0    345   22.2   21.0     93  16.50    180      7  298.3  346.4  301.2
+  953.0    462   21.4   20.7     96  16.42    184     16  298.6  346.6  301.6
+  936.9    610   20.8   20.5     98  16.52    190     28  299.5  347.9  302.5
+  925.0    720   20.4   20.4    100  16.61    200     33  300.2  349.0  303.1
+  904.5    914   19.3   19.3    100  15.81    205     36  300.9  347.6  303.8
+  896.0    995   18.8   18.8    100  15.49    209     38  301.3  347.0  304.1
+  890.0   1054   20.0   20.0    100  16.84    212     40  303.1  353.2  306.1
+  886.0   1093   22.2   19.0     82  15.87    214     41  305.7  353.5  308.6
+  873.3   1219   23.2   13.3     54  11.12    220     45  308.0  342.0  310.1
+  873.0   1222   23.2   13.2     53  11.04    220     45  308.1  341.8  310.1
+  850.0   1454   22.0    6.0     35   6.94    210     37  309.2  330.8  310.5
+  846.0   1495   21.8    3.8     31   5.97    210     37  309.4  328.1  310.5
+  813.8   1829   19.2   -1.7     24   4.16    210     34  310.1  323.4  310.8
+  802.0   1955   18.2   -3.8     22   3.62    214     32  310.3  322.0  311.0
+  785.0   2134   16.5   -4.7     23   3.46    220     29  310.4  321.6  311.1
+  757.1   2438   13.7   -6.2     25   3.20    230     30  310.6  321.0  311.2
+  730.1   2743   10.9   -7.7     26   2.96    240     31  310.7  320.4  311.3
+  700.0   3096    7.6   -9.4     29   2.69    245     30  310.9  319.7  311.4
+  653.3   3658    2.3  -10.9     37   2.56    250     26  311.1  319.5  311.6
+  639.0   3839    0.6  -11.4     40   2.52    251     31  311.1  319.4  311.6
+  606.0   4262   -2.9  -12.9     46   2.35    255     42  311.8  319.6  312.3
+  605.6   4267   -2.9  -12.9     46   2.35    255     42  311.9  319.6  312.3
+  584.0   4555   -4.5  -14.5     46   2.14    255     54  313.3  320.4  313.7
+  582.7   4572   -4.4  -14.4     46   2.17    255     55  313.6  320.9  314.0
+  582.0   4582   -4.3  -14.3     46   2.18    255     55  313.8  321.1  314.2
+  577.0   4650   -3.7  -24.7     18   0.90    256     52  315.3  318.5  315.5
+  571.0   4733   -3.3  -36.3      6   0.30    258     49  316.7  317.9  316.8
+  561.0   4873   -3.9  -31.9      9   0.47    260     44  317.6  319.4  317.7
+  560.7   4877   -3.9  -31.8      9   0.48    260     44  317.6  319.4  317.7
+  539.4   5182   -6.3  -27.4     17   0.75    255     41  318.4  321.1  318.5
+  539.0   5187   -6.3  -27.3     17   0.76    255     41  318.4  321.1  318.5
+  500.0   5770  -11.1  -29.1     21   0.69    260     48  319.4  322.0  319.6
+  478.9   6096  -13.7  -31.3     21   0.59    265     46  320.2  322.4  320.3
+  453.0   6515  -17.1  -34.1     21   0.47    263     44  321.1  322.8  321.1
+  443.0   6681  -18.3  -32.3     28   0.58    263     43  321.6  323.8  321.7
+  406.3   7315  -23.9  -37.0     29   0.39    260     40  322.4  323.9  322.5
+  400.0   7430  -24.9  -37.9     29   0.37    255     38  322.5  323.9  322.6
+  389.3   7620  -26.6  -39.3     29   0.33    250     35  322.8  324.0  322.8
+  327.3   8839  -37.9  -48.1     34   0.15    230     28  323.7  324.3  323.8
+  313.4   9144  -40.7  -50.3     35   0.12    235     26  323.9  324.4  323.9
+  300.0   9449  -43.5  -52.5     36   0.10    230     24  323.9  324.4  324.0
+  286.0   9769  -46.3  -55.3     35   0.07    240     28  324.4  324.7  324.4
+  250.0  10650  -52.1  -62.1     29   0.04    255     41  328.5  328.6  328.5
+  249.0  10676  -52.3  -62.3     29   0.04    255     41  328.6  328.7  328.6
+  220.0  11473  -54.1  -64.1     28   0.03    261     54  337.6  337.8  337.6
+  210.0  11770  -55.9  -65.9     27   0.03    263     58  339.3  339.4  339.3
+  200.0  12080  -56.5  -66.5     27   0.02    265     63  343.1  343.2  343.1
+  197.0  12176  -56.5  -66.5     27   0.03    265     64  344.6  344.7  344.6
+  196.5  12192  -56.5  -66.5     27   0.03    265     64  344.9  345.0  344.9
+  190.0  12405  -56.5  -66.5     27   0.03    265     63  348.2  348.3  348.2
+  181.0  12711  -57.9  -67.9     27   0.02    265     61  350.8  350.9  350.8
+  173.0  12996  -56.9  -66.9     27   0.03    264     60  357.0  357.1  357.0
+  159.0  13526  -59.3  -69.3     26   0.02    261     58  361.6  361.8  361.6
+  155.0  13685  -58.9  -68.9     26   0.02    260     57  365.0  365.1  365.0
+  154.2  13716  -59.0  -69.0     26   0.02    260     57  365.3  365.4  365.3
+  150.0  13890  -59.5  -69.5     26   0.02    260     51  367.4  367.5  367.4
+  148.0  13974  -59.9  -69.9     26   0.02    254     48  368.1  368.2  368.1
+  146.9  14021  -59.2  -69.2     26   0.02    250     47  370.0  370.2  370.0
+  146.0  14059  -58.7  -68.7     26   0.03    252     46  371.6  371.7  371.6
+  142.0  14233  -57.9  -67.9     27   0.03    262     40  376.0  376.1  376.0
+  140.0  14323  -56.7  -66.7     27   0.03    267     38  379.6  379.8  379.6
+  137.0  14460  -56.5  -66.5     27   0.04    275     33  382.3  382.5  382.3
+  133.3  14630  -57.7  -67.7     27   0.03    285     28  383.2  383.4  383.2
+  127.0  14935  -59.8  -69.8     26   0.02    275     17  384.8  384.9  384.8
+  126.0  14986  -60.1  -70.1     26   0.02    273     17  385.1  385.2  385.1
+  120.9  15240  -61.0  -71.0     25   0.02    260     16  387.9  388.1  387.9
+  111.0  15771  -62.9  -72.9     25   0.02    233     18  394.0  394.1  394.0
+  109.0  15882  -64.3  -73.3     28   0.02    227     18  393.4  393.5  393.4
+  104.0  16170  -63.3  -73.3     25   0.02    212     19  400.6  400.8  400.6
+  100.0  16410  -64.3  -74.3     24   0.02    200     20  403.2  403.3  403.2


Test patch:
diff --git a/tests/calc/test_indices.py b/tests/calc/test_indices.py
index 132e03d3dc..744edeb63e 100644
--- a/tests/calc/test_indices.py
+++ b/tests/calc/test_indices.py
@@ -11,7 +11,7 @@ import xarray as xr
 
 from metpy.calc import (bulk_shear, bunkers_storm_motion, critical_angle,
                         mean_pressure_weighted, precipitable_water, significant_tornado,
-                        supercell_composite)
+                        supercell_composite, weighted_continuous_average)
 from metpy.testing import assert_almost_equal, assert_array_almost_equal, get_upper_air_data
 from metpy.units import concatenate, units
 
@@ -118,6 +118,43 @@ def test_mean_pressure_weighted_elevated():
     assert_almost_equal(v, 1.7392601775853547 * units('m/s'), 7)
 
 
+def test_weighted_continuous_average():
+    """Test pressure-weighted mean wind function with vertical interpolation."""
+    data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
+    u, v = weighted_continuous_average(data['pressure'],
+                                       data['u_wind'],
+                                       data['v_wind'],
+                                       height=data['height'],
+                                       depth=6000 * units('meter'))
+    assert_almost_equal(u, 6.644137766806087 * units('m/s'), 7)
+    assert_almost_equal(v, 6.900543760612305 * units('m/s'), 7)
+
+
+@pytest.mark.xfail(reason='hgrecco/pint#1593')
+def test_weighted_continuous_average_temperature():
+    """Test pressure-weighted mean temperature function with vertical interpolation."""
+    data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
+    t, = weighted_continuous_average(data['pressure'],
+                                     data['temperature'],
+                                     height=data['height'],
+                                     depth=6000 * units('meter'))
+    # Commenting out since it won't run until the above can run without error
+    # assert_almost_equal(t, 279.3275828240889 * units('kelvin'), 7)
+
+
+def test_weighted_continuous_average_elevated():
+    """Test pressure-weighted mean wind function with a base above the surface."""
+    data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
+    u, v = weighted_continuous_average(data['pressure'],
+                                       data['u_wind'],
+                                       data['v_wind'],
+                                       height=data['height'],
+                                       depth=3000 * units('meter'),
+                                       bottom=data['height'][0] + 3000 * units('meter'))
+    assert_almost_equal(u, 8.279561625494285 * units('m/s'), 7)
+    assert_almost_equal(v, 1.616638856115755 * units('m/s'), 7)
+
+
 def test_precipitable_water_xarray():
     """Test precipitable water with xarray input."""
     data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py
index 1766b91bb2..d3249935f2 100644
--- a/tests/calc/test_thermo.py
+++ b/tests/calc/test_thermo.py
@@ -10,7 +10,7 @@ import pytest
 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