diff --git a/src/core_ocean/mode_init/Registry_global_ocean.xml b/src/core_ocean/mode_init/Registry_global_ocean.xml
index 03f4418034..c6422db842 100644
--- a/src/core_ocean/mode_init/Registry_global_ocean.xml
+++ b/src/core_ocean/mode_init/Registry_global_ocean.xml
@@ -11,15 +11,15 @@
description="Minimum depth where column contains all water-filled cells. The first layer with refBottomDepth greater than this value is included in whole, i.e. no partial bottom cells are used in this layer."
possible_values="Any positive real value greater than 0, but typically greater than 10 m."
/>
-
-
-
diff --git a/testing_and_setup/compass/ocean/global_ocean/ARM60to10/init/config_initial_state.xml b/testing_and_setup/compass/ocean/global_ocean/ARM60to10/init/config_initial_state.xml
index 708f1da2f9..5c7c5154d7 120000
--- a/testing_and_setup/compass/ocean/global_ocean/ARM60to10/init/config_initial_state.xml
+++ b/testing_and_setup/compass/ocean/global_ocean/ARM60to10/init/config_initial_state.xml
@@ -1 +1 @@
-../../config_files/config_initial_state_60_layers.xml
\ No newline at end of file
+../../config_files/config_initial_state.xml
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/ARM60to10/init/config_initial_state_60_layers.xml b/testing_and_setup/compass/ocean/global_ocean/ARM60to10/init/config_initial_state_60_layers.xml
new file mode 120000
index 0000000000..708f1da2f9
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/ARM60to10/init/config_initial_state_60_layers.xml
@@ -0,0 +1 @@
+../../config_files/config_initial_state_60_layers.xml
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/ARM60to10/init/define_vertical_grid.py b/testing_and_setup/compass/ocean/global_ocean/ARM60to10/init/define_vertical_grid.py
new file mode 120000
index 0000000000..4e1cac8e2b
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/ARM60to10/init/define_vertical_grid.py
@@ -0,0 +1 @@
+../../../scripts/vertical_grid/define_vertical_grid_64_layers.py
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/ARM60to6/init/config_initial_state.xml b/testing_and_setup/compass/ocean/global_ocean/ARM60to6/init/config_initial_state.xml
index 708f1da2f9..5c7c5154d7 120000
--- a/testing_and_setup/compass/ocean/global_ocean/ARM60to6/init/config_initial_state.xml
+++ b/testing_and_setup/compass/ocean/global_ocean/ARM60to6/init/config_initial_state.xml
@@ -1 +1 @@
-../../config_files/config_initial_state_60_layers.xml
\ No newline at end of file
+../../config_files/config_initial_state.xml
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/ARM60to6/init/config_initial_state_60_layers.xml b/testing_and_setup/compass/ocean/global_ocean/ARM60to6/init/config_initial_state_60_layers.xml
new file mode 120000
index 0000000000..708f1da2f9
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/ARM60to6/init/config_initial_state_60_layers.xml
@@ -0,0 +1 @@
+../../config_files/config_initial_state_60_layers.xml
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/ARM60to6/init/define_vertical_grid.py b/testing_and_setup/compass/ocean/global_ocean/ARM60to6/init/define_vertical_grid.py
new file mode 120000
index 0000000000..4e1cac8e2b
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/ARM60to6/init/define_vertical_grid.py
@@ -0,0 +1 @@
+../../../scripts/vertical_grid/define_vertical_grid_64_layers.py
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/CUSP12/init/config_initial_state_60_layers.xml b/testing_and_setup/compass/ocean/global_ocean/CUSP12/init/config_initial_state_60_layers.xml
new file mode 120000
index 0000000000..708f1da2f9
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/CUSP12/init/config_initial_state_60_layers.xml
@@ -0,0 +1 @@
+../../config_files/config_initial_state_60_layers.xml
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/CUSP12/init/define_vertical_grid.py b/testing_and_setup/compass/ocean/global_ocean/CUSP12/init/define_vertical_grid.py
new file mode 120000
index 0000000000..4e1cac8e2b
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/CUSP12/init/define_vertical_grid.py
@@ -0,0 +1 @@
+../../../scripts/vertical_grid/define_vertical_grid_64_layers.py
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/CUSP8/init/config_initial_state_60_layers.xml b/testing_and_setup/compass/ocean/global_ocean/CUSP8/init/config_initial_state_60_layers.xml
new file mode 120000
index 0000000000..708f1da2f9
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/CUSP8/init/config_initial_state_60_layers.xml
@@ -0,0 +1 @@
+../../config_files/config_initial_state_60_layers.xml
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/CUSP8/init/define_vertical_grid.py b/testing_and_setup/compass/ocean/global_ocean/CUSP8/init/define_vertical_grid.py
new file mode 120000
index 0000000000..4e1cac8e2b
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/CUSP8/init/define_vertical_grid.py
@@ -0,0 +1 @@
+../../../scripts/vertical_grid/define_vertical_grid_64_layers.py
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/EC60to30/init/config_initial_state.xml b/testing_and_setup/compass/ocean/global_ocean/EC60to30/init/config_initial_state.xml
deleted file mode 100644
index 87c01a0c1c..0000000000
--- a/testing_and_setup/compass/ocean/global_ocean/EC60to30/init/config_initial_state.xml
+++ /dev/null
@@ -1,53 +0,0 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- 128
-
-
-
-
diff --git a/testing_and_setup/compass/ocean/global_ocean/EC60to30/init/config_initial_state.xml b/testing_and_setup/compass/ocean/global_ocean/EC60to30/init/config_initial_state.xml
new file mode 120000
index 0000000000..5c7c5154d7
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/EC60to30/init/config_initial_state.xml
@@ -0,0 +1 @@
+../../config_files/config_initial_state.xml
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/EC60to30/init/config_initial_state_60_layers.xml b/testing_and_setup/compass/ocean/global_ocean/EC60to30/init/config_initial_state_60_layers.xml
new file mode 120000
index 0000000000..708f1da2f9
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/EC60to30/init/config_initial_state_60_layers.xml
@@ -0,0 +1 @@
+../../config_files/config_initial_state_60_layers.xml
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/EC60to30/init/define_vertical_grid.py b/testing_and_setup/compass/ocean/global_ocean/EC60to30/init/define_vertical_grid.py
new file mode 120000
index 0000000000..4e1cac8e2b
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/EC60to30/init/define_vertical_grid.py
@@ -0,0 +1 @@
+../../../scripts/vertical_grid/define_vertical_grid_64_layers.py
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/EC60to30wISC/init/config_initial_state.xml b/testing_and_setup/compass/ocean/global_ocean/EC60to30wISC/init/config_initial_state.xml
deleted file mode 100644
index ae9be2a62c..0000000000
--- a/testing_and_setup/compass/ocean/global_ocean/EC60to30wISC/init/config_initial_state.xml
+++ /dev/null
@@ -1,60 +0,0 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- 4
-
-
-
-
-
diff --git a/testing_and_setup/compass/ocean/global_ocean/EC60to30wISC/init/config_initial_state.xml b/testing_and_setup/compass/ocean/global_ocean/EC60to30wISC/init/config_initial_state.xml
new file mode 120000
index 0000000000..b42fcfc63a
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/EC60to30wISC/init/config_initial_state.xml
@@ -0,0 +1 @@
+../../config_files_ISC/config_initial_state.xml
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/EC60to30wISC/init/define_vertical_grid.py b/testing_and_setup/compass/ocean/global_ocean/EC60to30wISC/init/define_vertical_grid.py
new file mode 120000
index 0000000000..4e1cac8e2b
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/EC60to30wISC/init/define_vertical_grid.py
@@ -0,0 +1 @@
+../../../scripts/vertical_grid/define_vertical_grid_64_layers.py
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/QU240/init/config_initial_state_60_layers.xml b/testing_and_setup/compass/ocean/global_ocean/QU240/init/config_initial_state_60_layers.xml
new file mode 120000
index 0000000000..708f1da2f9
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/QU240/init/config_initial_state_60_layers.xml
@@ -0,0 +1 @@
+../../config_files/config_initial_state_60_layers.xml
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/QU240/init/define_vertical_grid.py b/testing_and_setup/compass/ocean/global_ocean/QU240/init/define_vertical_grid.py
new file mode 120000
index 0000000000..e92ea74842
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/QU240/init/define_vertical_grid.py
@@ -0,0 +1 @@
+../../../scripts/vertical_grid/define_vertical_grid_16_layers.py
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/QU240wISC/init/define_vertical_grid.py b/testing_and_setup/compass/ocean/global_ocean/QU240wISC/init/define_vertical_grid.py
new file mode 120000
index 0000000000..4e1cac8e2b
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/QU240wISC/init/define_vertical_grid.py
@@ -0,0 +1 @@
+../../../scripts/vertical_grid/define_vertical_grid_64_layers.py
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/SO60to10wISC/init/define_vertical_grid.py b/testing_and_setup/compass/ocean/global_ocean/SO60to10wISC/init/define_vertical_grid.py
new file mode 120000
index 0000000000..4e1cac8e2b
--- /dev/null
+++ b/testing_and_setup/compass/ocean/global_ocean/SO60to10wISC/init/define_vertical_grid.py
@@ -0,0 +1 @@
+../../../scripts/vertical_grid/define_vertical_grid_64_layers.py
\ No newline at end of file
diff --git a/testing_and_setup/compass/ocean/global_ocean/config_files/config_base_mesh.xml b/testing_and_setup/compass/ocean/global_ocean/config_files/config_base_mesh.xml
index a4685355bf..2f67262259 100644
--- a/testing_and_setup/compass/ocean/global_ocean/config_files/config_base_mesh.xml
+++ b/testing_and_setup/compass/ocean/global_ocean/config_files/config_base_mesh.xml
@@ -2,7 +2,7 @@
-
+
diff --git a/testing_and_setup/compass/ocean/global_ocean/config_files/config_initial_state.xml b/testing_and_setup/compass/ocean/global_ocean/config_files/config_initial_state.xml
index aa2bb895e0..76bcfaacb5 100644
--- a/testing_and_setup/compass/ocean/global_ocean/config_files/config_initial_state.xml
+++ b/testing_and_setup/compass/ocean/global_ocean/config_files/config_initial_state.xml
@@ -28,7 +28,9 @@
-
+
+
+
@@ -38,7 +40,9 @@
-
+
+
+
@@ -50,8 +54,13 @@
- 8
+ 4
-
-
+
+
+
+
+
+
+
diff --git a/testing_and_setup/compass/ocean/global_ocean/config_files/config_initial_state_60_layers.xml b/testing_and_setup/compass/ocean/global_ocean/config_files/config_initial_state_60_layers.xml
index ce7a561a89..d0706bfdce 100644
--- a/testing_and_setup/compass/ocean/global_ocean/config_files/config_initial_state_60_layers.xml
+++ b/testing_and_setup/compass/ocean/global_ocean/config_files/config_initial_state_60_layers.xml
@@ -1,5 +1,5 @@
-
+
diff --git a/testing_and_setup/compass/ocean/global_ocean/config_files_ISC/config_initial_state.xml b/testing_and_setup/compass/ocean/global_ocean/config_files_ISC/config_initial_state.xml
index 862247bf67..df9c15b3cd 100644
--- a/testing_and_setup/compass/ocean/global_ocean/config_files_ISC/config_initial_state.xml
+++ b/testing_and_setup/compass/ocean/global_ocean/config_files_ISC/config_initial_state.xml
@@ -1,9 +1,5 @@
-
-
-
-
@@ -28,8 +24,10 @@
+
-
+
+
@@ -41,7 +39,9 @@
-
+
+
+
@@ -56,7 +56,11 @@
4
+
+
-
+
+
+
diff --git a/testing_and_setup/compass/ocean/global_ocean/config_files_ISC/config_ssh_adjustment.xml b/testing_and_setup/compass/ocean/global_ocean/config_files_ISC/config_ssh_adjustment.xml
index 8bdf73a8d6..48af3fd9ac 100644
--- a/testing_and_setup/compass/ocean/global_ocean/config_files_ISC/config_ssh_adjustment.xml
+++ b/testing_and_setup/compass/ocean/global_ocean/config_files_ISC/config_ssh_adjustment.xml
@@ -4,7 +4,7 @@
-
+
diff --git a/testing_and_setup/compass/ocean/isomip_plus/2km/Ocean0/config_adjust_ssh.xml b/testing_and_setup/compass/ocean/isomip_plus/2km/Ocean0/config_adjust_ssh.xml
index 08ffa1304a..a4d4ec93d2 100644
--- a/testing_and_setup/compass/ocean/isomip_plus/2km/Ocean0/config_adjust_ssh.xml
+++ b/testing_and_setup/compass/ocean/isomip_plus/2km/Ocean0/config_adjust_ssh.xml
@@ -7,7 +7,7 @@
-
+
diff --git a/testing_and_setup/compass/ocean/isomip_plus/2km/Ocean1/config_adjust_ssh.xml b/testing_and_setup/compass/ocean/isomip_plus/2km/Ocean1/config_adjust_ssh.xml
index 08ffa1304a..a4d4ec93d2 100644
--- a/testing_and_setup/compass/ocean/isomip_plus/2km/Ocean1/config_adjust_ssh.xml
+++ b/testing_and_setup/compass/ocean/isomip_plus/2km/Ocean1/config_adjust_ssh.xml
@@ -7,7 +7,7 @@
-
+
diff --git a/testing_and_setup/compass/ocean/isomip_plus/2km/Ocean2/config_adjust_ssh.xml b/testing_and_setup/compass/ocean/isomip_plus/2km/Ocean2/config_adjust_ssh.xml
index 08ffa1304a..a4d4ec93d2 100644
--- a/testing_and_setup/compass/ocean/isomip_plus/2km/Ocean2/config_adjust_ssh.xml
+++ b/testing_and_setup/compass/ocean/isomip_plus/2km/Ocean2/config_adjust_ssh.xml
@@ -7,7 +7,7 @@
-
+
diff --git a/testing_and_setup/compass/ocean/isomip_plus/2km/time_varying_Ocean0/config_adjust_ssh.xml b/testing_and_setup/compass/ocean/isomip_plus/2km/time_varying_Ocean0/config_adjust_ssh.xml
index bc11fa2202..ea1cf884a9 100644
--- a/testing_and_setup/compass/ocean/isomip_plus/2km/time_varying_Ocean0/config_adjust_ssh.xml
+++ b/testing_and_setup/compass/ocean/isomip_plus/2km/time_varying_Ocean0/config_adjust_ssh.xml
@@ -7,7 +7,7 @@
-
+
diff --git a/testing_and_setup/compass/ocean/isomip_plus/5km/Ocean0/config_adjust_ssh.xml b/testing_and_setup/compass/ocean/isomip_plus/5km/Ocean0/config_adjust_ssh.xml
index fae3beafbf..8fca776ce0 100644
--- a/testing_and_setup/compass/ocean/isomip_plus/5km/Ocean0/config_adjust_ssh.xml
+++ b/testing_and_setup/compass/ocean/isomip_plus/5km/Ocean0/config_adjust_ssh.xml
@@ -7,7 +7,7 @@
-
+
diff --git a/testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py b/testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py
index 3ee781ad8b..26d3410361 100755
--- a/testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py
+++ b/testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py
@@ -28,8 +28,7 @@
from jigsaw_to_MPAS.inject_preserve_floodplain import \
inject_preserve_floodplain
-from define_base_mesh import define_base_mesh
-
+import define_base_mesh
def build_mesh(
preserve_floodplain=False,
diff --git a/testing_and_setup/compass/ocean/scripts/plot_globalStats.py b/testing_and_setup/compass/ocean/scripts/plots/plot_globalStats.py
similarity index 100%
rename from testing_and_setup/compass/ocean/scripts/plot_globalStats.py
rename to testing_and_setup/compass/ocean/scripts/plots/plot_globalStats.py
diff --git a/testing_and_setup/compass/ocean/scripts/plots/plot_initial_state.py b/testing_and_setup/compass/ocean/scripts/plots/plot_initial_state.py
new file mode 100755
index 0000000000..41756a2d4c
--- /dev/null
+++ b/testing_and_setup/compass/ocean/scripts/plots/plot_initial_state.py
@@ -0,0 +1,80 @@
+#!/usr/bin/env python
+"""
+This script creates plots of the initial condition.
+"""
+# import modules
+# {{{
+import matplotlib.pyplot as plt
+from netCDF4 import Dataset
+import numpy as np
+import argparse
+import matplotlib
+matplotlib.use('Agg')
+# }}}
+
+
+def main():
+ # parser
+ # {{{
+ parser = argparse.ArgumentParser(
+ description=__doc__,
+ formatter_class=argparse.RawTextHelpFormatter)
+ parser.add_argument(
+ '-i', '--input_file_name', dest='input_file_name',
+ default='initial_state.nc',
+ help='MPAS file name for input of initial state.')
+ parser.add_argument(
+ '-o', '--output_file_name', dest='output_file_name',
+ default='initial_state.png',
+ help='File name for output image.')
+ args = parser.parse_args()
+ # }}}
+
+ # load mesh variables
+ ncfile = Dataset(args.input_file_name, 'r')
+ nCells = ncfile.dimensions['nCells'].size
+ nEdges = ncfile.dimensions['nEdges'].size
+ nVertLevels = ncfile.dimensions['nVertLevels'].size
+
+ fig = plt.figure()
+ fig.set_size_inches(16.0, 8.0)
+ plt.clf()
+
+ print('plotting histograms of the initial condition')
+ txt = \
+ 'number cells: {}\n'.format(nCells) + \
+ 'number cells, millions: {:6.3f}\n'.format(nCells / 1.e6) + \
+ 'number layers: {}\n'.format(nVertLevels)
+ print(txt)
+ plt.subplot(2, 2, 1)
+ plt.text(0, 0, txt, fontsize=12)
+ plt.axis('off')
+
+ maxLevelCell = ncfile.variables['maxLevelCell']
+ plt.subplot(2, 2, 2)
+ plt.hist(maxLevelCell, bins=nVertLevels - 4)
+ plt.xlabel('maxLevelCell')
+ plt.ylabel('frequency')
+
+ bottomDepth = ncfile.variables['bottomDepth']
+ plt.subplot(2, 2, 4)
+ plt.hist(bottomDepth, bins=nVertLevels - 4)
+ plt.xlabel('bottomDepth [m]')
+ plt.ylabel('frequency')
+
+ rx1Edge = ncfile.variables['rx1Edge']
+ plt.subplot(2, 2, 3)
+ plt.hist(np.squeeze(rx1Edge[0, :]))
+ plt.xlabel('Haney Number')
+ plt.ylabel('frequency')
+ plt.title('Haney Number, max={:4.2f}'.format(
+ np.max(np.squeeze(rx1Edge[0, :, :]))))
+
+ plt.savefig(args.output_file_name)
+
+
+if __name__ == '__main__':
+ # If called as a primary module, run main
+ main()
+
+# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python
diff --git a/testing_and_setup/compass/ocean/scripts/vertical_grid/define_vertical_grid_16_layers.py b/testing_and_setup/compass/ocean/scripts/vertical_grid/define_vertical_grid_16_layers.py
new file mode 100755
index 0000000000..335c8a473b
--- /dev/null
+++ b/testing_and_setup/compass/ocean/scripts/vertical_grid/define_vertical_grid_16_layers.py
@@ -0,0 +1,12 @@
+#!/usr/bin/env python
+
+from make_vertical_grid import create_vertical_grid
+
+# Standard 16 layer vertical grid
+create_vertical_grid(
+ num_vert_levels=16,
+ bottom_depth=3000,
+ min_layer_thickness=3.,
+ max_layer_thickness=500.,
+ plot_vertical_grid=True,
+ outfile='vertical_grid.nc')
diff --git a/testing_and_setup/compass/ocean/scripts/vertical_grid/define_vertical_grid_64_layers.py b/testing_and_setup/compass/ocean/scripts/vertical_grid/define_vertical_grid_64_layers.py
new file mode 100755
index 0000000000..36f9997a91
--- /dev/null
+++ b/testing_and_setup/compass/ocean/scripts/vertical_grid/define_vertical_grid_64_layers.py
@@ -0,0 +1,12 @@
+#!/usr/bin/env python
+
+from make_vertical_grid import create_vertical_grid
+
+# Standard 64 layer vertical grid
+create_vertical_grid(
+ num_vert_levels=64,
+ bottom_depth=6000,
+ min_layer_thickness=2.,
+ max_layer_thickness=210.,
+ plot_vertical_grid=True,
+ outfile='vertical_grid.nc')
diff --git a/testing_and_setup/compass/ocean/scripts/vertical_grid/define_vertical_grid_80_layers.py b/testing_and_setup/compass/ocean/scripts/vertical_grid/define_vertical_grid_80_layers.py
new file mode 100755
index 0000000000..ebf3bc8ef7
--- /dev/null
+++ b/testing_and_setup/compass/ocean/scripts/vertical_grid/define_vertical_grid_80_layers.py
@@ -0,0 +1,12 @@
+#!/usr/bin/env python
+
+from make_vertical_grid import create_vertical_grid
+
+# Standard 80 layer vertical grid
+create_vertical_grid(
+ num_vert_levels=80,
+ bottom_depth=6000,
+ min_layer_thickness=1.,
+ max_layer_thickness=200.,
+ plot_vertical_grid=True,
+ outfile='vertical_grid.nc')
diff --git a/testing_and_setup/compass/ocean/scripts/vertical_grid/make_vertical_grid.py b/testing_and_setup/compass/ocean/scripts/vertical_grid/make_vertical_grid.py
new file mode 100755
index 0000000000..e9a59d0295
--- /dev/null
+++ b/testing_and_setup/compass/ocean/scripts/vertical_grid/make_vertical_grid.py
@@ -0,0 +1,278 @@
+#!/usr/bin/env python
+"""
+This script creates the vertical grid for MPAS-Ocean and writes it to a netcdf
+file.
+"""
+# import modules
+# {{{
+import matplotlib.pyplot as plt
+from netCDF4 import Dataset
+import numpy as np
+import argparse
+from scipy.optimize import root_scalar
+import matplotlib
+matplotlib.use('Agg')
+# }}}
+
+
+def main():
+ # parser
+ # {{{
+ parser = argparse.ArgumentParser(
+ description=__doc__,
+ formatter_class=argparse.RawTextHelpFormatter)
+ parser.add_argument(
+ '-nz', '--num_vert_levels', dest='nz',
+ default=64,
+ help='Number of vertical levels for the grid',
+ type=int)
+ parser.add_argument(
+ '-bd', '--bottom_depth', dest='bottom_depth',
+ default=5000.0,
+ help='bottom depth for the chosen vertical coordinate [m]',
+ type=float)
+ parser.add_argument(
+ '-dz1', '--min_layer_thickness', dest='dz1',
+ default=2.0,
+ help='Target thickness of the first layer [m]',
+ type=float)
+ parser.add_argument(
+ '-dz2', '--maxLayer_thickness', dest='dz2',
+ default=250.0,
+ help='Target maximum thickness in column [m]',
+ type=float)
+ parser.add_argument(
+ '-p', '--plot_vertical_grid', dest='plot_vertical_grid',
+ action='store_true')
+ parser.add_argument(
+ '-o', '--output_file_name', dest='output_file_name',
+ default='MPAS-Ocean_vertical_grid.nc',
+ help='MPAS file name for output of vertical grid.',
+ metavar='NAME')
+ args = parser.parse_args()
+
+ create_vertical_grid(args.bottom_depth, args.nz, args.dz1,
+ args.dz2, args.plot_vertical_grid,
+ outfile=args.output_file_name)
+# }}}
+
+
+def create_vertical_grid(
+ num_vert_levels=64,
+ bottom_depth=5000.0,
+ min_layer_thickness=2.0,
+ max_layer_thickness=250.0,
+ plot_vertical_grid=False,
+ outfile='MPAS-Ocean_vertical_grid.nc'):
+ # {{{
+ """
+ This function creates the vertical grid for MPAS-Ocean and writes it to a
+ NetCDF file.
+
+ Parameters
+ ----------
+
+ bottom_depth : float, optional
+ bottom depth for the chosen vertical coordinate [m]
+
+ num_vert_levels : int, optional
+ Number of vertical levels for the grid
+
+ min_layer_thickness : float, optional
+ Target thickness of the first layer [m]
+
+ max_layer_thickness : float, optional
+ Target maximum thickness in column [m]
+
+ plot_vertical_grid : bool, optional
+ Whether to plot the vertical grid
+
+ outfile : str, optional
+ MPAS file name for output of vertical grid
+ """
+
+ print('Creating mesh with ', num_vert_levels, ' layers...')
+ nz = num_vert_levels
+ dz1 = min_layer_thickness
+ dz2 = max_layer_thickness
+ # open a new netCDF file for writing.
+ ncfile = Dataset(outfile, 'w')
+ # create the depth_t dimension.
+ ncfile.createDimension('nVertLevels', nz)
+
+ refBottomDepth = ncfile.createVariable(
+ 'refBottomDepth', np.dtype('float64').char, ('nVertLevels',))
+ refMidDepth = ncfile.createVariable(
+ 'refMidDepth', np.dtype('float64').char, ('nVertLevels',))
+ refLayerThickness = ncfile.createVariable(
+ 'refLayerThickness', np.dtype('float64').char, ('nVertLevels',))
+
+ # the bracket here is large enough that it should hopefully encompass any
+ # reasonable value of delta, the characteristic length scale over which
+ # dz varies. The args are passed on to the match_bottom function below,
+ # and the root finder will determine a value of delta (sol.root) such that
+ # match_bottom is within a tolerance of zero, meaning the bottom of the
+ # coordinate computed by cumsum_z hits bottom_depth almost exactly
+ sol = root_scalar(match_bottom, method='brentq',
+ bracket=[dz1, 10 * bottom_depth],
+ args=(nz, dz1, dz2, bottom_depth))
+
+ delta = sol.root
+ layerThickness, z = cumsum_z(delta, nz, dz1, dz2)
+ nVertLevels = nz
+ botDepth = -z[1:]
+ midDepth = -0.5 * (z[0:-1] + z[1:])
+
+ refBottomDepth[:] = botDepth
+ refMidDepth[:] = midDepth
+ refLayerThickness[:] = layerThickness[:nVertLevels]
+ ncfile.close()
+
+ if plot_vertical_grid:
+ fig = plt.figure()
+ fig.set_size_inches(16.0, 8.0)
+ zInd = np.arange(1, nVertLevels + 1)
+ plt.clf()
+
+ plt.subplot(2, 2, 1)
+ plt.plot(zInd, midDepth, '.')
+ plt.gca().invert_yaxis()
+ plt.xlabel('vertical index (one-based)')
+ plt.ylabel('layer mid-depth [m]')
+ plt.grid()
+
+ plt.subplot(2, 2, 2)
+ plt.plot(layerThickness, midDepth, '.')
+ plt.gca().invert_yaxis()
+ plt.xlabel('layer thickness [m]')
+ plt.ylabel('layer mid-depth [m]')
+ plt.grid()
+
+ plt.subplot(2, 2, 3)
+ plt.plot(zInd, layerThickness, '.')
+ plt.xlabel('vertical index (one-based)')
+ plt.ylabel('layer thickness [m]')
+ plt.grid()
+
+ txt = \
+ 'number layers: {}\n'.format(nz) + \
+ 'bottom depth requested: {:8.2f}\n'.format(bottom_depth) + \
+ 'bottom depth actual: {:8.2f}\n'.format(np.amax(botDepth)) + \
+ 'min thickness reqeusted: {:8.2f}\n'.format(min_layer_thickness) + \
+ 'min thickness actual: {:8.2f}\n'.format(np.amin(layerThickness[:])) + \
+ 'max thickness reqeusted: {:8.2f}\n'.format(max_layer_thickness) + \
+ 'max thickness actual: {:8.2f}'.format(np.amax(layerThickness[:]))
+ print(txt)
+ plt.subplot(2, 2, 4)
+ plt.text(0, 0, txt, fontsize=12)
+ plt.axis('off')
+ plt.savefig('vertical_grid.png')
+
+# }}}
+
+
+def match_bottom(delta, nz, dz1, dz2, bottom_depth):
+ """
+ Compute the difference between the bottom depth computed with the given
+ parameters and the target ``bottom_depth``, used in the root finding
+ algorithm to determine which value of ``delta`` to use.
+
+ Parameters
+ ----------
+ delta : float
+ The characteristic length scale over which dz varies (this parameter
+ will be optimized to hit a target depth in a target number of layers)
+
+ nz : int
+ The number of layers
+
+ dz1 : float
+ The layer thickness at the top of the ocean (z = 0)
+
+ dz2 : float
+ The layer thickness at z --> -infinity
+
+ bottom_depth: float
+ depth of the bottom of the ocean that should match the bottom layer
+ interface. Note: the bottom_depth is positive, whereas the layer
+ interfaces are negative.
+
+ Returns
+ -------
+ diff : float
+ The computed bottom depth minus the target ``bottom_depth``. ``diff``
+ should be zero when we have found the desired ``delta``.
+ """
+ _, z = cumsum_z(delta, nz, dz1, dz2)
+ diff = -bottom_depth - z[-1]
+ return diff
+
+
+def cumsum_z(delta, nz, dz1, dz2):
+ """
+ Compute layer interface depths and layer thicknesses over ``nz`` layers
+
+ Parameters
+ ----------
+ delta : float
+ The characteristic length scale over which dz varies (this parameter
+ will be optimized to hit a target depth in a target number of layers)
+
+ nz : int
+ The number of layers
+
+ dz1 : float
+ The layer thickness at the top of the ocean (z = 0)
+
+ dz2 : float
+ The layer thickness at z --> -infinity
+
+ Returns
+ -------
+ dz : numpy.ndarray
+ The layer thicknesses for each layer
+
+ z : numpy.ndarray
+ The depth (positive up) of each layer interface (``nz + 1`` total
+ elements)
+ """
+ dz = np.zeros(nz)
+ z = np.zeros(nz + 1)
+ for zindex in range(nz):
+ dz[zindex] = dz_z(z[zindex], dz1, dz2, delta)
+ z[zindex + 1] = z[zindex] - dz[zindex]
+ return dz, z
+
+
+def dz_z(z, dz1, dz2, delta):
+ """
+ layer thickness as a funciton of depth
+
+ Parameters
+ ----------
+ z : float
+ Depth coordinate (positive up) at which to find the layer thickness
+
+ dz1 : float
+ The layer thickness at the top of the ocean (z = 0)
+
+ dz2 : float
+ The layer thickness at z --> -infinity
+
+ delta : float
+ The characteristic length scale over which dz varies (this parameter
+ will be optimized to hit a target depth in a target numer of layers)
+
+ Returns
+ -------
+ dz : float
+ The layer thickness
+ """
+ return (dz2 - dz1) * np.tanh(-z * np.pi / delta) + dz1
+
+
+if __name__ == '__main__':
+ # If called as a primary module, run main
+ main()
+
+# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python