diff --git a/.gitignore b/.gitignore
index 10e0c0d4f6..ecfe483e8b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,9 +8,6 @@ local.config
*.o
*.mod
-# Ignore all PNGs
-*.png
-
# Each tool should also have it's own .gitignore file that ignores the build files for that tool.
# Byte-compiled / optimized / DLL files
diff --git a/compass/ocean/__init__.py b/compass/ocean/__init__.py
index 032b659b10..045e7891b2 100644
--- a/compass/ocean/__init__.py
+++ b/compass/ocean/__init__.py
@@ -3,6 +3,7 @@
from compass.ocean.tests.global_convergence import GlobalConvergence
from compass.ocean.tests.global_ocean import GlobalOcean
from compass.ocean.tests.gotm import Gotm
+from compass.ocean.tests.internal_wave import InternalWave
from compass.ocean.tests.ice_shelf_2d import IceShelf2d
from compass.ocean.tests.isomip_plus import IsomipPlus
from compass.ocean.tests.ziso import Ziso
@@ -23,6 +24,7 @@ def __init__(self):
self.add_test_group(GlobalConvergence(mpas_core=self))
self.add_test_group(GlobalOcean(mpas_core=self))
self.add_test_group(Gotm(mpas_core=self))
+ self.add_test_group(InternalWave(mpas_core=self))
self.add_test_group(IceShelf2d(mpas_core=self))
self.add_test_group(IsomipPlus(mpas_core=self))
self.add_test_group(Ziso(mpas_core=self))
diff --git a/compass/ocean/tests/internal_wave/__init__.py b/compass/ocean/tests/internal_wave/__init__.py
new file mode 100644
index 0000000000..2941460799
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/__init__.py
@@ -0,0 +1,21 @@
+from compass.testgroup import TestGroup
+from compass.ocean.tests.internal_wave.default import Default
+from compass.ocean.tests.internal_wave.rpe_test import RpeTest
+from compass.ocean.tests.internal_wave.ten_day_test import TenDayTest
+
+
+class InternalWave(TestGroup):
+ """
+ A test group for General Ocean Turbulence Model (GOTM) 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='internal_wave')
+
+ self.add_test_case(Default(test_group=self))
+ self.add_test_case(RpeTest(test_group=self))
+ self.add_test_case(TenDayTest(test_group=self))
diff --git a/compass/ocean/tests/internal_wave/default/__init__.py b/compass/ocean/tests/internal_wave/default/__init__.py
new file mode 100644
index 0000000000..6b16442d46
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/default/__init__.py
@@ -0,0 +1,33 @@
+from compass.testcase import TestCase
+from compass.ocean.tests.internal_wave.initial_state import InitialState
+from compass.ocean.tests.internal_wave.forward import Forward
+from compass.ocean.tests.internal_wave.viz import Viz
+from compass.validate import compare_variables
+
+
+class Default(TestCase):
+ """
+ The default test case for the internal wave test
+ """
+
+ def __init__(self, test_group):
+ """
+ Create the test case
+
+ Parameters
+ ----------
+ test_group : compass.ocean.tests.internal_wave.InternalWave
+ The test group that this test case belongs to
+ """
+ super().__init__(test_group=test_group, name='default')
+ self.add_step(InitialState(test_case=self))
+ self.add_step(Forward(test_case=self, cores=4, threads=1))
+ self.add_step(Viz(test_case=self), run_by_default=False)
+
+ def validate(self):
+ """
+ Validate variables against a baseline
+ """
+ compare_variables(test_case=self,
+ variables=['layerThickness', 'normalVelocity'],
+ filename1='forward/output.nc')
diff --git a/compass/ocean/tests/internal_wave/forward.py b/compass/ocean/tests/internal_wave/forward.py
new file mode 100644
index 0000000000..ff09bb94b3
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/forward.py
@@ -0,0 +1,72 @@
+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 internal wave
+ test cases.
+ """
+ def __init__(self, test_case, name='forward', subdir=None, cores=1,
+ min_cores=None, threads=1, nu=None):
+ """
+ Create a new test case
+
+ Parameters
+ ----------
+ test_case : compass.TestCase
+ The test case this step belongs to
+
+ name : str
+ the name of the test case
+
+ subdir : str, optional
+ the subdirectory for the step. The default is ``name``
+
+ cores : int, optional
+ the number of cores the step would ideally use. If fewer cores
+ are available on the system, the step will run on all available
+ cores as long as this is not below ``min_cores``
+
+ min_cores : int, optional
+ the number of cores the step requires. If the system has fewer
+ than this number of cores, the step will fail
+
+ threads : int, optional
+ the number of threads the step will use
+
+ nu : float, optional
+ the viscosity (if different from the default for the test group)
+ """
+ if min_cores is None:
+ min_cores = cores
+ super().__init__(test_case=test_case, name=name, subdir=subdir,
+ cores=cores, min_cores=min_cores, threads=threads)
+ self.add_namelist_file('compass.ocean.tests.internal_wave',
+ 'namelist.forward')
+ if nu is not None:
+ # update the viscosity to the requested value
+ options = {'config_mom_del2': '{}'.format(nu)}
+ self.add_namelist_options(options)
+
+ self.add_streams_file('compass.ocean.tests.internal_wave',
+ 'streams.forward')
+
+ self.add_input_file(filename='init.nc',
+ target='../initial_state/ocean.nc')
+ self.add_input_file(filename='mesh.nc',
+ target='../initial_state/culled_mesh.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
+ """
+ run_model(self)
diff --git a/compass/ocean/tests/internal_wave/initial_state.py b/compass/ocean/tests/internal_wave/initial_state.py
new file mode 100644
index 0000000000..018756dd2d
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/initial_state.py
@@ -0,0 +1,132 @@
+import xarray
+import numpy
+
+from mpas_tools.planar_hex import make_planar_hex_mesh
+from mpas_tools.io import write_netcdf
+from mpas_tools.mesh.conversion import convert, cull
+
+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 internal wave test
+ cases
+ """
+ def __init__(self, test_case):
+ """
+ Create the step
+
+ Parameters
+ ----------
+ test_case : compass.ocean.tests.internal_wave.default.Default
+ The test case this step belongs to
+ """
+ super().__init__(test_case=test_case, name='initial_state', cores=1,
+ min_cores=1, threads=1)
+
+ self.add_namelist_file('compass.ocean.tests.internal_wave',
+ 'namelist.init')
+
+ self.add_streams_file('compass.ocean.tests.internal_wave',
+ 'streams.init')
+
+ for file in ['base_mesh.nc', 'culled_mesh.nc', 'culled_graph.info',
+ 'ocean.nc']:
+ self.add_output_file(file)
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+ config = self.config
+ logger = self.logger
+
+ replacements = dict()
+ replacements['config_periodic_planar_vert_levels'] = \
+ config.getfloat('vertical_grid', 'vert_levels')
+ replacements['config_periodic_planar_bottom_depth'] = \
+ config.getfloat('vertical_grid', 'bottom_depth')
+ self.update_namelist_at_runtime(options=replacements)
+
+ section = config['vertical_grid']
+ vert_levels = section.getint('vert_levels')
+ bottom_depth = section.getfloat('bottom_depth')
+
+ section = config['internal_wave']
+ nx = section.getint('nx')
+ ny = section.getint('ny')
+ dc = section.getfloat('dc')
+ use_distances = section.getboolean('use_distances')
+ amplitude_width_dist = section.getfloat('amplitude_width_dist')
+ amplitude_width_frac = section.getfloat('amplitude_width_frac')
+ bottom_temperature = section.getfloat('bottom_temperature')
+ surface_temperature = section.getfloat('surface_temperature')
+ temperature_difference = section.getfloat('temperature_difference')
+ salinity = section.getfloat('salinity')
+
+ logger.info(' * Make planar hex mesh')
+ dsMesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc, nonperiodic_x=False,
+ nonperiodic_y=True)
+ logger.info(' * Completed Make planar hex mesh')
+ write_netcdf(dsMesh, 'base_mesh.nc')
+
+ logger.info(' * Cull mesh')
+ dsMesh = cull(dsMesh, logger=logger)
+ logger.info(' * Convert mesh')
+ dsMesh = convert(dsMesh, graphInfoFileName='culled_graph.info',
+ logger=logger)
+ logger.info(' * Completed Convert mesh')
+ write_netcdf(dsMesh, 'culled_mesh.nc')
+
+ ds = dsMesh.copy()
+ yCell = ds.yCell
+
+ ds['bottomDepth'] = bottom_depth * xarray.ones_like(yCell)
+ ds['ssh'] = xarray.zeros_like(yCell)
+
+ init_vertical_coord(config, ds)
+
+ yMin = yCell.min().values
+ yMax = yCell.max().values
+
+ yMid = 0.5*(yMin + yMax)
+
+ if use_distances:
+ perturbation_width = amplitude_width_dist
+ else:
+ perturbation_width = (yMax - yMin) * amplitude_width_frac
+
+ # Set stratified temperature
+ temp_vert = (bottom_temperature
+ + (surface_temperature - bottom_temperature) *
+ ((ds.refZMid + bottom_depth) / bottom_depth))
+
+ depth_frac = xarray.zeros_like(temp_vert)
+ refBottomDepth = ds['refBottomDepth']
+ for k in range(1, vert_levels):
+ depth_frac[k] = refBottomDepth[k-1] / refBottomDepth[vert_levels-1]
+
+ # If cell is in the southern half, outside the sin width, subtract
+ # temperature difference
+ frac = xarray.where(numpy.abs(yCell - yMid) < perturbation_width,
+ numpy.cos(0.5 * numpy.pi * (yCell - yMid) /
+ perturbation_width) *
+ numpy.sin(numpy.pi * depth_frac),
+ 0.)
+
+ temperature = temp_vert - temperature_difference * frac
+ temperature = temperature.transpose('nCells', 'nVertLevels')
+ temperature = temperature.expand_dims(dim='Time', axis=0)
+
+ normalVelocity = xarray.zeros_like(ds.xEdge)
+ normalVelocity, _ = xarray.broadcast(normalVelocity, ds.refBottomDepth)
+ normalVelocity = normalVelocity.transpose('nEdges', 'nVertLevels')
+ normalVelocity = normalVelocity.expand_dims(dim='Time', axis=0)
+
+ ds['temperature'] = temperature
+ ds['salinity'] = salinity * xarray.ones_like(temperature)
+ ds['normalVelocity'] = normalVelocity
+
+ write_netcdf(ds, 'ocean.nc')
diff --git a/compass/ocean/tests/internal_wave/internal_wave.cfg b/compass/ocean/tests/internal_wave/internal_wave.cfg
new file mode 100644
index 0000000000..fe73eda65c
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/internal_wave.cfg
@@ -0,0 +1,56 @@
+# Options related to the vertical grid
+[vertical_grid]
+
+# the type of vertical grid
+grid_type = uniform
+
+# Number of vertical levels
+vert_levels = 20
+
+# Depth of the bottom of the ocean
+bottom_depth = 500.0
+
+# The type of vertical coordinate (e.g. z-level, z-star)
+coord_type = z-level
+
+# 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 internal wave test cases
+[internal_wave]
+
+# the number of grid cells in x and y
+nx = 4
+ny = 50
+
+# the size of grid cells (m)
+dc = 5000.0
+
+# Logical flag that determines if locations of features are defined by distance
+# or fractions. False means fractions.
+use_distances = False
+
+# Temperature of the surface in the northern half of the domain.
+surface_temperature = 20.1
+
+# Temperature of the bottom in the northern half of the domain.
+bottom_temperature = 10.1
+
+# Difference in the temperature field between top and bottom
+temperature_difference = 2.0
+
+# Fraction of domain in Y direction the temperature gradient should be linear
+# over.
+amplitude_width_frac = 0.33
+
+# Width of the temperature gradient around the center sin wave. Default value
+# is relative to a 500km domain in Y.
+amplitude_width_dist = 50e3
+
+# Salinity of the water in the entire domain.
+salinity = 35.0
+
+# isopycnal displacement = 125.0 Isopycnal configuration not implemented
diff --git a/compass/ocean/tests/internal_wave/namelist.forward b/compass/ocean/tests/internal_wave/namelist.forward
new file mode 100644
index 0000000000..5b4c080ab2
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/namelist.forward
@@ -0,0 +1,10 @@
+config_dt = '00:05:00'
+config_btr_dt = '00:00:15'
+config_time_integrator = 'split_explicit'
+config_run_duration = '0000_00:15:00'
+config_use_cvmix = .true.
+config_pio_stride = 4
+config_use_mom_del2 = .true.
+config_mom_del2 = 10.0
+config_implicit_bottom_drag_coeff = 1.0e-2
+config_use_cvmix_convection = .true.
diff --git a/compass/ocean/tests/internal_wave/namelist.init b/compass/ocean/tests/internal_wave/namelist.init
new file mode 100644
index 0000000000..1097b37834
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/namelist.init
@@ -0,0 +1,3 @@
+config_init_configuration = 'internal_waves'
+config_vert_levels = 1
+config_write_cull_cell_mask = .true.
diff --git a/compass/ocean/tests/internal_wave/rpe_test/__init__.py b/compass/ocean/tests/internal_wave/rpe_test/__init__.py
new file mode 100644
index 0000000000..a4e5287d82
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/rpe_test/__init__.py
@@ -0,0 +1,51 @@
+from compass.testcase import TestCase
+from compass.ocean.tests.internal_wave.initial_state import InitialState
+from compass.ocean.tests.internal_wave.forward import Forward
+from compass.ocean.tests.internal_wave.rpe_test.analysis import Analysis
+from compass.ocean.tests import internal_wave
+
+
+class RpeTest(TestCase):
+ """
+ The reference potential energy (RPE) test case for the internal wave
+ test group performs a 20-day integration of the model forward in time at
+ 5 different values of the viscosity at the given resolution.
+
+ Attributes
+ ----------
+ resolution : str
+ The resolution of the test case
+ """
+
+ def __init__(self, test_group):
+ """
+ Create the test case
+
+ Parameters
+ ----------
+ test_group : compass.ocean.tests.internal_wave.InternalWave
+ The test group that this test case belongs to
+ """
+ name = 'rpe_test'
+ super().__init__(test_group=test_group, name=name)
+
+ nus = [0.01, 1, 15, 150]
+
+ self.add_step(InitialState(test_case=self))
+
+ for index, nu in enumerate(nus):
+ name = 'rpe_test_{}_nu_{:g}'.format(index + 1, nu)
+ step = Forward(
+ test_case=self, name=name, subdir=name, cores=4,
+ threads=1, nu=float(nu))
+
+ step.add_namelist_file(
+ 'compass.ocean.tests.internal_wave.rpe_test',
+ 'namelist.forward')
+ step.add_streams_file(
+ 'compass.ocean.tests.internal_wave.rpe_test',
+ 'streams.forward')
+ self.add_step(step)
+
+ self.add_step(
+ Analysis(test_case=self, nus=nus))
diff --git a/compass/ocean/tests/internal_wave/rpe_test/analysis.py b/compass/ocean/tests/internal_wave/rpe_test/analysis.py
new file mode 100644
index 0000000000..8f96c739aa
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/rpe_test/analysis.py
@@ -0,0 +1,107 @@
+import numpy as np
+from netCDF4 import Dataset
+import matplotlib.pyplot as plt
+import cmocean
+
+from compass.step import Step
+
+
+class Analysis(Step):
+ """
+ A step for plotting the results of a series of RPE runs in the internal
+ wave test group
+
+ Attributes
+ ----------
+ resolution : str
+ The resolution of the test case
+
+ nus : list of float
+ A list of viscosities
+ """
+ def __init__(self, test_case, nus):
+ """
+ Create the step
+
+ Parameters
+ ----------
+ test_case : compass.TestCase
+ The test case this step belongs to
+
+ nus : list of float
+ A list of viscosities
+ """
+ super().__init__(test_case=test_case, name='analysis')
+ self.nus = nus
+
+ for index, nu in enumerate(nus):
+ self.add_input_file(
+ filename='output_{}.nc'.format(index+1),
+ target='../rpe_test_{}_nu_{}/output.nc'.format(index+1, nu))
+
+ self.add_output_file(
+ filename='sections_internal_wave.png')
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+ section = self.config['internal_wave']
+ nx = section.getint('nx')
+ ny = section.getint('ny')
+ _plot(nx, ny, self.outputs[0], self.nus)
+
+
+def _plot(nx, ny, filename, nus):
+ """
+ Plot section of the internal wave at different viscosities
+
+ Parameters
+ ----------
+ nx : int
+ The number of cells in the x direction
+
+ ny : int
+ The number of cells in the y direction (before culling)
+
+ filename : str
+ The output file name
+
+ nus : list of float
+ The viscosity values
+ """
+
+ plt.switch_backend('Agg')
+
+ fig = plt.gcf()
+ nRow = len(nus)
+ nCol = 2
+ iTime = [0, 1]
+ time = ['1', '21']
+
+ fig, axs = plt.subplots(nRow, nCol, figsize=(
+ 4.0 * nCol, 3.7 * nRow), constrained_layout=True)
+
+ for iRow in range(nRow):
+ ncfile = Dataset('output_' + str(iRow + 1) + '.nc', 'r')
+ var = ncfile.variables['temperature']
+ xtime = ncfile.variables['xtime']
+ for iCol in range(nCol):
+ ax = axs[iRow, iCol]
+ dis = ax.imshow(
+ var[iTime[iCol], 0::4, :].T,
+ extent=[0, 250, 500, 0],
+ aspect='0.5',
+ cmap='jet',
+ vmin=10,
+ vmax=20)
+ if iRow == nRow - 1:
+ ax.set_xlabel('x, km')
+ if iCol == 0:
+ ax.set_ylabel('depth, m')
+ if iCol == nCol - 1:
+ fig.colorbar(dis, ax=axs[iRow, iCol], aspect=10)
+ ax.set_title("day {}, $\\nu_h=${}".format(time[iCol], nus[iRow]))
+ ncfile.close()
+
+ plt.savefig(filename)
diff --git a/compass/ocean/tests/internal_wave/rpe_test/namelist.forward b/compass/ocean/tests/internal_wave/rpe_test/namelist.forward
new file mode 100644
index 0000000000..e3e9d43636
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/rpe_test/namelist.forward
@@ -0,0 +1 @@
+config_run_duration = '20_00:00:00'
diff --git a/compass/ocean/tests/internal_wave/rpe_test/streams.forward b/compass/ocean/tests/internal_wave/rpe_test/streams.forward
new file mode 100644
index 0000000000..1e3455b13c
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/rpe_test/streams.forward
@@ -0,0 +1,16 @@
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/compass/ocean/tests/internal_wave/streams.forward b/compass/ocean/tests/internal_wave/streams.forward
new file mode 100644
index 0000000000..c959aa1455
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/streams.forward
@@ -0,0 +1,24 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/compass/ocean/tests/internal_wave/streams.init b/compass/ocean/tests/internal_wave/streams.init
new file mode 100644
index 0000000000..8b25ef6a4b
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/streams.init
@@ -0,0 +1,23 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/compass/ocean/tests/internal_wave/ten_day_test/__init__.py b/compass/ocean/tests/internal_wave/ten_day_test/__init__.py
new file mode 100644
index 0000000000..13848ba5ce
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/ten_day_test/__init__.py
@@ -0,0 +1,43 @@
+from compass.testcase import TestCase
+from compass.ocean.tests.internal_wave.initial_state import InitialState
+from compass.ocean.tests.internal_wave.forward import Forward
+from compass.ocean.tests.internal_wave.viz import Viz
+from compass.validate import compare_variables
+
+
+class TenDayTest(TestCase):
+ """
+ The default test case for the internal wave test
+ """
+
+ def __init__(self, test_group):
+ """
+ Create the test case
+
+ Parameters
+ ----------
+ test_group : compass.ocean.tests.internal_wave.InternalWave
+ The test group that this test case belongs to
+ """
+ name = 'ten_day_test'
+ super().__init__(test_group=test_group, name=name)
+ self.add_step(InitialState(test_case=self))
+
+ step = Forward(test_case=self, cores=4, threads=1)
+ step.add_namelist_file(
+ 'compass.ocean.tests.internal_wave.ten_day_test',
+ 'namelist.forward')
+ step.add_streams_file(
+ 'compass.ocean.tests.internal_wave.ten_day_test',
+ 'streams.forward')
+ self.add_step(step)
+
+ self.add_step(Viz(test_case=self), run_by_default=False)
+
+ def validate(self):
+ """
+ Validate variables against a baseline
+ """
+ compare_variables(test_case=self,
+ variables=['layerThickness', 'normalVelocity'],
+ filename1='forward/output.nc')
diff --git a/compass/ocean/tests/internal_wave/ten_day_test/namelist.forward b/compass/ocean/tests/internal_wave/ten_day_test/namelist.forward
new file mode 100644
index 0000000000..8e95c74fe4
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/ten_day_test/namelist.forward
@@ -0,0 +1,14 @@
+config_dt = '00:05:00'
+config_btr_dt = '00:00:15'
+config_time_integrator = 'split_explicit'
+config_run_duration = '0010_00:00:00'
+config_use_cvmix = .true.
+config_pio_stride = 4
+config_use_mom_del2 = .true.
+config_mom_del2 = 10.0
+config_implicit_bottom_drag_coeff = 1.0e-2
+config_use_cvmix_convection = .true.
+config_AM_globalStats_enable = .true.
+config_AM_globalStats_compute_on_startup = .true.
+config_AM_globalStats_write_on_startup = .true.
+config_AM_globalStats_text_file = .true.
diff --git a/compass/ocean/tests/internal_wave/ten_day_test/streams.forward b/compass/ocean/tests/internal_wave/ten_day_test/streams.forward
new file mode 100644
index 0000000000..ff5e120dc4
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/ten_day_test/streams.forward
@@ -0,0 +1,27 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/compass/ocean/tests/internal_wave/viz.py b/compass/ocean/tests/internal_wave/viz.py
new file mode 100644
index 0000000000..5e5544220f
--- /dev/null
+++ b/compass/ocean/tests/internal_wave/viz.py
@@ -0,0 +1,204 @@
+import xarray
+import numpy
+import matplotlib.pyplot as plt
+
+from compass.step import Step
+
+
+class Viz(Step):
+ """
+ A step for visualizing a cross-section through the internal wave
+ """
+ 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='viz')
+
+ self.add_input_file(filename='output.nc',
+ target='../forward/output.nc')
+ self.add_output_file('uNormal_depth_section_t0.png')
+ self.add_output_file('pt_depth_section_t0.png')
+ self.add_output_file('sa_depth_section_t0.png')
+ self.add_output_file('layerThickness_depth_section_t0.png')
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+
+ ds = xarray.open_dataset('output.nc')
+ figsize = [6.4, 4.8]
+ markersize = 20
+ if 'Time' not in ds.dims:
+ print('Dataset missing time dimension')
+ return
+ nSteps = ds.sizes['Time'] # number of timesteps
+ tend = nSteps - 1
+
+ for j in [0, tend]:
+ ds1 = ds.isel(Time=j)
+
+ # prep all variables for uNormal plot
+ ds1 = ds1.sortby('yEdge')
+
+ nCells = ds1.sizes['nCells']
+ nEdges = ds1.sizes['nEdges']
+ nVertLevels = ds1.sizes['nVertLevels']
+
+ xEdge = numpy.zeros((nEdges))
+ xEdge = ds1.xEdge
+ xCell = numpy.zeros((nCells))
+ xCell = ds1.xCell
+ yCell = numpy.zeros((nCells))
+ yCell = ds1.yCell
+
+ xEdge_mid = numpy.median(xEdge)
+ edgeMask_x = numpy.equal(xEdge, xEdge_mid)
+
+ zIndex = xarray.DataArray(data=numpy.arange(nVertLevels),
+ dims='nVertLevels')
+
+ zInterface = numpy.zeros((nCells, nVertLevels + 1))
+ zInterface[:, 0] = ds1.ssh.values
+ for zIndex in range(nVertLevels):
+ thickness = ds1.layerThickness.isel(nVertLevels=zIndex)
+ thickness = thickness.fillna(0.)
+ zInterface[:, zIndex + 1] = \
+ zInterface[:, zIndex] - thickness.values
+
+ zMid = numpy.zeros((nCells, nVertLevels))
+ for zIndex in range(nVertLevels):
+ zMid[:, zIndex] = (zInterface[:, zIndex]
+ + numpy.divide(zInterface[:, zIndex + 1]
+ - zInterface[:, zIndex], 2.))
+
+ # Solve for lateral boundaries of uNormal at cell centers for
+ # x-section
+ cellsOnEdge = ds1.cellsOnEdge
+ cellsOnEdge_x = cellsOnEdge[edgeMask_x, :]
+ yEdges = numpy.zeros((len(cellsOnEdge_x)+1))
+ for i in range(len(cellsOnEdge_x)):
+ if cellsOnEdge[i, 1] == 0:
+ yEdges[i] = yCell[cellsOnEdge_x[i, 0] - 1]
+ yEdges[i+1] = yCell[cellsOnEdge_x[i, 0] - 1]
+ elif cellsOnEdge[i, 1] == 0:
+ yEdges[i] = yCell[cellsOnEdge_x[i, 1] - 1]
+ yEdges[i+1] = yCell[cellsOnEdge_x[i, 1] - 1]
+ else:
+ yEdges[i] = min(yCell[cellsOnEdge_x[i, 0] - 1],
+ yCell[cellsOnEdge_x[i, 1] - 1])
+ yEdges[i+1] = max(yCell[cellsOnEdge_x[i, 0] - 1],
+ yCell[cellsOnEdge_x[i, 1] - 1])
+
+ zInterfaces_mesh, yEdges_mesh = numpy.meshgrid(zInterface[0, :],
+ yEdges)
+
+ normalVelocity = numpy.zeros((nCells, nVertLevels))
+ normalVelocity = ds1.normalVelocity
+ normalVelocity_xmesh = normalVelocity[edgeMask_x, :]
+
+ # Figures
+ plt.figure(figsize=figsize, dpi=100)
+ cmax = numpy.max(numpy.abs(normalVelocity_xmesh.values))
+ plt.pcolormesh(numpy.divide(yEdges_mesh, 1e3),
+ zInterfaces_mesh,
+ normalVelocity_xmesh.values,
+ cmap='RdBu', vmin=-1.*cmax, vmax=cmax)
+ plt.xlabel('y (km)')
+ plt.ylabel('z (m)')
+ cbar = plt.colorbar()
+ cbar.ax.set_title('uNormal (m/s)')
+ plt.savefig('uNormal_depth_section_t{}.png'.format(j),
+ bbox_inches='tight', dpi=200)
+ plt.close()
+
+ # ------------------------------------------------------------------
+ # Plot cell-centered variables
+ # ------------------------------------------------------------------
+ # Prep variables for cell quantities
+ cellIndex = numpy.subtract(cellsOnEdge_x[1:, 0], 1)
+ yEdge = numpy.zeros((nEdges))
+ yEdge = ds1.yEdge
+ yEdge_x = yEdge[edgeMask_x]
+
+ zInterfaces_mesh, yCells_mesh = numpy.meshgrid(zInterface[0, :],
+ yEdge_x)
+
+ # Import cell quantities
+ layerThickness = numpy.zeros((nCells, nVertLevels))
+ layerThickness = ds1.layerThickness
+ layerThickness_x = layerThickness[cellIndex, :]
+ temperature = numpy.zeros((nCells, nVertLevels))
+ temperature = ds1.temperature
+ temperature_z = temperature.mean(dim='nCells')
+ zMid_z = zMid.mean(axis=0)
+ temperature_x = temperature[cellIndex, :]
+ salinity = numpy.zeros((nCells, nVertLevels))
+ salinity = ds1.salinity
+ salinity_x = salinity[cellIndex, :]
+ w = numpy.zeros((nCells, nVertLevels))
+ w = ds1.vertVelocityTop
+ w_x = w[cellIndex, :]
+
+ # Figures
+ plt.figure(figsize=figsize, dpi=100)
+ plt.plot(temperature_z.values, zMid_z)
+ plt.xlabel('PT (C)')
+ plt.ylabel('z (m)')
+ plt.savefig('pt_depth_t{}.png'.format(j),
+ bbox_inches='tight', dpi=200)
+ plt.close()
+
+ plt.figure(figsize=figsize, dpi=100)
+ plt.pcolormesh(numpy.divide(yCells_mesh, 1e3),
+ zInterfaces_mesh,
+ temperature_x.values, cmap='viridis')
+ plt.xlabel('y (km)')
+ plt.ylabel('z (m)')
+ cbar = plt.colorbar()
+ cbar.ax.set_title('PT (C)')
+ plt.savefig('pt_depth_section_t{}.png'.format(j),
+ bbox_inches='tight', dpi=200)
+ plt.close()
+
+ plt.figure(figsize=figsize, dpi=100)
+ plt.pcolormesh(numpy.divide(yCells_mesh, 1e3),
+ zInterfaces_mesh,
+ salinity_x.values, cmap='viridis')
+ plt.xlabel('y (km)')
+ plt.ylabel('z (m)')
+ cbar = plt.colorbar()
+ cbar.ax.set_title('SA (g/kg)')
+ plt.savefig('sa_depth_section_t{}.png'.format(j),
+ bbox_inches='tight', dpi=200)
+ plt.close()
+
+ plt.figure(figsize=figsize, dpi=100)
+ plt.pcolormesh(numpy.divide(yCells_mesh, 1e3),
+ zInterfaces_mesh,
+ w_x.values, cmap='viridis')
+ plt.xlabel('y (km)')
+ plt.ylabel('z (m)')
+ cbar = plt.colorbar()
+ cbar.ax.set_title('h (m)')
+ plt.savefig('w_depth_section_t{}.png'.format(j),
+ bbox_inches='tight', dpi=200)
+ plt.close()
+
+ plt.figure(figsize=figsize, dpi=100)
+ plt.pcolormesh(numpy.divide(yCells_mesh, 1e3),
+ zInterfaces_mesh,
+ layerThickness_x.values, cmap='viridis')
+ plt.xlabel('y (km)')
+ plt.ylabel('z (m)')
+ cbar = plt.colorbar()
+ cbar.ax.set_title('h (m)')
+ plt.savefig('layerThickness_depth_section_t{}.png'.format(j),
+ bbox_inches='tight', dpi=200)
+ plt.close()
diff --git a/docs/developers_guide/ocean/test_groups/internal_wave.rst b/docs/developers_guide/ocean/test_groups/internal_wave.rst
new file mode 100644
index 0000000000..60e417b42d
--- /dev/null
+++ b/docs/developers_guide/ocean/test_groups/internal_wave.rst
@@ -0,0 +1,89 @@
+.. _dev_ocean_internal_wave:
+
+internal_wave
+=============
+
+The ``internal_wave`` test group
+(:py:class:`compass.ocean.tests.internal_wave.InternalWave`)
+implements variants of the internal wave test case. Here,
+we describe the shared framework for this test group and the 3 test cases.
+
+.. _dev_ocean_internal_wave_framework:
+
+framework
+---------
+
+The shared config options for the ``internal_wave`` test group
+are described in :ref:`ocean_internal_wave` in the User's Guide.
+
+Additionally, the test group has a shared ``namelist.forward`` file with
+a few common namelist options related to run duration and default horizontal
+and vertical momentum and tracer diffusion, as well as a shared
+``streams.forward`` file that defines ``mesh``, ``input``, ``restart``, and
+``output`` streams.
+
+initial_state
+~~~~~~~~~~~~~
+
+The class :py:class:`compass.ocean.tests.internal_wave.initial_state.InitialState`
+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 y direction. A vertical grid is generated,
+with 20 layers of 25-m thickness each by default. Finally, the initial
+temperature field is computed along with uniform salinity and zero initial
+velocity.
+
+forward
+~~~~~~~
+
+The class :py:class:`compass.ocean.tests.internal_wave.forward.Forward`
+defines a step for running MPAS-Ocean from the initial condition produced in
+the ``initial_state`` step. If ``nu`` is provided as an argument to the
+constructor, the associate namelist option (``config_mom_del2``) will be given
+this value. Namelist and streams files are generate during ``setup()`` and
+MPAS-Ocean is run (including updating PIO namelist options and generating a
+graph partition) in ``run()``.
+
+.. _dev_ocean_internal_wave_default:
+
+default
+-------
+
+The :py:class:`compass.ocean.tests.internal_wave.default.Default`
+test performs a 15-minute run on 4 cores. It doesn't contain any
+:ref:`dev_validation`.
+
+.. _dev_ocean_internal_wave_rpe_test:
+
+rpe_test
+--------
+
+The :py:class:`compass.ocean.tests.internal_wave.rpe_test.RpeTest`
+performs a longer (20 day) integration of the model forward in time at 5
+different values of the viscosity. These ``nu`` values are later added as arguments to the ``Forward`` steps'
+constructors when they are added to the test case:
+
+.. code-block:: python
+
+ step = Forward(
+ test_case=self, name=name, subdir=name, cores=cores,
+ min_cores=min_cores, nu=float(nu))
+ ...
+ self.add_step(step)
+
+The ``analysis`` step defined by
+:py:class:`compass.ocean.tests.internal_wave.rpe_test.analysis.Analysis`
+makes plots of the final results with each value of the viscosity.
+
+This test is resource intensive enough that it is not used in regression
+testing.
+
+.. _dev_ocean_internal_wave_ten_day_test:
+
+ten_day_test
+------------
+
+The :py:class:`compass.ocean.tests.internal_wave.ten_day_test.TenDayTest`
+performs a longer (10 day) integration of the model forward in time.
diff --git a/docs/users_guide/ocean/test_groups/images/internal_wave.png b/docs/users_guide/ocean/test_groups/images/internal_wave.png
new file mode 100644
index 0000000000..7464e92768
Binary files /dev/null and b/docs/users_guide/ocean/test_groups/images/internal_wave.png differ
diff --git a/docs/users_guide/ocean/test_groups/internal_wave.rst b/docs/users_guide/ocean/test_groups/internal_wave.rst
new file mode 100644
index 0000000000..79f79cd161
--- /dev/null
+++ b/docs/users_guide/ocean/test_groups/internal_wave.rst
@@ -0,0 +1,145 @@
+.. _ocean_internal_wave:
+
+internal_wave
+=============
+
+The ``ocean/internal_wave`` test group induces internal wave propagation and is documented in
+`Ilicak et al. (2012) `_.
+
+The domain is periodic on the zonal boundaries and solid on the meridional boundaries.
+Salinity is constant throughout the domain (at 35 PSU). The initial
+temperature has a linear background stratification with a sinusoidal perturbation in the center
+of the domain. This perturbation initiates symmetric waves that propage out from the center and
+then back to the center after meeting the solid boundaries.
+with a gradient between the two halves that is sinusoidally perturbed in the
+meridional direction. The surface temperature is also warmer than at depth.
+
+.. image:: images/internal_wave.png
+ :width: 500 px
+ :align: center
+
+By default, the 20 vertical layers each have 25-m uniform
+thickness. By default, all cases have a horizontal resolution of 5 km.
+
+The test group includes 3 test cases. All test cases have 2 steps,
+``initial_state``, which defines the mesh and initial conditions for the model,
+and ``forward`` (given another name in many test cases to distinguish multiple
+forward runs), which performs time integration of the model. There is an optional ``viz``
+step which performs visualization of vertical cross-sections through the center of the domain.
+
+config options
+--------------
+
+All 3 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 = uniform
+
+ # Number of vertical levels
+ vert_levels = 20
+
+ # Depth of the bottom of the ocean
+ bottom_depth = 1000.0
+
+ # The type of vertical coordinate (e.g. z-level, z-star)
+ coord_type = z-level
+
+ # 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
+
+
+ # namelist options for internal wave testcases
+ [internal_wave]
+
+ # the number of grid cells in x and y
+ nx = 4
+ ny = 50
+
+ # the size of grid cells (m)
+ dc = 5000.0
+
+ # Logical flag that determines if locations of features are defined by distance
+ # or fractions. False means fractions.
+ use_distances = False
+
+ # Temperature of the surface in the northern half of the domain.
+ surface_temperature = 20.1
+
+ # Temperature of the bottom in the northern half of the domain.
+ bottom_temperature = 10.1
+
+ # Difference in the temperature field between top and bottom
+ temperature_difference = 2.0
+
+ # Fraction of domain in Y direction the temperature gradient should be linear
+ # over.
+ amplitude_width_frac = 0.33
+
+ # Width of the temperature gradient around the center sin wave. Default value
+ # is relative to a 500km domain in Y.
+ amplitude_width_dist = 50e3
+
+ # Salinity of the water in the entire domain.
+ salinity = 35.0
+
+ # Logical flag that determines if locations of features are defined by distance
+ # or fractions. False means fractions.
+ use_distances = False
+
+ # Temperature of the surface in the northern half of the domain.
+ surface_temperature = 13.1
+
+ # Temperature of the bottom in the northern half of the domain.
+ bottom_temperature = 10.1
+
+ # Difference in the temperature field between the northern and southern halves
+ # of the domain.
+ temperature_difference = 1.2
+
+ # Fraction of domain in Y direction the temperature gradient should be linear
+ # over.
+ gradient_width_frac = 0.08
+
+ # Width of the temperature gradient around the center sin wave. Default value
+ # is relative to a 500km domain in Y.
+ gradient_width_dist = 40e3
+
+ # Salinity of the water in the entire domain.
+ salinity = 35.0
+
+ # Coriolis parameter for entire domain.
+ coriolis_parameter = -1.2e-4
+
+All units are mks, with temperature in degrees Celsius and salinity in PSU.
+
+default
+-------
+
+``ocean/internal_wave/default`` is the default version of the
+internal wave test case for a short (15 min) test run and validation of
+prognostic variables for regression testing.
+
+rpe_test
+--------
+
+Since mixing is a strong function of horizontal viscosity, this test case performs
+20-day integrations of the model forward in time at 5 different values of the viscosity (with steps
+named ``rpe_test_1_nu_1``, ``rpe_test_2_nu_5``, etc.)
+``ocean/internal_wave/rpe_test``,
+Results of these tests have been used
+to show that MPAS-Ocean has lower spurious dissipation of reference potential
+energy (RPE) than POP, MOM and MITgcm models
+(`Petersen et al. 2015 `_).
+
+ten_day_test
+------------
+
+This test is identical to ``ocean/internal_wave/default`` except that the duration is 10 days.