diff --git a/eurec4a_environment/__init__.py b/eurec4a_environment/__init__.py
index c2d5f4d..52f9b57 100644
--- a/eurec4a_environment/__init__.py
+++ b/eurec4a_environment/__init__.py
@@ -37,7 +37,9 @@ def get_fields(ds, *names, **names_and_units):
     dataarrays = []
     coords = {}
 
-    all_fields = itertools.chain(names_and_units.items(), zip(names, [None] * len(names)))
+    all_fields = itertools.chain(
+        names_and_units.items(), zip(names, [None] * len(names))
+    )
 
     for (name, units) in all_fields:
         da = get_field(ds=ds, name=name, units=units)
diff --git a/eurec4a_environment/plot/profile.py b/eurec4a_environment/plot/profile.py
index 5139ba1..baa3b70 100644
--- a/eurec4a_environment/plot/profile.py
+++ b/eurec4a_environment/plot/profile.py
@@ -1,9 +1,6 @@
-import numpy as np
 import matplotlib.pyplot as plt
-import xarray as xr
 import matplotlib.dates as mdates
 import seaborn as sns
-import typhon
 
 sns.set(
     context="notebook",
@@ -93,8 +90,10 @@ def plot_profile_1D(
     height_labels=None,
     **kwargs
 ):
-    """Plot a 1D profile plot of two variables and a height level 
-    input a 1D dataset (that is profiles e.g. of one timestep or sounding or a mean profile)
+    """
+    Plot a 1D profile plot of two variables and a height level input a 1D
+    dataset (that is profiles e.g. of one timestep or sounding or a mean
+    profile)
     """
 
     fig, ax1 = plt.subplots(figsize=(12, 12))
diff --git a/eurec4a_environment/unit.py b/eurec4a_environment/unit.py
index 9afff24..dcee7e8 100644
--- a/eurec4a_environment/unit.py
+++ b/eurec4a_environment/unit.py
@@ -5,8 +5,6 @@
 
     HAS_UDUNITS2 = True
 except FileNotFoundError:
-    import warnings
-
     HAS_UDUNITS2 = False
 
 
diff --git a/eurec4a_environment/variables/boundary_layer/__init__.py b/eurec4a_environment/variables/boundary_layer/__init__.py
index 5ed4712..f5ea189 100644
--- a/eurec4a_environment/variables/boundary_layer/__init__.py
+++ b/eurec4a_environment/variables/boundary_layer/__init__.py
@@ -1,2 +1,2 @@
 from . import lcl  # noqa
-from . import mixed_layer_height
+from . import mixed_layer_height  # noqa
diff --git a/eurec4a_environment/variables/boundary_layer/inversion_height.py b/eurec4a_environment/variables/boundary_layer/inversion_height.py
index 2d87863..cafa0f5 100644
--- a/eurec4a_environment/variables/boundary_layer/inversion_height.py
+++ b/eurec4a_environment/variables/boundary_layer/inversion_height.py
@@ -1,5 +1,3 @@
-import xarray as xr
-
 from ... import get_field
 from ... import nomenclature as nom
 
diff --git a/eurec4a_environment/variables/boundary_layer/lcl.py b/eurec4a_environment/variables/boundary_layer/lcl.py
index ec10d6b..2405b97 100644
--- a/eurec4a_environment/variables/boundary_layer/lcl.py
+++ b/eurec4a_environment/variables/boundary_layer/lcl.py
@@ -7,7 +7,10 @@
 
 
 def find_LCL_Bolton(
-    ds, temperature=nom.TEMPERATURE, rh=nom.RELATIVE_HUMIDITY, altitude=nom.ALTITUDE,
+    ds,
+    temperature=nom.TEMPERATURE,
+    rh=nom.RELATIVE_HUMIDITY,
+    altitude=nom.ALTITUDE,
     vertical_coord=nom.ALTITUDE,
 ):
     """
diff --git a/eurec4a_environment/variables/boundary_layer/mixed_layer_height.py b/eurec4a_environment/variables/boundary_layer/mixed_layer_height.py
index 143e100..aa6e16e 100644
--- a/eurec4a_environment/variables/boundary_layer/mixed_layer_height.py
+++ b/eurec4a_environment/variables/boundary_layer/mixed_layer_height.py
@@ -3,8 +3,10 @@
 from scipy.signal import find_peaks
 import statsmodels.api as sm
 
-from ... import get_field
-from ... import nomenclature as nom
+from ... import get_field, get_fields, nomenclature as nom
+from ..utils import apply_by_column
+
+DEBUG = False
 
 
 def calc_peak_RH(
@@ -33,39 +35,50 @@ def calc_peakRH_linearize(
     ds,
     altitude=nom.ALTITUDE,
     rh=nom.RELATIVE_HUMIDITY,
-    time_dim="sounding",
+    time="sounding",
     z_min=200.0,
     z_max=1500.0,
     z_min_lin=50.0,
     z_max_lin=400.0,
+    vertical_coord=":altitude",
 ):
     """
-     Calculate maximum in relative humidity (RH) profile that minimizes difference between
-     observed RH profile and RH profile linearized around lower levels (i.e. 200-400m).
-     Assume moisture is well-mixed at these lower levels, such that RH would increase and
-     assume linear profile as an idealization
-     Inputs:
-         -- ds: dataset
-         -- altitude, rh, time_dim: variable names
-         -- z_min and z_max: lower and upper bounds for mixed layer height
-         -- z_min_lin and z_max_lin: bounds for linearization of observed RH profile
-     Outputs:
-         -- da: datarray containing h_peakRH_linfit
-     """
-    da_rh = get_field(ds=ds, name=rh, units="%")
-    da_alt = get_field(ds=ds, name=altitude, units="m")
-
-    h_peakRH_linfit = np.zeros(len(da_rh))
-
-    dz = int(da_alt.diff(dim=altitude)[1])  # dz=10m
-
-    mixed_layer = np.logical_and(da_alt >= z_min, da_alt <= z_max)
-    ml_linfit = np.logical_and(
-        da_alt >= z_min_lin, da_alt <= z_max_lin
-    )  # for linearized RH profile
-
-    for i in range(len(da_rh)):
-        rh_profile = da_rh.isel({time_dim: i}).interpolate_na(dim=altitude)
+    Calculate maximum in relative humidity (RH) profile that minimizes
+    difference between observed RH profile and RH profile linearized around
+    lower levels (i.e. 200-400m).  Assume moisture is well-mixed at these lower
+    levels, such that RH would increase and assume linear profile as an
+    idealization.
+
+    Inputs:
+        - ds: dataset
+        - altitude, rh, time: variable names
+        - z_min and z_max: lower and upper bounds for mixed layer height
+        - z_min_lin and z_max_lin: bounds for linearization of observed RH profile
+
+    Outputs:
+        - da: datarray containing h_peakRH_linfit
+    """
+    # if `vertical_coord` hasn't been changed we'll assume that the vertical
+    # coordinate has name of the altitude variable
+    if vertical_coord == ":altitude":
+        vertical_coord = altitude
+
+    import matplotlib.pyplot as plt
+
+    def _calc_mixed_layer_depth(ds_column):
+        da_alt = ds_column[altitude]
+
+        mixed_layer = np.logical_and(da_alt >= z_min, da_alt <= z_max)
+        ml_linfit = np.logical_and(
+            da_alt >= z_min_lin, da_alt <= z_max_lin
+        )  # for linearized RH profile
+
+        dz = int(da_alt.diff(dim=altitude)[1])  # dz=10m
+
+        rh_profile = ds_column[rh].interpolate_na(dim=altitude)
+
+        # fit a straight line to the relative humidity values in the
+        # "ml-linfit" region
         X_height = rh_profile[ml_linfit][altitude]
         X = sm.add_constant(X_height.values)  # add intercept
         model = sm.OLS(
@@ -78,6 +91,7 @@ def calc_peakRH_linearize(
         # add +X index from mixed layer (ml) lower bound (lb), divided by dz=10m
         add_to_index_ml_lb = int(z_min / dz)
         idx_peaks_RH = idx_peaks_RH_raw + add_to_index_ml_lb
+
         obs_RH_peaks = rh_profile[idx_peaks_RH]  # RH values at the local maxima
         linearized_RH_values = linearized_RH[
             idx_peaks_RH
@@ -85,19 +99,44 @@ def calc_peakRH_linearize(
         min_diff_RH = min(
             abs(obs_RH_peaks - linearized_RH_values)
         )  # find minimum RH difference between observed peaks and linearized profile
-        h_peakRH_linfit[i] = min_diff_RH[altitude]
 
-    # set 'anomalously' high values > z_max to nan
-    # ... somewhat arbitrary choice
-    h_peakRH_linfit[h_peakRH_linfit > z_max] = np.nan
+        h_peakRH_linfit = min_diff_RH[altitude]
 
-    dims = list(ds.dims.keys())
-    # drop the height coord
-    del dims[dims.index(altitude)]
-    da = xr.DataArray(h_peakRH_linfit, dims=dims, coords={d: ds[d] for d in dims})
-    da.attrs["long_name"] = "mixed layer height (from RH peak and linearization)"
-    da.attrs["units"] = "m"
-    return da
+        if DEBUG:
+            ds_column[rh].plot(y=altitude, label="raw profile")
+            rh_profile.plot(y=altitude, label="raw profile interpolated")
+            plt.plot(linearized_RH, da_alt, label="RH interpolation in linfit region")
+
+            rh_profile.isel({altitude: idx_peaks_RH}).plot(
+                marker="x", linestyle="--", label="interpolated profile peaks", y=altitude
+            )
+
+            plt.axhline(h_peakRH_linfit, label="h_peakRH_linfit")
+
+            plt.legend()
+            plt.show()
+
+        # set 'anomalously' high values > z_max to nan
+        # ... somewhat arbitrary choice
+        if h_peakRH_linfit > z_max:
+            h_peakRH_linfit = np.nan
+
+        import ipdb
+
+        ipdb.set_trace()
+
+        return h_peakRH_linfit
+
+    ds_subset = get_fields(ds=ds, **{rh: "%", altitude: "m"})
+    da_ml_depth = apply_by_column(
+        ds=ds_subset, vertical_coord=vertical_coord, fn=_calc_mixed_layer_depth
+    )
+
+    da_ml_depth.attrs[
+        "long_name"
+    ] = "mixed layer height (from RH peak and linearization)"
+    da_ml_depth.attrs["units"] = "m"
+    return da_ml_depth
 
 
 def calc_from_gradient(
@@ -115,7 +154,6 @@ def calc_from_gradient(
     Outputs: DataArray containing mixed layer height from gradient method
 
     Note that function is quite slow and should be optimized
-
     """
 
     def calculateHmix_var(
diff --git a/eurec4a_environment/variables/calc_static_stability.py b/eurec4a_environment/variables/calc_static_stability.py
new file mode 100755
index 0000000..e62002e
--- /dev/null
+++ b/eurec4a_environment/variables/calc_static_stability.py
@@ -0,0 +1,45 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import xarray as xr
+
+
+def calc_static_stability(ds, pres="p", temp="T", theta="theta", altitude="height"):
+    """
+    Equation from:
+    https://www.ncl.ucar.edu/Document/Functions/Contributed/static_stability.shtml
+
+    S = -(T/theta)*d(theta)/dp. From Bluestein (1992), pg 197, eqn 4.3.8
+
+    Static stability measures the gravitational resistance of an atmosphere to
+    vertical displacements.  It results from fundamental buoyant adjustments,
+    and so it is determined by the vertical stratification of density or
+    potential temperature.
+    """
+
+    # !! edit these lines based on units function !!
+    # Celsius to Kelvin
+    if ds[temp].max().values < 100:
+        temp_K = ds[temp] + 273.15
+    else:
+        temp_K = ds[temp]
+
+    # convert Pa to hPa
+    if ds[pres].max().values > 1400:
+        ds[pres] = ds[pres] / 100
+
+    static_stability = -(temp_K / ds[theta]) * (
+        ds[theta].diff(dim=altitude) / ds[pres].diff(dim=altitude)
+    )
+    static_stability = static_stability.transpose(
+        transpose_coords=True
+    )  # same dimension order as original dataset
+
+    dims = list(static_stability.dims)
+    da = xr.DataArray(
+        static_stability, dims=dims, coords={d: static_stability[d] for d in dims}
+    )
+    da.attrs["long_name"] = "static stability"
+    da.attrs["units"] = "K/hPa"
+
+    return da
diff --git a/eurec4a_environment/variables/inversion_layers/inversion_layer_height.py b/eurec4a_environment/variables/inversion_layers/inversion_layer_height.py
new file mode 100755
index 0000000..154e390
--- /dev/null
+++ b/eurec4a_environment/variables/inversion_layers/inversion_layer_height.py
@@ -0,0 +1,89 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import xarray as xr
+
+
+def calc_inversion_static_stability(
+    ds, time="sounding", altitude="height", z_min=1000, z_max=3000
+):
+
+    """
+    Returns inversion height defined as the height where static stability profiles
+    first exceed 0.1 K/hPa between given heights
+    Adapted from Bony and Stevens 2018, Measuring Area-Averaged Vertical Motions with Dropsondes, p 772
+    To do: combine static stability criterion with one for moisture (i.e. hydrolapse)
+    """
+
+    from eurec4a_environment.variables.calc_static_stability import (
+        calc_static_stability,
+    )
+
+    ds_inversion = ds.sel({altitude: slice(z_min, z_max)})
+    da_static_stability = calc_static_stability(ds_inversion)
+
+    inversion_height_stability = np.zeros(len(da_static_stability[time]))
+    for i in range(len(da_static_stability[time])):
+        stability_sounding = da_static_stability.isel({time: i}).values
+        idx_stability = np.argmax(stability_sounding > 0.1)
+        inversion_height_stability[i] = da_static_stability[altitude][idx_stability]
+
+    dims = list(da_static_stability.dims)
+    # drop the height coord
+    del dims[dims.index(altitude)]
+
+    da = xr.DataArray(
+        inversion_height_stability,
+        dims=dims,
+        coords={d: da_static_stability[d] for d in dims},
+    )
+    da.attrs["long_name"] = "inversion base height, static stability > 0.1 K/hPa"
+    da.attrs["units"] = "m"
+
+    return da
+
+
+def calc_inversion_grad_RH_T(
+    ds,
+    time="sounding",
+    altitude="height",
+    temperature="T",
+    rh="rh",
+    z_min=1000,
+    z_max=2500,
+):
+
+    """
+    Inversion base following Grindinger et al, 1992; Cao et al 2007
+    dT/dz > 0 and dRH/dz < 0
+    To do: find all inversion bases and heights, not just first inversion base
+    """
+
+    ds_inversion = ds.sel({altitude: slice(z_min, z_max)})
+    gradient_RH = ds_inversion[rh].differentiate(coord=altitude)
+    gradient_T = ds_inversion[temperature].differentiate(coord=altitude)
+    inversion_levels = gradient_T[altitude].where(
+        (gradient_RH <= 0) & (gradient_T >= 0)
+    )
+
+    num_iters = len(inversion_levels[time])
+    inv_base_height_gradients = np.zeros(num_iters)
+    for i in range(num_iters):
+        alt_inv_da = inversion_levels.isel({time: i})
+        inversion_alts = alt_inv_da[~np.isnan(alt_inv_da)]
+        if inversion_alts[altitude].size == 0:
+            print(f"sounding {i} has no inversion below {z_max} m")
+            continue
+        inv_base_height_gradients[i] = inversion_alts[0]  # find first inversion base
+    inv_base_height_gradients[inv_base_height_gradients == 0] = np.nan
+
+    dims = list(ds.dims.keys())
+    # drop the height coord
+    del dims[dims.index(altitude)]
+    da = xr.DataArray(
+        inv_base_height_gradients, dims=dims, coords={d: ds[d] for d in dims}
+    )
+    da.attrs["long_name"] = "inversion base height from RH and T gradient"
+    da.attrs["units"] = "m"
+    return da
diff --git a/eurec4a_environment/variables/utils.py b/eurec4a_environment/variables/utils.py
index f37ed4c..1d06528 100644
--- a/eurec4a_environment/variables/utils.py
+++ b/eurec4a_environment/variables/utils.py
@@ -1,5 +1,5 @@
 def apply_by_column(ds, vertical_coord, fn):
-    if not vertical_coord in ds.dims:
+    if vertical_coord not in ds.dims:
         raise Exception(
             f"`{vertical_coord}` is not a dimension of the provided " "dataset"
         )
diff --git a/tests/test_boundary_layer.py b/tests/test_boundary_layer.py
index 45f8cc7..a88b239 100644
--- a/tests/test_boundary_layer.py
+++ b/tests/test_boundary_layer.py
@@ -3,6 +3,7 @@
 from eurec4a_environment.variables import boundary_layer
 from eurec4a_environment.variables.boundary_layer import inversion_height
 from eurec4a_environment import nomenclature as nom
+import eurec4a_environment.source_data
 
 
 def test_LCL_Bolton(ds_isentropic_test_profiles):
@@ -25,6 +26,26 @@ def test_mixed_layer_height_RHmax(ds_isentropic_test_profiles):
     assert np.allclose(da_rh_peak, z0)
 
 
+def test_mixed_layer_height_RH_lin():
+    ds = eurec4a_environment.source_data.open_joanne_dataset()
+    ds = ds.isel(sounding=slice(0, 10))
+    da_lin = boundary_layer.mixed_layer_height.calc_peakRH_linearize(ds)
+    assert len(da_lin) == 10
+    assert np.all(da_lin < 1500.0)
+
+
+def test_mixed_layer_height_gradient():
+    ds = eurec4a_environment.source_data.open_joanne_dataset()
+    ds = ds.isel(sounding=slice(0, 10))
+    # adjust threshold based on units 0.4 g/kg/10m or 0.0004 kg/kg/10m
+    da_gradient = boundary_layer.mixed_layer_height.calc_from_gradient(
+        ds, var=nom.SPECIFIC_HUMIDITY, threshold=0.4 / 1000
+    )
+    assert len(da_gradient) == 10
+    assert np.all(da_gradient < 1500.0)
+    assert da_gradient.units == "m"
+
+
 def test_inversion_height_gradient_RH(ds_isentropic_test_profiles):
     ds = ds_isentropic_test_profiles.copy()
     z_INV = 2000.0
diff --git a/tests/test_inversion_layers.py b/tests/test_inversion_layers.py
new file mode 100755
index 0000000..5620210
--- /dev/null
+++ b/tests/test_inversion_layers.py
@@ -0,0 +1,25 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import numpy as np
+from eurec4a_environment.variables.inversion_layers import inversion_layer_height
+import eurec4a_environment.source_data
+
+
+def test_inversion_static_stability():
+    ds = eurec4a_environment.source_data.open_joanne_dataset()
+    ds = ds.isel(sounding=slice(0, 10))
+    da_ss = inversion_layer_height.calc_inversion_static_stability(ds)
+    assert len(da_ss) == 10
+    assert np.all(da_ss < 5000.0)
+    assert da_ss.units == "m"
+
+
+def test_inversion_gradient_RH_T():
+    ds = eurec4a_environment.source_data.open_joanne_dataset()
+    ds = ds.isel(sounding=slice(0, 10))
+    # adjust threshold based on units 0.4 g/kg/10m or 0.0004 kg/kg/10m
+    da_gradRH_T = inversion_layer_height.calc_inversion_grad_RH_T(ds)
+    assert len(da_gradRH_T) == 10
+    assert np.all(da_gradRH_T < 5000.0)
+    assert da_gradRH_T.units == "m"
diff --git a/tests/test_static_stability.py b/tests/test_static_stability.py
new file mode 100755
index 0000000..2a6e6f1
--- /dev/null
+++ b/tests/test_static_stability.py
@@ -0,0 +1,15 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+
+from eurec4a_environment.variables.calc_static_stability import calc_static_stability
+import eurec4a_environment.source_data
+
+
+def test_static_stability_calculation():
+    ds = eurec4a_environment.source_data.open_joanne_dataset()
+    da_static_stability = calc_static_stability(ds)
+    # check that we get some sensible numbers
+    assert da_static_stability.units == "K/hPa"
+    assert da_static_stability.max() < 3
+    assert da_static_stability.min() > -3
diff --git a/tests/test_tropical.py b/tests/test_tropical.py
index ed070a7..d4966d2 100644
--- a/tests/test_tropical.py
+++ b/tests/test_tropical.py
@@ -99,7 +99,9 @@ def test_eis_joanne_profile():
     )
     ds[LTS_name] = da_lts
 
-    da_eis = tropical_variables.estimated_inversion_strength(ds=ds, LCL=LCL_name, LTS=LTS_name)
+    da_eis = tropical_variables.estimated_inversion_strength(
+        ds=ds, LCL=LCL_name, LTS=LTS_name
+    )
 
     assert da_eis.units == "K"
     assert np.all(da_eis > 0.0)