diff --git a/compass/ocean/__init__.py b/compass/ocean/__init__.py
index 3c50d5125e..5d86c11962 100644
--- a/compass/ocean/__init__.py
+++ b/compass/ocean/__init__.py
@@ -1,4 +1,5 @@
 from compass.mpas_core import MpasCore
+from compass.ocean.tests.turbulence_closure import TurbulenceClosure
 from compass.ocean.tests.baroclinic_channel import BaroclinicChannel
 from compass.ocean.tests.dam_break import DamBreak
 from compass.ocean.tests.drying_slope import DryingSlope
@@ -32,6 +33,7 @@ def __init__(self):
         """
         super().__init__(name='ocean')
 
+        self.add_test_group(TurbulenceClosure(mpas_core=self))
         self.add_test_group(BaroclinicChannel(mpas_core=self))
         self.add_test_group(DamBreak(mpas_core=self))
         self.add_test_group(DryingSlope(mpas_core=self))
diff --git a/compass/ocean/tests/turbulence_closure/__init__.py b/compass/ocean/tests/turbulence_closure/__init__.py
new file mode 100644
index 0000000000..21cd27881d
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/__init__.py
@@ -0,0 +1,85 @@
+from compass.ocean.tests.turbulence_closure.decomp_test import DecompTest
+from compass.ocean.tests.turbulence_closure.default import Default
+from compass.ocean.tests.turbulence_closure.restart_test import RestartTest
+from compass.testgroup import TestGroup
+
+
+class TurbulenceClosure(TestGroup):
+    """
+    A test group for turbulence closure test cases
+    """
+    def __init__(self, mpas_core):
+        """
+        mpas_core : compass.MpasCore
+            the MPAS core that this test group belongs to
+        """
+        super().__init__(mpas_core=mpas_core, name='turbulence_closure')
+
+        for resolution in [1e4]:
+            self.add_test_case(
+                DecompTest(test_group=self, resolution=resolution))
+            self.add_test_case(
+                RestartTest(test_group=self, resolution=resolution))
+        for resolution in [1, 2, 1e4]:
+            for forcing in ['cooling', 'evaporation']:
+                self.add_test_case(
+                    Default(test_group=self, resolution=resolution,
+                            forcing=forcing))
+
+
+def configure(resolution, forcing, config):
+    """
+    Modify the configuration options for one of the turbulence closure test
+    cases
+
+    Parameters
+    ----------
+    resolution : float
+        The resolution of the test case in meters
+
+    config : configparser.ConfigParser
+        Configuration options for this test case
+    """
+    # The resolution parameters are different for different resolutions
+    # to match existing simulations
+    if resolution > 1e3:
+        nx = 16
+        ny = 50
+        vert_levels = 20
+        bottom_depth = 1e3
+    elif resolution <= 1e3 and resolution > 5:
+        nx = 50
+        ny = 50
+        vert_levels = 50
+        bottom_depth = 100.0
+    elif resolution <= 5 and resolution > 1:
+        nx = 150
+        ny = 150
+        vert_levels = 50
+        bottom_depth = 100.0
+    elif resolution <= 1:
+        nx = 128
+        ny = 128
+        vert_levels = 128
+        bottom_depth = 128.0
+
+    config.set('turbulence_closure', 'nx', f'{nx}')
+    config.set('turbulence_closure', 'ny', f'{ny}')
+    config.set('turbulence_closure', 'dc', f'{resolution}')
+    config.set('vertical_grid', 'vert_levels', f'{vert_levels}')
+    config.set('vertical_grid', 'bottom_depth', f'{bottom_depth}')
+
+    if forcing == 'cooling':
+        config.set('turbulence_closure', 'surface_heat_flux', '-100')
+        config.set('turbulence_closure', 'surface_freshwater_flux', '0')
+        config.set('turbulence_closure', 'interior_temperature_gradient',
+                   '0.1')
+        config.set('turbulence_closure', 'interior_salinity_gradient', '0')
+        config.set('turbulence_closure', 'wind_stress_zonal', '0')
+    if forcing == 'evaporation':
+        config.set('turbulence_closure', 'surface_heat_flux', '0')
+        config.set('turbulence_closure', 'surface_freshwater_flux', '0.429')
+        config.set('turbulence_closure', 'interior_temperature_gradient', '0')
+        config.set('turbulence_closure', 'interior_salinity_gradient',
+                   '-0.025')
+        config.set('turbulence_closure', 'wind_stress_zonal', '0')
diff --git a/compass/ocean/tests/turbulence_closure/boundary_values.py b/compass/ocean/tests/turbulence_closure/boundary_values.py
new file mode 100644
index 0000000000..c44fb3c0f5
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/boundary_values.py
@@ -0,0 +1,210 @@
+import numpy
+
+
+def add_boundary_arrays(mesh, x_nonperiodic, y_nonperiodic):  # noqa: C901
+    """
+    Computes boundary arrays necessary for LES test cases
+    Values are appended to a MPAS mesh dataset
+
+    if the domain is doubly periodic, returns array of all zeros
+    as arrays are meaningless in that case.
+
+    Parameters
+    ----------
+    mesh : xarray dataset containing an MPAS mesh
+
+    x_nonperiodic : logical if x direction is non periodic
+
+    y_nonperiodic : logical if y direction is non periodic
+    """
+
+    nCells = mesh.dims['nCells']
+    neonc = mesh.nEdgesOnCell.values
+    eonc = mesh.edgesOnCell.values - 1
+    conc = mesh.cellsOnCell.values - 1
+    nVertLevels = mesh.dims['nVertLevels']
+    mlc = mesh.maxLevelCell.values - 1
+    cone = mesh.cellsOnEdge.values - 1
+
+    if (not x_nonperiodic) and (not y_nonperiodic):
+        mesh['distanceToBoundary'] = (
+            ['nCells', 'nVertLevels'],
+            numpy.zeros((nCells, nVertLevels)))
+        mesh['boundaryNormalAngle'] = (
+            ['nCells', 'nVertLevels'],
+            numpy.zeros((nCells, nVertLevels)))
+        mesh['boundaryCellMask'] = (
+            ['nCells', 'nVertLevels'],
+            numpy.zeros((nCells, nVertLevels)).astype(int))
+        return mesh
+
+    # surface first and save until later
+    indVals = []
+    angle = mesh.angleEdge.values
+    xc = mesh.xEdge.values
+    yc = mesh.yEdge.values
+    xCell = mesh.xCell.values
+    yCell = mesh.yCell.values
+    aCell = numpy.zeros(nCells)
+    bCell = numpy.zeros(nCells)
+    xEdge = []
+    yEdge = []
+    for i in range(nCells):
+        if conc[i, :].min() < 0:
+            aEdge = 0
+            count = 0
+            for j in range(neonc[i]):
+                if conc[i, j] < 0:
+                    indVals.append(i)
+                    xEdge.append(xc[eonc[i, j]])
+                    yEdge.append(yc[eonc[i, j]])
+                    # compute dx, dy
+                    # check angle edge
+                    # add pi if needed
+                    # check if over 2pi correct if needed
+                    dx = xc[eonc[i, j]] - xCell[i]
+                    dy = yc[eonc[i, j]] - yCell[i]
+                    if dx <= 0 and dy <= 0:  # first quadrant
+                        if angle[eonc[i, j]] >= numpy.pi / 2.:
+                            aEdge += angle[eonc[i, j]] - numpy.pi
+                            count += 1
+                        else:
+                            aEdge += angle[eonc[i, j]]
+                            count += 1
+
+                    # fourth quadrant, but reference to 0 not 2pi
+                    elif dx <= 0 and dy >= 0:
+                        if angle[eonc[i, j]] < 3.0 * numpy.pi / 2.0:
+                            aEdge += angle[eonc[i, j]] + numpy.pi
+                            count += 1
+                        else:
+                            aEdge += angle[eonc[i, j]]
+                            count += 1
+
+                    elif dx >= 0 and dy >= 0:  # third quadrant
+                        if angle[eonc[i, j]] < numpy.pi / 2.0:
+                            # not in correct place
+                            aEdge += angle[eonc[i, j]] + numpy.pi
+                            count += 1
+                        else:
+                            aEdge += angle[eonc[i, j]]
+                            count += 1
+
+                    else:  # quadrant 2
+                        if angle[eonc[i, j]] > 3.0 * numpy.pi / 2.0:
+                            # wrong place
+                            aEdge += angle[eonc[i, j]] - numpy.pi
+                            count += 1
+                        else:
+                            aEdge += angle[eonc[i, j]]
+                            count += 1
+            if count > 0:
+                aCellTemp = aEdge / count
+                if aCellTemp > numpy.pi:
+                    aCellTemp = 2.0 * numpy.pi - aCellTemp
+                    aCell[i] = aCellTemp
+                    bCell[i] = 1
+
+    # with angle and index arrays, now can do distance
+    dist = numpy.zeros((nCells, nVertLevels))
+    angleCells = numpy.zeros((nCells, nVertLevels))
+    boundaryCells = numpy.zeros((nCells, nVertLevels))
+    for i in range(nCells):
+        d = numpy.sqrt((xCell[i] - xEdge)**2 + (yCell[i] - yEdge)**2)
+        ind = d.argmin()
+        angleCells[i, 0] = aCell[indVals[ind]]
+        boundaryCells[i, 0] = bCell[i]
+        dist[i, 0] = d.min()
+
+    # now to the rest
+    for k in range(1, nVertLevels):
+        inds = numpy.where(k == mlc)[0]
+        edgeList = []
+        cellList = []
+        if len(inds) > 0:
+            # First step is forming list of edges bordering maxLC
+            for i in range(len(inds)):
+                for j in range(neonc[inds[i]]):
+                    if conc[inds[i], j] not in inds:
+                        edgeList.append(eonc[inds[i], j])
+            for i in range(len(edgeList)):
+                c1 = cone[edgeList[i], 0]
+                c2 = cone[edgeList[i], 1]
+                if c1 in inds:
+                    if c2 not in cellList:
+                        cellList.append(c2)
+                if c2 in inds:
+                    if c1 not in cellList:
+                        cellList.append(c1)
+
+            for i in range(len(cellList)):
+                aEdge = 0
+                count = 0
+                for j in range(neonc[cellList[i]]):
+                    if eonc[cellList[i], j] in edgeList:
+                        indVals.append(cellList[i])
+                        xEdge.append(xc[eonc[cellList[i], j]])
+                        yEdge.append(yc[eonc[cellList[i], j]])
+                        # compute dx, dy
+                        # check angle edge
+                        # add pi if needed
+                        # check if over 2pi correct if needed
+                        dx = xc[eonc[cellList[i], j]] - xCell[cellList[i]]
+                        dy = yc[eonc[cellList[i], j]] - yCell[cellList[i]]
+                        if dx <= 0 and dy <= 0:  # first quadrant
+                            if angle[eonc[cellList[i], j]] >= numpy.pi / 2.:
+                                aEdge += angle[eonc[cellList[i], j]] - numpy.pi
+                                count += 1
+                            else:
+                                aEdge += angle[eonc[cellList[i], j]]
+                                count += 1
+                        # fourth quadrant, but reference to 0 not 2pi
+                        elif dx <= 0 and dy >= 0:
+                            if {(angle[eonc[cellList[i], j]] <
+                                 3.0 * numpy.pi / 2.0)}:
+                                aEdge += angle[eonc[cellList[i], j]] + numpy.pi
+                                count += 1
+                            else:
+                                aEdge += angle[eonc[cellList[i], j]]
+                                count += 1
+                        elif dx >= 0 and dy >= 0:  # third quadrant
+                            # not in correct place
+                            if angle[eonc[cellList[i], j]] < numpy.pi / 2.0:
+                                aEdge += angle[eonc[cellList[i], j]] + numpy.pi
+                                count += 1
+                            else:
+                                aEdge += angle[eonc[cellList[i], j]]
+                                count += 1
+                        else:  # quadrant 2
+                            # wrong place
+                            if {(angle[eonc[cellList[i], j]] >
+                                 3.0 * numpy.pi / 2.0)}:
+                                aEdge += angle[eonc[cellList[i], j]] - numpy.pi
+                                count += 1
+                            else:
+                                aEdge += angle[eonc[cellList[i], j]]
+                                count += 1
+                if count > 0:
+                    aCellTemp = aEdge / count
+                    if aCellTemp > numpy.pi:
+                        aCellTemp = 2.0 * numpy.pi - aCellTemp
+                        aCell[cellList[i]] = aCellTemp
+                        bCell[cellList[i]] = 1
+
+                for ii in range(nCells):
+                    d = numpy.sqrt((xCell[ii] - xEdge)**2 +
+                                   (yCell[ii] - yEdge)**2)
+                    ind = d.argmin()
+                    angleCells[ii, k] = aCell[indVals[ind]]
+                    boundaryCells[ii, k] = bCell[ii]
+                    dist[ii, k] = d.min()
+        else:
+            dist[:, k] = dist[:, 0]
+            angleCells[:, k] = angleCells[:, 0]
+            boundaryCells[:, k] = boundaryCells[:, 0]
+    mesh['distanceToBoundary'] = (['nCells', 'nVertLevels'], dist)
+    mesh['boundaryNormalAngle'] = (['nCells', 'nVertLevels'], angleCells)
+    mesh['boundaryCellMask'] = (['nCells', 'nVertLevels'],
+                                boundaryCells.astype(int))
+
+    return mesh
diff --git a/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py b/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py
new file mode 100644
index 0000000000..e0f7ae809b
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py
@@ -0,0 +1,67 @@
+from compass.testcase import TestCase
+from compass.ocean.tests.turbulence_closure.initial_state import InitialState
+from compass.ocean.tests.turbulence_closure.forward import Forward
+from compass.ocean.tests import turbulence_closure
+from compass.validate import compare_variables
+
+
+class DecompTest(TestCase):
+    """
+    A decomposition test case for the turbulence closure test group, which
+    makes sure the model produces identical results on 1 and 4 cores.
+
+    Attributes
+    ----------
+    resolution : float
+        The resolution of the test case in meters
+    """
+
+    def __init__(self, test_group, resolution, forcing='cooling'):
+        """
+        Create the test case
+
+        Parameters
+        ----------
+        test_group : compass.ocean.tests.turbulence_closure.TurbulenceClosure
+            The test group that this test case belongs to
+
+        resolution : float
+            The resolution of the test case in meters
+        """
+        name = 'decomp_test'
+        self.resolution = resolution
+        if resolution >= 1e3:
+            res_name = f'{int(resolution/1e3)}km'
+        else:
+            res_name = f'{int(resolution)}m'
+        subdir = f'{res_name}/{forcing}/{name}'
+        super().__init__(test_group=test_group, name=name,
+                         subdir=subdir)
+
+        self.add_step(
+            InitialState(test_case=self, resolution=resolution))
+
+        for procs in [4, 8]:
+            name = '{}proc'.format(procs)
+            self.add_step(
+                Forward(test_case=self, name=name, subdir=name, ntasks=procs,
+                        openmp_threads=1, resolution=resolution))
+
+    def configure(self):
+        """
+        Modify the configuration options for this test case.
+        """
+        turbulence_closure.configure(self.resolution, self.config)
+
+    # no run() method is needed
+
+    def validate(self):
+        """
+        Test cases can override this method to perform validation of variables
+        and timers
+        """
+        variables = ['temperature', 'salinity', 'layerThickness',
+                     'normalVelocity']
+        compare_variables(test_case=self, variables=variables,
+                          filename1='4proc/output.nc',
+                          filename2='8proc/output.nc')
diff --git a/compass/ocean/tests/turbulence_closure/default/__init__.py b/compass/ocean/tests/turbulence_closure/default/__init__.py
new file mode 100644
index 0000000000..70483f0418
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/default/__init__.py
@@ -0,0 +1,68 @@
+from compass.ocean.tests import turbulence_closure
+from compass.ocean.tests.turbulence_closure.forward import Forward
+from compass.ocean.tests.turbulence_closure.initial_state import InitialState
+from compass.ocean.tests.turbulence_closure.viz import Viz
+from compass.testcase import TestCase
+
+
+class Default(TestCase):
+    """
+    The default test case for the turbulence closure test group simply creates
+    the mesh and initial condition, then performs a short forward run on 4
+    cores.
+
+    Attributes
+    ----------
+    resolution : float
+        The resolution of the test case in meters
+    """
+
+    def __init__(self, test_group, resolution, forcing='cooling'):
+        """
+        Create the test case
+
+        Parameters
+        ----------
+        test_group : compass.ocean.tests.turbulence_closure.TurbulenceClosure
+            The test group that this test case belongs to
+
+        resolution : float
+            The resolution of the test case in meters
+
+        forcing: str
+            The forcing applied to the test case
+        """
+        name = 'default'
+        self.resolution = resolution
+        self.forcing = forcing
+        if resolution >= 1e3:
+            res_name = f'{int(resolution/1e3)}km'
+        else:
+            res_name = f'{int(resolution)}m'
+        subdir = f'{res_name}/{forcing}/{name}'
+        super().__init__(test_group=test_group, name=name,
+                         subdir=subdir)
+
+        if resolution <= 5:
+            ntasks = 128
+            plot_comparison = True
+        else:
+            ntasks = 4
+            plot_comparison = False
+
+        self.add_step(
+            InitialState(test_case=self, resolution=resolution))
+        self.add_step(
+            Forward(test_case=self, ntasks=ntasks, openmp_threads=1,
+                    resolution=resolution))
+        self.add_step(Viz(test_case=self, resolution=resolution,
+                          forcing=forcing, do_comparison=plot_comparison))
+
+    def configure(self):
+        """
+        Modify the configuration options for this test case.
+        """
+        turbulence_closure.configure(self.resolution, self.forcing,
+                                     self.config)
+
+    # no run() is needed because we're doing the default: running all steps
diff --git a/compass/ocean/tests/turbulence_closure/forward.py b/compass/ocean/tests/turbulence_closure/forward.py
new file mode 100644
index 0000000000..aac73f1af6
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/forward.py
@@ -0,0 +1,97 @@
+from compass.model import run_model
+from compass.step import Step
+
+
+class Forward(Step):
+    """
+    A step for performing forward MPAS-Ocean runs as part of turbulence
+    closure test cases.
+
+    Attributes
+    ----------
+    resolution : str
+        The resolution of the test case
+    """
+    def __init__(self, test_case, resolution, name='forward', subdir=None,
+                 ntasks=1, min_tasks=None, openmp_threads=1, nu=None):
+        """
+        Create a new test case
+
+        Parameters
+        ----------
+        test_case : compass.TestCase
+            The test case this step belongs to
+
+        resolution : str
+            The resolution of the test case
+
+        name : str
+            the name of the test case
+
+        subdir : str, optional
+            the subdirectory for the step.  The default is ``name``
+
+        ntasks: int, optional
+            the number of tasks the step would ideally use.  If fewer tasks
+            are available on the system, the step will run on all available
+            cores as long as this is not below ``min_tasks``
+
+        min_tasks : int, optional
+            the number of tasks the step requires.  If the system has fewer
+            than this number of tasks, the step will fail
+
+        openmp_threads : int, optional
+            the number of OpenMP threads the step will use
+        """
+        self.resolution = resolution
+        if min_tasks is None:
+            min_tasks = ntasks
+        super().__init__(test_case=test_case, name=name, subdir=subdir,
+                         ntasks=ntasks, min_tasks=min_tasks,
+                         openmp_threads=openmp_threads)
+        self.add_namelist_file('compass.ocean.tests.turbulence_closure',
+                               'namelist.forward')
+        if resolution < 500:
+            self.add_namelist_file('compass.ocean.tests.turbulence_closure',
+                                   'namelist.les.forward')
+        if resolution < 100:
+            self.add_namelist_file('compass.ocean.tests.turbulence_closure',
+                                   'namelist.nonhydro.forward')
+
+        # make sure output is double precision
+        self.add_streams_file('compass.ocean.streams', 'streams.output')
+
+        self.add_streams_file('compass.ocean.tests.turbulence_closure',
+                              'streams.forward')
+
+        self.add_input_file(filename='init.nc',
+                            target='../initial_state/ocean.nc')
+        self.add_input_file(
+            filename='forcing.nc',
+            target='../initial_state/init_mode_forcing_data.nc')
+        self.add_input_file(filename='graph.info',
+                            target='../initial_state/culled_graph.info')
+
+        self.add_model_as_input()
+
+        self.add_output_file(filename='output.nc')
+
+    # no setup() is needed
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+        # update the time step
+        resolution = self.resolution
+        if resolution == '10km':
+            self.update_namelist_at_runtime({'config_dt':
+                                             "'0000_00:00:01'"})
+        elif resolution == '2m':
+            self.update_namelist_at_runtime({'config_dt':
+                                             "'0000_00:00:00.1'"})
+        elif resolution == '1m':
+            self.update_namelist_at_runtime({'config_dt':
+                                             "'0000_00:00:00.1'"})
+
+        run_model(self)
diff --git a/compass/ocean/tests/turbulence_closure/initial_state.py b/compass/ocean/tests/turbulence_closure/initial_state.py
new file mode 100644
index 0000000000..d02f6da0e5
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/initial_state.py
@@ -0,0 +1,205 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import xarray as xr
+from mpas_tools.io import write_netcdf
+from mpas_tools.mesh.conversion import convert, cull
+from mpas_tools.planar_hex import make_planar_hex_mesh
+
+from compass.ocean.tests.turbulence_closure.boundary_values import (
+    add_boundary_arrays,
+)
+from compass.ocean.vertical import init_vertical_coord
+from compass.step import Step
+
+
+class InitialState(Step):
+    """
+    A step for creating a mesh and initial condition for turbulence closure
+    test cases
+
+    Attributes
+    ----------
+    resolution : str
+        The resolution of the test case
+    """
+    def __init__(self, test_case, resolution):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        test_case : compass.TestCase
+            The test case this step belongs to
+
+        resolution : str
+            The resolution of the test case
+        """
+        super().__init__(test_case=test_case, name='initial_state')
+        self.resolution = resolution
+
+        for file in ['base_mesh.nc', 'culled_graph.info', 'culled_mesh.nc',
+                     'init_mode_forcing_data.nc', 'ocean.nc']:
+            self.add_output_file(file)
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+        config = self.config
+        logger = self.logger
+
+        section = config['turbulence_closure']
+        nx = section.getint('nx')
+        ny = section.getint('ny')
+        dc = section.getfloat('dc')
+
+        x_nonperiodic = section.getboolean('x_nonperiodic')
+        y_nonperiodic = section.getboolean('y_nonperiodic')
+        disturbance_amplitude = section.getfloat('disturbance_amplitude')
+        surface_temperature = section.getfloat('surface_temperature')
+        interior_temperature_shape = section.get('interior_temperature_shape')
+        interior_salinity_shape = section.get('interior_salinity_shape')
+        surface_salinity = section.getfloat('surface_salinity')
+        mixed_layer_salinity_gradient = \
+            section.getfloat('mixed_layer_salinity_gradient')
+        mixed_layer_temperature_gradient = \
+            section.getfloat('mixed_layer_temperature_gradient')
+        mixed_layer_depth_salinity = \
+            section.getfloat('mixed_layer_depth_salinity')
+        mixed_layer_depth_temperature = \
+            section.getfloat('mixed_layer_depth_temperature')
+        interior_salinity_gradient = \
+            section.getfloat('interior_salinity_gradient')
+        interior_temperature_gradient = \
+            section.getfloat('interior_temperature_gradient')
+        interior_temperature_c0 = section.getfloat('interior_temperature_c0')
+        interior_temperature_c1 = section.getfloat('interior_temperature_c1')
+        interior_temperature_c2 = section.getfloat('interior_temperature_c2')
+        interior_temperature_c3 = section.getfloat('interior_temperature_c3')
+        interior_temperature_c4 = section.getfloat('interior_temperature_c4')
+        interior_temperature_c5 = section.getfloat('interior_temperature_c4')
+        interior_salinity_c0 = section.getfloat('interior_salinity_c0')
+        interior_salinity_c1 = section.getfloat('interior_salinity_c1')
+        interior_salinity_c2 = section.getfloat('interior_salinity_c2')
+        interior_salinity_c3 = section.getfloat('interior_salinity_c3')
+        interior_salinity_c4 = section.getfloat('interior_salinity_c4')
+        interior_salinity_c5 = section.getfloat('interior_salinity_c4')
+        surface_heat_flux = section.getfloat('surface_heat_flux')
+        surface_freshwater_flux = section.getfloat('surface_freshwater_flux')
+        wind_stress_zonal = section.getfloat('wind_stress_zonal')
+        wind_stress_meridional = section.getfloat('wind_stress_meridional')
+        coriolis_parameter = section.getfloat('coriolis_parameter')
+        bottom_depth = config.getfloat('vertical_grid', 'bottom_depth')
+
+        dsMesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc,
+                                      nonperiodic_x=x_nonperiodic,
+                                      nonperiodic_y=y_nonperiodic)
+
+        write_netcdf(dsMesh, 'base_mesh.nc')
+
+        dsMesh = cull(dsMesh, logger=logger)
+        dsMesh = convert(dsMesh, graphInfoFileName='culled_graph.info',
+                         logger=logger)
+        write_netcdf(dsMesh, 'culled_mesh.nc')
+
+        ds = dsMesh.copy()
+        dsForcing = dsMesh.copy()
+        xCell = ds.xCell
+
+        ds['bottomDepth'] = bottom_depth * xr.ones_like(xCell)
+        ds['ssh'] = xr.zeros_like(xCell)
+
+        init_vertical_coord(config, ds)
+
+        dsForcing['windStressZonal'] = wind_stress_zonal * xr.ones_like(xCell)
+        dsForcing['windStressMeridional'] = \
+            wind_stress_meridional * xr.ones_like(xCell)
+        dsForcing['sensibleHeatFlux'] = surface_heat_flux * xr.ones_like(xCell)
+        dsForcing['rainFlux'] = surface_freshwater_flux * xr.ones_like(xCell)
+
+        write_netcdf(dsForcing, 'init_mode_forcing_data.nc')
+
+        normalVelocity = xr.zeros_like(ds.xEdge)
+        normalVelocity, _ = xr.broadcast(normalVelocity, ds.refBottomDepth)
+        normalVelocity = normalVelocity.transpose('nEdges', 'nVertLevels')
+
+        temperature = xr.zeros_like(ds.layerThickness)
+
+        salinity = xr.zeros_like(ds.layerThickness)
+
+        zMid = ds.refZMid.values
+
+        temperature_at_mld = (surface_temperature -
+                              mixed_layer_temperature_gradient *
+                              mixed_layer_depth_temperature)
+        if interior_temperature_shape == 'linear':
+            initial_temperature = np.where(
+                zMid >= -mixed_layer_depth_temperature,
+                surface_temperature + mixed_layer_temperature_gradient * zMid,
+                temperature_at_mld + interior_temperature_gradient * zMid)
+        elif interior_temperature_shape == 'polynomial':
+            initial_temperature = np.where(
+                zMid >= -mixed_layer_depth_temperature,
+                surface_temperature + mixed_layer_temperature_gradient * zMid,
+                interior_temperature_c0 +
+                interior_temperature_c1 * zMid +
+                interior_temperature_c2 * zMid**2. +
+                interior_temperature_c3 * zMid**3. +
+                interior_temperature_c4 * zMid**4. +
+                interior_temperature_c5 * zMid**5.)
+        else:
+            print('interior_temperature_shape is not supported')
+
+        salinity_at_mld = (surface_salinity -
+                           mixed_layer_salinity_gradient *
+                           mixed_layer_depth_salinity)
+        if interior_salinity_shape == 'linear':
+            initial_salinity = np.where(
+                zMid >= -mixed_layer_depth_salinity,
+                surface_salinity + mixed_layer_salinity_gradient * zMid,
+                salinity_at_mld + interior_salinity_gradient * zMid)
+        elif interior_salinity_shape == 'polynomial':
+            initial_salinity = np.where(
+                zMid >= -mixed_layer_depth_salinity,
+                surface_salinity + mixed_layer_salinity_gradient * zMid,
+                interior_salinity_c0 +
+                interior_salinity_c1 * zMid +
+                interior_salinity_c2 * zMid**2. +
+                interior_salinity_c3 * zMid**3. +
+                interior_salinity_c4 * zMid**4. +
+                interior_salinity_c5 * zMid**5.)
+        else:
+            print('interior_salinity_shape is not supported')
+
+        normalVelocity = normalVelocity.expand_dims(dim='Time', axis=0)
+        ds['normalVelocity'] = normalVelocity
+
+        temperature[0, :, :] = initial_temperature
+        temperature[0, :, 0] += disturbance_amplitude * 2. * \
+            (np.random.rand(len(ds.xCell.values)) - 0.5)
+        salinity[0, :, :] = initial_salinity
+        ds['temperature'] = temperature
+        ds['salinity'] = salinity
+        ds['fCell'] = coriolis_parameter * xr.ones_like(xCell)
+        ds['fEdge'] = coriolis_parameter * xr.ones_like(ds.xEdge)
+        ds['fVertex'] = coriolis_parameter * xr.ones_like(ds.xVertex)
+
+        ds = add_boundary_arrays(ds, x_nonperiodic, y_nonperiodic)
+
+        write_netcdf(ds, 'ocean.nc')
+
+        plt.figure(dpi=100)
+        plt.plot(initial_temperature, zMid, 'k-')
+        plt.xlabel('PT (C)')
+        plt.ylabel('z (m)')
+        plt.savefig('pt_depth_t0h.png',
+                    bbox_inches='tight', dpi=200)
+        plt.close()
+
+        plt.figure(dpi=100)
+        plt.plot(initial_salinity, zMid, 'k-')
+        plt.xlabel('SA (g/kg)')
+        plt.ylabel('z (m)')
+        plt.savefig('sa_depth_t0h.png',
+                    bbox_inches='tight', dpi=200)
+        plt.close()
diff --git a/compass/ocean/tests/turbulence_closure/namelist.forward b/compass/ocean/tests/turbulence_closure/namelist.forward
new file mode 100644
index 0000000000..d1b9c76235
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/namelist.forward
@@ -0,0 +1,18 @@
+config_write_output_on_startup = .false.
+config_time_integrator = 'split_implicit'
+config_btr_si_partition_match_mode = .false.
+config_run_duration = '0004_00:00:00'
+config_use_implicit_bottom_drag = .false.
+config_zonal_ssh_grad = 0.0
+config_meridional_ssh_grad = 0.0
+config_use_mom_div2 = .true.
+config_mom_div2 = .00002
+config_use_mom_div4 = .false.
+config_mom_div4 = 1.0e12
+config_use_activeTracers_surface_bulk_forcing = .true.
+config_vertMom_del2 = 2e-6
+config_tracer_del2 = 2e-6
+config_mom_del2 = .true.
+config_mom_del2 = 2e-6
+config_cvmix_background_diffusion = 1.0e-2
+config_cvmix_background_viscosity = 1.0e-2
diff --git a/compass/ocean/tests/turbulence_closure/namelist.les.forward b/compass/ocean/tests/turbulence_closure/namelist.les.forward
new file mode 100644
index 0000000000..f26b031187
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/namelist.les.forward
@@ -0,0 +1,13 @@
+config_LES_mode = .true.
+config_use_two_equation_turbulence_model = .true.
+config_two_equation_model_choice = 'omega'
+config_use_cvmix = .false.
+config_two_equation_length_min = 0.1
+config_two_equation_tke_min = 1e-10
+config_two_equation_psi_min = 1e-10
+config_tke_ini = 1.0e-30
+config_LES_noise_enable = .true.
+config_LES_noise_max_time = 150
+config_LES_noise_min_level = 5
+config_LES_noise_max_level = 45
+config_LES_noise_perturbation_magnitude = 1e-4
diff --git a/compass/ocean/tests/turbulence_closure/namelist.nonhydro.forward b/compass/ocean/tests/turbulence_closure/namelist.nonhydro.forward
new file mode 100644
index 0000000000..c5b9b0c098
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/namelist.nonhydro.forward
@@ -0,0 +1,5 @@
+config_enable_nonhydrostatic_mode = .true.
+config_nonhydrostatic_preconditioner = 'bjacobi'
+config_nonhydrostatic_solver_type = 'cg'
+config_nonhydrostatic_solve_surface_boundary_condition = 'pressureTopGradientBottom'
+config_thickness_flux_type = 'upwind'
diff --git a/compass/ocean/tests/turbulence_closure/restart_test/__init__.py b/compass/ocean/tests/turbulence_closure/restart_test/__init__.py
new file mode 100644
index 0000000000..25b1ec4f46
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/restart_test/__init__.py
@@ -0,0 +1,75 @@
+from compass.testcase import TestCase
+from compass.ocean.tests.baroclinic_channel.initial_state import InitialState
+from compass.ocean.tests.baroclinic_channel.forward import Forward
+from compass.ocean.tests import baroclinic_channel
+from compass.validate import compare_variables
+
+
+class RestartTest(TestCase):
+    """
+    A restart test case for the baroclinic channel test group, which makes sure
+    the model produces identical results with one longer run and two shorter
+    runs with a restart in between.
+
+    Attributes
+    ----------
+    resolution : float
+        The resolution of the test case in meters
+    """
+
+    def __init__(self, test_group, resolution, forcing='cooling'):
+        """
+        Create the test case
+
+        Parameters
+        ----------
+        test_group : compass.ocean.tests.baroclinic_channel.BaroclinicChannel
+            The test group that this test case belongs to
+
+        resolution : float
+            The resolution of the test case in meters
+        """
+        name = 'restart_test'
+        self.resolution = resolution
+        if resolution >= 1e3:
+            res_name = f'{int(resolution/1e3)}km'
+        else:
+            res_name = f'{int(resolution)}m'
+        subdir = f'{res_name}/{forcing}/{name}'
+        super().__init__(test_group=test_group, name=name,
+                         subdir=subdir)
+
+        self.add_step(
+            InitialState(test_case=self, resolution=resolution))
+
+        for part in ['full', 'restart']:
+            name = '{}_run'.format(part)
+            step = Forward(test_case=self, name=name, subdir=name, ntasks=4,
+                           openmp_threads=1, resolution=resolution)
+
+            step.add_namelist_file(
+                'compass.ocean.tests.baroclinic_channel.restart_test',
+                'namelist.{}'.format(part))
+            step.add_streams_file(
+                'compass.ocean.tests.baroclinic_channel.restart_test',
+                'streams.{}'.format(part))
+            self.add_step(step)
+
+    def configure(self):
+        """
+        Modify the configuration options for this test case.
+        """
+        baroclinic_channel.configure(self.resolution, self.config)
+
+    # no run() method is needed
+
+    def validate(self):
+        """
+        Test cases can override this method to perform validation of variables
+        and timers
+        """
+        variables = ['temperature', 'salinity', 'layerThickness',
+                     'normalVelocity']
+        compare_variables(test_case=self, variables=variables,
+                          filename1='full_run/output.nc',
+                          filename2='restart_run/output.nc')
diff --git a/compass/ocean/tests/turbulence_closure/restart_test/namelist.full b/compass/ocean/tests/turbulence_closure/restart_test/namelist.full
new file mode 100644
index 0000000000..d99b25960d
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/restart_test/namelist.full
@@ -0,0 +1,3 @@
+config_start_time = '0001-01-01_00:00:00'
+config_run_duration = '0000_00:10:00'
+config_write_output_on_startup = .false.
diff --git a/compass/ocean/tests/turbulence_closure/restart_test/namelist.restart b/compass/ocean/tests/turbulence_closure/restart_test/namelist.restart
new file mode 100644
index 0000000000..620cd94498
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/restart_test/namelist.restart
@@ -0,0 +1,4 @@
+config_start_time = '0001-01-01_00:05:00'
+config_run_duration = '0000_00:05:00'
+config_write_output_on_startup = .false.
+config_do_restart = .true.
diff --git a/compass/ocean/tests/turbulence_closure/restart_test/streams.full b/compass/ocean/tests/turbulence_closure/restart_test/streams.full
new file mode 100644
index 0000000000..5d80d1abb5
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/restart_test/streams.full
@@ -0,0 +1,12 @@
+<streams>
+
+<immutable_stream name="restart"
+                  filename_template="../restarts/rst.$Y-$M-$D_$h.$m.$s.nc"
+                  filename_interval="output_interval"
+                  output_interval="0000-00-00_00:00:01"/>
+
+<stream name="output"
+        output_interval="0000-00-00_00:10:00">
+</stream>
+
+</streams>
diff --git a/compass/ocean/tests/turbulence_closure/restart_test/streams.restart b/compass/ocean/tests/turbulence_closure/restart_test/streams.restart
new file mode 100644
index 0000000000..c99ec09db0
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/restart_test/streams.restart
@@ -0,0 +1,12 @@
+<streams>
+
+<immutable_stream name="restart"
+                  filename_template="../restarts/rst.$Y-$M-$D_$h.$m.$s.nc"
+                  filename_interval="output_interval"
+                  output_interval="0000-00-00_00:00:01"/>
+
+<stream name="output"
+        output_interval="0000-00-00_00:05:00">
+</stream>
+
+</streams>
diff --git a/compass/ocean/tests/turbulence_closure/streams.forward b/compass/ocean/tests/turbulence_closure/streams.forward
new file mode 100644
index 0000000000..dc6a77fe5f
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/streams.forward
@@ -0,0 +1,69 @@
+<streams>
+
+<immutable_stream name="mesh"
+                  filename_template="init.nc"/>
+
+<immutable_stream name="input"
+                  filename_template="init.nc"/>
+
+<immutable_stream name="restart"
+                  output_interval="0000-00-00_00:10:00"/>
+
+<stream name="forcing"
+                  filename_template="forcing.nc"
+                  input_interval="initial_only"
+                  type="input">
+    <var name="sensibleHeatFlux"/>
+</stream>
+
+<stream name="output"
+        type="output"
+        filename_template="output.nc"
+        output_interval="0000_01:00:00"
+        clobber_mode="truncate">
+
+    <stream name="mesh"/>
+    <var_struct name="tracers"/>
+    <var name="xtime"/>
+    <var name="daysSinceStartOfSim"/>
+    <var name="vertDiffTopOfCell"/>
+    <var name="vertViscTopOfCell"/>
+    <var name="velocityZonal"/>
+    <var name="BruntVaisalaFreqTop"/>
+    <var name="tgrad"/>
+    <var name="ugrad"/>
+    <var name="vgrad"/>
+    <var name="Sh"/>
+    <var name="Sm"/>
+    <var name="RiTopOfCell"/>
+    <var name="wgrad"/>
+    <var name="velocityMeridional"/>
+    <var name="normalVelocity"/>
+    <var name="ssh"/>
+    <var name="layerThickness"/>
+    <var name="vertVelocityTop"/>
+    <var name="velocityZonal"/>
+    <var name="velocityMeridional"/>
+    <var name="relativeVorticity"/>
+    <var name="relativeVorticityCell"/>
+    <var name="kineticEnergyCell"/>
+    <var_array name="activeTracerHorMixTendency"/>
+    <var_array name="activeTracerHorizontalAdvectionTendency"/>
+    <var_array name="activeTracerVerticalAdvectionTendency"/>
+    <var_array name="activeTracerVertMixTendency"/>
+    <var name="turbLength"/>
+    <var name="glsPsi"/>
+    <var name="tke"/>
+    <var name="buoyancyProduction"/>
+    <var name="shearProduction"/>
+    <var name="glsAdvTend"/>
+    <var name="glsDissTend"/>
+    <var name="glsTransTend"/>
+    <var name="tkeAdvTend"/>
+    <var name="tkeDissTend"/>
+    <var name="tkeTransTend"/>
+    <var name="verticalVelocity"/>
+    <var name="nonhydrostaticPressure"/>
+</stream>
+
+</streams>
diff --git a/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg b/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg
new file mode 100644
index 0000000000..c755c8422d
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg
@@ -0,0 +1,87 @@
+# Options related to the vertical grid
+[vertical_grid]
+
+# the type of vertical grid
+grid_type = uniform
+
+# The type of vertical coordinate (e.g. z-level, z-star)
+coord_type = z-star
+
+# Whether to use "partial" or "full", or "None" to not alter the topography
+partial_cell_type = None
+
+# The minimum fraction of a layer for partial cells
+min_pc_fraction = 0.1
+
+# config options for baroclinic channel testcases
+[turbulence_closure]
+
+# Logical determining if grid is not periodic in x-direction
+x_nonperiodic = False
+
+# Logical determining if grid is not periodic in y-direction
+y_nonperiodic = False
+
+# Temperature at surface (^oC)
+surface_temperature = 20
+
+# Salinity at surface (PSU)
+surface_salinity = 35
+
+# Depth of mixed layer in temperature (m)
+mixed_layer_depth_temperature = 0
+
+# Depth of mixed layer in salinity (m)
+mixed_layer_depth_salinity = 0
+
+# Temperature gradient in the mixed layer (^oC/m)
+mixed_layer_temperature_gradient = 0
+
+# Salinity gradient in the mixed layer (PSU/m)
+mixed_layer_salinity_gradient = 0
+
+# Interior temperature shape, linear or polynomial
+interior_temperature_shape = linear
+
+# Temperature gradient in the interior (^oC/m)
+interior_temperature_gradient = 0.1
+
+# Temperature polynomial coefficients
+interior_temperature_c0 = 271.65
+interior_temperature_c1 = -0.42
+interior_temperature_c2 = -8.0e-3
+interior_temperature_c3 = -6.75e-5
+interior_temperature_c4 = -2.66e-7
+interior_temperature_c5 = -4.02e-10
+
+# Interior salinity shape, linear or polynomial
+interior_salinity_shape = linear
+
+# Salinity gradient in the interior (PSU/m)
+interior_salinity_gradient = 0
+
+# Salinity polynomial coefficients
+interior_salinity_c0 = 22.36
+interior_salinity_c1 = -0.3
+interior_salinity_c2 = -3.83e-3
+interior_salinity_c3 = -2.54e-5
+interior_salinity_c4 = -8.53e-8
+interior_salinity_c5 = -1.18e-10
+
+# Disturbance amplitude (added to initiate convection faster; ^oC)
+disturbance_amplitude = 0.005
+
+# Surface heat flux (W/m^2)
+surface_heat_flux = -100
+
+# Surface freshwater flux (kg / (m^2 s))
+surface_freshwater_flux = 0
+
+# Surface zonal windstress (N/m^2)
+wind_stress_zonal = 0
+
+# Surface meiridional windstress (N/m^2)
+wind_stress_meridional = 0
+
+# Coriolis parameter for entire domain.
+coriolis_parameter = 1.0e-4
diff --git a/compass/ocean/tests/turbulence_closure/viz.py b/compass/ocean/tests/turbulence_closure/viz.py
new file mode 100644
index 0000000000..405182ad9b
--- /dev/null
+++ b/compass/ocean/tests/turbulence_closure/viz.py
@@ -0,0 +1,195 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import xarray as xr
+
+from compass.step import Step
+
+
+class Viz(Step):
+    """
+    A step for visualizing a cross-section through the domain
+    """
+    def __init__(self, test_case, resolution, forcing, do_comparison=False):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        test_case : compass.TestCase
+            The test case this step belongs to
+        """
+        super().__init__(test_case=test_case, name='viz')
+
+        self.add_input_file(filename='output.nc',
+                            target='../forward/output.nc')
+        self.add_input_file(filename='mesh.nc',
+                            target='../initial_state/culled_mesh.nc')
+
+        if do_comparison:
+            if forcing == 'cooling':
+                forcing_name = 'c02'
+            elif forcing == 'evaporation':
+                forcing_name = 'e04'
+            else:
+                do_comparison = False
+
+        if do_comparison:
+            # Compare all cases with 1m resolution PALM output
+            suffix = 'g128_l128'
+            filename = f'case_{forcing_name}_{suffix}.nc'
+            self.add_input_file(filename='palm.nc', target=filename,
+                                database='turbulence_closure')
+
+        self.do_comparison = do_comparison
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+
+        figsize = [6.4, 4.8]
+        markersize = 5
+
+        dsInit = xr.open_dataset('../forward/init.nc')
+        ds = xr.open_dataset('output.nc')
+        if self.do_comparison:
+            ds_palm = xr.open_dataset('palm.nc')
+
+        if 'Time' not in ds.dims:
+            print('Dataset missing time dimension')
+            return
+
+        nt = ds.sizes['Time']  # number of timesteps
+        tidx = nt - 1
+        seconds_per_day = 24.0 * 3600.0
+        if 'daysSinceStartOfSim' in ds.keys():
+            time = ds.daysSinceStartOfSim
+        else:
+            # This routine is not generic but should not be used as
+            # daysSinceStartOfSim is included in the streams file
+            time0 = 2.0 + (7.0 / 24.0)
+            dt = 0.5 / 24.0
+            time = np.linspace(time0 + dt, time0 + nt * dt, num=nt)
+
+        if self.do_comparison:
+            time_palm = ds_palm.time
+            time_palm_day = (time_palm.astype('float64') /
+                             (seconds_per_day * 1e9))
+            tidx_palm = np.argmin(np.abs(np.subtract(
+                time_palm_day.values, time[tidx])))
+
+        ds = ds.isel(Time=tidx)
+        ds_palm = ds_palm.isel(time=tidx_palm)
+
+        if 'yEdge' not in ds.keys():
+            ds['yEdge'] = dsInit.yEdge
+        ds = ds.sortby('yEdge')
+
+        # Get mesh variables
+        xCell = dsInit.xCell
+        yCell = dsInit.yCell
+        zMid = dsInit.refZMid
+
+        # Import cell quantities
+        temperature = ds.temperature
+        temperature_z = temperature.mean(dim='nCells')
+        salinity = ds.salinity
+        salinity_z = salinity.mean(dim='nCells')
+        w = ds.verticalVelocity
+        w_zTop = w[:, 0]
+
+        velocityZonal = ds.velocityZonal
+        velocityZonal_z = velocityZonal.mean(dim='nCells')
+        velocityMeridional = ds.velocityMeridional
+        velocityMeridional_z = velocityMeridional.mean(dim='nCells')
+        buoyancyProduction = ds.buoyancyProduction
+        buoyancyProduction_z = buoyancyProduction.mean(dim='nCells')
+        wpt = ds.temperatureVerticalAdvectionTendency
+        wpt_z = wpt.mean(dim='nCells')
+
+        if self.do_comparison:
+            alpha_T = ds_palm.alpha_T
+            if 'beta_S' in ds_palm.keys():
+                beta_S = ds_palm.beta_S
+            else:
+                beta_S = 7.8e-4
+            pt_palm = ds_palm.pt - 273.15
+            sa_palm = ds_palm.sa
+            u_palm = ds_palm.u
+            v_palm = ds_palm.v
+            wpt_palm = np.add(ds_palm['w*pt*'].values, ds_palm['w"pt"'].values)
+            wsa_palm = np.add(ds_palm['w*sa*'].values, ds_palm['w"sa"'].values)
+            temp1 = np.multiply(alpha_T, wpt_palm)
+            temp2 = np.multiply(beta_S, wsa_palm)
+            buoyancyProduction_palm = np.multiply(9.81,
+                                                  np.subtract(temp1, temp2))
+            zu_palm = ds_palm.zu
+            z_palm = ds_palm.zprho
+
+        # Figures
+
+        plt.figure(figsize=figsize, dpi=100)
+        plt.plot(temperature_z.values, zMid, 'k-')
+        if self.do_comparison:
+            plt.plot(pt_palm.values, z_palm.values, 'b-')
+        plt.xlabel('PT (C)')
+        plt.ylabel('z (m)')
+        plt.savefig(f'pt_depth_t{int(time[tidx]*24.0)}h.png',
+                    bbox_inches='tight', dpi=200)
+        plt.close()
+
+        plt.figure(figsize=figsize, dpi=100)
+        plt.plot(wpt_z, zMid, 'k-')
+        if self.do_comparison:
+            plt.plot(wpt_palm, z_palm.values, 'b-')
+        plt.xlabel('wpt (C m/s)')
+        plt.ylabel('z (m)')
+        plt.savefig(f'wpt_depth_t{int(time[tidx]*24.0)}h.png',
+                    bbox_inches='tight', dpi=200)
+        plt.close()
+
+        plt.figure(figsize=figsize, dpi=100)
+        plt.plot(buoyancyProduction_z, zMid, 'k-')
+        if self.do_comparison:
+            plt.plot(buoyancyProduction_palm, z_palm.values, 'b-')
+        plt.xlabel('bouyancy production')
+        plt.ylabel('z (m)')
+        plt.savefig(f'buoy_depth_t{int(time[tidx]*24.0)}h.png',
+                    bbox_inches='tight', dpi=200)
+        plt.close()
+
+        plt.figure(figsize=figsize, dpi=100)
+        plt.plot(salinity_z.values, zMid, 'k-')
+        if self.do_comparison:
+            plt.plot(sa_palm.values, z_palm.values, 'b-')
+        plt.xlabel('SA (g/kg)')
+        plt.ylabel('z (m)')
+        plt.savefig(f'sa_depth_t{int(time[tidx]*24.0)}h.png',
+                    bbox_inches='tight', dpi=200)
+        plt.close()
+
+        plt.figure(figsize=figsize, dpi=100)
+        plt.plot(velocityZonal_z.values, zMid, 'k-')
+        plt.plot(velocityMeridional_z.values, zMid, 'k--')
+        if self.do_comparison:
+            plt.plot(u_palm.values, zu_palm.values, 'b-')
+            plt.plot(v_palm.values, zu_palm.values, 'b--')
+        plt.xlabel('u,v (m/s)')
+        plt.ylabel('z (m)')
+        plt.savefig(f'uv_depth_t{int(time[tidx]*24.0)}h.png',
+                    bbox_inches='tight', dpi=200)
+        plt.close()
+
+        plt.figure(figsize=figsize, dpi=100)
+        cmax = np.max(np.abs(w_zTop.values))
+        plt.scatter(np.divide(xCell, 1e3),
+                    np.divide(yCell, 1e3),
+                    s=markersize, c=w_zTop.values,
+                    cmap='cmo.balance', vmin=-1. * cmax, vmax=cmax)
+        plt.xlabel('x (km)')
+        plt.ylabel('y (km)')
+        cbar = plt.colorbar()
+        cbar.ax.set_title('w (m/s)')
+        plt.savefig(f'w_top_section_t{int(time[tidx]*24.0)}h.png',
+                    bbox_inches='tight', dpi=200)
+        plt.close()
diff --git a/docs/users_guide/ocean/test_groups/index.rst b/docs/users_guide/ocean/test_groups/index.rst
index 318d0b92ef..b80488f9f7 100644
--- a/docs/users_guide/ocean/test_groups/index.rst
+++ b/docs/users_guide/ocean/test_groups/index.rst
@@ -29,4 +29,5 @@ coming months.
    sphere_transport
    spherical_harmonic_transform
    tides
-   ziso
\ No newline at end of file
+   turbulence_closure
+   ziso
diff --git a/docs/users_guide/ocean/test_groups/turbulence_closure.rst b/docs/users_guide/ocean/test_groups/turbulence_closure.rst
new file mode 100644
index 0000000000..a4979a5269
--- /dev/null
+++ b/docs/users_guide/ocean/test_groups/turbulence_closure.rst
@@ -0,0 +1,115 @@
+.. _ocean_turbulence_closure:
+
+turbulence_closure
+==================
+
+The ``turbulence_closure`` test group implements studies of the surface layer
+evolution under different kinds of surface forcing for different turbulence
+closures.
+
+The domain is periodic in the horizontal dimensions and no-slip at both
+vertical boundaries.
+
+Initial temperature is 20 degC and salinity is 35 PSU. A linear gradient can be
+given over the mixed layer or the interior. Furthermore, a polynomial function
+can be used to describe the interior temperature and salinity. The default
+parameters are based on a fit to a time average over the whole 1 year
+deployment of ITP #118, located in the Beaufort Gyre (western Arctic).
+
+Variants of the test case are available at 1-m, 2-m and 10-km horizontal
+resolution. The O(1)m resolution test cases are designed to be run in LES mode
+whereas the O(1)km test cases are designed to be used to test turbulence
+closures in "standard" MPAS-Ocean.
+
+The vertical resolution for the LES test cases is the same as the horizontal
+resolution (1m or 2m) whereas the standard resolution test cases have a
+vertical resolution of 50m.
+
+There are currently two surface forcing variants available, cooling and
+evaporation.
+
+All test cases have 2 steps,
+``initial_state``, which defines the mesh and initial conditions for the model,
+and ``forward``, which performs time integration of the model. The ``default``
+test also has a viz step which performs visualization.
+
+ 
+config options
+--------------
+
+namelist options
+----------------
+When the horizontal resolution is less than 500m, LES mode is active with the
+following config options:
+
+.. code-block:: cfg
+
+    config_LES_mode = .true.
+    config_use_two_equation_turbulence_model = .true.
+
+    # TODO Explain what 'omega' means
+    config_two_equation_model_choice = 'omega'
+    config_use_cvmix = .false.
+
+    config_two_equation_length_min = 0.1
+    config_two_equation_tke_min = 1e-10
+    config_two_equation_psi_min = 1e-10
+    config_tke_ini = 1.0e-30
+
+    config_LES_noise_enable = .true.
+    config_LES_noise_max_time = 150
+    config_LES_noise_min_level = 5
+    config_LES_noise_max_level = 45
+
+    # TODO Can you provide any guidance here about what appropriate choices
+    # might be?
+    config_LES_noise_perturbation_magnitude = 1e-4
+
+When the horizontal resolution is less than 100m, nonhydrostatic mode is
+active with the following config options:
+
+.. code-block:: cfg
+
+    # turns the nonhydro model on or off
+    config_enable_nonhydrostatic_mode = .true.  
+
+    # preconditioner for the linear solver. Other options can be used, like
+    # jacobi, sor and asm, but they were found to be slower than block jacobi.
+    config_nonhydrostatic_preconditioner = 'bjacobi'
+
+    # linear solver. Other options can be used, like gmres
+    config_nonhydrostatic_solver_type = 'cg'
+
+    # do not change!
+    config_nonhydrostatic_solve_surface_boundary_condition = 'pressureTopGradientBottom'
+
+    # Nonhydro will work with either 'centered' or 'upwind'.
+    config_thickness_flux_type = 'upwind'
+
+default
+-------
+
+The default version of the turbulence closure test case is available for
+multiple resolutions: ``ocean/baroclinic_channel/1m/default``,
+``ocean/baroclinic_channel/2m/default``, and
+``ocean/baroclinic_channel/10km/default``. The default simulation time is 4 days.
+
+decomp_test
+-----------
+
+``ocean/turbulence_closure/10km/decomp_test`` runs a short (15 min) integration
+of the model forward in time on 4 (``4proc`` step) and then on 8 processors
+(``8proc`` step) to make sure the resulting prognostic variables are
+bit-for-bit identical between the two runs. Currently, only 10-km horizontal
+resolution is supported.
+
+restart_test
+------------
+
+``ocean/turbulence_closure/10km/restart_test`` runs a short (10 min)
+integration of the model forward in time (``full_run`` step), saving a restart
+file every 5 minutes.  Then, a second run (``restart_run`` step) is performed
+from the restart file 5 minutes into the simulation and prognostic variables
+are compared between the "full" and "restart" runs at minute 10 to make sure
+they are bit-for-bit identical. Currently, only 10-km horizontal resolution is
+supported.