diff --git a/.github/scripts/pr_mod_file_tests.py b/.github/scripts/pr_mod_file_tests.py
index 6b1ba1375..150febff8 100755
--- a/.github/scripts/pr_mod_file_tests.py
+++ b/.github/scripts/pr_mod_file_tests.py
@@ -140,15 +140,18 @@ def _main_prog():
 
     print("Generating list of modified files...")
 
-    #This should eventually be passed in via a command-line
-    #argument, and include everything inside the "lib" directory -JN:
-    testable_files = {"lib/adf_base.py",
-                      "lib/adf_config.py",
-                      "lib/adf_file_utils.py",
-                      "lib/adf_info.py",
-                      "lib/adf_obs.py",
-                      "lib/adf_web.py",
-                      "lib/adf_diag.py"}
+    # This should eventually be passed in via a command-line
+    # argument, and include everything inside the "lib" directory -JN:
+    testable_files = {
+        "lib/adf_base.py",
+        "lib/adf_config.py",
+        "lib/adf_file_utils.py",
+        "lib/adf_info.py",
+        "lib/adf_obs.py",
+        "lib/adf_units.py",
+        "lib/adf_web.py",
+        "lib/adf_diag.py",
+    }
 
     #+++++++++++++++++++++++
     #Read in input arguments
diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py
index 6caedee6f..227e3d0c5 100644
--- a/lib/adf_dataset.py
+++ b/lib/adf_dataset.py
@@ -5,6 +5,7 @@
 import xarray as xr
 
 import adf_utils as utils
+from adf_units import units_equivalent
 warnings.formatwarning = utils.my_formatwarning
 
 # "reference data"
@@ -294,6 +295,16 @@ def load_climo_ds(self, case, variablename):
         ds = self.load_dataset(fils)
         if ds is None:
             return None
+        if (scale_factor != 1 or add_offset != 0) and self.already_converted(
+            ds[variablename].attrs, variablename
+        ):
+            self.adf.debug_log(
+                f"\t    INFO: '{variablename}' climo file is already"
+                " in converted units; not converting again."
+            )
+            scale_factor = 1
+            add_offset = 0
+        # End if
         # xarray arithmetic drops attrs, so carry them across by hand -- otherwise
         # the regridded files lose 'units' and the plotting scripts KeyError on it.
         attrs = ds[variablename].attrs.copy()
@@ -368,10 +379,10 @@ def load_reference_climo_ds(self, case, variablename, apply_scaling=True):
         if ds is None:
             return None
         vname = self.ref_var_nam[variablename]  # name of variable in the reference data
-        # Check if already transformed (via attribute or units)
-        new_unit = self.adf.variable_defaults.get(variablename, {}).get('new_unit')
-        unit_match = new_unit is not None and ds[vname].attrs.get('units') == new_unit
-        if ds[vname].attrs.get('transformed', False) or unit_match:
+        # Check if already transformed.  The units comparison was a literal
+        # string match, which almost never fired: an observation file says
+        # "W/m2" where the defaults say "Wm$^{-2}$".
+        if self.already_converted(ds[vname].attrs, variablename):
             apply_scaling = False
         if not apply_scaling:
             add_offset = 0
@@ -381,6 +392,13 @@ def load_reference_climo_ds(self, case, variablename, apply_scaling=True):
         ds[vname] = ds[vname] * scale_factor + add_offset
         ds[vname].attrs = attrs
         if scale_factor != 1 or add_offset != 0:
+            # Rename the units as load_climo_ds does for the test cases.  Without
+            # this the reference kept the units it arrived with while holding
+            # converted values, so the file said 'm/s' over data in mm/day.
+            new_unit = self.adf.variable_defaults.get(variablename, {}).get("new_unit")
+            if new_unit:
+                ds[vname].attrs["units"] = new_unit
+            # End if
             # int, not bool: netCDF4 cannot store a Python bool as an attribute
             ds[vname].attrs['transformed'] = 1
         return ds
@@ -394,7 +412,13 @@ def load_reference_climo_da(self, case, variablename, apply_scaling=True):
             scale_factor = 1
         else:
             add_offset, scale_factor = self.get_value_converters(case, variablename)
-        return self.load_da(fils, vname, add_offset=add_offset, scale_factor=scale_factor)
+        return self.load_da(
+            fils,
+            vname,
+            field=variablename,
+            add_offset=add_offset,
+            scale_factor=scale_factor,
+        )
 
     def get_reference_climo_file(self, var):
         """Return a list of files to be used as reference (aka baseline) for variable var."""
@@ -510,7 +534,15 @@ def load_reference_regrid_da(self, case, field, apply_scaling=None):
         file_field = self.ref_var_nam[field] if self.adf.compare_obs else field
         add_offset, scale_factor = self._regrid_converters(fils, file_field, case, field,
                                                            apply_scaling)
-        return self.load_da(fils, file_field, add_offset=add_offset, scale_factor=scale_factor)
+        # field, not file_field: an observation file names the variable its own
+        # way, and the variable defaults are keyed by the ADF name
+        return self.load_da(
+            fils,
+            file_field,
+            field=field,
+            add_offset=add_offset,
+            scale_factor=scale_factor,
+        )
 
     def _regrid_converters(self, fils, file_field, case, field, apply_scaling):
         """Return the (add_offset, scale_factor) to use for a regridded file.
@@ -532,7 +564,7 @@ def _regrid_converters(self, fils, file_field, case, field, apply_scaling):
         if apply_scaling or (scale_factor == 1 and add_offset == 0):
             return add_offset, scale_factor
         ds = self.load_dataset(fils)
-        if ds is not None and ds[file_field].attrs.get('transformed', 0):
+        if ds is not None and self.already_converted(ds[file_field].attrs, field):
             return 0, 1
         return add_offset, scale_factor
 
@@ -569,10 +601,101 @@ def load_dataset(self, fils, use_time_bounds=False):
         # End if
         return ds
 
-    def load_da(self, fils, variablename, use_time_bounds=False, **kwargs):
+    def already_converted(self, attrs, field):
+        """Report whether the conversion for `field` has already been applied.
+
+        Parameters
+        ----------
+        attrs : dict
+            attributes of the variable as it was read from the file
+        field : str
+            ADF name of the variable, i.e. the key into the variable defaults
+
+        Returns
+        -------
+        bool
+            ``True`` when the file already holds converted values.
+
+        Notes
+        -----
+        Two things say so.  A file the ADF wrote carries ``transformed``,
+        stamped by whichever load applied the conversion.  A file it did not
+        write -- an observation file, or a regridded file from an older ADF --
+        carries no stamp, and then the only evidence is that its units are
+        already the units the defaults are converting to.  That comparison is
+        made with :func:`units_equivalent`, because the two strings come from
+        different hands: CAM writes ``W/m2`` where the defaults say
+        ``Wm$^{-2}$``, and comparing them literally answers the wrong
+        question.
+        """
+        if attrs.get("transformed", False):
+            return True
+        # End if
+        new_unit = self.adf.variable_defaults.get(field, {}).get("new_unit")
+        return units_equivalent(attrs.get("units"), new_unit)
+
+    def apply_conversion(self, data, field, case=None):
+        """Apply the variable-defaults conversion to `data`, once.
+
+        For a script that opens a file itself rather than through one of the
+        loaders here -- the TEM files and the vector plots do -- so that it
+        makes the same decision they do.
+
+        Parameters
+        ----------
+        data : xarray.DataArray
+            the variable as it was read from the file
+        field : str
+            ADF name of the variable, i.e. the key into the variable defaults
+        case : str, optional
+            case the data belongs to, which decides whether the observation
+            converters are used.  Defaults to `field`'s test-case converters.
+
+        Returns
+        -------
+        xarray.DataArray
+            The converted data, carrying its attributes, with ``units`` set
+            to the defaults' ``new_unit`` and ``transformed`` stamped on when
+            a conversion was applied.  Returned unchanged when the file
+            already holds converted values.
+        """
+        if data is None:
+            return None
+        # End if
+        add_offset, scale_factor = self.get_value_converters(
+            case if case is not None else field, field
+        )
+        if scale_factor == 1 and add_offset == 0:
+            return data
+        # End if
+        if self.already_converted(data.attrs, field):
+            dmsg = f"\t    INFO: '{field}' is already in converted units in the"
+            dmsg += " file, so the conversion is not applied again."
+            self.adf.debug_log(dmsg)
+            return data
+        # End if
+        attrs = data.attrs.copy()
+        converted = data * scale_factor + add_offset
+        converted.attrs = attrs
+        new_unit = self.adf.variable_defaults.get(field, {}).get("new_unit")
+        if new_unit:
+            converted.attrs["units"] = new_unit
+        # End if
+        # int, not bool: netCDF4 cannot store a Python bool as an attribute
+        converted.attrs["transformed"] = 1
+        return converted
+
+    def load_da(self, fils, variablename, use_time_bounds=False, field=None, **kwargs):
         """Return xarray DataArray from file(s) w/ optional scale factor, offset, new units.
 
         `use_time_bounds` is passed to `load_dataset`; see there.
+
+        `field` is the ADF name of the variable when it differs from its name
+        in the file, which is the case for observations.  The variable
+        defaults are keyed by the ADF name.
+
+        A conversion that has already been applied to the file is not applied
+        again; see :meth:`already_converted`.
         """
         ds = self.load_dataset(fils, use_time_bounds=use_time_bounds)
         if ds is None:
@@ -581,12 +704,25 @@ def load_da(self, fils, variablename, use_time_bounds=False, **kwargs):
         da = ds[variablename].squeeze()
         scale_factor = kwargs.get('scale_factor', 1)
         add_offset = kwargs.get('add_offset', 0)
+
+        if (scale_factor != 1 or add_offset != 0) and self.already_converted(
+            da.attrs, field if field is not None else variablename
+        ):
+            dmsg = f"\t    INFO: '{variablename}' is already in converted units"
+            dmsg += " in the file, so the conversion is not applied again."
+            self.adf.debug_log(dmsg)
+            scale_factor = 1
+            add_offset = 0
+        # End if
+
         attrs = da.attrs.copy()
         da = da * scale_factor + add_offset
         da.attrs = attrs
 
         if scale_factor != 1 or add_offset != 0:
-            new_unit = self.adf.variable_defaults.get(variablename, {}).get("new_unit")
+            new_unit = self.adf.variable_defaults.get(
+                field if field is not None else variablename, {}
+            ).get("new_unit")
             if new_unit:
                 da.attrs['units'] = new_unit
             # Stamp on any conversion, not only one that renames the units --
diff --git a/lib/adf_units.py b/lib/adf_units.py
new file mode 100644
index 000000000..46308d474
--- /dev/null
+++ b/lib/adf_units.py
@@ -0,0 +1,236 @@
+"""
+Unit string handling.
+
+Kept apart from adf_utils so it can be unit tested without the scientific
+stack: the ADF unit test workflow installs only PyYAML and pytest, so anything
+that imports xarray/geocat at module level cannot be exercised in CI.  This
+module imports nothing but re.
+
+Functions
+---------
+normalize_units(units)
+    Reduce a unit string to a canonical form for comparison.
+units_equivalent(first, second)
+    Report whether two unit strings mean the same thing.
+
+Notes
+-----
+The ADF compares units to decide whether the conversion named in the variable
+defaults has already been applied to a file.  The two strings being compared
+come from different places -- one written by CAM, the other by whoever edited
+`adf_variable_defaults.yaml` -- so they agree on the physics far more often
+than they agree character for character.  `W/m2`, `W m-2`, `Wm^-2` and the
+LaTeX `W m$^{-2}$` are one unit written four ways, and the shipped defaults
+alone contain both `Wm$^{-2}$` and `W m$^{-2}$`.  Comparing the raw strings
+answers "were these typed the same way", which is not the question being
+asked, and getting it wrong scales the data twice.
+
+Rendering a unit for somewhere with no LaTeX renderer, such as a table cell,
+is `adf_utils.plain_text_units`; this module only compares.
+"""
+
+import re
+
+# Unit names that mean the same thing.  Keys and values are compared after
+# normalization, so only real synonyms belong here, not spelling variants:
+_ALIASES = {
+    "fraction": "1",
+    "frac": "1",
+    "unitless": "1",
+    "dimensionless": "1",
+    "none": "1",
+    "-": "1",
+    "[-]": "1",
+    "percent": "%",
+    "pct": "%",
+    "degrees_kelvin": "k",
+    "deg_k": "k",
+    "kelvin": "k",
+    "degrees_celsius": "degc",
+    "celsius": "degc",
+    "deg_c": "degc",
+    "degrees_east": "degrees_east",
+    "meters": "m",
+    "meter": "m",
+    "metre": "m",
+    "seconds": "s",
+    "second": "s",
+    "sec": "s",
+    "days": "d",
+    "day": "d",
+    "grams": "g",
+    "gram": "g",
+    "micron": "um",
+    "microns": "um",
+}
+
+# Factors written without a space between them.  Case matters and is still
+# intact here, which is the point: "Nm" is newton-metre while "nm" is
+# nanometre, and lower-casing first makes them the same string.  Only forms
+# that actually occur are listed -- a general splitter guessing where to cut a
+# run of letters reads "Sv" as siemens-volt and "nm" as newton-metre, the
+# second of which reports two different units as equal.  Note that "ms" is
+# read as metre-second: no CAM field is reported in milliseconds, and
+# "ms$^{-1}$" for wind speed is in the shipped variable defaults today.
+_RUN_TOGETHER = {
+    "Wm": "W m",
+    "Nm": "N m",
+    "Jm": "J m",
+    "Km": "K m",
+    "kgm": "kg m",
+    "gm": "g m",
+    "ms": "m s",
+    "Pas": "Pa s",
+}
+
+# LaTeX and unicode fragments that carry no meaning for a comparison:
+_LATEX = (
+    (r"\mathrm", ""),
+    (r"\text", ""),
+    (r"\mu", "u"),
+    (r"\,", " "),
+    (r"\;", " "),
+    (r"\ ", " "),
+    ("µ", "u"),  # micro sign
+    ("μ", "u"),  # greek small letter mu
+    ("·", " "),  # middle dot, used as a multiplication sign
+    ("**", "^"),
+)
+
+# Superscript digits, which appear in units copied out of documents:
+_SUPERSCRIPTS = str.maketrans(
+    {
+        "⁰": "0",
+        "¹": "1",
+        "²": "2",
+        "³": "3",
+        "⁴": "4",
+        "⁵": "5",
+        "⁶": "6",
+        "⁷": "7",
+        "⁸": "8",
+        "⁹": "9",
+        "⁻": "-",
+    }
+)
+
+
+def _strip_markup(units):
+    """Return `units` with LaTeX, superscripts and stray braces removed."""
+    text = units.translate(_SUPERSCRIPTS)
+    for old, new in _LATEX:
+        text = text.replace(old, new)
+    # A LaTeX exponent is written ^{-2}; the braces are noise once the
+    # superscript digits are gone:
+    text = re.sub(r"\^\s*\{([^}]*)\}", r"^\1", text)
+    text = text.replace("$", "").replace("{", "").replace("}", "")
+    # Split the run-together factors before anything lower-cases the string,
+    # because which factors they are depends on their case:
+    for joined, apart in _RUN_TOGETHER.items():
+        text = re.sub(rf"(?<![A-Za-z]){joined}(?=[-^0-9\s]|$)", apart, text)
+    # End for
+    return text.strip()
+
+
+def _tokenize(text, sign):
+    """Return (name, exponent) pairs for one side of a solidus."""
+    tokens = []
+    for chunk in re.split(r"[\s*.]+", text):
+        if not chunk:
+            continue
+        # A trailing exponent, written either "m2", "m^2", "m-2" or "m^-2":
+        match = re.fullmatch(r"([a-z%_]+)\^?([+-]?\d+)?", chunk)
+        if match is None:
+            # Not something this understands; keep it whole so that two
+            # identical odd strings still compare equal:
+            tokens.append((chunk, sign))
+            continue
+        name, exponent = match.group(1), match.group(2)
+        tokens.append((name, sign * int(exponent if exponent else 1)))
+    return tokens
+
+
+def normalize_units(units):
+    """
+    Reduce a unit string to a canonical form for comparison.
+
+    Parameters
+    ----------
+    units : str or None
+        a unit string as it appears in a file or in the variable defaults
+
+    Returns
+    -------
+    str
+        A canonical form: lower case, no LaTeX, every factor written
+        ``name^exponent`` and sorted, so that the many spellings of one unit
+        reduce to a single string.  An empty string for ``None`` or for a
+        string that holds nothing.
+
+    Examples
+    --------
+    ``W/m2``, ``W m-2``, ``Wm^-2`` and ``W m$^{-2}$`` all give ``m^-2 w^1``.
+    """
+    if units is None:
+        return ""
+    text = _strip_markup(str(units)).lower()
+    if not text:
+        return ""
+    text = _ALIASES.get(text, text)
+    # Split on the solidus: everything after the first one is a denominator.
+    # "a/b/c" is read the way it is written, left to right.
+    parts = text.split("/")
+    tokens = _tokenize(parts[0], 1)
+    for part in parts[1:]:
+        tokens += _tokenize(part, -1)
+
+    # Fold the aliases once more, now that the factors are separated, so that
+    # "m/s" and "meters/second" agree:
+    folded = {}
+    for name, exponent in tokens:
+        name = _ALIASES.get(name, name)
+        folded[name] = folded.get(name, 0) + exponent
+    # A factor that cancels out carries no dimension, but it does carry an
+    # identity: "kg/kg" and "kg kg-1" are the same thing, and neither is
+    # "mol/mol".  A mass mixing ratio and a volume mixing ratio are both
+    # dimensionless and are not the same number, so what cancelled is kept.
+    remaining = {
+        name: exponent
+        for name, exponent in folded.items()
+        if exponent != 0 and name != "1"
+    }
+    if not remaining:
+        cancelled = sorted(name for name in folded if name != "1")
+        return f"1[{' '.join(cancelled)}]" if cancelled else "1"
+    # End if
+    return " ".join(
+        f"{name}^{exponent}" for name, exponent in sorted(remaining.items())
+    )
+
+
+def units_equivalent(first, second):
+    """
+    Report whether two unit strings mean the same thing.
+
+    Parameters
+    ----------
+    first, second : str or None
+        the unit strings to compare
+
+    Returns
+    -------
+    bool
+        ``True`` when the two describe the same unit, however they are
+        spelled.  Two strings that are both empty or ``None`` are not
+        equivalent to anything, including each other: nothing is known about
+        a variable that does not say what its units are.
+    """
+    left = normalize_units(first)
+    right = normalize_units(second)
+    if not left or not right:
+        return False
+    return left == right
+
+
+##############
+# END OF FILE
diff --git a/lib/adf_utils.py b/lib/adf_utils.py
index 858124578..7a01c50f3 100644
--- a/lib/adf_utils.py
+++ b/lib/adf_utils.py
@@ -7,6 +7,9 @@
 ts_file_span(),
 as_hist_str_list(), pick_hist_str()
     re-exported from adf_file_utils; time series file discovery
+normalize_units(), units_equivalent()
+    re-exported from adf_units; compare unit strings that are spelled
+    differently but mean the same thing
 plain_text_units()
     render a unit string from the variable defaults as plain text
 use_time_bounds_midpoint()
@@ -63,6 +66,7 @@
 #tested without importing the scientific stack (see adf_file_utils).  Re-export
 #here so `utils.find_ts_files(...)` keeps working for every existing caller:
 # pylint: disable=unused-import
+from adf_units import normalize_units, units_equivalent
 from adf_file_utils import (
     as_hist_str_list,
     describe_dir_problem,
diff --git a/lib/plotting_functions.py b/lib/plotting_functions.py
index 86911759e..88b9987e2 100644
--- a/lib/plotting_functions.py
+++ b/lib/plotting_functions.py
@@ -1210,8 +1210,12 @@ def plot_meridional_mean_and_save(wks, case_nickname, base_nickname,
         latweight = np.cos(np.radians(adata.lat))
         adata = adata.weighted(latweight).mean(dim='lat', keep_attrs=True)
     if 'time' in bdata.dims:
-        adata = bdata.mean(dim='time', keep_attrs=True)
-    if 'lat' in bdata.dims:
+        # bdata, not adata: this took the reference's time mean and put it in
+        # the test case, throwing away the latitude weighting just applied to
+        # it.  Latent, because seasonal_mean(is_climo=True) has already
+        # removed 'time' by the time the ADF calls this.
+        bdata = bdata.mean(dim="time", keep_attrs=True)
+    if "lat" in bdata.dims:
         latweight = np.cos(np.radians(bdata.lat))
         bdata = bdata.weighted(latweight).mean(dim='lat', keep_attrs=True)
     # If there are other dimensions, they are still going to be there:
@@ -1511,4 +1515,4 @@ def square_contour_difference(fld1, fld2, **kwargs):
     return fig
 
 #####################
-#END HELPER FUNCTIONS
\ No newline at end of file
+#END HELPER FUNCTIONS
diff --git a/lib/plotting_utils.py b/lib/plotting_utils.py
index e83836429..18e08ceb9 100644
--- a/lib/plotting_utils.py
+++ b/lib/plotting_utils.py
@@ -1,8 +1,9 @@
-"""                                                                    .
-Generic plotting helper functions
+"""Generic plotting helper functions.
 
 Functions
 ---------
+first_complete_plot_set(layouts)
+    pick the set of plots that is already on disk, if there is one
 use_this_norm()
     switches matplotlib color normalization method
 get_difference_colors(values)
@@ -27,6 +28,8 @@
 """
 
 #import statements:
+from pathlib import Path
+
 import numpy as np
 import xarray as xr
 import matplotlib as mpl
@@ -72,6 +75,40 @@ def load_dataset(fils):
 
 #######
 
+
+def first_complete_plot_set(layouts):
+    """
+    Return the first set of plots whose files all exist.
+
+    A plotting script can usually tell what it would produce from the
+    variable defaults and the configured seasons and pressure levels, without
+    opening any data -- except that a 2-D variable and a 3-D one are named
+    differently, and which one it has is a property of the data.  Both
+    spellings can be checked instead: if either set is complete on disk, the
+    script knows both that nothing needs drawing and which kind of variable it
+    was looking at, and never has to open the file.
+
+    Parameters
+    ----------
+    layouts : list
+        candidate sets, most specific first.  Each is a list of entries whose
+        first element is the plot's path; the rest is whatever the caller
+        needs to register the plot on the website.
+
+    Returns
+    -------
+    list
+        The first complete set, or an empty list when none is complete, which
+        means the data has to be opened to find out what is missing.
+    """
+    for entries in layouts:
+        if entries and all(Path(entry[0]).is_file() for entry in entries):
+            return entries
+        # End if
+    # End for
+    return []
+
+
 def use_this_norm():
     """Just use the right normalization; avoids a deprecation warning."""
 
diff --git a/lib/test/unit_tests/test_adf_units.py b/lib/test/unit_tests/test_adf_units.py
new file mode 100644
index 000000000..982e21d43
--- /dev/null
+++ b/lib/test/unit_tests/test_adf_units.py
@@ -0,0 +1,215 @@
+"""
+Collection of python unit tests for the "adf_units" unit string helpers.
+"""
+
+# +++++++++++++++++++++++
+# Import required modules
+# +++++++++++++++++++++++
+
+import unittest
+import sys
+import os
+import os.path
+
+# Set relevant path variables:
+_CURRDIR = os.path.abspath(os.path.dirname(__file__))
+_ADF_LIB_DIR = os.path.join(_CURRDIR, os.pardir, os.pardir)
+
+# Add ADF "lib" directory to python path:
+sys.path.append(_ADF_LIB_DIR)
+
+# adf_units imports nothing but re, so these run in CI, where only PyYAML and
+# pytest are installed:
+from adf_units import normalize_units, units_equivalent
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# Main adf_units testing routine, used when script is run directly
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+class AdfUnitsTestRoutine(unittest.TestCase):
+    """
+    Unit tests for comparing unit strings, which the ADF does to decide
+    whether a conversion has already been applied to a file.
+    """
+
+    def test_one_unit_written_many_ways(self):
+        """
+        Check that the spellings of watts per square metre all agree.
+
+        These are not hypothetical: CAM writes 'W/m2', and
+        adf_variable_defaults.yaml contains both 'Wm$^{-2}$' and
+        'W m$^{-2}$'.
+        """
+
+        spellings = [
+            "W/m2",
+            "W m-2",
+            "Wm-2",
+            "Wm^-2",
+            "W m^-2",
+            "W/m^2",
+            "W m**-2",
+            "W m$^{-2}$",
+            "Wm$^{-2}$",
+            "W m⁻²",
+        ]
+
+        canonical = {normalize_units(u) for u in spellings}
+
+        self.assertEqual(
+            len(canonical), 1, f"expected one canonical form, got {canonical}"
+        )
+        for unit in spellings:
+            self.assertTrue(units_equivalent(unit, "W/m2"), unit)
+
+    def test_defaults_spellings_match_file_spellings(self):
+        """
+        Check the pairs that actually occur: a unit as CAM writes it against
+        the same unit as the variable defaults write it.
+        """
+
+        pairs = [
+            ("W/m2", "Wm$^{-2}$"),
+            ("m/s", "ms$^{-1}$"),
+            ("mm/day", "mm d$^{-1}$"),
+            ("ug/m3", "$\\mu$g/m3"),
+            ("mol/mol", "mol mol$^{-1}$"),
+            ("fraction", "Fraction"),
+            ("%", "Percent"),
+        ]
+
+        for from_file, from_defaults in pairs:
+            with self.subTest(units=from_file):
+                self.assertTrue(units_equivalent(from_file, from_defaults))
+
+    def test_different_units_stay_different(self):
+        """
+        Check that units which differ are not reported as equivalent.
+
+        This is the direction that matters for correctness: a false match
+        means a conversion is skipped and the data is plotted in the wrong
+        units.
+        """
+
+        pairs = [
+            ("W/m2", "W/m3"),
+            ("m/s", "s/m"),
+            ("K", "W/m2"),
+            ("mm/day", "mm/s"),
+            ("Pa", "hPa"),
+            ("kg/m2", "kg/m3"),
+            ("ppbv", "ppmv"),
+            ("degrees_east", "degrees_north"),
+        ]
+
+        for first, second in pairs:
+            with self.subTest(units=(first, second)):
+                self.assertFalse(units_equivalent(first, second))
+
+    def test_dimensionless_forms_agree(self):
+        """
+        Check that the words for "no units" all mean the same thing.
+
+        A ratio that cancels is dimensionless too, but keeps what cancelled --
+        see test_dimensionless_ratios_keep_what_cancelled.
+        """
+
+        for unit in ["fraction", "Fraction", "1", "unitless", "none"]:
+            with self.subTest(units=unit):
+                self.assertEqual(normalize_units(unit), "1")
+
+    def test_case_decides_a_run_together_factor(self):
+        """
+        Check that "Nm" and "nm" are not read as the same unit.
+
+        Newton-metre and nanometre differ only in case, so the run-together
+        factors have to be split before anything lower-cases the string.
+        Reporting these two as equal would let a conversion be skipped on a
+        file whose units are nothing like the ones being converted to.
+        """
+
+        self.assertTrue(units_equivalent("N/m2", "Nm-2"))
+        self.assertFalse(units_equivalent("Nm-2", "nm-2"))
+        self.assertNotEqual(normalize_units("nm"), normalize_units("N m"))
+
+    def test_unit_names_are_not_chopped_into_letters(self):
+        """
+        Check that a unit name is left whole rather than read as a product of
+        whatever symbols happen to spell it.
+
+        "Sv" is sieverts, not siemens-volt; "cal" is calories.
+        """
+
+        for unit in ["Sv", "cal", "dam", "DU", "molec"]:
+            with self.subTest(units=unit):
+                self.assertEqual(normalize_units(unit), unit.lower() + "^1")
+        self.assertFalse(units_equivalent("Sv", "vs"))
+
+    def test_dimensionless_ratios_keep_what_cancelled(self):
+        """
+        Check that two dimensionless ratios of different things differ.
+
+        A mass mixing ratio and a volume mixing ratio are both dimensionless
+        and are not the same number.
+        """
+
+        self.assertTrue(units_equivalent("kg/kg", "kg kg-1"))
+        self.assertFalse(units_equivalent("kg/kg", "mol/mol"))
+        self.assertFalse(units_equivalent("kg/kg", "m3/m3"))
+
+    def test_units_with_several_factors(self):
+        """
+        Check units written with more than one solidus or exponent.
+        """
+
+        self.assertTrue(units_equivalent("kg/m2/s", "kg m-2 s-1"))
+        self.assertTrue(units_equivalent("W/m2/K", "W m-2 K-1"))
+        self.assertTrue(units_equivalent("1/s", "s-1"))
+        self.assertFalse(units_equivalent("kg/m2/s", "kg/m2"))
+
+    def test_missing_units_match_nothing(self):
+        """
+        Check that an absent unit is not equivalent to anything.
+
+        Nothing is known about a variable that does not say what its units
+        are, so the ADF must not conclude a conversion was already applied.
+        """
+
+        self.assertEqual(normalize_units(None), "")
+        for missing in [None, "", "   "]:
+            with self.subTest(units=missing):
+                self.assertFalse(units_equivalent(missing, "K"))
+                self.assertFalse(units_equivalent("K", missing))
+                self.assertFalse(units_equivalent(missing, missing))
+
+    def test_unfamiliar_unit_matches_itself(self):
+        """
+        Check that a unit this does not understand still compares equal to
+        itself, and unequal to something else.
+        """
+
+        self.assertTrue(units_equivalent("ppbv", "ppbv"))
+        self.assertTrue(units_equivalent("kg m-2 s-1", "kg/m2/s"))
+        self.assertFalse(units_equivalent("frobnicate", "ppbv"))
+
+    def test_hpa_is_not_hecto_pascal_split(self):
+        """
+        Check that a known symbol is not split into smaller ones.
+
+        'hPa' would become hour-pascal if the run-together splitting were
+        applied to a symbol that stands on its own.
+        """
+
+        self.assertEqual(normalize_units("hPa"), "hpa^1")
+        self.assertFalse(units_equivalent("hPa", "Pa"))
+
+
+# ++++++++++++++++++
+
+# Run unit tests if this script is called directly:
+if __name__ == "__main__":
+    unittest.main()
+
+#############
+# End of file
diff --git a/lib/test/unit_tests/test_regrid_scaling.py b/lib/test/unit_tests/test_regrid_scaling.py
index 08f85860c..576f9f9ee 100644
--- a/lib/test/unit_tests/test_regrid_scaling.py
+++ b/lib/test/unit_tests/test_regrid_scaling.py
@@ -25,11 +25,19 @@
 NO_CONVERSION = (0, 1)
 
 
-def _data(stamped, converters=CONVERTED):
+def _data(stamped, converters=CONVERTED, units=None, new_unit=None):
     """An AdfData stand-in holding one variable, with or without the stamp."""
     attrs = {'transformed': 1} if stamped else {}
+    if units is not None:
+        attrs["units"] = units
+    defaults = {"TAUX": {"new_unit": new_unit}} if new_unit else {}
     obj = SimpleNamespace(
         get_value_converters=lambda case, field: converters,
+        adf=SimpleNamespace(variable_defaults=defaults),
+    )
+    # The real decision, so that this pins down what the loaders actually do:
+    obj.already_converted = lambda attrs_, field: AdfData.already_converted(
+        obj, attrs_, field
     )
     obj.reads = 0
 
@@ -63,3 +71,18 @@ def test_variable_without_a_conversion_never_opens_the_file():
 def test_explicit_override_wins_over_the_stamp():
     assert _call(_data(stamped=False), apply_scaling=False) == NO_CONVERSION
     assert _call(_data(stamped=True), apply_scaling=True) == CONVERTED
+
+
+def test_units_already_the_converted_ones_are_not_converted_again():
+    """A file an older ADF wrote carries no stamp, but its units give it away.
+
+    The comparison has to survive the two being spelled differently: CAM
+    writes "W/m2" where the variable defaults say "Wm$^{-2}$".
+    """
+    obj = _data(stamped=False, units="W/m2", new_unit="Wm$^{-2}$")
+    assert _call(obj) == NO_CONVERSION
+
+
+def test_units_that_differ_still_convert():
+    obj = _data(stamped=False, units="m/s", new_unit="Wm$^{-2}$")
+    assert _call(obj) == CONVERTED
diff --git a/scripts/plotting/aod_latlon.py b/scripts/plotting/aod_latlon.py
index 926dc884f..00e055516 100644
--- a/scripts/plotting/aod_latlon.py
+++ b/scripts/plotting/aod_latlon.py
@@ -32,7 +32,18 @@ class AODPlotConfig:
 def aod_latlon(adfobj):
     """Generate AOD comparison plots."""
     config = AODPlotConfig()
-    
+
+    # Every panel is named from the observation source and the season, both
+    # of which are known here.  Settling a re-run first means neither the
+    # observation files nor any case's climatologies are opened, which is
+    # what this script spends its time on:
+    if _register_existing_panels(adfobj, config):
+        print(
+            "\t    INFO: All AOD panels exist. Existing plots added to" " website data."
+        )
+        return
+    # End if
+
     # Load observations
     obs_data = load_observations(adfobj)
     if not obs_data:
@@ -490,4 +501,57 @@ def monthly_to_seasonal(ds, obs=False):
     ds_seasonal['season'] = seasons
     ds_seasonal = ds_seasonal.transpose('lat', 'lon', 'season')
 
-    return ds_seasonal
\ No newline at end of file
+    return ds_seasonal
+
+
+def _register_existing_panels(adfobj, config):
+    """
+    Put the AOD panels that are already on disk onto the website.
+
+    Parameters
+    ----------
+    adfobj : AdfDiag
+        The diagnostics object containing configuration
+    config : AODPlotConfig
+        the observation sources and seasons being plotted
+
+    Returns
+    -------
+    bool
+        ``True`` when every panel exists and ``redo_plot`` is false, so there
+        is nothing left to draw.  ``False`` leaves the plots alone: the run
+        goes on to make them, and registers them as it does.
+    """
+    if adfobj.get_basic_info("redo_plot"):
+        return False
+    # End if
+    file_type = adfobj.read_config_var("diag_basic_info").get("plot_type", "png")
+    plot_dir = Path(adfobj.plot_location[0])
+    panels = []
+    for obs_name in config.obs_sources:
+        for season in config.seasons:
+            label = obs_name.replace(" ", "_")
+            panels.append(
+                (
+                    plot_dir / f"AOD_diff_{label}_{season}_LatLon_Mean.{file_type}",
+                    f"AOD_diff_{label}",
+                    season,
+                )
+            )
+        # End for
+    # End for
+    if not all(path.is_file() for path, _, _ in panels):
+        return False
+    # End if
+    for path, web_name, season in panels:
+        adfobj.add_website_data(
+            path,
+            web_name,
+            None,
+            season=season,
+            multi_case=True,
+            plot_type="LatLon",
+            category="4-Panel AOD Diags",
+        )
+    # End for
+    return True
diff --git a/scripts/plotting/global_latlon_vect_map.py b/scripts/plotting/global_latlon_vect_map.py
index 403a0d8a8..c38004aa4 100644
--- a/scripts/plotting/global_latlon_vect_map.py
+++ b/scripts/plotting/global_latlon_vect_map.py
@@ -185,6 +185,32 @@ def global_latlon_vect_map(adfobj):
         # otherwise defaults to 180
         vres['central_longitude'] = plot_utils.get_central_longitude(adfobj)
 
+        # A complete set of plots can be recognised from the file names alone,
+        # because a 2-D vector and a 3-D one are named differently.  Doing that
+        # first means a re-run never opens the four regridded files a vector
+        # pair needs just to find out which of the two it is looking at.
+        settled = _existing_plot_set(
+            plot_locations, case_names, var_name, seasons, pres_levs, plot_type
+        )
+        if settled and not redo_plot:
+            for path, web_name, case_name, season in settled:
+                adfobj.debug_log(f"'{path}' exists and clobber is false.")
+                adfobj.add_website_data(
+                    path,
+                    web_name,
+                    case_name,
+                    category=web_category,
+                    season=season,
+                    plot_type="LatLon_Vector",
+                )
+            # End for
+            print(
+                f"\t    INFO: All plots exist for {var_name}. "
+                f"Redo is {redo_plot}. Existing plots added to website data."
+            )
+            continue
+        # End if
+
         #Determine observations to compare against:
         if adfobj.compare_obs:
             if var not in adfobj.data.ref_var_nam:
@@ -272,9 +298,12 @@ def global_latlon_vect_map(adfobj):
             uodata = uoclim_ds[data_var[0]].squeeze()  # squeeze in case of degenerate dimensions
             vodata = voclim_ds[data_var[1]].squeeze()  # squeeze in case of degenerate dimensions
 
-            #Convert units if requested (assumes units between model and data are the same):
-            uodata = uodata * vres.get("scale_factor",1) + vres.get("add_offset", 0)
-            vodata = vodata * vres.get("scale_factor",1) + vres.get("add_offset", 0)
+            # Convert units if requested (assumes units between model and data are the same).
+            # Through the ADF's data layer, so that a file already holding
+            # converted values -- which is what the regridding stage writes --
+            # is not scaled a second time:
+            uodata = adfobj.data.apply_conversion(uodata, var)
+            vodata = adfobj.data.apply_conversion(vodata, var_pair)
 
             #Check zonal mean dimensions
             has_lat_ref, has_lev_ref = utils.zm_validate_dims(uodata)
@@ -331,9 +360,9 @@ def global_latlon_vect_map(adfobj):
                 umdata = umclim_ds[var].squeeze()
                 vmdata = vmclim_ds[var_pair].squeeze()
 
-                #Convert units if requested:
-                umdata = umdata * vres.get("scale_factor",1) + vres.get("add_offset", 0)
-                vmdata = vmdata * vres.get("scale_factor",1) + vres.get("add_offset", 0)
+                # Convert units if requested, once -- see above:
+                umdata = adfobj.data.apply_conversion(umdata, var)
+                vmdata = adfobj.data.apply_conversion(vmdata, var_pair)
 
                 #Check dimensions:
                 has_lat, has_lev = utils.zm_validate_dims(umdata)
@@ -346,12 +375,9 @@ def global_latlon_vect_map(adfobj):
                     continue
                 # End if
 
-                # update units
-                # NOTE: looks like our climo files don't have all their metadata
-                uodata.attrs['units'] = vres.get("new_unit", uodata.attrs.get('units', 'none'))
-                vodata.attrs['units'] = vres.get("new_unit", vodata.attrs.get('units', 'none'))
-                umdata.attrs['units'] = vres.get("new_unit", umdata.attrs.get('units', 'none'))
-                vmdata.attrs['units'] = vres.get("new_unit", vmdata.attrs.get('units', 'none'))
+                # Units are set by apply_conversion above, which renames them
+                # only when it actually converts.  Relabelling here regardless
+                # claimed the defaults' units for data still in the file's own.
 
                 #Determine if observations/baseline have the correct dimensions:
                 if has_lev:
@@ -519,5 +545,71 @@ def global_latlon_vect_map(adfobj):
     #Notify user that script has ended:
     print("  ...lat/lon vector maps have been generated successfully.")
 
+
+def _existing_plot_set(
+    plot_locations, case_names, var_name, seasons, pres_levs, plot_type
+):
+    """
+    Return the complete set of vector plots for `var_name`, if one is on disk.
+
+    A 2-D vector is plotted as ``{name}_{season}_LatLon_Vector_Mean`` and a
+    3-D one as ``{name}_{pressure}hpa_{season}_LatLon_Vector_Mean``, so which
+    of the two a vector pair is can be read off the file names.
+
+    Parameters
+    ----------
+    plot_locations : list
+        output plot directory for each case
+    case_names : list
+        names of the test cases
+    var_name : str
+        the vector's name from the variable defaults, e.g. "Wind"
+    seasons : list or dict
+        the seasons being plotted
+    pres_levs : list
+        the configured pressure levels, empty when none are set
+    plot_type : str
+        file extension the plots are written with, e.g. "png"
+
+    Returns
+    -------
+    list
+        ``(path, website name, case, season)`` for every plot, when a
+        complete set exists; an empty list when it does not.
+    """
+    import plotting_utils as plot_utils
+
+    flat = []
+    levelled = []
+    for case_idx, case_name in enumerate(case_names):
+        plot_loc = Path(plot_locations[case_idx])
+        for season in seasons:
+            flat.append(
+                (
+                    plot_loc / f"{var_name}_{season}_LatLon_Vector_Mean.{plot_type}",
+                    var_name,
+                    case_name,
+                    season,
+                )
+            )
+            for lev in pres_levs:
+                levelled.append(
+                    (
+                        plot_loc
+                        / (
+                            f"{var_name}_{lev}hpa_{season}"
+                            f"_LatLon_Vector_Mean.{plot_type}"
+                        ),
+                        f"{var_name}_{lev}hpa",
+                        case_name,
+                        season,
+                    )
+                )
+            # End for
+        # End for
+    # End for
+    return plot_utils.first_complete_plot_set([levelled, flat])
+
+
 ##############
-#END OF SCRIPT
\ No newline at end of file
+#END OF SCRIPT
diff --git a/scripts/plotting/global_mean_timeseries.py b/scripts/plotting/global_mean_timeseries.py
index 5ea5e9731..54ff1c9d3 100644
--- a/scripts/plotting/global_mean_timeseries.py
+++ b/scripts/plotting/global_mean_timeseries.py
@@ -34,6 +34,8 @@ def global_mean_timeseries(adfobj):
     # Gather ADF configurations
     plot_loc = get_plot_loc(adfobj)
     plot_type = adfobj.read_config_var("diag_basic_info").get("plot_type", "png")
+    redo_plot = adfobj.get_basic_info("redo_plot")
+    print(f"\t NOTE: redo_plot is set to {redo_plot}")
     res = adfobj.variable_defaults # will be dict of variable-specific plot preferences
     # or an empty dictionary if use_defaults was not specified in YAML.
 
@@ -42,6 +44,24 @@ def global_mean_timeseries(adfobj):
         #Notify user of variable being plotted:
         print(f"\t - time series plot for {field}")
 
+        # This field's plot is named the same way whatever the data turns out
+        # to hold, so a re-run can settle it here -- before opening the
+        # reference time series and every case's, and before taking the
+        # global means of them, which is the whole cost of this script:
+        plot_name = plot_loc / f"{field}_GlobalMean_ANN_TimeSeries_Mean.{plot_type}"
+        if (not redo_plot) and plot_name.is_file():
+            adfobj.debug_log(f"'{plot_name}' exists and clobber is false.")
+            adfobj.add_website_data(
+                plot_name,
+                f"{field}_GlobalMean",
+                None,
+                season="ANN",
+                multi_case=True,
+                plot_type="TimeSeries",
+            )
+            continue
+        # End if
+
         # Check res for any variable specific options that need to be used BEFORE going to the plot:
         if field in res:
             vres = res[field]
@@ -162,7 +182,6 @@ def global_mean_timeseries(adfobj):
 
         unit = vres.get("new_unit","[-]")
         ax.set_ylabel(getattr(ref_ts_da,"unit", unit)) # add units
-        plot_name = plot_loc / f"{field}_GlobalMean_ANN_TimeSeries_Mean.{plot_type}"
 
         conditional_save(adfobj, plot_name, fig)
 
@@ -313,4 +332,4 @@ def make_plot(case_ts, lens2, label=None, ref_ts_da=None):
 
 
 ##############
-#END OF SCRIPT
\ No newline at end of file
+#END OF SCRIPT
diff --git a/scripts/plotting/meridional_mean.py b/scripts/plotting/meridional_mean.py
index dc47854e2..9fa9df3e9 100644
--- a/scripts/plotting/meridional_mean.py
+++ b/scripts/plotting/meridional_mean.py
@@ -1,76 +1,67 @@
 from pathlib import Path
 import numpy as np
-import xarray as xr
 import plotting_functions as pf
 
 import adf_utils as utils
 import warnings  # use to warn user about missing files.
 warnings.formatwarning = utils.my_formatwarning
 
+
 def meridional_mean(adfobj):
 
     """
-    This script plots meridional averages.
-    Follows the old AMWG convention of plotting 5S to 5N.
-    **Note:** the constraint of 5S to 5N is easily changed;
-    the function that calculates the average can take any range of latitudes.
-    Compare CAM climatologies against
-    other climatological data (observations or baseline runs).
+    Plots meridional average from climatological files (annual and seasonal).
+
+    Follows the old AMWG convention of averaging 5S to 5N.  **Note:** that
+    constraint is easily changed; the function that calculates the average
+    takes any range of latitudes.  Compares CAM climatologies against other
+    climatological data (observations or baseline runs).
+
+    Parameters
+    ----------
+    adfobj : AdfDiag
+        The diagnostics object that contains all the configuration information
+
+    Returns
+    -------
+    None
+        Does not return value, produces files.
+
+    Notes
+    -----
+    Uses AdfData for loading data described by adfobj.
+
+    Directly uses adfobj for the following:
+    plot_var_list, plot_location, climo_yrs, variable_defaults,
+    read_config_var, get_basic_info, add_website_data, debug_log,
+    compare_obs
+
+    Every plot this makes is named the same way whether or not the variable
+    has a `lev` dimension, so a case whose plots are all present can be
+    settled from the file names alone -- see the pre-flight pass below, which
+    is what keeps re-running the ADF cheap.
     """
 
     #Notify user that script has started:
     msg = "\n  Generating meridional mean plots..."
     print(f"{msg}\n  {'-' * (len(msg)-3)}")
 
-    #Extract needed quantities from ADF object:
-    #-----------------------------------------
     var_list = adfobj.plot_var_list
-    model_rgrid_loc = adfobj.get_basic_info("cam_regrid_loc", required=True)
 
     #Special ADF variable which contains the output paths for
     #all generated plots and tables:
     plot_locations = adfobj.plot_location
 
-    #CAM simulation variables (this is always assumed to be a list):
-    case_names = adfobj.get_cam_info("cam_case_name", required=True)
-
     #Grab case years
     syear_cases = adfobj.climo_yrs["syears"]
     eyear_cases = adfobj.climo_yrs["eyears"]
 
-    # CAUTION:
-    # "data" here refers to either obs or a baseline simulation,
-    # Until those are both treated the same (via intake-esm or similar)
-    # we will do a simple check and switch options as needed:
-    if adfobj.get_basic_info("compare_obs"):
-        #Set obs call for observation details for plot titles
-        obs = True
-
-        #Extract variable-obs dictionary:
-        var_obs_dict = adfobj.var_obs_dict
-
-        #If dictionary is empty, then  there are no observations to regrid to,
-        #so quit here:
-        if not var_obs_dict:
-            print("\t No observations found to plot against, so no meridional-mean maps will be generated.")
-            return
-    else:
-        obs = False
-        data_name = adfobj.get_baseline_info("cam_case_name", required=True) # does not get used, is just here as a placemarker
-        data_list = [data_name] # gets used as just the name to search for climo files HAS TO BE LIST
-        data_loc  = model_rgrid_loc #Just use the re-gridded model data path
-    #End if
-
     #Grab baseline years (which may be empty strings if using Obs):
     syear_baseline = adfobj.climo_yrs["syear_baseline"]
     eyear_baseline = adfobj.climo_yrs["eyear_baseline"]
 
-    #Grab all case nickname(s)
-    test_nicknames = adfobj.case_nicknames["test_nicknames"]
-    base_nickname = adfobj.case_nicknames["base_nickname"]
-
     res = adfobj.variable_defaults # will be dict of variable-specific plot preferences
-    # or an empty dictionary if use_defaults was not specified in the config YAML file.
+    # or an empty dictionary if use_defaults was not specified in YAML.
 
     #Set plot file type:
     # -- this should be set in basic_info_dict, but is not required
@@ -84,13 +75,6 @@ def meridional_mean(adfobj):
     print(f"\t NOTE: redo_plot is set to {redo_plot}")
     #-----------------------------------------
 
-    #Set data path variables:
-    #-----------------------
-    mclimo_rg_loc = Path(model_rgrid_loc)
-    if not adfobj.compare_obs:
-        dclimo_loc  = Path(data_loc)
-    #-----------------------
-
     #Set seasonal ranges:
     seasons = {"ANN": np.arange(1,13,1),
                "DJF": [12, 1, 2],
@@ -98,152 +82,211 @@ def meridional_mean(adfobj):
                "MAM": [3, 4, 5],
                "SON": [9, 10, 11]}
 
+    # Check if plots already exist and redo_plot boolean
+    # If redo_plot is false and file exists, keep track and skip the calculation
+    # entirely, which is what makes re-running the ADF fast.  A set, because
+    # every season of every variable is looked up in it once per case.
+    merid_skip = set()
+
+    # Loop over model cases:
+    for case_idx, case_name in enumerate(adfobj.data.case_names):
+        # Set output plot location:
+        plot_loc = Path(plot_locations[case_idx])
+
+        # Check if plot output directory exists, and if not, then create it:
+        if not plot_loc.is_dir():
+            print(f"    {plot_loc} not found, making new directory")
+            plot_loc.mkdir(parents=True)
+        # End if
+
+        # Loop over the variables for each season
+        for var in var_list:
+            for s in seasons:
+                plot_name = plot_loc / f"{var}_{s}_Meridional_Mean.{plot_type}"
+
+                # Check redo_plot. If set to True: remove old plot, if it already exists:
+                if (not redo_plot) and plot_name.is_file():
+                    merid_skip.add(plot_name)
+                    # Add already-existing plot to website (if enabled):
+                    adfobj.add_website_data(
+                        plot_name, var, case_name, season=s, plot_type="Meridional"
+                    )
+                    continue
+                elif (redo_plot) and plot_name.is_file():
+                    plot_name.unlink()
+                # End if
+            # End for (seasons)
+        # End for (variables)
+    # End for (cases)
+    #
+    # End redo plots check
+    #
+
+    #
+    # Setup Plotting
+    #
     #Loop over variables:
     for var in var_list:
         #Notify user of variable being plotted:
         print(f"\t - meridional mean plots for {var}")
 
-        if adfobj.compare_obs:
-            #Check if obs exist for the variable:
-            if var in var_obs_dict:
-                #Note: In the future these may all be lists, but for
-                #now just convert the target_list.
-                #Extract target file:
-                dclimo_loc = var_obs_dict[var]["obs_file"]
-                #Extract target list (eventually will be a list, for now need to convert):
-                data_list = [var_obs_dict[var]["obs_name"]]
-                #Extract target variable name:
-                data_var = var_obs_dict[var]["obs_var"]
-            else:
-                dmsg = f"No obs found for variable `{var}`, meridional mean plotting skipped."
-                adfobj.debug_log(dmsg)
-                continue
-            #End if
-        else:
-            #Set "data_var" for consistent use below:
-            data_var = var
+        # Nothing to draw for this variable at all: settle it from the file
+        # names, before opening any data
+        if _all_plots_present(
+            plot_locations, adfobj.data.case_names, var, seasons, plot_type, merid_skip
+        ):
+            continue
+        # End if
+
+        if var not in adfobj.data.ref_var_nam:
+            dmsg = f"\t    WARNING: No reference data found for variable `{var}`, meridional mean plotting skipped."
+            adfobj.debug_log(dmsg)
+            print(dmsg)
+            continue
         #End if
 
         # Check res for any variable specific options that need to be used BEFORE going to the plot:
         if var in res:
             vres = res[var]
-            #If found then notify user, assuming debug log is enabled:
-            adfobj.debug_log(f"meridional_mean: Found variable defaults for {var}")
-
+            # If found then notify user, assuming debug log is enabled:
+            adfobj.debug_log(
+                f"\t    INFO: meridional_mean: Found variable defaults for {var}"
+            )
         else:
             vres = {}
         #End if
 
-        #loop over different data sets to plot model against:
-        for data_src in data_list:
-            # load data (observational) comparison files
-            # (we should explore intake as an alternative to having this kind of repeated code):
-            if adfobj.compare_obs:
-                #For now, only grab one file (but convert to list for use below)
-                oclim_fils = [dclimo_loc]
-            else:
-                oclim_fils = sorted(dclimo_loc.glob(f"{data_src}_{var}_baseline.nc"))
+        # load reference data (observational or baseline)
+        if not adfobj.compare_obs:
+            base_name = adfobj.data.ref_case_label
+        else:
+            base_name = adfobj.data.ref_labels[var]
+        # End if
+
+        # Gather reference variable data
+        odata = adfobj.data.load_reference_regrid_da(base_name, var)
+
+        # Check if regridded file exists, if not skip meridional plot for this var
+        if odata is None:
+            dmsg = f"\t    WARNING: No regridded baseline file for {base_name} for variable `{var}`, meridional mean plotting skipped."
+            adfobj.debug_log(dmsg)
+            continue
+        # End if
+
+        # Check meridional mean dimensions
+        has_dims_ref = utils.validate_dims(odata, ["lat", "lev"])
+
+        # check if there is a lat dimension:
+        # if not, skip test cases and move to next variable
+        if not has_dims_ref["has_lat"]:
+            print(
+                f"\t    WARNING: Variable {var} is missing a lat dimension for '{base_name}', cannot continue to plot."
+            )
+            continue
+        # End if
+
+        # Loop over model cases:
+        for case_idx, case_name in enumerate(adfobj.data.case_names):
+
+            # Set case nickname:
+            case_nickname = adfobj.data.test_nicknames[case_idx]
+
+            # Set output plot location:
+            plot_loc = Path(plot_locations[case_idx])
+
+            # Nothing to draw for this case: settle it before opening the file
+            if all(
+                (plot_loc / f"{var}_{s}_Meridional_Mean.{plot_type}") in merid_skip
+                for s in seasons
+            ):
+                continue
             #End if
-            oclim_ds = utils.load_dataset(oclim_fils)
 
-            #Loop over model cases:
-            for case_idx, case_name in enumerate(case_names):
+            # load re-gridded model files:
+            mdata = adfobj.data.load_regrid_da(case_name, var)
+
+            if mdata is None:
+                dmsg = f"\t    WARNING: No regridded test file for {case_name} for variable `{var}`, meridional mean plotting skipped."
+                adfobj.debug_log(dmsg)
+                continue
+            # End if
 
-                #Set case nickname:
-                case_nickname = test_nicknames[case_idx]
+            # determine whether it's 2D or 3D
+            # 3D triggers search for surface pressure
+            has_dims = utils.validate_dims(mdata, ["lat", "lev"])
 
-                #Set output plot location:
-                plot_loc = Path(plot_locations[case_idx])
+            # check if there is a lat dimension:
+            if not has_dims["has_lat"]:
+                print(
+                    f"\t    WARNING: Variable {var} is missing a lat dimension for '{case_name}', cannot continue to plot."
+                )
+                continue
+            # End if
 
-                #Check if plot output directory exists, and if not, then create it:
-                if not plot_loc.is_dir():
-                    print(f"    {plot_loc} not found, making new directory")
-                    plot_loc.mkdir(parents=True)
+            has_lev = has_dims["has_lev"]
 
-                # load re-gridded model files:
-                mclim_fils = sorted(mclimo_rg_loc.glob(f"{data_src}_{case_name}_{var}_*.nc"))
-                mclim_ds = utils.load_dataset(mclim_fils)
+            # Notify user of level dimension:
+            if has_lev:
+                print(f"\t    INFO: {var} has lev dimension.")
+            # End if
 
-                # stop if data is invalid:
-                if (oclim_ds is None) or (mclim_ds is None):
-                    warnings.warn(f"invalid data, skipping meridional mean plot of {var}")
+            # Check to make sure each case has vertical levels if one of the cases does
+            if has_lev != has_dims_ref["has_lev"]:
+                print(
+                    f"\t    WARNING: expecting lev boolean for both case: {has_lev} and ref: {has_dims_ref['has_lev']}"
+                )
+                continue
+            # End if
+
+            #
+            # Seasonal Averages
+            # Note: xarray can do seasonal averaging, but depends on having time accessor,
+            # which these prototype climo files don't.
+            #
+
+            # Create new dictionaries:
+            mseasons = {}
+            oseasons = {}
+
+            # Loop over season dictionary:
+            for s in seasons:
+                # Set the file name
+                plot_name = plot_loc / f"{var}_{s}_Meridional_Mean.{plot_type}"
+
+                # Found above and redo_plot is false, so there is nothing to
+                # draw.  The seasonal averages below are the expensive part of
+                # a re-run, so work that out before doing them.
+                if plot_name in merid_skip:
                     continue
+                # End if
 
-                #Extract variable of interest
-                odata = oclim_ds[data_var].squeeze()  # squeeze in case of degenerate dimensions
-                mdata = mclim_ds[var].squeeze()
-
-                # APPLY UNITS TRANSFORMATION IF SPECIFIED:
-                # NOTE: looks like our climo files don't have all their metadata
-                mdata = mdata * vres.get("scale_factor",1) + vres.get("add_offset", 0)
-                # update units
-                mdata.attrs['units'] = vres.get("new_unit", mdata.attrs.get('units', 'none'))
-
-                # Do the same for the baseline case if need be:
-                if not adfobj.compare_obs:
-                    odata = odata * vres.get("scale_factor",1) + vres.get("add_offset", 0)
-                    # update units
-                    odata.attrs['units'] = vres.get("new_unit", odata.attrs.get('units', 'none'))
-                # Or for observations
-                else:
-                    odata = odata * vres.get("obs_scale_factor",1) + vres.get("obs_add_offset", 0)
-                    # Note: we are going to assume that the specification ensures the conversion makes the units the same.
-                    #       Doesn't make sense to add a different unit.
-
-                # determine whether it's 2D or 3D
-                # 3D triggers search for surface pressure
-                validate_lat_lev = utils.validate_dims(mdata, ['lat', 'lev']) # keys=> 'has_lat', 'has_lev', with T/F values
-
-                #Notify user of level dimension:
-                if validate_lat_lev['has_lev']:
-                    print(f"\t    INFO: {var} has lev dimension.")
-                    has_lev = True
-                else:
-                    has_lev = False
-
-                #
                 # Seasonal Averages
-                # Note: xarray can do seasonal averaging, but depends on having time accessor,
-                # which these prototype climo files don't.
-                #
-
-                #Create new dictionaries:
-                mseasons = {}
-                oseasons = {}
-
-                #Loop over season dictionary:
-                for s in seasons:
-                    plot_name = plot_loc / f"{var}_{s}_Meridional_Mean.{plot_type}"
-
-                    # Check redo_plot. If set to True: remove old plot, if it already exists:
-                    if (not redo_plot) and plot_name.is_file():
-                        #Add already-existing plot to website (if enabled):
-                        adfobj.debug_log(f"'{plot_name}' exists and clobber is false.")
-                        adfobj.add_website_data(plot_name, var, case_name, season=s,
-                                                plot_type="Meridional")
-                        #Continue to next iteration:
-                        continue
-                    elif (redo_plot) and plot_name.is_file():
-                        plot_name.unlink()
-
-                    mseasons[s] = utils.seasonal_mean(mdata, season=s, is_climo=True)
-                    oseasons[s] = utils.seasonal_mean(odata, season=s, is_climo=True)
-
-
-                    #Create new plot:
-                    pf.plot_meridional_mean_and_save(plot_name, case_nickname, base_nickname,
-                                                [syear_cases[case_idx],eyear_cases[case_idx]],
-                                                [syear_baseline,eyear_baseline],
-                                                mseasons[s], oseasons[s], has_lev, latbounds=slice(-5,5), obs=obs, **vres)
-
-                    #Add plot to website (if enabled):
-                    adfobj.add_website_data(plot_name, var, case_name, season=s,
-                                            plot_type="Meridional")
-
-                #End for (seasons loop)
-            #End for (case names loop)
-        #End for (obs/baseline loop)
-    #End for (variables loop)
+                mseasons[s] = utils.seasonal_mean(mdata, season=s, is_climo=True)
+                oseasons[s] = utils.seasonal_mean(odata, season=s, is_climo=True)
+
+                # Create new plot:
+                pf.plot_meridional_mean_and_save(
+                    plot_name,
+                    case_nickname,
+                    adfobj.data.ref_nickname,
+                    [syear_cases[case_idx], eyear_cases[case_idx]],
+                    [syear_baseline, eyear_baseline],
+                    mseasons[s],
+                    oseasons[s],
+                    has_lev,
+                    latbounds=slice(-5, 5),
+                    obs=adfobj.compare_obs,
+                    **vres,
+                )
+
+                # Add plot to website (if enabled):
+                adfobj.add_website_data(
+                    plot_name, var, case_name, season=s, plot_type="Meridional"
+                )
+
+            # End for (seasons loop)
+        # End for (case names loop)
+    # End for (variables loop)
 
     #Notify user that script has ended:
     print("  ...Meridional mean plots have been generated successfully.")
@@ -254,5 +297,41 @@ def meridional_mean(adfobj):
 #########
 
 
+def _all_plots_present(plot_locations, case_names, var, seasons, plot_type, merid_skip):
+    """
+    Report whether every plot this variable would produce is already there.
+
+    Parameters
+    ----------
+    plot_locations : list
+        output plot directory for each case
+    case_names : list
+        names of the test cases
+    var : str
+        ADF name of the variable being plotted
+    seasons : dict
+        the seasons being plotted, keyed by name
+    plot_type : str
+        file extension the plots are written with, e.g. "png"
+    merid_skip : set
+        plots found by the pre-flight pass above
+
+    Returns
+    -------
+    bool
+        ``True`` when nothing is left to draw for this variable, so the
+        reference and test data never have to be opened.
+    """
+    for case_idx in range(len(case_names)):
+        plot_loc = Path(plot_locations[case_idx])
+        for s in seasons:
+            if (plot_loc / f"{var}_{s}_Meridional_Mean.{plot_type}") not in merid_skip:
+                return False
+            # End if
+        # End for
+    # End for
+    return True
+
+
 ##############
-#END OF SCRIPT
\ No newline at end of file
+# END OF FILE
diff --git a/scripts/plotting/polar_map.py b/scripts/plotting/polar_map.py
index 7893fee8e..5d634b49a 100644
--- a/scripts/plotting/polar_map.py
+++ b/scripts/plotting/polar_map.py
@@ -4,6 +4,7 @@
 
 # ADF library
 import plotting_functions as pf
+import plotting_utils as plot_utils
 import adf_utils as utils
 
 def get_hemisphere(hemi_type):
@@ -138,10 +139,32 @@ def polar_map(adfobj):
         vres = res.get(var, {})
         web_category = vres.get("category", None)
 
+        # A complete set of plots can be recognised from the file names alone,
+        # because a 2-D variable and a 3-D one are named differently.  Doing
+        # that first means a re-run never opens the regridded files just to
+        # find out which of the two it is looking at.
+        settled = _existing_plot_set(
+            plot_locations, case_names, var, seasons, pres_levs, plot_type
+        )
+        if settled and not redo_plot:
+            for path, web_name, case_name, season, hemi_type in settled:
+                adfobj.add_website_data(
+                    path,
+                    web_name,
+                    case_name,
+                    category=web_category,
+                    season=season,
+                    plot_type=hemi_type,
+                )
+            # End for
+            print(f"\t    Skipping {var} - all plots already exist")
+            continue
+        # End if
+
         # Get all plot info and check existence
         plot_info = []
         all_plots_exist = True
-        
+
         for case_idx, case_name in enumerate(case_names):
             plot_loc = Path(plot_locations[case_idx])
 
@@ -266,5 +289,73 @@ def polar_map(adfobj):
 ##############
 #END OF `polar_map` function
 
+
+def _existing_plot_set(plot_locations, case_names, var, seasons, pres_levs, plot_type):
+    """
+    Return the complete set of polar plots for `var`, if one is on disk.
+
+    A 2-D variable is plotted as ``{var}_{season}_{hemisphere}_Mean`` and a
+    3-D one as ``{var}_{pressure}hpa_{season}_{hemisphere}_Mean``, so which
+    of the two a variable is can be read off the file names.  That is what
+    lets a re-run decide there is nothing to do without opening the regridded
+    file to look for a ``lev`` dimension.
+
+    Parameters
+    ----------
+    plot_locations : list
+        output plot directory for each case
+    case_names : list
+        names of the test cases
+    var : str
+        ADF name of the variable being plotted
+    seasons : list or dict
+        the seasons being plotted
+    pres_levs : list
+        the configured pressure levels, empty when none are set
+    plot_type : str
+        file extension the plots are written with, e.g. "png"
+
+    Returns
+    -------
+    list
+        ``(path, website name, case, season, hemisphere)`` for every plot,
+        when a complete set exists; an empty list when it does not.
+    """
+    flat = []
+    levelled = []
+    for case_idx, case_name in enumerate(case_names):
+        plot_loc = Path(plot_locations[case_idx])
+        for season in seasons:
+            for hemi_type in ["NHPolar", "SHPolar"]:
+                flat.append(
+                    (
+                        plot_loc / f"{var}_{season}_{hemi_type}_Mean.{plot_type}",
+                        var,
+                        case_name,
+                        season,
+                        hemi_type,
+                    )
+                )
+                for pres in pres_levs:
+                    levelled.append(
+                        (
+                            plot_loc
+                            / (
+                                f"{var}_{pres}hpa_{season}_{hemi_type}"
+                                f"_Mean.{plot_type}"
+                            ),
+                            f"{var}_{pres}hpa",
+                            case_name,
+                            season,
+                            hemi_type,
+                        )
+                    )
+                # End for
+            # End for
+        # End for
+    # End for
+    return plot_utils.first_complete_plot_set([levelled, flat])
+
+
 ##############
-# END OF FILE
\ No newline at end of file
+# END OF FILE
diff --git a/scripts/plotting/tem.py b/scripts/plotting/tem.py
index 13f9c49d4..a72c153ec 100644
--- a/scripts/plotting/tem.py
+++ b/scripts/plotting/tem.py
@@ -408,19 +408,17 @@ def tem(adf):
                 mdata = ds[var].squeeze()
                 odata = ds_base[var].squeeze()
 
-                #Apply the unit conversion from the variable defaults. TEM
-                #files carry little metadata, so the new unit is taken from the
-                #defaults too. Observations have their own scaling, which is
-                #assumed to bring them to the same units, so they keep the unit
-                #string they arrived with.
-                mdata = mdata * vres.get("scale_factor", 1) + vres.get("add_offset", 0)
-                mdata.attrs['units'] = vres.get("new_unit", mdata.attrs.get('units', 'none'))
+                # Apply the unit conversion from the variable defaults, through
+                # the ADF's data layer so that a TEM file already holding
+                # converted values is not scaled a second time. Observations
+                # have their own scaling, which is assumed to bring them to the
+                # same units, so they keep the unit string they arrived with.
+                mdata = adf.data.apply_conversion(mdata, var)
                 if obs:
                     odata = (odata * vres.get("obs_scale_factor", 1)
                              + vres.get("obs_add_offset", 0))
                 else:
-                    odata = odata * vres.get("scale_factor", 1) + vres.get("add_offset", 0)
-                    odata.attrs['units'] = vres.get("new_unit", odata.attrs.get('units', 'none'))
+                    odata = adf.data.apply_conversion(odata, var)
                 #End if
 
                 #Month-length weighted seasonal (or annual) mean. The weighted
@@ -571,4 +569,4 @@ def _mesh(da):
     print("  ...TEM plots have been generated successfully.")
 
 # Helper functions
-##################
\ No newline at end of file
+##################
diff --git a/scripts/plotting/zonal_mean.py b/scripts/plotting/zonal_mean.py
index 2f8630fb3..97a19928f 100644
--- a/scripts/plotting/zonal_mean.py
+++ b/scripts/plotting/zonal_mean.py
@@ -100,7 +100,7 @@ def zonal_mean(adfobj):
         for var in var_list:
             for s in seasons:
                 #Check zonal log-p:
-                plot_name_log = plot_loc / f"{var}_{s}_Zonal_logp_Mean.{plot_type}"
+                plot_name_log = plot_loc / f"{var}_logp_{s}_Zonal_Mean.{plot_type}"
 
                 # Check redo_plot. If set to True: remove old plot, if it already exists:
                 if (not redo_plot) and plot_name_log.is_file():
@@ -252,10 +252,6 @@ def zonal_mean(adfobj):
                 # because we can let any pressure-level interpolation happen there
                 # This could be re-visited for efficiency or improved code structure.
 
-                #Seasonal Averages
-                mseasons[s] = utils.seasonal_mean(mdata, season=s, is_climo=True)
-                oseasons[s] = utils.seasonal_mean(odata, season=s, is_climo=True)
-
                 #Set the file name
                 plot_name = plot_loc / f"{var}_{s}_Zonal_Mean.{plot_type}"
                 plot_name_log = None
@@ -265,6 +261,19 @@ def zonal_mean(adfobj):
                     plot_name_log = plot_loc / f"{var}_logp_{s}_Zonal_Mean.{plot_type}"
                 #End if
 
+                # Nothing to draw for this season: both plots were found above
+                # and redo_plot is false.  The seasonal averages below are the
+                # expensive part of a re-run, so work that out before doing them.
+                if (plot_name in zonal_skip) and (
+                    (plot_name_log is None) or (plot_name_log in logp_zonal_skip)
+                ):
+                    continue
+                # End if
+
+                # Seasonal Averages
+                mseasons[s] = utils.seasonal_mean(mdata, season=s, is_climo=True)
+                oseasons[s] = utils.seasonal_mean(odata, season=s, is_climo=True)
+
                 #Create plots
                 if plot_name not in zonal_skip:
 
@@ -299,4 +308,4 @@ def zonal_mean(adfobj):
 
 
 ##############
-#END OF SCRIPT
\ No newline at end of file
+#END OF SCRIPT