diff --git a/ci/recipe/meta.yaml b/ci/recipe/meta.yaml
index be418c4bf..34d02fee8 100644
--- a/ci/recipe/meta.yaml
+++ b/ci/recipe/meta.yaml
@@ -1,5 +1,5 @@
 {% set name = "MPAS-Analysis" %}
-{% set version = "1.9.0" %}
+{% set version = "1.9.1" %}
 
 package:
   name: {{ name|lower }}
@@ -35,7 +35,7 @@ requirements:
     - lxml
     - mache >=1.11.0
     - matplotlib-base >=3.6.0,!=3.7.2
-    - mpas_tools >=0.16.0
+    - mpas_tools >=0.30.0
     - nco >=4.8.1
     - netcdf4
     - numpy
diff --git a/dev-spec.txt b/dev-spec.txt
index 27fbe93b3..af35884ae 100644
--- a/dev-spec.txt
+++ b/dev-spec.txt
@@ -18,7 +18,7 @@ mache >=1.11.0
 # 3.7.2 contains a bug with tight layouts and insets
 # https://github.com/matplotlib/matplotlib/pull/26291
 matplotlib-base>=3.6.0,!=3.7.2
-mpas_tools>=0.16.0
+mpas_tools>=0.30.0
 nco>=4.8.1
 netcdf4
 numpy
diff --git a/mpas_analysis/__init__.py b/mpas_analysis/__init__.py
index c8a904509..0e791a096 100644
--- a/mpas_analysis/__init__.py
+++ b/mpas_analysis/__init__.py
@@ -3,5 +3,5 @@
 import matplotlib as mpl
 mpl.use('Agg')
 
-__version_info__ = (1, 9, 0)
+__version_info__ = (1, 9, 1)
 __version__ = '.'.join(str(vi) for vi in __version_info__)
diff --git a/mpas_analysis/shared/html/pages.py b/mpas_analysis/shared/html/pages.py
index d8f8acbaf..a2292597e 100644
--- a/mpas_analysis/shared/html/pages.py
+++ b/mpas_analysis/shared/html/pages.py
@@ -92,6 +92,9 @@ def generate_html(config, analyses, controlConfig, customConfigFiles):
             html_dir = config.get('output', 'htmlSubdirectory')
             if html_dir.startswith(base_path):
                 url = base_url + html_dir[len(base_path):]
+                if not url.endswith('/'):
+                    # lack of trailing '/' is causing issues on NERSC's portal
+                    url = f'{url}/'
                 print(f'Web page: {url}')
     if url is None:
         print("Done.")
diff --git a/mpas_analysis/shared/plot/climatology_map.py b/mpas_analysis/shared/plot/climatology_map.py
index 1a81ca4b3..1f9a16d40 100644
--- a/mpas_analysis/shared/plot/climatology_map.py
+++ b/mpas_analysis/shared/plot/climatology_map.py
@@ -702,7 +702,15 @@ def plot_panel(ax, title, array, colormap, norm, levels, ticks, contours,
 def _add_stats(modelArray, refArray, diffArray, Lats, axes):
     """ compute the means, std devs. and Pearson correlation """
     weights = np.cos(np.deg2rad(Lats))
-    modelMean = np.average(modelArray, weights=weights)
+
+    model_weights = weights
+    model_mask = None
+    if isinstance(modelArray, np.ma.MaskedArray):
+        # make sure we're using the MPAS land mask for all 3 sets of stats
+        model_mask = modelArray.mask
+        model_weights = np.ma.array(weights, mask=model_mask)
+
+    modelMean = np.average(modelArray, weights=model_weights)
 
     _add_stats_text(
         names=['Min', 'Mean', 'Max'],
@@ -710,26 +718,28 @@ def _add_stats(modelArray, refArray, diffArray, Lats, axes):
         ax=axes[0], loc='upper')
 
     if refArray is not None:
+        ref_weights = weights
+        ref_mask = None
         if isinstance(modelArray, np.ma.MaskedArray):
             # make sure we're using the MPAS land mask for all 3 sets of stats
-            mask = modelArray.mask
             if isinstance(refArray, np.ma.MaskedArray):
                 # mask invalid where either model or ref array is invalid
-                mask = np.logical_or(mask, refArray.mask)
-            refArray = np.ma.array(refArray, mask=mask)
+                ref_mask = np.logical_or(model_mask, refArray.mask)
+                ref_weights = np.ma.array(weights, mask=ref_mask)
+            refArray = np.ma.array(refArray, mask=ref_mask)
         modelAnom = modelArray - modelMean
-        modelVar = np.average(modelAnom ** 2, weights=weights)
-        refMean = np.average(refArray, weights=weights)
+        modelVar = np.average(modelAnom ** 2, weights=ref_weights)
+        refMean = np.average(refArray, weights=ref_weights)
         refAnom = refArray - refMean
-        refVar = np.average(refAnom**2, weights=weights)
+        refVar = np.average(refAnom**2, weights=ref_weights)
 
         _add_stats_text(
             names=['Min', 'Mean', 'Max'],
             values=[np.amin(refArray), refMean, np.amax(refArray)],
             ax=axes[1], loc='upper')
 
-        diffMean = np.average(diffArray, weights=weights)
-        diffVar = np.average((diffArray - diffMean)**2, weights=weights)
+        diffMean = np.average(diffArray, weights=ref_weights)
+        diffVar = np.average((diffArray - diffMean)**2, weights=ref_weights)
         diffRMSE = np.sqrt(diffVar)
 
         _add_stats_text(
@@ -737,7 +747,7 @@ def _add_stats(modelArray, refArray, diffArray, Lats, axes):
             values=[np.amin(diffArray), diffMean, np.amax(diffArray)],
             ax=axes[2], loc='upper')
 
-        covar = np.average(modelAnom*refAnom, weights=weights)
+        covar = np.average(modelAnom*refAnom, weights=ref_weights)
 
         corr = covar/np.sqrt(modelVar*refVar)