diff --git a/compass/ocean/tests/tides/__init__.py b/compass/ocean/tests/tides/__init__.py
index ec029d4db5..0288850c34 100644
--- a/compass/ocean/tests/tides/__init__.py
+++ b/compass/ocean/tests/tides/__init__.py
@@ -1,8 +1,7 @@
-from compass.testgroup import TestGroup
-
-from compass.ocean.tests.tides.mesh import Mesh
-from compass.ocean.tests.tides.init import Init
 from compass.ocean.tests.tides.forward import Forward
+from compass.ocean.tests.tides.init import Init
+from compass.ocean.tests.tides.mesh import Mesh
+from compass.testgroup import TestGroup
 
 
 class Tides(TestGroup):
@@ -17,7 +16,7 @@ def __init__(self, mpas_core):
         super().__init__(mpas_core=mpas_core,
                          name='tides')
 
-        for mesh_name in ['Icos7']:
+        for mesh_name in ['Icos7', 'VR45to5']:
 
             mesh = Mesh(test_group=self, mesh_name=mesh_name)
             self.add_test_case(mesh)
diff --git a/compass/ocean/tests/tides/analysis/__init__.py b/compass/ocean/tests/tides/analysis/__init__.py
index 789e85dbc7..cda4f59676 100644
--- a/compass/ocean/tests/tides/analysis/__init__.py
+++ b/compass/ocean/tests/tides/analysis/__init__.py
@@ -1,14 +1,15 @@
-from compass.step import Step
+import os
 
-import netCDF4
-import matplotlib.pyplot as plt
-import matplotlib.cm as cm
-import numpy as np
 import cartopy.crs as ccrs
 import cartopy.feature as cfeature
-import os
+import matplotlib.cm as cm
+import matplotlib.pyplot as plt
+import netCDF4
+import numpy as np
 from mpas_tools.logging import check_call
 
+from compass.step import Step
+
 
 class Analysis(Step):
     """
@@ -103,7 +104,7 @@ def write_coordinate_file(self, idx):
         # Write coordinate file for OTPS2
         f = open('lat_lon', 'w')
         for i in range(nCells):
-            f.write(str(lat_grid[i])+'  '+str(lon_grid[i])+'\n')
+            f.write(str(lat_grid[i]) + '  ' + str(lon_grid[i]) + '\n')
         f.close()
 
     def setup_otps2(self):
@@ -121,7 +122,7 @@ def setup_otps2(self):
                       'comment': '! 2. latitude/longitude/<time> file'},
                      {'inp': 'z',
                       'comment': '! 3. z/U/V/u/v'},
-                     {'inp':  con,
+                     {'inp': con,
                       'comment': '! 4. tidal constituents to include'},
                      {'inp': 'AP',
                       'comment': '! 5. AP/RI'},
@@ -137,10 +138,10 @@ def setup_otps2(self):
                 os.mkdir('inputs')
 
             # Write the setup_con file
-            f = open('inputs/'+con+'_setup', 'w')
+            f = open(f'inputs/{con}_setup', 'w')
             for line in lines:
                 spaces = 28 - len(line['inp'])
-                f.write(line['inp'] + spaces*' ' + line['comment'] + '\n')
+                f.write(line['inp'] + spaces * ' ' + line['comment'] + '\n')
             f.close()
 
             # Write the Model_atlas_con file
@@ -180,18 +181,18 @@ def read_otps2_output(self, idx):
                 line_sp = line.split()
                 if line_sp[2] != '*************':
                     val = float(line_sp[2])
-                    self.mesh_AP[con]['amp'][start+i] = val
+                    self.mesh_AP[con]['amp'][start + i] = val
                 else:
-                    self.mesh_AP[con]['amp'][start+i] = -9999
+                    self.mesh_AP[con]['amp'][start + i] = -9999
 
                 if line_sp[3] != 'Site':
                     val = float(line_sp[3])
                     if val < 0:
                         val = val + 360.0
-                    self.mesh_AP[con]['phase'][start+i] = val
+                    self.mesh_AP[con]['phase'][start + i] = val
 
                 else:
-                    self.mesh_AP[con]['phase'][start+i] = -9999
+                    self.mesh_AP[con]['phase'][start + i] = -9999
 
     def append_tpxo_data(self):
         """
@@ -268,8 +269,8 @@ def plot(self):
 
         lon_grid = np.mod(data_nc.variables['lonCell'][:] + np.pi,
                           2.0 * np.pi) - np.pi
-        lon_grid = lon_grid*180.0/np.pi
-        lat_grid = data_nc.variables['latCell'][:]*180.0/np.pi
+        lon_grid = lon_grid * 180.0 / np.pi
+        lat_grid = data_nc.variables['latCell'][:] * 180.0 / np.pi
 
         nCells = lon_grid.size
         data1 = np.zeros((nCells))
@@ -287,14 +288,14 @@ def plot(self):
         # Use these to fix up the plots
         subplot_ticks = [[np.linspace(0, 0.65, 10), np.linspace(0, 0.65, 10),
                           np.linspace(0, 0.13, 10), np.linspace(0, 0.13, 10)],
-                         [np.linspace(0, 1.4,  10), np.linspace(0, 1.4,  10),
+                         [np.linspace(0, 1.40, 10), np.linspace(0, 1.40, 10),
                           np.linspace(0, 0.22, 10), np.linspace(0, 0.25, 10)],
                          [np.linspace(0, 0.22, 10), np.linspace(0, 0.22, 10),
                           np.linspace(0, 0.05, 10), np.linspace(0, 0.05, 10)],
-                         [np.linspace(0, 0.5,  10), np.linspace(0, 0.5,  10),
+                         [np.linspace(0, 0.50, 10), np.linspace(0, 0.50, 10),
                           np.linspace(0, 0.08, 10), np.linspace(0, 0.08, 10)],
-                         [np.linspace(0, 0.7,  10), np.linspace(0, 0.7,  10),
-                          np.linspace(0, 0.5,  10), np.linspace(0, 0.5,  10)]]
+                         [np.linspace(0, 0.70, 10), np.linspace(0, 0.70, 10),
+                          np.linspace(0, 0.50, 10), np.linspace(0, 0.50, 10)]]
 
         for i, con in enumerate(constituent_list):
 
@@ -311,39 +312,45 @@ def plot(self):
             data2_phase[:] = data_nc.variables[
                 f'{con}Phase{self.tpxo_version}'][:]
 
-            data1_phase = data1_phase*np.pi/180.0
-            data2_phase = data2_phase*np.pi/180.0
+            data1_phase = data1_phase * np.pi / 180.0
+            data2_phase = data2_phase * np.pi / 180.0
 
             # Calculate RMSE values
-            rmse_amp = 0.5*(data1 - data2)**2
-            rmse_com = 0.5*(data2**2 + data1**2) \
-                - data1*data2*np.cos(data2_phase - data1_phase)
+            rmse_amp = 0.5 * (data1 - data2)**2
+            rmse_com = 0.5 * (data2**2 + data1**2) \
+                - data1 * data2 * np.cos(data2_phase - data1_phase)
 
             # Calculate mean (global) values
-            idx = np.where((depth > 20)
-                           & (rmse_com < 1000) & (rmse_amp < 1000))
+            idx = np.where((depth > 20) &
+                           (rmse_com < 1000) & (rmse_amp < 1000))
             area_tot = np.sum(area[idx])
-            glo_rmse_amp = np.sqrt(np.sum(rmse_amp[idx]*area[idx])/area_tot)
-            glo_rmse_com = np.sqrt(np.sum(rmse_com[idx]*area[idx])/area_tot)
+            glo_rmse_amp = np.sqrt(np.sum(rmse_amp[idx] * area[idx]) /
+                                   area_tot)
+            glo_rmse_com = np.sqrt(np.sum(rmse_com[idx] * area[idx]) /
+                                   area_tot)
             print('Global RMSE (Amp) = ', glo_rmse_amp)
             print('Global RMSE (Com) = ', glo_rmse_com)
 
             # Calculate shallow RMSE (<=1000m)
-            idx = np.where((depth > 20) & (depth < 1000)
-                           & (np.abs(lat_grid) < 66)
-                           & (rmse_com < 1000) & (rmse_amp < 1000))
+            idx = np.where((depth > 20) & (depth < 1000) &
+                           (np.abs(lat_grid) < 66) &
+                           (rmse_com < 1000) & (rmse_amp < 1000))
             area_tot = np.sum(area[idx])
-            shal_rmse_amp = np.sqrt(np.sum(rmse_amp[idx]*area[idx])/area_tot)
-            shal_rmse_com = np.sqrt(np.sum(rmse_com[idx]*area[idx])/area_tot)
+            shal_rmse_amp = np.sqrt(np.sum(rmse_amp[idx] * area[idx]) /
+                                    area_tot)
+            shal_rmse_com = np.sqrt(np.sum(rmse_com[idx] * area[idx]) /
+                                    area_tot)
             print('Shallow RMSE (Amp) = ', shal_rmse_amp)
             print('Shallow RMSE (Com) = ', shal_rmse_com)
 
             # Calculate deep RMSE (>1000m)
-            idx = np.where((depth >= 1000) & (np.abs(lat_grid) < 66)
-                           & (rmse_com < 1000) & (rmse_amp < 1000))
+            idx = np.where((depth >= 1000) & (np.abs(lat_grid) < 66) &
+                           (rmse_com < 1000) & (rmse_amp < 1000))
             area_tot = np.sum(area[idx])
-            deep_rmse_amp = np.sqrt(np.sum(rmse_amp[idx]*area[idx])/area_tot)
-            deep_rmse_com = np.sqrt(np.sum(rmse_com[idx]*area[idx])/area_tot)
+            deep_rmse_amp = np.sqrt(np.sum(rmse_amp[idx] * area[idx]) /
+                                    area_tot)
+            deep_rmse_com = np.sqrt(np.sum(rmse_com[idx] * area[idx]) /
+                                    area_tot)
             print('Deep RMSE (Amp) = ', deep_rmse_amp)
             print('Deep RMSE (Com) = ', deep_rmse_com)
 
@@ -353,12 +360,12 @@ def plot(self):
             # Plot data
             fig = plt.figure(figsize=(18, 12))
             subplot_title = [f'{con} Amplitude (simulation) [m]',
-                             f'{con} Amplitude (TPXO8) [m]',
+                             f'{con} Amplitude ({self.tpxo_version}) [m]',
                              f'{con} RMSE (Amplitude) [m]',
                              f'{con} RMSE (Complex) [m]']
 
             for subplot in range(0, 4):
-                ax = fig.add_subplot(2, 2, subplot+1,
+                ax = fig.add_subplot(2, 2, subplot + 1,
                                      projection=ccrs.PlateCarree())
                 ax.set_title(subplot_title[subplot], fontsize=20)
                 levels = subplot_ticks[i][subplot][:]
@@ -407,9 +414,9 @@ def plot(self):
                 cbar.ax.tick_params(labelsize=16)
 
             fig.tight_layout()
-            global_err = str(round(glo_rmse_com*100, 3))
-            deep_err = str(round(deep_rmse_com*100, 3))
-            shallow_err = str(round(shal_rmse_com*100, 3))
+            global_err = str(round(glo_rmse_com * 100, 3))
+            deep_err = str(round(deep_rmse_com * 100, 3))
+            shallow_err = str(round(shal_rmse_com * 100, 3))
             fig.suptitle(f'Complex RMSE: Global = {global_err} cm; '
                          f'Deep = {deep_err} cm; '
                          f'Shallow = {shallow_err} cm',
@@ -430,7 +437,7 @@ def run(self):
 
         # Setup chunking for TPXO extraction with large meshes
         indices = np.arange(self.nCells)
-        nchunks = np.ceil(self.nCells/200000)
+        nchunks = np.ceil(self.nCells / 200000)
         index_chunks = np.array_split(indices, nchunks)
 
         # Initialize data structure for TPXO values
diff --git a/compass/ocean/tests/tides/dem/__init__.py b/compass/ocean/tests/tides/dem/__init__.py
new file mode 100644
index 0000000000..8db0eceeba
--- /dev/null
+++ b/compass/ocean/tests/tides/dem/__init__.py
@@ -0,0 +1,59 @@
+import os
+
+import compass.ocean.tests.tides.dem.dem_pixel as dem_pixel
+from compass.step import Step
+
+
+class CreatePixelFile(Step):
+    """
+    A step for creating a pixel file for creating MPAS meshes
+
+    Attributes
+    ----------
+    """
+
+    def __init__(self, test_case):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        test_case : compass.ocean.tests.tides.init.Init
+            The test case this step belongs to
+        """
+        super().__init__(test_case=test_case, name='pixel',
+                         ntasks=1, min_tasks=1, openmp_threads=1)
+
+        self.add_input_file(
+            filename='GEBCO_2023_sub_ice_topo.nc',
+            target='GEBCO_2023_sub_ice_topo.nc',
+            database='bathymetry_database')
+
+        self.add_input_file(
+            filename='RTopo-2.0.4_30sec_bedrock_topography.nc',
+            target='RTopo-2.0.4_30sec_bedrock_topography.nc',
+            database='bathymetry_database')
+
+        self.add_input_file(
+            filename='RTopo-2.0.4_30sec_surface_elevation.nc',
+            target='RTopo-2.0.4_30sec_surface_elevation.nc',
+            database='bathymetry_database')
+
+        self.add_input_file(
+            filename='RTopo-2.0.4_30sec_ice_base_topography.nc',
+            target='RTopo-2.0.4_30sec_ice_base_topography.nc',
+            database='bathymetry_database')
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+
+        self.init_path = './'
+
+        if not os.path.exists('RTopo_2_0_4_30sec_pixel.nc'):
+            dem_pixel.rtopo_30sec(self.init_path, self.init_path)
+        if not os.path.exists('GEBCO_v2023_30sec_pixel.nc'):
+            dem_pixel.gebco_30sec(self.init_path, self.init_path)
+        if not os.path.exists('RTopo_2_0_4_GEBCO_v2023_30sec_pixel.nc'):
+            dem_pixel.rtopo_gebco_30sec(self.init_path, self.init_path)
diff --git a/compass/ocean/tests/tides/dem/dem_names.py b/compass/ocean/tests/tides/dem/dem_names.py
new file mode 100644
index 0000000000..a30957d570
--- /dev/null
+++ b/compass/ocean/tests/tides/dem/dem_names.py
@@ -0,0 +1,21 @@
+
+# -- long_name strings for NetCDF variables
+
+class base:
+    pass
+
+
+names = base()
+names.bed_elevation = "elevation of bed"
+names.ocn_thickness = "thickness of ocn"
+names.ice_thickness = "thickness of ice"
+names.ocn_cover = "fractional cover of ocn, 0-1"
+names.ice_cover = "fractional cover of ice, 0-1"
+names.bed_slope = "RMS magnitude of bed slopes"
+names.bed_slope_deg = "arc-tangent of RMS bed slopes"
+names.bed_dz_dx = "derivative of bed elevation along lon.-axis"
+names.bed_dz_dy = "derivative of bed elevation along lat.-axis"
+names.bed_elevation_profile = "sub-grid percentiles of bed elevation"
+names.bed_slope_profile = "sub-grid percentiles of RMS bed slope"
+names.ocn_thickness_profile = "sub-grid percentiles of ocn thickness"
+names.ice_thickness_profile = "sub-grid percentiles of ice thickness"
diff --git a/compass/ocean/tests/tides/init/dem_pixel.py b/compass/ocean/tests/tides/dem/dem_pixel.py
similarity index 56%
rename from compass/ocean/tests/tides/init/dem_pixel.py
rename to compass/ocean/tests/tides/dem/dem_pixel.py
index 337f5eacda..31c226f55b 100644
--- a/compass/ocean/tests/tides/init/dem_pixel.py
+++ b/compass/ocean/tests/tides/dem/dem_pixel.py
@@ -6,8 +6,229 @@
 import numpy as np
 from scipy.ndimage import gaussian_filter
 
+from compass.ocean.tests.tides.dem.dem_names import names
+
 # Authors: Darren Engwirda
 
+RSPH = 6371220.
+
+
+def circ_dist(xa, ya, xb, yb, rs=1.):
+    """
+    Calculate geodesic-length of great circles [PA, PB] on a
+    sphere of radius RS.
+
+    """
+
+    dlon = .5 * (xa - xb)
+    dlat = .5 * (ya - yb)
+
+    dist = 2. * rs * np.arcsin(np.sqrt(
+        np.sin(dlat) ** 2 +
+        np.sin(dlon) ** 2 * np.cos(ya) * np.cos(yb)
+    ))
+
+    return dist
+
+
+def cell_dzdx(xlon, ylat, vals, rsph):
+    """
+    Eval. the lon.- and lat.-aligned elevation gradients and
+    their magnitude via a D8 stencil.
+
+    """
+
+    print("Building slopes...")
+
+    cols = xlon.size - 2
+    rows = ylat.size - 2
+
+    xmid = .5 * (xlon[:-1:] + xlon[1::])
+    ymid = .5 * (ylat[:-1:] + ylat[1::])
+
+    ridx, icol = np.ogrid[:rows + 1, :cols + 1]
+
+    # dz/dx is poorly conditioned at poles
+    beta = (np.tanh((ymid + 87.5) * 2.5) -
+            np.tanh((ymid - 87.5) * 2.5)) * 0.5
+    beta = np.reshape(beta, (beta.size, 1))
+
+    filt = np.asarray(vals, dtype=np.float32)
+
+    part = rows // 20
+    fbot = gaussian_filter(filt[+0:part * +1, :],
+                           sigma=(4., cols / 512.), mode=("reflect", "wrap"))
+
+    ftop = gaussian_filter(filt[19 * part:-1, :],
+                           sigma=(4., cols / 512.), mode=("reflect", "wrap"))
+
+    filt[+0:part * +1, :] = fbot
+    filt[19 * part:-1, :] = ftop
+
+    vals *= beta  # careful with mem. alloc.
+    beta = (+1. - beta)
+    vals += beta * filt
+    # vals = beta * vals + (1. - beta) * filt
+
+    del filt
+    del ftop
+    del fbot
+
+    xmid = xmid * np.pi / 180.
+    ymid = ymid * np.pi / 180.
+
+    dzds = np.zeros((
+        rows + 1, cols + 1), dtype=np.float32)
+
+    dzdx = np.zeros((
+        rows + 1, cols + 1), dtype=np.float32)
+    dzdy = np.zeros((
+        rows + 1, cols + 1), dtype=np.float32)
+
+    indx = np.asarray(np.round(
+        np.linspace(-1, rows, 32)), dtype=np.int64)
+
+    print("* process tiles:")
+
+    for tile in range(0, indx.size - 1):
+
+        head = indx[tile + 0] + 1
+        tail = indx[tile + 1] + 1
+
+        slab = tail - head + 0
+
+        irow = ridx[head:tail]
+
+        zdel = np.zeros((
+            slab + 0, cols + 1, 8), dtype=np.float32)
+
+        xdel = np.zeros((
+            slab + 0, cols + 1, 1), dtype=np.float32)
+        ydel = np.zeros((
+            slab + 0, cols + 1, 1), dtype=np.float32)
+
+        zero = 0
+        wcol = icol - 1
+        wcol[wcol < zero] = cols
+        ecol = icol + 1
+        ecol[ecol > cols] = zero
+        nrow = irow + 1
+        nrow[nrow > rows] = rows
+        srow = irow - 1
+        srow[srow < zero] = zero
+
+        # -- index D4 neighbours
+
+        xdel[:, :, 0] = \
+            vals[irow, ecol] - vals[irow, wcol]
+
+        dist = circ_dist(xmid[ecol], ymid[irow],
+                         xmid[wcol], ymid[irow])
+
+        dist += 1.E-12
+        xdel[:, :, 0] /= dist * rsph
+
+        ydel[:, :, 0] = \
+            vals[nrow, icol] - vals[srow, icol]
+
+        dist = circ_dist(xmid[icol], ymid[nrow],
+                         xmid[icol], ymid[srow])
+
+        dist += 1.E-12
+        ydel[:, :, 0] /= dist * rsph
+
+        # -- index D8 neighbours =
+        # -- NN, EE, SS, WW, NE, SE, SW, NW
+
+        zdel[:, :, 0] = np.abs(
+            vals[nrow, icol] - vals[irow, icol])
+        zdel[:, :, 1] = np.abs(
+            vals[irow, ecol] - vals[irow, icol])
+        zdel[:, :, 2] = np.abs(
+            vals[srow, icol] - vals[irow, icol])
+        zdel[:, :, 3] = np.abs(
+            vals[irow, wcol] - vals[irow, icol])
+
+        zdel[:, :, 4] = np.abs(
+            vals[nrow, ecol] - vals[irow, icol])
+        zdel[:, :, 5] = np.abs(
+            vals[srow, ecol] - vals[irow, icol])
+        zdel[:, :, 6] = np.abs(
+            vals[srow, wcol] - vals[irow, icol])
+        zdel[:, :, 7] = np.abs(
+            vals[nrow, wcol] - vals[irow, icol])
+
+        dist = circ_dist(xmid[icol], ymid[irow],
+                         xmid[icol], ymid[nrow])
+
+        dist += 1.E-12
+        zdel[:, :, 0] /= dist * rsph
+
+        dist = circ_dist(xmid[icol], ymid[irow],
+                         xmid[ecol], ymid[irow])
+
+        dist += 1.E-12
+        zdel[:, :, 1] /= dist * rsph
+
+        dist = circ_dist(xmid[icol], ymid[irow],
+                         xmid[icol], ymid[srow])
+
+        dist += 1.E-12
+        zdel[:, :, 2] /= dist * rsph
+
+        dist = circ_dist(xmid[icol], ymid[irow],
+                         xmid[wcol], ymid[irow])
+
+        dist += 1.E-12
+        zdel[:, :, 3] /= dist * rsph
+
+        dist = circ_dist(xmid[icol], ymid[irow],
+                         xmid[ecol], ymid[nrow])
+
+        dist += 1.E-12
+        zdel[:, :, 4] /= dist * rsph
+
+        dist = circ_dist(xmid[icol], ymid[irow],
+                         xmid[ecol], ymid[srow])
+
+        dist += 1.E-12
+        zdel[:, :, 5] /= dist * rsph
+
+        dist = circ_dist(xmid[icol], ymid[irow],
+                         xmid[wcol], ymid[srow])
+
+        dist += 1.E-12
+        zdel[:, :, 6] /= dist * rsph
+
+        dist = circ_dist(xmid[icol], ymid[irow],
+                         xmid[wcol], ymid[nrow])
+
+        dist += 1.E-12
+        zdel[:, :, 7] /= dist * rsph
+
+        dzds[irow, icol] = \
+            np.sqrt(np.mean(zdel**2, axis=2))
+
+        dzdx[irow, icol] = xdel[:, :, 0]
+        dzdy[irow, icol] = ydel[:, :, 0]
+
+        del zdel
+        del xdel
+        del ydel
+
+        print("* compute local D8 slope:",
+              tile, "of", indx.size - 1)
+
+    dzds[+1, :] = dzds[+2, :]
+    dzdx[+1, :] = dzdx[+2, :]
+    dzdy[+1, :] = dzdy[+2, :]
+
+    dzds[-1, :] = dzds[-2, :]
+    dzdx[-1, :] = dzdx[-2, :]
+    dzdy[-1, :] = dzdy[-2, :]
+
+    return dzds, dzdx, dzdy
+
 
 def blend_front(e1st, i1st, e2nd, halo, sdev):
     """
@@ -59,13 +280,14 @@ def blend_front(e1st, i1st, e2nd, halo, sdev):
 
     for inum in range(halo):
 
-        print("* blending sweep:", inum)
-
         part = np.minimum(part, part[npos, :] + 1.)
         part = np.minimum(part, part[spos, :] + 1.)
         part = np.minimum(part, part[:, epos] + 1.)
         part = np.minimum(part, part[:, wpos] + 1.)
 
+        print("* compute local blending:",
+              inum, "of", halo)
+
     mask[sidx, :] = part
 
     print("* blending NORTH.")
@@ -89,13 +311,14 @@ def blend_front(e1st, i1st, e2nd, halo, sdev):
 
     for inum in range(halo):
 
-        print("* blending sweep:", inum)
-
         part = np.minimum(part, part[npos, :] + 1.)
         part = np.minimum(part, part[spos, :] + 1.)
         part = np.minimum(part, part[:, epos] + 1.)
         part = np.minimum(part, part[:, wpos] + 1.)
 
+        print("* compute local blending:",
+              inum, "of", halo)
+
     mask[nidx, :] = part
 
     print("* blending final.")
@@ -106,10 +329,11 @@ def blend_front(e1st, i1st, e2nd, halo, sdev):
     mask = mask ** 1.50
     mask[i1st >= 1] = 0.
 
-    return mask
+    return np.asarray(mask, dtype=np.float32)
 
 
 def rtopo_60sec(elev_path, save_path):
+
     """
     Create a zipped and pixel centred version of RTopo 2.0.4
     (60 arc-sec) to support remapping of elevation data.
@@ -121,7 +345,7 @@ def rtopo_60sec(elev_path, save_path):
     data = nc.Dataset(os.path.join(
         elev_path, "RTopo-2.0.4_1min_data.nc"), "r")
 
-    data.set_auto_maskandscale(False)  # quite valid_min/max
+    data.set_auto_maskandscale(False)  # quiet valid_min/max
 
     xpos = np.asarray(data["lon"][:], dtype=np.float64)
     ypos = np.asarray(data["lat"][:], dtype=np.float64)
@@ -199,7 +423,7 @@ def rtopo_30sec(elev_path, save_path):
         elev_path,
         "RTopo-2.0.4_30sec_bedrock_topography.nc"), "r")
 
-    data.set_auto_maskandscale(False)  # quite valid_min/max
+    data.set_auto_maskandscale(False)  # quiet valid_min/max
 
     xpos = np.asarray(data["lon"][:], dtype=np.float64)
     ypos = np.asarray(data["lat"][:], dtype=np.float64)
@@ -211,7 +435,7 @@ def rtopo_30sec(elev_path, save_path):
         elev_path,
         "RTopo-2.0.4_30sec_surface_elevation.nc"), "r")
 
-    data.set_auto_maskandscale(False)  # quite valid_min/max
+    data.set_auto_maskandscale(False)  # quiet valid_min/max
 
     surf = np.asarray(
         data["surface_elevation"][:], dtype=np.float32)
@@ -220,7 +444,7 @@ def rtopo_30sec(elev_path, save_path):
         elev_path,
         "RTopo-2.0.4_30sec_ice_base_topography.nc"), "r")
 
-    data.set_auto_maskandscale(False)  # quite valid_min/max
+    data.set_auto_maskandscale(False)  # quiet valid_min/max
 
     base = np.asarray(
         data["ice_base_topography"][:], dtype=np.float32)
@@ -278,19 +502,183 @@ def rtopo_30sec(elev_path, save_path):
     root.close()
 
 
-def srtmp_60sec(elev_path, save_path):
+def rtopo_15sec(elev_path, save_path):
+    """
+    Create a zipped and pixel centred version of RTopo 2.0.4
+    (15 arc-sec) to support remapping of elevation data.
+
+    """
+
+    print("Making RTopo-2.0.4 (15 arc-sec) pixel...")
+
+    data = nc.Dataset(os.path.join(
+        elev_path,
+        "RTopo-2.0.4_30sec_bedrock_topography.nc"), "r")
+
+    data.set_auto_maskandscale(False)  # quiet valid_min/max
+
+    elev = np.asarray(
+        data["bedrock_topography"][:], dtype=np.float32)
+
+    zmid = (elev[:-1:, :-1:] + elev[+1::, :-1:] +
+            elev[:-1:, +1::] + elev[+1::, +1::]) / 4.
+
+    zlhs = (elev[:-1:, :-1:] + elev[+1::, :-1:]) / 2.
+    zrhs = (elev[:-1:, +1::] + elev[+1::, +1::]) / 2.
+    zbot = (elev[:-1:, :-1:] + elev[:-1:, +1::]) / 2.
+    ztop = (elev[+1::, :-1:] + elev[+1::, +1::]) / 2.
+
+    znew = np.zeros((43200, 86400), dtype=np.float32)
+    znew[0::2, 0::2] = (
+        elev[:-1:, :-1:] + zbot + zmid + zlhs) / 4.
+    znew[0::2, 1::2] = (
+        zbot + elev[:-1:, +1::] + zrhs + zmid) / 4.
+    znew[1::2, 1::2] = (
+        zmid + zrhs + elev[+1::, +1::] + ztop) / 4.
+    znew[1::2, 0::2] = (
+        zlhs + zmid + ztop + elev[+1::, :-1:]) / 4.
+
+    del zmid
+    del zlhs
+    del zrhs
+    del zbot
+    del ztop
+
+    elev = np.asarray(np.round(znew), dtype=np.int16)
+    data.close()
+
+    data = nc.Dataset(os.path.join(
+        elev_path,
+        "RTopo-2.0.4_30sec_surface_elevation.nc"), "r")
+
+    data.set_auto_maskandscale(False)  # quiet valid_min/max
+
+    surf = np.asarray(
+        data["surface_elevation"][:], dtype=np.float32)
+
+    zmid = (surf[:-1:, :-1:] + surf[+1::, :-1:] +
+            surf[:-1:, +1::] + surf[+1::, +1::]) / 4.
+
+    zlhs = (surf[:-1:, :-1:] + surf[+1::, :-1:]) / 2.
+    zrhs = (surf[:-1:, +1::] + surf[+1::, +1::]) / 2.
+    zbot = (surf[:-1:, :-1:] + surf[:-1:, +1::]) / 2.
+    ztop = (surf[+1::, :-1:] + surf[+1::, +1::]) / 2.
+
+    znew = np.zeros((43200, 86400), dtype=np.float32)
+    znew[0::2, 0::2] = (
+        surf[:-1:, :-1:] + zbot + zmid + zlhs) / 4.
+    znew[0::2, 1::2] = (
+        zbot + surf[:-1:, +1::] + zrhs + zmid) / 4.
+    znew[1::2, 1::2] = (
+        zmid + zrhs + surf[+1::, +1::] + ztop) / 4.
+    znew[1::2, 0::2] = (
+        zlhs + zmid + ztop + surf[+1::, :-1:]) / 4.
+
+    del zmid
+    del zlhs
+    del zrhs
+    del zbot
+    del ztop
+
+    surf = np.asarray(np.round(znew), dtype=np.int16)
+    data.close()
+
+    data = nc.Dataset(os.path.join(
+        elev_path,
+        "RTopo-2.0.4_30sec_ice_base_topography.nc"), "r")
+
+    data.set_auto_maskandscale(False)  # quiet valid_min/max
+
+    base = np.asarray(
+        data["ice_base_topography"][:], dtype=np.float32)
+
+    zmid = (base[:-1:, :-1:] + base[+1::, :-1:] +
+            base[:-1:, +1::] + base[+1::, +1::]) / 4.
+
+    zlhs = (base[:-1:, :-1:] + base[+1::, :-1:]) / 2.
+    zrhs = (base[:-1:, +1::] + base[+1::, +1::]) / 2.
+    zbot = (base[:-1:, :-1:] + base[:-1:, +1::]) / 2.
+    ztop = (base[+1::, :-1:] + base[+1::, +1::]) / 2.
+
+    znew = np.zeros((43200, 86400), dtype=np.float32)
+    znew[0::2, 0::2] = (
+        base[:-1:, :-1:] + zbot + zmid + zlhs) / 4.
+    znew[0::2, 1::2] = (
+        zbot + base[:-1:, +1::] + zrhs + zmid) / 4.
+    znew[1::2, 1::2] = (
+        zmid + zrhs + base[+1::, +1::] + ztop) / 4.
+    znew[1::2, 0::2] = (
+        zlhs + zmid + ztop + base[+1::, :-1:]) / 4.
+
+    del zmid
+    del zlhs
+    del zrhs
+    del zbot
+    del ztop
+
+    base = np.asarray(np.round(znew), dtype=np.int16)
+    data.close()
+
+    iceh = surf - base
+    iceh[base == 0] = 0
+
+    ocnh = np.maximum(0, base - elev)
+
+    xpos = np.linspace(
+        -180., +180., elev.shape[1] + 1, dtype=np.float64)
+    ypos = np.linspace(
+        -90.0, +90.0, elev.shape[0] + 1, dtype=np.float64)
+
+    root = nc.Dataset(
+        os.path.join(
+            save_path, "RTopo_2_0_4_15sec_pixel.nc"),
+        "w", format="NETCDF4")
+    root.description = "A zipped RTopo-2.0.4 (15 arc-sec) " + \
+        "data-set, pixel centred and compressed to int16_t."
+    root.source = "RTopo-2.0.4_30sec_data.nc"
+    root.references = "doi.pangaea.de/10.1594/PANGAEA.905295"
+    root.createDimension("num_lon", elev.shape[1] + 1)
+    root.createDimension("num_col", elev.shape[1])
+    root.createDimension("num_lat", elev.shape[0] + 1)
+    root.createDimension("num_row", elev.shape[0])
+
+    data = root.createVariable("lon", "f8", ("num_lon"))
+    data.units = "degrees_east"
+    data[:] = xpos
+    data = root.createVariable("lat", "f8", ("num_lat"))
+    data.units = "degrees_north"
+    data[:] = ypos
+    data = root.createVariable(
+        "bed_elevation", "i2", ("num_row", "num_col"))
+    data.units = "m"
+    data[:, :] = elev
+
+    data = root.createVariable(
+        "ice_thickness", "i2", ("num_row", "num_col"))
+    data.units = "m"
+    data[:, :] = iceh
+    data = root.createVariable(
+        "ocn_thickness", "i2", ("num_row", "num_col"))
+    data.units = "m"
+    data[:, :] = ocnh
+
+    root.close()
+
+
+def gebco_60sec(elev_path, save_path):
     """
-    Create a zipped and pixel centred version of SRTM15+V2.1
+    Create a zipped and pixel centred version of GEBCO[2023]
     (15 arc-sec) at 60 arc-sec.
 
     """
 
-    print("Making SRTM15+V2.1 (60 arc-sec) pixel...")
+    print("Making GEBCO[2023] (60 arc-sec) pixel...")
 
     data = nc.Dataset(os.path.join(
-        elev_path, "SRTM15+V2.1.nc"), "r")
+        elev_path, "GEBCO_2023_sub_ice_topo.nc"), "r")
 
-    elev = np.asarray(data["z"][:], dtype=np.float32)
+    elev = np.asarray(
+        data["elevation"][:], dtype=np.float32)
 
     halo = 4
     z_60 = np.zeros((10800, 21600), dtype=np.float32)
@@ -311,7 +699,7 @@ def srtmp_60sec(elev_path, save_path):
 
     root = nc.Dataset(
         os.path.join(
-            save_path, "SRTM15+V2.1_60sec_pixel.nc"),
+            save_path, "GEBCO_v2023_60sec_pixel.nc"),
         "w", format="NETCDF4")
     root.createDimension("num_lon", elev.shape[1] + 1)
     root.createDimension("num_col", elev.shape[1])
@@ -325,26 +713,27 @@ def srtmp_60sec(elev_path, save_path):
     data.units = "degrees_north"
     data[:] = ypos
     data = root.createVariable(
-        "top_elevation", "i2", ("num_row", "num_col"))
+        "bed_elevation", "i2", ("num_row", "num_col"))
     data.units = "m"
     data[:, :] = elev
 
     root.close()
 
 
-def srtmp_30sec(elev_path, save_path):
+def gebco_30sec(elev_path, save_path):
     """
-    Create a zipped and pixel centred version of SRTM15+V2.1
+    Create a zipped and pixel centred version of GEBCO[2023]
     (15 arc-sec) at 30 arc-sec.
 
     """
 
-    print("Making SRTM15+V2.1 (30 arc-sec) pixel...")
+    print("Making GEBCO[2023] (30 arc-sec) pixel...")
 
     data = nc.Dataset(os.path.join(
-        elev_path, "SRTM15+V2.1.nc"), "r")
+        elev_path, "GEBCO_2023_sub_ice_topo.nc"), "r")
 
-    elev = np.asarray(data["z"][:], dtype=np.float32)
+    elev = np.asarray(
+        data["elevation"][:], dtype=np.float32)
 
     halo = 2
     z_30 = np.zeros((21600, 43200), dtype=np.float32)
@@ -365,7 +754,7 @@ def srtmp_30sec(elev_path, save_path):
 
     root = nc.Dataset(
         os.path.join(
-            save_path, "SRTM15+V2.1_30sec_pixel.nc"),
+            save_path, "GEBCO_v2023_30sec_pixel.nc"),
         "w", format="NETCDF4")
     root.createDimension("num_lon", elev.shape[1] + 1)
     root.createDimension("num_col", elev.shape[1])
@@ -379,48 +768,27 @@ def srtmp_30sec(elev_path, save_path):
     data.units = "degrees_north"
     data[:] = ypos
     data = root.createVariable(
-        "top_elevation", "i2", ("num_row", "num_col"))
+        "bed_elevation", "i2", ("num_row", "num_col"))
     data.units = "m"
     data[:, :] = elev
 
     root.close()
 
 
-def rtopo_srtmp_60sec(elev_path, save_path):
+def gebco_15sec(elev_path, save_path):
     """
-    Create a zipped and pixel centred 'blend' of RTopo 2.0.4
-    and SRTM15+V2.1 at 60 arc-sec.
+    Create a zipped and pixel centred version of GEBCO[2023]
+    (15 arc-sec) at 15 arc-sec.
 
     """
 
-    print("Making RTopo-SRTM+ (60 arc-sec) blend...")
-
-    data = nc.Dataset(os.path.join(
-        elev_path, "RTopo_2_0_4_60sec_pixel.nc"), "r")
-
-    e1st = np.asarray(
-        data["bed_elevation"][:], dtype=np.int16)
-
-    i1st = np.asarray(
-        data["ice_thickness"][:], dtype=np.int16)
-
-    o1st = np.asarray(
-        data["ocn_thickness"][:], dtype=np.int16)
+    print("Making GEBCO[2023] (15 arc-sec) pixel...")
 
     data = nc.Dataset(os.path.join(
-        elev_path, "SRTM15+V2.1_60sec_pixel.nc"), "r")
-
-    e2nd = np.asarray(
-        data["top_elevation"][:], dtype=np.int16)
+        elev_path, "GEBCO_2023_sub_ice_topo.nc"), "r")
 
-    mask = blend_front(e1st, i1st, e2nd, halo=20, sdev=2.0)
-
-    elev = np.asarray(np.round(
-        (1. - mask) * e1st + mask * e2nd), dtype=np.int16)
-
-    iceh = i1st
-    ocnh = o1st
-    ocnh[i1st == 0] = np.maximum(0, -elev[i1st == 0])
+    elev = np.asarray(
+        data["elevation"][:], dtype=np.int16)
 
     xpos = np.linspace(
         -180., +180., elev.shape[1] + 1, dtype=np.float64)
@@ -428,18 +796,9 @@ def rtopo_srtmp_60sec(elev_path, save_path):
         -90.0, +90.0, elev.shape[0] + 1, dtype=np.float64)
 
     root = nc.Dataset(
-        os.path.join(save_path,
-                     "RTopo_2_0_4_SRTM15+V2_1_60sec_pixel.nc"),
+        os.path.join(
+            save_path, "GEBCO_v2023_15sec_pixel.nc"),
         "w", format="NETCDF4")
-    root.description = "Blend of RTopo-2.0.4 (60 arc-sec) " + \
-        "and SRTM15+V2.1 (60 arc-sec) - pixel centred and " + \
-        "compressed to int16_t. RTopo data used under ice " + \
-        "sheets/shelves."
-    root.source = \
-        "RTopo-2.0.4_1min_data.nc and SRTM15+V2.1.nc"
-    root.references = \
-        "doi.pangaea.de/10.1594/PANGAEA.905295 and " + \
-        "doi.org/10.1029/2019EA000658"
     root.createDimension("num_lon", elev.shape[1] + 1)
     root.createDimension("num_col", elev.shape[1])
     root.createDimension("num_lat", elev.shape[0] + 1)
@@ -455,29 +814,21 @@ def rtopo_srtmp_60sec(elev_path, save_path):
         "bed_elevation", "i2", ("num_row", "num_col"))
     data.units = "m"
     data[:, :] = elev
-    data = root.createVariable(
-        "ice_thickness", "i2", ("num_row", "num_col"))
-    data.units = "m"
-    data[:, :] = iceh
-    data = root.createVariable(
-        "ocn_thickness", "i2", ("num_row", "num_col"))
-    data.units = "m"
-    data[:, :] = ocnh
 
     root.close()
 
 
-def rtopo_srtmp_30sec(elev_path, save_path):
+def rtopo_gebco_60sec(elev_path, save_path):
     """
     Create a zipped and pixel centred 'blend' of RTopo 2.0.4
-    and SRTM15+V2.1 at 30 arc-sec.
+    and GEBCO[2023] at 60 arc-sec.
 
     """
 
-    print("Making RTopo-SRTM+ (30 arc-sec) blend...")
+    print("Making RTopo-GEBCO (60 arc-sec) blend...")
 
     data = nc.Dataset(os.path.join(
-        elev_path, "RTopo_2_0_4_30sec_pixel.nc"), "r")
+        elev_path, "RTopo_2_0_4_60sec_pixel.nc"), "r")
 
     e1st = np.asarray(
         data["bed_elevation"][:], dtype=np.int16)
@@ -489,12 +840,12 @@ def rtopo_srtmp_30sec(elev_path, save_path):
         data["ocn_thickness"][:], dtype=np.int16)
 
     data = nc.Dataset(os.path.join(
-        elev_path, "SRTM15+V2.1_30sec_pixel.nc"), "r")
+        elev_path, "GEBCO_v2023_60sec_pixel.nc"), "r")
 
     e2nd = np.asarray(
-        data["top_elevation"][:], dtype=np.int16)
+        data["bed_elevation"][:], dtype=np.int16)
 
-    mask = blend_front(e1st, i1st, e2nd, halo=40, sdev=4.0)
+    mask = blend_front(e1st, i1st, e2nd, halo=10, sdev=1.0)
 
     elev = np.asarray(np.round(
         (1. - mask) * e1st + mask * e2nd), dtype=np.int16)
@@ -510,17 +861,17 @@ def rtopo_srtmp_30sec(elev_path, save_path):
 
     root = nc.Dataset(
         os.path.join(save_path,
-                     "RTopo_2_0_4_SRTM15+V2_1_30sec_pixel.nc"),
+                     "RTopo_2_0_4_GEBCO_v2023_60sec_pixel.nc"),
         "w", format="NETCDF4")
-    root.description = "Blend of RTopo-2.0.4 (30 arc-sec) " + \
-        "and SRTM15+V2.1 (30 arc-sec) - pixel centred and " + \
+    root.description = "Blend of RTopo-2.0.4 (60 arc-sec) " + \
+        "and GEBCO[2023] (15 arc-sec) - pixel centred and " + \
         "compressed to int16_t. RTopo data used under ice " + \
-        "sheets/shelves."
+        "sheets/shelves. Remapped to 60 arc-sec."
     root.source = \
-        "RTopo-2.0.4_30sec_data.nc and SRTM15+V2.1.nc"
+        "RTopo-2.0.4_1min_data.nc and GEBCO_2023.nc"
     root.references = \
         "doi.pangaea.de/10.1594/PANGAEA.905295 and " + \
-        "doi.org/10.1029/2019EA000658"
+        "doi.org/10.5285/a29c5465-b138-234d-e053-6c86abc040b9"
     root.createDimension("num_lon", elev.shape[1] + 1)
     root.createDimension("num_col", elev.shape[1])
     root.createDimension("num_lat", elev.shape[0] + 1)
@@ -535,140 +886,53 @@ def rtopo_srtmp_30sec(elev_path, save_path):
     data = root.createVariable(
         "bed_elevation", "i2", ("num_row", "num_col"))
     data.units = "m"
+    data.long_name = names.bed_elevation
     data[:, :] = elev
     data = root.createVariable(
         "ice_thickness", "i2", ("num_row", "num_col"))
     data.units = "m"
+    data.long_name = names.ocn_thickness
     data[:, :] = iceh
     data = root.createVariable(
         "ocn_thickness", "i2", ("num_row", "num_col"))
     data.units = "m"
+    data.long_name = names.ice_thickness
     data[:, :] = ocnh
 
-    root.close()
+    # filt. grid-scale noise that imprints on dz/dx...
+    filt = gaussian_filter(np.asarray(
+        elev, dtype=np.float32), sigma=.625, mode="wrap")
 
+    zslp, dzdx, dzdy = \
+        cell_dzdx(xpos, ypos, filt, RSPH)
 
-def gebco_60sec(elev_path, save_path):
-    """
-    Create a zipped and pixel centred version of GEBCO[2020]
-    (15 arc-sec) at 60 arc-sec.
-
-    """
-
-    print("Making GEBCO[2020] (60 arc-sec) pixel...")
-
-    data = nc.Dataset(os.path.join(
-        elev_path, "GEBCO_2020.nc"), "r")
-
-    elev = np.asarray(
-        data["elevation"][:], dtype=np.float32)
-
-    halo = 4
-    z_60 = np.zeros((10800, 21600), dtype=np.float32)
-
-    for ipos in range(halo):
-        for jpos in range(halo):
-            iend = elev.shape[0] - halo + ipos + 1
-            jend = elev.shape[1] - halo + jpos + 1
-            z_60 += elev[ipos:iend:halo, jpos:jend:halo]
-
-    elev = np.asarray(
-        np.round(z_60 / float(halo ** 2)), dtype=np.int16)
-
-    xpos = np.linspace(
-        -180., +180., elev.shape[1] + 1, dtype=np.float64)
-    ypos = np.linspace(
-        -90.0, +90.0, elev.shape[0] + 1, dtype=np.float64)
-
-    root = nc.Dataset(
-        os.path.join(
-            save_path, "GEBCO_v2020_60sec_pixel.nc"),
-        "w", format="NETCDF4")
-    root.createDimension("num_lon", elev.shape[1] + 1)
-    root.createDimension("num_col", elev.shape[1])
-    root.createDimension("num_lat", elev.shape[0] + 1)
-    root.createDimension("num_row", elev.shape[0])
-
-    data = root.createVariable("lon", "f8", ("num_lon"))
-    data.units = "degrees_east"
-    data[:] = xpos
-    data = root.createVariable("lat", "f8", ("num_lat"))
-    data.units = "degrees_north"
-    data[:] = ypos
     data = root.createVariable(
-        "top_elevation", "i2", ("num_row", "num_col"))
-    data.units = "m"
-    data[:, :] = elev
-
-    root.close()
-
-
-def gebco_30sec(elev_path, save_path):
-    """
-    Create a zipped and pixel centred version of GEBCO[2020]
-    (15 arc-sec) at 30 arc-sec.
-
-    """
-
-    print("Making GEBCO[2020] (30 arc-sec) pixel...")
-
-    data = nc.Dataset(os.path.join(
-        elev_path, "GEBCO_2020.nc"), "r")
-
-    elev = np.asarray(
-        data["elevation"][:], dtype=np.float32)
-
-    halo = 2
-    z_30 = np.zeros((21600, 43200), dtype=np.float32)
-
-    for ipos in range(halo):
-        for jpos in range(halo):
-            iend = elev.shape[0] - halo + ipos + 1
-            jend = elev.shape[1] - halo + jpos + 1
-            z_30 += elev[ipos:iend:halo, jpos:jend:halo]
-
-    elev = np.asarray(
-        np.round(z_30 / float(halo ** 2)), dtype=np.int16)
-
-    xpos = np.linspace(
-        -180., +180., elev.shape[1] + 1, dtype=np.float64)
-    ypos = np.linspace(
-        -90.0, +90.0, elev.shape[0] + 1, dtype=np.float64)
-
-    root = nc.Dataset(
-        os.path.join(
-            save_path, "GEBCO_v2020_30sec_pixel.nc"),
-        "w", format="NETCDF4")
-    root.createDimension("num_lon", elev.shape[1] + 1)
-    root.createDimension("num_col", elev.shape[1])
-    root.createDimension("num_lat", elev.shape[0] + 1)
-    root.createDimension("num_row", elev.shape[0])
-
-    data = root.createVariable("lon", "f8", ("num_lon"))
-    data.units = "degrees_east"
-    data[:] = xpos
-    data = root.createVariable("lat", "f8", ("num_lat"))
-    data.units = "degrees_north"
-    data[:] = ypos
+        "bed_slope", "f4", ("num_row", "num_col"))
+    data.long_name = names.bed_slope
+    data[:, :] = zslp[:, :]
     data = root.createVariable(
-        "top_elevation", "i2", ("num_row", "num_col"))
-    data.units = "m"
-    data[:, :] = elev
+        "bed_dz_dx", "f4", ("num_row", "num_col"))
+    data.long_name = names.bed_dz_dx
+    data[:, :] = dzdx[:, :]
+    data = root.createVariable(
+        "bed_dz_dy", "f4", ("num_row", "num_col"))
+    data.long_name = names.bed_dz_dy
+    data[:, :] = dzdy[:, :]
 
     root.close()
 
 
-def rtopo_gebco_60sec(elev_path, save_path):
+def rtopo_gebco_30sec(elev_path, save_path):
     """
     Create a zipped and pixel centred 'blend' of RTopo 2.0.4
-    and GEBCO[2020] at 60 arc-sec.
+    and GEBCO[2023] at 30 arc-sec.
 
     """
 
-    print("Making RTopo-GEBCO (60 arc-sec) blend...")
+    print("Making RTopo-GEBCO (30 arc-sec) blend...")
 
     data = nc.Dataset(os.path.join(
-        elev_path, "RTopo_2_0_4_60sec_pixel.nc"), "r")
+        elev_path, "RTopo_2_0_4_30sec_pixel.nc"), "r")
 
     e1st = np.asarray(
         data["bed_elevation"][:], dtype=np.int16)
@@ -680,12 +944,12 @@ def rtopo_gebco_60sec(elev_path, save_path):
         data["ocn_thickness"][:], dtype=np.int16)
 
     data = nc.Dataset(os.path.join(
-        elev_path, "GEBCO_v2020_60sec_pixel.nc"), "r")
+        elev_path, "GEBCO_v2023_30sec_pixel.nc"), "r")
 
     e2nd = np.asarray(
-        data["top_elevation"][:], dtype=np.int16)
+        data["bed_elevation"][:], dtype=np.int16)
 
-    mask = blend_front(e1st, i1st, e2nd, halo=20, sdev=2.0)
+    mask = blend_front(e1st, i1st, e2nd, halo=10, sdev=1.0)
 
     elev = np.asarray(np.round(
         (1. - mask) * e1st + mask * e2nd), dtype=np.int16)
@@ -701,14 +965,14 @@ def rtopo_gebco_60sec(elev_path, save_path):
 
     root = nc.Dataset(
         os.path.join(save_path,
-                     "RTopo_2_0_4_GEBCO_v2020_60sec_pixel.nc"),
+                     "RTopo_2_0_4_GEBCO_v2023_30sec_pixel.nc"),
         "w", format="NETCDF4")
-    root.description = "Blend of RTopo-2.0.4 (60 arc-sec) " + \
-        "and GEBCO[2020] (15 arc-sec) - pixel centred and " + \
+    root.description = "Blend of RTopo-2.0.4 (30 arc-sec) " + \
+        "and GEBCO[2023] (15 arc-sec) - pixel centred and " + \
         "compressed to int16_t. RTopo data used under ice " + \
-        "sheets/shelves. Remapped to 60 arc-sec."
+        "sheets/shelves. Remapped to 30 arc-sec."
     root.source = \
-        "RTopo-2.0.4_1min_data.nc and GEBCO_2020.nc"
+        "RTopo-2.0.4_30sec_data.nc and GEBCO_2023.nc"
     root.references = \
         "doi.pangaea.de/10.1594/PANGAEA.905295 and " + \
         "doi.org/10.5285/a29c5465-b138-234d-e053-6c86abc040b9"
@@ -726,54 +990,88 @@ def rtopo_gebco_60sec(elev_path, save_path):
     data = root.createVariable(
         "bed_elevation", "i2", ("num_row", "num_col"))
     data.units = "m"
+    data.long_name = names.bed_elevation
     data[:, :] = elev
     data = root.createVariable(
         "ice_thickness", "i2", ("num_row", "num_col"))
     data.units = "m"
+    data.long_name = names.ocn_thickness
     data[:, :] = iceh
     data = root.createVariable(
         "ocn_thickness", "i2", ("num_row", "num_col"))
     data.units = "m"
+    data.long_name = names.ice_thickness
     data[:, :] = ocnh
 
+    # filt. grid-scale noise that imprints on dz/dx...
+    filt = gaussian_filter(np.asarray(
+        elev, dtype=np.float32), sigma=1.25, mode="wrap")
+
+    zslp, dzdx, dzdy = \
+        cell_dzdx(xpos, ypos, filt, RSPH)
+
+    data = root.createVariable(
+        "bed_slope", "f4", ("num_row", "num_col"))
+    data.long_name = names.bed_slope
+    data[:, :] = zslp[:, :]
+    data = root.createVariable(
+        "bed_dz_dx", "f4", ("num_row", "num_col"))
+    data.long_name = names.bed_dz_dx
+    data[:, :] = dzdx[:, :]
+    data = root.createVariable(
+        "bed_dz_dy", "f4", ("num_row", "num_col"))
+    data.long_name = names.bed_dz_dy
+    data[:, :] = dzdy[:, :]
+
     root.close()
 
 
-def rtopo_gebco_30sec(elev_path, save_path):
+def rtopo_gebco_15sec(elev_path, save_path):
     """
     Create a zipped and pixel centred 'blend' of RTopo 2.0.4
-    and GEBCO[2020] at 30 arc-sec.
+    and GEBCO[2023] at 15 arc-sec.
 
     """
 
-    print("Making RTopo-GEBCO (30 arc-sec) blend...")
+    print("Making RTopo-GEBCO (15 arc-sec) blend...")
 
     data = nc.Dataset(os.path.join(
-        elev_path, "RTopo_2_0_4_30sec_pixel.nc"), "r")
+        elev_path, "RTopo_2_0_4_15sec_pixel.nc"), "r")
 
-    e1st = np.asarray(
+    elev = np.asarray(
         data["bed_elevation"][:], dtype=np.int16)
 
-    i1st = np.asarray(
+    iceh = np.asarray(
         data["ice_thickness"][:], dtype=np.int16)
 
-    o1st = np.asarray(
-        data["ocn_thickness"][:], dtype=np.int16)
-
     data = nc.Dataset(os.path.join(
-        elev_path, "GEBCO_v2020_30sec_pixel.nc"), "r")
+        elev_path, "GEBCO_v2023_15sec_pixel.nc"), "r")
 
     e2nd = np.asarray(
-        data["top_elevation"][:], dtype=np.int16)
+        data["bed_elevation"][:], dtype=np.int16)
 
-    mask = blend_front(e1st, i1st, e2nd, halo=40, sdev=4.0)
+    mask = blend_front(elev, iceh, e2nd, halo=20, sdev=2.0)
 
-    elev = np.asarray(np.round(
-        (1. - mask) * e1st + mask * e2nd), dtype=np.int16)
+    # careful w mem. alloc.
+    del iceh
+    elev = np.asarray(elev, dtype=np.float32)
+    elev -= mask * elev
+    elev += mask * e2nd
+    del e2nd
+    del mask
 
-    iceh = i1st
-    ocnh = o1st
-    ocnh[i1st == 0] = np.maximum(0, -elev[i1st == 0])
+    elev = np.asarray(np.round(elev), dtype=np.int16)
+
+    data = nc.Dataset(os.path.join(
+        elev_path, "RTopo_2_0_4_15sec_pixel.nc"), "r")
+
+    iceh = np.asarray(
+        data["ice_thickness"][:], dtype=np.int16)
+
+    ocnh = np.asarray(
+        data["ocn_thickness"][:], dtype=np.int16)
+
+    ocnh[iceh == 0] = np.maximum(0, -elev[iceh == 0])
 
     xpos = np.linspace(
         -180., +180., elev.shape[1] + 1, dtype=np.float64)
@@ -782,14 +1080,14 @@ def rtopo_gebco_30sec(elev_path, save_path):
 
     root = nc.Dataset(
         os.path.join(save_path,
-                     "RTopo_2_0_4_GEBCO_v2020_30sec_pixel.nc"),
+                     "RTopo_2_0_4_GEBCO_v2023_15sec_pixel.nc"),
         "w", format="NETCDF4")
     root.description = "Blend of RTopo-2.0.4 (30 arc-sec) " + \
-        "and GEBCO[2020] (15 arc-sec) - pixel centred and " + \
+        "and GEBCO[2023] (15 arc-sec) - pixel centred and " + \
         "compressed to int16_t. RTopo data used under ice " + \
-        "sheets/shelves. Remapped to 30 arc-sec."
+        "sheets/shelves. Remapped to 15 arc-sec."
     root.source = \
-        "RTopo-2.0.4_30sec_data.nc and GEBCO_2020.nc"
+        "RTopo-2.0.4_30sec_data.nc and GEBCO_2023.nc"
     root.references = \
         "doi.pangaea.de/10.1594/PANGAEA.905295 and " + \
         "doi.org/10.5285/a29c5465-b138-234d-e053-6c86abc040b9"
@@ -807,16 +1105,43 @@ def rtopo_gebco_30sec(elev_path, save_path):
     data = root.createVariable(
         "bed_elevation", "i2", ("num_row", "num_col"))
     data.units = "m"
+    data.long_name = names.bed_elevation
     data[:, :] = elev
     data = root.createVariable(
         "ice_thickness", "i2", ("num_row", "num_col"))
     data.units = "m"
+    data.long_name = names.ice_thickness
     data[:, :] = iceh
     data = root.createVariable(
         "ocn_thickness", "i2", ("num_row", "num_col"))
     data.units = "m"
+    data.long_name = names.ocn_thickness
     data[:, :] = ocnh
 
+    # filt. grid-scale noise that imprints on dz/dx...
+    filt = gaussian_filter(np.asarray(
+        elev, dtype=np.float32), sigma=2.50, mode="wrap")
+
+    del elev
+    del ocnh
+    del iceh
+
+    zslp, dzdx, dzdy = \
+        cell_dzdx(xpos, ypos, filt, RSPH)
+
+    data = root.createVariable(
+        "bed_slope", "f4", ("num_row", "num_col"))
+    data.long_name = names.bed_slope
+    data[:, :] = zslp
+    data = root.createVariable(
+        "bed_dz_dx", "f4", ("num_row", "num_col"))
+    data.long_name = names.bed_dz_dx
+    data[:, :] = dzdx
+    data = root.createVariable(
+        "bed_dz_dy", "f4", ("num_row", "num_col"))
+    data.long_name = names.bed_dz_dy
+    data[:, :] = dzdy
+
     root.close()
 
 
@@ -841,20 +1166,12 @@ def rtopo_gebco_30sec(elev_path, save_path):
 
     rtopo_60sec(elev_path, save_path)
     rtopo_30sec(elev_path, save_path)
-
-    """
-    SRTM15+V2.1 data seems to include high-freq 'noise'
-    near coastlines, so don't use for now...
-
-    srtmp_60sec(elev_path, save_path)
-    srtmp_30sec(elev_path, save_path)
-
-    rtopo_srtmp_60sec(elev_path, save_path)
-    rtopo_srtmp_30sec(elev_path, save_path)
-    """
+    rtopo_15sec(elev_path, save_path)
 
     gebco_60sec(elev_path, save_path)
     gebco_30sec(elev_path, save_path)
+    gebco_15sec(elev_path, save_path)
 
     rtopo_gebco_60sec(elev_path, save_path)
     rtopo_gebco_30sec(elev_path, save_path)
+    rtopo_gebco_15sec(elev_path, save_path)
diff --git a/compass/ocean/tests/tides/dem/dem_remap.py b/compass/ocean/tests/tides/dem/dem_remap.py
new file mode 100644
index 0000000000..de96144a08
--- /dev/null
+++ b/compass/ocean/tests/tides/dem/dem_remap.py
@@ -0,0 +1,818 @@
+
+import argparse
+import time
+
+import netCDF4 as nc
+import numpy as np
+from scipy import spatial
+from scipy.sparse import csr_matrix
+
+from compass.ocean.tests.tides.dem.dem_names import names
+
+# Authors: Darren Engwirda
+
+
+def map_to_r3(mesh, xlon, ylat, head, tail):
+    """
+    Map lon-lat coordinates to XYZ points. Restricted to the
+    panel LAT[HEAD:TAIL] to manage memory use.
+
+    """
+
+    sinx = np.sin(xlon * np.pi / 180.)
+    cosx = np.cos(xlon * np.pi / 180.)
+    siny = np.sin(ylat * np.pi / 180.)
+    cosy = np.cos(ylat * np.pi / 180.)
+
+    sinu, sinv = np.meshgrid(sinx, siny[head:tail])
+    cosu, cosv = np.meshgrid(cosx, cosy[head:tail])
+
+    rsph = mesh.sphere_radius
+
+    xpos = rsph * cosu * cosv
+    ypos = rsph * sinu * cosv
+    zpos = rsph * sinv
+
+    return np.vstack(
+        (xpos.ravel(), ypos.ravel(), zpos.ravel())).T
+
+
+def tria_area(pa, pb, pc, rs=1.):
+    """
+    Calculate areas of spherical triangles [PA, PB, PC] on a
+    sphere of radius RS.
+
+    """
+
+    lena = circ_dist(
+        pa[:, 0], pa[:, 1], pb[:, 0], pb[:, 1], 1.)
+    lenb = circ_dist(
+        pb[:, 0], pb[:, 1], pc[:, 0], pc[:, 1], 1.)
+    lenc = circ_dist(
+        pc[:, 0], pc[:, 1], pa[:, 0], pa[:, 1], 1.)
+
+    slen = 0.5 * (lena + lenb + lenc)
+
+    tana = np.tan(0.5 * (slen - lena))
+    tanb = np.tan(0.5 * (slen - lenb))
+    tanc = np.tan(0.5 * (slen - lenc))
+
+    edel = 4.0 * np.arctan(np.sqrt(
+        np.tan(0.5 * slen) * tana * tanb * tanc))
+
+    return edel * rs ** 2
+
+
+def circ_dist(xa, ya, xb, yb, rs=1.):
+    """
+    Calculate geodesic-length of great circles [PA, PB] on a
+    sphere of radius RS.
+
+    """
+
+    dlon = .5 * (xa - xb)
+    dlat = .5 * (ya - yb)
+
+    dist = 2. * rs * np.arcsin(np.sqrt(
+        np.sin(dlat) ** 2 +
+        np.sin(dlon) ** 2 * np.cos(ya) * np.cos(yb)
+    ))
+
+    return dist
+
+
+def sample_1d(xpos, ffun, xnew):
+    """
+    A fast(er) 1-dim. interpolation routine for DEM sampling.
+
+    """
+
+    ipos = np.searchsorted(xpos, xnew, side="left")
+
+    ipos = np.maximum(ipos, 1)
+    ipos = np.minimum(ipos, xpos.size - 1)
+
+    spac = np.diff(xpos)
+    scal = (xnew - xpos[ipos - 1]) / spac[ipos - 1]
+
+    fnew = (
+        (1.0 - scal) * ffun[ipos - 1] +
+        (0.0 + scal) * ffun[ipos - 0]
+    )
+
+    return fnew
+
+
+def linear_2d(xxww, yyss, xxee, yynn,
+              xpos, ypos,
+              ffsw, ffse, ffne, ffnw):
+    """
+    Standard bilinear interpolation on aligned quadrilateral.
+
+    """
+
+    aane = (xxee - xpos) * (yynn - ypos)
+    aanw = (xpos - xxww) * (yynn - ypos)
+    aase = (xxee - xpos) * (ypos - yyss)
+    aasw = (xpos - xxww) * (ypos - yyss)
+
+    asum = (aane + aanw + aase + aasw)
+
+    return (ffnw * aase + ffne * aasw +
+            ffsw * aane + ffse * aanw) / asum
+
+
+def sample_2d(xlon, ylat, vals, xpos, ypos):
+    """
+    A fast(er) 2-dim. interpolation routine for DEM sampling.
+
+    """
+
+    cols = xlon.size - 1
+    rows = ylat.size - 1
+
+    dlon = (xlon[-1] - xlon[+0]) / cols
+    dlat = (ylat[-1] - ylat[+0]) / rows
+
+    icol = (xpos - np.min(xlon)) / dlon
+    irow = (ypos - np.min(ylat)) / dlat
+
+    icol = np.asarray(icol, dtype=int)
+    irow = np.asarray(irow, dtype=int)
+
+    xmid = .5 * (xlon[:-1:] + xlon[1::])
+    ymid = .5 * (ylat[:-1:] + ylat[1::])
+
+    cols = cols - 1
+    rows = rows - 1
+    zero = 0
+    wcol = icol - 1
+    wcol[wcol < zero] = cols
+    ecol = icol + 1
+    ecol[ecol > cols] = zero
+    nrow = irow + 1
+    nrow[nrow > rows] = rows
+    srow = irow - 1
+    srow[srow < zero] = zero
+
+    # -- Sub-pixel bilinear interpolation - pixel DEM values are
+    # -- considered to be located at pixel centres; edge middle
+    # -- and corner values are reconstructed using simple linear
+    # -- averaging. Each pixel is thus split into 4 "sub-pixels"
+    # -- with standard bilinear interpolation applied for each.
+
+    # -- Sub-pixel scheme preserves min/max DEM values - without
+    # -- exessive smoothing.
+
+    vmid = vals[irow, icol]
+
+    vvnw = 0.25 * (
+        vals[irow, icol] + vals[irow, wcol] +
+        vals[nrow, icol] + vals[nrow, wcol]
+    )
+    vvne = 0.25 * (
+        vals[irow, icol] + vals[irow, ecol] +
+        vals[nrow, icol] + vals[nrow, ecol]
+    )
+    vvsw = 0.25 * (
+        vals[irow, icol] + vals[irow, wcol] +
+        vals[srow, icol] + vals[srow, wcol]
+    )
+    vvse = 0.25 * (
+        vals[irow, icol] + vals[irow, ecol] +
+        vals[srow, icol] + vals[srow, ecol]
+    )
+
+    vvnn = 0.50 * (
+        vals[irow, icol] + vals[nrow, icol]
+    )
+    vvee = 0.50 * (
+        vals[irow, icol] + vals[irow, ecol]
+    )
+    vvss = 0.50 * (
+        vals[irow, icol] + vals[srow, icol]
+    )
+    vvww = 0.50 * (
+        vals[irow, icol] + vals[irow, wcol]
+    )
+
+    isnw = np.logical_and(
+        xpos <= xmid[icol + 0], ypos >= ymid[irow + 0]
+    )
+    isne = np.logical_and(
+        xpos >= xmid[icol + 0], ypos >= ymid[irow + 0]
+    )
+    issw = np.logical_and(
+        xpos <= xmid[icol + 0], ypos <= ymid[irow + 0]
+    )
+    isse = np.logical_and(
+        xpos >= xmid[icol + 0], ypos <= ymid[irow + 0]
+    )
+
+    vnew = np.zeros(vmid.shape, dtype=vals.dtype)
+
+    vnew[isnw] = linear_2d(
+        xlon[icol[isnw] + 0], ymid[irow[isnw] + 0],
+        xmid[icol[isnw] + 0], ylat[irow[isnw] + 1],
+        xpos[isnw], ypos[isnw],
+        vvww[isnw], vmid[isnw], vvnn[isnw], vvnw[isnw]
+    )
+
+    vnew[isne] = linear_2d(
+        xmid[icol[isne] + 0], ymid[irow[isne] + 0],
+        xlon[icol[isne] + 1], ylat[irow[isne] + 1],
+        xpos[isne], ypos[isne],
+        vmid[isne], vvee[isne], vvne[isne], vvnn[isne]
+    )
+
+    vnew[issw] = linear_2d(
+        xlon[icol[issw] + 0], ylat[irow[issw] + 0],
+        xmid[icol[issw] + 0], ymid[irow[issw] + 0],
+        xpos[issw], ypos[issw],
+        vvsw[issw], vvss[issw], vmid[issw], vvww[issw]
+    )
+
+    vnew[isse] = linear_2d(
+        xmid[icol[isse] + 0], ylat[irow[isse] + 0],
+        xlon[icol[isse] + 1], ymid[irow[isse] + 0],
+        xpos[isse], ypos[isse],
+        vvss[isse], vvse[isse], vvee[isse], vmid[isse]
+    )
+
+    return vnew
+
+
+def cell_quad(mesh, xlon, ylat, vals):
+    """
+    Eval. finite-volume integrals for mesh cells - splitting
+    each into a series of triangles and employing an order-2
+    quadrature rule.
+
+    """
+
+    class base:
+        pass
+
+    ncel = mesh.dimensions["nCells"].size
+    pcel = np.zeros((ncel, 2), dtype=np.float64)
+    pcel[:, 0] = mesh.variables["lonCell"][:]
+    pcel[:, 1] = mesh.variables["latCell"][:]
+
+    pcel = pcel * 180. / np.pi
+    pcel[pcel[:, 0] > 180., 0] -= 360.
+
+    fcel = sample_2d(xlon, ylat, vals, pcel[:, 0], pcel[:, 1])
+
+    pcel = pcel * np.pi / 180.
+
+    nvrt = mesh.dimensions["nVertices"].size
+    pvrt = np.zeros((nvrt, 2), dtype=np.float64)
+    pvrt[:, 0] = mesh.variables["lonVertex"][:]
+    pvrt[:, 1] = mesh.variables["latVertex"][:]
+
+    pvrt = pvrt * 180. / np.pi
+    pvrt[pvrt[:, 0] > 180., 0] -= 360.
+
+    fvrt = sample_2d(xlon, ylat, vals, pvrt[:, 0], pvrt[:, 1])
+
+    pvrt = pvrt * np.pi / 180.
+
+    cell = base()
+    cell.edge = mesh.variables["edgesOnCell"][:, :]
+    cell.topo = mesh.variables["nEdgesOnCell"][:]
+
+    edge = base()
+    edge.vert = mesh.variables["verticesOnEdge"][:]
+
+    abar = np.zeros((ncel, 1), dtype=np.float64)
+    fbar = np.zeros((ncel, 1), dtype=np.float64)
+
+    for epos in range(np.max(cell.topo)):
+
+        mask = cell.topo > epos
+
+        icel = np.argwhere(mask).ravel()
+
+        ifac = cell.edge[mask, epos] - 1
+
+        ivrt = edge.vert[ifac, 0] - 1
+        jvrt = edge.vert[ifac, 1] - 1
+
+        rsph = mesh.sphere_radius
+
+        atri = tria_area(
+            pcel[icel], pvrt[ivrt], pvrt[jvrt], rsph)
+
+        atri = np.reshape(atri, (atri.size, 1))
+
+        ftri = (fcel[icel] + fvrt[ivrt] + fvrt[jvrt])
+
+        ftri = np.reshape(ftri, (ftri.size, 1))
+
+        abar[icel] += atri
+        fbar[icel] += atri * ftri / 3.
+
+    fbar = np.asarray(fbar / abar, dtype=np.float32)
+
+    return fvrt, fcel, fbar
+
+
+def cell_prfl(mesh, smat,
+              nlev, zdem, zvrt, zcel, zbar):
+    """
+    Build elev. profiles for each cell in the mesh - sorting
+    the DEM pixel values assigned to each cell and assigning
+    to NLEV bins.
+
+    """
+
+    prfl = np.tile(zbar, (1, nlev))
+
+    nvrt = mesh.variables["nEdgesOnCell"][:]
+    ivrt = mesh.variables["verticesOnCell"][:, :] - 1
+
+    for cell in range(0, smat.shape[0]):
+
+        # -- extract set of DEM pixels per MPAS cell. Sparse SMAT
+        # -- contains cell-to-DEM mapping: ith row is the list of
+        # -- pixels overlapping with the ith cell
+
+        head = smat.indptr[cell + 0] + 0
+        tail = smat.indptr[cell + 1] + 0
+        idem = smat.indices[head:tail]
+
+        # -- build the cell elev. profiles: sort pixels by height
+        # -- and (linearly) interpolate to profile band endpoints
+
+        if (idem.size > nvrt[cell] + 1):
+
+            # -- list of contained DEM pixel values, for smooth cases
+            nvec = idem.shape[+0]
+            zvec = zdem[idem]
+
+            prfl[cell, :] = np.interp(
+                np.linspace(0., nvec - 1., nlev),
+                np.arange(0., nvec),
+                np.sort(zvec))
+
+        else:
+
+            # -- cell vert. + centre interpolations, for coarse cases
+            nvec = nvrt[cell] + 1
+            zvec = zvrt[ivrt[cell, :nvrt[cell]]]
+
+            prfl[cell, :] = np.interp(
+                np.linspace(0., nvec - 1., nlev),
+                np.arange(0., nvec),
+                np.sort(np.append(zvec, zcel[cell])))
+
+    return prfl
+
+
+def dem_remap(elev_file, mpas_file, elev_band=0):  # noqa: C901
+    """
+    Map elevation and ice+ocn-thickness data from a "zipped"
+    RTopo data-set onto the cells in an MPAS mesh.
+
+    Cell values are a blending of an approx. integral remap
+    and a local quadrature rule.
+
+    """
+
+    NLEV = elev_band + 1   # no. of evel. profile bands
+
+    print("Loading assests...")
+
+    elev = nc.Dataset(elev_file, "r+")
+    mesh = nc.Dataset(mpas_file, "r+")
+
+    xlon = np.asarray(elev["lon"][:], dtype=np.float64)
+    ylat = np.asarray(elev["lat"][:], dtype=np.float64)
+
+    # -- add dummy data to elev file if missing
+
+    if ("bed_elevation" not in elev.variables.keys()):
+        print("*bed_elevation variable not found")
+        elev.createVariable("bed_elevation",
+                            "i2", ("num_row", "num_col"))
+
+    if ("bed_slope" not in elev.variables.keys()):
+        print("*bed_slope variable not found")
+        elev.createVariable("bed_slope",
+                            "f4", ("num_row", "num_col"))
+
+    if ("bed_dz_dx" not in elev.variables.keys()):
+        print("*bed_dz_dx variable not found")
+        elev.createVariable("bed_dz_dx",
+                            "f4", ("num_row", "num_col"))
+
+    if ("bed_dz_dy" not in elev.variables.keys()):
+        print("*bed_dz_dy variable not found")
+        elev.createVariable("bed_dz_dy",
+                            "f4", ("num_row", "num_col"))
+
+    if ("ocn_thickness" not in elev.variables.keys()):
+        print("*ocn_thickness variable not found")
+        elev.createVariable("ocn_thickness",
+                            "i2", ("num_row", "num_col"))
+
+    if ("ice_thickness" not in elev.variables.keys()):
+        print("*ice_thickness variable not found")
+        elev.createVariable("ice_thickness",
+                            "i2", ("num_row", "num_col"))
+
+    # -- Compute an approximate remapping, associating pixels in
+    # -- the DEM with cells in the MPAS mesh. Since polygons are
+    # -- Voronoi, the point-in-cell query can be computed by
+    # -- finding nearest neighbours. This remapping is an approx.
+    # -- as no partial pixel-cell intersection is computed.
+
+    print("Building KDtree...")
+
+    ppos = np.zeros(
+        (mesh.dimensions["nCells"].size, 3), dtype=np.float64)
+    ppos[:, 0] = mesh["xCell"][:]
+    ppos[:, 1] = mesh["yCell"][:]
+    ppos[:, 2] = mesh["zCell"][:]
+
+    tree = spatial.cKDTree(ppos, leafsize=32)
+
+    print("Remap elevation...")
+
+    xmid = .5 * (xlon[:-1:] + xlon[1::])
+    ymid = .5 * (ylat[:-1:] + ylat[1::])
+
+    indx = np.asarray(np.round(
+        np.linspace(-1, ymid.size, 32)), dtype=int)
+
+    print("* process tiles:")
+
+    nset = []
+    for tile in range(indx.size - 1):
+
+        head = indx[tile + 0] + 1
+        tail = indx[tile + 1] + 1
+
+        qpos = map_to_r3(mesh, xmid, ymid, head, tail)
+
+        ttic = time.time()
+        try:  # ridiculous argument renaming...
+            __, cell = tree.query(qpos, n_jobs=-1)
+        except:  # noqa: E722
+            __, cell = tree.query(qpos, workers=-1)
+        ttoc = time.time()
+        print("* built node-to-cell map:",
+              tile, "of", indx.size - 1)
+
+        nset.append(
+            np.asarray(cell, dtype=np.uint32))
+
+    del tree
+    del ppos
+    del qpos
+
+    near = np.concatenate(nset)
+
+    del nset
+    del cell
+
+    # -- Build cell-to-pixel map as a sparse matrix, and average
+    # -- RTopo pixel values within each cell.
+
+    print("Form map matrix...")
+
+    ttic = time.time()
+
+    ncel = mesh.dimensions["nCells"].size
+
+    cols = np.arange(0, near.size, dtype=np.uint32)
+    vals = np.ones(near.size, dtype=np.int8)
+
+    smat = csr_matrix((vals, (near, cols)),
+                      shape=(ncel, near.size), dtype=np.int8)
+
+    del near
+    del cols
+    del vals
+
+    nmap = np.asarray(
+        smat.sum(axis=1), dtype=np.float32)
+
+    vals = np.asarray(
+        elev["bed_elevation"][:], dtype=np.float32)
+    vals = np.reshape(vals, (vals.size, 1))
+
+    emap = (smat * vals) / np.maximum(1., nmap)
+
+    vals = np.asarray(
+        elev["bed_slope"][:], dtype=np.float32)
+    vals = np.reshape(vals, (vals.size, 1))
+
+    smap = (smat * vals) / np.maximum(1., nmap)
+
+    vals = np.asarray(
+        elev["bed_dz_dx"][:], dtype=np.float32)
+    vals = np.reshape(vals, (vals.size, 1))
+
+    xmap = (smat * vals) / np.maximum(1., nmap)
+
+    vals = np.asarray(
+        elev["bed_dz_dy"][:], dtype=np.float32)
+    vals = np.reshape(vals, (vals.size, 1))
+
+    ymap = (smat * vals) / np.maximum(1., nmap)
+
+    vals = np.asarray(
+        elev["ocn_thickness"][:], dtype=np.float32)
+    vals = np.reshape(vals, (vals.size, 1))
+
+    omap = (smat * vals) / np.maximum(1., nmap)
+
+    frac = np.zeros(vals.shape, dtype=vals.dtype)
+    frac[vals > 0.0] = 1.
+    ofrc = (smat * frac) / np.maximum(1., nmap)
+
+    vals = np.asarray(
+        elev["ice_thickness"][:], dtype=np.float32)
+    vals = np.reshape(vals, (vals.size, 1))
+
+    imap = (smat * vals) / np.maximum(1., nmap)
+
+    frac = np.zeros(vals.shape, dtype=vals.dtype)
+    frac[vals > 0.0] = 1.
+    ifrc = (smat * frac) / np.maximum(1., nmap)
+
+    ttoc = time.time()
+
+    del vals
+    del frac
+
+    print("* built remapping matrix:",
+          np.round(ttoc - ttic, decimals=1), "sec")
+
+    # -- If the resolution of the mesh is greater, or comparable
+    # -- to that of the DEM, the approx. remapping (above) will
+    # -- result in a low order interpolation.
+    # -- Thus, blend with a local order-2 quadrature formulation
+
+    print("Eval. elevation...")
+
+    vals = np.asarray(
+        elev["bed_elevation"][:], dtype=np.float32)
+
+    ttic = time.time()
+    evrt, ecel, eint = cell_quad(mesh, xlon, ylat, vals)
+    ttoc = time.time()
+    print("* compute cell integrals:",
+          np.round(ttoc - ttic, decimals=1), "sec")
+
+    vals = np.asarray(
+        elev["bed_slope"][:], dtype=np.float32)
+
+    ttic = time.time()
+    svrt, scel, sint = cell_quad(mesh, xlon, ylat, vals)
+    ttoc = time.time()
+    print("* compute cell integrals:",
+          np.round(ttoc - ttic, decimals=1), "sec")
+
+    vals = np.asarray(
+        elev["bed_dz_dx"][:], dtype=np.float32)
+
+    ttic = time.time()
+    xvrt, xcel, xint = cell_quad(mesh, xlon, ylat, vals)
+    ttoc = time.time()
+    print("* compute cell integrals:",
+          np.round(ttoc - ttic, decimals=1), "sec")
+
+    vals = np.asarray(
+        elev["bed_dz_dy"][:], dtype=np.float32)
+
+    ttic = time.time()
+    yvrt, ycel, yint = cell_quad(mesh, xlon, ylat, vals)
+    ttoc = time.time()
+    print("* compute cell integrals:",
+          np.round(ttoc - ttic, decimals=1), "sec")
+
+    vals = np.asarray(
+        elev["ocn_thickness"][:], dtype=np.float32)
+
+    ttic = time.time()
+    ovrt, ocel, oint = cell_quad(mesh, xlon, ylat, vals)
+    ttoc = time.time()
+    print("* compute cell integrals:",
+          np.round(ttoc - ttic, decimals=1), "sec")
+
+    vals = np.asarray(
+        elev["ice_thickness"][:], dtype=np.float32)
+
+    ttic = time.time()
+    ivrt, icel, iint = cell_quad(mesh, xlon, ylat, vals)
+    ttoc = time.time()
+    print("* compute cell integrals:",
+          np.round(ttoc - ttic, decimals=1), "sec")
+
+    print("Save to dataset...")
+
+    ebar = (np.multiply(nmap, emap) + 4 * eint) / (4 + nmap)
+    sbar = (np.multiply(nmap, smap) + 4 * sint) / (4 + nmap)
+    xbar = (np.multiply(nmap, xmap) + 4 * xint) / (4 + nmap)
+    ybar = (np.multiply(nmap, ymap) + 4 * yint) / (4 + nmap)
+    obar = (np.multiply(nmap, omap) + 4 * oint) / (4 + nmap)
+    ibar = (np.multiply(nmap, imap) + 4 * iint) / (4 + nmap)
+
+    obar = np.maximum(0., np.minimum(obar, -ebar))
+
+    tfrc = np.zeros(ofrc.shape, dtype=np.float32)
+    tfrc[ebar <= +0.] = -ebar[ebar <= +0.]
+    tfrc[ebar <= -1.] = 1.
+    ofrc = (np.multiply(nmap, ofrc) + 4 * tfrc) / (4 + nmap)
+
+    tfrc = np.zeros(ofrc.shape, dtype=np.float32)
+    tfrc[ibar >= +0.] = +ibar[ibar >= +0.]
+    tfrc[ibar >= +1.] = 1.
+    ifrc = (np.multiply(nmap, ifrc) + 4 * tfrc) / (4 + nmap)
+
+    if ("bed_elevation" not in mesh.variables.keys()):
+        mesh.createVariable("bed_elevation", "f4", ("nCells"))
+
+    if ("bed_slope" not in mesh.variables.keys()):
+        mesh.createVariable("bed_slope", "f4", ("nCells"))
+
+    if ("bed_dz_dx" not in mesh.variables.keys()):
+        mesh.createVariable("bed_dz_dx", "f4", ("nCells"))
+
+    if ("bed_dz_dy" not in mesh.variables.keys()):
+        mesh.createVariable("bed_dz_dy", "f4", ("nCells"))
+
+    if ("ocn_thickness" not in mesh.variables.keys()):
+        mesh.createVariable("ocn_thickness", "f4", ("nCells"))
+
+    if ("ice_thickness" not in mesh.variables.keys()):
+        mesh.createVariable("ice_thickness", "f4", ("nCells"))
+
+    mesh["bed_elevation"].units = "m"
+    mesh["bed_elevation"][:] = ebar
+    mesh["bed_elevation"].long_name = names.bed_elevation
+
+    mesh["ocn_thickness"].units = "m"
+    mesh["ocn_thickness"][:] = obar
+    mesh["ocn_thickness"].long_name = names.ocn_thickness
+
+    mesh["ice_thickness"].units = "m"
+    mesh["ice_thickness"][:] = ibar
+    mesh["ice_thickness"].long_name = names.ice_thickness
+
+    mesh["bed_dz_dx"][:] = xbar
+    mesh["bed_dz_dx"].long_name = names.bed_dz_dx
+    mesh["bed_dz_dy"][:] = ybar
+    mesh["bed_dz_dy"].long_name = names.bed_dz_dy
+
+    mesh["bed_slope"].units = "deg"
+    mesh["bed_slope"][:] = \
+        np.arctan(sbar) * 180. / np.pi  # degrees, for ELM
+    mesh["bed_slope"].long_name = names.bed_slope_deg
+
+    if ("ocn_cover" not in mesh.variables.keys()):
+        mesh.createVariable("ocn_cover", "f4", ("nCells"))
+
+    if ("ice_cover" not in mesh.variables.keys()):
+        mesh.createVariable("ice_cover", "f4", ("nCells"))
+
+    mesh["ocn_cover"][:] = ofrc
+    mesh["ocn_cover"].long_name = names.ocn_cover
+    mesh["ice_cover"][:] = ifrc
+    mesh["ice_cover"].long_name = names.ice_cover
+
+    del emap
+    del eint
+    del smap
+    del sint
+    del xmap
+    del xint
+    del ymap
+    del yint
+    del omap
+    del oint
+    del imap
+    del iint
+    del ofrc
+    del ifrc
+
+    # -- Also compute profiles (ie. histograms) of elev. outputs
+    # -- per cell, dividing distributions of DEM pixel
+    # -- values into NLEV-1 bands. Write band endpoints to file.
+
+    if (elev_band <= 0):
+        return
+
+    print("Eval. histogram...")
+
+    vals = np.asarray(
+        elev["bed_elevation"][:], dtype=np.float32)
+    vals = np.reshape(vals, (vals.size))
+
+    ttic = time.time()
+    eprf = cell_prfl(
+        mesh, smat, NLEV, vals, evrt, ecel, ebar)
+    ttoc = time.time()
+    print("* compute elev. profiles:",
+          np.round(ttoc - ttic, decimals=1), "sec")
+
+    vals = np.asarray(
+        elev["bed_slope"][:], dtype=np.float32)
+    vals = np.reshape(vals, (vals.size))
+
+    ttic = time.time()
+    sprf = cell_prfl(
+        mesh, smat, NLEV, vals, svrt, scel, sbar)
+    ttoc = time.time()
+    print("* compute elev. profiles:",
+          np.round(ttoc - ttic, decimals=1), "sec")
+
+    vals = np.asarray(
+        elev["ocn_thickness"][:], dtype=np.float32)
+    vals = np.reshape(vals, (vals.size))
+
+    ttic = time.time()
+    oprf = cell_prfl(
+        mesh, smat, NLEV, vals, ovrt, ocel, obar)
+    ttoc = time.time()
+    print("* compute elev. profiles:",
+          np.round(ttoc - ttic, decimals=1), "sec")
+
+    vals = np.asarray(
+        elev["ice_thickness"][:], dtype=np.float32)
+    vals = np.reshape(vals, (vals.size))
+
+    ttic = time.time()
+    iprf = cell_prfl(
+        mesh, smat, NLEV, vals, ivrt, icel, ibar)
+    ttoc = time.time()
+    print("* compute elev. profiles:",
+          np.round(ttoc - ttic, decimals=1), "sec")
+
+    print("Save to dataset...")
+
+    if ("nProfiles" not in mesh.dimensions.keys()):
+        mesh.createDimension("nProfiles", NLEV)
+
+    if ("bed_elevation_profile" not in mesh.variables.keys()):
+        mesh.createVariable("bed_elevation_profile",
+                            "f4", ("nCells", "nProfiles"))
+
+    if ("bed_slope_profile" not in mesh.variables.keys()):
+        mesh.createVariable("bed_slope_profile",
+                            "f4", ("nCells", "nProfiles"))
+
+    if ("ocn_thickness_profile" not in mesh.variables.keys()):
+        mesh.createVariable("ocn_thickness_profile",
+                            "f4", ("nCells", "nProfiles"))
+
+    if ("ice_thickness_profile" not in mesh.variables.keys()):
+        mesh.createVariable("ice_thickness_profile",
+                            "f4", ("nCells", "nProfiles"))
+
+    mesh["bed_elevation_profile"].units = "m"
+    mesh["bed_elevation_profile"][:, :] = eprf
+    mesh["bed_elevation_profile"].long_name = names.bed_elevation_profile
+
+    mesh["bed_slope_profile"][:, :] = sprf
+    mesh["bed_slope_profile"].long_name = names.bed_slope_profile
+
+    mesh["ocn_thickness_profile"].units = "m"
+    mesh["ocn_thickness_profile"][:, :] = oprf
+    mesh["ocn_thickness_profile"].long_name = names.ocn_thickness_profile
+
+    mesh["ice_thickness_profile"].units = "m"
+    mesh["ice_thickness_profile"][:, :] = iprf
+    mesh["ice_thickness_profile"].long_name = names.ice_thickness_profile
+
+    elev.close()
+    mesh.close()
+
+
+if (__name__ == "__main__"):
+    parser = argparse.ArgumentParser(
+        description=__doc__,
+        formatter_class=argparse.RawTextHelpFormatter)
+
+    parser.add_argument(
+        "--mpas-file", dest="mpas_file", type=str,
+        required=True, help="Name of user MPAS mesh.")
+
+    parser.add_argument(
+        "--elev-file", dest="elev_file", type=str,
+        required=True, help="Name of DEM pixel file.")
+
+    parser.add_argument(
+        "--elev-band", dest="elev_band", type=int,
+        default=0,
+        required=False, help="Elev. profile band(s).")
+
+    dem_remap(parser.parse_args())
+
+    args = parser.parse_args()
+    elev_file = args.elev_file
+    mpas_file = args.mpas_file
+    elev_band = args.elev_band
+    dem_remap(elev_file, mpas_file, elev_band)
diff --git a/compass/ocean/tests/tides/dem/dem_trnsf.py b/compass/ocean/tests/tides/dem/dem_trnsf.py
new file mode 100644
index 0000000000..418be28391
--- /dev/null
+++ b/compass/ocean/tests/tides/dem/dem_trnsf.py
@@ -0,0 +1,248 @@
+
+import argparse
+
+import netCDF4 as nc
+import numpy as np
+
+from compass.ocean.tests.tides.dem.dem_names import names
+
+# Authors: Darren Engwirda
+
+
+def dem_trnsf(base_mesh, part_mesh):  # noqa: C901
+    """
+    Transfer elevation & ice+ocn-thickness data from a full
+    sphere MPAS mesh onto a partial sub-mesh resulting from
+    a culling operation, or equiv.
+
+    The sub-mesh is expected to be a piece of the base mesh
+    such that (some) cells match 1:1.
+
+    """
+
+    base = nc.Dataset(base_mesh, "r")
+    part = nc.Dataset(part_mesh, "r+")
+
+    # -- make a vector of cell-centre positions to match against
+
+    xpos = np.vstack((
+        np.vstack((
+            base["xCell"][:],
+            base["yCell"][:],
+            base["zCell"][:])).T,
+        np.vstack((
+            part["xCell"][:],
+            part["yCell"][:],
+            part["zCell"][:])).T
+    ))
+
+    xidx = np.hstack((
+        np.arange(
+            0, base.dimensions["nCells"].size),
+        np.arange(
+            0, part.dimensions["nCells"].size)
+    ))
+
+    # -- culling shouldn't introduce fp round-off - but truncate
+    # -- anyway...
+
+    xpos = np.round(xpos, decimals=9)
+
+    # -- use stable sorting to bring matching cell xyz (and idx)
+    # -- into "ascending" order
+
+    imap = np.argsort(xpos[:, 2], kind="stable")
+    xpos = xpos[imap, :]
+    xidx = xidx[imap]
+    imap = np.argsort(xpos[:, 1], kind="stable")
+    xpos = xpos[imap, :]
+    xidx = xidx[imap]
+    imap = np.argsort(xpos[:, 0], kind="stable")
+    xpos = xpos[imap, :]
+    xidx = xidx[imap]
+
+    diff = xpos[1:, :] - xpos[:-1, :]
+
+    same = np.argwhere(np.logical_and.reduce((
+        diff[:, 0] == 0,
+        diff[:, 1] == 0,
+        diff[:, 2] == 0))).ravel()
+
+    # -- cell inew in part matches iold in base - re-index elev.
+    # -- data-sets
+
+    inew = xidx[same + 1]
+    iold = xidx[same + 0]
+
+    if ("bed_elevation" in base.variables.keys()):
+        if ("bed_elevation" not in part.variables.keys()):
+            part.createVariable("bed_elevation", "f4", ("nCells"))
+
+        btmp = np.asarray(
+            base["bed_elevation"][:], dtype=np.float32)
+
+        ncel = part.dimensions["nCells"].size
+
+        ptmp = np.zeros(ncel, dtype=np.float32)
+        ptmp[inew] = btmp[iold]
+
+        part["bed_elevation"].units = "m"
+        part["bed_elevation"].long_name = names.bed_elevation
+        part["bed_elevation"][:] = ptmp
+
+    if ("bed_slope" in base.variables.keys()):
+        if ("bed_slope" not in part.variables.keys()):
+            part.createVariable("bed_slope", "f4", ("nCells"))
+
+        btmp = np.asarray(
+            base["bed_slope"][:], dtype=np.float32)
+
+        ncel = part.dimensions["nCells"].size
+
+        ptmp = np.zeros(ncel, dtype=np.float32)
+        ptmp[inew] = btmp[iold]
+
+        part["bed_slope"].units = "deg"
+        part["bed_slope"].long_name = names.bed_slope_deg
+        part["bed_slope"][:] = ptmp
+
+    if ("bed_dz_dx" in base.variables.keys()):
+        if ("bed_dz_dx" not in part.variables.keys()):
+            part.createVariable("bed_dz_dx", "f4", ("nCells"))
+
+        btmp = np.asarray(
+            base["bed_dz_dx"][:], dtype=np.float32)
+
+        ncel = part.dimensions["nCells"].size
+
+        ptmp = np.zeros(ncel, dtype=np.float32)
+        ptmp[inew] = btmp[iold]
+
+        part["bed_dz_dx"].long_name = names.bed_dz_dx
+        part["bed_dz_dx"][:] = ptmp
+
+    if ("bed_dz_dy" in base.variables.keys()):
+        if ("bed_dz_dy" not in part.variables.keys()):
+            part.createVariable("bed_dz_dy", "f4", ("nCells"))
+
+        btmp = np.asarray(
+            base["bed_dz_dy"][:], dtype=np.float32)
+
+        ncel = part.dimensions["nCells"].size
+
+        ptmp = np.zeros(ncel, dtype=np.float32)
+        ptmp[inew] = btmp[iold]
+
+        part["bed_dz_dy"].long_name = names.bed_dz_dy
+        part["bed_dz_dy"][:] = ptmp
+
+    if ("ocn_thickness" not in part.variables.keys()):
+        if ("ocn_thickness" not in part.variables.keys()):
+            part.createVariable("ocn_thickness", "f4", ("nCells"))
+
+        btmp = np.asarray(
+            base["ocn_thickness"][:], dtype=np.float32)
+
+        ncel = part.dimensions["nCells"].size
+
+        ptmp = np.zeros(ncel, dtype=np.float32)
+        ptmp[inew] = btmp[iold]
+
+        part["ocn_thickness"].units = "m"
+        part["ocn_thickness"].long_name = names.ocn_thickness
+        part["ocn_thickness"][:] = ptmp
+
+    if ("ice_thickness" not in part.variables.keys()):
+        if ("ice_thickness" not in part.variables.keys()):
+            part.createVariable("ice_thickness", "f4", ("nCells"))
+
+        btmp = np.asarray(
+            base["ice_thickness"][:], dtype=np.float32)
+
+        ncel = part.dimensions["nCells"].size
+
+        ptmp = np.zeros(ncel, dtype=np.float32)
+        ptmp[inew] = btmp[iold]
+
+        part["ice_thickness"].units = "m"
+        part["ice_thickness"].long_name = names.ice_thickness
+        part["ice_thickness"][:] = ptmp
+
+    if ("bed_elevation_profile" in base.variables.keys()):
+        if ("bed_elevation_profile" not in part.variables.keys()):
+            part.createVariable("bed_elevation_profile",
+                                "f4", ("nCells", "nProfiles"))
+
+        btmp = np.asarray(
+            base["bed_elevation_profile"][:, :], dtype=np.float32)
+
+        ncel = part.dimensions["nCells"].size
+        nprf = part.dimensions["nProfiles"].size
+
+        ptmp = np.zeros((ncel, nprf), dtype=np.float32)
+        ptmp[inew, :] = btmp[iold, :]
+
+        part["bed_elevation_profile"].units = "m"
+        part["bed_elevation_profile"].long_name = names.bed_elevation_profile
+        part["bed_elevation_profile"][:, :] = ptmp
+
+    if ("ocn_thickness_profile" in base.variables.keys()):
+        if ("ocn_thickness_profile" not in part.variables.keys()):
+            part.createVariable("ocn_thickness_profile",
+                                "f4", ("nCells", "nProfiles"))
+
+        btmp = np.asarray(
+            base["ocn_thickness_profile"][:, :], dtype=np.float32)
+
+        ncel = part.dimensions["nCells"].size
+        nprf = part.dimensions["nProfiles"].size
+
+        ptmp = np.zeros((ncel, nprf), dtype=np.float32)
+        ptmp[inew, :] = btmp[iold, :]
+
+        part["ocn_thickness_profile"].units = "m"
+        part["ocn_thickness_profile"].long_name = names.ocn_thickness_profile
+        part["ocn_thickness_profile"][:, :] = ptmp
+
+    if ("ice_thickness_profile" in base.variables.keys()):
+        if ("ice_thickness_profile" not in part.variables.keys()):
+            part.createVariable("ice_thickness_profile",
+                                "f4", ("nCells", "nProfiles"))
+
+        btmp = np.asarray(
+            base["ice_thickness_profile"][:, :], dtype=np.float32)
+
+        ncel = part.dimensions["nCells"].size
+        nprf = part.dimensions["nProfiles"].size
+
+        ptmp = np.zeros((ncel, nprf), dtype=np.float32)
+        ptmp[inew, :] = btmp[iold, :]
+
+        part["ice_thickness_profile"].units = "m"
+        part["ice_thickness_profile"].long_name = names.ice_thickness_profile
+        part["ice_thickness_profile"][:, :] = ptmp
+
+    base.close()
+    part.close()
+
+    return
+
+
+if (__name__ == "__main__"):
+    parser = argparse.ArgumentParser(
+        description=__doc__,
+        formatter_class=argparse.RawTextHelpFormatter)
+
+    parser.add_argument(
+        "--base-mesh", dest="base_mesh", type=str,
+        required=True, help="Name of (full) MPAS mesh.")
+
+    parser.add_argument(
+        "--part-mesh", dest="part_mesh", type=str,
+        required=True, help="Name of culled MPAS mesh.")
+
+    args = parser.parse_args()
+    base_mesh = args.base_mesh
+    part_mesh = args.part_mesh
+
+    dem_trnsf(base_mesh, part_mesh)
diff --git a/compass/ocean/tests/tides/forward/forward.py b/compass/ocean/tests/tides/forward/forward.py
index 9522e54715..c5bb4c5798 100644
--- a/compass/ocean/tests/tides/forward/forward.py
+++ b/compass/ocean/tests/tides/forward/forward.py
@@ -49,7 +49,7 @@ def __init__(self, test_case, mesh, init, name='forward', subdir=None):
         self.add_namelist_file(mesh_package, 'namelist.ocean')
 
         initial_state_path = f'{init.path}/initial_state'
-        interpolate_path = f'{init.path}/interpolate'
+        wave_drag_path = f'{init.path}/wave_drag'
         initial_state_target = f'{initial_state_path}/initial_state.nc'
 
         self.add_input_file(
@@ -60,7 +60,7 @@ def __init__(self, test_case, mesh, init, name='forward', subdir=None):
             work_dir_target=f'{initial_state_path}/init_mode_forcing_data.nc')
         self.add_input_file(
             filename='topographic_wave_drag.nc',
-            work_dir_target=f'{interpolate_path}/topographic_wave_drag.nc')
+            work_dir_target=f'{wave_drag_path}/topographic_wave_drag.nc')
         self.add_input_file(
             filename='graph.info',
             work_dir_target=f'{initial_state_path}/graph.info')
diff --git a/compass/ocean/tests/tides/forward/namelist.ocean b/compass/ocean/tests/tides/forward/namelist.ocean
index 9c5f58645c..9d3004f66f 100644
--- a/compass/ocean/tests/tides/forward/namelist.ocean
+++ b/compass/ocean/tests/tides/forward/namelist.ocean
@@ -4,10 +4,14 @@ config_calendar_type = 'noleap'
 config_pio_num_iotasks = 0
 config_pio_stride = 1
 config_time_integrator = 'RK4'
+config_use_mom_del2 = .true.
+config_mom_del2 = 3.333
+config_use_mom_del4 = .true.
+config_mom_del4 = 3.333e10
 config_hmix_scaleWithMesh = .true.
 config_hmix_use_ref_cell_width = .true.
 config_use_cvmix = .false.
-config_use_bulk_wind_stress = .true.
+config_use_bulk_wind_stress = .false.
 config_use_self_attraction_loading = .true.
 config_self_attraction_loading_compute_interval = '0000_00:30:00'
 config_use_parallel_self_attraction_loading = .true.
@@ -21,8 +25,10 @@ config_use_tidal_potential_forcing_P1 = .false.
 config_tidal_potential_ramp = 15.0
 config_implicit_bottom_drag_type = 'loglaw'
 config_use_topographic_wave_drag = .true.
+config_Rayleigh_damping_coeff = 1.0e-5
 config_use_wetting_drying = .true.
-config_prevent_drying = .true
+config_prevent_drying = .true.
+config_land_ice_flux_mode = 'pressure_only'
 config_drying_min_cell_height = 0.1
 config_zero_drying_velocity = .true.
 config_pressure_gradient_type = 'ssh_gradient'
@@ -36,7 +42,7 @@ config_AM_globalStats_compute_on_startup = .true.
 config_AM_globalStats_write_on_startup = .true.
 config_AM_globalStats_text_file = .true.
 config_AM_harmonicAnalysis_enable = .true.
-config_AM_harmonicAnalysis_compute_interval = 00:30:00
+config_AM_harmonicAnalysis_compute_interval = '00:30:00'
 config_AM_harmonicAnalysis_start_delay = 30
 config_AM_harmonicAnalysis_duration = 90
 config_AM_harmonicAnalysis_use_K2 = .false.
diff --git a/compass/ocean/tests/tides/forward/streams.ocean b/compass/ocean/tests/tides/forward/streams.ocean
index 6a8bd3a27c..7201e7154e 100755
--- a/compass/ocean/tests/tides/forward/streams.ocean
+++ b/compass/ocean/tests/tides/forward/streams.ocean
@@ -29,7 +29,6 @@
 
     <stream name="mesh"/>
     <var name="ssh"/>
-    <var name="ssh_sal"/>
     <var name="layerThickness"/>
     <var name="xtime"/>
     <var name="kineticEnergyCell"/>
@@ -61,9 +60,14 @@
 <stream name="topographic_wave_drag"
         type="input"
         filename_template="topographic_wave_drag.nc"
-        input_interval="initial_only" >
-
-        <var name="topographic_wave_drag"/>
+        input_interval="initial_only">
+
+        <var name="topo_buoyancy_N1V"/>
+        <var name="topo_buoyancy_N1B"/>
+        <var name="bathy_stddev"/>
+        <var name="bed_slope_edges"/>
+        <var name="latGradEdge"/>
+        <var name="lonGradEdge"/>
 </stream>
 
 <stream name="globalStatsOutput"
diff --git a/compass/ocean/tests/tides/init/__init__.py b/compass/ocean/tests/tides/init/__init__.py
index be94f9b622..cbaa2d2b73 100644
--- a/compass/ocean/tests/tides/init/__init__.py
+++ b/compass/ocean/tests/tides/init/__init__.py
@@ -1,12 +1,15 @@
 import os
 
-from compass.testcase import TestCase
-from compass.ocean.tests.tides.init.interpolate_wave_drag \
-    import InterpolateWaveDrag
-from compass.ocean.tests.tides.init.remap_bathymetry \
-    import RemapBathymetry
-from compass.ocean.tests.tides.init.initial_state import InitialState
 from compass.ocean.tests.tides.configure import configure_tides
+from compass.ocean.tests.tides.init.calculate_wave_drag import (
+    CalculateWaveDrag,
+)
+from compass.ocean.tests.tides.init.initial_state import InitialState
+from compass.ocean.tests.tides.init.interpolate_wave_drag import (
+    InterpolateWaveDrag,
+)
+from compass.ocean.tests.tides.init.remap_bathymetry import RemapBathymetry
+from compass.testcase import TestCase
 
 
 class Init(TestCase):
@@ -37,7 +40,8 @@ def __init__(self, test_group, mesh):
 
         self.mesh = mesh
 
-        self.add_step(InterpolateWaveDrag(test_case=self, mesh=mesh))
+        # self.add_step(InterpolateWaveDrag(test_case=self, mesh=mesh))
+        self.add_step(CalculateWaveDrag(test_case=self, mesh=mesh))
         self.add_step(RemapBathymetry(test_case=self, mesh=mesh))
         self.add_step(InitialState(test_case=self, mesh=mesh))
 
diff --git a/compass/ocean/tests/tides/init/calculate_wave_drag.py b/compass/ocean/tests/tides/init/calculate_wave_drag.py
new file mode 100644
index 0000000000..f2dadb855b
--- /dev/null
+++ b/compass/ocean/tests/tides/init/calculate_wave_drag.py
@@ -0,0 +1,335 @@
+# Creates topographic_wave_drag.nc file containing all data
+# needed for each type of wave drag scheme.
+
+import netCDF4
+import numpy as np
+import xarray as xr
+from scipy import interpolate
+from scipy.spatial import cKDTree
+
+from compass.step import Step
+
+
+class CalculateWaveDrag(Step):
+    """
+    A step for calculating wave drag information
+
+    Attributes
+    ----------
+    bathy_file : str
+       File name for the blended RTopo/GEBCO pixel file
+
+    mesh_file : str
+       File name for the culled mesh from the mesh test case
+
+    bouy_file : str
+       File name for the WOA bouyancy data
+
+    output_file : str
+       File name for the output file containing wave drag information
+    """
+
+    def __init__(self, test_case, mesh):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        test_case : compass.ocean.tests.tides.init.Init
+            The test case this step belongs to
+
+         mesh : compass.ocean.tests.tides.mesh.Mesh
+            The test case the produces the mesh for this case
+
+        """
+
+        super().__init__(test_case=test_case, name='wave_drag',
+                         ntasks=1, min_tasks=1, openmp_threads=1)
+
+        self.bathy_file = 'bathy.nc'
+        self.mesh_file = 'mesh.nc'
+        self.buoy_file = 'buoy.nc'
+
+        self.add_input_file(
+            filename=self.buoy_file,
+            target='WOA2013_buoyancies.nc',
+            database='initial_condition_database')
+
+        mesh_path = mesh.steps['cull_mesh'].path
+        self.add_input_file(
+            filename=self.mesh_file,
+            work_dir_target=f'{mesh_path}/culled_mesh.nc')
+
+        bathy_path = mesh.steps['pixel'].path
+        self.add_input_file(
+            filename=self.bathy_file,
+            work_dir_target=f'{bathy_path}/'
+                            'RTopo_2_0_4_GEBCO_v2023_30sec_pixel.nc')
+
+        self.output_file = 'topographic_wave_drag.nc'
+        self.add_output_file(filename=self.output_file)
+
+    def interpolate_data_to_grid(self, grid_file, data_file, var):
+        """
+        Interpolate a variable from a data file to an the MPAS grid
+        """
+
+        # Open files
+        data_nc = netCDF4.Dataset(data_file, 'r')
+        grid_nc = netCDF4.Dataset(grid_file, 'r')
+
+        # Get grid from data file
+        lon_data = data_nc.variables['longitude'][:]
+        lat_data = data_nc.variables['latitude'][:]
+        nsnaps = 1
+
+        # Get grid from grid file
+        lon_grid = np.mod(grid_nc.variables['lonEdge'][:] + np.pi,
+                          2.0 * np.pi) - np.pi
+        lon_grid = lon_grid * 180.0 / np.pi
+        lat_grid = grid_nc.variables['latEdge'][:] * 180.0 / np.pi
+
+        grid_points = np.column_stack((lon_grid, lat_grid))
+        nEdges = lon_grid.size
+        interp_data = np.zeros((nsnaps, nEdges))
+
+        # Interpolate timesnaps
+        print(f'Interpolating {var}')
+
+        # Get data to interpolate
+        data = data_nc.variables[var][:]
+        # bottom_index = data_nc.variables['bottom_index'][:]
+        data[data == np.nan] = 0
+        data[data > 0.015] = 0.015
+        data[data < 1.405189e-4] = 1.405189e-4
+        # Interpolate data onto new grid
+        interpolator = interpolate.RegularGridInterpolator(
+            (lon_data, lat_data), data,
+            bounds_error=False, fill_value=None, method='nearest')
+        interp_data[0, :] = interpolator(grid_points)
+
+        return interp_data
+
+    def run(self):
+        """
+        Run the step to calculate wave drag information
+        """
+
+        # Open Datasets
+        print("Opening Datasets")
+        bd = xr.open_dataset(self.bathy_file)
+        elev = np.asarray(bd.bed_elevation)  # (lat, lon)
+        ylat = np.asarray(bd.lat.values[:], dtype=np.float64)
+        xlon = np.asarray(bd.lon.values[:], dtype=np.float64)
+        bed_slope = np.asarray(bd.bed_slope.values)
+        xgrad = np.asarray(bd.bed_dz_dx.values)
+        ygrad = np.asarray(bd.bed_dz_dy.values)
+
+        md = xr.open_dataset(self.mesh_file)
+        nEdges = np.size(md.nEdges.values)
+        xEdges = np.asarray(md.xEdge.values)
+        yEdges = np.asarray(md.yEdge.values)
+        zEdges = np.asarray(md.zEdge.values)
+        xyzVertex = np.vstack((xEdges, yEdges, zEdges)).T
+
+        print("Getting nearest neighbors")
+        tree = cKDTree(xyzVertex)
+        xmid = .5 * (xlon[:-1:] + xlon[1::])
+        ymid = .5 * (ylat[:-1:] + ylat[1::])
+        indx = np.asarray(np.round(
+            np.linspace(-1, ymid.size, 17)), dtype=int)
+        nset = []
+
+        for tile in range(indx.size - 1):
+
+            head = indx[tile + 0] + 1
+            tail = indx[tile + 1] + 1
+
+            qpos = self.map_to_r3(xmid, ymid, head, tail)
+
+            __, nloc = tree.query(qpos, workers=-1)
+
+            nset.append(nloc)
+
+        near = np.concatenate(nset)
+        near = near.reshape((np.size(ymid), np.size(xmid)))
+        del tree
+        del qpos
+        del nset
+        del nloc
+
+        print("Calculating Standard Deviation and Interpolating Bed Slope")
+        stddev = np.zeros((nEdges))
+        bed_slope_edges = np.zeros((nEdges))
+        xGradEdges = np.zeros((nEdges))
+        yGradEdges = np.zeros((nEdges))
+        N1V_interp = np.zeros((nEdges))
+        N1B_interp = np.zeros((nEdges))
+        xlon_nn = [[] for _ in range(nEdges)]
+        ylat_nn = [[] for _ in range(nEdges)]
+        print("Making lists")
+        for i in range(np.size(ymid)):
+            for j in range(np.size(xmid)):
+                ylat_nn, xlon_nn = self.make_nn_lists(ylat_nn, xlon_nn,
+                                                      near, i, j)
+
+        print("Calculating stats")
+        for edge in range(nEdges):
+            bed_slope_edge, sd, xGradEdge, yGradEdge = self.calc_stats(
+                xlon_nn, ylat_nn, near, bed_slope,
+                edge, xmid, ymid, elev, xgrad, ygrad)
+
+            bed_slope_edges[edge] = bed_slope_edge
+            stddev[edge] = sd
+            xGradEdges[edge] = xGradEdge
+            yGradEdges[edge] = yGradEdge
+
+        # Fix some nn issues
+        print("Fixing nans in data")
+        while (np.sum(np.isnan(stddev)) > 0):
+            print("stddev", np.sum(np.isnan(stddev)))
+            stddev = self.fix_nans(stddev, md)
+        while (np.sum(np.isnan(bed_slope_edges)) > 0):
+            print("bed_slope_edges", np.sum(np.isnan(stddev)))
+            bed_slope_edges = self.fix_nans(bed_slope_edges, md)
+        while (np.sum(np.isnan(xGradEdges)) > 0):
+            print("xGradEdges", np.sum(np.isnan(stddev)))
+            xGradEdges = self.fix_nans(xGradEdges, md)
+        while (np.sum(np.isnan(yGradEdges)) > 0):
+            print("yGradEdges", np.sum(np.isnan(stddev)))
+            yGradEdge = self.fix_nans(yGradEdges, md)
+
+        # Remove any 0s from the bottom Depth
+        # bottomDepthEdges[bottomDepthEdges==0.0] = 0.01
+
+        N1V_interp = self.interpolate_data_to_grid(self.mesh_file,
+                                                   self.buoy_file, 'N1V')
+        N1B_interp = self.interpolate_data_to_grid(self.mesh_file,
+                                                   self.buoy_file, 'N1B')
+        twd = xr.Dataset()
+        twd["topo_buoyancy_N1V"] = xr.DataArray(N1V_interp[0, :],
+                                                dims={"nEdges"})
+        twd["topo_buoyancy_N1B"] = xr.DataArray(N1B_interp[0, :],
+                                                dims={"nEdges"})
+        twd["bathy_stddev"] = xr.DataArray(stddev, dims={"nEdges"})
+        twd["bed_slope_edges"] = xr.DataArray(bed_slope_edges, dims={"nEdges"})
+        twd["lonGradEdge"] = xr.DataArray(xGradEdges, dims={"nEdges"})
+        twd["latGradEdge"] = xr.DataArray(yGradEdges, dims={"nEdges"})
+        twd.to_netcdf(self.output_file)
+
+        del xlon_nn
+        del ylat_nn
+
+    def fix_nans(self, data, mesh_nc):
+        """
+        Replace NaN values with average of surrounding (non-NaN) values
+        """
+        if (np.sum(np.isnan(data))):
+            edgesOnEdge = mesh_nc.edgesOnEdge.data
+
+            nanEdges = np.isnan(data)
+            for edge in range(mesh_nc.nEdges.data[-1]):
+                if nanEdges[edge]:
+                    surrounding = edgesOnEdge[edge]
+                    surrounding = surrounding[surrounding > 0]
+                    surrounding = surrounding[~np.isnan(data[surrounding])]
+
+                    data[edge] = np.mean(data[surrounding])
+
+        return data
+
+    def make_nn_lists(self, ylat_nn, xlon_nn, near, i, j):
+        """
+        Create the nearest neighbor lists
+        """
+        ylat_nn[near[i, j]].append(i)
+        xlon_nn[near[i, j]].append(j)
+        return ylat_nn, xlon_nn
+
+    def calc_stats(self, xlon_nn, ylat_nn, near, bed_slope,
+                   edge, xlon, ylat, elev, xgrad, ygrad):
+        """
+        Calculate the mean bathymetry slopes/gradients and stadard deviation
+        """
+
+        # Prepare for calculating area integral
+        # dlat = ylat[1] - ylat[0]
+        # dlon = xlon[1] - xlon[0]
+        # R = 6378206.4
+
+        # Standard deviation
+        xlon_idx = xlon_nn[edge]
+        ylat_idx = ylat_nn[edge]
+        Y = ylat[ylat_idx]
+        X = xlon[xlon_idx]
+        values = elev[ylat_idx, xlon_idx]
+
+        # Fit: H = a+bx+cy+dxy
+        coeffs = self.polyfit2d(X, Y, values)
+        values = values - (coeffs[0] + coeffs[1] * X +
+                           coeffs[2] * Y + coeffs[3] * Y * X)
+        stddev = np.std(values)
+
+        # Bed Slope
+        nnlat = ylat_nn[edge]
+        nnlon = xlon_nn[edge]
+
+        # nndA = dlat * dlon * R * R * np.cos(nnlat)
+        # dA = np.sum(nndA)
+
+        nnbs = bed_slope[nnlat, nnlon]
+        # nnel = elev[nnlat, nnlon]
+        nnxg = xgrad[nnlat, nnlon]
+        nnyg = ygrad[nnlat, nnlon]
+
+        mean_slope = np.mean(nnbs)
+        mean_xgrad = np.mean(nnxg)
+        mean_ygrad = np.mean(nnyg)
+        # mean_slope= np.sum(nnbs*nndA)/dA
+        # mean_elev = np.sum(nnbs*nndA)/dA
+        # mean_xgrad =np.sum(nnxg*nndA)/dA
+        # mean_ygrad = np.sum(nnyg*nndA)/dA
+
+        # return mean_slope, stddev, mean_elev, mean_xgrad, mean_ygrad
+        return mean_slope, stddev, mean_xgrad, mean_ygrad
+
+    def polyfit2d(self, X, Y, Z):
+        """
+        Calculate linear least squares fit for standard deviation calculation
+        """
+        # X, Y = np.meshgrid(x, y, copy=False)
+        # X = X.flatten()
+        # Y = Y.flatten()
+
+        # Fit: H = a+bx+cy+dxy (from Jayne and St. Laurent)
+        A = np.array([X * 0 + 1, X, Y, X * Y]).T
+        # B = Z.flatten()
+        B = Z
+
+        coeffs, __, __, __ = np.linalg.lstsq(A, B, rcond=None)
+
+        return coeffs
+
+    def map_to_r3(self, xlon, ylat, head, tail):
+        """
+        Map lon-lat coordinates to XYZ points. Restricted to the
+        panel LAT[HEAD:TAIL] to manage memory use.
+
+        """
+
+        sinx = np.sin(xlon * np.pi / 180.)
+        cosx = np.cos(xlon * np.pi / 180.)
+        siny = np.sin(ylat * np.pi / 180.)
+        cosy = np.cos(ylat * np.pi / 180.)
+
+        sinu, sinv = np.meshgrid(sinx, siny[head:tail])
+        cosu, cosv = np.meshgrid(cosx, cosy[head:tail])
+
+        rsph = 6378206.4
+
+        xpos = rsph * cosu * cosv
+        ypos = rsph * sinu * cosv
+        zpos = rsph * sinv
+
+        return np.vstack(
+            (xpos.ravel(), ypos.ravel(), zpos.ravel())).T
diff --git a/compass/ocean/tests/tides/init/dem_remap.py b/compass/ocean/tests/tides/init/dem_remap.py
deleted file mode 100644
index 7c67fe7065..0000000000
--- a/compass/ocean/tests/tides/init/dem_remap.py
+++ /dev/null
@@ -1,344 +0,0 @@
-
-import time
-import numpy as np
-import netCDF4 as nc
-from scipy import spatial
-from scipy.sparse import csr_matrix
-import argparse
-
-# Authors: Darren Engwirda
-
-
-def map_to_r3(mesh, xlon, ylat, head, tail):
-    """
-    Map lon-lat coordinates to XYZ points. Restricted to the
-    panel LAT[HEAD:TAIL] to manage memory use.
-
-    """
-
-    sinx = np.sin(xlon * np.pi / 180.)
-    cosx = np.cos(xlon * np.pi / 180.)
-    siny = np.sin(ylat * np.pi / 180.)
-    cosy = np.cos(ylat * np.pi / 180.)
-
-    sinu, sinv = np.meshgrid(sinx, siny[head:tail])
-    cosu, cosv = np.meshgrid(cosx, cosy[head:tail])
-
-    rsph = mesh.sphere_radius
-
-    xpos = rsph * cosu * cosv
-    ypos = rsph * sinu * cosv
-    zpos = rsph * sinv
-
-    return np.vstack(
-        (xpos.ravel(), ypos.ravel(), zpos.ravel())).T
-
-
-def tria_area(rs, pa, pb, pc):
-    """
-    Calculate areas of spherical triangles [PA, PB, PC] on a
-    sphere of radius RS.
-
-    """
-
-    lena = circ_dist(1., pa, pb)
-    lenb = circ_dist(1., pb, pc)
-    lenc = circ_dist(1., pc, pa)
-
-    slen = 0.5 * (lena + lenb + lenc)
-
-    tana = np.tan(0.5 * (slen - lena))
-    tanb = np.tan(0.5 * (slen - lenb))
-    tanc = np.tan(0.5 * (slen - lenc))
-
-    edel = 4.0 * np.arctan(np.sqrt(
-        np.tan(0.5 * slen) * tana * tanb * tanc))
-
-    return edel * rs ** 2
-
-
-def circ_dist(rs, pa, pb):
-    """
-    Calculate geodesic-length of great circles [PA, PB] on a
-    sphere of radius RS.
-
-    """
-
-    dlon = .5 * (pa[:, 0] - pb[:, 0])
-    dlat = .5 * (pa[:, 1] - pb[:, 1])
-
-    dist = 2. * rs * np.arcsin(np.sqrt(
-        np.sin(dlat) ** 2 +
-        np.sin(dlon) ** 2 * np.cos(pa[:, 1]) * np.cos(pb[:, 1])
-    ))
-
-    return dist
-
-
-def find_vals(xlon, ylat, vals, xpos, ypos):
-
-    cols = xlon.size - 1
-    rows = ylat.size - 1
-
-    dlon = (xlon[-1] - xlon[+0]) / cols
-    dlat = (ylat[-1] - ylat[+0]) / rows
-
-    icol = (xpos - np.min(xlon)) / dlon
-    irow = (ypos - np.min(ylat)) / dlat
-
-    icol = np.asarray(icol, dtype=np.uint32)
-    irow = np.asarray(irow, dtype=np.uint32)
-
-    return vals[irow, icol]
-
-
-def cell_quad(mesh, xlon, ylat, vals):
-    """
-    Eval. finite-volume integrals for mesh cells - splitting
-    each into a series of triangles and employing an order-2
-    quadrature rule.
-
-    """
-
-    class base:
-        pass
-
-    abar = np.zeros(
-        (mesh.dimensions["nCells"].size, 1), dtype=np.float32)
-    fbar = np.zeros(
-        (mesh.dimensions["nCells"].size, 1), dtype=np.float32)
-
-    pcel = np.zeros(
-        (mesh.dimensions["nCells"].size, 2), dtype=np.float64)
-    pcel[:, 0] = mesh.variables["lonCell"][:]
-    pcel[:, 1] = mesh.variables["latCell"][:]
-
-    pcel = pcel * 180. / np.pi
-    pcel[pcel[:, 0] > 180., 0] -= 360.
-
-    fcel = find_vals(xlon, ylat, vals, pcel[:, 0], pcel[:, 1])
-
-    pcel = pcel * np.pi / 180.
-
-    pvrt = np.zeros(
-        (mesh.dimensions["nVertices"].size, 2),
-        dtype=np.float64)
-    pvrt[:, 0] = mesh.variables["lonVertex"][:]
-    pvrt[:, 1] = mesh.variables["latVertex"][:]
-
-    pvrt = pvrt * 180. / np.pi
-    pvrt[pvrt[:, 0] > 180., 0] -= 360.
-
-    fvrt = find_vals(xlon, ylat, vals, pvrt[:, 0], pvrt[:, 1])
-
-    pvrt = pvrt * np.pi / 180.
-
-    cell = base()
-    cell.edge = np.asarray(
-        mesh.variables["edgesOnCell"][:], dtype=np.int32)
-    cell.topo = np.asarray(
-        mesh.variables["nEdgesOnCell"][:], dtype=np.int32)
-
-    edge = base()
-    edge.vert = np.asarray(
-        mesh.variables["verticesOnEdge"][:], dtype=np.int32)
-
-    for epos in range(np.max(cell.topo)):
-
-        mask = cell.topo > epos
-
-        icel = np.argwhere(mask).ravel()
-
-        ifac = cell.edge[mask, epos] - 1
-
-        ivrt = edge.vert[ifac, 0] - 1
-        jvrt = edge.vert[ifac, 1] - 1
-
-        rsph = mesh.sphere_radius
-
-        atri = tria_area(
-            rsph, pcel[icel], pvrt[ivrt], pvrt[jvrt])
-
-        atri = np.reshape(atri, (atri.size, 1))
-
-        ftri = (fcel[icel] + fvrt[ivrt] + fvrt[jvrt])
-
-        ftri = np.reshape(ftri, (ftri.size, 1))
-
-        abar[icel] += atri
-        fbar[icel] += atri * ftri / 3.
-
-    return fbar / abar
-
-
-def dem_remap(elev_file, mpas_file):
-    """
-    Map elevation and ice+ocn-thickness data from a "zipped"
-    RTopo data-set onto the cells in an MPAS mesh.
-
-    Cell values are a blending of an approx. integral remap
-    and a local quadrature rule.
-
-    """
-
-    print("Loading assests...")
-
-    elev = nc.Dataset(elev_file, "r")
-    mesh = nc.Dataset(mpas_file, "r+")
-
-# -- Compute an approximate remapping, associating pixels in
-# -- the DEM with cells in the MPAS mesh. Since polygons are
-# -- Voronoi, the point-in-cell query can be computed by
-# -- finding nearest neighbours. This remapping is an approx.
-# -- as no partial pixel-cell intersection is computed.
-
-    print("Building KDtree...")
-
-    ppos = np.zeros(
-        (mesh.dimensions["nCells"].size, 3), dtype=np.float64)
-    ppos[:, 0] = mesh["xCell"][:]
-    ppos[:, 1] = mesh["yCell"][:]
-    ppos[:, 2] = mesh["zCell"][:]
-
-    tree = spatial.cKDTree(ppos, leafsize=8)
-
-    print("Remap elevation...")
-
-    xlon = np.asarray(elev["lon"][:], dtype=np.float64)
-    ylat = np.asarray(elev["lat"][:], dtype=np.float64)
-
-    xmid = .5 * (xlon[:-1:] + xlon[1::])
-    ymid = .5 * (ylat[:-1:] + ylat[1::])
-
-    indx = np.asarray(np.round(
-        np.linspace(-1, ymid.size, 9)), dtype=np.int32)
-
-    print("* process tiles:")
-
-    nset = []
-    for tile in range(indx.size - 1):
-
-        head = indx[tile + 0] + 1
-        tail = indx[tile + 1] + 1
-
-        qpos = map_to_r3(mesh, xmid, ymid, head, tail)
-
-        ttic = time.time()
-        __, nloc = tree.query(qpos, workers=-1)
-        ttoc = time.time()
-        print("* built node-to-cell map:", ttoc - ttic)
-
-        nset.append(nloc)
-
-    near = np.concatenate(nset)
-
-    del tree
-    del qpos
-    del nset
-    del nloc
-
-# -- Build cell-to-pixel map as a sparse matrix, and average
-# -- RTopo pixel values within each cell.
-
-    cols = np.arange(0, near.size)
-    vals = np.ones(near.size, dtype=np.int8)
-
-    smat = csr_matrix((vals, (near, cols)))
-
-    nmap = np.asarray(
-        np.sum(smat, axis=1), dtype=np.float32)
-
-    vals = np.asarray(
-        elev["bed_elevation"][:], dtype=np.float32)
-    vals = np.reshape(vals, (vals.size, 1))
-
-    emap = (smat * vals) / np.maximum(1., nmap)
-
-    vals = np.asarray(
-        elev["ocn_thickness"][:], dtype=np.float32)
-    vals = np.reshape(vals, (vals.size, 1))
-
-    omap = (smat * vals) / np.maximum(1., nmap)
-
-    vals = np.asarray(
-        elev["ice_thickness"][:], dtype=np.float32)
-    vals = np.reshape(vals, (vals.size, 1))
-
-    imap = (smat * vals) / np.maximum(1., nmap)
-
-    del smat
-    del cols
-    del vals
-    del near
-
-# -- If the resolution of the mesh is greater, or comparable
-# -- to that of the DEM, the approx. remapping (above) will
-# -- result in a low order interpolation.
-# -- Thus, blend with a local order-2 quadrature formulation
-
-    print("Eval. elevation...")
-
-    vals = np.asarray(
-        elev["bed_elevation"][:], dtype=np.float32)
-
-    ttic = time.time()
-    eint = cell_quad(mesh, xlon, ylat, vals)
-    ttoc = time.time()
-    print("* compute cell integrals:", ttoc - ttic)
-
-    vals = np.asarray(
-        elev["ocn_thickness"][:], dtype=np.float32)
-
-    ttic = time.time()
-    oint = cell_quad(mesh, xlon, ylat, vals)
-    ttoc = time.time()
-    print("* compute cell integrals:", ttoc - ttic)
-
-    vals = np.asarray(
-        elev["ice_thickness"][:], dtype=np.float32)
-
-    ttic = time.time()
-    iint = cell_quad(mesh, xlon, ylat, vals)
-    ttoc = time.time()
-    print("* compute cell integrals:", ttoc - ttic)
-
-    print("Save to dataset...")
-
-    ebar = (np.multiply(nmap, emap) + 6 * eint) / (6 + nmap)
-    obar = (np.multiply(nmap, omap) + 6 * oint) / (6 + nmap)
-    ibar = (np.multiply(nmap, imap) + 6 * iint) / (6 + nmap)
-
-    if ("bed_elevation" not in mesh.variables.keys()):
-        mesh.createVariable("bed_elevation", "f8", ("nCells"))
-
-    if ("ocn_thickness" not in mesh.variables.keys()):
-        mesh.createVariable("ocn_thickness", "f8", ("nCells"))
-
-    if ("ice_thickness" not in mesh.variables.keys()):
-        mesh.createVariable("ice_thickness", "f8", ("nCells"))
-
-    mesh["bed_elevation"][:] = ebar
-    mesh["ocn_thickness"][:] = obar
-    mesh["ice_thickness"][:] = ibar
-
-    elev.close()
-    mesh.close()
-
-
-if (__name__ == "__main__"):
-    parser = argparse.ArgumentParser(
-        description=__doc__,
-        formatter_class=argparse.RawTextHelpFormatter)
-
-    parser.add_argument(
-        "--mpas-file", dest="mpas_file", type=str,
-        required=True, help="Name of user MPAS mesh.")
-
-    parser.add_argument(
-        "--elev-file", dest="elev_file", type=str,
-        required=True, help="Name of DEM pixel file.")
-
-    args = parser.parse_args()
-    elev_file = args.elev_file
-    mpas_file = args.mpas_file
-    dem_remap(elev_file, mpas_file)
diff --git a/compass/ocean/tests/tides/init/dem_trnsf.py b/compass/ocean/tests/tides/init/dem_trnsf.py
deleted file mode 100644
index 1c41a2bacb..0000000000
--- a/compass/ocean/tests/tides/init/dem_trnsf.py
+++ /dev/null
@@ -1,114 +0,0 @@
-
-import numpy as np
-import netCDF4 as nc
-import argparse
-
-# Authors: Darren Engwirda
-
-
-def dem_trnsf(base_mesh, part_mesh):
-    """
-    Transfer elevation & ice+ocn-thickness data from a full
-    sphere MPAS mesh onto a partial sub-mesh resulting from
-    a culling operation, or equiv.
-
-    The sub-mesh is expected to be a piece of the base mesh
-    such that (some) cells match 1:1.
-
-    """
-
-    base = nc.Dataset(base_mesh, "r")
-    part = nc.Dataset(part_mesh, "r+")
-
-# -- make a vector of cell-centre positions to match against
-
-    xpos = np.vstack((
-        np.vstack((
-            base["xCell"][:],
-            base["yCell"][:],
-            base["zCell"][:])).T,
-        np.vstack((
-            part["xCell"][:],
-            part["yCell"][:],
-            part["zCell"][:])).T
-    ))
-
-    xidx = np.hstack((
-        np.arange(
-            0, base.dimensions["nCells"].size),
-        np.arange(
-            0, part.dimensions["nCells"].size)
-    ))
-
-# -- culling shouldn't introduce fp round-off - but truncate
-# -- anyway...
-
-    xpos = np.round(xpos, decimals=8)
-
-# -- use stable sorting to bring matching cell xyz (and idx)
-# -- into "ascending" order
-
-    imap = np.argsort(xpos[:, 2], kind="stable")
-    xpos = xpos[imap, :]
-    xidx = xidx[imap]
-    imap = np.argsort(xpos[:, 1], kind="stable")
-    xpos = xpos[imap, :]
-    xidx = xidx[imap]
-    imap = np.argsort(xpos[:, 0], kind="stable")
-    xpos = xpos[imap, :]
-    xidx = xidx[imap]
-
-    diff = xpos[1:, :] - xpos[:-1, :]
-
-    same = np.argwhere(np.logical_and.reduce((
-        diff[:, 0] == 0,
-        diff[:, 1] == 0,
-        diff[:, 2] == 0))).ravel()
-
-# -- cell inew in part matches iold in base - re-index elev.
-# -- data-sets
-
-    inew = xidx[same + 1]
-    iold = xidx[same + 0]
-
-    if ("bed_elevation" not in base.variables.keys() or
-            "ocn_thickness" not in base.variables.keys() or
-            "ice_thickness" not in base.variables.keys()):
-        raise Exception("Base does not contain elev. data!")
-
-    if ("bed_elevation" not in part.variables.keys()):
-        part.createVariable("bed_elevation", "f8", ("nCells"))
-
-    if ("ocn_thickness" not in part.variables.keys()):
-        part.createVariable("ocn_thickness", "f8", ("nCells"))
-
-    if ("ice_thickness" not in part.variables.keys()):
-        part.createVariable("ice_thickness", "f8", ("nCells"))
-
-    part["bed_elevation"][inew] = base["bed_elevation"][iold]
-    part["ocn_thickness"][inew] = base["ocn_thickness"][iold]
-    part["ice_thickness"][inew] = base["ice_thickness"][iold]
-
-    base.close()
-    part.close()
-
-    return
-
-
-if (__name__ == "__main__"):
-    parser = argparse.ArgumentParser(
-        description=__doc__,
-        formatter_class=argparse.RawTextHelpFormatter)
-
-    parser.add_argument(
-        "--base-mesh", dest="base_mesh", type=str,
-        required=True, help="Name of (full) MPAS mesh.")
-
-    parser.add_argument(
-        "--part-mesh", dest="part_mesh", type=str,
-        required=True, help="Name of culled MPAS mesh.")
-
-    args = parser.parse_args()
-    base_mesh = args.base_mesh
-    part_mesh = args.part_mesh
-    dem_trnsf(base_mesh, part_mesh)
diff --git a/compass/ocean/tests/tides/init/initial_state.py b/compass/ocean/tests/tides/init/initial_state.py
index 2e744f8776..a604bfed1f 100644
--- a/compass/ocean/tests/tides/init/initial_state.py
+++ b/compass/ocean/tests/tides/init/initial_state.py
@@ -1,4 +1,4 @@
-import netCDF4
+import netCDF4 as nc
 import numpy as np
 
 from compass.model import run_model
@@ -113,13 +113,89 @@ def run(self):
         max_depth = config.getfloat('vertical_grid', 'bottom_depth')
         min_depth = config.getfloat('vertical_grid', 'min_depth')
 
-        init = netCDF4.Dataset("initial_state.nc", "r+")
+        init = nc.Dataset("initial_state.nc", "r+")
         init["bottomDepth"][:] = \
             np.minimum(max_depth,
                        np.maximum(min_depth, init["bottomDepthObserved"][:]))
 
         init["layerThickness"][0, :, 0] = init["bottomDepth"][:]
         init["restingThickness"][:, 0] = init["bottomDepth"][:]
+
+        # -- Estimate vert. grid for ice-shelves, min.-thicknesses, etc
+        # -- Darren Engwirda
+
+        print("Est. layering to account for ice-shelves")
+
+        mesh = nc.Dataset("mesh.nc", "r")
+
+        botd = np.asarray(init["bottomDepthObserved"][:], dtype=np.float64)
+        # ossh = np.asarray(init["ssh"][0,:], dtype=np.float64)
+        ossh = 0. * botd  # assume ssh is zero
+
+        grav = 9.80665  # gravitational accel.
+        irho = float(init.config_land_ice_flux_rho_ice)
+        orho = float(init.config_density0)
+        minh = float(init.config_drying_min_cell_height) / 2.
+
+        print("ice-shelf density:", irho)
+        print("ocn-const density:", orho)
+        print("min-layer thickness:", minh)
+
+        iceh = np.asarray(mesh["ice_thickness"][:], dtype=np.float64)
+        # icef = np.asarray(mesh["ice_cover"][:], dtype=np.float64)
+
+        icep = irho * grav * iceh  # ice pressure
+        iced = icep / grav / orho  # ice draft
+
+        # ensure thin-layer beneath ice-shelves
+        iced = np.minimum(iced, +botd - minh)
+        iced = np.maximum(iced, +0.0)
+        ossh = ossh - iced
+
+        icep[iced <= 0.] = 0.
+
+        # allow thin-layer in partially flooded zone
+        ossh = np.maximum(ossh, -botd + minh)
+
+        print("max ice-draft:", np.max(iced))
+        print("max ice-pressure:", np.max(icep))
+
+        if ("ssh" not in init.variables.keys()):
+            init.createVariable("ssh", "f8", ("Time", "nCells"))
+
+        if ("landIceDraft" not in init.variables.keys()):
+            init.createVariable("landIceDraft", "f8", ("Time", "nCells"))
+
+        if ("landIcePressure" not in init.variables.keys()):
+            init.createVariable("landIcePressure", "f8", ("Time", "nCells"))
+
+        if ("landIceMask" not in init.variables.keys()):
+            init.createVariable("landIceMask", "i4", ("Time", "nCells"))
+
+        if ("landIceFloatingMask" not in init.variables.keys()):
+            init.createVariable("landIceFloatingMask", "i4",
+                                ("Time", "nCells"))
+
+        if ("landIceFraction" not in init.variables.keys()):
+            init.createVariable("landIceFraction", "f8", ("Time", "nCells"))
+
+        if ("landIceFloatingFraction" not in init.variables.keys()):
+            init.createVariable("landIceFloatingFraction", "f8",
+                                ("Time", "nCells"))
+
+        init["landIceDraft"][0, :] = -iced  # NB. sign
+        init["landIcePressure"][0, :] = icep
+        init["landIceMask"][0, :] = (icep > 0.)
+        init["landIceFloatingMask"][0, :] = (icep > 0.)
+        init["landIceFraction"][0, :] = (icep > 0.)
+        init["landIceFloatingFraction"][0, :] = (icep > 0.)
+
+        init["bottomDepth"][:] = botd
+        init["ssh"][0, :] = ossh
+
+        init["layerThickness"][0, :, 0] = ossh + botd
+        # init["restingThickness"][:,0] = ossh + botd
+
         init.close()
 
     def _get_resources(self):
diff --git a/compass/ocean/tests/tides/init/namelist.init b/compass/ocean/tests/tides/init/namelist.init
index 853428f1c0..d36a53859d 100644
--- a/compass/ocean/tests/tides/init/namelist.init
+++ b/compass/ocean/tests/tides/init/namelist.init
@@ -4,7 +4,7 @@ config_vert_levels = -1
 config_write_cull_cell_mask = .false.
 config_expand_sphere = .true.
 config_realistic_coriolis_parameter = .true.
-config_global_ocean_minimum_depth = 0.1 
+config_global_ocean_minimum_depth = 0.01
 config_global_ocean_topography_source = 'mpas_variable'
 config_global_ocean_depth_file = 'vertical_grid.nc'
 config_global_ocean_temperature_file = 'temperature.nc'
@@ -20,7 +20,7 @@ config_global_ocean_salinity_varname = 'SALT'
 config_global_ocean_tracer_lat_varname = 't_lat'
 config_global_ocean_tracer_lon_varname = 't_lon'
 config_global_ocean_tracer_depth_varname = 'depth_t'
-config_global_ocean_smooth_TS_iterations = 5 
+config_global_ocean_smooth_TS_iterations = 2
 config_global_ocean_swData_nlat_dimname = 't_lat'
 config_global_ocean_swData_nlon_dimname = 't_lon'
 config_global_ocean_swData_lat_varname = 't_lat'
@@ -49,3 +49,8 @@ config_global_ocean_ecosys_ndepth_dimname = 'depth_t'
 config_global_ocean_ecosys_lat_varname = 't_lat'
 config_global_ocean_ecosys_lon_varname = 't_lon'
 config_global_ocean_ecosys_depth_varname = 'depth_t'
+config_global_ocean_piston_velocity = 1.5856e-05
+config_use_wetting_drying = .true.
+config_prevent_drying = .true.
+config_drying_min_cell_height = 0.01
+config_land_ice_flux_mode = "pressure_only"
diff --git a/compass/ocean/tests/tides/init/remap_bathymetry.py b/compass/ocean/tests/tides/init/remap_bathymetry.py
index 6f3a077507..f6d19ee8fc 100644
--- a/compass/ocean/tests/tides/init/remap_bathymetry.py
+++ b/compass/ocean/tests/tides/init/remap_bathymetry.py
@@ -1,11 +1,7 @@
+import compass.ocean.tests.tides.dem.dem_remap as dem_remap
+import compass.ocean.tests.tides.dem.dem_trnsf as dem_trnsf
 from compass.step import Step
 
-import os
-
-import compass.ocean.tests.tides.init.dem_pixel as dem_pixel
-import compass.ocean.tests.tides.init.dem_remap as dem_remap
-import compass.ocean.tests.tides.init.dem_trnsf as dem_trnsf
-
 
 class RemapBathymetry(Step):
     """
@@ -33,28 +29,9 @@ def __init__(self, test_case, mesh):
 
         self.plot = True
 
-        self.add_input_file(
-            filename='GEBCO_2020.nc',
-            target='GEBCO_2020.nc',
-            database='bathymetry_database')
-
-        self.add_input_file(
-            filename='RTopo-2.0.4_30sec_bedrock_topography.nc',
-            target='RTopo-2.0.4_30sec_bedrock_topography.nc',
-            database='bathymetry_database')
-
-        self.add_input_file(
-            filename='RTopo-2.0.4_30sec_surface_elevation.nc',
-            target='RTopo-2.0.4_30sec_surface_elevation.nc',
-            database='bathymetry_database')
-
-        self.add_input_file(
-            filename='RTopo-2.0.4_30sec_ice_base_topography.nc',
-            target='RTopo-2.0.4_30sec_ice_base_topography.nc',
-            database='bathymetry_database')
-
         cull_mesh_path = mesh.steps['cull_mesh'].path
         base_mesh_path = mesh.steps['base_mesh'].path
+        pixel_path = mesh.steps['pixel'].path
 
         self.add_input_file(
             filename='mesh.nc',
@@ -62,21 +39,15 @@ def __init__(self, test_case, mesh):
         self.add_input_file(
             filename='base_mesh.nc',
             work_dir_target=f'{base_mesh_path}/base_mesh.nc')
+        self.add_input_file(
+            filename='bathy.nc',
+            work_dir_target=f'{pixel_path}/'
+                            'RTopo_2_0_4_GEBCO_v2023_30sec_pixel.nc')
 
     def run(self):
         """
         Run this step of the test case
         """
-
-        self.init_path = './'
-
-        if not os.path.exists('RTopo_2_0_4_30sec_pixel.nc'):
-            dem_pixel.rtopo_30sec(self.init_path, self.init_path)
-        if not os.path.exists('GEBCO_v2020_30sec_pixel.nc'):
-            dem_pixel.gebco_30sec(self.init_path, self.init_path)
-        if not os.path.exists('RTopo_2_0_4_GEBCO_v2020_30sec_pixel.nc'):
-            dem_pixel.rtopo_gebco_30sec(self.init_path, self.init_path)
-
-        dem_remap.dem_remap('RTopo_2_0_4_GEBCO_v2020_30sec_pixel.nc',
+        dem_remap.dem_remap('bathy.nc',
                             'base_mesh.nc')
         dem_trnsf.dem_trnsf('base_mesh.nc', 'mesh.nc')
diff --git a/compass/ocean/tests/tides/mesh/__init__.py b/compass/ocean/tests/tides/mesh/__init__.py
index 96d0663877..f4fe523197 100644
--- a/compass/ocean/tests/tides/mesh/__init__.py
+++ b/compass/ocean/tests/tides/mesh/__init__.py
@@ -1,7 +1,9 @@
-from compass.testcase import TestCase
 from compass.mesh.spherical import IcosahedralMeshStep
-from compass.ocean.tests.tides.configure import configure_tides
 from compass.ocean.mesh.cull import CullMeshStep
+from compass.ocean.tests.tides.configure import configure_tides
+from compass.ocean.tests.tides.dem import CreatePixelFile
+from compass.ocean.tests.tides.mesh.vr45to5 import VRTidesMesh
+from compass.testcase import TestCase
 
 
 class Mesh(TestCase):
@@ -25,11 +27,23 @@ def __init__(self, test_group, mesh_name):
         subdir = '{}/{}'.format(mesh_name, name)
         super().__init__(test_group=test_group, name=name, subdir=subdir)
 
+        pixel_step = CreatePixelFile(self)
+        self.add_step(pixel_step)
+
         name = 'base_mesh'
         if mesh_name == 'Icos7':
             base_mesh_step = IcosahedralMeshStep(
                 self, name=name, subdivisions=7)
             mesh_lower = 'icos7'
+        elif mesh_name == 'VR45to5':
+            base_mesh_step = VRTidesMesh(
+                self, pixel_step,
+                name='base_mesh', subdir=None,
+                elev_file='RTopo_2_0_4_GEBCO_v2023_30sec_pixel.nc',
+                spac_dhdx=0.125, spac_hmin=5, spac_hmax=45, spac_hbar=45,
+                ncell_nwav=80, ncell_nslp=4,
+                filt_sdev=0.5, filt_halo=50, filt_plev=0.325)
+            mesh_lower = 'vr45to5'
         else:
             raise ValueError(f'Unexpected mesh name {mesh_name}')
 
diff --git a/compass/ocean/tests/tides/mesh/vr45to5/__init__.py b/compass/ocean/tests/tides/mesh/vr45to5/__init__.py
new file mode 100644
index 0000000000..40d9701f7d
--- /dev/null
+++ b/compass/ocean/tests/tides/mesh/vr45to5/__init__.py
@@ -0,0 +1,233 @@
+import os
+
+import jigsawpy
+import netCDF4 as nc
+import numpy as np
+from mpas_tools.cime.constants import constants
+from mpas_tools.logging import check_call
+from skimage.filters import gaussian, scharr
+
+from compass.mesh.spherical import QuasiUniformSphericalMeshStep
+
+# from skimage.filters import farid
+# from skimage.filters.rank import median, percentile
+
+
+class VRTidesMesh(QuasiUniformSphericalMeshStep):
+    """
+    A step for creating a variable resolution tides mesh
+    """
+
+    def __init__(self, test_case, pixel, name='base_mesh', subdir=None,
+                 elev_file='RTopo_2_0_4_GEBCO_v2023_30sec_pixel.nc',
+                 spac_dhdx=0.1, spac_hmin=10.0, spac_hmax=75.0, spac_hbar=60.0,
+                 ncell_nwav=60, ncell_nslp=0,
+                 filt_sdev=3.0, filt_halo=50, filt_plev=0.325):
+        """
+
+        """
+
+        super().__init__(test_case=test_case, name=name, subdir=subdir)
+
+        self.spac_dhdx = spac_dhdx
+        self.spac_hmin = spac_hmin
+        self.spac_hmax = spac_hmax
+        self.spac_hbar = spac_hbar
+
+        self.filt_halo = filt_halo
+        self.filt_sdev = filt_sdev
+        self.filt_plev = filt_plev
+
+        self.ncell_wav = ncell_nwav
+        self.ncell_slp = ncell_nslp
+
+        self.elev_file = elev_file
+
+        pixel_path = pixel.path
+
+        self.add_input_file(
+            filename='bathy.nc',
+            work_dir_target=f'{pixel_path}/{elev_file}')
+
+    def make_jigsaw_mesh(self, lon, lat, cell_width):
+        """
+        Build the JIGSAW mesh.
+        """
+        self.opts.optm_kern = "cvt+dqdx"
+        self.opts.optm_iter = 32                 # tighter opt. tol
+        self.opts.optm_qtol = +1.00E-05
+
+        super().make_jigsaw_mesh(lon, lat, cell_width)
+
+    def build_cell_width_lat_lon(self):
+        """
+
+        """
+        spac = jigsawpy.jigsaw_msh_t()
+
+        # ------------------------------------ define spacing pattern
+
+        dhdx = self.spac_dhdx   # max allowable slope in spacing
+        hmin = self.spac_hmin   # min mesh-spacing value (km)
+        hmax = self.spac_hmax   # max mesh-spacing value (km)
+        hbar = self.spac_hbar   # constant spacing value (km)
+
+        halo = self.filt_halo   # DEM pixels per filtering radii
+        sdev = self.filt_sdev   # std-dev for gaussian filter
+        plev = self.filt_plev   # median-style filter percentile
+
+        nwav = self.ncell_wav   # number of cells per wavelength
+        nslp = self.ncell_slp   # number of cells per grad(elev)
+
+        print("Loading elevation assets...")
+
+        data = nc.Dataset("bathy.nc", "r")
+
+        ocnh = np.asarray(
+            data["ocn_thickness"][:])
+
+        print("Computing background h(x)...")
+
+        hmat = np.full(
+            (ocnh.shape[:]), (hbar))
+
+        if (nwav >= 1):
+            hmat = np.minimum(
+                hmat, self.swe_wavelength_spacing(
+                    ocnh, nwav, hmin, hmax, halo, plev))
+
+        if (nslp >= 1):
+            hmat = np.minimum(
+                hmat, self.elev_sharpness_spacing(
+                    ocnh, nslp, hmin, hmax, halo, plev,
+                    sdev))
+
+        hmat[ocnh <= 0.] = hmax
+        hmat = np.maximum(hmat, hmin)
+        hmat = np.minimum(hmat, hmax)
+
+        # -- pack h(x) data into jigsaw data-type: average pixel-to-
+        # -- node, careful with periodic BC's.
+
+        hmat = self.coarsen_spacing_pixels(hmat, down=4)
+
+        FULL_SPHERE_RADIUS = constants["SHR_CONST_REARTH"] / 1.E+003
+
+        spac.mshID = "ellipsoid-grid"       # use the elv. grid
+        spac.radii = np.full(
+            3, FULL_SPHERE_RADIUS, dtype=spac.REALS_t)
+
+        spac.xgrid = np.linspace(
+            -1. * np.pi, +1. * np.pi, hmat.shape[1] + 1)
+
+        spac.ygrid = np.linspace(
+            -.5 * np.pi, +.5 * np.pi, hmat.shape[0] + 1)
+
+        R = hmat.shape[0]
+        C = hmat.shape[1]
+
+        spac.value = np.zeros(
+            (R + 1, C + 1))
+
+        npos = np.arange(+0, hmat.shape[0] + 1)
+        epos = np.arange(-1, hmat.shape[1] - 0)
+        spos = np.arange(-1, hmat.shape[0] - 0)
+        wpos = np.arange(+0, hmat.shape[1] + 1)
+
+        npos[npos >= +R] = R - 1
+        spos[spos <= -1] = +0
+        epos[epos <= -1] = C - 1
+        wpos[wpos >= +C] = +0
+
+        npos, epos = np.meshgrid(
+            npos, epos, sparse=True, indexing="ij")
+        spos, wpos = np.meshgrid(
+            spos, wpos, sparse=True, indexing="ij")
+
+        spac.value += hmat[npos, epos] * (+1. / 4.)
+        spac.value += hmat[npos, wpos] * (+1. / 4.)
+        spac.value += hmat[spos, epos] * (+1. / 4.)
+        spac.value += hmat[spos, wpos] * (+1. / 4.)
+
+        spac = self.limit_spacing_gradient(spac, dhdx=dhdx)
+
+        return spac.value, np.degrees(spac.xgrid), np.degrees(spac.ygrid)
+
+    def limit_spacing_gradient(self, spac, dhdx):
+
+        print("Smoothing h(x) via |dh/dx| limits...")
+
+        opts = jigsawpy.jigsaw_jig_t()
+
+        spac.slope = np.full(spac.value.shape, dhdx)
+        opts.hfun_file = os.path.join('.', "spac_pre_smooth.msh")
+        opts.jcfg_file = os.path.join('.', "opts_pre_smooth.jig")
+        jigsawpy.savemsh(opts.hfun_file, spac)
+
+        opts.verbosity = +1
+
+        jigsawpy.cmd.marche(opts, spac)
+
+        return spac
+
+    def swe_wavelength_spacing(self,
+                               ocnh, nwav, hmin, hmax, halo, plev,
+                               T_M2=12. * 3600., grav=9.806):
+
+        print("Computing wavelength heuristic...")
+
+        vals = T_M2 * np.sqrt(
+            grav * np.maximum(5, ocnh)) / nwav / 1000.
+
+        vals[ocnh <= 0.] = hmax
+        vals = np.maximum(vals, hmin)
+        vals = np.minimum(vals, hmax)
+
+        vals = np.asarray(vals, dtype=np.uint16)
+        # vals = percentile(
+        #     vals, footprint=disk(halo), mask=(ocnh>0.), p0=plev)
+
+        return vals
+
+    def elev_sharpness_spacing(self,
+                               ocnh, nslp, hmin, hmax, halo, plev, sdev):
+
+        print("Computing GRAD(elev) heuristic...")
+
+        dzdx = scharr(gaussian(np.asarray(
+            ocnh, dtype=np.float32), sigma=sdev, mode="wrap"))
+
+        dzdx = np.maximum(1.E-08, dzdx)  # no divide-by-zero
+
+        vals = np.maximum(
+            5., np.abs(ocnh)) / dzdx * 2. * np.pi / nslp
+
+        vals = np.maximum(vals, hmin)
+        vals = np.minimum(vals, hmax)
+
+        vals = np.asarray(vals, dtype=np.uint16)
+        # vals = percentile(
+        #     vals, footprint=disk(halo), mask=(ocnh>0.), p0=plev)
+
+        return vals
+
+    def coarsen_spacing_pixels(self, hmat, down):
+
+        print("Coarsening mesh-spacing pixels...")
+
+        rows = hmat.shape[0] // down
+        cols = hmat.shape[1] // down
+
+        htmp = np.full(
+            (rows, cols), (np.amax(hmat)), dtype=hmat.dtype)
+
+        for jpos in range(down):
+            for ipos in range(down):
+
+                iend = hmat.shape[0] - down + ipos + 1
+                jend = hmat.shape[1] - down + jpos + 1
+
+                htmp = np.minimum(htmp,
+                                  hmat[ipos:iend:down, jpos:jend:down])
+
+        return htmp
diff --git a/compass/ocean/tests/tides/mesh/vr45to5/namelist.ocean b/compass/ocean/tests/tides/mesh/vr45to5/namelist.ocean
new file mode 100644
index 0000000000..a4d3c261c8
--- /dev/null
+++ b/compass/ocean/tests/tides/mesh/vr45to5/namelist.ocean
@@ -0,0 +1,6 @@
+config_dt = '00:00:20'
+config_topographic_wave_drag_scheme = 'ZAE'
+config_topographic_wave_drag_coeff = 0.3
+config_topographic_wave_drag_cutoff_depth = 500
+config_topographic_wave_drag_cutoff_width = 10
+config_drying_min_cell_height = 0.01
diff --git a/compass/ocean/tests/tides/mesh/vr45to5/vr45to5.cfg b/compass/ocean/tests/tides/mesh/vr45to5/vr45to5.cfg
new file mode 100644
index 0000000000..6ddc17404d
--- /dev/null
+++ b/compass/ocean/tests/tides/mesh/vr45to5/vr45to5.cfg
@@ -0,0 +1,62 @@
+# options for spherical meshes
+[spherical_mesh]
+
+## config options related to the step for culling land from the mesh
+# number of cores to use
+cull_mesh_cpus_per_task = 18
+# minimum of cores, below which the step fails
+cull_mesh_min_cpus_per_task = 1
+# maximum memory usage allowed (in MB)
+cull_mesh_max_memory = 1000
+
+# for icosahedral meshes, whether to use cell_width to determine the number of
+# subdivisions or to use subdivisions directly
+icosahedral_method = subdivisions
+# Whether to convert the culled mesh file to CDF5 format
+convert_culled_mesh_to_cdf5 = False
+# Minimum latitude, in degrees, for masking land-locked cells
+latitude_threshold = 43.0
+# Maximum number of sweeps to search for land-locked cells
+sweep_count = 20
+
+
+# options for global ocean testcases
+[global_ocean]
+
+# The following options are detected from .gitconfig if not explicitly entered
+author = autodetect
+email = autodetect
+
+[vertical_grid]
+
+grid_type = uniform
+vert_levels = 1
+bottom_depth = 12000
+min_depth = 0.1
+
+
+# options for tides testcases
+[tides]
+
+## config options related to the initial_state step
+# number of cores to use
+init_ntasks = 36
+# minimum of cores, below which the step fails
+init_min_tasks = 8
+# maximum memory usage allowed (in MB)
+init_max_memory = 1000
+# number of threads
+init_threads = 1
+
+## config options related to the forward steps
+# number of cores to use
+forward_ntasks = 1024
+# minimum of cores, below which the step fails
+forward_min_tasks = 1024
+# maximum memory usage allowed (in MB)
+forward_max_memory = 1000
+# number of threads
+forward_threads = 1
+
+# TPXO version for validation
+tpxo_version = TPXO9
diff --git a/docs/developers_guide/ocean/api.rst b/docs/developers_guide/ocean/api.rst
index edb7822257..39a6a2ee96 100644
--- a/docs/developers_guide/ocean/api.rst
+++ b/docs/developers_guide/ocean/api.rst
@@ -864,6 +864,14 @@ test cases and steps
    mesh.Mesh.configure
    mesh.Mesh.run
 
+   mesh.vr45to5.VRTidesMesh
+   mesh.vr45to5.VRTidesMesh.make_jigsaw_mesh
+   mesh.vr45to5.VRTidesMesh.build_cell_width_lat_lon
+   mesh.vr45to5.VRTidesMesh.limit_spacing_gradient
+   mesh.vr45to5.VRTidesMesh.swe_wavelength_spacing
+   mesh.vr45to5.VRTidesMesh.elev_sharpness_spacing
+   mesh.vr45to5.VRTidesMesh.coarsen_spacing_pixels
+
    init.Init
    init.Init.configure
    init.Init.run
@@ -872,12 +880,41 @@ test cases and steps
    init.initial_state.InitialState
    init.initial_state.InitialState.setup
    init.initial_state.InitialState.run
-   init.interpolate_wave_drag.InterpolateWaveDrag
-   init.interpolate_wave_drag.InterpolateWaveDrag.interpolate_data_to_grid
-   init.interpolate_wave_drag.InterpolateWaveDrag.plot_interp_data
-   init.interpolate_wave_drag.InterpolateWaveDrag.write_to_file
-   init.interpolate_wave_drag.InterpolateWaveDrag.run
-
+   init.calculate_wave_drag.CalculateWaveDrag
+   init.calculate_wave_drag.CalculateWaveDrag.interpolate_data_to_grid
+   init.calculate_wave_drag.CalculateWaveDrag.fix_nans
+   init.calculate_wave_drag.CalculateWaveDrag.make_nn_lists
+   init.calculate_wave_drag.CalculateWaveDrag.calc_stats
+   init.calculate_wave_drag.CalculateWaveDrag.polyfit2d
+   init.calculate_wave_drag.CalculateWaveDrag.map_to_r3
+   init.calculate_wave_drag.CalculateWaveDrag.run
+
+   dem.CreatePixelFile
+
+   dem.dem_pixel.circ_dist
+   dem.dem_pixel.cell_dzdx
+   dem.dem_pixel.blend_front
+   dem.dem_pixel.rtopo_60sec
+   dem.dem_pixel.rtopo_30sec
+   dem.dem_pixel.rtopo_15sec
+   dem.dem_pixel.gebco_60sec
+   dem.dem_pixel.gebco_30sec
+   dem.dem_pixel.gebco_15sec
+   dem.dem_pixel.rtopo_gebco_60sec
+   dem.dem_pixel.rtopo_gebco_30sec
+   dem.dem_pixel.rtopo_gebco_15sec
+
+   dem.dem_remap.map_to_r3
+   dem.dem_remap.tria_area
+   dem.dem_remap.circ_dist
+   dem.dem_remap.sample_1d
+   dem.dem_remap.linear_2d
+   dem.dem_remap.sample_2d
+   dem.dem_remap.cell_quad
+   dem.dem_remap.cell_prfl
+   dem.dem_remap.dem_remap
+
+   dem.dem_trnsf.dem_trnsf
 
    forward.Forward
    forward.Forward.configure
diff --git a/docs/developers_guide/ocean/test_groups/tides.rst b/docs/developers_guide/ocean/test_groups/tides.rst
index f5418f5384..ed41a078ad 100644
--- a/docs/developers_guide/ocean/test_groups/tides.rst
+++ b/docs/developers_guide/ocean/test_groups/tides.rst
@@ -1,36 +1,78 @@
 .. _dev_ocean_tides:
 
 tides
-=========
+=====
 
 The ``tides`` test group implements single layer
 barotropic, tidal cases as described in :ref:`ocean_tides` in
 the User's Guide.
 
+mesh test case
+--------------
+This test case generates the bathymetric dataset and horizontal mesh for tidal simulations.
+A :py:class:`compass.ocean.tests.tides.mesh.Mesh` object is created with ``mesh_name`` as one of its
+arguments. Based on this argument, it determines the appropriate child class of 
+:py:class:`compass.mesh.spherical.SphericalBaseStep` to create the base mesh and adds a
+:py:class:`compass.ocean.mesh.cull.CullMeshStep`. Prior to creation of the base mesh, a
+:py:class:`compass.ocean.tests.tides.dem.CreatePixelFile` step is added to create the 
+"pixel" files used to remap bathymetry data onto the mesh in the initial state test case.
+
+meshes
+~~~~~~
+``tides`` currently defines 2 meshes, with more to come.
+
+Icos7
+^^^^^
+This is a uniform resolution mesh based on icosahedral subdivision. It has approximately 60 km
+resolution globally. It is defined by :py:class:`compass.mesh.spherical.IcosahedralMeshStep`.
+
+vr45to5
+^^^^^^^
+This is a variable resolution mesh which has refinement based on bathymetric depth and slope
+criteria. It has a maximum resolution of 45 km and a minimum resolution of 5 km along coastlines
+and steep bathymetric gradients. It is defined by :py:class:`compass.ocean.tests.tides.mesh.vr45to5.VRTidesMesh`
+which inherits from :py:class:`compass.mesh.spherical.QuasiUniformSphericalMeshStep` and overrides
+the ``build_cell_width_lat_lon`` method with the resolution specification previously described.
+
+
+initial state test case
+-----------------------
+This test case contains steps to calculate parameters necessary for the 
+wave drag schemes in MPAS-Ocean, remaps the bathymetry onto the mesh, and
+generates the initial state. These steps are added in the construction of 
+the :py:class:`compass.ocean.tests.tides.init.Init` object. 
+
 initial_state
--------------
+~~~~~~~~~~~~~
 The class :py:class:`compass.ocean.tests.tides.init.initial_state.InitialState`
 defines a step for running MPAS-Ocean in init mode. The vertical mesh is
 set up with a single layer. 
 
-interpolate_wave_drag
----------------------
-The class :py:class:`compass.ocean.tests.tides.init.interpolate_wave_drag.InterpolateWaveDrag`
-defines a step for interpolating HYCOM data onto the MPAS-O mesh
-for the topographic wave drag scheme.
+calculate_wave_drag
+~~~~~~~~~~~~~~~~~~~
+The class :py:class:`compass.ocean.tests.tides.init.calculate_wave_drag.CalculateWaveDrag`
+defines a step for calculating bathymetric slopes and interpolating buoyancy frequency data onto
+the MPAS-O mesh for the topographic wave drag parameterizations.
 
 remap_bathymetry
-----------------
+~~~~~~~~~~~~~~~~
 The class :py:class:`compass.ocean.tests.tides.init.remap_bathymetry.RemapBathymetry`
-defines a step to perform an integral remap of bathyetry data onto the MPAS-O mesh.
+defines a step to perform an integral remap of bathymetry data onto the MPAS-O mesh.
+
+forward test case
+-----------------
+The forward test case contains steps to run the forward simulation and 
+compare the harmonic constituents with the TPXO database. These steps
+are added in the creation of the :py:class:`compass.ocean.tests.tides.forward.Forward`
+object.
 
 forward
--------
+~~~~~~~
 The class :py:class:`compass.ocean.tests.tides.forward.forward.ForwardStep`
 defines a step to run MPAS-Ocean in forward mode.
 
 analysis
---------
+~~~~~~~~
 The class :py:class:`compass.ocean.tests.tides.analysis.Analysis`
 defines a step to extract harmonic constituent data from the TPXO database.
 These values are used to compute and plot errors.
diff --git a/docs/users_guide/ocean/test_groups/images/M2_plot.png b/docs/users_guide/ocean/test_groups/images/M2_plot_icos7.png
similarity index 100%
rename from docs/users_guide/ocean/test_groups/images/M2_plot.png
rename to docs/users_guide/ocean/test_groups/images/M2_plot_icos7.png
diff --git a/docs/users_guide/ocean/test_groups/images/M2_plot_vr45to5.png b/docs/users_guide/ocean/test_groups/images/M2_plot_vr45to5.png
new file mode 100644
index 0000000000..35c86e18ec
Binary files /dev/null and b/docs/users_guide/ocean/test_groups/images/M2_plot_vr45to5.png differ
diff --git a/docs/users_guide/ocean/test_groups/tides.rst b/docs/users_guide/ocean/test_groups/tides.rst
index c1389b293f..247018bab5 100644
--- a/docs/users_guide/ocean/test_groups/tides.rst
+++ b/docs/users_guide/ocean/test_groups/tides.rst
@@ -8,11 +8,15 @@ initial conditions, forward simulations, and validation for global,
 realistic ocean domains. These are single layer, barotropic simulations
 forced with tidal potential. Self-attraction and loading effects are included
 as well as a parameterization for topographic wave drag.
+Currently the `Zarron and Egbert (2006) <https://doi.org/10.1175/JPO2878.1>`_
+wave drag parameterization is the default option. But MPAS-Ocean also supports
+the `local generation formula <https://doi.org/10.5194/gmd-14-1125-2021>`_ and
+`Jayne and St.Laurent (2001) <https://doi.org/10.1029/2000GL012044>`_.
 Wetting and drying is turned on in these configurations to prevent minimum
 depth issues.
 The tidal solution is decomposed into harmonic constituents during the 
 simulation and are compared with the TPXO database.
-Currently, the icosaheral 7 mesh is supported with
+Currently, the icosaheral 7 and vr45to5 meshes are supported with
 more mesh resolutions to be added in the future.
 
 Shared config options
@@ -87,6 +91,20 @@ mesh test case
 The mesh test case produces the horizontal mesh. The base mesh has global coverage
 and is culled to remove land cells. Cells beneath ice shelves are retained in the mesh.
 
+Icos7
+^^^^^
+This mesh is an Icosahedral mesh with 7 refienemnt steps, resulting in a globally uniform
+resolution of about 60km. This mesh is meant for efficient testing of tidal physics and
+does not produce very accurate global tides (Deep RMSE M2 ~11cm).
+
+VR45to5
+^^^^^^^
+This is a variable resolution mesh that ranges between 45km to 5km
+at the coasts. It uses refiniment criteria based on depth and bathymetric slope.
+More details can be found in `Barton et al. (2022) <https://doi.org/10.1029/2022MS003207>`_.
+This mesh achieves a more accurate tidal result competitive with other non-data assimilative
+models (Deep RMSE M2 ~3.3cm).
+
 .. _tides_init:
 
 init test case
@@ -98,17 +116,21 @@ remap bathymetry step
 ^^^^^^^^^^^^^^^^^^^^^
 This step performs an integral remap of bathymetric data onto the MPAS-O mesh
 
-interpolate wave drag step
+calculate wave drag step
 ^^^^^^^^^^^^^^^^^^^^^^^^^^
-In this step, a timescale factor from HYCOM data is interpolated onto
+In this step, a several values are calculated and interpolated onto
 the MPAS-O mesh to be used in the topographic wave drag parameterization.
-A plot is also produced to verify the interpolation.
+The climatological depth-averaged and bottom bouancy frequency values are
+interpolated from WOA are interpolated onto the MPAS mesh. The bathymetric
+gradients are also computed along with the standard deviaion of the subgrid
+bathymetry.
 
 initial state step
 ^^^^^^^^^^^^^^^^^^
 The initial state step runs MPAS-Ocean in init mode to create the initial
 condition file for the forward run. The vertical mesh is setup for a
-single layer.
+single layer. This step also computes the adjusted sea surface height,
+accounting for land ice pressure.
 
 .. _tides_forward:
 
@@ -132,6 +154,14 @@ global, shallow, and deep RMS errors. The global MPAS-O and TPXO solutions are
 plotted along with the spatial error fields. Below is an example of the type 
 of plots produced for the 5 major constituents.
 
-.. image:: images/M2_plot.png
+For the Icos7 mesh:
+
+.. image:: images/M2_plot_icos7.png
+   :width: 500px
+   :align: center
+
+And for the VR45to5 case:
+
+.. image:: images/M2_plot_vr45to5.png
    :width: 500px
    :align: center