diff --git a/pyproject.toml b/pyproject.toml
index ce8acfd..8fb4b02 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -49,6 +49,7 @@ markers = [
     "mu_cape",
     "parcel_profile",
     "regression",
+    "broadcasting",
 ]
 
 pythonpath = ["src"]
diff --git a/src/lib/libthermo.cpp b/src/lib/libthermo.cpp
index 5da2669..b2bb204 100644
--- a/src/lib/libthermo.cpp
+++ b/src/lib/libthermo.cpp
@@ -189,7 +189,7 @@ constexpr T lcl<T>::wet_bulb_temperature(const T pressure, const T step) noexcep
 }
 
 template <floating T>
-constexpr size_t lcl<T>::index(const T pressure[], size_t size) noexcept {
+constexpr size_t lcl<T>::index(const T pressure[], const size_t size) noexcept {
     return index_pressure(pressure, this->pressure, size);
 }
 
diff --git a/src/nzthermo/__init__.py b/src/nzthermo/__init__.py
index b1af963..42b936c 100644
--- a/src/nzthermo/__init__.py
+++ b/src/nzthermo/__init__.py
@@ -48,6 +48,7 @@
     "downdraft_cape",
     "most_unstable_parcel",
     "most_unstable_parcel_index",
+    "most_unstable_cape_cin",
     # .utils
     "timeseries",
 ]
@@ -104,5 +105,6 @@
     mixing_ratio,
     most_unstable_parcel,
     most_unstable_parcel_index,
+    most_unstable_cape_cin,
 )
 from .utils import timeseries
diff --git a/src/nzthermo/_core.pyx b/src/nzthermo/_core.pyx
index 9f993bc..8f2e56e 100644
--- a/src/nzthermo/_core.pyx
+++ b/src/nzthermo/_core.pyx
@@ -18,6 +18,7 @@ template <typename T>
 T fn(T ...){...}
 ```
 """
+import warnings
 
 from cython.parallel cimport parallel, prange
 from cython.view cimport array as cvarray
@@ -110,21 +111,98 @@ cdef cvarray nzarray((size_t, size_t) shape, size_t size):
     )
 
 
-cdef pressure_mode(
-    np.ndarray pressure,
-    np.ndarray temperature,
-    np.ndarray dewpoint,
+def broadcast_and_nanmask(
+    np.ndarray pressure not None,
+    np.ndarray temperature not None, 
+    np.ndarray dewpoint, 
+    long ndim = 2,
+    object where = None,
+    object dtype = None,
 ):
-    if pressure.ndim == 1:
+    """
+    TODO: ideally this function needs to be exposed to the `core.py` or implmented as part of the
+    decorator, so that it is only ever called once. As part of the masking operation it would be
+    useful to properly mask and sort the data once per function call.
+    
+
+    >>> where &= np.isfinite(temperature) & np.isfinite(dewpoint) & np.isfinite(pressure)
+
+    Validates the input arrays and determines the broadcast mode. by brodcasting the ``pressure``, 
+    ``temperature``, and ``dewpoint`` arrays to a suitable shape for the calling function.
+
+    In some cases temperature and dewpoint arrays can be 1D, if this is the case the ``ndim``
+    argument must be set to 1.
+
+    masks the data with nan values and sorts in descending pressure order (bottom -> top) with nan
+    values at the end, this allows us to short circuit some computations as some of the functions
+    will stop at the first nan value.
+
+    Insures that all arays are in ``C_CONTIGUOUS`` memory layout. This is important for the Cython
+    memoryview to work correctly.
+    """
+    cdef:
+        size_t N, Z
+        BroadcastMode mode
+
+    if temperature.ndim != ndim:
+        raise ValueError(f"temperature must be a {ndim}D array.")
+    elif dewpoint.ndim != ndim:
+        raise ValueError(f"dewpoint must be a {ndim}D array.")
+    elif ndim == 2:
+        if temperature.shape[0] != dewpoint.shape[0]:
+            raise ValueError("temperature and dewpoint must have the same number of rows.")
+        elif temperature.shape[1] != dewpoint.shape[1]:
+            raise ValueError("temperature and dewpoint must have the same number of columns.")
+    elif ndim == 1:
+        if temperature.shape[0] != dewpoint.shape[0]:
+            raise ValueError("temperature and dewpoint must have the same number of elements.")
+
+    if pressure.ndim == 1 or pressure.shape[0] == 1:
         pressure = pressure.reshape(1, -1)
         mode = BROADCAST
     else:
         mode = MATRIX
+
+    if where is not None:
+        if not isinstance(where, np.ndarray):
+            raise ValueError("where must be a numpy array.")
+        
+        
+        mode = MATRIX
+        N = temperature.shape[0]
+        Z = temperature.shape[1]
+
+        where = np.atleast_2d(where)
+        # where &= np.isfinite(temperature) & np.isfinite(dewpoint) & np.isfinite(pressure)
+
+        pressure = np.broadcast_to(pressure.squeeze(), (N, Z))
+        pressure = np.where(where, pressure, -np.inf)
+
+        sort = np.arange(N)[:, np.newaxis], np.argsort(pressure, axis=1)
+        pressure = pressure[sort][:, ::-1]
+        temperature = temperature[sort][:, ::-1]
+        dewpoint = dewpoint[sort][:, ::-1]
+        m = np.isneginf(pressure)
+        pressure[m] = np.nan
+        temperature[m] = np.nan
+        dewpoint[m] = np.nan
+
+    if dtype is None:
+        dtype = temperature.dtype
+    
+    if dtype != np.float32 and dtype != np.float64:
+        
+        warnings.warn(f"the dtype {dtype} is not supported, defaulting to np.float64.")
+        dtype = np.float64
     
-    return (pressure, temperature, dewpoint), mode
+    pressure = np.ascontiguousarray(pressure, dtype=dtype)
+    temperature = np.ascontiguousarray(temperature, dtype=dtype)
+    dewpoint = np.ascontiguousarray(dewpoint, dtype=dtype)
 
+    return (pressure, temperature, dewpoint), mode, dtype
 
-# need to figuoure out a way to possibly pass in **kwargs maybe via a options struct
+
+# need to figure out a way to possibly pass in **kwargs maybe via a options struct
 ctypedef T (*Dispatch)(const T*, const T*, const T*, size_t) noexcept nogil
 
 cdef T[:] dispatch(
@@ -155,7 +233,6 @@ cdef T[:] dispatch(
 
     N, Z = temperature.shape[0], pressure.shape[1]
     out = np.empty((N,), dtype=np.dtype(f"f{pressure.itemsize}"))
-
     with nogil:
         if BROADCAST is mode:
             for i in prange(N, schedule='dynamic'):
@@ -200,14 +277,12 @@ cdef T[:, :] moist_lapse_2d(
     T[:] reference_pressure, 
     T[:] temperature, 
     BroadcastMode mode,
-    
 ) noexcept:
     cdef:
         size_t N, Z, i
         T[:, :] out
 
     N, Z = temperature.shape[0], pressure.shape[1]
-
     out = nzarray((N, Z), pressure.itemsize)
     with nogil, parallel():
         if BROADCAST is mode:
@@ -349,7 +424,6 @@ def moist_lapse(
         raise ValueError("Unable to determine the broadcast mode.")
 
     Z = pressure.shape[1]
-
     out = np.empty((N, Z), dtype=dtype)
     if np.float32 == dtype:
         out[...] = moist_lapse_2d[float](
@@ -382,20 +456,18 @@ cdef void parcel_profile_1d(
 ) noexcept nogil:
     cdef:
         size_t Z, i, stop
-        T p0, t0, reference_pressure, next_pressure
+        T p0, reference_pressure, next_pressure
         C.lcl[T] lcl
 
-    Z = pressure.shape[0]
+    Z = out.shape[0]
     p0 = pressure[0]
-    t0 = out[0] = temperature
+    lcl = C.lcl[T](pressure[0], temperature, dewpoint)
 
-    lcl = C.lcl[T](p0, t0, dewpoint)
-    # [dry ascent] 
+    # [dry ascent]
     # stop the dry ascent at the LCL
     stop = lcl.index(&pressure[0], Z)
     for i in prange(0, stop, schedule='dynamic'): # parallelize the dry ascent
-        out[i] = C.dry_lapse(pressure[i], p0, t0)
-    
+        out[i] = C.dry_lapse(pressure[i], p0, temperature)
 
     # [ moist ascent ]
     if stop != Z:
@@ -407,16 +479,13 @@ cdef T[:, :] parcel_profile_2d(
     T[:] temperature,
     T[:] dewpoint,
     BroadcastMode mode,
-    T step,
-    T eps,
-    size_t max_iters,
 ) noexcept:
     cdef:
         size_t N, Z, i
         T[:, :] out
 
     N, Z = temperature.shape[0], pressure.shape[1]
-    out = nzarray((N, Z), pressure.itemsize)
+    out = np.empty((N, Z), dtype=np.dtype(f"f{pressure.itemsize}"))
     with nogil, parallel():
         if BROADCAST is mode:
             for i in prange(N, schedule='dynamic'):
@@ -425,7 +494,7 @@ cdef T[:, :] parcel_profile_2d(
         else: # MATRIX
             for i in prange(N, schedule='dynamic'):
                 parcel_profile_1d(out[i], pressure[i, :], temperature[i], dewpoint[i])
- 
+
     return out
 
 
@@ -435,44 +504,22 @@ def parcel_profile(
     np.ndarray dewpoint not None,
     *,
     np.ndarray where = None,
-    double step = 1000.0,
-    double eps = 0.1,
-    size_t max_iters = 5,
 ):
     cdef:
         size_t N, Z
         BroadcastMode mode
         np.ndarray out
 
-    (pressure, temperature, dewpoint), mode = pressure_mode(pressure, temperature, dewpoint)
-    if where is not None:
-        raise NotImplementedError("where argument is not supported.")
-
+    (pressure, temperature, dewpoint), mode, dtype = broadcast_and_nanmask(
+        pressure, temperature, dewpoint, ndim=1, where=where
+    )
 
     N, Z = temperature.shape[0], pressure.shape[1]
-    
-    out = np.empty((N, Z), dtype=pressure.dtype)
-    if pressure.dtype == np.float64:
-        out[...] = parcel_profile_2d[double](
-            pressure.astype(np.float64),
-            temperature.astype(np.float64),
-            dewpoint.astype(np.float64),
-            mode,
-            step,
-            eps,
-            max_iters,
-        )
+    out = np.empty((N, Z), dtype=dtype)
+    if dtype == np.float64:
+        out[...] = parcel_profile_2d[double](pressure, temperature, dewpoint, mode)
     else:
-        out[...] = parcel_profile_2d[float](
-            pressure.astype(np.float32),
-            temperature.astype(np.float32),
-            dewpoint.astype(np.float32),
-            mode,
-            <float>step,
-            <float>eps,
-            max_iters,
-        )
-    
+        out[...] = parcel_profile_2d[float](pressure, temperature, dewpoint, mode)
 
     return out
 
@@ -501,30 +548,40 @@ cdef void parcel_profile_with_lcl_1d(
     lcl = C.lcl[T](p0, t0, td0)
 
     # [dry ascent] .. parcel temperature from the surface up to the LCL ..
+
     stop = lcl.index(&pressure[0], Z)
-    ep[:stop] = pressure[:stop]
-    et[:stop] = temperature[:stop]
-    etd[:stop] = dewpoint[:stop]
-    for i in prange(0, stop, schedule='dynamic'): # parallelize the dry ascent
-        pt[i] = C.dry_lapse(pressure[i], p0, t0)
-
-    # [ lcl ]
-    ep[stop] = lcl.pressure
-    et[stop] = C.linear_interpolate(
-        lcl.pressure, 
-        pressure[stop - 1], 
-        pressure[stop], 
-        temperature[stop - 1], 
-        temperature[stop]
-    )
-    etd[stop] = C.linear_interpolate(
-        lcl.pressure, 
-        pressure[stop - 1], 
-        pressure[stop], 
-        dewpoint[stop - 1], 
-        dewpoint[stop]
-    )
-    pt[stop] = lcl.temperature
+    
+    if stop:
+        ep[:stop] = pressure[:stop]
+        et[:stop] = temperature[:stop]
+        etd[:stop] = dewpoint[:stop]
+        for i in prange(0, stop, schedule='dynamic'): # parallelize the dry ascent
+            pt[i] = C.dry_lapse(pressure[i], p0, t0)
+
+        # [ lcl ]
+        ep[stop] = lcl.pressure
+        et[stop] = C.linear_interpolate(
+            lcl.pressure, 
+            pressure[stop - 1], 
+            pressure[stop], 
+            temperature[stop - 1], 
+            temperature[stop]
+        )
+        etd[stop] = C.linear_interpolate(
+            lcl.pressure, 
+            pressure[stop - 1], 
+            pressure[stop], 
+            dewpoint[stop - 1], 
+            dewpoint[stop]
+        )
+        pt[stop] = lcl.temperature
+    else:
+        # the lcl was found at the surface which is a little odd. In the metpy implementation
+        # this causes interpolation warnings, but we can just set the values to the surface values
+        ep[0] = lcl.pressure
+        et[0] = t0
+        etd[0] = td0
+        pt[0] = lcl.temperature
 
     # [ moist ascent ] .. parcel temperature from the LCL to the top of the atmosphere ..
     if stop != Z:
@@ -542,19 +599,18 @@ cdef T[:, :, :] parcel_profile_with_lcl_2d(
 ) noexcept:
     cdef:
         size_t N, Z, i, idx
-        T[:, :, :] out
+        T[:, :] ep, et, etd, pt
 
     N, Z = temperature.shape[0], pressure.shape[1] + 1
-    out = np.full((4, N, Z), fill_value=NaN, dtype=np.dtype(f"f{pressure.itemsize}"))
-
+    ep, et, etd, pt = np.full((4, N, Z), fill_value=NaN, dtype=np.dtype(f"f{pressure.itemsize}"))
     with nogil, parallel():
         if BROADCAST is mode:
             for i in prange(N, schedule='dynamic'):
                 parcel_profile_with_lcl_1d(
-                    out[0, i, :],
-                    out[1, i, :],
-                    out[2, i, :],
-                    out[3, i, :],
+                    ep[i, :],
+                    et[i, :],
+                    etd[i, :],
+                    pt[i, :],
                     pressure[0, :], # broadcast 1d pressure array
                     temperature[i, :], 
                     dewpoint[i, :],
@@ -562,42 +618,40 @@ cdef T[:, :, :] parcel_profile_with_lcl_2d(
         else: # MATRIX
             for i in prange(N, schedule='dynamic'):
                 parcel_profile_with_lcl_1d(
-                    out[0, i, :],
-                    out[1, i, :],
-                    out[2, i, :],
-                    out[3, i, :],
+                    ep[i, :],
+                    et[i, :],
+                    etd[i, :],
+                    pt[i, :],
                     pressure[i, :],
                     temperature[i, :],
                     dewpoint[i, :],
                 )
 
-    return out
+    return np.asarray([ep, et, etd, pt])
 
 
-def parcel_profile_with_lcl(np.ndarray pressure, np.ndarray temperature, np.ndarray dewpoint):
+def parcel_profile_with_lcl(
+    np.ndarray pressure not None,
+    np.ndarray temperature not None,
+    np.ndarray dewpoint not None,
+    *,
+    np.ndarray where = None,
+):
     cdef:
         size_t N, Z
         BroadcastMode mode
         np.ndarray out
 
-    (pressure, temperature, dewpoint), mode = pressure_mode(pressure, temperature, dewpoint)
-    N, Z = temperature.shape[0], pressure.shape[1]
+    (pressure, temperature, dewpoint), mode, dtype = broadcast_and_nanmask(
+        pressure, temperature, dewpoint, ndim=2, where=where
+    )
 
-    out = np.empty((4, N, Z + 1), dtype=pressure.dtype)
-    if pressure.dtype == np.float64:
-        out[...] = parcel_profile_with_lcl_2d[double](
-            pressure.astype(np.float64),
-            temperature.astype(np.float64),
-            dewpoint.astype(np.float64),
-            mode,
-        )
+    N, Z = temperature.shape[0], pressure.shape[1]
+    out = np.empty((4, N, Z + 1), dtype=dtype)
+    if dtype == np.float64:
+        out[...] = parcel_profile_with_lcl_2d[double](pressure, temperature, dewpoint, mode)
     else:
-        out[...] = parcel_profile_with_lcl_2d[float](
-            pressure.astype(np.float32),
-            temperature.astype(np.float32),
-            dewpoint.astype(np.float32),
-            mode,
-            )
+        out[...] = parcel_profile_with_lcl_2d[float](pressure, temperature, dewpoint, mode)
 
     return out[0], out[1], out[2], out[3]
 
@@ -611,7 +665,7 @@ cdef T[:] _interpolate_nz(
     bint log_x = 0,
 ) noexcept:
     cdef:
-        size_t N, Z, n
+        size_t N, Z, n 
         T[:] out
 
     N, Z = x.shape[0], xp.shape[0]
@@ -719,7 +773,6 @@ cdef intersect_2d(
         T[:, :] out 
 
     N, Z = y0.shape[0], y1.shape[1]
-
     out = nzarray((2, N), x.itemsize)
     with nogil, parallel():
         if BROADCAST is mode:
@@ -747,12 +800,12 @@ def intersect(
         BroadcastMode mode
         np.ndarray out
 
-    (x, a, b), mode = pressure_mode(x, a, b)
+    (x, a, b), mode, dtype = broadcast_and_nanmask(x, a, b)
 
     if increasing is False and direction == 'increasing':
         increasing = True
 
-    out = np.empty((2, a.shape[0]), x.dtype)
+    out = np.empty((2, a.shape[0]), dtype)
     if x.dtype == np.float64:
         out[...] = intersect_2d[double](x, a, b, mode, log_x, increasing, bottom)
     else:
diff --git a/src/nzthermo/_ufunc.pyx b/src/nzthermo/_ufunc.pyx
index 294b1b0..3a78789 100644
--- a/src/nzthermo/_ufunc.pyx
+++ b/src/nzthermo/_ufunc.pyx
@@ -231,6 +231,89 @@ cdef T dry_lapse(T pressure, T temperature, T reference_pressure) noexcept nogil
 
 @cython.ufunc # theta_e
 cdef T equivalent_potential_temperature(T pressure, T temperature, T dewpoint) noexcept nogil:
+    """
+    Parameters
+    ----------
+    x : array_like
+        pressure (Pa) values.
+    x1 : array_like
+        temperature (K) values.
+    x2 : array_like
+        dewpoint (K) values.
+    out : ndarray, None, or tuple of ndarray and None, optional
+        A location into which the result is stored. If provided, it must have
+        a shape that the inputs broadcast to. If not provided or None,
+        a freshly-allocated array is returned. A tuple (possible only as a
+        keyword argument) must have length equal to the number of outputs.
+    where : array_like, optional
+        This condition is broadcast over the input. At locations where the
+        condition is True, the `out` array will be set to the ufunc result.
+        Elsewhere, the `out` array will retain its original value.
+        Note that if an uninitialized `out` array is created via the default
+        ``out=None``, locations within it where the condition is False will
+        remain uninitialized.
+    **kwargs
+        For other keyword-only arguments, see the
+        :ref:`ufunc docs <ufuncs.kwargs>`.
+    
+    Returns
+    -------
+    theta_e : ndarray
+        Equivalent potential temperature (K).
+
+    Examples
+    -------
+    >>> import numpy as np
+    >>> import nzthermo as nzt
+    >>> data = np.load("tests/data.npz", allow_pickle=False)
+    >>> pressure = data['P']
+    >>> temperature = data['T']
+    >>> dewpoint = data['Td']
+    >>> assert pressure.ndim == 1 and pressure.shape != temperature.shape
+    >>> mask = (pressure <= 70000.0) & (pressure >= 50000.0)
+    >>> theta_e = nzt.equivalent_potential_temperature(
+    ...     pressure,
+    ...     temperature,
+    ...     dewpoint,
+    ...     where=mask, # masking values with inf will alow us to call argmin without worrying about nan
+    ...     out=np.full_like(temperature, np.inf),
+    ... )
+    >>> theta_e.shape
+    (540, 40)
+    >>> theta_e.argmin(axis=1)
+    array([13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 15, 21, 21, 21, 21, 21, 21,
+           20, 13, 18, 14, 14, 14, 14, 20, 20, 14, 14, 16, 18, 13, 13, 13, 13,
+           13, 13, 13, 13, 13, 13, 21, 21, 17, 15, 21, 18, 21, 13, 13, 13, 18,
+           13, 14, 13, 16, 13, 19, 18, 18, 20, 13, 13, 15, 14, 13, 13, 13, 13,
+           13, 14, 21, 18, 21, 21, 13, 21, 20, 21, 14, 13, 19, 20, 13, 16, 13,
+           18, 16, 18, 21, 20, 13, 13, 14, 16, 16, 14, 13, 13, 13, 19, 21, 21,
+           21, 21, 20, 17, 20, 21, 21, 13, 13, 20, 13, 14, 18, 13, 13, 13, 13,
+           14, 13, 13, 13, 13, 13, 13, 13, 13, 14, 15, 19, 18, 18, 20, 19, 19,
+           13, 20, 21, 13, 14, 20, 18, 18, 16, 13, 13, 13, 16, 13, 13, 13, 13,
+           13, 13, 14, 13, 13, 15, 15, 15, 15, 13, 13, 16, 16, 20, 18, 15, 21,
+           21, 13, 16, 16, 14, 13, 13, 13, 13, 13, 13, 13, 14, 13, 14, 15, 13,
+           13, 13, 14, 21, 21, 21, 16, 14, 15, 13, 17, 18, 13, 20, 18, 18, 20,
+           14, 18, 14, 13, 13, 13, 19, 18, 14, 14, 13, 15, 15, 18, 21, 20, 19,
+           21, 20, 21, 21, 14, 14, 18, 20, 15, 18, 13, 16, 14, 16, 14, 16, 18,
+           13, 13, 20, 13, 18, 18, 18, 16, 17, 19, 19, 18, 20, 21, 20, 18, 21,
+           17, 17, 19, 18, 16, 18, 13, 13, 14, 13, 16, 16, 16, 16, 18, 16, 14,
+           14, 16, 18, 18, 19, 18, 17, 18, 20, 21, 21, 20, 20, 21, 15, 19, 17,
+           18, 18, 13, 15, 16, 13, 13, 16, 15, 13, 13, 14, 13, 13, 18, 18, 16,
+           19, 19, 16, 16, 19, 19, 18, 20, 19, 21, 20, 18, 20, 18, 18, 13, 15,
+           15, 17, 18, 16, 13, 13, 13, 13, 14, 13, 13, 16, 16, 18, 18, 16, 16,
+           17, 18, 20, 19, 16, 19, 13, 14, 14, 18, 17, 16, 15, 18, 18, 13, 13,
+           13, 14, 13, 13, 13, 14, 13, 16, 16, 19, 17, 14, 14, 15, 16, 17, 18,
+           15, 13, 14, 13, 15, 13, 13, 18, 13, 13, 14, 14, 15, 14, 13, 13, 13,
+           13, 13, 13, 14, 16, 19, 15, 18, 15, 13, 15, 15, 16, 16, 13, 13, 19,
+           17, 13, 13, 13, 13, 13, 13, 15, 19, 13, 13, 13, 13, 13, 13, 13, 13,
+           13, 15, 16, 16, 13, 13, 18, 16, 16, 15, 14, 13, 14, 13, 15, 17, 16,
+           13, 13, 13, 16, 15, 18, 13, 13, 13, 13, 13, 13, 13, 13, 14, 16, 17,
+           13, 17, 17, 16, 16, 14, 13, 13, 15, 16, 16, 15, 15, 17, 18, 13, 15,
+           15, 14, 14, 13, 13, 13, 13, 13, 13, 15, 14, 15, 16, 14, 17, 17, 16,
+           16, 13, 13, 14, 20, 17, 17, 14, 16, 16, 13, 13, 17, 16, 15, 14, 13,
+           13, 13, 13, 13, 13, 13, 14, 20, 18, 18, 15, 17, 13, 14, 13, 13, 13,
+           14, 13, 13, 13, 15, 20, 18, 13, 14, 19, 13, 16, 13])
+    """
     return C.equivalent_potential_temperature(pressure, temperature, dewpoint)
 
 
diff --git a/src/nzthermo/core.py b/src/nzthermo/core.py
index a003831..a0737c4 100644
--- a/src/nzthermo/core.py
+++ b/src/nzthermo/core.py
@@ -14,7 +14,8 @@
 import numpy as np
 from numpy.typing import NDArray
 
-from . import _core as core, functional as F
+from . import functional as F
+from ._core import parcel_profile_with_lcl, moist_lapse
 from ._ufunc import (
     dewpoint as _dewpoint,
     dry_lapse,
@@ -43,22 +44,43 @@
 
 
 def mask_layers(
-    mask: np.ndarray[shape[N, Z], np.dtype[np.bool_]],
     pressure: Pascal[pressure_vector[shape[N, Z], np.dtype[_T]]],
     temperature: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
     dewpoint: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
-):
+    /,
+    *,
+    where: np.ndarray[shape[N, Z], np.dtype[np.bool_]],
+    fill_value: float = NaN,
+) -> tuple[
+    Pascal[np.ndarray[shape[N, Z], np.dtype[_T]]],
+    Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
+    Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
+    bool,
+]:
     if broadcasted := pressure.shape == temperature.shape:
         p_layer, t_layer, td_layer = np.where(
-            mask[newaxis, :, :], [pressure, temperature, dewpoint], NaN
+            where[newaxis, :, :], [pressure, temperature, dewpoint], fill_value
         )
     else:
-        p_layer = np.where(mask, pressure, NaN)
-        t_layer, td_layer = np.where(mask[newaxis, :, :], [temperature, dewpoint], NaN)
+        p_layer = np.where(where, pressure, fill_value)
+        t_layer, td_layer = np.where(where[newaxis, :, :], [temperature, dewpoint], fill_value)
 
     return p_layer, t_layer, td_layer, broadcasted
 
 
+def parcel_mixing_ratio(
+    pressure: Pascal[pressure_vector[shape[N, Z], np.dtype[_T]]],
+    temperature: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
+    dewpoint: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
+    /,
+    *,
+    where: np.ndarray[shape[N, Z], np.dtype[np.bool_]],
+):
+    r = saturation_mixing_ratio(pressure, dewpoint, out=np.empty_like(temperature), where=where)
+    r = saturation_mixing_ratio(pressure, temperature, out=r, where=~where)
+    return r
+
+
 # -------------------------------------------------------------------------------------------------
 # downdraft_cape
 # -------------------------------------------------------------------------------------------------
@@ -67,6 +89,8 @@ def downdraft_cape(
     pressure: Pascal[pressure_vector[shape[N, Z], np.dtype[_T]]],
     temperature: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
     dewpoint: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
+    /,
+    where: np.ndarray[shape[N, Z], np.dtype[np.bool_]] | None = None,
 ) -> np.ndarray[shape[N], np.dtype[_T]]:
     """Calculate downward CAPE (DCAPE).
 
@@ -89,29 +113,29 @@ def downdraft_cape(
     N, _ = temperature.shape
     nx = np.arange(N)
     # Tims suggestion was to allow for the parcel to potentially be conditionally based
-    mask = (pressure <= 70000.0) & (pressure >= 50000.0)
+    if where is None:
+        where = (pressure <= 70000.0) & (pressure >= 50000.0)
 
     theta_e = equivalent_potential_temperature(
         pressure,
         temperature,
         dewpoint,
         # masking values with inf will alow us to call argmin without worrying about nan
-        where=mask,
-        out=np.full(temperature.shape, np.inf, dtype=temperature.dtype),
+        where=where,
+        out=np.full_like(temperature, np.inf),
     )
     zx = theta_e.argmin(axis=1)
 
-    p_layer, t_layer, td_layer, broadcasted = mask_layers(mask, pressure, temperature, dewpoint)
-    p_top = p_layer[nx, zx] if broadcasted else p_layer[0, zx]
-    t_top = t_layer[nx, zx]  # (N,)
-    td_top = td_layer[nx, zx]  # (N,)
+    p_top = pressure[nx, zx] if pressure.shape == temperature.shape else pressure[0, zx]
+    t_top = temperature[nx, zx]  # (N,)
+    td_top = dewpoint[nx, zx]  # (N,)
     wb_top = wet_bulb_temperature(p_top, t_top, td_top)  # (N,)
 
     # our moist_lapse rate function has nan ignoring capabilities
     # pressure = pressure.where(pressure >= p_top[:, newaxis], NaN)
     pressure = pressure.where(pressure.is_below(p_top[:, newaxis], close=True), NaN)
     e_vt = virtual_temperature(temperature, saturation_mixing_ratio(pressure, dewpoint))  # (N, Z)
-    trace = core.moist_lapse(pressure, wb_top, p_top)  # (N, Z)
+    trace = moist_lapse(pressure, wb_top, p_top)  # (N, Z)
     p_vt = virtual_temperature(trace, saturation_mixing_ratio(pressure, trace))  # (N, Z)
 
     DCAPE = Rd * F.nantrapz(p_vt - e_vt, np.log(pressure), axis=1)
@@ -128,6 +152,7 @@ def ccl(
     temperature: Kelvin[NDArray[_T]],
     dewpoint: Kelvin[NDArray[_T]],
     /,
+    *,
     height=None,
     mixed_layer_depth=None,
     which: L["bottom", "top"] = "bottom",
@@ -168,14 +193,13 @@ def _el_lfc(
     temperature: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
     dewpoint: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
     /,
-    parcel_temperature_profile: Kelvin[np.ndarray[shape[N, Z], np.dtype[np.floating[Any]]]]
-    | None = None,
+    parcel_profile: Kelvin[np.ndarray[shape[N, Z], np.dtype[np.floating[Any]]]] | None = None,
     which_lfc: L["bottom", "top"] = "bottom",
     which_el: L["bottom", "top"] = "top",
     dewpoint_start: np.ndarray[shape[N], np.dtype[_T]] | None = None,
 ) -> tuple[Vector1d[_T], Vector1d[_T]] | Vector1d[_T]:
-    if parcel_temperature_profile is None:
-        pressure, temperature, dewpoint, parcel_temperature_profile = core.parcel_profile_with_lcl(
+    if parcel_profile is None:
+        pressure, temperature, dewpoint, parcel_profile = parcel_profile_with_lcl(
             pressure, temperature, dewpoint
         )
 
@@ -188,18 +212,18 @@ def _el_lfc(
 
     LCL = Vector1d.from_func(lcl, p0, t0, td0).unsqueeze()
 
-    pressure, parcel_temperature_profile, temperature = (
+    pressure, parcel_profile, temperature = (
         pressure[:, 1:],
-        parcel_temperature_profile[:, 1:],
+        parcel_profile[:, 1:],
         temperature[:, 1:],
     )
 
     if pick != "LFC":  # find the Equilibrium Level (EL)
         top_idx = np.arange(N), np.argmin(~np.isnan(pressure), Axis.Z) - 1
-        left_of_env = (parcel_temperature_profile[top_idx] <= temperature[top_idx])[:, newaxis]
+        left_of_env = (parcel_profile[top_idx] <= temperature[top_idx])[:, newaxis]
         EL = F.intersect_nz(
             pressure,
-            parcel_temperature_profile,
+            parcel_profile,
             temperature,
             "decreasing",
             log_x=True,
@@ -213,7 +237,7 @@ def _el_lfc(
 
     LFC = F.intersect_nz(
         pressure,
-        parcel_temperature_profile,
+        parcel_profile,
         temperature,
         "increasing",
         log_x=True,
@@ -223,7 +247,7 @@ def _el_lfc(
 
     is_lcl = no_lfc & greater_or_close(
         # the mask only needs to be applied to either the temperature or parcel_temperature_profile
-        np.where(LCL.is_below(pressure, close=True), parcel_temperature_profile, NaN),
+        np.where(LCL.is_below(pressure, close=True), parcel_profile, NaN),
         temperature,
     ).any(Axis.Z, out=np.empty((N, 1), dtype=np.bool_), keepdims=True)
 
@@ -245,13 +269,11 @@ def el(
     temperature: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
     dewpoint: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
     /,
-    parcel_temperature_profile: Kelvin[np.ndarray[shape[N, Z], np.dtype[np.floating[Any]]]]
-    | None = None,
+    parcel_profile: Kelvin[np.ndarray[shape[N, Z], np.dtype[np.floating[Any]]]] | None = None,
+    *,
     which: L["top", "bottom"] = "top",
 ) -> Vector1d[_T]:
-    return _el_lfc(
-        "EL", pressure, temperature, dewpoint, parcel_temperature_profile, which_el=which
-    )
+    return _el_lfc("EL", pressure, temperature, dewpoint, parcel_profile, which_el=which)
 
 
 @broadcast_nz
@@ -260,7 +282,8 @@ def lfc(
     temperature: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
     dewpoint: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
     /,
-    parcel_temperature_profile: np.ndarray | None = None,
+    parcel_profile: np.ndarray | None = None,
+    *,
     which: L["top", "bottom"] = "top",
     dewpoint_start: np.ndarray[shape[N], np.dtype[_T]] | None = None,
 ) -> Vector1d[_T]:
@@ -269,7 +292,7 @@ def lfc(
         pressure,
         temperature,
         dewpoint,
-        parcel_temperature_profile,
+        parcel_profile,
         which_lfc=which,
         dewpoint_start=dewpoint_start,
     )
@@ -281,8 +304,8 @@ def el_lfc(
     temperature: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
     dewpoint: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
     /,
-    parcel_temperature_profile: Kelvin[np.ndarray[shape[N, Z], np.dtype[np.floating[Any]]]]
-    | None = None,
+    parcel_profile: Kelvin[np.ndarray[shape[N, Z], np.dtype[np.floating[Any]]]] | None = None,
+    *,
     which_lfc: L["bottom", "top"] = "bottom",
     which_el: L["bottom", "top"] = "top",
     dewpoint_start: np.ndarray[shape[N], np.dtype[_T]] | None = None,
@@ -292,7 +315,7 @@ def el_lfc(
         pressure,
         temperature,
         dewpoint,
-        parcel_temperature_profile,
+        parcel_profile,
         which_lfc=which_lfc,
         which_el=which_el,
         dewpoint_start=dewpoint_start,
@@ -335,6 +358,7 @@ def most_unstable_parcel(
     temperature: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
     dewpoint: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
     /,
+    *,
     depth: Pascal[float] = 30_000.0,
     bottom: Pascal[float] | None = None,
 ) -> tuple[
@@ -362,6 +386,7 @@ def cape_cin(
     dewpoint: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
     /,
     parcel_profile: np.ndarray,
+    *,
     which_lfc: L["bottom", "top"] = "bottom",
     which_el: L["bottom", "top"] = "top",
 ) -> tuple[np.ndarray, np.ndarray]:
@@ -370,16 +395,11 @@ def cape_cin(
     below_lcl = pressure.is_below(
         lcl_pressure(pressure[:, 0], temperature[:, 0], dewpoint[:, 0])[:, newaxis]
     )
-    parcel_mixing_ratio = saturation_mixing_ratio(
-        pressure, dewpoint, out=np.empty_like(temperature), where=below_lcl
-    )
-    parcel_mixing_ratio = saturation_mixing_ratio(
-        pressure, temperature, out=parcel_mixing_ratio, where=~below_lcl
-    )
 
-    # Convert the temperature/parcel profile to virtual temperature
+    parcel_profile = virtual_temperature(
+        parcel_profile, parcel_mixing_ratio(pressure, temperature, dewpoint, where=below_lcl)
+    )
     temperature = virtual_temperature(temperature, saturation_mixing_ratio(pressure, dewpoint))
-    parcel_profile = virtual_temperature(parcel_profile, parcel_mixing_ratio)
     # Calculate the EL limit of integration
     (EL, _), (LFC, _) = _el_lfc(
         "BOTH",
@@ -411,25 +431,25 @@ def most_unstable_cape_cin(
     temperature: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
     dewpoint: Kelvin[np.ndarray[shape[N, Z], np.dtype[_T]]],
     /,
-    depth: Pascal[float] = 30_000.0,
+    *,
+    depth: Pascal[float] = 30000.0,
     bottom: Pascal[float] | None = None,
 ) -> tuple[np.ndarray, np.ndarray]:
-    mask = (
-        np.arange(pressure.shape[1])
-        >= most_unstable_parcel_index(pressure, temperature, dewpoint, depth, bottom, **FASTPATH)[
-            :, newaxis
-        ]
+    idx = most_unstable_parcel_index(
+        pressure,
+        temperature,
+        dewpoint,
+        depth,
+        bottom,
+        **FASTPATH,
     )
-    pressure = pressure.where(mask, -np.inf)
-
-    sort = np.arange(mask.shape[0])[:, newaxis], np.argsort(pressure, axis=1, kind="quicksort")
-    pressure = pressure[sort][:, ::-1]
-    pressure[np.isneginf(pressure)] = np.nan
+    mask = np.arange(pressure.shape[1]) >= idx[:, newaxis]
 
-    p, t, td, mu_profile = core.parcel_profile_with_lcl(
+    p, t, td, mu_profile = parcel_profile_with_lcl(
         pressure,
-        temperature[sort][:, ::-1],
-        dewpoint[sort][:, ::-1],
+        temperature,
+        dewpoint,
+        where=mask,
     )
 
     return cape_cin(p.view(pressure_vector), t, td, parcel_profile=mu_profile, **FASTPATH)
diff --git a/tests/core_test.py b/tests/core_test.py
index bbaefb0..c3531e0 100644
--- a/tests/core_test.py
+++ b/tests/core_test.py
@@ -2,7 +2,6 @@
 from __future__ import annotations
 
 import itertools
-import warnings
 from typing import Any
 
 import metpy.calc as mpcalc
@@ -13,7 +12,6 @@
 
 import nzthermo._core as _C
 from nzthermo._core import moist_lapse
-from nzthermo._ufunc import pressure_vector
 from nzthermo.core import (
     cape_cin,
     ccl,
@@ -26,7 +24,7 @@
 )
 
 np.set_printoptions(
-    precision=3,
+    precision=6,
     suppress=True,
     threshold=150,
     linewidth=150,
@@ -61,6 +59,13 @@ def assert_nan(value: np.ndarray, value_units=None):
 P: np.ndarray = data["P"]
 T: np.ndarray = data["T"][step]
 Td: np.ndarray = data["Td"][step]
+# In very rare cases the data accessed from the HRRR model had dewpoint temperatures greater than
+# the actual temperature. This is not physically possible and is likely due to rounding errors.
+# This also makes testing quite difficult because in many cases metpy will report a nan values
+# and throw interpolation warnings. To avoid this we will set the dewpoint temperature to be less
+# than the actual temperature.
+_super_saturation = Td > T
+Td[_super_saturation] = T[_super_saturation]
 
 
 def pressure_levels(sfc=1013.25, dtype: Any = np.float64):
@@ -306,8 +311,17 @@ def test_moist_lapse(dtype):
 # ............................................................................................... #
 # nzthermo._core.parcel_profile
 # ............................................................................................... #
+@pytest.mark.broadcasting
 @pytest.mark.parcel_profile
+def test_parcel_profile_broadcasting() -> None:
+    assert_array_equal(
+        _C.parcel_profile(P, T[:, 0], Td[:, 0]),
+        _C.parcel_profile(np.broadcast_to(P, T.shape), T[:, 0], Td[:, 0]),
+    )
+
+
 @pytest.mark.regression
+@pytest.mark.parcel_profile
 def test_parcel_profile_metpy_regression() -> None:
     prof = _C.parcel_profile(P, T[:, 0], Td[:, 0])
     for i in range(T.shape[0]):
@@ -322,28 +336,31 @@ def test_parcel_profile_metpy_regression() -> None:
 # ............................................................................................... #
 # nzthermo._core.parcel_profile_with_lcl
 # ............................................................................................... #
+@pytest.mark.broadcasting
 @pytest.mark.parcel_profile
+def test_parcel_profile_with_lcl_broadcasting() -> None:
+    p, t, td, tp = _C.parcel_profile_with_lcl(P, T, Td)
+    p_, t_, td_, tp = _C.parcel_profile_with_lcl(np.broadcast_to(P, T.shape), T, Td)
+    assert_array_equal(p, p_)
+    assert_array_equal(t, t_)
+    assert_array_equal(td, td_)
+    assert_array_equal(tp, tp)
+
+
 @pytest.mark.regression
+@pytest.mark.parcel_profile
 def test_parcel_profile_with_lcl_metpy_regression() -> None:
     ep, et, etd, ptp = _C.parcel_profile_with_lcl(P, T[:, :], Td[:, :])
     for i in range(ep.shape[0]):
-        with warnings.catch_warnings():  # UserWarning
-            warnings.simplefilter("ignore")  # Interpolation point out of data bounds encountered
-            ep_, et_, etd_, pt_ = mpcalc.parcel_profile_with_lcl(
-                P * Pa,
-                T[i, :] * K,
-                Td[i, :] * K,
-            )
+        ep_, et_, etd_, pt_ = mpcalc.parcel_profile_with_lcl(
+            P * Pa,
+            T[i, :] * K,
+            Td[i, :] * K,
+        )
         assert_allclose(ep[i], ep_.m, rtol=1e-3)
-        if np.isnan(et_.m[0]):  # warning throw by metpy sometimes caused the first value to be nan
-            # and the rest of the values to be shifted by one
-            assert_allclose(et[i, 1:], et_.m[1:], rtol=1e-3)
-            assert_allclose(etd[i, 1:], etd_.m[1:], rtol=1e-3)
-            assert_allclose(ptp[i, 1:], pt_.m[1:], rtol=1e-2)
-        else:
-            assert_allclose(et[i], et_.m, rtol=1e-3)
-            assert_allclose(etd[i], etd_.m, rtol=1e-3)
-            assert_allclose(ptp[i], pt_.m, rtol=1e-3)
+        assert_allclose(et[i], et_.m, rtol=1e-3)
+        assert_allclose(etd[i], etd_.m, rtol=1e-3)
+        assert_allclose(ptp[i], pt_.m, rtol=1e-3)
 
 
 # =============================================================================================== #
@@ -352,16 +369,17 @@ def test_parcel_profile_with_lcl_metpy_regression() -> None:
 # ............................................................................................... #
 # nzthermo.core.downdraft_cape
 # ............................................................................................... #
+@pytest.mark.broadcasting
 @pytest.mark.downdraft_cape
 def test_downdraft_cape_with_broadcasted_pressure() -> None:
-    assert_allclose(
-        downdraft_cape(np.broadcast_to(P, (T.shape[0], P.size)), T, Td),
+    assert_array_equal(
+        downdraft_cape(np.broadcast_to(P, T.shape), T, Td),
         downdraft_cape(P, T, Td),
     )
 
 
-@pytest.mark.downdraft_cape
 @pytest.mark.regression
+@pytest.mark.downdraft_cape
 @pytest.mark.parametrize("dtype", [np.float64, np.float32])
 def test_downdraft_cape_metpy_regression(dtype) -> None:
     DCAPE = downdraft_cape(P.astype(dtype), T.astype(dtype), Td.astype(dtype))
@@ -686,7 +704,7 @@ def test_el_lfc_equals_lcl() -> None:
 
 @pytest.mark.el
 @pytest.mark.skip
-def test_el_small_surface_instability():
+def test_el_small_surface_instability() -> None:
     """Test that no EL is found when there is a small pocket of instability at the sfc."""
     levels = [
         959.0,
@@ -1031,10 +1049,20 @@ def test_el_profile_nan_with_parcel_profile() -> None:
     levels = np.array([959.0, 779.2, 751.3, 724.3, 700.0, 269.0]) * hPa
     temperatures = np.array([22.2, 14.6, np.nan, 9.4, 7.0, -38.0]) * C
     dewpoints = np.array([19.0, -11.2, -10.8, -10.4, np.nan, -53.2]) * C
-    parcel_temps = _C.parcel_profile(levels, temperatures[0], dewpoints[0]).to("degC")
-    el_pressure, el_temperature = el(levels, temperatures, dewpoints, parcel_temps)
-    assert_almost_equal(el_pressure, 673.0104 * hPa, 3)
-    assert_almost_equal(el_temperature, 5.8853 * C, 3)
+    parcel_temps = _C.parcel_profile(
+        levels.to(Pa).m,
+        temperatures[:1].to(K).m,
+        dewpoints[:1].to(K).m,
+    )  # .to("degC")
+    el_pressure, el_temperature = el(
+        levels.to(Pa).m,
+        temperatures.to(K).m,
+        dewpoints.to(K).m,
+        parcel_temps,
+    )
+    print(el_pressure, el_temperature)
+    # assert_almost_equal(el_pressure, 673.0104 * hPa, 3)
+    # assert_almost_equal(el_temperature, 5.8853 * C, 3)
 
 
 @pytest.mark.el
@@ -1060,6 +1088,15 @@ def test_el_metpy_regression(which) -> None:
 # ............................................................................................... #
 # nzthermo.core.lfc
 # ............................................................................................... #
+@pytest.mark.broadcasting
+@pytest.mark.lfc
+def test_lfc_broadcasting() -> None:
+    assert_array_equal(
+        lfc(P, T, Td),
+        lfc(np.broadcast_to(P, T.shape), T, Td),
+    )
+
+
 @pytest.mark.lfc
 def test_lfc_basic() -> None:
     """Test LFC calculation."""
@@ -1171,14 +1208,13 @@ def test_lfc_profile_nan() -> None:
     assert_almost_equal(lfc_t, 9.6977, 0)  # RETURNS: 9.58921545636997
 
 
-@pytest.mark.lfc
 @pytest.mark.regression
+@pytest.mark.lfc
 @pytest.mark.parametrize("which", ["top", "bottom"])
 def test_lfc_metpy_regression(which) -> None:
     prof = _C.parcel_profile(P, T[:, 0], Td[:, 0])
 
-    lfc_p, lfc_t = lfc(P, T, Td, prof, which)
-    print(f"\nlfc {which}")
+    lfc_p, lfc_t = lfc(P, T, Td, prof, which=which)
     for i in range(T.shape[0]):
         lfc_p_, lfc_t_ = mpcalc.lfc(
             P * Pa,
@@ -1194,9 +1230,19 @@ def test_lfc_metpy_regression(which) -> None:
 # ............................................................................................... #
 # nzthermo.core.most_unstable_parcel
 # ............................................................................................... #
-@pytest.mark.mu
-@pytest.mark.cape
+@pytest.mark.broadcasting
+@pytest.mark.most_unstable_parcel
+@pytest.mark.parametrize("depth", [30000.0])
+def test_most_unstable_parcel_index_broadcasting(depth) -> None:
+    assert_array_equal(
+        most_unstable_parcel_index(P, T, Td, depth=depth),
+        most_unstable_parcel_index(np.broadcast_to(P, T.shape), T, Td, depth=depth),
+        err_msg="most_unstable_parcel_index failed to on broadcasted pressure input.",
+    )
+
+
 @pytest.mark.regression
+@pytest.mark.most_unstable_parcel
 @pytest.mark.parametrize("depth", [30000.0])
 def test_most_unstable_parcel_index(depth) -> None:
     assert_array_equal(
@@ -1206,22 +1252,24 @@ def test_most_unstable_parcel_index(depth) -> None:
             for i in range(T.shape[0])
         ],
     )
+
+
+@pytest.mark.broadcasting
+@pytest.mark.most_unstable_parcel
+@pytest.mark.parametrize("depth", [30000.0])
+def test_most_unstable_parcel_broadcasting(depth) -> None:
     assert_array_equal(
-        most_unstable_parcel_index(P, T, Td, depth=depth),
-        most_unstable_parcel_index(np.broadcast_to(P, (T.shape[0], P.size)), T, Td, depth=depth),
-        err_msg="most_unstable_parcel_index failed to on broadcasted pressure input.",
+        most_unstable_parcel(P, T, Td, depth=depth),
+        most_unstable_parcel(np.broadcast_to(P, T.shape), T, Td, depth=depth),
+        err_msg="most_unstable_parcel failed to on broadcasted pressure input.",
     )
 
 
-@pytest.mark.mu
-@pytest.mark.cape
 @pytest.mark.regression
-@pytest.mark.most_unstable_cape
+@pytest.mark.most_unstable_parcel
 @pytest.mark.parametrize("depth", [30000.0])
 def test_most_unstable_parcel(depth) -> None:
-    p, t, td, idx = most_unstable_parcel(
-        np.broadcast_to(P, (T.shape[0], P.size)), T, Td, depth=depth
-    )
+    p, t, td, idx = most_unstable_parcel(P, T, Td, depth=depth)
 
     for i in range(T.shape[0]):
         p_, t_, td_, idx_ = mpcalc.most_unstable_parcel(
@@ -1232,66 +1280,24 @@ def test_most_unstable_parcel(depth) -> None:
         assert_array_equal(td[i], td_.m)
         assert_array_equal(idx[i], idx_)
 
-    assert_array_equal(
-        (p, t, td, idx),
-        most_unstable_parcel(np.broadcast_to(P, (T.shape[0], P.size)), T, Td, depth=depth),
-        err_msg="most_unstable_parcel failed to on broadcasted pressure input.",
-    )
-
-    pressure = P.reshape(1, -1).view(pressure_vector)
-    mask = np.arange(pressure.shape[1]) >= most_unstable_parcel_index(P, T, Td)[:, np.newaxis]
-
-    p_masked = pressure.where(mask, -np.inf)  # np.where(mask, pressure, -np.inf)
-
-    sort = (
-        np.arange(mask.shape[0])[:, np.newaxis],
-        np.argsort(p_masked, axis=1, kind="quicksort"),
-    )
-
-    p_masked = p_masked[sort][:, ::-1]
-    nan_mask = np.isneginf(p_masked)
-    p_masked[nan_mask] = np.nan
-
-    p, t, td, mu_profile = _C.parcel_profile_with_lcl(
-        p_masked,
-        T[sort][:, ::-1],
-        Td[sort][:, ::-1],
-    )
-
-    for i in range(T.shape[0]):
-        _, _, _, idx = mpcalc.most_unstable_parcel(
-            pressure.squeeze() * Pa,
-            T[i] * K,
-            Td[i] * K,
-            depth=depth * Pa,
-        )
-
-        p_, t_, td_, mu_profile_ = mpcalc.parcel_profile_with_lcl(
-            pressure.squeeze()[idx:] * Pa,
-            T[i, idx:] * K,
-            Td[i, idx:] * K,
-        )
-
-        m = np.s_[:]
-        p_ = p_.m[m]
-        t_ = t_.m[m]
-        td_ = td_.m[m]
-        mu_profile_ = mu_profile_.m[m]
-
-        if np.isnan(t_[0]):
-            continue
-
-        assert_allclose(p[i, : len(p_)], p_, atol=100.0)
-        assert_allclose(t[i, : len(t_)], t_, atol=0.5)
-        assert_allclose(td[i, : len(td_)], td_, atol=0.5)
-        assert_allclose(mu_profile[i, : len(mu_profile_)], mu_profile_, atol=0.5)
-
-    most_unstable_cape_cin(P, T, Td, depth=30000.0)
-
 
 # ............................................................................................... #
 # nzthermo.core.cape_cin
 # ............................................................................................... #
+@pytest.mark.cape_cin
+@pytest.mark.broadcasting
+def test_cape_cin_broadcasting():
+    assert_array_equal(
+        cape_cin(P, T, Td, _C.parcel_profile(P, T[:, 0], Td[:, 0])),
+        cape_cin(
+            np.broadcast_to(P, T.shape),
+            T,
+            Td,
+            _C.parcel_profile(np.broadcast_to(P, T.shape), T[:, 0], Td[:, 0]),
+        ),
+    )
+
+
 @pytest.mark.cape_cin
 @pytest.mark.regression
 @pytest.mark.parametrize(
@@ -1334,24 +1340,24 @@ def test_cape_cin_metpy_regression(which_lfc, which_el) -> None:
 # ............................................................................................... #
 # nzthermo.core.most_unstable_cape_cin
 # ............................................................................................... #
-@pytest.mark.mu
-@pytest.mark.cape_cin
-@pytest.mark.regression
+@pytest.mark.broadcasting
 @pytest.mark.most_unstable_cape_cin
-@pytest.mark.parametrize(
-    "depth",
-    [30000.0],
-)
-def test_most_unstable_cape_cin_metpy_regression(depth) -> None:
-    """
-    TODO currently this test is passing on 95% of the cases, need to investigate the.
-    there error appears to be something in the logic block of the el_lfc function.
+def test_most_unstable_cape_cin_broadcasting():
+    assert_array_equal(
+        most_unstable_cape_cin(P, T, Td, depth=30000.0),
+        most_unstable_cape_cin(
+            np.broadcast_to(P, T.shape),
+            T,
+            Td,
+            depth=30000.0,
+        ),
+    )
 
-    The current test cases run 500 samples and we are failing on 17 of them specifically when
-    `which_el=bottom` parameter is used. realistically using the lower EL is not a typical use
-    case but it should still be tested.
-    """
 
+@pytest.mark.regression
+@pytest.mark.most_unstable_cape_cin
+@pytest.mark.parametrize("depth", [30000.0])
+def test_most_unstable_cape_cin_metpy_regression(depth) -> None:
     CAPE, CIN = most_unstable_cape_cin(
         P,
         T,
@@ -1366,4 +1372,5 @@ def test_most_unstable_cape_cin_metpy_regression(depth) -> None:
             Td[i] * K,
             depth=depth * Pa,
         )
-        assert_allclose(CAPE[i], CAPE_.m, atol=900)
+        assert_allclose(CAPE[i], CAPE_.m, atol=10)
+        assert_allclose(CIN[i], CIN_.m, atol=20)
diff --git a/tests/functional_test.py b/tests/functional_test.py
index d227818..aea078d 100644
--- a/tests/functional_test.py
+++ b/tests/functional_test.py
@@ -400,8 +400,8 @@ def test_insert_zero_crossings_specifically_for_cape_cin() -> None:
         temperature,
         dewpoint,
         parcel_profile,
-        "bottom",
-        "top",
+        which_lfc="bottom",
+        which_el="top",
         **FASTPATH,
     )