diff --git a/docs/developers_guide/ocean/api.md b/docs/developers_guide/ocean/api.md
index 20bb6c0982..96d4ca8623 100644
--- a/docs/developers_guide/ocean/api.md
+++ b/docs/developers_guide/ocean/api.md
@@ -72,10 +72,11 @@
    analysis.Analysis.exact_solution
    analysis.Analysis.run
 
-   forward.compute_max_time_step
    forward.Forward
    forward.Forward.compute_cell_count
+   forward.Forward.compute_max_time_step
    forward.Forward.dynamic_model_config
+   forward.Forward.setup
 
    init.Init
    init.Init.setup
diff --git a/docs/developers_guide/ocean/tasks/barotropic_gyre.md b/docs/developers_guide/ocean/tasks/barotropic_gyre.md
index 6b22472e41..9860ada55f 100644
--- a/docs/developers_guide/ocean/tasks/barotropic_gyre.md
+++ b/docs/developers_guide/ocean/tasks/barotropic_gyre.md
@@ -22,7 +22,7 @@ defines a step for setting up the initial state for each test case.
 First, a mesh appropriate for the resolution is generated using
 {py:func}`mpas_tools.planar_hex.make_planar_hex_mesh()`.  Then, the mesh is
 culled to remove periodicity in the x and y directions.  A vertical grid is
-generated, with 1 layer by default.  Next, the wind stress forcing field is
+generated, with 3 layers by default.  Next, the wind stress forcing field is
 generated.
 
 ### forward
@@ -39,9 +39,10 @@ tasks have as their target and minimum (if the resources are not explicitly
 prescribed).  For MPAS-Ocean, PIO namelist options are modified and a
 graph partition is generated as part of `runtime_setup()`.  Next, the ocean
 model is run. If `run_time_steps` is provided then this determines the run
-duration, otherwise the duration is 3 years. Finally, validation of
-`layerThickness` and `normalVelocity` in the `output.nc` file are performed
-against a baseline if one is provided when calling {ref}`dev-polaris-setup`.
+duration, otherwise the duration is given by the config option `run_duration`.
+Finally, validation of `layerThickness` and `normalVelocity` in the `output.nc`
+file are performed against a baseline if one is provided when calling
+{ref}`dev-polaris-setup`.
 
 ### analysis
 
@@ -51,9 +52,8 @@ the analytical solution for the linearized dynamics. This step also produces a
 figure with the model solution, the analytical solution, and the difference
 between the two.
 
-(dev-ocean-baroclinic-gyre-default)=
+## barotropic_gyre
 
-## default
-
-The {py:class}`polaris.tasks.ocean.baroclinic_channel.default.Default`
-test performs a test of the linearized dynamics.
+The {py:class}`polaris.tasks.ocean.barotropic_gyre.BarotropicGyre`
+test performs either the Munk free-slip or the Munk no-slip cases depending
+on the `test_name` and `boundary_condition` parameters provided.
diff --git a/docs/developers_guide/quick_start.md b/docs/developers_guide/quick_start.md
index d7b6123c40..0970fee55b 100644
--- a/docs/developers_guide/quick_start.md
+++ b/docs/developers_guide/quick_start.md
@@ -493,8 +493,12 @@ Otherwise, to build Omega,
 source ./load_<env_name>_<machine>_<compiler>_<mpi>.sh
 git submodule update --init e3sm_submodules/Omega
 cd e3sm_submodules/Omega
-git submodule update --init --recursive externals/YAKL externals/ekat \
-    externals/scorpio cime
+    externals/YAKL \
+    externals/ekat \
+    externals/scorpio \
+    externals/cpptrace \
+    components/omega/external \
+    cime
 cd components/omega
 mkdir build
 cd build
diff --git a/docs/users_guide/ocean/tasks/barotropic_gyre.md b/docs/users_guide/ocean/tasks/barotropic_gyre.md
index 14a154f1b0..d1a3cb1d86 100644
--- a/docs/users_guide/ocean/tasks/barotropic_gyre.md
+++ b/docs/users_guide/ocean/tasks/barotropic_gyre.md
@@ -23,11 +23,19 @@ clockwise circulation. The comparison with the analytical solution is performed
 via a streamfunction. The L2 norm is printed and the following figure is
 produced:
 
-```{image} images/barotropic_gyre_solution.png
+```{image} images/barotropic_gyre_free-slip_solution.png
 :align: center
 :width: 800 px
 ```
 
+The analytic solution for the free slip condition is derived from
+Henderschott's General Circulation of the Ocean doi:10.1007/978-1-4612-4636-7
+p.223 unnumbered equation after 2.3.9.
+
+The analytic solution for the no slip condition is derived from
+Vallis's Atmospheric and Oceanic Fluid Dynamics doi:10.1017/9781107588417
+p.743 equation 19.49.
+
 ## mesh
 
 The mesh is planar with dimensions given by the config options `lx` and `ly`
@@ -92,11 +100,14 @@ where `tau_0` is given by a config option.
 :width: 500 px
 ```
 
-### config options
+## config options
 
 ```cfg
 [barotropic_gyre]
 
+# time integrator
+time_integrator = split_explicit_ab2
+
 # distance in kilometers between cell centers
 resolution = 20
 
@@ -109,24 +120,21 @@ ly = 1200
 # Maximum amplitude of the zonal wind stress [N m-2]
 tau_0 = 0.1
 
-# Horizontal visocity [m2 s-1]
-nu_2 = 4e2
-
-# [s-1 m-1]
-beta = 10e-11
-
-# [s-1]
-f_0 = 10e-4
+# Horizontal gradient in coriolis parameter [s-1 m-1]
+beta = 1.0e-10
 
 # homogenous fluid density [kg m-3]
 rho_0 = 1000
+
+# Duration of forward step if run_time_steps not given [years]
+run_duration = 2.
 ```
 
 The config option `nu_2` specifies the del2 horizontal viscosity. This value
 will be compared against the resolution to check for stability at set-up. All
 other config options are explained further in previous sections.
 
-### cores
+## cores
 
 The number of cores is determined by `goal_cells_per_core` and
 `max_cells_per_core` in the `ocean` section of the config file.
@@ -135,15 +143,7 @@ The number of cores is determined by `goal_cells_per_core` and
 
 ## supported models
 
-These tasks support only MPAS-Ocean.
-
-## default
-
-The default case is designed such that only a short forward run is performed.
-However, when users set up the default case, a long forward run is also set
-up and may be run along with the analysis step to compare the numerical
-solution with an analytical solution. The user should change the
-`barotropic_gyre: steps_to_run` config option to include desired steps.
+These tasks support both MPAS-Ocean and Omega.
 
 ### time step and run duration
 
@@ -152,5 +152,47 @@ condition assuming a maximum velocity of 1 m/s and a function of the coriolis
 parameter `f_0`. The stability parameter is set to 0.25 based on the stability
 of MPAS-Ocean for this test case.
 
-The run duration is 3 time steps for the `short_forward` step and 3 years for
+The run duration is 3 time steps for the `short_forward` step and 2 years for
 the `long_forward` step.
+
+## munk free-slip
+
+The Munk free-slip test case uses its own set of default horizontal viscosity
+and coriolis parameter values. It also automatically compares against the
+analytic solution for the free-slip case during the `analysis` step.
+
+The case is designed such that only a short forward run is performed.
+However, when users set up the default case, a long forward run and analysis
+step are also set up and may be run together to compare the numerical
+solution with an analytical solution. The user should change the
+`barotropic_gyre_munk_free-slip: steps_to_run` config option to include desired
+steps.
+
+### config options
+
+```cfg
+[barotropic_gyre_munk_free-slip]
+
+# Horizontal visocity [m2 s-1]
+nu_2 = 1e2
+
+# Coriolis parameter [s-1]
+f_0 = 1.0e-4
+```
+
+## munk no-slip
+
+This case has the same attributes described in munk free-slip except the no-
+slip solution and a separate cfg section are used.
+
+### config options
+
+```cfg
+[barotropic_gyre_munk_no-slip]
+
+# Horizontal visocity [m2 s-1]
+nu_2 = 4e2
+
+# Coriolis parameter [s-1]
+f_0 = 1.0e-3
+```
diff --git a/docs/users_guide/ocean/tasks/images/barotropic_gyre_solution.png b/docs/users_guide/ocean/tasks/images/barotropic_gyre_free-slip_solution.png
similarity index 100%
rename from docs/users_guide/ocean/tasks/images/barotropic_gyre_solution.png
rename to docs/users_guide/ocean/tasks/images/barotropic_gyre_free-slip_solution.png
diff --git a/polaris/ocean/model/mpaso_to_omega.yaml b/polaris/ocean/model/mpaso_to_omega.yaml
index 0eea1c772d..23ce54f885 100644
--- a/polaris/ocean/model/mpaso_to_omega.yaml
+++ b/polaris/ocean/model/mpaso_to_omega.yaml
@@ -3,6 +3,7 @@ dimensions:
   nCells: NCells
   nEdges: NEdges
   nVertices: NVertices
+  nTracers: NTracers
   maxEdges: MaxEdges
   maxEdges2: MaxEdges2
   TWO: MaxCellsOnEdge
@@ -11,6 +12,16 @@ dimensions:
   nVertLevelsP1: NVertLevelsP1
 
 variables:
+  # mesh
+  cellsOnCell: CellsOnCell
+  edgesOnCell: EdgesOnCell
+  verticesOnCell: VerticesOnCell
+  cellsOnEdge: CellsOnEdge
+  edgesOnEdge: EdgesOnEdge
+  verticesOnEdge: VerticesOnEdge
+  cellsOnVertex: CellsOnVertex
+  edgesOnVertex: EdgesOnVertex
+
   # tracers
   temperature: Temperature
   salinity: Salinity
@@ -23,7 +34,10 @@ variables:
   normalVelocity: NormalVelocity
 
   # auxiliary state
-  ssh: SshCellDefault
+  ssh: SshCell
+  windStressZonal: WindStressZonal
+  windStressMeridional: WindStressMeridional
+  relativeVorticity: RelVortVertex
 
 config:
 - section:
@@ -67,6 +81,21 @@ config:
     config_use_mom_del4: VelHyperDiffTendencyEnable
     config_mom_del4: ViscDel4
 
+- section:
+    forcing: Tendencies
+  options:
+    config_use_bulk_wind_stress: WindForcingTendencyEnable
+
+- section:
+    debug: Tendencies
+  options:
+    config_disable_vel_explicit_bottom_drag: BottomDragTendencyEnable
+
+- section:
+    bottom_drag: Tendencies
+  options:
+    config_explicit_bottom_drag_coeff: BottomDragCoeff
+
 - section:
     manufactured_solution: ManufacturedSolution
   options:
diff --git a/polaris/suites/ocean/barotropic_gyre.txt b/polaris/suites/ocean/barotropic_gyre.txt
new file mode 100644
index 0000000000..b31537816d
--- /dev/null
+++ b/polaris/suites/ocean/barotropic_gyre.txt
@@ -0,0 +1,2 @@
+ocean/planar/barotropic_gyre/munk/free-slip
+ocean/planar/barotropic_gyre/munk/no-slip
diff --git a/polaris/suites/ocean/omega_pr.txt b/polaris/suites/ocean/omega_pr.txt
index 358c185d21..e60602b168 100644
--- a/polaris/suites/ocean/omega_pr.txt
+++ b/polaris/suites/ocean/omega_pr.txt
@@ -1,4 +1,5 @@
 ocean/planar/manufactured_solution/convergence_both/default
+ocean/planar/barotropic_gyre/munk/free-slip
 ocean/spherical/icos/rotation_2d
 ocean/spherical/icos/cosine_bell/decomp
 ocean/spherical/icos/cosine_bell/restart
diff --git a/polaris/tasks/ocean/barotropic_gyre/__init__.py b/polaris/tasks/ocean/barotropic_gyre/__init__.py
index 871766cf9f..7f08ad60dd 100644
--- a/polaris/tasks/ocean/barotropic_gyre/__init__.py
+++ b/polaris/tasks/ocean/barotropic_gyre/__init__.py
@@ -19,18 +19,42 @@ def add_barotropic_gyre_tasks(component):
     component : polaris.tasks.ocean.Ocean
         the ocean component that the task will be added to
     """
-    test_name = 'default'
-    component.add_task(
-        BarotropicGyre(component=component, test_name=test_name)
+    group_name = 'barotropic_gyre'
+    group_dir = os.path.join('planar', group_name)
+    config_filename = f'{group_name}.cfg'
+    config_filepath = os.path.join(component.name, group_dir, config_filename)
+    config = PolarisConfigParser(filepath=config_filepath)
+    config.add_from_package(
+        f'polaris.tasks.ocean.{group_name}', config_filename
     )
 
+    for boundary_condition in ['free-slip', 'no-slip']:
+        component.add_task(
+            BarotropicGyre(
+                component=component,
+                subdir=group_dir,
+                test_name='munk',
+                boundary_condition=boundary_condition,
+                config=config,
+                config_filename=config_filename,
+            )
+        )
+
 
 class BarotropicGyre(Task):
     """
     The convergence test case for inertial gravity waves
     """
 
-    def __init__(self, component, test_name):
+    def __init__(
+        self,
+        component,
+        subdir,
+        test_name,
+        boundary_condition,
+        config,
+        config_filename,
+    ):
         """
         Create the test case
 
@@ -38,53 +62,62 @@ def __init__(self, component, test_name):
         ----------
         component : polaris.tasks.ocean.Ocean
             The ocean component that this task belongs to
+        subdir : str
+            The subdirectory for the task
+        test_name : str
+            The name of the test (e.g., 'munk')
+        boundary_condition : str
+            The type of boundary condition ('free-slip' or 'no-slip')
+        config : PolarisConfigParser
+            The configuration parser for the task
+        config_filename : str
+            The name of the configuration file
         """
-        group_name = 'barotropic_gyre'
-        name = f'{group_name}_{test_name}'
-        subdir = os.path.join('planar', group_name, test_name)
-        super().__init__(component=component, name=name, subdir=subdir)
-
-        config = self.config
-        config_filename = f'{group_name}.cfg'
-        config.filepath = os.path.join(component.name, subdir, config_filename)
-        config.add_from_package(
-            f'polaris.tasks.ocean.{group_name}', config_filename
-        )
+        name = f'{test_name}/{boundary_condition}'
+        indir = f'{subdir}/{name}'
+        super().__init__(component=component, name=name, subdir=indir)
         self.set_shared_config(config, link=config_filename)
 
-        init = Init(component=component, subdir=subdir)
-        init.set_shared_config(config, link=config_filename)
-        self.add_step(init)
+        init_step = Init(
+            component=component,
+            indir=indir,
+            boundary_condition=boundary_condition,
+            test_name=test_name,
+        )
+        self.add_step(init_step)
 
-        forward = Forward(
+        forward_step = Forward(
             component=component,
-            indir=self.subdir,
+            indir=indir,
             ntasks=None,
             min_tasks=None,
             openmp_threads=1,
+            boundary_condition=boundary_condition,
             name='short_forward',
             run_time_steps=3,
-            graph_target=os.path.join(init.path, 'culled_graph.info'),
+            graph_target=os.path.join(init_step.path, 'culled_graph.info'),
         )
-        forward.set_shared_config(config, link=config_filename)
-        self.add_step(forward)
+        forward_step.set_shared_config(config, link=config_filename)
+        self.add_step(forward_step)
 
         forward = Forward(
             component=component,
-            indir=self.subdir,
+            indir=indir,
             ntasks=None,
             min_tasks=None,
             openmp_threads=1,
+            boundary_condition=boundary_condition,
             name='long_forward',
-            graph_target=os.path.join(init.path, 'culled_graph.info'),
+            graph_target=os.path.join(init_step.path, 'culled_graph.info'),
         )
         forward.set_shared_config(config, link=config_filename)
         self.add_step(forward, run_by_default=False)
 
         analysis = Analysis(
             component=component,
-            indir=self.subdir,
-            boundary_condition='free slip',
+            indir=indir,
+            test_name=test_name,
+            boundary_condition=boundary_condition,
         )
         analysis.set_shared_config(config, link=config_filename)
         self.add_step(analysis, run_by_default=False)
diff --git a/polaris/tasks/ocean/barotropic_gyre/analysis.py b/polaris/tasks/ocean/barotropic_gyre/analysis.py
index 4e0122916e..59827afc3e 100644
--- a/polaris/tasks/ocean/barotropic_gyre/analysis.py
+++ b/polaris/tasks/ocean/barotropic_gyre/analysis.py
@@ -4,7 +4,6 @@
 import matplotlib.pyplot as plt
 import mosaic
 import numpy as np
-import xarray as xr
 from matplotlib import colors as mcolors
 from mpas_tools.ocean import compute_barotropic_streamfunction
 
@@ -15,21 +14,40 @@
 
 class Analysis(OceanIOStep):
     """
-    A step for analysing the output from the barotropic gyre
-    test case
+    A step for analyzing the output from the barotropic gyre test case.
+
+    Attributes
+    ----------
+    boundary_condition : str
+        The type of boundary condition to use ('free-slip' or 'no-slip').
+
+    test_name : str
+        The name of the test case (e.g., 'munk').
     """
 
-    def __init__(self, component, indir, boundary_condition='free slip'):
+    def __init__(
+        self,
+        component,
+        indir,
+        test_name='munk',
+        boundary_condition='free-slip',
+    ):
         """
-        Create the step
+        Create the analysis step.
 
         Parameters
         ----------
         component : polaris.Component
-            The component the step belongs to
+            The component the step belongs to.
 
         indir : str
-            the directory the step is in, to which ``name`` will be appended
+            The directory the step is in, to which ``name`` will be appended.
+
+        test_name : str, optional
+            The name of the test case (default is 'munk').
+
+        boundary_condition : str, optional
+            The type of boundary condition to use (default is 'free-slip').
         """
         super().__init__(component=component, name='analysis', indir=indir)
         self.add_input_file(
@@ -40,15 +58,19 @@ def __init__(self, component, indir, boundary_condition='free slip'):
             filename='output.nc', target='../long_forward/output.nc'
         )
         self.boundary_condition = boundary_condition
+        self.test_name = test_name
 
     def run(self):
-        ds_mesh = xr.open_dataset('mesh.nc')
-        ds_init = xr.open_dataset('init.nc')
-        ds = xr.open_dataset('output.nc')
+        logger = self.logger
+        ds_mesh = self.open_model_dataset('mesh.nc')
+        ds_init = self.open_model_dataset('init.nc')
+        ds = self.open_model_dataset('output.nc')
 
         field_mpas = compute_barotropic_streamfunction(
             ds_init.isel(Time=0), ds, prefix='', time_index=-1
         )
+        x_maxpsi = ds_mesh.xVertex.isel(nVertices=np.argmax(field_mpas.values))
+        logger.info(f'Streamfunction reaches maximum at x = {x_maxpsi.values}')
         field_exact = self.exact_solution(
             ds_mesh,
             self.config,
@@ -66,13 +88,14 @@ def run(self):
             variable_name='psi',
             boundary_condition=self.boundary_condition,
         )
-        print(f'L2 error norm for {self.boundary_condition} bsf: {error:1.2e}')
+        logger.info(
+            f'L2 error norm for {self.boundary_condition} bsf: {error:1.2e}'
+        )
 
         descriptor = mosaic.Descriptor(ds_mesh)
 
         use_mplstyle()
         pad = 20
-        fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(12, 2))
         x0 = ds_mesh.xEdge.min().values
         y0 = ds_mesh.yEdge.min().values
 
@@ -82,6 +105,8 @@ def run(self):
         # convert to km
         descriptor.vertex_patches *= 1.0e-3
 
+        fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(12, 2))
+
         eta0 = max(
             np.max(np.abs(field_exact.values)),
             np.max(np.abs(field_mpas.values)),
@@ -205,15 +230,28 @@ def exact_solution(
         self, ds_mesh, config, loc='Cell', boundary_condition='free slip'
     ):
         """
-        Exact solution to the sea surface height for the linearized Munk layer
-        experiments.
+        Exact solution to the barotropic streamfunction for the linearized Munk
+        layer experiments.
 
         Parameters
         ----------
         ds_mesh : xarray.Dataset
-            Must contain the fields: f'x{loc}', f'y{loc}'
+            The mesh dataset. Must contain the fields: f'x{loc}', f'y{loc}'.
+
+        config : polaris.config.PolarisConfigParser
+            The configuration options for the test case.
+
+        loc : str, optional
+            The location type ('Cell', 'Vertex', etc.) for which to compute
+            the solution.
+
+        boundary_condition : str, optional
+            The type of boundary condition to use ('free-slip' or 'no-slip').
         """
 
+        logger = self.logger
+        test_name = self.test_name
+        boundary_condition = self.boundary_condition
         x = ds_mesh[f'x{loc}']
         x = x - ds_mesh.xEdge.min()
         y = ds_mesh[f'y{loc}']
@@ -224,13 +262,17 @@ def exact_solution(
         # df/dy where f is coriolis parameter
         beta = config.getfloat('barotropic_gyre', 'beta')
         # Laplacian viscosity
-        nu = config.getfloat('barotropic_gyre', 'nu_2')
+        nu = config.getfloat(
+            f'barotropic_gyre_{test_name}_{boundary_condition}', 'nu_2'
+        )
 
         # Compute some non-dimensional numbers
         delta_m = (nu / (beta * L_y**3.0)) ** (1.0 / 3.0)
         gamma = (np.sqrt(3.0) * x) / (2.0 * delta_m * L_x)
+        x_maxpsi = 2 * delta_m / np.sqrt(3.0)
+        logger.info(f'Streamfunction should reach maximum at x = {x_maxpsi}')
 
-        if boundary_condition == 'no slip':
+        if boundary_condition == 'no-slip':
             psi = (
                 pi
                 * np.sin(pi * y / L_y)
@@ -246,15 +288,18 @@ def exact_solution(
                 )
             )
 
-        elif boundary_condition == 'free slip':
+        elif boundary_condition == 'free-slip':
             psi = (
                 pi
-                * (1.0 - x / L_x)
-                * np.sin(pi * y / L_y)
+                * np.sin(pi * (y / L_y))
                 * (
-                    1.0
-                    - np.exp(-x / (2.0 * delta_m * L_x))
-                    * (np.cos(gamma) + (1.0 / np.sqrt(3.0)) * np.sin(gamma))
+                    (1.0 - (x / L_x) - delta_m)
+                    + (np.exp((-(x / L_x)) / (2.0 * delta_m)))
+                    * (
+                        (-2 / 3) * (1 - delta_m) * np.cos(gamma - (pi / 6))
+                        + ((2.0 / np.sqrt(3.0)) * np.sin(gamma))
+                    )
+                    + delta_m * np.exp((((x / L_x) - 1) / delta_m))
                 )
             )
         return psi
diff --git a/polaris/tasks/ocean/barotropic_gyre/barotropic_gyre.cfg b/polaris/tasks/ocean/barotropic_gyre/barotropic_gyre.cfg
index 4d3388ee54..af569cea35 100644
--- a/polaris/tasks/ocean/barotropic_gyre/barotropic_gyre.cfg
+++ b/polaris/tasks/ocean/barotropic_gyre/barotropic_gyre.cfg
@@ -1,28 +1,44 @@
 [barotropic_gyre]
 
+# time integrator
+time_integrator = RK4
+
 # distance in kilometers between cell centers
-resolution = 20
+resolution = 10
 
 # Longituidinal domain length in kilometers
-lx = 1200
+lx = 1000
 
 # Latitudinal domain length in kilometers
-ly = 1200
+ly = 1000
 
 # Maximum amplitude of the zonal wind stress [N m-2]
 tau_0 = 0.1
 
+# Horizontal gradient in coriolis parameter [s-1 m-1]
+beta = 1.0e-10
+
+# homogenous fluid density [kg m-3]
+rho_0 = 1000
+
+# Duration of forward step if run_time_steps not given [years]
+run_duration = 2.
+
+[barotropic_gyre_munk_free-slip]
+
 # Horizontal visocity [m2 s-1]
-nu_2 = 4e2
+nu_2 = 1e2
 
-# [s-1 m-1]
-beta = 10e-11
+# Coriolis parameter [s-1]
+f_0 = 1.0e-4
 
-# [s-1]
-f_0 = 10e-4
+[barotropic_gyre_munk_no-slip]
 
-# homogenous fluid density [kg m-3]
-rho_0 = 1000
+# Horizontal visocity [m2 s-1]
+nu_2 = 4e2
+
+# Coriolis parameter [s-1]
+f_0 = 1.0e-3
 
 [vertical_grid]
 
diff --git a/polaris/tasks/ocean/barotropic_gyre/forward.py b/polaris/tasks/ocean/barotropic_gyre/forward.py
index 53e6e997b1..67a59b91f3 100644
--- a/polaris/tasks/ocean/barotropic_gyre/forward.py
+++ b/polaris/tasks/ocean/barotropic_gyre/forward.py
@@ -1,5 +1,7 @@
 import time
-from math import ceil, floor
+from math import ceil, floor, pi
+
+import numpy as np
 
 from polaris.mesh.planar import compute_planar_hex_nx_ny
 from polaris.ocean.model import OceanModelStep, get_time_interval_string
@@ -12,14 +14,14 @@ class Forward(OceanModelStep):
 
     Attributes
     ----------
-    resolution : float
-        The resolution of the task in km
-
-    dt : float
-        The model time step in seconds
-
     run_time_steps : int or None
         Number of time steps to run for
+
+    test_name : str
+        The name of the test case (e.g., 'munk')
+
+    boundary_condition : str
+        The type of boundary condition ('free-slip' or 'no-slip')
     """
 
     def __init__(
@@ -28,51 +30,61 @@ def __init__(
         name='forward',
         subdir=None,
         indir=None,
+        test_name='munk',
+        boundary_condition='free-slip',
         ntasks=None,
         min_tasks=None,
         openmp_threads=1,
-        nu=None,
         run_time_steps=None,
         graph_target='graph.info',
     ):
         """
-        Create a new task
+        Create a new Forward step for the barotropic gyre task.
 
         Parameters
         ----------
         component : polaris.Component
-            The component the step belongs to
+            The component the step belongs to.
 
-        name : str
-            the name of the task
+        name : str, optional
+            The name of the step. Default is 'forward'.
 
         subdir : str, optional
-            the subdirectory for the step.  If neither this nor ``indir``
-             are provided, the directory is the ``name``
+            The subdirectory for the step. If neither this nor ``indir``
+            are provided, the directory is the ``name``.
 
         indir : str, optional
-            the directory the step is in, to which ``name`` will be appended
+            The directory the step is in, to which ``name`` will be appended.
+
+        test_name : str, optional
+            The name of the test case (e.g., 'munk'). Default is 'munk'.
+
+        boundary_condition : str, optional
+            The type of boundary condition ('free-slip' or 'no-slip').
+            Default is 'free-slip'.
 
         ntasks : int, optional
-            the number of tasks the step would ideally use.  If fewer tasks
+            The number of tasks the step would ideally use. If fewer tasks
             are available on the system, the step will run on all available
-            tasks as long as this is not below ``min_tasks``
+            tasks 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
+            The minimum number of tasks required. 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
+            The number of OpenMP threads the step will use. Default is 1.
 
-        run_time_steps : int or None
-            Number of time steps to run for
+        run_time_steps : int or None, optional
+            Number of time steps to run for. If None, uses config default.
 
         graph_target : str, optional
             The graph file name (relative to the base work directory).
-            If none, it will be created.
+            Default is 'graph.info'.
         """
         self.run_time_steps = run_time_steps
+        self.test_name = test_name
+        self.boundary_condition = boundary_condition
         super().__init__(
             component=component,
             name=name,
@@ -88,7 +100,6 @@ def __init__(
         self.add_yaml_file('polaris.ocean.config', 'output.yaml')
 
         self.add_input_file(filename='init.nc', target='../init/init.nc')
-        self.add_input_file(filename='forcing.nc', target='../init/forcing.nc')
 
         self.add_output_file(
             filename='output.nc',
@@ -131,40 +142,119 @@ def dynamic_model_config(self, at_setup):
         super().dynamic_model_config(at_setup)
 
         config = self.config
+        logger = self.logger
 
         model = config.get('ocean', 'model')
-        if model == 'mpas-ocean':
+        vert_levels = config.getfloat('vertical_grid', 'vert_levels')
+        if model == 'mpas-ocean' and vert_levels == 1:
             self.add_yaml_file('polaris.ocean.config', 'single_layer.yaml')
 
-        nu = config.getfloat('barotropic_gyre', 'nu_2')
+        resolution = config.getfloat('barotropic_gyre', 'resolution')
+        # Laplacian viscosity
+        if self.test_name == 'munk':
+            nu = config.getfloat(
+                f'barotropic_gyre_{self.test_name}_{self.boundary_condition}',
+                'nu_2',
+            )
+        else:
+            nu = 0.0
         rho_0 = config.getfloat('barotropic_gyre', 'rho_0')
+        # beta = df/dy where f is coriolis parameter
+        beta = config.getfloat('barotropic_gyre', 'beta')
+
+        # calculate the boundary layer thickness for specified parameters
+        m = (pi * 2) / np.sqrt(3) * (nu / beta) ** (1.0 / 3.0)
+        # ensure the boundary layer is at least 3 gridcells wide
+        logger.info(
+            'Lateral boundary layer has an anticipated width of '
+            f'{(m * 1e-3):03g} km'
+        )
+        if m <= 3.0 * resolution * 1.0e3:
+            logger.warn(
+                f'Resolution {resolution} km is too coarse to '
+                'properly resolve the lateral boundary layer '
+                f'with anticipated width of {(m * 1e-3):03g} km'
+            )
+
+        # check whether viscosity suitable for stability
+        stability_parameter_max = 0.6
+        dt_max = self.compute_max_time_step(config)
+        nu_max = (
+            stability_parameter_max
+            * (resolution * 1.0e3) ** 2.0
+            / (8 * dt_max)
+        )
+        if nu > nu_max:
+            raise ValueError(
+                f'Laplacian viscosity cannot be set to {nu}; '
+                f'maximum value is {nu_max}'
+            )
 
-        dt_max = compute_max_time_step(config)
-        dt = floor(dt_max / 3.0)
+        dt = floor(dt_max / 5.0)
         dt_str = get_time_interval_string(seconds=dt)
+        dt_btr_str = get_time_interval_string(seconds=dt / 20.0)
+
+        nu_max = stability_parameter_max * (resolution * 1.0e3) ** 2.0 / dt
+        if nu > nu_max:
+            raise ValueError(
+                f'Laplacian viscosity cannot be set to {nu}; '
+                f'maximum value is {nu_max} or decrease the time step'
+            )
 
+        model = config.get('ocean', 'model')
         options = {'config_dt': dt_str, 'config_density0': rho_0}
         self.add_model_config_options(
             options=options, config_model='mpas-ocean'
         )
 
         if self.run_time_steps is not None:
+            output_interval_units = 'Seconds'
             run_duration = ceil(self.run_time_steps * dt)
             stop_time_str = time.strftime(
                 '0001-01-01_%H:%M:%S', time.gmtime(run_duration)
             )
-            output_interval_str = time.strftime(
-                '0000_%H:%M:%S', time.gmtime(run_duration)
-            )
+            if model == 'omega':
+                output_interval_str = str(run_duration)
+            else:
+                output_interval_str = get_time_interval_string(
+                    seconds=run_duration
+                )
         else:
-            stop_time_str = time.strftime('0004-01-01_00:00:00')
-            output_interval_str = time.strftime('0000-01-00_00:00:00')
+            output_interval = 1
+            output_interval_units = 'Months'
+            run_duration = config.getfloat('barotropic_gyre', 'run_duration')
+            stop_time_str = time.strftime(
+                f'{run_duration + 1.0:04g}-01-01_00:00:00'
+            )
+            if model == 'omega':
+                output_interval_str = str(output_interval)
+            else:
+                output_interval_str = get_time_interval_string(
+                    days=output_interval * 30.0
+                )
+
+        slip_factor_dict = {'no-slip': 0.0, 'free-slip': 1.0}
+        time_integrator = config.get('barotropic_gyre', 'time_integrator')
+        time_integrator_map = dict([('RK4', 'RungeKutta4')])
+        if model == 'omega':
+            if time_integrator in time_integrator_map.keys():
+                time_integrator = time_integrator_map[time_integrator]
+            else:
+                print(
+                    'Warning: mapping from time integrator '
+                    f'{time_integrator} to omega not found, '
+                    'retaining name given in config'
+                )
 
         replacements = dict(
             dt=dt_str,
+            dt_btr=dt_btr_str,
             stop_time=stop_time_str,
             output_interval=output_interval_str,
+            output_interval_units=output_interval_units,
+            time_integrator=time_integrator,
             nu=f'{nu:02g}',
+            slip_factor=f'{slip_factor_dict[self.boundary_condition]:02g}',
         )
 
         # make sure output is double precision
@@ -174,27 +264,39 @@ def dynamic_model_config(self, at_setup):
             template_replacements=replacements,
         )
 
+    def setup(self):
+        """
+        TEMP: symlink initial condition to name hard-coded in Omega
+        """
+        super().setup()
+        model = self.config.get('ocean', 'model')
+        # TODO: remove as soon as Omega no longer hard-codes this file
+        if model == 'omega':
+            self.add_input_file(filename='OmegaMesh.nc', target='init.nc')
 
-def compute_max_time_step(config):
-    """
-    Compute the approximate maximum time step for stability
+    def compute_max_time_step(self, config):
+        """
+        Compute the approximate maximum time step for stability
 
-    Parameters
-    ----------
-    config : polaris.config.PolarisConfigParser
-        Config options for test case
+        Parameters
+        ----------
+        config : polaris.config.PolarisConfigParser
+            Config options for test case
 
-    Returns
-    -------
-    dt_max : float
-        The approximate maximum time step for stability
-    """
-    u_max = 1.0  # m/s
-    stability_parameter_max = 0.25
-    resolution = config.getfloat('barotropic_gyre', 'resolution')
-    f_0 = config.getfloat('barotropic_gyre', 'f_0')
-    dt_max = min(
-        stability_parameter_max * resolution * 1e3 / (2 * u_max),
-        stability_parameter_max / f_0,
-    )
-    return dt_max
+        Returns
+        -------
+        dt_max : float
+            The approximate maximum time step for stability
+        """
+        u_max = 1.0  # m/s
+        stability_parameter_max = 0.25
+        resolution = config.getfloat('barotropic_gyre', 'resolution')
+        f_0 = config.getfloat(
+            f'barotropic_gyre_{self.test_name}_{self.boundary_condition}',
+            'f_0',
+        )
+        dt_max = min(
+            stability_parameter_max * resolution * 1e3 / (2 * u_max),
+            stability_parameter_max / f_0,
+        )
+        return dt_max
diff --git a/polaris/tasks/ocean/barotropic_gyre/forward.yaml b/polaris/tasks/ocean/barotropic_gyre/forward.yaml
index e592baca4c..b68985b8b3 100644
--- a/polaris/tasks/ocean/barotropic_gyre/forward.yaml
+++ b/polaris/tasks/ocean/barotropic_gyre/forward.yaml
@@ -4,17 +4,39 @@ ocean:
     config_run_duration: none
   time_integration:
     config_dt: {{ dt }}
-    config_time_integrator: RK4
+    config_time_integrator: {{ time_integrator }}
   hmix_del2:
     config_use_mom_del2: true
     config_mom_del2: {{ nu }}
+  hmix_del4:
+    config_use_mom_del4: false
+  forcing:
+    config_use_bulk_wind_stress: true
+  debug:
+    config_disable_vel_explicit_bottom_drag: true
+
+mpas-ocean:
+  run_modes:
+    config_ocean_run_mode: forward
+  decomposition:
+    config_block_decomp_file_prefix: graph.info.part.
+  advection:
+    config_thickness_flux_type: constant
+  lateral_walls:
+    config_wall_slip_factor: {{ slip_factor }}
+  debug:
+    config_disable_thick_vadv: true
+    config_disable_vel_hadv: true
+    config_disable_vel_hmix: false
+    config_disable_tr_all_tend: true
+
   streams:
     mesh:
       filename_template: init.nc
     input:
       filename_template: init.nc
     forcing:
-      filename_template: forcing.nc
+      filename_template: init.nc
       input_interval: initial_only
       type: input
       contents:
@@ -31,22 +53,41 @@ ocean:
       - normalVelocity
       - layerThickness
       - ssh
-
-mpas-ocean:
-  advection:
-    config_thickness_flux_type: constant
-  bottom_drag:
-    config_implicit_constant_bottom_drag_coeff: 0.
-  forcing:
-    config_use_bulk_wind_stress: true
-  debug:
-    config_disable_vel_hadv: true
-    config_disable_vel_hmix: false
-    config_disable_tr_all_tend: true
+      - relativeVorticity
 
 Omega:
   Tendencies:
-    VelDiffTendencyEnable: true
-    VelHyperDiffTendencyEnable: false
+    TracerHorzAdvTendencyEnable : false
+    TracerDiffTendencyEnable : false
+    TracerHyperDiffTendencyEnable : false
+  Advection:
+    FluxThicknessType: Center
   Dimension:
     NVertLevels: 1
+  IOStreams:
+    InitialState:
+      UsePointerFile: false
+      Filename: init.nc
+      Mode: read
+      Precision: double
+      Freq: 1
+      FreqUnits: OnStartup
+      UseStartEnd: false
+      Contents:
+      - State
+      - WindStressZonal
+      - WindStressMeridional
+    History:
+      UsePointerFile: false
+      Filename: output.nc
+      Mode: write
+      IfExists: append
+      Precision: double
+      Freq: {{ output_interval }}
+      FreqUnits: {{ output_interval_units }}
+      UseStartEnd: false
+      Contents:
+      - State
+      - SshCell
+      FileFreq: 9999
+      FileFreqUnits: years
diff --git a/polaris/tasks/ocean/barotropic_gyre/init.py b/polaris/tasks/ocean/barotropic_gyre/init.py
index 459907fa32..c3ad1d8f7b 100644
--- a/polaris/tasks/ocean/barotropic_gyre/init.py
+++ b/polaris/tasks/ocean/barotropic_gyre/init.py
@@ -1,28 +1,36 @@
 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 polaris import Step
 from polaris.mesh.planar import compute_planar_hex_nx_ny
+from polaris.ocean.model import OceanIOStep
 from polaris.ocean.vertical import init_vertical_coord
-from polaris.tasks.ocean.barotropic_gyre.forward import compute_max_time_step
 from polaris.viz import plot_horiz_field
 
 
-class Init(Step):
+class Init(OceanIOStep):
     """
-    A step for creating a mesh and initial condition for baroclinic channel
-    tasks
+    A step for creating a mesh and initial condition for barotropic gyre tasks
 
     Attributes
     ----------
-    resolution : float
-        The resolution of the task in km
+    name : str
+        The name of the step
+    test_name : str
+        The name of the test case (e.g., 'munk')
+    boundary_condition : str
+        The type of boundary condition (e.g., 'free-slip')
     """
 
-    def __init__(self, component, subdir):
+    def __init__(
+        self,
+        component,
+        indir,
+        name='init',
+        test_name='munk',
+        boundary_condition='free-slip',
+    ):
         """
         Create the step
 
@@ -30,52 +38,31 @@ def __init__(self, component, subdir):
         ----------
         component : polaris.Component
             The component the step belongs to
+        indir : str
+            The input directory for the step
+        name : str, optional
+            The name of the step (default is 'init')
+        test_name : str, optional
+            The name of the test case (default is 'munk')
+        boundary_condition : str, optional
+            The type of boundary condition (default is 'free-slip')
         """
-        super().__init__(component=component, name='init', indir=subdir)
+        super().__init__(component=component, name=name, indir=indir)
 
         for file in [
             'base_mesh.nc',
             'culled_mesh.nc',
             'culled_graph.info',
-            'forcing.nc',
         ]:
             self.add_output_file(file)
+        self.name = name
+        self.test_name = test_name
+        self.boundary_condition = boundary_condition
         self.add_output_file('init.nc', validate_vars=['layerThickness'])
 
     def setup(self):
         super().setup()
 
-        config = self.config
-
-        resolution = config.getfloat('barotropic_gyre', 'resolution')
-        # coriolis parameter
-        beta = config.getfloat('barotropic_gyre', 'beta')
-        # laplacian horizontal viscosity
-        nu_2 = config.getfloat('barotropic_gyre', 'nu_2')
-
-        # calculate the boundary layer thickness for specified parameters
-        m = (np.pi * 2) / np.sqrt(3) * (nu_2 / beta) ** (1.0 / 3.0)
-        # ensure the boundary layer is at least 3 gridcells wide
-        if m <= 3.0 * resolution * 1.0e3:
-            raise ValueError(
-                f'Resolution {resolution} km is too coarse to '
-                'properly resolve the lateral boundary layer '
-                f'with anticipated width of {(m * 1e-3):03g} km'
-            )
-
-        # check whether viscosity suitable for stability
-        stability_parameter_max = 0.6
-        dt_max = compute_max_time_step(config)
-        dt = dt_max / 3.0
-        nu_2_max = (
-            stability_parameter_max * (resolution * 1.0e3) ** 2.0 / (8 * dt)
-        )
-        if nu_2 > nu_2_max:
-            raise ValueError(
-                f'Laplacian viscosity cannot be set to {nu_2}; '
-                f'maximum value is {nu_2_max}'
-            )
-
     def run(self):
         """
         Create the at rest inital condition for the barotropic gyre testcase
@@ -94,18 +81,21 @@ def run(self):
         ds_mesh = make_planar_hex_mesh(
             nx=nx, ny=ny, dc=dc, nonperiodic_x=True, nonperiodic_y=True
         )
-        write_netcdf(ds_mesh, 'base_mesh.nc')
+        self.write_model_dataset(ds_mesh, 'base_mesh.nc')
 
         ds_mesh = cull(ds_mesh, logger=logger)
         ds_mesh = convert(
             ds_mesh, graphInfoFileName='culled_graph.info', logger=logger
         )
-        write_netcdf(ds_mesh, 'culled_mesh.nc')
+        self.write_model_dataset(ds_mesh, 'culled_mesh.nc')
 
         # vertical coordinate parameters
         bottom_depth = config.getfloat('vertical_grid', 'bottom_depth')
         # coriolis parameters
-        f0 = config.getfloat('barotropic_gyre', 'f_0')
+        f0 = config.getfloat(
+            f'barotropic_gyre_{self.test_name}_{self.boundary_condition}',
+            'f_0',
+        )
         beta = config.getfloat('barotropic_gyre', 'beta')
         # surface (wind) forcing parameters
         tau_0 = config.getfloat('barotropic_gyre', 'tau_0')
@@ -135,30 +125,28 @@ def run(self):
         normal_velocity = normal_velocity.expand_dims(dim='Time', axis=0)
         ds['normalVelocity'] = normal_velocity
 
-        # write the initial condition file
-        write_netcdf(ds, 'init.nc')
-
         # set the wind stress forcing
-        ds_forcing = xr.Dataset()
         # Convert from km to m
         ly = ly * 1e3
         wind_stress_zonal = -tau_0 * np.cos(
             np.pi * (ds.yCell - ds.yCell.min()) / ly
         )
         wind_stress_meridional = xr.zeros_like(ds.xCell)
-        ds_forcing['windStressZonal'] = wind_stress_zonal.expand_dims(
+        ds['windStressZonal'] = wind_stress_zonal.expand_dims(
             dim='Time', axis=0
         )
-        ds_forcing['windStressMeridional'] = (
-            wind_stress_meridional.expand_dims(dim='Time', axis=0)
+        ds['windStressMeridional'] = wind_stress_meridional.expand_dims(
+            dim='Time', axis=0
         )
-        write_netcdf(ds_forcing, 'forcing.nc')
+
+        # write the initial condition file
+        self.write_model_dataset(ds, 'init.nc')
 
         cell_mask = ds.maxLevelCell >= 1
 
         plot_horiz_field(
             ds_mesh,
-            ds_forcing['windStressZonal'],
+            ds['windStressZonal'],
             'forcing_wind_stress_zonal.png',
             cmap='cmo.balance',
             show_patch_edges=True,
diff --git a/polaris/tasks/ocean/cosine_bell/forward.yaml b/polaris/tasks/ocean/cosine_bell/forward.yaml
index 042de693fa..0b951925d4 100644
--- a/polaris/tasks/ocean/cosine_bell/forward.yaml
+++ b/polaris/tasks/ocean/cosine_bell/forward.yaml
@@ -79,14 +79,15 @@ Omega:
     TracerHorzAdvTendencyEnable : true
     TracerDiffTendencyEnable : false
     TracerHyperDiffTendencyEnable : false
-  Dimension:
-    NVertLevels: 1
   IOStreams:
+    InitialVertCoord:
+      Filename: init.nc
     InitialState:
       Filename: init.nc
       Contents:
         - State
         - Tracers
+    RestartRead: {}
     History:
       Filename: output.nc
       Freq: {{ output_freq }}
diff --git a/polaris/tasks/ocean/cosine_bell/restart/forward.yaml b/polaris/tasks/ocean/cosine_bell/restart/forward.yaml
index 2547ab6bec..ece4f0f6a4 100644
--- a/polaris/tasks/ocean/cosine_bell/restart/forward.yaml
+++ b/polaris/tasks/ocean/cosine_bell/restart/forward.yaml
@@ -19,6 +19,8 @@ mpas-ocean:
 
 Omega:
   IOStreams:
+    InitialVertCoord:
+      Filename: init.nc
     InitialState:
       FreqUnits: {{ init_freq_units }}
     RestartRead:
diff --git a/polaris/tasks/ocean/manufactured_solution/forward.yaml b/polaris/tasks/ocean/manufactured_solution/forward.yaml
index 025b9eca2c..faf3bff028 100644
--- a/polaris/tasks/ocean/manufactured_solution/forward.yaml
+++ b/polaris/tasks/ocean/manufactured_solution/forward.yaml
@@ -44,15 +44,16 @@ Omega:
     VelHyperDiffTendencyEnable: false
     UseCustomTendency: true
     ManufacturedSolutionTendency: true
-  Dimension:
-    NVertLevels: 1
   IOStreams:
+    InitialVertCoord:
+      Filename: init.nc
     InitialState:
       UsePointerFile: false
       Filename: init.nc
       Contents:
       - NormalVelocity
       - LayerThickness
+    RestartRead: {}
     History:
       Filename: output.nc
       Freq: {{ output_freq }}
@@ -64,4 +65,4 @@ Omega:
       Contents:
       - NormalVelocity
       - LayerThickness
-      - SshCellDefault
+      - SshCell
diff --git a/polaris/tasks/ocean/sphere_transport/forward.yaml b/polaris/tasks/ocean/sphere_transport/forward.yaml
index a1626e6459..24ef51e799 100644
--- a/polaris/tasks/ocean/sphere_transport/forward.yaml
+++ b/polaris/tasks/ocean/sphere_transport/forward.yaml
@@ -71,14 +71,15 @@ Omega:
     TracerHorzAdvTendencyEnable : true
     TracerDiffTendencyEnable : false
     TracerHyperDiffTendencyEnable : false
-  Dimension:
-    NVertLevels: 1
   IOStreams:
+    InitialVertCoord:
+      Filename: init.nc
     InitialState:
       Filename: init.nc
       Contents:
         - State
         - Tracers
+    RestartRead: {}
     History:
       Filename: output.nc
       Freq: {{ output_freq }}