diff --git a/compass/ocean/__init__.py b/compass/ocean/__init__.py
index 7bf9296c18..c9e8619193 100644
--- a/compass/ocean/__init__.py
+++ b/compass/ocean/__init__.py
@@ -1,5 +1,6 @@
from compass.mpas_core import MpasCore
from compass.ocean.tests.baroclinic_channel import BaroclinicChannel
+from compass.ocean.tests.baroclinic_gyre import BaroclinicGyre
from compass.ocean.tests.buttermilk_bay import ButtermilkBay
from compass.ocean.tests.dam_break import DamBreak
from compass.ocean.tests.drying_slope import DryingSlope
@@ -38,6 +39,7 @@ def __init__(self):
super().__init__(name='ocean')
self.add_test_group(BaroclinicChannel(mpas_core=self))
+ self.add_test_group(BaroclinicGyre(mpas_core=self))
self.add_test_group(ButtermilkBay(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/baroclinic_gyre/__init__.py b/compass/ocean/tests/baroclinic_gyre/__init__.py
new file mode 100644
index 0000000000..b39f70d505
--- /dev/null
+++ b/compass/ocean/tests/baroclinic_gyre/__init__.py
@@ -0,0 +1,22 @@
+from compass.ocean.tests.baroclinic_gyre.gyre_test_case import GyreTestCase
+from compass.testgroup import TestGroup
+
+
+class BaroclinicGyre(TestGroup):
+ """
+ A test group for baroclinic gyre 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='baroclinic_gyre')
+
+ for resolution in [20000., 80000.]:
+ self.add_test_case(
+ GyreTestCase(test_group=self, resolution=resolution,
+ long=False))
+ self.add_test_case(
+ GyreTestCase(test_group=self, resolution=resolution,
+ long=True))
diff --git a/compass/ocean/tests/baroclinic_gyre/baroclinic_gyre.cfg b/compass/ocean/tests/baroclinic_gyre/baroclinic_gyre.cfg
new file mode 100644
index 0000000000..757918f28c
--- /dev/null
+++ b/compass/ocean/tests/baroclinic_gyre/baroclinic_gyre.cfg
@@ -0,0 +1,57 @@
+
+# Options related to the vertical grid
+[vertical_grid]
+
+# the type of vertical grid
+grid_type = linear_dz
+
+# the linear rate of thickness (m) increase for linear_dz
+linear_dz_rate = 10.
+
+# Number of vertical levels
+vert_levels = 15
+
+# Total water column depth in m
+bottom_depth = 1800.
+
+# 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 the baroclinic gyre
+[baroclinic_gyre]
+# Basin dimensions
+lat_min = 15
+lat_max = 75
+lon_min = 0
+lon_max = 60
+
+# Initial vertical temperature profile (C)
+initial_temp_top = 33.
+initial_temp_bot = 1.
+
+# Constant salinity value (also used in restoring)
+initial_salinity = 34.
+
+# Maximum zonal wind stress value (N m-2)
+wind_stress_max = 0.1
+
+# Surface temperature restoring profile
+restoring_temp_min = 0.
+restoring_temp_max = 30.
+
+# Restoring timescale for surface temperature (in days)
+restoring_temp_timescale = 30.
+
+# config options for the post processing (moc and viz)
+[baroclinic_gyre_post]
+# latitude bin increment for the moc calculation
+dlat = 0.25
+# number of years to average over for the mean state plots
+time_averaging_length = 1
diff --git a/compass/ocean/tests/baroclinic_gyre/cull_mesh.py b/compass/ocean/tests/baroclinic_gyre/cull_mesh.py
new file mode 100644
index 0000000000..36cccc5f61
--- /dev/null
+++ b/compass/ocean/tests/baroclinic_gyre/cull_mesh.py
@@ -0,0 +1,61 @@
+import numpy as np
+import xarray as xr
+from mpas_tools.io import write_netcdf
+from mpas_tools.mesh.conversion import convert, cull
+
+from compass.step import Step
+
+
+class CullMesh(Step):
+ """
+ Cull a global mesh to only a single basin
+ Attributes
+ ----------
+ """
+ def __init__(self, test_case):
+ """
+ Create the step
+ Parameters
+ ----------
+ test_case : compass.TestCase
+ The test case this step belongs to
+ """
+ super().__init__(test_case=test_case, name='cull_mesh')
+
+ self.add_input_file(
+ filename='base_mesh.nc',
+ target='../base_mesh/base_mesh.nc')
+
+ for file in ['culled_mesh.nc', 'culled_graph.info']:
+ self.add_output_file(file)
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+ config = self.config
+ section = config['baroclinic_gyre']
+ logger = self.logger
+ ds_mesh = xr.open_dataset('base_mesh.nc')
+ ds_mask = xr.Dataset()
+
+ lon = np.rad2deg(ds_mesh.lonCell.values)
+ lat = np.rad2deg(ds_mesh.latCell.values)
+ lon_min = section.getfloat('lon_min')
+ lon_max = section.getfloat('lon_max')
+ lat_min = section.getfloat('lat_min')
+ lat_max = section.getfloat('lat_max')
+
+ mask = np.logical_and(
+ np.logical_and(lon >= lon_min, lon <= lon_max),
+ np.logical_and(lat >= lat_min, lat <= lat_max))
+
+ n_cells = ds_mesh.sizes['nCells']
+ ds_mask['regionCellMasks'] = (('nCells', 'nRegions'),
+ mask.astype(int).reshape(n_cells, 1))
+
+ ds_mesh = cull(ds_mesh, dsInverse=ds_mask, logger=logger)
+ ds_mesh = convert(ds_mesh, graphInfoFileName='culled_graph.info',
+ logger=logger)
+
+ write_netcdf(ds_mesh, 'culled_mesh.nc')
diff --git a/compass/ocean/tests/baroclinic_gyre/forward.py b/compass/ocean/tests/baroclinic_gyre/forward.py
new file mode 100644
index 0000000000..0b31277a3e
--- /dev/null
+++ b/compass/ocean/tests/baroclinic_gyre/forward.py
@@ -0,0 +1,120 @@
+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
+ the baroclinic gyre test cases.
+
+ Attributes
+ ----------
+ resolution : float
+ The resolution of the test case (m)
+
+ """
+ def __init__(self, test_case, resolution, name='forward', subdir=None,
+ long=False):
+ """
+ Create a new test case
+
+ Parameters
+ ----------
+ test_case : compass.TestCase
+ The test case this step belongs to
+
+ resolution : float
+ The resolution of the test case (m)
+
+ name : str
+ the name of the test case
+
+ subdir : str, optional
+ the subdirectory for the step. The default is ``name``
+
+ long : bool, optional
+ Whether to run a long (3-year) simulation to quasi-equilibrium
+ """
+ self.long = long
+ self.resolution = resolution
+ if resolution >= 1e3:
+ res_name = f'{int(resolution/1e3)}km'
+ else:
+ res_name = f'{int(resolution)}m'
+
+ res_params = {'80km': {'ntasks': 32,
+ 'min_tasks': 3,
+ 'dt': "'01:00:00'",
+ 'btr_dt': "'00:03:00'",
+ 'mom_del4': "5.0e11",
+ 'run_duration': "'0000_03:00:00'"},
+ '20km': {'ntasks': 384,
+ 'min_tasks': 32,
+ 'dt': "'00:20:00'",
+ 'btr_dt': "'0000_00:00:20'",
+ 'mom_del4': "2.0e10 ",
+ 'run_duration': "'0000_01:00:00'"}}
+
+ if res_name not in res_params:
+ raise ValueError(
+ f'Unsupported resolution {res_name}. Supported values are: '
+ f'{list(res_params)}')
+
+ res_params = res_params[res_name]
+
+ ntasks = res_params['ntasks']
+ min_tasks = res_params['min_tasks']
+
+ super().__init__(test_case=test_case, name=name, subdir=subdir,
+ ntasks=ntasks, min_tasks=min_tasks, openmp_threads=1)
+
+ # make sure output is double precision
+ self.add_streams_file('compass.ocean.streams', 'streams.output')
+
+ self.add_namelist_file('compass.ocean.tests.baroclinic_gyre',
+ 'namelist.forward')
+ if long:
+ output_interval = "0000-01-00_00:00:00"
+ restart_interval = "0001-00-00_00:00:00"
+ else:
+ output_interval = res_params['run_duration'].replace("'", "")
+ restart_interval = "0030_00:00:00"
+ replacements = dict(
+ output_interval=output_interval, restart_interval=restart_interval)
+ self.add_streams_file(package='compass.ocean.tests.baroclinic_gyre',
+ streams='streams.forward',
+ template_replacements=replacements)
+ options = dict()
+ for option in ['dt', 'btr_dt', 'mom_del4', 'run_duration']:
+ options[f'config_{option}'] = res_params[option]
+ if long:
+ # run for 3 years instead of 3 time steps
+ options['config_run_duration'] = "'0003-00-00_00:00:00'"
+
+ self.add_namelist_options(options=options)
+ self.add_input_file(filename='init.nc',
+ target='../initial_state/initial_state.nc')
+ self.add_input_file(filename='forcing.nc',
+ target='../initial_state/forcing.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/output.0001-01-01_00.00.00.nc')
+ if long:
+ self.add_output_file(filename='output/'
+ 'timeSeriesStatsMonthly.0001-01-01.nc')
+
+ # no setup() is needed
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+ run_model(self, partition_graph=True)
+
+ if self.long:
+ replacements = {'config_do_restart': '.true.',
+ 'config_start_time': "'file'"}
+ self.update_namelist_at_runtime(replacements)
diff --git a/compass/ocean/tests/baroclinic_gyre/gyre_test_case.py b/compass/ocean/tests/baroclinic_gyre/gyre_test_case.py
new file mode 100644
index 0000000000..b8baafab0c
--- /dev/null
+++ b/compass/ocean/tests/baroclinic_gyre/gyre_test_case.py
@@ -0,0 +1,81 @@
+from compass.mesh import QuasiUniformSphericalMeshStep
+from compass.ocean.tests.baroclinic_gyre.cull_mesh import CullMesh
+from compass.ocean.tests.baroclinic_gyre.forward import Forward
+from compass.ocean.tests.baroclinic_gyre.initial_state import InitialState
+from compass.ocean.tests.baroclinic_gyre.moc import Moc
+from compass.ocean.tests.baroclinic_gyre.viz import Viz
+from compass.testcase import TestCase
+from compass.validate import compare_variables
+
+
+class GyreTestCase(TestCase):
+ """
+ A class to define the baroclinic gyre test cases
+
+ Attributes
+ ----------
+ resolution : float
+ The resolution of the test case (m)
+ """
+
+ def __init__(self, test_group, resolution, long):
+ """
+ Create the test case
+
+ Parameters
+ ----------
+ test_group :
+ compass.ocean.tests.baroclinic_gyre.BaroclinicGyre
+ The test group that this test case belongs to
+
+ resolution : float
+ The resolution of the test case (m)
+
+ long : bool
+ Whether to run a long (3-year) simulation to quasi-equilibrium
+ """
+ name = 'performance_test'
+ self.resolution = resolution
+ self.long = long
+
+ if long:
+ name = '3_year_test'
+
+ if resolution >= 1e3:
+ res_name = f'{int(resolution/1e3)}km'
+ else:
+ res_name = f'{int(resolution)}m'
+ subdir = f'{res_name}/{name}'
+ super().__init__(test_group=test_group, name=name,
+ subdir=subdir)
+
+ self.add_step(QuasiUniformSphericalMeshStep(
+ test_case=self, cell_width=int(resolution / 1e3)))
+ self.add_step(CullMesh(test_case=self))
+ self.add_step(
+ InitialState(test_case=self, resolution=resolution))
+ self.add_step(
+ Forward(test_case=self, resolution=resolution,
+ long=long))
+ if long:
+ self.add_step(
+ Moc(test_case=self, resolution=resolution))
+ self.add_step(
+ Viz(test_case=self, resolution=resolution))
+
+ def configure(self):
+ """
+ Set config options for the test case
+ """
+ config = self.config
+ config.add_from_package('compass.mesh', 'mesh.cfg')
+
+ def validate(self):
+ """
+ Validate variables against a baseline
+ """
+ compare_variables(test_case=self,
+ variables=['layerThickness', 'temperature',
+ 'ssh'],
+ filename1='forward/output/'
+ 'output.0001-01-01_00.00.00.nc')
diff --git a/compass/ocean/tests/baroclinic_gyre/initial_state.py b/compass/ocean/tests/baroclinic_gyre/initial_state.py
new file mode 100644
index 0000000000..0b5bdaac1a
--- /dev/null
+++ b/compass/ocean/tests/baroclinic_gyre/initial_state.py
@@ -0,0 +1,145 @@
+import numpy as np
+import xarray as xr
+from mpas_tools.cime.constants import constants
+from mpas_tools.io import write_netcdf
+
+from compass.ocean.vertical import init_vertical_coord
+from compass.step import Step
+
+
+class InitialState(Step):
+ """
+ A step for creating the initial condition for the baroclinic gyre
+ test cases
+
+ Attributes
+ ----------
+ resolution : float
+ The resolution of the test case (m)
+ """
+ def __init__(self, test_case, resolution):
+ """
+ Create the step
+
+ Parameters
+ ----------
+ test_case : compass.TestCase
+ The test case this step belongs to
+
+ resolution : float
+ The resolution of the test case (m)
+ """
+ super().__init__(test_case=test_case, name='initial_state')
+ self.resolution = resolution
+ self.add_input_file(
+ filename='culled_mesh.nc',
+ target='../cull_mesh/culled_mesh.nc')
+ self.add_input_file(
+ filename='culled_graph.info',
+ target='../cull_mesh/culled_graph.info')
+
+ self.add_output_file('initial_state.nc')
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+ config = self.config
+
+ dsMesh = xr.open_dataset('culled_mesh.nc')
+
+ ds = _write_initial_state(config, dsMesh)
+ _write_forcing(config, ds.latCell, ds.refBottomDepth)
+
+
+def _write_initial_state(config, dsMesh):
+
+ ds = dsMesh.copy()
+
+ bottom_depth = config.getfloat('vertical_grid', 'bottom_depth')
+ ds['bottomDepth'] = bottom_depth * xr.ones_like(ds.xCell)
+ ds['ssh'] = xr.zeros_like(ds.xCell)
+
+ init_vertical_coord(config, ds)
+
+ # setting the initial conditions
+ section = config['baroclinic_gyre']
+ initial_salinity = section.getfloat('initial_salinity')
+ temp_top = section.getfloat('initial_temp_top')
+ temp_bot = section.getfloat('initial_temp_bot')
+
+ cc = 0.049 # 0.049 value from optimization on 1800m
+ aa = temp_top / np.log(cc)
+ bb = (cc ** (temp_bot / temp_top) - cc)
+ temperature = aa * np.log(bb * -1.0 * ds.zMid / bottom_depth + cc)
+ val_bot = aa * np.log(bb + cc)
+ val_top = aa * np.log(cc)
+ print(f'analytical bottom T: {val_bot} at '
+ f'depth : {bottom_depth}')
+ print(f'analytical surface T: {val_top} at '
+ f'depth = 0')
+
+ print(f'bottom layer T: {temperature[0, 0, -1]} and '
+ f'surface layer T: {temperature[0, 0, 0]}')
+ salinity = initial_salinity * xr.ones_like(temperature)
+
+ normalVelocity = xr.zeros_like(ds.xEdge)
+ normalVelocity, _ = xr.broadcast(normalVelocity, ds.refBottomDepth)
+ normalVelocity = normalVelocity.transpose('nEdges', 'nVertLevels')
+ normalVelocity = normalVelocity.expand_dims(dim='Time', axis=0)
+
+ ds['temperature'] = temperature
+ ds['salinity'] = salinity
+ ds['normalVelocity'] = normalVelocity
+
+ omega = 2. * np.pi / constants['SHR_CONST_SDAY']
+ ds['fCell'] = 2. * omega * np.sin(ds.latCell)
+ ds['fEdge'] = 2. * omega * np.sin(ds.latEdge)
+ ds['fVertex'] = 2. * omega * np.sin(ds.latVertex)
+
+ write_netcdf(ds, 'initial_state.nc')
+ return ds
+
+
+def _write_forcing(config, lat, refBottomDepth):
+ section = config['baroclinic_gyre']
+ latMin = section.getfloat('lat_min')
+ latMax = section.getfloat('lat_max')
+ tauMax = section.getfloat('wind_stress_max')
+ tempMin = section.getfloat('restoring_temp_min')
+ tempMax = section.getfloat('restoring_temp_max')
+ restoring_temp_timescale = section.getfloat('restoring_temp_timescale')
+ initial_salinity = section.getfloat('initial_salinity')
+ lat = np.rad2deg(lat)
+ # set wind stress
+ windStressZonal = - tauMax * np.cos(2 * np.pi * (lat - latMin) /
+ (latMax - latMin))
+
+ windStressZonal = windStressZonal.expand_dims(dim='Time', axis=0)
+
+ windStressMeridional = xr.zeros_like(windStressZonal)
+
+ # surface restoring
+ temperatureSurfaceRestoringValue = \
+ (tempMax - tempMin) * (latMax - lat) / (latMax - latMin) + tempMin
+ temperatureSurfaceRestoringValue = \
+ temperatureSurfaceRestoringValue.expand_dims(dim='Time', axis=0)
+
+ temperaturePistonVelocity = \
+ (refBottomDepth[0] * xr.ones_like(temperatureSurfaceRestoringValue) /
+ (restoring_temp_timescale * 24. * 3600.))
+
+ salinitySurfaceRestoringValue = \
+ initial_salinity * xr.ones_like(temperatureSurfaceRestoringValue)
+ salinityPistonVelocity = xr.zeros_like(temperaturePistonVelocity)
+
+ dsForcing = xr.Dataset()
+ dsForcing['windStressZonal'] = windStressZonal
+ dsForcing['windStressMeridional'] = windStressMeridional
+ dsForcing['temperaturePistonVelocity'] = temperaturePistonVelocity
+ dsForcing['salinityPistonVelocity'] = salinityPistonVelocity
+ dsForcing['temperatureSurfaceRestoringValue'] = \
+ temperatureSurfaceRestoringValue
+ dsForcing['salinitySurfaceRestoringValue'] = salinitySurfaceRestoringValue
+
+ write_netcdf(dsForcing, 'forcing.nc')
diff --git a/compass/ocean/tests/baroclinic_gyre/moc.py b/compass/ocean/tests/baroclinic_gyre/moc.py
new file mode 100644
index 0000000000..261d4435bd
--- /dev/null
+++ b/compass/ocean/tests/baroclinic_gyre/moc.py
@@ -0,0 +1,94 @@
+import numpy as np
+import xarray
+from mpas_tools.io import write_netcdf
+
+from compass.step import Step
+
+
+class Moc(Step):
+ """
+ A step for computing the zonally-averaged meridional overturning
+ streamfunction in the single basin of the baroclinic gyre
+
+ Attributes
+ ----------
+ resolution : float
+ The horizontal resolution (km) 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 : float
+ The horizontal resolution (km) of the test case
+ """
+ super().__init__(test_case=test_case, name='moc')
+ self.resolution = resolution
+
+ self.add_input_file('../initial_state/initial_state.nc')
+ self.add_input_file(
+ '../forward/output/timeSeriesStatsMonthly.0001-01-01.nc')
+ self.add_output_file('moc.nc')
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+
+ in_dir = '../forward/output'
+ out_dir = '.'
+
+ lat_min = self.config.getfloat('baroclinic_gyre', 'lat_min')
+ lat_max = self.config.getfloat('baroclinic_gyre', 'lat_max')
+ dlat = self.config.getfloat('baroclinic_gyre_post', 'dlat')
+ latBins = np.arange(lat_min + 2 * dlat,
+ lat_max + 2 * dlat, dlat)
+ nz = self.config.getint('vertical_grid', 'vert_levels')
+
+ dsMesh = xarray.open_dataset('../initial_state/initial_state.nc')
+
+ ds = xarray.open_mfdataset(
+ f'{in_dir}/timeSeriesStatsMonthly*.nc',
+ concat_dim='Time', combine='nested')
+
+ moc = self._compute_amoc(dsMesh, ds, latBins, nz)
+
+ dsMOC = xarray.Dataset()
+ dsMOC['xtime_startMonthly'] = ds.xtime_startMonthly
+ dsMOC['xtime_endMonthly'] = ds.xtime_endMonthly
+ dsMOC['moc'] = (["Time", "latBins", "nVertLevelsP1"], moc)
+ dsMOC.coords['latBins'] = latBins
+ dsMOC.coords['nVertLevelsP1'] = np.arange(nz + 1)
+ dsMOC.moc.attrs['units'] = 'Sv'
+ dsMOC.moc.attrs['description'] = \
+ 'zonally-averaged meridional overturning streamfunction'
+
+ outputFileName = f'{out_dir}/moc.nc'
+ write_netcdf(dsMOC, outputFileName)
+
+ def _compute_amoc(self, dsMesh, ds, latBins, nz):
+ """
+ compute the overturning streamfunction for the given mesh
+ """
+
+ latCell = 180. / np.pi * dsMesh.variables['latCell'][:]
+ nt = np.shape(ds.timeMonthly_avg_vertVelocityTop)[0]
+ mocTop = np.zeros([nt, np.size(latBins), nz + 1])
+ indlat_all = [np.logical_and(
+ latCell >= latBins[iLat - 1], latCell < latBins[iLat])
+ for iLat in range(1, np.size(latBins))]
+ for tt in range(nt):
+ for iLat in range(1, np.size(latBins)):
+ indlat = indlat_all[iLat - 1]
+ velArea = (ds.timeMonthly_avg_vertVelocityTop[tt, :, :] *
+ dsMesh.areaCell[:])
+ mocTop[tt, iLat, :] = (mocTop[tt, iLat - 1, :] +
+ np.nansum(velArea[indlat, :], axis=0))
+ # convert m^3/s to Sverdrup
+ mocTop = mocTop * 1e-6
+ return mocTop
diff --git a/compass/ocean/tests/baroclinic_gyre/namelist.forward b/compass/ocean/tests/baroclinic_gyre/namelist.forward
new file mode 100644
index 0000000000..8c37456dbe
--- /dev/null
+++ b/compass/ocean/tests/baroclinic_gyre/namelist.forward
@@ -0,0 +1,25 @@
+config_dt = '00:30:00'
+config_use_mom_del2 = .true.
+config_mom_del2 = 5.0e3
+config_use_tracer_del2 = .true.
+config_tracer_del2 = 1000.0
+config_use_mom_del4 = .false.
+config_use_cvmix_convection = .true.
+config_cvmix_background_diffusion = 1.0e-5
+config_cvmix_background_viscosity = 1.0e-2
+config_use_bulk_wind_stress = .true.
+config_use_activeTracers_surface_restoring = .true.
+config_eos_linear_alpha = 0.2
+config_eos_linear_beta = 0.8
+config_eos_linear_Tref = 5.0
+config_eos_linear_Sref = 35.0
+config_eos_linear_densityref = 1000.0
+config_bottom_drag_mode = 'implicit'
+config_implicit_bottom_drag_type = 'constant'
+config_cvmix_convective_basedOnBVF = .true.
+config_use_cvmix_shear = .true.
+config_cvmix_shear_mixing_scheme = 'KPP'
+config_use_cvmix_kpp = .true.
+config_use_GM = .true.
+config_use_Redi = .true.
+config_AM_timeSeriesStatsMonthly_enable = .true.
diff --git a/compass/ocean/tests/baroclinic_gyre/streams.forward b/compass/ocean/tests/baroclinic_gyre/streams.forward
new file mode 100644
index 0000000000..2f34120ff8
--- /dev/null
+++ b/compass/ocean/tests/baroclinic_gyre/streams.forward
@@ -0,0 +1,113 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/compass/ocean/tests/baroclinic_gyre/viz.py b/compass/ocean/tests/baroclinic_gyre/viz.py
new file mode 100644
index 0000000000..7ae61f4a20
--- /dev/null
+++ b/compass/ocean/tests/baroclinic_gyre/viz.py
@@ -0,0 +1,147 @@
+import cmocean
+import matplotlib.pyplot as plt
+import numpy as np
+import xarray
+from mpas_tools.cime.constants import constants
+
+from compass.step import Step
+
+
+class Viz(Step):
+ """
+ A step for visualizing output from baroclinic gyre
+
+ Attributes
+ ----------
+ resolution : float
+ The horizontal resolution (km) 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 : float
+ The horizontal resolution (km) of the test case
+ """
+ super().__init__(test_case=test_case, name='viz')
+ self.resolution = resolution
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+
+ out_dir = '.'
+ moc_dir = '../moc'
+ mon_dir = '../forward/output'
+ dsMesh = xarray.open_dataset('../initial_state/initial_state.nc')
+ ds = xarray.open_dataset(f'{moc_dir}/moc.nc')
+ # Insert plots here
+ self._plot_moc(ds, dsMesh, out_dir)
+ self._plot_mean_surface_state(mon_dir, dsMesh, out_dir)
+ self._plot_spinup(mon_dir, dsMesh, out_dir)
+
+ def _plot_moc(self, ds, dsMesh, out_dir):
+ """
+ Plot the time-mean moc state for the test case
+ """
+ avg_len = self.config.getint('baroclinic_gyre_post',
+ 'time_averaging_length')
+ moc = ds.moc[-12 * avg_len:, :, :].mean(axis=0).T.values
+ latbins = ds.latBins
+ plt.contourf(
+ latbins, dsMesh.refInterfaces, moc,
+ cmap="RdBu_r", vmin=-12, vmax=12)
+ plt.gca().invert_yaxis()
+ plt.ylabel('Depth (m)')
+ plt.xlabel('Latitude')
+ idx = np.unravel_index(np.argmax(moc), moc.shape)
+ max_moc = round(np.max(moc), 1)
+ max_moc_loc = [latbins[idx[-1]].values,
+ int(dsMesh.refInterfaces[idx[0]].values)]
+ ref_moc_loc = 175 # to be improved
+ max_moc_ref = round(np.max(moc[:, ref_moc_loc]), 1)
+ ref_moc_lat = [latbins[175 - 1].values, latbins[175].values]
+
+ amoc = f"max MOC = {max_moc:.1e}"
+ maxloc = f'at lat = {max_moc_loc[0]} and z = {max_moc_loc[1]}m'
+ maxval = (f'max MOC = {max_moc_ref:.1e} at '
+ f'lat={ref_moc_lat[0]}-{ref_moc_lat[1]}')
+ plt.annotate(amoc + '\n' + maxloc + '\n' + maxval,
+ xy=(0.01, 0.05), xycoords='axes fraction')
+ plt.colorbar()
+ plt.savefig(f'{out_dir}/time_avg_moc_last{avg_len}years.png')
+
+ def _plot_spinup(self, mon_dir, dsMesh, out_dir):
+ """
+ Plot the timeseries of monthy layer-mean
+ kinetic energy and temperature for the test case
+ """
+
+ ds = xarray.open_mfdataset(
+ f'{mon_dir}/timeSeriesStatsMonthly*.nc',
+ concat_dim='Time', combine='nested')
+ KE = ds.timeMonthly_avg_kineticEnergyCell[:, :, :].mean(axis=1)
+ T = ds.timeMonthly_avg_activeTracers_temperature[:, :, :].mean(axis=1)
+ midlayer = (dsMesh.refInterfaces +
+ 0.5 *
+ (np.roll(dsMesh.refInterfaces, -1) - dsMesh.refInterfaces)
+ ).values[:-1]
+
+ fig, ax = plt.subplots(2, 1, figsize=(6, 8))
+ for ll in [0, 3, 6, 10, 14]:
+ ax[0].plot(KE[:, ll], label=f'{int(midlayer[ll])}m')
+ ax[1].plot(T[:, ll], label=f'{int(midlayer[ll])}m')
+ ax[0].legend()
+ ax[1].legend()
+ ax[0].set_xlabel('Time (months)')
+ ax[1].set_xlabel('Time (months)')
+ ax[0].set_ylabel('Layer Mean Kinetic Energy ($m^2 s^{-2}$)')
+ ax[1].set_ylabel(r'Layer Mean Temperature ($^{\circ}$C)')
+ plt.savefig(f'{out_dir}/spinup_ft.png', bbox_inches='tight')
+
+ def _plot_mean_surface_state(self, mon_dir, dsMesh, out_dir):
+
+ lon = 180. / np.pi * dsMesh.variables['lonCell'][:]
+ lat = 180. / np.pi * dsMesh.variables['latCell'][:]
+ ds = xarray.open_mfdataset(
+ f'{mon_dir}/timeSeriesStatsMonthly*.nc',
+ concat_dim='Time', combine='nested')
+ heatflux = (
+ ds.timeMonthly_avg_activeTracersSurfaceFlux_temperatureSurfaceFlux[:, :] * # noqa: E501
+ constants['SHR_CONST_CPSW'] * constants['SHR_CONST_RHOSW'])
+ avg_len = self.config.getint('baroclinic_gyre_post',
+ 'time_averaging_length')
+ absmax = np.max(np.abs(np.mean(heatflux[-12 * avg_len:, :].values, axis=0))) # noqa: E501
+ fig, ax = plt.subplots(1, 3, figsize=[18, 5])
+ ax[0].tricontour(lon, lat, np.mean(ds.timeMonthly_avg_ssh[-12 * avg_len:, :], axis=0), # noqa: E501
+ levels=14, linewidths=0.5, colors='k')
+ ssh = ax[0].tricontourf(lon, lat, np.mean(ds.timeMonthly_avg_ssh[-12 * avg_len:, :], axis=0), # noqa: E501
+ levels=14, cmap="RdBu_r")
+ plt.colorbar(ssh, ax=ax[0])
+ ax[1].tricontour(lon, lat, np.mean(ds.timeMonthly_avg_ssh[-12 * avg_len:, :], axis=0), # noqa: E501
+ levels=14, linewidths=.8, colors='k')
+ temp = ax[1].tricontourf(lon, lat,
+ np.mean(ds.timeMonthly_avg_activeTracers_temperature[-12 * avg_len:, :, 0], axis=0), # noqa: E501
+ levels=15, cmap=cmocean.cm.thermal)
+ plt.colorbar(temp, ax=ax[1])
+ ax[2].tricontour(lon, lat, np.mean(ds.timeMonthly_avg_ssh[-12 * avg_len:, :], axis=0), # noqa: E501
+ levels=14, linewidths=.8, colors='k')
+ hf = ax[2].tricontourf(lon, lat, np.mean(heatflux[-12 * avg_len:, :], axis=0), # noqa: E501
+ levels=np.linspace(- absmax, absmax, 27), cmap="RdBu_r") # noqa: E501
+ plt.colorbar(hf, ax=ax[2])
+
+ ax[0].set_title('SSH (m)')
+ ax[1].set_title(r'SST ($^\circ$C)')
+ ax[2].set_title('Heat Flux (W/m$^{2}$)')
+
+ ax[0].set_ylabel(r'Latitude ($^\circ$)')
+ for axis in ax:
+ axis.set_xlabel(r'Longitude ($^\circ$)')
+ plt.savefig(f'{out_dir}/meansurfacestate_last{avg_len}years.png',
+ bbox_inches='tight')
diff --git a/compass/ocean/vertical/grid_1d/__init__.py b/compass/ocean/vertical/grid_1d/__init__.py
index fab57b5f23..89064fb8ec 100644
--- a/compass/ocean/vertical/grid_1d/__init__.py
+++ b/compass/ocean/vertical/grid_1d/__init__.py
@@ -8,6 +8,7 @@
from compass.ocean.vertical.grid_1d.index_tanh_dz import (
create_index_tanh_dz_grid,
)
+from compass.ocean.vertical.grid_1d.linear_dz import create_linear_dz_grid
from compass.ocean.vertical.grid_1d.tanh_dz import create_tanh_dz_grid
@@ -32,6 +33,13 @@ def generate_1d_grid(config):
if grid_type == 'uniform':
vert_levels = section.getint('vert_levels')
interfaces = _generate_uniform(vert_levels)
+ elif grid_type == 'linear_dz':
+ vert_levels = section.getint('vert_levels')
+ bottom_depth = section.getfloat('bottom_depth')
+ linear_dz_rate = section.getfloat('linear_dz_rate')
+ interfaces = create_linear_dz_grid(vert_levels,
+ bottom_depth,
+ linear_dz_rate)
elif grid_type == 'tanh_dz':
vert_levels = section.getint('vert_levels')
min_layer_thickness = section.getfloat('min_layer_thickness')
diff --git a/compass/ocean/vertical/grid_1d/linear_dz.py b/compass/ocean/vertical/grid_1d/linear_dz.py
new file mode 100644
index 0000000000..465fd8bd9a
--- /dev/null
+++ b/compass/ocean/vertical/grid_1d/linear_dz.py
@@ -0,0 +1,38 @@
+import numpy
+import numpy as np
+
+
+def create_linear_dz_grid(num_vert_levels, bottom_depth,
+ linear_dz_rate):
+ """
+ Creates the linear vertical grid for MPAS-Ocean and
+ writes it to a NetCDF file
+
+ Parameters
+ ----------
+ num_vert_levels : int
+ Number of vertical levels for the grid
+
+ bottom_depth : float
+ bottom depth for the chosen vertical coordinate [m]
+
+ linear_dz_rate : float
+ rate of layer thickness increase (for linear_dz) [m]
+
+ Returns
+ -------
+ interfaces : numpy.ndarray
+ A 1D array of positive depths for layer interfaces in meters
+ """
+
+ nz = num_vert_levels
+ layerThickness = [(bottom_depth / nz) - (np.floor(nz / 2) - k) *
+ linear_dz_rate for k in np.arange(0, nz)]
+ min_layer_thickness = layerThickness[0]
+ max_layer_thickness = layerThickness[-1]
+ print('Linear dz vertical grid')
+ print(f'min layer thickness: {min_layer_thickness}; '
+ f'max layer thickness {max_layer_thickness} in m;')
+ interfaces = - np.append([0], np.cumsum(layerThickness))
+
+ return interfaces
diff --git a/docs/developers_guide/ocean/api.rst b/docs/developers_guide/ocean/api.rst
index 710c81c4fe..bd94c3cac1 100644
--- a/docs/developers_guide/ocean/api.rst
+++ b/docs/developers_guide/ocean/api.rst
@@ -83,6 +83,34 @@ buttermilk_bay
viz.Viz.contour_plots
+baroclinic_gyre
+~~~~~~~~~~~~~~~
+
+.. currentmodule:: compass.ocean.tests.baroclinic_gyre
+
+.. autosummary::
+ :toctree: generated/
+
+ BaroclinicGyre
+
+ GyreTestCase
+ GyreTestCase.configure
+
+ cull_mesh.CullMesh
+ cull_mesh.CullMesh.run
+
+ forward.Forward
+ forward.Forward.run
+
+ initial_state.InitialState
+ initial_state.InitialState.run
+
+ moc.Moc
+ moc.Moc.run
+
+ viz.Viz
+ viz.Viz.run
+
dam_break
~~~~~~~~~
diff --git a/docs/developers_guide/ocean/test_groups/baroclinic_gyre.rst b/docs/developers_guide/ocean/test_groups/baroclinic_gyre.rst
new file mode 100644
index 0000000000..6eaa8c6325
--- /dev/null
+++ b/docs/developers_guide/ocean/test_groups/baroclinic_gyre.rst
@@ -0,0 +1,125 @@
+.. _dev_baroclinic_gyre:
+
+baroclinic_gyre
+===============
+
+The ``baroclinic_gyre`` test group implements variants of the
+Baroclinic ocean gyre set-up from the
+`MITgcm test case `_.
+
+This test case is described in detail in the User guide
+(see :ref:`baroclinic_gyre`). Here, we describe the shared framework
+for this test group.
+
+Framework
+---------
+
+At this stage, the test case is available at 80-km and 20-km horizontal
+resolution. By default, the 15 vertical layers vary linearly in thickness
+with depth, from 50m at the surface to 190m at depth (full depth: 1800m).
+
+The test group includes 2 test cases, called ``performance_test`` for a short
+(3-time-step) run, and ``3_year_test`` for a 3-year simulation. Note that
+3 years are insufficient to bring the standard test case to full equilibrium.
+Both test cases have 2 steps, ``initial_state``, which defines the mesh and
+initial conditions for the model, and ``forward``, which performs the time
+integration of the model.
+
+Additionally, the test group has a shared ``namelist.forward`` file with
+a few common namelist options related to horizontal
+and vertical momentum and tracer diffusion, as well as a shared
+``streams.forward`` file that defines ``mesh``, ``input``, ``restart``, and
+``output`` streams.
+
+cull_mesh
+~~~~~~~~~
+
+The class :py:class:`compass.ocean.tests.baroclinic_gyre.cull_mesh.CullMesh`
+defines a step for setting up the mesh for the baroclinic gyre test case.
+
+First, a mesh appropriate for the resolution is generated using
+:py:class:`compass.mesh.QuasiUniformSphericalMeshStep`. Then, the mesh is
+culled to keep a single ocean basin (with lat, lon bounds set in `.cfg` file).
+
+initial_state
+~~~~~~~~~~~~~
+
+The class :py:class:`compass.ocean.tests.baroclinic_gyre.initial_state.InitialState`
+defines a step for setting up the initial state for each test case.
+
+A vertical grid is generated, with 15 layers of thickness that increases
+linearly with depth (10m increase by default), from 50m at the surface to 190m
+at depth (full depth: 1800m). Finally, the initial temperature field is
+initialized with a vertical profile to approximate the discrete values set in
+the `MITgcm test case `_,
+uniform in horizontal space. The surface and bottom values are set in the
+`.cfg` file. Salinity is set to a constant value (34 psu, set in the `.cfg`
+file) and initial velocity is set to zero.
+
+The ``initial_state`` step also generates the forcing, defined as zonal wind
+stress that varies with latitude, surface temperature restoring that varies
+with latitutde, and writes it to `forcing.nc`.
+
+forward
+~~~~~~~
+
+The class :py:class:`compass.ocean.tests.baroclinic_gyre.forward.Forward`
+defines a step for running MPAS-Ocean from the initial condition produced in
+the ``initial_state`` step.
+
+performance_test
+----------------
+
+``ocean/baroclinic_gyre/80km/performance_test`` is the default version of the
+baroclinic eddies test case for a short (3 time steps) test run and validation of
+prognostic variables for regression testing. Currently, only the 80-km and 20-km horizontal
+resolutions are supported.
+
+3_year_test
+-----------
+
+``ocean/baroclinic_gyre/80km/3_year_test`` performs a longer (3 year) integration
+of the model forward in time. The point is to (ultimately) compare the quasi-
+steady state with theroretical scaling and results from other models.
+Currently, only the 80-km and 20-km horizontal resolutions are supported.
+Note that 3 years is not long enough to reach steady state.
+For the 80km configuration, the estimated time to equilibration is roughly 50 years.
+For a detailed comparison of the mean state against theory and results from other models,
+the mean state at 100 years may be most appropriate to be aligned with the MITgcm results.
+
+moc
+---
+
+The class :py:class:`compass.ocean.tests.baroclinic_gyre.moc.Moc`
+defines a step for computing the zonally-averaged meridional overturning
+in the baroclinic gyre test case.
+
+It calculates the moc by summing the monthly mean vertical velocities
+in latitude bins (the bin width ``dlat`` is set in the ``.cfg`` file).
+It outputs a netcdf file ``moc.nc`` in the local folder, with the monthly
+value of zonally-averaged meriodional overturning in Sv.
+
+By default, this step is only called after the ``3_year_test`` step since
+the ``performance_test`` is too short to output monthly mean output.
+
+viz
+---
+
+The class :py:class:`compass.ocean.tests.baroclinic_gyre.viz.Viz`
+defines a step for visualizing output from baroclinic gyre.
+Currently, it includes 3 separate plots.
+
+First, it plots 2 timeseries to assess the test case spin-up: the top plot is
+layer-mean kinetic energy over 5 different layers (to show dynamical spin-up),
+while the bottom plot is layer mean temperature (to show thermal adjustment).
+The values plotted are all basin-averaged monthly output.
+
+It also plots the time-averaged final state of the simulation. The averaging
+period is set up in the ``.cfg`` file (by default, 1 year). Plots include
+time-mean sea surface height (SSH), sea surface temperature, and surface heat
+flux (in W/m2, due to surface restoring) with contours of SSH superimposed.
+A separate figure is generated showing the time-mean overturning streamfunction,
+based on the ``moc.nc`` file generated in the ``moc`` step described above.
+
+By default, this step is only called after the ``3_year_test`` step since
+the ``performance_test`` is too short to output monthly mean output.
diff --git a/docs/developers_guide/ocean/test_groups/index.rst b/docs/developers_guide/ocean/test_groups/index.rst
index 258bda2911..19fe85b3af 100644
--- a/docs/developers_guide/ocean/test_groups/index.rst
+++ b/docs/developers_guide/ocean/test_groups/index.rst
@@ -8,6 +8,7 @@ Test groups
:titlesonly:
baroclinic_channel
+ baroclinic_gyre
buttermilk_bay
dam_break
drying_slope
diff --git a/docs/users_guide/ocean/test_groups/baroclinic_gyre.rst b/docs/users_guide/ocean/test_groups/baroclinic_gyre.rst
new file mode 100644
index 0000000000..8990555e26
--- /dev/null
+++ b/docs/users_guide/ocean/test_groups/baroclinic_gyre.rst
@@ -0,0 +1,227 @@
+.. _baroclinic_gyre:
+
+baroclinic_gyre
+===============
+
+The ``baroclinic_gyre`` test group implements variants of the
+baroclinic ocean gyre set-up from the `MITgcm test case`_.
+
+This test simulates a baroclinic, wind and buoyancy-forced, double-gyre ocean
+circulation. The grid employs spherical polar coordinates with 15 vertical
+layers. The configuration is similar to the double-gyre setup first solved
+numerically in `Cox and Bryan (1984) `_:
+the model is configured to represent an enclosed sector of fluid on a sphere,
+spanning the tropics to mid-latitudes,
+:math:`60^{\circ} \times 60^{\circ}` in lateral extent.
+The fluid is :math:`1.8`\ km deep and is forced by a zonal wind stress which
+is constant in time, :math:`\tau_{\lambda}`, varying sinusoidally in the
+north-south direction.
+
+.. figure:: ./images/baroclinic_gyre_config.png
+ :width: 95%
+ :align: center
+ :alt: baroclinic gyre configuration
+ :name: baroclinic_gyre_config
+
+ Schematic of simulation domain and wind-stress forcing function for
+ baroclinic gyre numerical experiment. The domain is enclosed by solid walls.
+ From `MITgcm test case`_.
+
+Forcing
+-------
+
+The Coriolis parameter, :math:`f`, is defined
+according to latitude :math:`\varphi`
+
+.. math::
+ f(\varphi) = 2 \Omega \sin( \varphi )
+
+with the rotation rate, :math:`\Omega` set to
+:math:`\frac{2 \pi}{86164} \text{s}^{-1}` (i.e., corresponding to the
+standard Earth rotation rate, using the CIME constant to ensure consistency).
+The sinusoidal wind-stress variations are defined according to
+
+.. math::
+ \tau_{\lambda}(\varphi) = -\tau_{0}\cos \left(2 \pi \frac{\varphi-\varphi_o}{L_{\varphi}} \right)
+
+where :math:`L_{\varphi}` is the lateral domain extent
+(:math:`60^{\circ}`), :math:`\varphi_o` is set to :math:`15^{\circ} \text{N}`
+and :math:`\tau_0` is :math:`0.1 \text{ N m}^{-2}`.
+:ref:`baroclinic_gyre.cfg` summarizes the configuration options used in this
+simulation.
+
+Temperature is restored in the surface layer to a linear profile:
+
+.. math::
+ {\cal F}_\theta = - U_{piston} (\theta-\theta^*), \phantom{WWW}
+ \theta^* = \frac{\theta_{\rm max} - \theta_{\rm min}}{L_\varphi} (\varphi_{\rm max} - \varphi) + \theta_{\rm min}
+ :label: baroc_restore_theta
+
+where the piston velocity :math:`U_{piston}` (in m.s^{-1}) is calculated by
+applying a relaxation timescale of 30 days (set in config file) over the
+thickness of the top layer (50 m by default) and
+:math:`\theta_{\rm max}=30^{\circ}` C, :math:`\theta_{\rm min}=0^{\circ}` C.
+
+Initial state
+-------------
+
+Initially the fluid is stratified with a reference potential temperature
+profile that varies from (approximately) :math:`\theta=30.6 \text{ } ^{\circ}`\ C
+in the surface layer to :math:`\theta=1.56 \text{ } ^{\circ}`\ C in the
+bottom layer. To ensure that the profile is independent of the vertical
+discretization, the profile is now set by a surface value (at the top
+interface) and a bottom value (at the bottom interface), set in the ``.cfg``
+file. The default values have been chosen for the layer values (calculated
+with ``zMid``) to approximate the discrete values presented in the MITgcm
+test case. The temperature functional form (and inner parameter ``cc`` was
+determined by fitting an analytical function to the MITgcm discrete layer
+values (originally ranging from 2 to :math:`30 \text{ } ^{\circ}`\ C. If the
+``bottom_depth`` is different from the default 1800m value, the temperature
+profile is stretched in the vertical to fit the surface and bottom
+temperature constraints, but the thermocline depth and the discrete layer
+values will move away from the MITgcm test case.
+The equation of state used in this experiment is linear:
+
+.. math::
+ \rho = \rho_{0} ( 1 - \alpha_{\theta}\theta^{\prime} )
+ :label: rho_lineareos
+
+with :math:`\rho_{0}=999.8\,{\rm kg\,m}^{-3}` and
+:math:`\alpha_{\theta}=2\times10^{-4}\,{\rm K}^{-1}`. The salinity is set to
+a uniform value of :math:`S=34`\ psu (set in the `.cfg` file)
+Given the linear equation of state, in this configuration the model state
+variable for temperature is equivalent to either in-situ temperature,
+:math:`T`, or potential temperature, :math:`\theta`. For simplicity, here we
+use the variable :math:`\theta` to represent temperature.
+
+Analysis
+--------
+
+For scientific validation, this test case is meant to be run to quasi-steady
+state and its mean state compared to the MITgcm test case and / or
+theoretical scaling. This is done through an analysis step in the
+``3_year_test`` case. Note that 3 years are likely insufficient to bring the
+test case to full equilibrium. Examples of qualitative plots include:
+i) equilibrated SSH contours on top of surface heat fluxes,
+ii) barotropic streamfunction (compared to MITgcm or a barotropic gyre test
+case).
+
+Examples of checks against theory include:
+iii) max of simulated barotropic streamfunction ~ Sverdrup transport,
+iv) simulated thermocline depth ~ scaling argument for penetration depth
+(Vallis (2017) or Cushman-Roisin and Beckers (2011)).
+
+Consider the Sverdrup transport:
+
+.. math:: \rho v_{\rm bt} = \hat{\boldsymbol{k}} \cdot \frac{ \nabla \times \vec{\boldsymbol{\tau}}}{\beta}
+
+If we plug in a typical mid-latitude value for :math:`\beta`
+(:math:`2 \times 10^{-11}` m\ :sup:`-1` s\ :sup:`-1`) and note that
+:math:`\tau` varies by :math:`0.1` Nm\ :sup:`2` over :math:`15^{\circ}`
+latitude, and multiply by the width of our ocean sector, we obtain an
+estimate of approximately 20 Sv.
+
+This scaling is obtained via thermal wind and the linearized barotropic
+vorticity equation), the depth of the thermocline :math:`h` should scale as:
+
+.. math:: h = \left( \frac{w_{\rm Ek} f^2 L_x}{\beta \Delta b} \right) ^2 = \left( \frac{(\tau / L_y) f L_x}{\beta \rho'} \right) ^2
+
+where :math:`w_{\rm Ek}` is a representive value for Ekman pumping,
+:math:`\Delta b = g \rho' / \rho_0`
+is the variation in buoyancy across the gyre,
+and :math:`L_x` and :math:`L_y` are length scales in the
+:math:`x` and :math:`y` directions, respectively.
+Plugging in applicable values at :math:`30^{\circ}`\ N,
+we obtain an estimate for :math:`h` of 200 m.
+
+.. _baroclinic_gyre.cfg:
+
+config options
+--------------
+
+All 2 test cases share the same set of config options:
+
+.. code-block:: cfg
+
+ # Options related to the vertical grid
+ [vertical_grid]
+
+ # the type of vertical grid
+ grid_type = linear_dz
+
+ # the linear rate of thickness (m) increase for linear_dz
+ linear_dz_rate = 10.
+
+ # Number of vertical levels
+ vert_levels = 15
+
+ # Total water column depth in m
+ bottom_depth = 1800.
+
+ # 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 the baroclinic gyre
+ [baroclinic_gyre]
+ # Basin dimensions
+ lat_min = 15
+ lat_max = 75
+ lon_min = 0
+ lon_max = 60
+
+ # Initial vertical temperature profile (C)
+ initial_temp_top = 33.
+ initial_temp_bot = 1.
+
+ # Constant salinity value (also used in restoring)
+ initial_salinity = 34.
+
+ # Maximum zonal wind stress value (N m-2)
+ wind_stress_max = 0.1
+
+ # Surface temperature restoring profile
+ restoring_temp_min = 0.
+ restoring_temp_max = 30.
+
+ # Restoring timescale for surface temperature (in days)
+ restoring_temp_timescale = 30.
+
+
+ # config options for the post processing (moc and viz)
+ [baroclinic_gyre_post]
+ # latitude bin increment for the moc calculation
+ dlat = 0.25
+ # number of years to average over for the mean state plots
+ time_averaging_length = 1
+
+performance_test
+----------------
+
+``ocean/baroclinic_gyre/performance_test`` is the default version of the
+baroclinic_gyre test case for a short (3 time steps) test run and validation
+of prognostic variables for regression testing.
+
+3_year_test
+-----------
+
+``ocean/baroclinic_gyre/3_year_test`` is an additional version of the
+baroclinic_gyre test case with a longer (3-year) spin-up. By default, it
+includes monthly mean output, and plots the mean state of the simulation for
+the last 1 year (option in the config file). Note that for the 80km
+configuration, the estimated time to equilibration is roughly 50 years
+(approx 3 hours of compute time on default layout). This can be done by
+running the ``forward`` step several times (adding 3 years each time), or by
+editing the ``forward/namelist.ocean`` file to set the ``config_run_duration``
+to the desired simulation duration.
+For a detailed comparison of the mean state against theory and results from
+other models, the mean state at 100 years may be most appropriate to be
+aligned with the MITgcm results.
+
+.. _MITgcm test case: https://mitgcm.readthedocs.io/en/latest/examples/baroclinic_gyre/baroclinic_gyre.html
diff --git a/docs/users_guide/ocean/test_groups/images/baroclinic_gyre_config.png b/docs/users_guide/ocean/test_groups/images/baroclinic_gyre_config.png
new file mode 100644
index 0000000000..64d18b0a7c
Binary files /dev/null and b/docs/users_guide/ocean/test_groups/images/baroclinic_gyre_config.png differ
diff --git a/docs/users_guide/ocean/test_groups/index.rst b/docs/users_guide/ocean/test_groups/index.rst
index 56e55dbbcc..60e9c5a7fd 100644
--- a/docs/users_guide/ocean/test_groups/index.rst
+++ b/docs/users_guide/ocean/test_groups/index.rst
@@ -12,6 +12,7 @@ coming months.
:titlesonly:
baroclinic_channel
+ baroclinic_gyre
buttermilk_bay
dam_break
drying_slope