diff --git a/setup.cfg b/setup.cfg
index f5466efc466..7a12ad5fed1 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -71,7 +71,7 @@ max-line-length = 95
 
 [flake8]
 max-line-length = 95
-application-import-names = metpy
+application-import-names = metpy flake8_metpy
 import-order-style = google
 copyright-check = True
 copyright-author = MetPy Developers
@@ -83,15 +83,20 @@ docstring-convention = numpy
 exclude = docs build src/metpy/io/metar_parser.py
 select = A B C D E F H I M Q RST S T W B902
 ignore = F405 W503 RST902 SIM106
-per-file-ignores = examples/*.py: D T003 T001
-                   tutorials/*.py: D T003 T001
+per-file-ignores = examples/*.py: D MPY001 T003 T001
+                   tutorials/*.py: D MPY001 T003 T001
                    src/metpy/xarray.py: RST304
                    src/metpy/deprecation.py: C801
                    src/metpy/calc/*.py: RST306
                    src/metpy/interpolate/*.py: RST306
                    src/metpy/future.py: RST307
-		   src/metpy/constants.py: RST306
+                   src/metpy/constants.py: RST306
                    docs/doc-server.py: T001
+                   tests/*.py: MPY001
+
+[flake8:local-plugins]
+extension = MPY = flake8_metpy:MetPyChecker
+paths = ./tools/flake8_metpy
 
 [tool:pytest]
 # https://github.com/matplotlib/pytest-mpl/issues/69
diff --git a/src/metpy/calc/basic.py b/src/metpy/calc/basic.py
index 4e8a131025f..71e37d6bd9d 100644
--- a/src/metpy/calc/basic.py
+++ b/src/metpy/calc/basic.py
@@ -24,9 +24,9 @@
 exporter = Exporter(globals())
 
 # The following variables are constants for a standard atmosphere
-t0 = 288. * units.kelvin
-p0 = 1013.25 * units.hPa
-gamma = 6.5 * units('K/km')
+t0 = units.Quantity(288., 'kelvin')
+p0 = units.Quantity(1013.25, 'hPa')
+gamma = units.Quantity(6.5, 'K/km')
 
 
 @exporter.export
@@ -89,24 +89,24 @@ def wind_direction(u, v, convention='from'):
     of 0.
 
     """
-    wdir = 90. * units.deg - np.arctan2(-v, -u)
+    wdir = units.Quantity(90., 'deg') - np.arctan2(-v, -u)
     origshape = wdir.shape
     wdir = np.atleast_1d(wdir)
 
     # Handle oceanographic convection
     if convention == 'to':
-        wdir -= 180 * units.deg
+        wdir -= units.Quantity(180., 'deg')
     elif convention not in ('to', 'from'):
         raise ValueError('Invalid kwarg for "convention". Valid options are "from" or "to".')
 
     mask = wdir <= 0
     if np.any(mask):
-        wdir[mask] += 360. * units.deg
+        wdir[mask] += units.Quantity(360., 'deg')
     # avoid unintended modification of `pint.Quantity` by direct use of magnitude
     calm_mask = (np.asarray(u.magnitude) == 0.) & (np.asarray(v.magnitude) == 0.)
     # np.any check required for legacy numpy which treats 0-d False boolean index as zero
     if np.any(calm_mask):
-        wdir[calm_mask] = 0. * units.deg
+        wdir[calm_mask] = units.Quantity(0., 'deg')
     return wdir.reshape(origshape).to('degrees')
 
 
@@ -199,7 +199,7 @@ def windchill(temperature, speed, face_level_winds=False, mask_undefined=True):
         # noinspection PyAugmentAssignment
         speed = speed * 1.5
 
-    temp_limit, speed_limit = 10. * units.degC, 3 * units.mph
+    temp_limit, speed_limit = units.Quantity(10., 'degC'), units.Quantity(3, 'mph')
     speed_factor = speed.to('km/hr').magnitude ** 0.16
     wcti = units.Quantity((0.6215 + 0.3965 * speed_factor) * temperature.to('degC').magnitude
                           - 11.37 * speed_factor + 13.12, units.degC).to(temperature.units)
@@ -257,38 +257,39 @@ def heat_index(temperature, relative_humidity, mask_undefined=True):
     relative_humidity = np.atleast_1d(relative_humidity)
     # assign units to relative_humidity if they currently are not present
     if not hasattr(relative_humidity, 'units'):
-        relative_humidity = relative_humidity * units.dimensionless
-    delta = temperature.to(units.degF) - 0. * units.degF
+        relative_humidity = units.Quantity(relative_humidity, 'dimensionless')
+    delta = temperature.to(units.degF) - units.Quantity(0., 'degF')
     rh2 = relative_humidity * relative_humidity
     delta2 = delta * delta
 
     # Simplifed Heat Index -- constants converted for relative_humidity in [0, 1]
-    a = -10.3 * units.degF + 1.1 * delta + 4.7 * units.delta_degF * relative_humidity
+    a = (units.Quantity(-10.3, 'degF') + 1.1 * delta
+         + units.Quantity(4.7, 'delta_degF') * relative_humidity)
 
     # More refined Heat Index -- constants converted for relative_humidity in [0, 1]
-    b = (-42.379 * units.degF
+    b = (units.Quantity(-42.379, 'degF')
          + 2.04901523 * delta
-         + 1014.333127 * units.delta_degF * relative_humidity
+         + units.Quantity(1014.333127, 'delta_degF') * relative_humidity
          - 22.475541 * delta * relative_humidity
-         - 6.83783e-3 / units.delta_degF * delta2
-         - 5.481717e2 * units.delta_degF * rh2
-         + 1.22874e-1 / units.delta_degF * delta2 * relative_humidity
+         - units.Quantity(6.83783e-3, '1/delta_degF') * delta2
+         - units.Quantity(5.481717e2, 'delta_degF') * rh2
+         + units.Quantity(1.22874e-1, '1/delta_degF') * delta2 * relative_humidity
          + 8.5282 * delta * rh2
-         - 1.99e-2 / units.delta_degF * delta2 * rh2)
+         - units.Quantity(1.99e-2, '1/delta_degF') * delta2 * rh2)
 
     # Create return heat index
-    hi = np.full(np.shape(temperature), np.nan) * units.degF
+    hi = units.Quantity(np.full(np.shape(temperature), np.nan), 'degF')
     # Retain masked status of temperature with resulting heat index
     if hasattr(temperature, 'mask'):
         hi = masked_array(hi)
 
     # If T <= 40F, Heat Index is T
-    sel = (temperature <= 40. * units.degF)
+    sel = (temperature <= units.Quantity(40., 'degF'))
     if np.any(sel):
         hi[sel] = temperature[sel].to(units.degF)
 
     # If a < 79F and hi is unset, Heat Index is a
-    sel = (a < 79. * units.degF) & np.isnan(hi)
+    sel = (a < units.Quantity(79., 'degF')) & np.isnan(hi)
     if np.any(sel):
         hi[sel] = a[sel]
 
@@ -298,26 +299,28 @@ def heat_index(temperature, relative_humidity, mask_undefined=True):
         hi[sel] = b[sel]
 
     # Adjustment for RH <= 13% and 80F <= T <= 112F
-    sel = ((relative_humidity <= 13. * units.percent) & (temperature >= 80. * units.degF)
-           & (temperature <= 112. * units.degF))
+    sel = ((relative_humidity <= units.Quantity(13., 'percent'))
+           & (temperature >= units.Quantity(80., 'degF'))
+           & (temperature <= units.Quantity(112., 'degF')))
     if np.any(sel):
         rh15adj = ((13. - relative_humidity[sel] * 100.) / 4.
-                   * np.sqrt((17. * units.delta_degF
-                              - np.abs(delta[sel] - 95. * units.delta_degF))
-                             / 17. * units.delta_degF))
-        hi[sel] = hi[sel] - rh15adj
+                   * np.sqrt((units.Quantity(17., 'delta_degF')
+                              - np.abs(delta[sel] - units.Quantity(95., 'delta_degF')))
+                             / units.Quantity(17., 'delta_degF')))
+        hi[sel] = hi[sel] - units.Quantity(rh15adj, 'delta_degF')
 
     # Adjustment for RH > 85% and 80F <= T <= 87F
-    sel = ((relative_humidity > 85. * units.percent) & (temperature >= 80. * units.degF)
-           & (temperature <= 87. * units.degF))
+    sel = ((relative_humidity > units.Quantity(85., 'percent'))
+           & (temperature >= units.Quantity(80., 'degF'))
+           & (temperature <= units.Quantity(87., 'degF')))
     if np.any(sel):
-        rh85adj = 0.02 * (relative_humidity[sel] * 100. - 85.) * (87. * units.delta_degF
-                                                                  - delta[sel])
+        rh85adj = (0.02 * (relative_humidity[sel] * 100. - 85.)
+                   * (units.Quantity(87., 'delta_degF') - delta[sel]))
         hi[sel] = hi[sel] + rh85adj
 
     # See if we need to mask any undefined values
     if mask_undefined:
-        mask = np.array(temperature < 80. * units.degF)
+        mask = np.array(temperature < units.Quantity(80., 'degF'))
         if mask.any():
             hi = masked_array(hi, mask=mask)
     return hi
@@ -398,7 +401,7 @@ def apparent_temperature(temperature, relative_humidity, speed, face_level_winds
     # If no values are masked and provided temperature does not have a mask
     # we should return a non-masked array
     if not np.any(app_temperature.mask) and not hasattr(temperature, 'mask'):
-        app_temperature = np.array(app_temperature.m) * temperature.units
+        app_temperature = units.Quantity(np.array(app_temperature.m), temperature.units)
 
     if is_not_scalar:
         return app_temperature
@@ -1103,7 +1106,7 @@ def altimeter_to_station_pressure(altimeter_value, height):
 
     return ((altimeter_value ** n
              - ((p0.to(altimeter_value.units) ** n * gamma * height) / t0)) ** (1 / n)
-            + 0.3 * units.hPa)
+            + units.Quantity(0.3, 'hPa'))
 
 
 @exporter.export
diff --git a/src/metpy/calc/cross_sections.py b/src/metpy/calc/cross_sections.py
index eb05dd830ad..6d5cbe40b29 100644
--- a/src/metpy/calc/cross_sections.py
+++ b/src/metpy/calc/cross_sections.py
@@ -48,8 +48,8 @@ def distances_from_cross_section(cross):
         y = distance * np.cos(np.deg2rad(forward_az))
 
         # Build into DataArrays
-        x = xr.DataArray(x * units.meter, coords=lon.coords, dims=lon.dims)
-        y = xr.DataArray(y * units.meter, coords=lat.coords, dims=lat.dims)
+        x = xr.DataArray(units.Quantity(x, 'meter'), coords=lon.coords, dims=lon.dims)
+        y = xr.DataArray(units.Quantity(y, 'meter'), coords=lat.coords, dims=lat.dims)
 
     elif check_axis(cross.metpy.x, 'x') and check_axis(cross.metpy.y, 'y'):
 
@@ -88,8 +88,8 @@ def latitude_from_cross_section(cross):
             inverse=True,
             radians=False
         )[1]
-        latitude = xr.DataArray(latitude * units.degrees_north, coords=y.coords, dims=y.dims)
-        return latitude
+        return xr.DataArray(units.Quantity(latitude, 'degrees_north'), coords=y.coords,
+                            dims=y.dims)
 
 
 @exporter.export
diff --git a/src/metpy/calc/indices.py b/src/metpy/calc/indices.py
index 7277485e67e..a1c975b0978 100644
--- a/src/metpy/calc/indices.py
+++ b/src/metpy/calc/indices.py
@@ -78,8 +78,7 @@ def precipitable_water(pressure, dewpoint, *, bottom=None, top=None):
     w = mixing_ratio(saturation_vapor_pressure(dewpoint_layer), pres_layer)
 
     # Since pressure is in decreasing order, pw will be the opposite sign of that expected.
-    pw = -1. * (np.trapz(w.magnitude, pres_layer.magnitude) * (w.units * pres_layer.units)
-                / (mpconsts.g * mpconsts.rho_l))
+    pw = -np.trapz(w, pres_layer) / (mpconsts.g * mpconsts.rho_l)
     return pw.to('millimeters')
 
 
@@ -187,26 +186,29 @@ def bunkers_storm_motion(pressure, u, v, height):
 
     """
     # mean wind from sfc-6km
-    _, u_mean, v_mean = get_layer(pressure, u, v, height=height, depth=6000 * units('meter'))
-    wind_mean = [np.mean(u_mean).m, np.mean(v_mean).m] * u_mean.units
+    _, u_mean, v_mean = get_layer(pressure, u, v, height=height,
+                                  depth=units.Quantity(6000, 'meter'))
+    wind_mean = units.Quantity([np.mean(u_mean).m, np.mean(v_mean).m], u_mean.units)
 
     # mean wind from sfc-500m
-    _, u_500m, v_500m = get_layer(pressure, u, v, height=height, depth=500 * units('meter'))
-    wind_500m = [np.mean(u_500m).m, np.mean(v_500m).m] * u_500m.units
+    _, u_500m, v_500m = get_layer(pressure, u, v, height=height,
+                                  depth=units.Quantity(500, 'meter'))
+    wind_500m = units.Quantity([np.mean(u_500m).m, np.mean(v_500m).m], u_500m.units)
 
     # mean wind from 5.5-6km
-    _, u_5500m, v_5500m = get_layer(pressure, u, v, height=height, depth=500 * units('meter'),
-                                    bottom=height[0] + 5500 * units('meter'))
-    wind_5500m = [np.mean(u_5500m).m, np.mean(v_5500m).m] * u_5500m.units
+    _, u_5500m, v_5500m = get_layer(pressure, u, v, height=height,
+                                    depth=units.Quantity(500, 'meter'),
+                                    bottom=height[0] + units.Quantity(5500, 'meter'))
+    wind_5500m = units.Quantity([np.mean(u_5500m).m, np.mean(v_5500m).m], u_5500m.units)
 
     # Calculate the shear vector from sfc-500m to 5.5-6km
     shear = wind_5500m - wind_500m
 
     # Take the cross product of the wind shear and k, and divide by the vector magnitude and
-    # multiply by the deviaton empirically calculated in Bunkers (2000) (7.5 m/s)
+    # multiply by the deviation empirically calculated in Bunkers (2000) (7.5 m/s)
     shear_cross = concatenate([shear[1], -shear[0]])
-    shear_mag = np.hypot(*(arg.magnitude for arg in shear)) * shear.units
-    rdev = shear_cross * (7.5 * units('m/s').to(u.units) / shear_mag)
+    shear_mag = units.Quantity(np.hypot(*(arg.magnitude for arg in shear)), shear.units)
+    rdev = shear_cross * (units.Quantity(7.5, 'm/s').to(u.units) / shear_mag)
 
     # Add the deviations to the layer average wind to get the RM motion
     right_mover = wind_mean + rdev
@@ -309,12 +311,12 @@ def supercell_composite(mucape, effective_storm_helicity, effective_shear):
         Supercell composite
 
     """
-    effective_shear = np.clip(np.atleast_1d(effective_shear), None, 20 * units('m/s'))
-    effective_shear[effective_shear < 10 * units('m/s')] = 0 * units('m/s')
-    effective_shear = effective_shear / (20 * units('m/s'))
+    effective_shear = np.clip(np.atleast_1d(effective_shear), None, units.Quantity(20, 'm/s'))
+    effective_shear[effective_shear < units.Quantity(10, 'm/s')] = units.Quantity(0, 'm/s')
+    effective_shear = effective_shear / units.Quantity(20, 'm/s')
 
-    return ((mucape / (1000 * units('J/kg')))
-            * (effective_storm_helicity / (50 * units('m^2/s^2')))
+    return ((mucape / units.Quantity(1000, 'J/kg'))
+            * (effective_storm_helicity / units.Quantity(50, 'm^2/s^2'))
             * effective_shear).to('dimensionless')
 
 
@@ -361,17 +363,18 @@ def significant_tornado(sbcape, surface_based_lcl_height, storm_helicity_1km, sh
 
     """
     surface_based_lcl_height = np.clip(np.atleast_1d(surface_based_lcl_height),
-                                       1000 * units.m, 2000 * units.m)
-    surface_based_lcl_height[surface_based_lcl_height > 2000 * units.m] = 0 * units.m
-    surface_based_lcl_height = ((2000. * units.m - surface_based_lcl_height)
-                                / (1000. * units.m))
-    shear_6km = np.clip(np.atleast_1d(shear_6km), None, 30 * units('m/s'))
-    shear_6km[shear_6km < 12.5 * units('m/s')] = 0 * units('m/s')
-    shear_6km /= 20 * units('m/s')
-
-    return ((sbcape / (1500. * units('J/kg')))
+                                       units.Quantity(1000., 'm'), units.Quantity(2000., 'm'))
+    surface_based_lcl_height[surface_based_lcl_height > units.Quantity(2000., 'm')] = \
+        units.Quantity(0, 'm')
+    surface_based_lcl_height = ((units.Quantity(2000., 'm') - surface_based_lcl_height)
+                                / units.Quantity(1000., 'm'))
+    shear_6km = np.clip(np.atleast_1d(shear_6km), None, units.Quantity(30, 'm/s'))
+    shear_6km[shear_6km < units.Quantity(12.5, 'm/s')] = units.Quantity(0, 'm/s')
+    shear_6km /= units.Quantity(20, 'm/s')
+
+    return ((sbcape / units.Quantity(1500., 'J/kg'))
             * surface_based_lcl_height
-            * (storm_helicity_1km / (150. * units('m^2/s^2')))
+            * (storm_helicity_1km / units.Quantity(150., 'm^2/s^2'))
             * shear_6km)
 
 
@@ -436,7 +439,7 @@ def critical_angle(pressure, u, v, height, u_storm, v_storm):
     v = v[sort_inds]
 
     # Calculate sfc-500m shear vector
-    shr5 = bulk_shear(pressure, u, v, height=height, depth=500 * units('meter'))
+    shr5 = bulk_shear(pressure, u, v, height=height, depth=units.Quantity(500, 'meter'))
 
     # Make everything relative to the sfc wind orientation
     umn = u_storm - u[0]
@@ -445,6 +448,6 @@ def critical_angle(pressure, u, v, height, u_storm, v_storm):
     vshr = np.asarray([shr5[0].magnitude, shr5[1].magnitude])
     vsm = np.asarray([umn.magnitude, vmn.magnitude])
     angle_c = np.dot(vshr, vsm) / (np.linalg.norm(vshr) * np.linalg.norm(vsm))
-    critical_angle = np.arccos(angle_c) * units('radian')
+    critical_angle = units.Quantity(np.arccos(angle_c), 'radian')
 
     return critical_angle.to('degrees')
diff --git a/src/metpy/calc/kinematics.py b/src/metpy/calc/kinematics.py
index 5e2281cdefc..5de5a8501a7 100644
--- a/src/metpy/calc/kinematics.py
+++ b/src/metpy/calc/kinematics.py
@@ -563,8 +563,7 @@ def montgomery_streamfunction(height, temperature):
 @preprocess_and_wrap()
 @check_units('[length]', '[speed]', '[speed]', '[length]',
              bottom='[length]', storm_u='[speed]', storm_v='[speed]')
-def storm_relative_helicity(height, u, v, depth, *, bottom=0 * units.m,
-                            storm_u=0 * units('m/s'), storm_v=0 * units('m/s')):
+def storm_relative_helicity(height, u, v, depth, *, bottom=None, storm_u=None, storm_v=None):
     # Partially adapted from similar SharpPy code
     r"""Calculate storm relative helicity.
 
@@ -622,6 +621,13 @@ def storm_relative_helicity(height, u, v, depth, *, bottom=0 * units.m,
        ``storm_v`` parameters to keyword-only arguments
 
     """
+    if bottom is None:
+        bottom = units.Quantity(0, 'm')
+    if storm_u is None:
+        storm_u = units.Quantity(0, 'm/s')
+    if storm_v is None:
+        storm_v = units.Quantity(0, 'm/s')
+
     _, u, v = get_layer_heights(height, depth, u, v, with_agl=True, bottom=bottom)
 
     storm_relative_u = u - storm_u
@@ -634,10 +640,10 @@ def storm_relative_helicity(height, u, v, depth, *, bottom=0 * units.m,
     # mask will return a masked value rather than 0. See numpy/numpy#11736
     positive_srh = int_layers[int_layers.magnitude > 0.].sum()
     if np.ma.is_masked(positive_srh):
-        positive_srh = 0.0 * units('meter**2 / second**2')
+        positive_srh = units.Quantity(0.0, 'meter**2 / second**2')
     negative_srh = int_layers[int_layers.magnitude < 0.].sum()
     if np.ma.is_masked(negative_srh):
-        negative_srh = 0.0 * units('meter**2 / second**2')
+        negative_srh = units.Quantity(0.0, 'meter**2 / second**2')
 
     return (positive_srh.to('meter ** 2 / second ** 2'),
             negative_srh.to('meter ** 2 / second ** 2'),
@@ -789,8 +795,8 @@ def potential_vorticity_baroclinic(
         (np.shape(potential_temperature)[y_dim] == 1)
         and (np.shape(potential_temperature)[x_dim] == 1)
     ):
-        dthtady = 0 * units.K / units.m  # axis=y_dim only has one dimension
-        dthtadx = 0 * units.K / units.m  # axis=x_dim only has one dimension
+        dthtady = units.Quantity(0, 'K/m')  # axis=y_dim only has one dimension
+        dthtadx = units.Quantity(0, 'K/m')  # axis=x_dim only has one dimension
     else:
         dthtady = first_derivative(potential_temperature, delta=dy, axis=y_dim)
         dthtadx = first_derivative(potential_temperature, delta=dx, axis=x_dim)
@@ -798,8 +804,7 @@ def potential_vorticity_baroclinic(
     dvdp = first_derivative(v, x=pressure, axis=vertical_dim)
 
     return (-mpconsts.g * (dudp * dthtady - dvdp * dthtadx
-                           + avor * dthtadp)).to(units.kelvin * units.meter**2
-                                                 / (units.second * units.kilogram))
+                           + avor * dthtadp)).to('K * m**2 / (s * kg)')
 
 
 @exporter.export
diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py
index 4357c580e9d..4bb2791eaba 100644
--- a/src/metpy/calc/thermo.py
+++ b/src/metpy/calc/thermo.py
@@ -21,7 +21,7 @@
 
 exporter = Exporter(globals())
 
-sat_pressure_0c = 6.112 * units.millibar
+sat_pressure_0c = units.Quantity(6.112, 'millibar')
 
 
 @exporter.export
@@ -408,7 +408,8 @@ def _lcl_iter(p, p0, w, t):
 
     # np.isclose needed if surface is LCL due to precision error with np.log in dewpoint.
     # Causes issues with parcel_profile_with_lcl if removed. Issue #1187
-    lcl_p = np.where(np.isclose(lcl_p, pressure.m), pressure.m, lcl_p) * pressure.units
+    lcl_p = units.Quantity(np.where(np.isclose(lcl_p, pressure.m), pressure.m, lcl_p),
+                           pressure.units)
 
     return lcl_p, globals()['dewpoint'](vapor_pressure(lcl_p, w)).to(temperature.units)
 
@@ -508,7 +509,8 @@ def lfc(pressure, temperature, dewpoint, parcel_temperature_profile=None, dewpoi
         mask = pressure < this_lcl[0]
         if np.all(_less_or_close(parcel_temperature_profile[mask], temperature[mask])):
             # LFC doesn't exist
-            x, y = np.nan * pressure.units, np.nan * temperature.units
+            x = units.Quantity(np.nan, pressure.units)
+            y = units.Quantity(np.nan, temperature.units)
         else:  # LFC = LCL
             x, y = this_lcl
         return x, y
@@ -522,7 +524,8 @@ def lfc(pressure, temperature, dewpoint, parcel_temperature_profile=None, dewpoi
                                             temperature[1:], direction='decreasing',
                                             log_x=True)
             if np.min(el_pres) > this_lcl[0]:
-                x, y = np.nan * pressure.units, np.nan * temperature.units
+                x = units.Quantity(np.nan, pressure.units)
+                y = units.Quantity(np.nan, temperature.units)
             else:
                 x, y = this_lcl
             return x, y
@@ -666,7 +669,8 @@ def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which='
 
     # If the top of the sounding parcel is warmer than the environment, there is no EL
     if parcel_temperature_profile[-1] > temperature[-1]:
-        return np.nan * pressure.units, np.nan * temperature.units
+        return (units.Quantity(np.nan, pressure.units),
+                units.Quantity(np.nan, temperature.units))
 
     # Interpolate in log space to find the appropriate pressure - units have to be stripped
     # and reassigned to allow np.log() to function properly.
@@ -679,7 +683,8 @@ def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which='
                                         parcel_temperature_profile, temperature, dewpoint,
                                         intersect_type='EL')
     else:
-        return np.nan * pressure.units, np.nan * temperature.units
+        return (units.Quantity(np.nan, pressure.units),
+                units.Quantity(np.nan, temperature.units))
 
 
 @exporter.export
@@ -896,7 +901,7 @@ def _insert_lcl_level(pressure, temperature, lcl_pressure):
     # Pressure needs to be increasing for searchsorted, so flip it and then convert
     # the index back to the original array
     loc = pressure.size - pressure[::-1].searchsorted(lcl_pressure)
-    return temperature.units * np.insert(temperature.m, loc, interp_temp.m)
+    return units.Quantity(np.insert(temperature.m, loc, interp_temp.m), temperature.units)
 
 
 @exporter.export
@@ -971,8 +976,8 @@ def saturation_vapor_pressure(temperature):
     """
     # Converted from original in terms of C to use kelvin. Using raw absolute values of C in
     # a formula plays havoc with units support.
-    return sat_pressure_0c * np.exp(17.67 * (temperature - 273.15 * units.kelvin)
-                                    / (temperature - 29.65 * units.kelvin))
+    return sat_pressure_0c * np.exp(17.67 * (temperature - units.Quantity(273.15, 'kelvin'))
+                                    / (temperature - units.Quantity(29.65, 'kelvin')))
 
 
 @exporter.export
@@ -1041,7 +1046,8 @@ def dewpoint(vapor_pressure):
 
     """
     val = np.log(vapor_pressure / sat_pressure_0c)
-    return 0. * units.degC + 243.5 * units.delta_degC * val / (17.67 - val)
+    return (units.Quantity(0., 'degC')
+            + units.Quantity(243.5, 'delta_degC') * val / (17.67 - val))
 
 
 @exporter.export
@@ -1380,7 +1386,7 @@ def density(pressure, temperature, mixing_ratio, molecular_weight_ratio=mpconsts
 
     """
     virttemp = virtual_temperature(temperature, mixing_ratio, molecular_weight_ratio)
-    return (pressure / (mpconsts.Rd * virttemp)).to(units.kilogram / units.meter ** 3)
+    return (pressure / (mpconsts.Rd * virttemp)).to('kg/m**3')
 
 
 @exporter.export
@@ -1441,7 +1447,7 @@ def relative_humidity_wet_psychrometric(pressure, dry_bulb_temperature, wet_bulb
 )
 @check_units('[pressure]', '[temperature]', '[temperature]')
 def psychrometric_vapor_pressure_wet(pressure, dry_bulb_temperature, wet_bulb_temperature,
-                                     psychrometer_coefficient=6.21e-4 / units.kelvin):
+                                     psychrometer_coefficient=None):
     r"""Calculate the vapor pressure with wet bulb and dry bulb temperatures.
 
     This uses a psychrometric relationship as outlined in [WMO8]_, with
@@ -1490,6 +1496,8 @@ def psychrometric_vapor_pressure_wet(pressure, dry_bulb_temperature, wet_bulb_te
     saturation_vapor_pressure
 
     """
+    if psychrometer_coefficient is None:
+        psychrometer_coefficient = units.Quantity(6.21e-4, '1/K')
     return (saturation_vapor_pressure(wet_bulb_temperature) - psychrometer_coefficient
             * pressure * (dry_bulb_temperature - wet_bulb_temperature).to('kelvin'))
 
@@ -1788,7 +1796,7 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom'
 
     # If there is no LFC, no need to proceed.
     if np.isnan(lfc_pressure):
-        return 0 * units('J/kg'), 0 * units('J/kg')
+        return units.Quantity(0, 'J/kg'), units.Quantity(0, 'J/kg')
     else:
         lfc_pressure = lfc_pressure.magnitude
 
@@ -1814,7 +1822,7 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom'
     x_clipped = x[p_mask].magnitude
     y_clipped = y[p_mask].magnitude
     cape = (mpconsts.Rd
-            * (np.trapz(y_clipped, np.log(x_clipped)) * units.degK)).to(units('J/kg'))
+            * units.Quantity(np.trapz(y_clipped, np.log(x_clipped)), 'K')).to(units('J/kg'))
 
     # CIN
     # Only use data between the surface and LFC for calculation
@@ -1822,11 +1830,11 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom'
     x_clipped = x[p_mask].magnitude
     y_clipped = y[p_mask].magnitude
     cin = (mpconsts.Rd
-           * (np.trapz(y_clipped, np.log(x_clipped)) * units.degK)).to(units('J/kg'))
+           * units.Quantity(np.trapz(y_clipped, np.log(x_clipped)), 'K')).to(units('J/kg'))
 
     # Set CIN to 0 if it's returned as a positive value (#1190)
-    if cin > 0 * units('J/kg'):
-        cin = 0 * units('J/kg')
+    if cin > units.Quantity(0, 'J/kg'):
+        cin = units.Quantity(0, 'J/kg')
     return cape, cin
 
 
@@ -1853,7 +1861,8 @@ def _find_append_zero_crossings(x, y):
         Y values of data
 
     """
-    crossings = find_intersections(x[1:], y[1:], y.units * np.zeros_like(y[1:]), log_x=True)
+    crossings = find_intersections(x[1:], y[1:],
+                                   units.Quantity(np.zeros_like(y[1:]), y.units), log_x=True)
     x = concatenate((x, crossings[0]))
     y = concatenate((y, crossings[1]))
 
@@ -1872,8 +1881,8 @@ def _find_append_zero_crossings(x, y):
 @exporter.export
 @preprocess_and_wrap()
 @check_units('[pressure]', '[temperature]', '[temperature]')
-def most_unstable_parcel(pressure, temperature, dewpoint, height=None,
-                         bottom=None, depth=300 * units.hPa):
+def most_unstable_parcel(pressure, temperature, dewpoint, height=None, bottom=None,
+                         depth=None):
     """
     Determine the most unstable parcel in a layer.
 
@@ -1924,6 +1933,8 @@ def most_unstable_parcel(pressure, temperature, dewpoint, height=None,
        Renamed ``heights`` parameter to ``height``
 
     """
+    if depth is None:
+        depth = units.Quantity(300, 'hPa')
     p_layer, t_layer, td_layer = get_layer(pressure, temperature, dewpoint, bottom=bottom,
                                            depth=depth, height=height, interpolate=False)
     theta_e = equivalent_potential_temperature(p_layer, t_layer, td_layer)
@@ -2016,7 +2027,8 @@ def _isen_iter(iter_log_p, isentlevs_nd, ka, a, b, pok):
     slices = [np.newaxis] * ndim
     slices[vertical_dim] = slice(None)
     slices = tuple(slices)
-    pres = np.broadcast_to(pres[slices].magnitude, temperature.shape) * pres.units
+    pres = units.Quantity(np.broadcast_to(pres[slices].magnitude, temperature.shape),
+                          pres.units)
 
     # Sort input data
     sort_pres = np.argsort(pres.m, axis=vertical_dim)
@@ -2077,11 +2089,11 @@ def _isen_iter(iter_log_p, isentlevs_nd, ka, a, b, pok):
     isentprs[~(good & _less_or_close(isentprs, np.max(pres.m)))] = np.nan
 
     # create list for storing output data
-    ret = [isentprs * units.hPa]
+    ret = [units.Quantity(isentprs, 'hPa')]
 
     # if temperature_out = true, calculate temperature and output as last item in list
     if temperature_out:
-        ret.append((isentlevs_nd / ((mpconsts.P0.m / isentprs) ** ka)) * units.kelvin)
+        ret.append(units.Quantity((isentlevs_nd / ((mpconsts.P0.m / isentprs) ** ka)), 'K'))
 
     # do an interpolation for each additional argument
     if args:
@@ -2350,7 +2362,7 @@ def mixed_layer_cape_cin(pressure, temperature, dewpoint, **kwargs):
     Quantities even when given xarray DataArray profiles.
 
     """
-    depth = kwargs.get('depth', 100 * units.hPa)
+    depth = kwargs.get('depth', units.Quantity(100, 'hPa'))
     parcel_pressure, parcel_temp, parcel_dewpoint = mixed_parcel(pressure, temperature,
                                                                  dewpoint, **kwargs)
 
@@ -2370,7 +2382,7 @@ def mixed_layer_cape_cin(pressure, temperature, dewpoint, **kwargs):
 @preprocess_and_wrap()
 @check_units('[pressure]', '[temperature]', '[temperature]')
 def mixed_parcel(pressure, temperature, dewpoint, parcel_start_pressure=None,
-                 height=None, bottom=None, depth=100 * units.hPa, interpolate=True):
+                 height=None, bottom=None, depth=None, interpolate=True):
     r"""Calculate the properties of a parcel mixed from a layer.
 
     Determines the properties of an air parcel that is the result of complete mixing of a
@@ -2429,6 +2441,9 @@ def mixed_parcel(pressure, temperature, dewpoint, parcel_start_pressure=None,
     if not parcel_start_pressure:
         parcel_start_pressure = pressure[0]
 
+    if depth is None:
+        depth = units.Quantity(100, 'hPa')
+
     # Calculate the potential temperature and mixing ratio over the layer
     theta = potential_temperature(pressure, temperature)
     mixing_ratio = saturation_mixing_ratio(pressure, dewpoint)
@@ -2455,8 +2470,7 @@ def mixed_parcel(pressure, temperature, dewpoint, parcel_start_pressure=None,
 @exporter.export
 @preprocess_and_wrap()
 @check_units('[pressure]')
-def mixed_layer(pressure, *args, height=None, bottom=None, depth=100 * units.hPa,
-                interpolate=True):
+def mixed_layer(pressure, *args, height=None, bottom=None, depth=None, interpolate=True):
     r"""Mix variable(s) over a layer, yielding a mass-weighted average.
 
     This function will integrate a data variable with respect to pressure and determine the
@@ -2499,6 +2513,8 @@ def mixed_layer(pressure, *args, height=None, bottom=None, depth=100 * units.hPa
        Renamed ``p``, ``heights`` parameters to ``pressure``, ``height``
 
     """
+    if depth is None:
+        depth = units.Quantity(100, 'hPa')
     layer = get_layer(pressure, *args, height=height, bottom=bottom,
                       depth=depth, interpolate=interpolate)
     p_layer = layer[0]
@@ -2507,8 +2523,8 @@ def mixed_layer(pressure, *args, height=None, bottom=None, depth=100 * units.hPa
     ret = []
     for datavar_layer in datavars_layer:
         actual_depth = abs(p_layer[0] - p_layer[-1])
-        ret.append((-1. / actual_depth.m) * np.trapz(datavar_layer.m, p_layer.m)
-                   * datavar_layer.units)
+        ret.append(units.Quantity(np.trapz(datavar_layer.m, p_layer.m) / -actual_depth.m,
+                   datavar_layer.units))
     return ret
 
 
@@ -2675,8 +2691,9 @@ def thickness_hydrostatic(pressure, temperature, mixing_ratio=None,
             layer_virttemp = virtual_temperature(layer_temp, layer_w, molecular_weight_ratio)
 
     # Take the integral (with unit handling) and return the result in meters
-    return (- mpconsts.Rd / mpconsts.g * np.trapz(
-        layer_virttemp.m_as('K'), x=np.log(layer_p.m_as('hPa'))) * units.K).to('m')
+    integral = units.Quantity(np.trapz(layer_virttemp.m_as('K'), np.log(layer_p.m_as('hPa'))),
+                              units.K)
+    return (-mpconsts.Rd / mpconsts.g * integral).to('m')
 
 
 @exporter.export
@@ -2949,17 +2966,16 @@ def wet_bulb_temperature(pressure, temperature, dewpoint):
                    flags=['buffered'])
 
     for press, lpress, ltemp, ret in it:
-        press = press * pressure.units
-        lpress = lpress * lcl_press.units
-        ltemp = ltemp * lcl_temp.units
-        moist_adiabat_temperatures = moist_lapse(press, ltemp, lpress)
+        moist_adiabat_temperatures = moist_lapse(units.Quantity(press, pressure.units),
+                                                 units.Quantity(ltemp, lcl_temp.units),
+                                                 units.Quantity(lpress, lcl_press.units))
         ret[...] = moist_adiabat_temperatures.magnitude
 
     # If we started with a scalar, return a scalar
     ret = it.operands[3]
     if ret.size == 1:
         ret = ret[0]
-    return ret * moist_adiabat_temperatures.units
+    return units.Quantity(ret, moist_adiabat_temperatures.units)
 
 
 @exporter.export
@@ -3210,7 +3226,7 @@ def lifted_index(pressure, temperature, parcel_profile):
 
     """
     # find the index for the 500 hPa pressure level.
-    idx = np.where(pressure == 500 * units.hPa)
+    idx = np.where(pressure == units.Quantity(500, 'hPa'))
     # find the measured temperature at 500 hPa.
     T500 = temperature[idx]
     # find the parcel profile temperature at 500 hPa.
diff --git a/src/metpy/calc/tools.py b/src/metpy/calc/tools.py
index 3e10d89d065..5b92f08432b 100644
--- a/src/metpy/calc/tools.py
+++ b/src/metpy/calc/tools.py
@@ -30,8 +30,8 @@
     UND
 )  # note the order matters!
 
-MAX_DEGREE_ANGLE = 360 * units.degree
-BASE_DEGREE_MULTIPLIER = 22.5 * units.degree
+MAX_DEGREE_ANGLE = units.Quantity(360, 'degree')
+BASE_DEGREE_MULTIPLIER = units.Quantity(22.5, 'degree')
 
 DIR_DICT = {dir_str: i * BASE_DEGREE_MULTIPLIER for i, dir_str in enumerate(DIR_STRS)}
 DIR_DICT[UND] = np.nan
@@ -388,9 +388,9 @@ def _get_bound_pressure_height(pressure, bound, height=None, interpolate=True):
                     # Need to cast back to the input type since interp (up to at least numpy
                     # 1.13 always returns float64. This can cause upstream users problems,
                     # resulting in something like np.append() to upcast.
-                    bound_pressure = (np.interp(np.atleast_1d(bound.m), height.m,
-                                                pressure.m).astype(np.result_type(bound))
-                                      * pressure.units)
+                    bound_pressure = units.Quantity(np.interp(np.atleast_1d(bound.m), height.m,
+                                                    pressure.m).astype(np.result_type(bound)),
+                                                    pressure.units)
                 else:
                     idx = (np.abs(height - bound)).argmin()
                     bound_pressure = pressure[idx]
@@ -521,8 +521,7 @@ def get_layer_heights(height, depth, *args, bottom=None, interpolate=True, with_
 @exporter.export
 @preprocess_and_wrap()
 @check_units('[pressure]')
-def get_layer(pressure, *args, height=None, bottom=None, depth=100 * units.hPa,
-              interpolate=True):
+def get_layer(pressure, *args, height=None, bottom=None, depth=None, interpolate=True):
     r"""Return an atmospheric layer from upper air data with the requested bottom and depth.
 
     This function will subset an upper air dataset to contain only the specified layer. The
@@ -566,7 +565,7 @@ def get_layer(pressure, *args, height=None, bottom=None, depth=100 * units.hPa,
     """
     # If we get the depth kwarg, but it's None, set it to the default as well
     if depth is None:
-        depth = 100 * units.hPa
+        depth = units.Quantity(100, 'hPa')
 
     # Make sure pressure and datavars are the same length
     for datavar in args:
@@ -607,9 +606,11 @@ def get_layer(pressure, *args, height=None, bottom=None, depth=100 * units.hPa,
     if interpolate:
         # If we don't have the bottom or top requested, append them
         if not np.any(np.isclose(top_pressure, p_interp)):
-            p_interp = np.sort(np.append(p_interp.m, top_pressure.m)) * pressure.units
+            p_interp = units.Quantity(np.sort(np.append(p_interp.m, top_pressure.m)),
+                                      pressure.units)
         if not np.any(np.isclose(bottom_pressure, p_interp)):
-            p_interp = np.sort(np.append(p_interp.m, bottom_pressure.m)) * pressure.units
+            p_interp = units.Quantity(np.sort(np.append(p_interp.m, bottom_pressure.m)),
+                                      pressure.units)
 
     ret.append(p_interp[::-1])
 
@@ -847,7 +848,7 @@ def lat_lon_grid_deltas(longitude, latitude, x_dim=-1, y_dim=-2, geod=None):
                               latitude[take_x(slice(1, None))])
     dx[(forward_az < 0.) | (forward_az > 180.)] *= -1
 
-    return dx * units.meter, dy * units.meter
+    return units.Quantity(dx, 'meter'), units.Quantity(dy, 'meter')
 
 
 @exporter.export
@@ -1290,7 +1291,7 @@ def _process_deriv_args(f, axis, x, delta):
             delta_units = getattr(delta, 'units', None)
             delta = np.broadcast_to(delta, diff_size, subok=True)
             if not hasattr(delta, 'units') and delta_units is not None:
-                delta = delta * delta_units
+                delta = units.Quantity(delta, delta_units)
         else:
             delta = _broadcast_to_axis(delta, axis, n)
     elif x is not None:
@@ -1393,13 +1394,9 @@ def angle_to_direction(input_angle, full=False, level=3):
     else:
         scalar = False
 
-    # clean any numeric strings, negatives, and None
-    # does not handle strings with alphabet
-    input_angle = np.array(input_angle).astype(float)
-    with np.errstate(invalid='ignore'):  # warns about the np.nan
-        input_angle[np.where(input_angle < 0)] = np.nan
-
-    input_angle = input_angle * origin_units
+    # clean any numeric strings, negatives, and None does not handle strings with alphabet
+    input_angle = units.Quantity(np.array(input_angle).astype(float), origin_units)
+    input_angle[input_angle < 0] = units.Quantity(np.nan, origin_units)
 
     # normalizer used for angles > 360 degree to normalize between 0 - 360
     normalizer = np.array(input_angle.m / MAX_DEGREE_ANGLE.m, dtype=int)
diff --git a/src/metpy/constants.py b/src/metpy/constants.py
index 8d46ae5ddf2..a4bc64c0bf3 100644
--- a/src/metpy/constants.py
+++ b/src/metpy/constants.py
@@ -78,49 +78,49 @@
 # Export all the variables defined in this block
 with exporter:
     # Earth
-    earth_gravity = g = 9.80665 * units('m / s^2')
-    Re = earth_avg_radius = 6371008.7714 * units('m')
-    G = gravitational_constant = 6.67430e-11 * units('m^3 / kg / s^2')
-    GM = geocentric_gravitational_constant = 3986005e8 * units('m^3 / s^2')
-    omega = earth_avg_angular_vel = 7292115e-11 * units('rad / s')
-    d = earth_sfc_avg_dist_sun = 149597870700. * units('m')
-    S = earth_solar_irradiance = 1360.8 * units('W / m^2')
-    delta = earth_max_declination = 23.45 * units('degrees')
-    earth_orbit_eccentricity = 0.0167 * units('dimensionless')
+    earth_gravity = g = units.Quantity(9.80665, 'm / s^2')
+    Re = earth_avg_radius = units.Quantity(6371008.7714, 'm')
+    G = gravitational_constant = units.Quantity(6.67430e-11, 'm^3 / kg / s^2')
+    GM = geocentric_gravitational_constant = units.Quantity(3986005e8, 'm^3 / s^2')
+    omega = earth_avg_angular_vel = units.Quantity(7292115e-11, 'rad / s')
+    d = earth_sfc_avg_dist_sun = units.Quantity(149597870700., 'm')
+    S = earth_solar_irradiance = units.Quantity(1360.8, 'W / m^2')
+    delta = earth_max_declination = units.Quantity(23.45, 'degrees')
+    earth_orbit_eccentricity = units.Quantity(0.0167, 'dimensionless')
     earth_mass = me = geocentric_gravitational_constant / gravitational_constant
 
     # molar gas constant
-    R = 8.314462618 * units('J / mol / K')
+    R = units.Quantity(8.314462618, 'J / mol / K')
 
     # Water
-    Mw = water_molecular_weight = 18.015268 * units('g / mol')
+    Mw = water_molecular_weight = units.Quantity(18.015268, 'g / mol')
     Rv = water_gas_constant = R / Mw
-    rho_l = density_water = 999.97495 * units('kg / m^3')
-    wv_specific_heat_ratio = 1.330 * units('dimensionless')
+    rho_l = density_water = units.Quantity(999.97495, 'kg / m^3')
+    wv_specific_heat_ratio = units.Quantity(1.330, 'dimensionless')
     Cp_v = wv_specific_heat_press = (
         wv_specific_heat_ratio * Rv / (wv_specific_heat_ratio - 1)
     )
     Cv_v = wv_specific_heat_vol = Cp_v / wv_specific_heat_ratio
-    Cp_l = water_specific_heat = 4.2194 * units('kJ / kg / K')
-    Lv = water_heat_vaporization = 2.50084e6 * units('J / kg')
-    Lf = water_heat_fusion = 3.337e5 * units('J / kg')
-    Cp_i = ice_specific_heat = 2090 * units('J / kg / K')
-    rho_i = density_ice = 917 * units('kg / m^3')
+    Cp_l = water_specific_heat = units.Quantity(4.2194, 'kJ / kg / K')
+    Lv = water_heat_vaporization = units.Quantity(2.50084e6, 'J / kg')
+    Lf = water_heat_fusion = units.Quantity(3.337e5, 'J / kg')
+    Cp_i = ice_specific_heat = units.Quantity(2090, 'J / kg / K')
+    rho_i = density_ice = units.Quantity(917, 'kg / m^3')
 
     # Dry air
-    Md = dry_air_molecular_weight = 28.96546e-3 * units('kg / mol')
+    Md = dry_air_molecular_weight = units.Quantity(28.96546e-3, 'kg / mol')
     Rd = dry_air_gas_constant = R / Md
-    dry_air_spec_heat_ratio = 1.4 * units('dimensionless')
+    dry_air_spec_heat_ratio = units.Quantity(1.4, 'dimensionless')
     Cp_d = dry_air_spec_heat_press = (
         dry_air_spec_heat_ratio * Rd / (dry_air_spec_heat_ratio - 1)
     )
     Cv_d = dry_air_spec_heat_vol = Cp_d / dry_air_spec_heat_ratio
     rho_d = dry_air_density_stp = (
-        (1000. * units('mbar')) / (Rd * 273.15 * units('K'))
+        units.Quantity(1000., 'mbar') / (Rd * units.Quantity(273.15, 'K'))
     ).to('kg / m^3')
 
     # General meteorology constants
-    P0 = pot_temp_ref_press = 1000. * units('mbar')
+    P0 = pot_temp_ref_press = units.Quantity(1000., 'mbar')
     kappa = poisson_exponent = (Rd / Cp_d).to('dimensionless')
     gamma_d = dry_adiabatic_lapse_rate = g / Cp_d
     epsilon = molecular_weight_ratio = (Mw / Md).to('dimensionless')
diff --git a/src/metpy/io/metar.py b/src/metpy/io/metar.py
index 67a60069431..da1680b4a42 100644
--- a/src/metpy/io/metar.py
+++ b/src/metpy/io/metar.py
@@ -158,17 +158,16 @@ def parse_metar_to_dataframe(metar_text, *, year=None, month=None):
     try:
         # Create a field for sea-level pressure and make sure it is a float
         df['air_pressure_at_sea_level'] = float(altimeter_to_sea_level_pressure(
-            df.altimeter.values * units('inHg'),
-            df.elevation.values * units('meters'),
-            df.temperature.values * units('degC')).to('hPa').magnitude)
+            units.Quantity(df.altimeter.values, 'inHg'),
+            units.Quantity(df.elevation.values, 'meters'),
+            units.Quantity(df.temperature.values, 'degC')).to('hPa').magnitude)
     except AttributeError:
         df['air_pressure_at_sea_level'] = [np.nan]
 
     # Use get wind components and assign them to u and v variables
-    df['eastward_wind'], df['northward_wind'] = wind_components((df.wind_speed.values
-                                                                 * units.kts),
-                                                                df.wind_direction.values
-                                                                * units.degree)
+    df['eastward_wind'], df['northward_wind'] = wind_components(
+        units.Quantity(df.wind_speed.values, 'kts'),
+        units.Quantity(df.wind_direction.values, 'degree'))
 
     # Round the altimeter and sea-level pressure values
     df['altimeter'] = df.altimeter.round(2)
@@ -424,7 +423,7 @@ def parse_metar_to_named_tuple(metar_text, station_metadata, year, month):
         if (float(tree.altim.text.strip()[1:5])) > 1100:
             altim = float(tree.altim.text.strip()[1:5]) / 100
         else:
-            altim = (int(tree.altim.text.strip()[1:5]) * units.hPa).to('inHg').magnitude
+            altim = units.Quantity(int(tree.altim.text.strip()[1:5]), 'hPa').to('inHg').m
 
     # Returns a named tuple with all the relevant variables
     return Metar(station_id, lat, lon, elev, date_time, wind_dir, wind_spd,
@@ -617,15 +616,14 @@ def merge(x, key='     '):
 
     # Calculate sea-level pressure from function in metpy.calc
     df['air_pressure_at_sea_level'] = altimeter_to_sea_level_pressure(
-        altim * units('inHg'),
-        elev * units('meters'),
-        temp * units('degC')).to('hPa').magnitude
+        units.Quantity(altim, 'inHg'),
+        units.Quantity(elev, 'meters'),
+        units.Quantity(temp, 'degC')).to('hPa').magnitude
 
     # Use get wind components and assign them to eastward and northward winds
-    df['eastward_wind'], df['northward_wind'] = wind_components((df.wind_speed.values
-                                                                 * units.kts),
-                                                                df.wind_direction.values
-                                                                * units.degree)
+    df['eastward_wind'], df['northward_wind'] = wind_components(
+        units.Quantity(df.wind_speed.values, 'kts'),
+        units.Quantity(df.wind_direction.values, 'degree'))
 
     # Drop duplicate values
     df = df.drop_duplicates(subset=['date_time', 'latitude', 'longitude'], keep='last')
diff --git a/src/metpy/io/station_data.py b/src/metpy/io/station_data.py
index 5d31a1c14ab..a22b2dbf498 100644
--- a/src/metpy/io/station_data.py
+++ b/src/metpy/io/station_data.py
@@ -113,7 +113,8 @@ def _read_airports_file(input_file=None):
     station_map = pd.DataFrame({'id': df.ident.values, 'synop_id': 99999,
                                 'latitude': df.latitude_deg.values,
                                 'longitude': df.longitude_deg.values,
-                                'altitude': ((df.elevation_ft.values * units.ft).to('m')).m,
+                                'altitude': units.Quantity(
+                                    df.elevation_ft.values, 'ft').to('m').m,
                                 'country': df.iso_region.str.split('-', n=1,
                                                                    expand=True)[1].values,
                                 'source': input_file
diff --git a/src/metpy/plots/declarative.py b/src/metpy/plots/declarative.py
index a2b49d41555..686da82bcb5 100644
--- a/src/metpy/plots/declarative.py
+++ b/src/metpy/plots/declarative.py
@@ -1570,7 +1570,8 @@ def _build(self):
             else:
                 field_kwargs['plot_units'] = self.plot_units[0]
             if hasattr(self.data, 'units') and (field_kwargs['plot_units'] is not None):
-                parameter = data[ob_type][subset].values * units(self.data.units[ob_type])
+                parameter = units.Quantity(data[ob_type][subset].values,
+                                           self.data.units[ob_type])
             else:
                 parameter = data[ob_type][subset]
             if field_kwargs['formatter'] is not None:
@@ -1592,10 +1593,10 @@ def _build(self):
             vector_kwargs['color'] = self.vector_field_color
             vector_kwargs['plot_units'] = self.vector_plot_units
             if hasattr(self.data, 'units') and (vector_kwargs['plot_units'] is not None):
-                u = (data[self.vector_field[0]][subset].values
-                     * units(self.data.units[self.vector_field[0]]))
-                v = (data[self.vector_field[1]][subset].values
-                     * units(self.data.units[self.vector_field[1]]))
+                u = units.Quantity(data[self.vector_field[0]][subset].values,
+                                   self.data.units[self.vector_field[0]])
+                v = units.Quantity(data[self.vector_field[1]][subset].values,
+                                   self.data.units[self.vector_field[1]])
             else:
                 vector_kwargs.pop('plot_units')
                 u = data[self.vector_field[0]][subset]
diff --git a/src/metpy/plots/skewt.py b/src/metpy/plots/skewt.py
index d08fbacb45a..304d6dd0e0d 100644
--- a/src/metpy/plots/skewt.py
+++ b/src/metpy/plots/skewt.py
@@ -485,14 +485,15 @@ def plot_dry_adiabats(self, t0=None, pressure=None, **kwargs):
         # Determine set of starting temps if necessary
         if t0 is None:
             xmin, xmax = self.ax.get_xlim()
-            t0 = np.arange(xmin, xmax + 1, 10) * units.degC
+            t0 = units.Quantity(np.arange(xmin, xmax + 1, 10), 'degC')
 
         # Get pressure levels based on ylims if necessary
         if pressure is None:
-            pressure = np.linspace(*self.ax.get_ylim()) * units.mbar
+            pressure = units.Quantity(np.linspace(*self.ax.get_ylim()), 'mbar')
 
         # Assemble into data for plotting
-        t = dry_lapse(pressure, t0[:, np.newaxis], 1000. * units.mbar).to(units.degC)
+        t = dry_lapse(pressure, t0[:, np.newaxis],
+                      units.Quantity(1000., 'mbar')).to(units.degC)
         linedata = [np.vstack((ti.m, pressure.m)).T for ti in t]
 
         # Add to plot
@@ -542,15 +543,16 @@ def plot_moist_adiabats(self, t0=None, pressure=None, **kwargs):
         # Determine set of starting temps if necessary
         if t0 is None:
             xmin, xmax = self.ax.get_xlim()
-            t0 = np.concatenate((np.arange(xmin, 0, 10),
-                                 np.arange(0, xmax + 1, 5))) * units.degC
+            t0 = units.Quantity(np.concatenate((np.arange(xmin, 0, 10),
+                                                np.arange(0, xmax + 1, 5))), 'degC')
 
         # Get pressure levels based on ylims if necessary
         if pressure is None:
-            pressure = np.linspace(*self.ax.get_ylim()) * units.mbar
+            pressure = units.Quantity(np.linspace(*self.ax.get_ylim()), 'mbar')
 
         # Assemble into data for plotting
-        t = moist_lapse(pressure, t0[:, np.newaxis], 1000. * units.mbar).to(units.degC)
+        t = moist_lapse(pressure, t0[:, np.newaxis],
+                        units.Quantity(1000., 'mbar')).to(units.degC)
         linedata = [np.vstack((ti.m, pressure.m)).T for ti in t]
 
         # Add to plot
@@ -600,7 +602,7 @@ def plot_mixing_lines(self, mixing_ratio=None, pressure=None, **kwargs):
 
         # Set pressure range if necessary
         if pressure is None:
-            pressure = np.linspace(600, max(self.ax.get_ylim())) * units.mbar
+            pressure = units.Quantity(np.linspace(600, max(self.ax.get_ylim())), 'mbar')
 
         # Assemble data for plotting
         td = dewpoint(vapor_pressure(pressure, mixing_ratio))
@@ -930,10 +932,10 @@ def plot_colormapped(self, u, v, c, intervals=None, colors=None, **kwargs):
             if intervals.dimensionality == {'[length]': 1.0}:
 
                 # Find any intervals not in the data and interpolate them
-                interpolation_heights = [bound.m for bound in intervals if bound not in c]
-                interpolation_heights = np.array(interpolation_heights) * intervals.units
-                interpolation_heights = (np.sort(interpolation_heights.magnitude)
-                                         * interpolation_heights.units)
+                interpolation_heights = np.array([bound.m for bound in intervals
+                                                  if bound not in c])
+                interpolation_heights = units.Quantity(np.sort(interpolation_heights),
+                                                       intervals.units)
                 (interpolated_heights, interpolated_u,
                  interpolated_v) = interpolate_1d(interpolation_heights, c, c, u, v)
 
@@ -952,7 +954,7 @@ def plot_colormapped(self, u, v, c, intervals=None, colors=None, **kwargs):
                 intervals = intervals.to_base_units()
             # If segmenting by anything else, do not interpolate, just use the data
             else:
-                intervals = np.asarray(intervals) * intervals.units
+                intervals = units.Quantity(np.asarray(intervals), intervals.units)
 
             norm = mcolors.BoundaryNorm(intervals.magnitude, cmap.N)
             cmap.set_over('none')
diff --git a/src/metpy/testing.py b/src/metpy/testing.py
index d1556e35411..ec56a1f2663 100644
--- a/src/metpy/testing.py
+++ b/src/metpy/testing.py
@@ -84,12 +84,12 @@ def to_float(s):
 
     p, z, t, td, direc, spd = np.array(arr_data).T
 
-    p = p * units.hPa
-    z = z * units.meters
-    t = t * units.degC
-    td = td * units.degC
-    direc = direc * units.degrees
-    spd = spd * units.knots
+    p = units.Quantity(p, 'hPa')
+    z = units.Quantity(z, 'meters')
+    t = units.Quantity(t, 'degC')
+    td = units.Quantity(td, 'degC')
+    direc = units.Quantity(direc, 'degrees')
+    spd = units.Quantity(spd, 'knots')
 
     u, v = wind_components(spd, direc)
 
@@ -158,9 +158,9 @@ def check_mask(actual, desired):
     np.testing.assert_array_equal(actual_mask, desired_mask)
 
 
-def assert_nan(value, units):
+def assert_nan(value, value_units):
     """Check for nan with proper units."""
-    value, _ = check_and_drop_units(value, np.nan * units)
+    value, _ = check_and_drop_units(value, units.Quantity(np.nan, value_units))
     assert np.isnan(value)
 
 
diff --git a/src/metpy/units.py b/src/metpy/units.py
index 8e172d18da2..5915bbd6259 100644
--- a/src/metpy/units.py
+++ b/src/metpy/units.py
@@ -92,7 +92,7 @@ def pandas_dataframe_to_unit_arrays(df, column_units=None):
     res = {}
     for column in df:
         if column in column_units and column_units[column]:
-            res[column] = df[column].values * units(column_units[column])
+            res[column] = units.Quantity(df[column].values, column_units[column])
         else:
             res[column] = df[column].values
     return res
diff --git a/src/metpy/xarray.py b/src/metpy/xarray.py
index 41e9dc01a0b..29b6409c917 100644
--- a/src/metpy/xarray.py
+++ b/src/metpy/xarray.py
@@ -444,8 +444,8 @@ def coordinates_identical(self, other):
     @property
     def time_deltas(self):
         """Return the time difference of the data in seconds (to microsecond precision)."""
-        return (np.diff(self._data_array.values).astype('timedelta64[us]').astype('int64')
-                / 1e6 * units.s)
+        us_diffs = np.diff(self._data_array.values).astype('timedelta64[us]').astype('int64')
+        return units.Quantity(us_diffs / 1e6, 's')
 
     def find_axis_name(self, axis):
         """Return the name of the axis corresponding to the given identifier.
@@ -790,7 +790,7 @@ def _rebuild_coords(self, var, crs):
                         new_coord_var = coord_var.copy(
                             data=(
                                 coord_var.metpy.unit_array
-                                * (height * units.meter)
+                                * units.Quantity(height, 'meter')
                             ).m_as('meter')
                         )
                         new_coord_var.attrs['units'] = 'meter'
@@ -1411,10 +1411,12 @@ def grid_deltas_from_dataarray(f, kind='default'):
         dy_units = units(y.attrs.get('units'))
 
     # Broadcast to input and attach units
-    dx = dx_var.set_dims(f.dims, shape=[dx_var.sizes[dim] if dim in dx_var.dims else 1
-                                        for dim in f.dims]).data * dx_units
-    dy = dy_var.set_dims(f.dims, shape=[dy_var.sizes[dim] if dim in dy_var.dims else 1
-                                        for dim in f.dims]).data * dy_units
+    dx_var = dx_var.set_dims(f.dims, shape=[dx_var.sizes[dim] if dim in dx_var.dims else 1
+                                            for dim in f.dims])
+    dx = units.Quantity(dx_var.data, dx_units)
+    dy_var = dy_var.set_dims(f.dims, shape=[dy_var.sizes[dim] if dim in dy_var.dims else 1
+                                            for dim in f.dims])
+    dy = units.Quantity(dy_var.data, dy_units)
 
     return dx, dy
 
diff --git a/tests/calc/test_calc_tools.py b/tests/calc/test_calc_tools.py
index d18346556ab..d2876a51ddc 100644
--- a/tests/calc/test_calc_tools.py
+++ b/tests/calc/test_calc_tools.py
@@ -358,6 +358,17 @@ def test_get_layer(pressure, variable, heights, bottom, depth, interp, expected)
     assert_array_almost_equal(y_layer, expected[1], 3)
 
 
+def test_get_layer_masked():
+    """Test get_layer with masked arrays as input."""
+    p = units.Quantity(np.ma.array([1000, 500, 400]), 'hPa')
+    u = units.Quantity(np.arange(3), 'm/s')
+    p_layer, u_layer = get_layer(p, u, depth=units.Quantity(6000, 'm'))
+    true_p_layer = units.Quantity([1000., 500., 464.4742], 'hPa')
+    true_u_layer = units.Quantity([0., 1., 1.3303], 'm/s')
+    assert_array_almost_equal(p_layer, true_p_layer, 4)
+    assert_array_almost_equal(u_layer, true_u_layer, 4)
+
+
 def test_greater_or_close():
     """Test floating point greater or close to."""
     x = np.array([0.0, 1.0, 1.49999, 1.5, 1.5000, 1.7])
diff --git a/tools/flake8-metpy/flake8_metpy.py b/tools/flake8-metpy/flake8_metpy.py
new file mode 100644
index 00000000000..c9a7337c87e
--- /dev/null
+++ b/tools/flake8-metpy/flake8_metpy.py
@@ -0,0 +1,72 @@
+# Copyright (c) 2021 MetPy Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+"""
+Custom flake8 plugin to catch MetPy-specific bad style/practice.
+
+Currently this only looks for multiplication or division by units, since that can break
+masked arrays and is slower than calling ``Quantity()``.
+"""
+
+import ast
+from collections import namedtuple
+
+Error = namedtuple('Error', 'lineno col code')
+
+
+class MetPyVisitor(ast.NodeVisitor):
+    """Visit nodes of the AST looking for violations."""
+
+    def __init__(self):
+        """Initialize the visitor."""
+        self.errors = []
+
+    @staticmethod
+    def _is_unit(node):
+        """Check whether a node should be considered to represent "units"."""
+        # Looking for a .units attribute, a units.myunit, or a call to units()
+        is_units_attr = isinstance(node, ast.Attribute) and node.attr == 'units'
+        is_reg_attr = (isinstance(node, ast.Attribute)
+                       and isinstance(node.value, ast.Name) and node.value.id == 'units')
+        is_reg_call = (isinstance(node, ast.Call)
+                       and isinstance(node.func, ast.Name) and node.func.id == 'units')
+        is_unit_alias = isinstance(node, ast.Name) and 'unit' in node.id
+
+        return is_units_attr or is_reg_attr or is_reg_call or is_unit_alias
+
+    def visit_BinOp(self, node):
+        """Visit binary operations."""
+        # Check whether this is multiplying or dividing by units
+        if (isinstance(node.op, (ast.Mult, ast.Div))
+                and (self._is_unit(node.right) or self._is_unit(node.left))):
+            self.error(node.lineno, node.col_offset, 1)
+
+        super().generic_visit(node)
+
+    def error(self, lineno, col, code):
+        """Add an error to our output list."""
+        self.errors.append(Error(lineno, col, code))
+
+
+class MetPyChecker:
+    """Flake8 plugin class to check MetPy style/best practice."""
+
+    name = __name__
+    version = '1.0'
+
+    def __init__(self, tree):
+        """Initialize the plugin."""
+        self.tree = tree
+
+    def run(self):
+        """Run the plugin and yield errors."""
+        visitor = MetPyVisitor()
+        visitor.visit(self.tree)
+        for err in visitor.errors:
+            yield self.error(err)
+
+    def error(self, err):
+        """Format errors into Flake8's required format."""
+        return (err.lineno, err.col,
+                f'MPY{err.code:03d}: Multiplying/dividing by units--use units.Quantity()',
+                type(self))
diff --git a/tools/flake8-metpy/test_flake8_metpy.py b/tools/flake8-metpy/test_flake8_metpy.py
new file mode 100644
index 00000000000..ce1d990c2b6
--- /dev/null
+++ b/tools/flake8-metpy/test_flake8_metpy.py
@@ -0,0 +1,28 @@
+# Copyright (c) 2021 MetPy Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+"""Tests for custom Flake8 plugin for MetPy."""
+
+import ast
+
+import pytest
+
+from flake8_metpy import MetPyChecker
+
+
+@pytest.mark.parametrize('source, errs', [
+    ('5 * pressure.units', 1),
+    ('pw = -1. * (np.trapz(w.magnitude, pres.magnitude) * (w.units * pres.units))', 1),
+    ("""def foo():
+    return ret * moist_adiabat_temperatures.units""", 1),
+    ('p_interp = np.sort(np.append(p_interp.m, top_pressure.m)) * pressure.units', 1),
+    ('parameter = data[ob_type][subset].values * units(self.data.units[ob_type])', 1),
+    ('np.nan * pressure.units', 1),
+    ('np.array([1, 2, 3]) * units.m', 1),
+    ('np.arange(4) * units.s', 1),
+    ('np.ma.array([1, 2, 3]) * units.hPa', 1)
+])
+def test_plugin(source, errs):
+    """Test that the flake8 checker works correctly."""
+    checker = MetPyChecker(ast.parse(source))
+    assert len(list(checker.run())) == errs