diff --git a/compass/ocean/__init__.py b/compass/ocean/__init__.py
index ea766c088b..605db80b34 100644
--- a/compass/ocean/__init__.py
+++ b/compass/ocean/__init__.py
@@ -15,6 +15,7 @@
 from compass.ocean.tests.sphere_transport import SphereTransport
 from compass.ocean.tests.spherical_harmonic_transform import \
     SphericalHarmonicTransform
+from compass.ocean.tests.tides import Tides
 from compass.ocean.tests.ziso import Ziso
 
 
@@ -44,4 +45,5 @@ def __init__(self):
         self.add_test_group(Soma(mpas_core=self))
         self.add_test_group(SphereTransport(mpas_core=self))
         self.add_test_group(SphericalHarmonicTransform(mpas_core=self))
+        self.add_test_group(Tides(mpas_core=self))
         self.add_test_group(Ziso(mpas_core=self))
diff --git a/compass/ocean/tests/tides/__init__.py b/compass/ocean/tests/tides/__init__.py
new file mode 100644
index 0000000000..ec029d4db5
--- /dev/null
+++ b/compass/ocean/tests/tides/__init__.py
@@ -0,0 +1,30 @@
+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
+
+
+class Tides(TestGroup):
+    """
+    A test group for tidal simulations with MPAS-Ocean
+    """
+    def __init__(self, mpas_core):
+        """
+        mpas_core : compass.ocean.Ocean
+            the MPAS core that this test group belongs to
+        """
+        super().__init__(mpas_core=mpas_core,
+                         name='tides')
+
+        for mesh_name in ['Icos7']:
+
+            mesh = Mesh(test_group=self, mesh_name=mesh_name)
+            self.add_test_case(mesh)
+
+            init = Init(test_group=self, mesh=mesh)
+            self.add_test_case(init)
+
+            self.add_test_case(Forward(test_group=self,
+                                       mesh=mesh,
+                                       init=init))
diff --git a/compass/ocean/tests/tides/analysis/__init__.py b/compass/ocean/tests/tides/analysis/__init__.py
new file mode 100644
index 0000000000..789e85dbc7
--- /dev/null
+++ b/compass/ocean/tests/tides/analysis/__init__.py
@@ -0,0 +1,452 @@
+from compass.step import Step
+
+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
+from mpas_tools.logging import check_call
+
+
+class Analysis(Step):
+    """
+    A step for producing harmonic constituent errors and validation plots
+
+    Attributes
+    ----------
+    harmonic_analysis_file : str
+        File containing MPAS-O constitents
+
+    grid_file : str
+        Name of file containing MPAS-O mesh information
+
+    constituents : list
+        List of constituents to extract from TPXO database
+
+    tpxo_version : str
+        Version of TPXO to use in validation
+    """
+    def __init__(self, test_case):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        test_case : compass.ocean.tests.tides.forward.Forward
+            The test case this step belongs to
+        """
+        super().__init__(test_case=test_case, name='analysis')
+
+        self.harmonic_analysis_file = 'harmonicAnalysis.nc'
+        self.grid_file = 'initial_state.nc'
+        self.constituents = ['k1', 'k2', 'm2', 'n2', 'o1', 'p1', 'q1', 's2']
+
+        self.add_input_file(
+            filename=self.harmonic_analysis_file,
+            target='../forward/analysis_members/harmonicAnalysis.nc')
+
+        self.add_input_file(
+            filename=self.grid_file,
+            target='../forward/initial_state.nc')
+
+    def setup(self):
+        """
+        Setup test case and download data
+        """
+
+        config = self.config
+        self.tpxo_version = config.get('tides', 'tpxo_version')
+
+        os.makedirs(f'{self.work_dir}/TPXO_data', exist_ok=True)
+        if self.tpxo_version == 'TPXO9':
+            for constituent in self.constituents:
+                self.add_input_file(
+                    filename=f'TPXO_data/h_{constituent}_tpxo',
+                    target=f'TPXO9/h_{constituent}_tpxo9_atlas_30_v5',
+                    database='tides')
+                self.add_input_file(
+                    filename=f'TPXO_data/u_{constituent}_tpxo',
+                    target=f'TPXO9/u_{constituent}_tpxo9_atlas_30_v5',
+                    database='tides')
+            self.add_input_file(
+                filename='TPXO_data/grid_tpxo',
+                target='TPXO9/grid_tpxo9_atlas_30_v5',
+                database='tides')
+        elif self.tpxo_version == 'TPXO8':
+            for constituent in self.constituents:
+                self.add_input_file(
+                    filename=f'TPXO_data/h_{constituent}_tpxo',
+                    target=f'TPXO8/hf.{constituent}_tpxo8_atlas_30c_v1.out',
+                    database='tides')
+                self.add_input_file(
+                    filename=f'TPXO_data/u_{constituent}_tpxo',
+                    target=f'TPXO8/uv.{constituent}_tpxo8_atlas_30c_v1.out',
+                    database='tides')
+            self.add_input_file(
+                filename='TPXO_data/grid_tpxo',
+                target='TPXO8/grid_tpxo8_atlas_30_v1',
+                database='tides')
+
+    def write_coordinate_file(self, idx):
+        """
+        Write mesh coordinates for TPXO extraction
+        """
+
+        # Read in mesh
+        grid_nc = netCDF4.Dataset(self.grid_file, 'r')
+        lon_grid = np.degrees(grid_nc.variables['lonCell'][idx])
+        lat_grid = np.degrees(grid_nc.variables['latCell'][idx])
+        nCells = len(lon_grid)
+
+        # 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.close()
+
+    def setup_otps2(self):
+        """
+        Write input files for TPXO extraction
+        """
+
+        for con in self.constituents:
+            print(f'setup {con}')
+
+            # Lines for the setup_con files
+            lines = [{'inp': f'inputs/Model_atlas_{con}',
+                      'comment': '! 1. tidal model control file'},
+                     {'inp': 'lat_lon',
+                      'comment': '! 2. latitude/longitude/<time> file'},
+                     {'inp': 'z',
+                      'comment': '! 3. z/U/V/u/v'},
+                     {'inp':  con,
+                      'comment': '! 4. tidal constituents to include'},
+                     {'inp': 'AP',
+                      'comment': '! 5. AP/RI'},
+                     {'inp': 'oce',
+                      'comment': '! 6. oce/geo'},
+                     {'inp': '1',
+                      'comment': '! 7. 1/0 correct for minor constituents'},
+                     {'inp': f'outputs/{con}.out',
+                      'comment': '! 8. output file (ASCII)'}]
+
+            # Create directory for setup_con and Model_atlas_con files
+            if not os.path.exists('inputs'):
+                os.mkdir('inputs')
+
+            # Write the setup_con file
+            f = open('inputs/'+con+'_setup', 'w')
+            for line in lines:
+                spaces = 28 - len(line['inp'])
+                f.write(line['inp'] + spaces*' ' + line['comment'] + '\n')
+            f.close()
+
+            # Write the Model_atlas_con file
+            f = open(f'inputs/Model_atlas_{con}', 'w')
+            f.write(f'TPXO_data/h_{con}_tpxo\n')
+            f.write(f'TPXO_data/u_{con}_tpxo\n')
+            f.write('TPXO_data/grid_tpxo')
+            f.close()
+
+            # Create directory for the con.out files
+            if not os.path.exists('outputs'):
+                os.mkdir('outputs')
+
+    def run_otps2(self):
+        """
+        Perform TPXO extraction
+        """
+
+        # Run the executable
+        for con in self.constituents:
+            print('')
+            print(f'run {con}')
+            check_call(f'extract_HC < inputs/{con}_setup',
+                       logger=self.logger, shell=True)
+
+    def read_otps2_output(self, idx):
+        """
+        Read TPXO extraction output
+        """
+
+        start = idx[0]
+        for con in self.constituents:
+
+            f = open(f'outputs/{con}.out', 'r')
+            lines = f.read().splitlines()
+            for i, line in enumerate(lines[3:]):
+                line_sp = line.split()
+                if line_sp[2] != '*************':
+                    val = float(line_sp[2])
+                    self.mesh_AP[con]['amp'][start+i] = val
+                else:
+                    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
+
+                else:
+                    self.mesh_AP[con]['phase'][start+i] = -9999
+
+    def append_tpxo_data(self):
+        """
+        Inject TPXO data into harmonic analysis file
+        """
+
+        data_nc = netCDF4.Dataset(self.harmonic_analysis_file, 'a',
+                                  format='NETCDF3_64BIT_OFFSET')
+        for con in self.constituents:
+
+            # Inject amplitude
+            amp_varname = f'{con.upper()}Amplitude{self.tpxo_version}'
+            amp_var = data_nc.createVariable(
+                amp_varname,
+                np.float64,
+                ('nCells'))
+            amp_var.units = 'm'
+            amp_var.long_name = f'Amplitude of {con.upper()} tidal ' \
+                'consitiuent at each cell center from the ' \
+                f'{self.tpxo_version} model'
+            amp_var[:] = self.mesh_AP[con]['amp'][:]
+
+            # Inject phase
+            phase_varname = f'{con.upper()}Phase{self.tpxo_version}'
+            phase_var = data_nc.createVariable(
+                phase_varname,
+                np.float64,
+                ('nCells'))
+            phase_var.units = 'deg'
+            phase_var.long_name = f'Phase of {con.upper()} tidal ' \
+                'consitiuent at each cell center from the ' \
+                f'{self.tpxo_version} model'
+            phase_var[:] = self.mesh_AP[con]['phase'][:]
+
+        data_nc.close()
+
+    def check_tpxo_data(self):
+        """
+        Check if TPXO data exists in harmonic analysis file
+        """
+
+        data_nc = netCDF4.Dataset(self.harmonic_analysis_file, 'r',
+                                  format='NETCDF3_64BIT_OFFSET')
+        self.nCells = len(data_nc.dimensions['nCells'])
+
+        for con in self.constituents[:]:
+            amp_var = f'{con.upper()}Amplitude{self.tpxo_version}'
+            phase_var = f'{con.upper()}Phase{self.tpxo_version}'
+            if (amp_var in data_nc.variables) \
+                    and (phase_var in data_nc.variables):
+
+                self.constituents.remove(con)
+                print(f'{con} TPXO Constituent already exists '
+                      f'in {self.harmonic_analysis_file}')
+
+        data_nc.close()
+
+    def plot(self):
+        """
+        Calculate errors and plot consitituents
+        """
+
+        plt.switch_backend('agg')
+        cmap_reversed = cm.get_cmap('Spectral_r')
+
+        # Initialize plotting variables
+        TW = 2                         # Tick width
+        TL = 2                         # Tick length
+        TF = 8                         # Tick label size
+
+        # Open data file
+        data_file = self.harmonic_analysis_file
+        data_nc = netCDF4.Dataset(data_file, 'r')
+
+        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
+
+        nCells = lon_grid.size
+        data1 = np.zeros((nCells))
+        data2 = np.zeros((nCells))
+        data1_phase = np.zeros((nCells))
+        data2_phase = np.zeros((nCells))
+        depth = np.zeros((nCells))
+        area = np.zeros((nCells))
+
+        depth[:] = data_nc.variables['bottomDepth'][:]
+        area[:] = data_nc.variables['areaCell'][:]
+
+        constituent_list = ['K1', 'M2', 'N2', 'O1', 'S2']
+
+        # 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, 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.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)]]
+
+        for i, con in enumerate(constituent_list):
+
+            print('')
+            print(f' ====== {con} Constituent ======')
+
+            # Get data
+            data1[:] = data_nc.variables[
+                f'{con}Amplitude'][:]
+            data1_phase[:] = data_nc.variables[
+                f'{con}Phase'][:]
+            data2[:] = data_nc.variables[
+                f'{con}Amplitude{self.tpxo_version}'][:]
+            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
+
+            # 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)
+
+            # Calculate mean (global) values
+            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)
+            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))
+            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)
+            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))
+            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)
+            print('Deep RMSE (Amp) = ', deep_rmse_amp)
+            print('Deep RMSE (Com) = ', deep_rmse_com)
+
+            rmse_amp = rmse_amp**0.5
+            rmse_com = rmse_com**0.5
+
+            # Plot data
+            fig = plt.figure(figsize=(18, 12))
+            subplot_title = [f'{con} Amplitude (simulation) [m]',
+                             f'{con} Amplitude (TPXO8) [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,
+                                     projection=ccrs.PlateCarree())
+                ax.set_title(subplot_title[subplot], fontsize=20)
+                levels = subplot_ticks[i][subplot][:]
+
+                # MPAS amplitude and phase
+                if subplot == 0:
+                    cf = ax.tricontourf(lon_grid, lat_grid, data1,
+                                        levels=levels,
+                                        transform=ccrs.PlateCarree(),
+                                        cmap=cmap_reversed)
+                    ax.tricontour(lon_grid, lat_grid, data1_phase,
+                                  levels=10, linewidths=0.5, colors='k')
+
+                # TPXO amplitude and phase
+                elif subplot == 1:
+                    ix = np.logical_and(data2_phase >= 0, data2_phase < 360)
+                    cf = ax.tricontourf(lon_grid, lat_grid, data2,
+                                        levels=levels,
+                                        transform=ccrs.PlateCarree(),
+                                        cmap=cmap_reversed)
+                    ax.tricontour(lon_grid[ix], lat_grid[ix], data2_phase[ix],
+                                  levels=10, linewidths=0.5, colors='k')
+
+                # Amplitude RMSE
+                elif subplot == 2:
+                    cf = ax.tricontourf(lon_grid, lat_grid, rmse_amp,
+                                        levels=levels,
+                                        transform=ccrs.PlateCarree(),
+                                        cmap='OrRd')
+
+                # Complex RMSE
+                elif subplot == 3:
+                    cf = ax.tricontourf(lon_grid, lat_grid, rmse_com,
+                                        levels=levels,
+                                        transform=ccrs.PlateCarree(),
+                                        cmap='OrRd')
+
+                ax.set_extent([-180, 180, -90, 90], crs=ccrs.PlateCarree())
+                ax.add_feature(cfeature.LAND, zorder=100)
+                ax.add_feature(cfeature.LAKES, alpha=0.5, zorder=101)
+                ax.add_feature(cfeature.COASTLINE, zorder=101)
+                ax.tick_params(axis='both', which='major',
+                               length=TL, width=TW, labelsize=TF)
+                cbar = fig.colorbar(cf, ax=ax,
+                                    ticks=levels.round(2), shrink=0.6)
+                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))
+            fig.suptitle(f'Complex RMSE: Global = {global_err} cm; '
+                         f'Deep = {deep_err} cm; '
+                         f'Shallow = {shallow_err} cm',
+                         fontsize=20)
+            plt.savefig(f'{con}_plot.png')
+            plt.close()
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+
+        # Check if TPXO values aleady exist in harmonic_analysis.nc
+        self.check_tpxo_data()
+
+        # Setup input files for TPXO extraction
+        self.setup_otps2()
+
+        # Setup chunking for TPXO extraction with large meshes
+        indices = np.arange(self.nCells)
+        nchunks = np.ceil(self.nCells/200000)
+        index_chunks = np.array_split(indices, nchunks)
+
+        # Initialize data structure for TPXO values
+        self.mesh_AP = {}
+        for con in self.constituents:
+            self.mesh_AP[con] = {'amp': np.zeros((self.nCells)),
+                                 'phase': np.zeros((self.nCells))}
+
+        # Extract TPXO values
+        for idx in index_chunks:
+            self.write_coordinate_file(idx)
+            self.run_otps2()
+            self.read_otps2_output(idx)
+
+        # Inject TPXO values
+        self.append_tpxo_data()
+
+        # Calulate and plot global errors
+        self.plot()
diff --git a/compass/ocean/tests/tides/configure.py b/compass/ocean/tests/tides/configure.py
new file mode 100644
index 0000000000..21c92db11d
--- /dev/null
+++ b/compass/ocean/tests/tides/configure.py
@@ -0,0 +1,26 @@
+from compass.ocean.tests.global_ocean.metadata import \
+    get_author_and_email_from_git
+
+
+def configure_tides(test_case, mesh):
+    """
+    Modify the configuration options for this test case
+
+    Parameters
+    ----------
+    test_case : compass.TestCase
+        The test case to configure
+
+    mesh : compass.ocean.tests.global_ocean.mesh.Mesh
+        The test case that produces the mesh for this run
+
+    init : compass.ocean.tests.global_ocean.init.Init, optional
+        The test case that produces the initial condition for this run
+    """
+    config = test_case.config
+
+    config.add_from_package('compass.mesh', 'mesh.cfg')
+    config.add_from_package(mesh.package, mesh.mesh_config_filename,
+                            exception=True)
+
+    get_author_and_email_from_git(config)
diff --git a/compass/ocean/tests/tides/forward/__init__.py b/compass/ocean/tests/tides/forward/__init__.py
new file mode 100644
index 0000000000..057e8e27bf
--- /dev/null
+++ b/compass/ocean/tests/tides/forward/__init__.py
@@ -0,0 +1,56 @@
+from compass.testcase import TestCase
+from compass.ocean.tests.tides.forward.forward import ForwardStep
+from compass.ocean.tests.tides.analysis import Analysis
+from compass.ocean.tests.tides.configure import configure_tides
+import os
+
+
+class Forward(TestCase):
+    """
+    A test case for performing a forward run for a tidal case
+
+    Attributes
+    ----------
+    mesh : compass.ocean.tests.global_ocean.mesh.Mesh
+        The test case that produces the mesh for this run
+
+    init : compass.ocean.tests.tides.init.Init
+        The test case that produces the initial condition for this run
+    """
+    def __init__(self, test_group, mesh, init):
+        """
+        Create test case
+
+        Parameters
+        ----------
+        test_group : compass.ocean.tests.tides.Hurricane
+            The test group that this test case belongs to
+
+        mesh : compass.ocean.tests.global_ocean.mesh.Mesh
+            The test case that produces the mesh for this run
+
+        init : compass.ocean.tests.tides.init.Init
+            The test case that produces the initial condition for this run
+        """
+        name = 'forward'
+        mesh_name = mesh.mesh_name
+        subdir = os.path.join(mesh_name, name)
+        super().__init__(test_group=test_group,
+                         subdir=subdir,
+                         name=name)
+        self.mesh = mesh
+
+        step = ForwardStep(test_case=self, mesh=mesh, init=init)
+
+        step.add_output_file(filename='output/output.nc')
+        step.add_output_file(filename='analysis_members/harmonicAnalysis.nc')
+        self.add_step(step)
+
+        step = Analysis(test_case=self)
+        self.add_step(step)
+
+    def configure(self):
+        """
+        Modify the configuration options for this test case
+        """
+        configure_tides(test_case=self, mesh=self.mesh)
diff --git a/compass/ocean/tests/tides/forward/forward.py b/compass/ocean/tests/tides/forward/forward.py
new file mode 100644
index 0000000000..f6ce890212
--- /dev/null
+++ b/compass/ocean/tests/tides/forward/forward.py
@@ -0,0 +1,83 @@
+from compass.model import run_model
+from compass.step import Step
+
+
+class ForwardStep(Step):
+    """
+    A step for performing forward MPAS-Ocean runs as part of tides test
+    cases.
+
+    Attributes
+    ----------
+    mesh : compass.ocean.tests.tides.mesh.Mesh
+        The test case that produces the mesh for this run
+
+    init : compass.ocean.tests.tides.init.Init
+        The test case that produces the initial condition for this run
+    """
+    def __init__(self, test_case, mesh, init, name='forward', subdir=None):
+        """
+        Create a new step
+
+        Parameters
+        ----------
+        test_case : compass.ocean.tests.tides.forward.Forward
+            The test case this step belongs to
+
+        mesh : compass.ocean.tests.global_ocean.mesh.Mesh
+            The test case that produces the mesh for this run
+
+        init : compass.ocean.tests.tides.init.Init
+            The test case that produces the initial condition for this run
+
+        name : str, optional
+            the name of the step
+
+        subdir : str, optional
+            the subdirectory for the step.  The default is ``name``
+        """
+        self.mesh = mesh
+        self.init = init
+        super().__init__(test_case=test_case, name=name)
+
+        self.add_namelist_file(
+            'compass.ocean.tests.tides.forward', 'namelist.ocean')
+        self.add_streams_file(
+            'compass.ocean.tests.tides.forward', 'streams.ocean')
+
+        mesh_package = mesh.package
+        self.add_namelist_file(mesh_package, 'namelist.ocean')
+
+        initial_state_path = f'{init.path}/initial_state'
+        interpolate_path = f'{init.path}/interpolate'
+        initial_state_target = f'{initial_state_path}/initial_state.nc'
+
+        self.add_input_file(
+            filename='initial_state.nc',
+            work_dir_target=initial_state_target)
+        self.add_input_file(
+            filename='forcing_data.nc',
+            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')
+        self.add_input_file(
+            filename='graph.info',
+            work_dir_target=f'{initial_state_path}/graph.info')
+
+        self.add_model_as_input()
+
+    def setup(self):
+        """
+        Set up the test case in the work directory, including downloading any
+        dependencies
+        """
+        self.ntasks = self.config.getint('tides', 'forward_ntasks')
+        self.min_tasks = self.config.getint('tides', 'forward_min_tasks')
+        self.threads = self.config.getint('tides', 'forward_threads')
+
+    def run(self):
+        """
+        Run this step of the testcase
+        """
+        run_model(self)
diff --git a/compass/ocean/tests/tides/forward/namelist.ocean b/compass/ocean/tests/tides/forward/namelist.ocean
new file mode 100644
index 0000000000..7144e04a64
--- /dev/null
+++ b/compass/ocean/tests/tides/forward/namelist.ocean
@@ -0,0 +1,47 @@
+config_start_time = '2012-10-01_00:00:00'
+config_run_duration = '00-00-125_00:00:00'
+config_calendar_type = 'noleap'
+config_pio_num_iotasks = 0
+config_pio_stride = 1
+config_time_integrator = 'RK4'
+config_hmix_scaleWithMesh = .true.
+config_hmix_use_ref_cell_width = .true.
+config_use_cvmix = .false.
+config_use_bulk_wind_stress = .true.
+config_use_self_attraction_loading = .true. 
+config_self_attraction_loading_compute_interval = '0000_00:30:00'
+config_use_parallel_self_attraction_loading = .true.
+config_parallel_self_attraction_loading_order = 50
+config_self_attraction_and_loading_beta = 0.0
+config_use_tidal_potential_forcing = .true.
+config_tidal_potential_reference_time = '2012-10-01_00:00:00'
+config_use_tidal_potential_forcing_K2 = .false.
+config_use_tidal_potential_forcing_Q1 = .false.
+config_use_tidal_potential_forcing_P1 = .false.
+config_tidal_potential_ramp = 15.0
+config_use_implicit_bottom_drag = .false.
+config_use_implicit_bottom_roughness = .true.
+config_use_topographic_wave_drag = .true.
+config_use_wetting_drying = .true.
+config_prevent_drying = .true
+config_drying_min_cell_height = 0.1
+config_zero_drying_velocity = .true.
+config_pressure_gradient_type = 'ssh_gradient'
+config_disable_thick_vadv = .true.
+config_disable_thick_sflux = .true.
+config_disable_vel_explicit_bottom_drag = .true.
+config_disable_vel_vadv = .true.
+config_disable_tr_all_tend = .true.
+config_AM_globalStats_enable = .true.
+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_start_delay = 30
+config_AM_harmonicAnalysis_duration = 90
+config_AM_harmonicAnalysis_use_K2 = .false.
+config_AM_harmonicAnalysis_use_Q1 = .false.
+config_AM_harmonicAnalysis_use_P1 = .false.
+
+
diff --git a/compass/ocean/tests/tides/forward/streams.ocean b/compass/ocean/tests/tides/forward/streams.ocean
new file mode 100755
index 0000000000..6a8bd3a27c
--- /dev/null
+++ b/compass/ocean/tests/tides/forward/streams.ocean
@@ -0,0 +1,151 @@
+<streams>
+
+<immutable_stream name="mesh"
+                  type="input"
+                  filename_template="initial_state.nc"
+                  input_interval="initial_only"/>
+
+<immutable_stream name="input"
+                  type="input"
+                  filename_template="initial_state.nc"
+                  input_interval="initial_only"/>
+
+<immutable_stream name="restart"
+                  type="input;output"
+                  filename_template="restarts/restart.$Y-$M-$D_$h.$m.$s.nc"
+                  filename_interval="output_interval"
+                  reference_time="0001-01-00_00:00:00"
+                  clobber_mode="truncate"
+                  input_interval="initial_only"
+                  output_interval="00-00-20_00:00:00"/>
+
+<stream name="output"
+        filename_template="output/output.nc"
+        filename_interval="01-00-00_00:00:00"
+        clobber_mode="truncate"
+        reference_time="0001-01-00_00:00:00"
+        type="output"
+        output_interval="00-00-10_00:00:00">
+
+    <stream name="mesh"/>
+    <var name="ssh"/>
+    <var name="ssh_sal"/>
+    <var name="layerThickness"/>
+    <var name="xtime"/>
+    <var name="kineticEnergyCell"/>
+    <var name="relativeVorticityCell"/>
+    <var name="areaCellGlobal"/>
+    <var name="areaEdgeGlobal"/>
+    <var name="areaTriangleGlobal"/>
+    <var name="volumeCellGlobal"/>
+    <var name="volumeEdgeGlobal"/>
+    <var name="CFLNumberGlobal"/>
+</stream>
+
+<stream name="block_debug_output"
+        filename_template="output_debug_block_$B.nc"
+        filename_interval="1000-00-00_00:00:00"
+        clobber_mode="truncate"
+        reference_time="0001-01-01_00:00:00"
+        type="output"
+        output_interval="1000-00-00_00:00:00">
+
+    <stream name="mesh"/>
+    <var_struct name="tracers"/>
+    <var name="xtime"/>
+    <var name="layerThickness"/>
+    <var name="normalVelocity"/>
+    <var name="ssh"/>
+</stream>
+
+<stream name="topographic_wave_drag"
+        type="input"
+        filename_template="topographic_wave_drag.nc"
+        input_interval="initial_only" >
+
+        <var name="topographic_wave_drag"/>
+</stream>
+
+<stream name="globalStatsOutput"
+        runtime_format="single_file"
+        filename_interval="01-00-00_00:00:00"
+        clobber_mode="truncate"
+        output_interval="0000_06:00:00"
+        reference_time="0001-01-01_00:00:00"
+        filename_template="analysis_members/globalStats.$Y-$M-$D_$h.$m.$s.nc"
+        packages="globalStatsAMPKG"
+	type="output">
+
+    <var_array name="minGlobalStats"/>
+    <var_array name="maxGlobalStats"/>
+    <var_array name="sumGlobalStats"/>
+    <var_array name="rmsGlobalStats"/>
+    <var_array name="avgGlobalStats"/>
+    <var_array name="vertSumMinGlobalStats"/>
+    <var_array name="vertSumMaxGlobalStats"/>
+    <var name="daysSinceStartOfSim"/>
+    <var name="xtime"/>
+    <var name="volumeCellGlobal"/>
+    <var name="CFLNumberGlobal"/>
+</stream>
+
+<stream name="forcing_data"
+	filename_template="forcing_data.nc"
+	input_interval="initial_only"
+	type="input">
+
+    <var_struct name="tracersSurfaceRestoringFields"/>
+    <var_struct name="tracersInteriorRestoringFields"/>
+    <var_struct name="tracersExponentialDecayFields"/>
+    <var_struct name="tracersIdealAgeFields"/>
+    <var_struct name="tracersTTDFields"/>
+    <var name="windStressZonal"/>
+    <var name="windStressMeridional"/>
+    <var name="landIceSurfaceTemperature"/>
+    <var name="atmosphericPressure"/>
+    <var name="latentHeatFlux"/>
+    <var name="sensibleHeatFlux"/>
+    <var name="shortWaveHeatFlux"/>
+    <var name="evaporationFlux"/>
+    <var name="rainFlux"/>
+</stream>
+
+<stream name="harmonicAnalysisOutput"
+        filename_template="analysis_members/harmonicAnalysis.nc"
+        type="output"
+        mode="forward;analysis"
+        output_interval="1000-00-00_00:00:00"
+        packages="harmonicAnalysisAMPKG"
+        clobber_mode="truncate"
+        runtime_format="single_file">
+
+    <stream name="mesh"/>
+    <var name="M2Amplitude"/>
+    <var name="M2Phase"/>
+    <var name="S2Amplitude"/>
+    <var name="S2Phase"/>
+    <var name="N2Amplitude"/>
+    <var name="N2Phase"/>
+    <var name="K1Amplitude"/>
+    <var name="K1Phase"/>
+    <var name="O1Amplitude"/>
+    <var name="O1Phase"/>
+</stream>
+
+<stream name="harmonicAnalysisRestart"
+        filename_template="restarts/mpaso.rst.am.harmonicAnalysisRestart.$Y-$M-$D_$h.$m.$s.nc"
+        filename_interval="output_interval"
+        type="input;output"
+        mode="forward;analysis"
+        output_interval="stream:restart:output_interval"
+        input_interval="initial_only"
+        packages="harmonicAnalysisAMPKG"
+        clobber_mode="truncate"
+        runtime_format="single_file"
+        reference_time="0001-01-01_00:00:00">
+
+    <var name="leastSquaresLHSMatrix"/>
+    <var name="leastSquaresRHSVector"/>
+</stream>
+
+</streams>
diff --git a/compass/ocean/tests/tides/init/__init__.py b/compass/ocean/tests/tides/init/__init__.py
new file mode 100644
index 0000000000..be94f9b622
--- /dev/null
+++ b/compass/ocean/tests/tides/init/__init__.py
@@ -0,0 +1,48 @@
+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
+
+
+class Init(TestCase):
+    """
+    A test case for creating initial conditions on a global MPAS-Ocean mesh
+
+    Attributes
+    ----------
+    mesh : compass.ocean.tests.tides.mesh.Mesh
+        The test case that creates the mesh used by this test case
+    """
+    def __init__(self, test_group, mesh):
+        """
+        Create the test case
+
+        Parameters
+        ----------
+        test_group : compass.ocean.tests.tides.Tides
+            The tides test group that this test case belongs to
+
+        mesh : compass.ocean.tests.tides.mesh.Mesh
+            The test case that creates the mesh used by this test case
+        """
+        name = 'init'
+        mesh_name = mesh.mesh_name
+        subdir = os.path.join(mesh_name, name)
+        super().__init__(test_group=test_group, name=name, subdir=subdir)
+
+        self.mesh = mesh
+
+        self.add_step(InterpolateWaveDrag(test_case=self, mesh=mesh))
+        self.add_step(RemapBathymetry(test_case=self, mesh=mesh))
+        self.add_step(InitialState(test_case=self, mesh=mesh))
+
+    def configure(self):
+        """
+        Modify the configuration options for this test case
+        """
+        configure_tides(test_case=self, mesh=self.mesh)
diff --git a/compass/ocean/tests/tides/init/dem_pixel.py b/compass/ocean/tests/tides/init/dem_pixel.py
new file mode 100644
index 0000000000..e2ac404da3
--- /dev/null
+++ b/compass/ocean/tests/tides/init/dem_pixel.py
@@ -0,0 +1,859 @@
+
+import os
+import numpy as np
+import netCDF4 as nc
+from scipy.ndimage import gaussian_filter
+import argparse
+
+# Authors: Darren Engwirda
+
+
+def blend_front(e1st, i1st, e2nd, halo, sdev):
+    """
+    Create a mask of linear weights to 'blend' two elev. fun
+    at the ice sheet/shelf front:
+
+        ELEV = (1.-MASK) * E1ST + (0.+MASK) * E2ND
+
+    Elev. data assoc. with the 1st data-set is preserved in
+    pixels where ice-thickness is non-zero. The two datasets
+    are blended over a distance of approx. HALO pixels. The
+    mask is additionally smoothed via a Gaussian filter with
+    standard-deviation sigma=SDEV.
+
+    """
+
+    mask = np.full(
+        e1st.shape, halo + 1, dtype=np.float32)
+
+    nidx = np.full(
+        e1st.shape[0], False, dtype=np.bool)
+    sidx = np.full(
+        e1st.shape[0], False, dtype=np.bool)
+
+    bnds = np.asarray(np.round(
+        np.linspace(0, e1st.shape[0], 5)), dtype=np.uint32)
+
+    sidx[bnds[0]:bnds[1]:] = True
+    nidx[bnds[3]:bnds[4]:] = True
+
+    print("* blending SOUTH.")
+
+    part = mask[sidx, :]
+
+    npos = np.arange(+1, part.shape[0] + 1)
+    epos = np.arange(-1, part.shape[1] - 1)
+    spos = np.arange(-1, part.shape[0] - 1)
+    wpos = np.arange(+1, part.shape[1] + 1)
+
+    Y = part.shape[0]
+    X = part.shape[1]
+
+    npos[npos >= +Y] = Y - 1
+    spos[spos <= -1] = 0
+    epos[epos <= -1] = X - 1
+    wpos[wpos >= +X] = 0
+
+    part[i1st[sidx, :] > 0] = 0.
+
+    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.)
+
+    mask[sidx, :] = part
+
+    print("* blending NORTH.")
+
+    part = mask[nidx, :]
+
+    npos = np.arange(+1, part.shape[0] + 1)
+    epos = np.arange(-1, part.shape[1] - 1)
+    spos = np.arange(-1, part.shape[0] - 1)
+    wpos = np.arange(+1, part.shape[1] + 1)
+
+    Y = part.shape[0]
+    X = part.shape[1]
+
+    npos[npos >= +Y] = Y - 1
+    spos[spos <= -1] = 0
+    epos[epos <= -1] = X - 1
+    wpos[wpos >= +X] = 0
+
+    part[i1st[nidx, :] > 0] = 0.
+
+    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.)
+
+    mask[nidx, :] = part
+
+    print("* blending final.")
+
+    mask /= float(halo + 1.00)
+
+    mask = gaussian_filter(mask, sigma=sdev, mode="wrap")
+    mask = mask ** 1.50
+    mask[i1st >= 1] = 0.
+
+    return mask
+
+
+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.
+
+    """
+
+    print("Making RTopo-2.0.4 (60 arc-sec) pixel...")
+
+    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
+
+    xpos = np.asarray(data["lon"][:], dtype=np.float64)
+    ypos = np.asarray(data["lat"][:], dtype=np.float64)
+
+    elev = np.asarray(
+        data["bedrock_topography"][:], dtype=np.float32)
+    surf = np.asarray(
+        data["surface_elevation"][:], dtype=np.float32)
+    base = np.asarray(
+        data["ice_base_topography"][:], dtype=np.float32)
+
+    elev = (elev[:-1:, :-1:] + elev[+1::, :-1:] +
+            elev[:-1:, +1::] + elev[+1::, +1::]) / 4.
+
+    surf = (surf[:-1:, :-1:] + surf[+1::, :-1:] +
+            surf[:-1:, +1::] + surf[+1::, +1::]) / 4.
+
+    base = (base[:-1:, :-1:] + base[+1::, :-1:] +
+            base[:-1:, +1::] + base[+1::, +1::]) / 4.
+
+    elev = np.asarray(np.round(elev), dtype=np.int16)
+    surf = np.asarray(np.round(surf), dtype=np.int16)
+    base = np.asarray(np.round(base), dtype=np.int16)
+
+    iceh = surf - base
+    iceh[base == 0] = 0
+
+    ocnh = np.maximum(0, base - elev)
+
+    root = nc.Dataset(
+        os.path.join(
+            save_path, "RTopo_2_0_4_60sec_pixel.nc"),
+        "w", format="NETCDF4")
+    root.description = "A zipped RTopo-2.0.4 (60 arc-sec) " + \
+        "data-set, pixel centred and compressed to int16_t."
+    root.source = "RTopo-2.0.4_1min_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 rtopo_30sec(elev_path, save_path):
+    """
+    Create a zipped and pixel centred version of RTopo 2.0.4
+    (30 arc-sec) to support remapping of elevation data.
+
+    """
+
+    print("Making RTopo-2.0.4 (30 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)  # quite valid_min/max
+
+    xpos = np.asarray(data["lon"][:], dtype=np.float64)
+    ypos = np.asarray(data["lat"][:], dtype=np.float64)
+
+    elev = np.asarray(
+        data["bedrock_topography"][:], dtype=np.float32)
+
+    data = nc.Dataset(os.path.join(
+        elev_path,
+        "RTopo-2.0.4_30sec_surface_elevation.nc"), "r")
+
+    data.set_auto_maskandscale(False)  # quite valid_min/max
+
+    surf = np.asarray(
+        data["surface_elevation"][:], dtype=np.float32)
+
+    data = nc.Dataset(os.path.join(
+        elev_path,
+        "RTopo-2.0.4_30sec_ice_base_topography.nc"), "r")
+
+    data.set_auto_maskandscale(False)  # quite valid_min/max
+
+    base = np.asarray(
+        data["ice_base_topography"][:], dtype=np.float32)
+
+    elev = (elev[:-1:, :-1:] + elev[+1::, :-1:] +
+            elev[:-1:, +1::] + elev[+1::, +1::]) / 4.
+
+    surf = (surf[:-1:, :-1:] + surf[+1::, :-1:] +
+            surf[:-1:, +1::] + surf[+1::, +1::]) / 4.
+
+    base = (base[:-1:, :-1:] + base[+1::, :-1:] +
+            base[:-1:, +1::] + base[+1::, +1::]) / 4.
+
+    elev = np.asarray(np.round(elev), dtype=np.int16)
+    surf = np.asarray(np.round(surf), dtype=np.int16)
+    base = np.asarray(np.round(base), dtype=np.int16)
+
+    iceh = surf - base
+    iceh[base == 0] = 0
+
+    ocnh = np.maximum(0, base - elev)
+
+    root = nc.Dataset(
+        os.path.join(
+            save_path, "RTopo_2_0_4_30sec_pixel.nc"),
+        "w", format="NETCDF4")
+    root.description = "A zipped RTopo-2.0.4 (30 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 srtmp_60sec(elev_path, save_path):
+    """
+    Create a zipped and pixel centred version of SRTM15+V2.1
+    (15 arc-sec) at 60 arc-sec.
+
+    """
+
+    print("Making SRTM15+V2.1 (60 arc-sec) pixel...")
+
+    data = nc.Dataset(os.path.join(
+        elev_path, "SRTM15+V2.1.nc"), "r")
+
+    elev = np.asarray(data["z"][:], 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, "SRTM15+V2.1_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 srtmp_30sec(elev_path, save_path):
+    """
+    Create a zipped and pixel centred version of SRTM15+V2.1
+    (15 arc-sec) at 30 arc-sec.
+
+    """
+
+    print("Making SRTM15+V2.1 (30 arc-sec) pixel...")
+
+    data = nc.Dataset(os.path.join(
+        elev_path, "SRTM15+V2.1.nc"), "r")
+
+    elev = np.asarray(data["z"][:], 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, "SRTM15+V2.1_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
+    data = root.createVariable(
+        "top_elevation", "i2", ("num_row", "num_col"))
+    data.units = "m"
+    data[:, :] = elev
+
+    root.close()
+
+
+def rtopo_srtmp_60sec(elev_path, save_path):
+    """
+    Create a zipped and pixel centred 'blend' of RTopo 2.0.4
+    and SRTM15+V2.1 at 60 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)
+
+    data = nc.Dataset(os.path.join(
+        elev_path, "SRTM15+V2.1_60sec_pixel.nc"), "r")
+
+    e2nd = np.asarray(
+        data["top_elevation"][:], dtype=np.int16)
+
+    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])
+
+    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_SRTM15+V2_1_60sec_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)
+    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 rtopo_srtmp_30sec(elev_path, save_path):
+    """
+    Create a zipped and pixel centred 'blend' of RTopo 2.0.4
+    and SRTM15+V2.1 at 30 arc-sec.
+
+    """
+
+    print("Making RTopo-SRTM+ (30 arc-sec) blend...")
+
+    data = nc.Dataset(os.path.join(
+        elev_path, "RTopo_2_0_4_30sec_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)
+
+    data = nc.Dataset(os.path.join(
+        elev_path, "SRTM15+V2.1_30sec_pixel.nc"), "r")
+
+    e2nd = np.asarray(
+        data["top_elevation"][:], dtype=np.int16)
+
+    mask = blend_front(e1st, i1st, e2nd, halo=40, sdev=4.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])
+
+    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_SRTM15+V2_1_30sec_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 " + \
+        "compressed to int16_t. RTopo data used under ice " + \
+        "sheets/shelves."
+    root.source = \
+        "RTopo-2.0.4_30sec_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)
+    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 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
+    data = root.createVariable(
+        "top_elevation", "i2", ("num_row", "num_col"))
+    data.units = "m"
+    data[:, :] = elev
+
+    root.close()
+
+
+def rtopo_gebco_60sec(elev_path, save_path):
+    """
+    Create a zipped and pixel centred 'blend' of RTopo 2.0.4
+    and GEBCO[2020] at 60 arc-sec.
+
+    """
+
+    print("Making RTopo-GEBCO (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)
+
+    data = nc.Dataset(os.path.join(
+        elev_path, "GEBCO_v2020_60sec_pixel.nc"), "r")
+
+    e2nd = np.asarray(
+        data["top_elevation"][:], dtype=np.int16)
+
+    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])
+
+    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_GEBCO_v2020_60sec_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 " + \
+        "compressed to int16_t. RTopo data used under ice " + \
+        "sheets/shelves. Remapped to 60 arc-sec."
+    root.source = \
+        "RTopo-2.0.4_1min_data.nc and GEBCO_2020.nc"
+    root.references = \
+        "doi.pangaea.de/10.1594/PANGAEA.905295 and " + \
+        "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)
+    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 rtopo_gebco_30sec(elev_path, save_path):
+    """
+    Create a zipped and pixel centred 'blend' of RTopo 2.0.4
+    and GEBCO[2020] at 30 arc-sec.
+
+    """
+
+    print("Making RTopo-GEBCO (30 arc-sec) blend...")
+
+    data = nc.Dataset(os.path.join(
+        elev_path, "RTopo_2_0_4_30sec_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)
+
+    data = nc.Dataset(os.path.join(
+        elev_path, "GEBCO_v2020_30sec_pixel.nc"), "r")
+
+    e2nd = np.asarray(
+        data["top_elevation"][:], dtype=np.int16)
+
+    mask = blend_front(e1st, i1st, e2nd, halo=40, sdev=4.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])
+
+    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_GEBCO_v2020_30sec_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 " + \
+        "compressed to int16_t. RTopo data used under ice " + \
+        "sheets/shelves. Remapped to 30 arc-sec."
+    root.source = \
+        "RTopo-2.0.4_30sec_data.nc and GEBCO_2020.nc"
+    root.references = \
+        "doi.pangaea.de/10.1594/PANGAEA.905295 and " + \
+        "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)
+    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()
+
+
+if (__name__ == "__main__"):
+    parser = argparse.ArgumentParser(
+        description=__doc__,
+        formatter_class=argparse.RawTextHelpFormatter)
+
+    parser.add_argument(
+        "--elev-path", dest="elev_path", type=str,
+        required=False,
+        default="", help="Path to raw DEM data-sets.")
+
+    parser.add_argument(
+        "--save-path", dest="save_path", type=str,
+        required=False,
+        default="", help="Path to store output data.")
+
+    parser.parse_args()
+    elev_path = parser.elev_path
+    save_path = parser.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)
+    """
+
+    gebco_60sec(elev_path, save_path)
+    gebco_30sec(elev_path, save_path)
+
+    rtopo_gebco_60sec(elev_path, save_path)
+    rtopo_gebco_30sec(elev_path, save_path)
diff --git a/compass/ocean/tests/tides/init/dem_remap.py b/compass/ocean/tests/tides/init/dem_remap.py
new file mode 100644
index 0000000000..7c67fe7065
--- /dev/null
+++ b/compass/ocean/tests/tides/init/dem_remap.py
@@ -0,0 +1,344 @@
+
+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
new file mode 100644
index 0000000000..1c41a2bacb
--- /dev/null
+++ b/compass/ocean/tests/tides/init/dem_trnsf.py
@@ -0,0 +1,114 @@
+
+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
new file mode 100644
index 0000000000..e88419dded
--- /dev/null
+++ b/compass/ocean/tests/tides/init/initial_state.py
@@ -0,0 +1,121 @@
+from compass.model import run_model
+from compass.step import Step
+
+from compass.ocean.vertical.grid_1d import generate_1d_grid, write_1d_grid
+from compass.ocean.plot import plot_vertical_grid
+import netCDF4
+import numpy as np
+
+
+class InitialState(Step):
+    """
+    A step for creating a mesh and initial condition for tidal
+    test cases
+
+    Attributes
+    ----------
+    mesh : compass.ocean.tests.tides.mesh.mesh.MeshStep
+        The step for creating the mesh
+
+    """
+    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 that creates the mesh used by this test case
+
+        """
+
+        super().__init__(test_case=test_case, name='initial_state')
+        self.mesh = mesh
+
+        package = 'compass.ocean.tests.tides.init'
+
+        # generate the namelist, replacing a few default options
+        self.add_namelist_file(package, 'namelist.init', mode='init')
+
+        # generate the streams file
+        self.add_streams_file(package, 'streams.init', mode='init')
+
+        mesh_path = mesh.steps['cull_mesh'].path
+
+        self.add_input_file(
+            filename='mesh.nc',
+            work_dir_target=f'{mesh_path}/culled_mesh.nc')
+
+        self.add_input_file(
+            filename='graph.info',
+            work_dir_target=f'{mesh_path}/culled_graph.info')
+
+        self.add_input_file(
+            target='PotentialTemperature.01.filled.60levels.PHC.151106.nc',
+            filename='temperature.nc',
+            database='initial_condition_database')
+
+        self.add_input_file(
+            target='Salinity.01.filled.60levels.PHC.151106.nc',
+            filename='salinity.nc',
+            database='initial_condition_database')
+
+        self.add_input_file(
+            target='windStress.ncep_1958-2000avg.interp3600x2431.151106.nc',
+            filename='wind_stress.nc',
+            database='initial_condition_database')
+
+        self.add_input_file(
+            target='chlorophyllA_monthly_averages_1deg.151201.nc',
+            filename='swData.nc',
+            database='initial_condition_database')
+
+        self.add_input_file(
+            target='BedMachineAntarctica_and_GEBCO_2019_0.05_degree.200128.nc',
+            filename='topography.nc',
+            database='bathymetry_database')
+
+        self.add_model_as_input()
+
+        for file in ['initial_state.nc', 'graph.info']:
+            self.add_output_file(filename=file)
+
+    def setup(self):
+        """
+        Set up the test case in the work directory, including downloading any
+        dependencies.
+        """
+        # get the these properties from the config options
+        config = self.config
+        self.ntasks = config.getint('tides', 'init_ntasks')
+        self.min_tasks = config.getint('tides', 'init_min_tasks')
+        self.threads = config.getint('tides', 'init_threads')
+
+    def run(self):
+        """
+        Run this step of the testcase
+        """
+
+        config = self.config
+        interfaces = generate_1d_grid(config=config)
+
+        write_1d_grid(interfaces=interfaces, out_filename='vertical_grid.nc')
+        plot_vertical_grid(grid_filename='vertical_grid.nc', config=config,
+                           out_filename='vertical_grid.png')
+
+        run_model(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["bottomDepth"][:] = \
+            np.minimum(max_depth,
+                       np.maximum(min_depth, init["bottomDepthObserved"][:]))
+
+        init["layerThickness"][0, :, 0] = init["bottomDepth"][:]
+        init["restingThickness"][:, 0] = init["bottomDepth"][:]
+        init.close()
diff --git a/compass/ocean/tests/tides/init/interpolate_wave_drag.py b/compass/ocean/tests/tides/init/interpolate_wave_drag.py
new file mode 100644
index 0000000000..86c7333ed8
--- /dev/null
+++ b/compass/ocean/tests/tides/init/interpolate_wave_drag.py
@@ -0,0 +1,185 @@
+from compass.step import Step
+from mpas_tools.mesh.interpolation import interp_bilin
+
+import netCDF4
+import matplotlib.pyplot as plt
+import numpy as np
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
+
+
+class InterpolateWaveDrag(Step):
+    """
+    A step for interpolating the topographic wave drag data
+    onto the MPAS-Ocean mesh
+
+    Attributes
+    ----------
+    plot : bool
+        Whether to produce plots of interpolated atmospheric data
+
+    rinv_file: str
+        Name of file for rinv data
+
+    wave_drag_file : str
+        Output file with interpolated rinv data
+
+    grid_file : str
+        Name of mesh file
+
+    """
+    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 that creates the mesh used by this test case
+        """
+        super().__init__(test_case=test_case, name='interpolate',
+                         ntasks=1, min_tasks=1, openmp_threads=1)
+
+        self.plot = True
+
+        self.rinv_file = 'jsl_lim24_inv_hrs.nc'
+        self.wave_drag_file = 'topographic_wave_drag.nc'
+        self.grid_file = 'mesh.nc'
+
+        self.add_input_file(
+            filename=self.rinv_file,
+            target=self.rinv_file,
+            database='initial_condition_database')
+
+        mesh_path = mesh.steps['cull_mesh'].path
+
+        self.add_input_file(
+            filename='mesh.nc',
+            work_dir_target=f'{mesh_path}/culled_mesh.nc')
+
+        self.add_output_file(filename=self.wave_drag_file)
+
+    def interpolate_data_to_grid(self, grid_file, data_file):
+        """
+        Interpolate time snaps of gridded data field to MPAS mesh
+        """
+
+        # 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'][:]
+        lon_data = np.append(lon_data, 180.0)
+        lon_data = np.insert(lon_data, 0, -180.0)
+        lat_data = data_nc.variables['Latitude'][:]
+
+        data = data_nc.variables['rinv'][:]
+        data = np.squeeze(data)
+        data = np.insert(data, 0, data[:, 0], axis=1)
+        data = np.append(data, data[:, -1][np.newaxis].T, axis=1)
+        data = 1.0/data[:, :]
+
+        # 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
+
+        # Interpolate
+        interp_data = interp_bilin(lon_data, lat_data,
+                                   data,
+                                   lon_grid, lat_grid)
+
+        return (lon_grid, lat_grid, interp_data), \
+            (lon_data, lat_data, data)
+
+    def plot_interp_data(self, orig_data, interp_data):
+        """
+        Plot original gridded data and interpolated fields
+        """
+
+        plt.switch_backend('agg')
+
+        if not self.plot:
+            return
+
+        lon_data = orig_data[0]
+        lat_data = orig_data[1]
+        lon_grid = interp_data[0]
+        lat_grid = interp_data[1]
+
+        data = orig_data[2][:, :]
+        interp = interp_data[2][:]
+
+        # Plot data
+        fig = plt.figure()
+        levels = np.linspace(np.amin(data), np.amax(data), 100)
+        ax0 = fig.add_subplot(2, 1, 1, projection=ccrs.PlateCarree())
+        cf = ax0.contourf(lon_data, lat_data, data, levels=levels,
+                          transform=ccrs.PlateCarree())
+        ax0.set_extent([0, 359.9, -90, 90], crs=ccrs.PlateCarree())
+        ax0.add_feature(cfeature.LAND, zorder=100)
+        ax0.add_feature(cfeature.LAKES, alpha=0.5, zorder=101)
+        ax0.add_feature(cfeature.COASTLINE, zorder=101)
+        ax0.set_title('data')
+        cbar = fig.colorbar(cf, ax=ax0)
+        cbar.set_label('rinv')
+
+        # Plot interpolated data
+        ax1 = fig.add_subplot(2, 1, 2, projection=ccrs.PlateCarree())
+        levels = np.linspace(np.amin(interp), np.amax(interp), 100)
+        cf = ax1.tricontourf(lon_grid, lat_grid, interp, levels=levels,
+                             transform=ccrs.PlateCarree())
+        ax1.set_extent([0, 359.9, -90, 90], crs=ccrs.PlateCarree())
+        ax1.add_feature(cfeature.LAND, zorder=100)
+        ax1.add_feature(cfeature.LAKES, alpha=0.5, zorder=101)
+        ax1.add_feature(cfeature.COASTLINE, zorder=101)
+        ax1.set_title('interpolated data')
+        cbar = fig.colorbar(cf, ax=ax1)
+        cbar.set_label('rinv')
+
+        # Save figure
+        fig.tight_layout()
+        fig.savefig('rinv.png',
+                    bbox_inches='tight')
+        plt.close()
+
+    def write_to_file(self, data):
+        """
+        Write data to netCDF file
+        """
+
+        data_nc = netCDF4.Dataset(self.wave_drag_file, 'w',
+                                  format='NETCDF3_64BIT_OFFSET')
+
+        # Find dimesions
+        nEdges = data.shape[0]
+
+        # Declare dimensions
+        data_nc.createDimension('nEdges', nEdges)
+
+        # Set variables
+        data_var = data_nc.createVariable('topographic_wave_drag',
+                                          np.float64,
+                                          ('nEdges'))
+        data_var[:] = data[:]
+        data_nc.close()
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+
+        # Interpolation of rinv
+        interp, data = self.interpolate_data_to_grid(
+            self.grid_file,
+            self.rinv_file)
+
+        self.write_to_file(interp[2])
+
+        # Plot rinv field
+        self.plot_interp_data(data, interp)
diff --git a/compass/ocean/tests/tides/init/namelist.init b/compass/ocean/tests/tides/init/namelist.init
new file mode 100644
index 0000000000..853428f1c0
--- /dev/null
+++ b/compass/ocean/tests/tides/init/namelist.init
@@ -0,0 +1,51 @@
+config_ocean_run_mode = 'init'
+config_init_configuration = 'global_ocean'
+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_topography_source = 'mpas_variable'
+config_global_ocean_depth_file = 'vertical_grid.nc'
+config_global_ocean_temperature_file = 'temperature.nc'
+config_global_ocean_salinity_file = 'salinity.nc'
+config_global_ocean_swData_file = 'swData.nc'
+config_global_ocean_topography_file = 'topography.nc'
+config_global_ocean_windstress_file = 'wind_stress.nc'
+config_global_ocean_tracer_nlat_dimname = 't_lat'
+config_global_ocean_tracer_nlon_dimname = 't_lon'
+config_global_ocean_tracer_ndepth_dimname = 'depth_t'
+config_global_ocean_temperature_varname = 'TEMP'
+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_swData_nlat_dimname = 't_lat'
+config_global_ocean_swData_nlon_dimname = 't_lon'
+config_global_ocean_swData_lat_varname = 't_lat'
+config_global_ocean_swData_lon_varname = 't_lon'
+config_global_ocean_chlorophyll_varname = 'Chlorophyll'
+config_global_ocean_zenithAngle_varname = 'zenithAngle'
+config_global_ocean_clearSky_varname = 'clearSky'
+config_global_ocean_topography_nlat_dimname = 'y'
+config_global_ocean_topography_nlon_dimname = 'x'
+config_global_ocean_topography_lat_varname = 'y'
+config_global_ocean_topography_lon_varname = 'x'
+config_global_ocean_topography_varname = 'z'
+config_global_ocean_deepen_critical_passages = .false.
+config_global_ocean_cull_inland_seas = .false.
+config_global_ocean_windstress_nlat_dimname = 'u_lat'
+config_global_ocean_windstress_nlon_dimname = 'u_lon'
+config_global_ocean_windstress_lat_varname = 'u_lat'
+config_global_ocean_windstress_lon_varname = 'u_lon'
+config_global_ocean_windstress_zonal_varname = 'TAUX'
+config_global_ocean_windstress_meridional_varname = 'TAUY'
+config_global_ocean_ecosys_file = 'ecosys.nc'
+config_global_ocean_ecosys_forcing_file = 'ecosys_forcing.nc'
+config_global_ocean_ecosys_nlat_dimname = 't_lat'
+config_global_ocean_ecosys_nlon_dimname = 't_lon'
+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'
diff --git a/compass/ocean/tests/tides/init/remap_bathymetry.py b/compass/ocean/tests/tides/init/remap_bathymetry.py
new file mode 100644
index 0000000000..6f3a077507
--- /dev/null
+++ b/compass/ocean/tests/tides/init/remap_bathymetry.py
@@ -0,0 +1,82 @@
+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):
+    """
+    A step for remapping bathymetric data onto the MPAS-Ocean mesh
+
+    Attributes
+    ----------
+    plot : bool
+        Whether to produce plots of remapped data
+    """
+    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 that creates the mesh used by this test case
+        """
+        super().__init__(test_case=test_case, name='remap',
+                         ntasks=1, min_tasks=1, openmp_threads=1)
+
+        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
+
+        self.add_input_file(
+            filename='mesh.nc',
+            work_dir_target=f'{cull_mesh_path}/culled_mesh.nc')
+        self.add_input_file(
+            filename='base_mesh.nc',
+            work_dir_target=f'{base_mesh_path}/base_mesh.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',
+                            'base_mesh.nc')
+        dem_trnsf.dem_trnsf('base_mesh.nc', 'mesh.nc')
diff --git a/compass/ocean/tests/tides/init/streams.init b/compass/ocean/tests/tides/init/streams.init
new file mode 100644
index 0000000000..a8e2d9e6a0
--- /dev/null
+++ b/compass/ocean/tests/tides/init/streams.init
@@ -0,0 +1,105 @@
+<streams>
+
+<immutable_stream name="input_init"
+                  type="input"
+                  filename_template="mesh.nc"
+                  input_interval="initial_only"/>
+
+<immutable_stream name="atmospheric_forcing"
+                  type="input"
+                  filename_template="atmospheric_forcing.nc"
+                  filename_interval="none"
+                  packages="timeVaryingAtmosphericForcingPKG"
+                  input_interval="none"/>
+
+<immutable_stream name="land_ice_forcing"
+                  type="input"
+                  filename_template="land_ice_forcing.nc"
+                  filename_interval="none"
+                  packages="timeVaryingLandIceForcingPKG"
+                  input_interval="none"/>
+
+<stream name="output_init"
+        type="output"
+        output_interval="0000_00:00:01"
+        clobber_mode="truncate"
+        filename_template="initial_state.nc">
+
+    <stream name="input_init"/>
+    <var_struct name="tracers"/>
+    <var name="refZMid"/>
+    <var name="normalVelocity"/>
+    <var name="layerThickness"/>
+    <var name="restingThickness"/>
+    <var name="seaIcePressure"/>
+    <var name="atmosphericPressure"/>
+    <var name="boundaryLayerDepth"/>
+    <var name="refBottomDepth"/>
+    <var name="bottomDepth"/>
+    <var name="bottomDepthObserved"/>
+    <var name="oceanFracObserved"/>
+    <var name="maxLevelCell"/>
+    <var name="vertCoordMovementWeights"/>
+    <var name="edgeMask"/>
+    <var name="cullCell"/>
+    <var name="landIceDraftObserved"/>
+    <var name="landIceThkObserved"/>
+    <var name="landIceFracObserved"/>
+    <var name="landIceGroundedFracObserved"/>
+    <var name="landIceFraction"/>
+    <var name="landIceMask"/>
+    <var name="landIcePressure"/>
+    <var name="landIceDraft"/>
+    <var name="ssh"/>
+    <var name="modifyLandIcePressureMask"/>
+    <var name="rx1Cell"/>
+    <var name="rx1Edge"/>
+    <var name="rx1MaxCell"/>
+    <var name="rx1MaxEdge"/>
+    <var name="globalRx1Max"/>
+    <var name="density"/>
+    <var name="rx1InitSmoothingMask"/>
+    <var name="verticalStretch"/>
+</stream>
+
+<stream name="forcing_data_init"
+        type="output"
+        output_interval="0000_00:00:01"
+        clobber_mode="truncate"
+        filename_template="init_mode_forcing_data.nc">
+
+    <var_struct name="tracersSurfaceRestoringFields"/>
+    <var_struct name="tracersInteriorRestoringFields"/>
+    <var_struct name="tracersExponentialDecayFields"/>
+    <var_struct name="tracersIdealAgeFields"/>
+    <var_struct name="tracersTTDFields"/>
+    <var name="windStressZonal"/>
+    <var name="windStressMeridional"/>
+    <var name="landIceSurfaceTemperature"/>
+    <var name="seaIcePressure"/>
+    <var name="atmosphericPressure"/>
+</stream>
+
+<stream name="shortwave_forcing_data_init"
+        type="output"
+        output_interval="0000-00-00_00:00:01"
+        clobber_mode="truncate"
+        filename_template="init_mode_shortwaveData.nc">
+
+    <var name="chlorophyllData"/>
+    <var name="zenithAngle"/>
+    <var name="clearSkyRadiation"/>
+    <var name="xtime"/>
+</stream>
+
+<stream name="critical_passages"
+        filename_template="critical_passages.nc"
+        input_interval="initial_only"
+        type="input">
+
+    <var name="transectCellMasks"/>
+    <var name="depthTransects"/>
+</stream>
+
+</streams>
+
diff --git a/compass/ocean/tests/tides/mesh/__init__.py b/compass/ocean/tests/tides/mesh/__init__.py
new file mode 100644
index 0000000000..96d0663877
--- /dev/null
+++ b/compass/ocean/tests/tides/mesh/__init__.py
@@ -0,0 +1,49 @@
+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
+
+
+class Mesh(TestCase):
+    """
+    A test case for creating a global MPAS-Ocean mesh
+    """
+    def __init__(self, test_group, mesh_name):
+        """
+        Create test case for creating a global MPAS-Ocean mesh
+
+        Parameters
+        ----------
+        test_group : compass.ocean.tests.tides.Tides
+            The test group that this test case belongs to
+
+        mesh_name : str
+            The name of the mesh
+        """
+        self.mesh_name = mesh_name
+        name = 'mesh'
+        subdir = '{}/{}'.format(mesh_name, name)
+        super().__init__(test_group=test_group, name=name, subdir=subdir)
+
+        name = 'base_mesh'
+        if mesh_name == 'Icos7':
+            base_mesh_step = IcosahedralMeshStep(
+                self, name=name, subdivisions=7)
+            mesh_lower = 'icos7'
+        else:
+            raise ValueError(f'Unexpected mesh name {mesh_name}')
+
+        self.package = f'compass.ocean.tests.tides.mesh.{mesh_lower}'
+        self.mesh_config_filename = f'{mesh_lower}.cfg'
+
+        self.add_step(base_mesh_step)
+
+        self.add_step(CullMeshStep(
+            test_case=self, base_mesh_step=base_mesh_step,
+            with_ice_shelf_cavities=True))
+
+    def configure(self):
+        """
+        Modify the configuration options for this test case
+        """
+        configure_tides(test_case=self, mesh=self)
diff --git a/compass/ocean/tests/tides/mesh/icos7/__init__.py b/compass/ocean/tests/tides/mesh/icos7/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/compass/ocean/tests/tides/mesh/icos7/icos7.cfg b/compass/ocean/tests/tides/mesh/icos7/icos7.cfg
new file mode 100644
index 0000000000..5eda375b38
--- /dev/null
+++ b/compass/ocean/tests/tides/mesh/icos7/icos7.cfg
@@ -0,0 +1,56 @@
+# 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
+
+
+# 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 = 6000
+min_depth = 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 = 180
+# minimum of cores, below which the step fails
+forward_min_tasks = 160
+# 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/compass/ocean/tests/tides/mesh/icos7/namelist.ocean b/compass/ocean/tests/tides/mesh/icos7/namelist.ocean
new file mode 100644
index 0000000000..429e61daf0
--- /dev/null
+++ b/compass/ocean/tests/tides/mesh/icos7/namelist.ocean
@@ -0,0 +1,2 @@
+config_dt = '00:01:00'
+config_topographic_wave_drag_coeff = 5e-4
diff --git a/conda/compass_env/spec-file.template b/conda/compass_env/spec-file.template
index 410920c796..f26d922c1e 100644
--- a/conda/compass_env/spec-file.template
+++ b/conda/compass_env/spec-file.template
@@ -28,7 +28,7 @@ pyamg >=4.2.2
 pyproj
 pyremap>=0.0.13,<0.1.0
 requests
-scipy
+scipy>=1.8.0
 shapely>=1.8.0,<2.0.0
 xarray
 
diff --git a/conda/recipe/meta.yaml b/conda/recipe/meta.yaml
index fe08d871f6..1715c4806c 100644
--- a/conda/recipe/meta.yaml
+++ b/conda/recipe/meta.yaml
@@ -65,7 +65,7 @@ requirements:
     - pyproj
     - pyremap >=0.0.13,<0.1.0
     - requests
-    - scipy
+    - scipy >=1.8.0
     - shapely >=1.8.0,<2.0.0
     - xarray
 
diff --git a/docs/developers_guide/ocean/api.rst b/docs/developers_guide/ocean/api.rst
index 6e53b49e59..368d156698 100644
--- a/docs/developers_guide/ocean/api.rst
+++ b/docs/developers_guide/ocean/api.rst
@@ -677,6 +677,57 @@ sphere_transport framework
    process_output.plot_filament
    process_output.plot_over_and_undershoot_errors
 
+tides
+~~~~~
+
+test cases and steps
+''''''''''''''''''''
+
+.. currentmodule:: compass.ocean.tests.tides
+
+.. autosummary::
+   :toctree: generated/
+
+   Tides
+   configure
+
+   mesh.Mesh
+   mesh.Mesh.configure
+   mesh.Mesh.run
+
+   init.Init
+   init.Init.configure
+   init.Init.run
+   init.remap_bathymetry.RemapBathymetry
+   init.remap_bathymetry.RemapBathymetry.run
+   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
+
+
+   forward.Forward
+   forward.Forward.configure
+   forward.Forward.run
+   forward.forward.ForwardStep
+   forward.forward.ForwardStep.setup
+   forward.forward.ForwardStep.run
+
+   analysis.Analysis
+   analysis.Analysis.setup
+   analysis.Analysis.write_coordinate_file
+   analysis.Analysis.setup_otps2
+   analysis.Analysis.run_otps2
+   analysis.Analysis.read_otps2_output
+   analysis.Analysis.append_tpxo_data
+   analysis.Analysis.check_tpxo_data
+   analysis.Analysis.plot
+   analysis.Analysis.run
+
 ziso
 ~~~~
 
diff --git a/docs/developers_guide/ocean/test_groups/index.rst b/docs/developers_guide/ocean/test_groups/index.rst
index b174ba71ee..e0c1c5c5b8 100644
--- a/docs/developers_guide/ocean/test_groups/index.rst
+++ b/docs/developers_guide/ocean/test_groups/index.rst
@@ -21,5 +21,6 @@ Test groups
    planar_convergence
    soma
    sphere_transport
-   ziso
    spherical_harmonic_transform
+   tides
+   ziso
diff --git a/docs/developers_guide/ocean/test_groups/tides.rst b/docs/developers_guide/ocean/test_groups/tides.rst
new file mode 100644
index 0000000000..f5418f5384
--- /dev/null
+++ b/docs/developers_guide/ocean/test_groups/tides.rst
@@ -0,0 +1,38 @@
+.. _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.
+
+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.
+
+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.
+
+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.png
new file mode 100644
index 0000000000..b457b1c0e2
Binary files /dev/null and b/docs/users_guide/ocean/test_groups/images/M2_plot.png differ
diff --git a/docs/users_guide/ocean/test_groups/index.rst b/docs/users_guide/ocean/test_groups/index.rst
index e434bb4f31..4787512eae 100644
--- a/docs/users_guide/ocean/test_groups/index.rst
+++ b/docs/users_guide/ocean/test_groups/index.rst
@@ -25,5 +25,6 @@ coming months.
    planar_convergence
    soma
    sphere_transport
-   ziso
    spherical_harmonic_transform
+   tides
+   ziso
diff --git a/docs/users_guide/ocean/test_groups/tides.rst b/docs/users_guide/ocean/test_groups/tides.rst
new file mode 100644
index 0000000000..c1389b293f
--- /dev/null
+++ b/docs/users_guide/ocean/test_groups/tides.rst
@@ -0,0 +1,137 @@
+.. _ocean_tides:
+
+tides
+=====
+
+The ``ocean/tides`` test group defines meshes,
+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.
+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
+more mesh resolutions to be added in the future.
+
+Shared config options
+---------------------
+
+All ``tides`` test cases start the following shared config options.
+Note that meshes and test cases may modify these options, as noted below.
+
+.. code-block:: cfg
+
+    # 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
+    
+    
+    # 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 = 6000
+    
+    
+    # options for tides testcases
+    [tides]
+    
+    ## config options related to the initial_state step
+    # number of cores to use
+    init_cores = 36
+    # minimum of cores, below which the step fails
+    init_min_cores = 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 = 180
+    # minimum of cores, below which the step fails
+    forward_min_tasks = 160
+    # maximum memory usage allowed (in MB)
+    forward_max_memory = 1000
+    # number of threads
+    forward_threads = 1 
+    
+    # TPXO version for validation 
+    tpxo_version = TPXO9
+
+
+.. _tides_mesh:
+
+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.
+
+.. _tides_init:
+
+init test case
+--------------
+The init test performs steps to set up the vertical mesh, initial conditions,
+atmospheric forcing, and prepares the station locations for timeseries output.
+
+remap bathymetry step
+^^^^^^^^^^^^^^^^^^^^^
+This step performs an integral remap of bathymetric data onto the MPAS-O mesh
+
+interpolate wave drag step
+^^^^^^^^^^^^^^^^^^^^^^^^^^
+In this step, a timescale factor from HYCOM data is interpolated onto
+the MPAS-O mesh to be used in the topographic wave drag parameterization.
+A plot is also produced to verify the interpolation.
+
+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.
+
+.. _tides_forward:
+
+forward test case
+-----------------
+The forward test case is responsible for the forward model simulation and
+analysis.
+
+forward step
+^^^^^^^^^^^^
+The forward step runs a 125 day model simulation. The simulation
+begins with a spinup period, where the tides 
+are ramped to their full value to avoid shocking the system.
+The harmonic analysis period begins following this spinup period.
+
+analysis step
+^^^^^^^^^^^^^
+The analysis step extracts the harmonic constituent amplitude and phase values
+at each cell center from the TPXO database. The TPXO values are used to compute
+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
+   :width: 500px
+   :align: center
diff --git a/setup.py b/setup.py
index ad0f22f66e..1b7cda35eb 100644
--- a/setup.py
+++ b/setup.py
@@ -32,7 +32,7 @@ def package_files(directory, prefixes, extensions):
      'pyamg',
      'pyproj',
      'requests',
-     'scipy',
+     'scipy>=1.8.0',
      'shapely',
      'xarray']