diff --git a/compliance_checker/cfutil.py b/compliance_checker/cfutil.py
index 8bcce331..c46976c9 100644
--- a/compliance_checker/cfutil.py
+++ b/compliance_checker/cfutil.py
@@ -8,9 +8,10 @@
 from collections import defaultdict
 from functools import lru_cache, partial
 
-from cf_units import Unit
 from importlib_resources import files
 
+from compliance_checker.units import UndefinedUnitError, units
+
 _UNITLESS_DB = None
 _SEA_NAMES = None
 
@@ -111,8 +112,8 @@ def is_dimensionless_standard_name(standard_name_table, standard_name):
         f".//entry[@id='{standard_name}']",
     )
     if found_standard_name is not None:
-        canonical_units = Unit(found_standard_name.find("canonical_units").text)
-        return canonical_units.is_dimensionless()
+        canonical_units = units(found_standard_name.find("canonical_units").text)
+        return canonical_units.dimensionless
     # if the standard name is not found, assume we need units for the time being
     else:
         return False
@@ -2037,8 +2038,8 @@ def units_convertible(units1, units2, reftimeistime=True):
     :param str units2: A string representing the units
     """
     try:
-        u1 = Unit(units1)
-        u2 = Unit(units2)
-    except ValueError:
+        u1 = units(units1)
+        u2 = units(units2)
+    except UndefinedUnitError:
         return False
-    return u1.is_convertible(u2)
+    return u1.is_compatible_with(u2)
diff --git a/compliance_checker/ioos.py b/compliance_checker/ioos.py
index 2deed991..19895c21 100644
--- a/compliance_checker/ioos.py
+++ b/compliance_checker/ioos.py
@@ -6,7 +6,6 @@
 from numbers import Number
 
 import validators
-from cf_units import Unit
 from lxml.etree import XPath
 from owslib.namespaces import Namespaces
 
@@ -28,6 +27,7 @@
     get_instrument_variables,
     get_z_variables,
 )
+from compliance_checker.units import units
 
 
 class IOOSBaseCheck(BaseCheck):
@@ -1379,12 +1379,12 @@ def check_vertical_coordinates(self, ds):
             )
 
             unit_def_set = {
-                Unit(unit_str).definition for unit_str in expected_unit_strs
+                str(units(unit_str).to_root_units()) for unit_str in expected_unit_strs
             }
 
             try:
-                units = Unit(units_str)
-                pass_stat = units.definition in unit_def_set
+                u = units(units_str)
+                pass_stat = str(u.to_root_units()) in unit_def_set
             # unknown unit not convertible to UDUNITS
             except ValueError:
                 pass_stat = False
diff --git a/compliance_checker/units.py b/compliance_checker/units.py
new file mode 100644
index 00000000..bcae5a5c
--- /dev/null
+++ b/compliance_checker/units.py
@@ -0,0 +1,124 @@
+"""Module to provide unit support via pint approximating UDUNITS/CF."""
+
+import functools
+import re
+
+import pint
+from pint import (  # noqa: F401
+    DimensionalityError,
+    UndefinedUnitError,
+    UnitStrippedWarning,
+)
+
+# from `xclim`'s unit support module with permission of the maintainers
+try:
+
+    @pint.register_unit_format("cf")
+    def short_formatter(unit, registry, **options):
+        """Return a CF-compliant unit string from a `pint` unit.
+
+        Parameters
+        ----------
+        unit : pint.UnitContainer
+            Input unit.
+        registry : pint.UnitRegistry
+            the associated registry
+        **options
+            Additional options (may be ignored)
+
+        Returns
+        -------
+        out : str
+            Units following CF-Convention, using symbols.
+        """
+        import re
+
+        # convert UnitContainer back to Unit
+        unit = registry.Unit(unit)
+        # Print units using abbreviations (millimeter -> mm)
+        s = f"{unit:~D}"
+
+        # Search and replace patterns
+        pat = r"(?P<inverse>(?:1 )?/ )?(?P<unit>\w+)(?: \*\* (?P<pow>\d))?"
+
+        def repl(m):
+            i, u, p = m.groups()
+            p = p or (1 if i else "")
+            neg = "-" if i else ""
+
+            return f"{u}{neg}{p}"
+
+        out, n = re.subn(pat, repl, s)
+
+        # Remove multiplications
+        out = out.replace(" * ", " ")
+        # Delta degrees:
+        out = out.replace("Δ°", "delta_deg")
+        return out.replace("percent", "%")
+
+except ImportError:
+    pass
+
+# ------
+# Reused with modification from MetPy under the terms of the BSD 3-Clause License.
+# Copyright (c) 2015,2017,2019 MetPy Developers.
+# Create registry, with preprocessors for UDUNITS-style powers (m2 s-2) and percent signs
+units = pint.UnitRegistry(
+    autoconvert_offset_to_baseunit=True,
+    preprocessors=[
+        lambda x: (
+            "count" if x == "1" else x
+        ),  # Should be salinity as well but all we can is that it is integer and dimensionless
+        lambda x: "S_K" if x.lower() in ["0.001", "1e-3"] else x,
+        functools.partial(
+            re.compile(
+                r"(?<=[A-Za-z])(?![A-Za-z])(?<![0-9\-][eE])(?<![0-9\-])(?=[0-9\-])",
+            ).sub,
+            "**",
+        ),
+        lambda string: string.replace("%", "percent"),
+    ],
+    force_ndarray_like=True,
+)
+# ----- end block copied from metpy
+
+# need to insert to make sure this is the first preprocessor
+# This ensures we convert integer `1` to string `"1"`, as needed by pint.
+units.preprocessors.insert(0, str)
+
+# -----
+units.define("percent = 0.01 = %")
+
+# Define commonly encountered units (both CF and non-CF) not defined by pint
+units.define("@alias meter = gpm")
+# ----- end block copied from metpy
+
+# -----
+# The following redefinitions were copied from xclim under the terms of their Apache-2 license
+# In pint, the default symbol for year is "a" which is not CF-compliant (stands for "are")
+units.define("year = 365.25 * day = yr")
+
+# Define commonly encountered units not defined by pint
+units.define("@alias degC = C = deg_C = Celsius = degrees_Celsius")
+units.define("@alias degK = deg_K")
+units.define("@alias day = d")
+units.define("@alias hour = h")  # Not the Planck constant...
+units.define(
+    "degrees_north = degree = degrees_north = degrees_N = degreesN = degree_north = degree_N = Degrees_N = degreeN",
+)
+units.define(
+    "degrees_east = degree = degrees_east = degrees_E = degreesE = degree_east = degree_E = Degrees_E = degreeE",
+)
+# degrees for grid_longitude / grid_latitude for grid_mappings
+units.define("degrees = degree = degrees")
+units.define("[speed] = [length] / [time]")
+# ----- end block copied from xclim
+
+# Add other specific aliases (by cf_xarray developers)
+units.define("practical_salinity_unit = [] = psu = PSU")
+units.define("parts_per_thousand = 0.001 = sea_water_knudsen_salinity = S_K")
+units.define("US_survey_foot = 0.304800609601219 m = US_survey_feet")
+# end of vendored code from MetPy
+
+# Set as application registry
+pint.set_application_registry(units)
diff --git a/requirements.txt b/requirements.txt
index 914faa1e..7de1880c 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -9,6 +9,7 @@ netcdf4>=1.6.4
 owslib>=0.8.3
 packaging
 pendulum>=1.2.4
+pint
 pygeoif>=0.6
 pyproj>=2.2.1
 regex>=2017.07.28