From 0d5e14f0c548120f64a90ba3259b3e823af39ea2 Mon Sep 17 00:00:00 2001
From: Mark Petersen <mpetersen@lanl.gov>
Date: Mon, 30 Mar 2020 06:12:00 -0600
Subject: [PATCH 1/9] original script from Luke

---
 .../ocean/scripts/make_vertical_grid.py       | 103 ++++++++++++++++++
 1 file changed, 103 insertions(+)
 create mode 100755 testing_and_setup/compass/ocean/scripts/make_vertical_grid.py

diff --git a/testing_and_setup/compass/ocean/scripts/make_vertical_grid.py b/testing_and_setup/compass/ocean/scripts/make_vertical_grid.py
new file mode 100755
index 0000000000..61459ece37
--- /dev/null
+++ b/testing_and_setup/compass/ocean/scripts/make_vertical_grid.py
@@ -0,0 +1,103 @@
+def dz_z(z,Hmax,epsilon,dzmax):
+    import numpy as np
+    return dzmax*np.tanh(-z*np.pi/Hmax) + epsilon
+
+from netCDF4 import Dataset
+import numpy as np
+"""
+Write vertical grid to a netcdf file
+"""
+def create_vertical_grid(bottom_depth, nz, dz1, dz2, maxit=1000,
+    epsilon=1e-2, outFile='MPAS-Ocean_vertical_grid.nc'):
+
+  # open a new netCDF file for writing.
+  ncfile = Dataset(outFile,'w')
+  # create the depth_t dimension.
+  ncfile.createDimension('dim_layer',nz)
+
+  refBottomDepth = ncfile.createVariable('refBottomDepth',np.dtype('float64').char,('dim_layer'))
+  refMidDepth = ncfile.createVariable('refMidDepth',np.dtype('float64').char,('dim_layer'))
+  refLayerThickness = ncfile.createVariable('refLayerThickness',np.dtype('float64').char,('dim_layer'))
+
+  Hmax = bottom_depth
+  nLayers = 0
+
+  layerThickness = np.zeros(nz)
+  dz = [epsilon]
+  z = [0]
+  count=0
+
+  while nLayers != nz and count < maxit:
+    zval = -epsilon
+    dz = [epsilon]
+    z = [0]
+    nLayers = 0
+    while zval > -Hmax:
+        difference = dz_z(zval,Hmax,epsilon,dz2) - zval
+        while abs(difference) > 0.3:
+            zval -= epsilon
+            difference = dz_z(zval,Hmax,epsilon,dz2) -z[nLayers] + zval
+        z.append(zval)
+        dz.append(dz_z(zval,Hmax,epsilon,dz2))
+        nLayers += 1
+        zval -= epsilon
+
+    dz_arr = np.asarray(dz)
+    ind = abs(dz_arr - dz1).argmin()
+
+    dztemp = dz_arr[ind:]
+    nLayers = len(dztemp)
+    change = nz-nLayers
+    dz2 -= float(change)
+    count += 1
+
+  if count >= maxit:
+    print("Error: grid did not converge, adjust parameters")
+  else:
+    layerThickness[:nLayers] = dztemp
+
+    nVertLevels = nz 
+    botDepth = np.zeros(nVertLevels)
+    midDepth = np.zeros(nVertLevels)
+    botDepth[0] = layerThickness[0]
+    midDepth[0] = 0.5*layerThickness[0]
+
+    for i in range(1,nVertLevels):
+      botDepth[i] = botDepth[i-1] + layerThickness[i]
+      midDepth[i] = midDepth[i-1] + 0.5*(layerThickness[i] + layerThickness[i-1])
+
+    refBottomDepth[:] = botDepth
+    refMidDepth[:] = midDepth
+    refLayerThickness[:] = layerThickness[:nVertLevels]
+
+    ncfile.close()
+
+if __name__ == "__main__":
+  import argparse
+  parser = argparse.ArgumentParser(description=__doc__,
+                                  formatter_class=argparse.RawTextHelpFormatter)
+  parser.add_argument("-o", "--output_file_pattern", dest="output_filename_pattern",
+      default='MPAS-Ocean_vertical_grid.nc',
+      help="MPAS Filename pattern for file output of vertical grid.", metavar="NAME")
+  parser.add_argument("-bd", "--bottom_depth", dest="bottom_depth",
+      help="bottom depth for the chosen vertical coordinate", 
+      type=float,required=True)
+  parser.add_argument("-nz", "--num_vert_levels", dest="nz",
+      help="Number of vertical levels for the grid", type=int,
+      required=True)
+  parser.add_argument("-dz1", "--layer1_thickness", dest="dz1",
+      help="Target thickness of the first layer", 
+      type=float, required=True)
+  parser.add_argument("-dz2", "--maxLayer_thickness", dest="dz2",
+      help="Target maximum thickness in column", 
+      type=float, required=True)
+  parser.add_argument("-eps", "--error_tolerance", dest="epsilon",
+      help="Threshold for iterations", type=float, default=1e-2)
+  parser.add_argument("-maxit", "--max_iterations", dest="maxit",
+      help="maximum number of iterations for grid convergences",
+      type=int, default=1000)
+  args = parser.parse_args()
+
+  create_vertical_grid(args.bottom_depth,args.nz,args.dz1,
+      args.dz2,maxit=args.maxit,epsilon=args.epsilon,
+      outFile=args.output_filename_pattern)

From c11d2c65e35c16e09cf0ad594353341035533b09 Mon Sep 17 00:00:00 2001
From: Mark Petersen <mpetersen@lanl.gov>
Date: Mon, 30 Mar 2020 07:50:26 -0600
Subject: [PATCH 2/9] Add 64 layer to domains, add plotting, set parameters

---
 .../mode_init/Registry_global_ocean.xml       |   6 +-
 .../ARM60to10/init/config_initial_state.xml   |   2 +-
 .../ARM60to6/init/config_initial_state.xml    |   2 +-
 .../EC60to30/init/config_initial_state.xml    |  54 +----
 .../EC60to30/init/define_vertical_grid.py     |   1 +
 .../EC60to30wISC/init/define_vertical_grid.py |   1 +
 .../QU240/init/define_vertical_grid.py        |   1 +
 .../QU240wISC/init/define_vertical_grid.py    |   1 +
 .../config_files/config_base_mesh.xml         |   2 +-
 .../config_files/config_initial_state.xml     |  12 +-
 .../config_files_ISC/config_initial_state.xml |  19 +-
 .../ocean/jigsaw_to_MPAS/build_mesh.py        |   3 +-
 .../ocean/scripts/make_vertical_grid.py       | 103 ---------
 .../define_vertical_grid_64_layers.py         |  14 ++
 .../vertical_grid/make_vertical_grid.py       | 198 ++++++++++++++++++
 15 files changed, 245 insertions(+), 174 deletions(-)
 mode change 100644 => 120000 testing_and_setup/compass/ocean/global_ocean/EC60to30/init/config_initial_state.xml
 create mode 120000 testing_and_setup/compass/ocean/global_ocean/EC60to30/init/define_vertical_grid.py
 create mode 120000 testing_and_setup/compass/ocean/global_ocean/EC60to30wISC/init/define_vertical_grid.py
 create mode 120000 testing_and_setup/compass/ocean/global_ocean/QU240/init/define_vertical_grid.py
 create mode 120000 testing_and_setup/compass/ocean/global_ocean/QU240wISC/init/define_vertical_grid.py
 delete mode 100755 testing_and_setup/compass/ocean/scripts/make_vertical_grid.py
 create mode 100755 testing_and_setup/compass/ocean/scripts/vertical_grid/define_vertical_grid_64_layers.py
 create mode 100755 testing_and_setup/compass/ocean/scripts/vertical_grid/make_vertical_grid.py

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."
 		/>
-		<nml_option name="config_global_ocean_depth_file" type="character" default_value="none" units="unitless"
+		<nml_option name="config_global_ocean_depth_file" type="character" default_value="vertical_grid.nc" units="unitless"
 			description="Path to the depth initial condition file."
 			possible_values="path/to/temperature/file.nc"
 		/>
-		<nml_option name="config_global_ocean_depth_dimname" type="character" default_value="none" units="unitless"
+		<nml_option name="config_global_ocean_depth_dimname" type="character" default_value="nVertLevels" units="unitless"
 			description="Name of the dimension defining the number of vertical levels in input files."
 			possible_values="Dim name from input files."
 		/>
-		<nml_option name="config_global_ocean_depth_varname" type="character" default_value="none" units="unitless"
+		<nml_option name="config_global_ocean_depth_varname" type="character" default_value="refMidDepth" units="unitless"
 			description="Name of the variable containing mid-depth of levels in temperature and salinity initial condition files."
 			possible_values="Variable name from input files."
 		/>
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/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/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 @@
-<?xml version="1.0"?>
-<config case="initial_state">
-
-	<get_file dest_path="initial_condition_database" file_name="PotentialTemperature.01.filled.60levels.PHC.151106.nc">
-		<mirror protocol="wget" url="https://web.lcrc.anl.gov/public/e3sm/mpas_standalonedata/mpas-ocean/initial_condition_database/"/>
-	</get_file>
-
-	<get_file dest_path="initial_condition_database" file_name="Salinity.01.filled.60levels.PHC.151106.nc">
-		<mirror protocol="wget" url="https://web.lcrc.anl.gov/public/e3sm/mpas_standalonedata/mpas-ocean/initial_condition_database/"/>
-	</get_file>
-
-	<get_file dest_path="initial_condition_database" file_name="windStress.ncep_1958-2000avg.interp3600x2431.151106.nc">
-		<mirror protocol="wget" url="https://web.lcrc.anl.gov/public/e3sm/mpas_standalonedata/mpas-ocean/initial_condition_database/"/>
-	</get_file>
-
-	<get_file dest_path="initial_condition_database" file_name="ETOPO2v2c_f4_151106.nc">
-		<mirror protocol="wget" url="https://web.lcrc.anl.gov/public/e3sm/mpas_standalonedata/mpas-ocean/initial_condition_database/"/>
-	</get_file>
-
-	<get_file dest_path="initial_condition_database" file_name="chlorophyllA_monthly_averages_1deg.151201.nc">
-		<mirror protocol="wget" url="https://web.lcrc.anl.gov/public/e3sm/mpas_standalonedata/mpas-ocean/initial_condition_database/"/>
-	</get_file>
-
-
-	<add_link source="../culled_mesh/culled_mesh.nc" dest="mesh.nc"/>
-	<add_link source="../culled_mesh/culled_graph.info" dest="graph.info"/>
-	<add_link source="../culled_mesh/critical_passages_mask_final.nc" dest="critical_passages.nc"/>
-	<add_link source_path="initial_condition_database" source="PotentialTemperature.01.filled.60levels.PHC.151106.nc" dest="temperature.nc"/>
-	<add_link source_path="initial_condition_database" source="Salinity.01.filled.60levels.PHC.151106.nc" dest="salinity.nc"/>
-    <add_link source_path="initial_condition_database" source="Salinity.01.filled.60levels.PHC.151106.nc" dest="layer_depth.nc"/>
-	<add_link source_path="initial_condition_database" source="windStress.ncep_1958-2000avg.interp3600x2431.151106.nc" dest="wind_stress.nc"/>
-	<add_link source_path="initial_condition_database" source="ETOPO2v2c_f4_151106.nc" dest="topography.nc"/>
-	<add_link source_path="initial_condition_database" source="chlorophyllA_monthly_averages_1deg.151201.nc" dest="swData.nc"/>
-
-	<namelist name="namelist.ocean" mode="init">
-		<template file="template_initial_state.xml" path_base="script_core_dir" path="global_ocean"/>
-		<template file="template_critical_passages.xml" path_base="script_core_dir" path="global_ocean"/>
-        <option name="config_global_ocean_depth_conversion_factor">0.01</option>
-        <option name="config_global_ocean_tracer_depth_conversion_factor">0.01</option>
-	</namelist>
-
-	<streams name="streams.ocean" keep="immutable" mode="init">
-		<template file="template_initial_state.xml" path_base="script_core_dir" path="global_ocean"/>
-		<template file="template_critical_passages.xml" path_base="script_core_dir" path="global_ocean"/>
-	</streams>
-
-	<run_script name="run.py">
-		<step executable="gpmetis">
-			<argument flag="graph.info">128</argument>
-		</step>
-		<model_run procs="128" threads="1" namelist="namelist.ocean" streams="streams.ocean"/>
-	</run_script>
-</config>
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/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/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/define_vertical_grid.py b/testing_and_setup/compass/ocean/global_ocean/QU240/init/define_vertical_grid.py
new file mode 120000
index 0000000000..4e1cac8e2b
--- /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_64_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/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 @@
 <config case="base_mesh">
 
 	<add_link source_path="mpas_model" source="testing_and_setup/compass/ocean/jigsaw_to_MPAS" dest="jigsaw_to_MPAS"/>
-	<add_link source_path="script_test_dir" source="." dest="define_base_mesh"/>
+	<copy_file source_path="script_test_dir" source="define_base_mesh.py" dest="define_base_mesh.py"/>
 
 	<run_script name="run.py">
 		<step executable="python">
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..20a429653e 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,8 @@
 	<add_link source="../culled_mesh/culled_mesh.nc" dest="mesh.nc"/>
 	<add_link source="../culled_mesh/culled_graph.info" dest="graph.info"/>
 	<add_link source="../culled_mesh/critical_passages_mask_final.nc" dest="critical_passages.nc"/>
-	<add_link source_path="initial_condition_database" source="layer_depth.80Layer.180619.nc" dest="layer_depth.nc"/>
+	<copy_file source_path="script_test_dir" source="define_vertical_grid.py" dest="define_vertical_grid.py"/>
+	<add_link source_path="script_core_dir" source="scripts/vertical_grid/make_vertical_grid.py" dest="make_vertical_grid.py"/>
 	<add_link source_path="initial_condition_database" source="PTemp.Jan_p3.filled.mpas100levs.160127.nc" dest="temperature.nc"/>
 	<add_link source_path="initial_condition_database" source="Salt.Jan_p3.noBlackCaspian.filled.mpas100levs.160127.nc" dest="salinity.nc"/>
 	<add_link source_path="initial_condition_database" source="windStress.ncep_1958-2000avg.interp3600x2431.151106.nc" dest="wind_stress.nc"/>
@@ -38,7 +39,9 @@
 	<namelist name="namelist.ocean" mode="init">
 		<template file="template_initial_state.xml" path_base="script_configuration_dir"/>
 		<template file="template_critical_passages.xml" path_base="script_core_dir" path="global_ocean"/>
-		<option name="config_global_ocean_depth_file">'layer_depth.nc'</option>
+		<option name="config_global_ocean_depth_file">'vertical_grid.nc'</option>
+		<option name="config_global_ocean_depth_dimname">'nVertLevels'</option>
+		<option name="config_global_ocean_depth_varname">'refMidDepth'</option>
 		<option name="config_global_ocean_minimum_depth">10</option>
 		<option name="config_global_ocean_deepen_critical_passages">.false.</option>
 	</namelist>
@@ -52,6 +55,9 @@
 		<step executable="gpmetis">
 			<argument flag="graph.info">8</argument>
 		</step>
+		<step executable="./define_vertical_grid.py">
+		</step>
 		<model_run procs="8" threads="1" namelist="namelist.ocean" streams="streams.ocean"/>
-    </run_script>
+	</run_script>
+
 </config>
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..f410a9920f 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
@@ -28,22 +28,25 @@
 	<add_link source="../culled_mesh/culled_mesh.nc" dest="mesh.nc"/>
 	<add_link source="../culled_mesh/culled_graph.info" dest="graph.info"/>
 	<add_link source="../culled_mesh/critical_passages_mask_final.nc" dest="critical_passages.nc"/>
-	<add_link source="../culled_mesh/land_ice_mask.nc" dest="land_ice_mask.nc"/>
-	<add_link source_path="initial_condition_database" source="layer_depth.80Layer.180619.nc" dest="layer_depth.nc"/>
+	<copy_file source_path="script_test_dir" source="define_vertical_grid.py" dest="define_vertical_grid.py"/>
+	<add_link source_path="script_core_dir" source="scripts/vertical_grid/make_vertical_grid.py" dest="make_vertical_grid.py"/>
 	<add_link source_path="initial_condition_database" source="PTemp.Jan_p3.filled.mpas100levs.160127.nc" dest="temperature.nc"/>
 	<add_link source_path="initial_condition_database" source="Salt.Jan_p3.noBlackCaspian.filled.mpas100levs.160127.nc" dest="salinity.nc"/>
 	<add_link source_path="initial_condition_database" source="windStress.ncep_1958-2000avg.interp3600x2431.151106.nc" dest="wind_stress.nc"/>
 	<add_link source_path="bathymetry_database" source="BedMachineAntarctica_and_GEBCO_2019_0.05_degree.200128.nc" dest="topography.nc"/>
 	<add_link source_path="initial_condition_database" source="chlorophyllA_monthly_averages_1deg.151201.nc" dest="swData.nc"/>
+	<add_link source="../culled_mesh/land_ice_mask.nc" dest="land_ice_mask.nc"/>
 
 	<namelist name="namelist.ocean" mode="init">
 		<template file="template_initial_state.xml" path_base="script_configuration_dir"/>
 		<template file="template_critical_passages.xml" path_base="script_core_dir" path="global_ocean"/>
 		<template file="template_init_with_land_ice.xml" path_base="script_configuration_dir"/>
-		<option name="config_iterative_init_variable">'landIcePressure_from_top_density'</option>
-		<option name="config_global_ocean_depth_file">'layer_depth.nc'</option>
+		<option name="config_global_ocean_depth_file">'vertical_grid.nc'</option>
+		<option name="config_global_ocean_depth_dimname">'nVertLevels'</option>
+		<option name="config_global_ocean_depth_varname">'refMidDepth'</option>
 		<option name="config_global_ocean_minimum_depth">10</option>
 		<option name="config_global_ocean_deepen_critical_passages">.false.</option>
+		<option name="config_iterative_init_variable">'landIcePressure_from_top_density'</option>
 	</namelist>
 
 	<streams name="streams.ocean" keep="immutable" mode="init">
@@ -54,9 +57,11 @@
 
 	<run_script name="run.py">
 		<step executable="gpmetis">
-			<argument flag="graph.info">4</argument>
+			<argument flag="graph.info">8</argument>
+		</step>
+		<step executable="./define_vertical_grid.py">
 		</step>
-		<model_run procs="4" threads="1" namelist="namelist.ocean" streams="streams.ocean"/>
-    </run_script>
+		<model_run procs="8" threads="1" namelist="namelist.ocean" streams="streams.ocean"/>
+	</run_script>
 
 </config>
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/make_vertical_grid.py b/testing_and_setup/compass/ocean/scripts/make_vertical_grid.py
deleted file mode 100755
index 61459ece37..0000000000
--- a/testing_and_setup/compass/ocean/scripts/make_vertical_grid.py
+++ /dev/null
@@ -1,103 +0,0 @@
-def dz_z(z,Hmax,epsilon,dzmax):
-    import numpy as np
-    return dzmax*np.tanh(-z*np.pi/Hmax) + epsilon
-
-from netCDF4 import Dataset
-import numpy as np
-"""
-Write vertical grid to a netcdf file
-"""
-def create_vertical_grid(bottom_depth, nz, dz1, dz2, maxit=1000,
-    epsilon=1e-2, outFile='MPAS-Ocean_vertical_grid.nc'):
-
-  # open a new netCDF file for writing.
-  ncfile = Dataset(outFile,'w')
-  # create the depth_t dimension.
-  ncfile.createDimension('dim_layer',nz)
-
-  refBottomDepth = ncfile.createVariable('refBottomDepth',np.dtype('float64').char,('dim_layer'))
-  refMidDepth = ncfile.createVariable('refMidDepth',np.dtype('float64').char,('dim_layer'))
-  refLayerThickness = ncfile.createVariable('refLayerThickness',np.dtype('float64').char,('dim_layer'))
-
-  Hmax = bottom_depth
-  nLayers = 0
-
-  layerThickness = np.zeros(nz)
-  dz = [epsilon]
-  z = [0]
-  count=0
-
-  while nLayers != nz and count < maxit:
-    zval = -epsilon
-    dz = [epsilon]
-    z = [0]
-    nLayers = 0
-    while zval > -Hmax:
-        difference = dz_z(zval,Hmax,epsilon,dz2) - zval
-        while abs(difference) > 0.3:
-            zval -= epsilon
-            difference = dz_z(zval,Hmax,epsilon,dz2) -z[nLayers] + zval
-        z.append(zval)
-        dz.append(dz_z(zval,Hmax,epsilon,dz2))
-        nLayers += 1
-        zval -= epsilon
-
-    dz_arr = np.asarray(dz)
-    ind = abs(dz_arr - dz1).argmin()
-
-    dztemp = dz_arr[ind:]
-    nLayers = len(dztemp)
-    change = nz-nLayers
-    dz2 -= float(change)
-    count += 1
-
-  if count >= maxit:
-    print("Error: grid did not converge, adjust parameters")
-  else:
-    layerThickness[:nLayers] = dztemp
-
-    nVertLevels = nz 
-    botDepth = np.zeros(nVertLevels)
-    midDepth = np.zeros(nVertLevels)
-    botDepth[0] = layerThickness[0]
-    midDepth[0] = 0.5*layerThickness[0]
-
-    for i in range(1,nVertLevels):
-      botDepth[i] = botDepth[i-1] + layerThickness[i]
-      midDepth[i] = midDepth[i-1] + 0.5*(layerThickness[i] + layerThickness[i-1])
-
-    refBottomDepth[:] = botDepth
-    refMidDepth[:] = midDepth
-    refLayerThickness[:] = layerThickness[:nVertLevels]
-
-    ncfile.close()
-
-if __name__ == "__main__":
-  import argparse
-  parser = argparse.ArgumentParser(description=__doc__,
-                                  formatter_class=argparse.RawTextHelpFormatter)
-  parser.add_argument("-o", "--output_file_pattern", dest="output_filename_pattern",
-      default='MPAS-Ocean_vertical_grid.nc',
-      help="MPAS Filename pattern for file output of vertical grid.", metavar="NAME")
-  parser.add_argument("-bd", "--bottom_depth", dest="bottom_depth",
-      help="bottom depth for the chosen vertical coordinate", 
-      type=float,required=True)
-  parser.add_argument("-nz", "--num_vert_levels", dest="nz",
-      help="Number of vertical levels for the grid", type=int,
-      required=True)
-  parser.add_argument("-dz1", "--layer1_thickness", dest="dz1",
-      help="Target thickness of the first layer", 
-      type=float, required=True)
-  parser.add_argument("-dz2", "--maxLayer_thickness", dest="dz2",
-      help="Target maximum thickness in column", 
-      type=float, required=True)
-  parser.add_argument("-eps", "--error_tolerance", dest="epsilon",
-      help="Threshold for iterations", type=float, default=1e-2)
-  parser.add_argument("-maxit", "--max_iterations", dest="maxit",
-      help="maximum number of iterations for grid convergences",
-      type=int, default=1000)
-  args = parser.parse_args()
-
-  create_vertical_grid(args.bottom_depth,args.nz,args.dz1,
-      args.dz2,maxit=args.maxit,epsilon=args.epsilon,
-      outFile=args.output_filename_pattern)
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..f3d1906933
--- /dev/null
+++ b/testing_and_setup/compass/ocean/scripts/vertical_grid/define_vertical_grid_64_layers.py
@@ -0,0 +1,14 @@
+#!/usr/bin/env python
+
+import make_vertical_grid
+
+# Standard 64 layer vertical grid
+make_vertical_grid.create_vertical_grid(
+     5500,            # bottom_depth
+     64,              # num_vert_levels
+     2,               # layer1_thickness
+     200,             # maxLayer_thickness
+     True,            # plot_vertical_grid
+     maxit=1000,
+     epsilon=1.0e-2,
+     outFile= 'vertical_grid.nc')# output_file_name  
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..832690a59c
--- /dev/null
+++ b/testing_and_setup/compass/ocean/scripts/vertical_grid/make_vertical_grid.py
@@ -0,0 +1,198 @@
+#!/usr/bin/env python
+'''
+This script creates the vertical grid for MPAS-Ocean and writes it to a netcdf file.
+'''
+# import modules
+# {{{
+from netCDF4 import Dataset
+import numpy as np
+import argparse
+# }}}
+
+
+def main():
+    # parser
+    # {{{
+    parser = argparse.ArgumentParser(
+        description=__doc__,
+        formatter_class=argparse.RawTextHelpFormatter)
+    parser.add_argument(
+        '-o', '--output_file_name', dest='output_filename_name',
+        default='MPAS-Ocean_vertical_grid.nc',
+        help='MPAS file name for output of vertical grid.',
+        metavar='NAME')
+    parser.add_argument(
+        '-p', '--plot_vertical_grid', dest='plot_vertical_grid',
+        action='store_true')
+    parser.add_argument(
+        '-bd', '--bottom_depth', dest='bottom_depth',
+        default=5000,
+        help='bottom depth for the chosen vertical coordinate [m]',
+        type=float)
+    parser.add_argument(
+        '-nz', '--num_vert_levels', dest='nz',
+        default=64,
+        help='Number of vertical levels for the grid',
+        type=int)
+    parser.add_argument(
+        '-dz1', '--layer1_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(
+        '-eps', '--error_tolerance', dest='epsilon',
+        default=1e-2,
+        help='Threshold for iterations', type=float)
+    parser.add_argument(
+        '-maxit',
+        '--max_iterations',
+        dest='maxit',
+        default=1000,
+        help='maximum number of iterations for grid convergences',
+        type=int)
+    args = parser.parse_args()
+
+    create_vertical_grid(args.bottom_depth, args.nz, args.dz1,
+                         args.dz2, args.plot_vertical_grid,
+                         maxit=args.maxit, epsilon=args.epsilon,
+                         outFile=args.output_filename_name)
+# }}}
+
+
+def create_vertical_grid(
+        bottom_depth,
+        nz,
+        dz1_in,
+        dz2_in,
+        plot_vertical_grid,
+        maxit,
+        epsilon,
+        outFile):
+    # {{{
+    print('Creating mesh with ', nz, ' layers...')
+    dz1 = dz1_in
+    dz2 = dz2_in
+    # 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'))
+
+    Hmax = bottom_depth
+    nLayers = 0
+
+    layerThickness = np.zeros(nz)
+    dz = [epsilon]
+    z = [0]
+    count = 0
+
+    while nLayers != nz and count < maxit:
+        zval = -epsilon
+        dz = [epsilon]
+        z = [0]
+        nLayers = 0
+        while zval > -Hmax:
+            difference = dz_z(zval, Hmax, epsilon, dz2) - zval
+            while abs(difference) > 0.3:
+                zval -= epsilon
+                difference = dz_z(zval, Hmax, epsilon, dz2) - z[nLayers] + zval
+            z.append(zval)
+            dz.append(dz_z(zval, Hmax, epsilon, dz2))
+            nLayers += 1
+            zval -= epsilon
+
+        dz_arr = np.asarray(dz)
+        ind = abs(dz_arr - dz1).argmin()
+
+        dztemp = dz_arr[ind:]
+        nLayers = len(dztemp)
+        change = nz - nLayers
+        dz2 -= float(change)
+        count += 1
+
+    layerThickness[:nLayers] = dztemp
+    nVertLevels = nz
+    botDepth = np.zeros(nVertLevels)
+    midDepth = np.zeros(nVertLevels)
+    botDepth[0] = layerThickness[0]
+    midDepth[0] = 0.5 * layerThickness[0]
+
+    for i in range(1, nVertLevels):
+        botDepth[i] = botDepth[i - 1] + layerThickness[i]
+        midDepth[i] = midDepth[i - 1] + 0.5 * \
+            (layerThickness[i] + layerThickness[i - 1])
+
+    if count >= maxit:
+        print('Error: grid did not converge, adjust parameters')
+    else:
+        refBottomDepth[:] = botDepth
+        refMidDepth[:] = midDepth
+        refLayerThickness[:] = layerThickness[:nVertLevels]
+        ncfile.close()
+
+    if plot_vertical_grid:
+        import matplotlib
+        import matplotlib.pyplot as plt
+        matplotlib.use('Agg')
+        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: ' + str(nz) + '\n' + \
+            'bottom depth requested:  ' + '{:8.2f}'.format(bottom_depth) + '\n' +  \
+            'bottom depth actual:     ' + '{:8.2f}'.format(np.amax(botDepth)) + '\n' +  \
+            'min thickness reqeusted: ' + '{:8.2f}'.format(dz1_in) + '\n' + \
+            'min thickness actual:    ' + '{:8.2f}'.format(np.amin(layerThickness[:])) + '\n' + \
+            'max thickness reqeusted: ' + '{:8.2f}'.format(dz2_in) + '\n' + \
+            '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 dz_z(z, Hmax, epsilon, dzmax):
+    return dzmax * np.tanh(-z * np.pi / Hmax) + epsilon
+
+
+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

From 888a13ec2996d5b5fe7f1dd17f33dd99f243b835 Mon Sep 17 00:00:00 2001
From: Xylar Asay-Davis <xylarstorm@gmail.com>
Date: Tue, 31 Mar 2020 16:40:02 +0200
Subject: [PATCH 3/9] Switch algorithm for finding layer thicknesses

This is a more straight-forward root finding approach using
scipy.  It finds a characteristic lengh scale `delta` over which
layer thickenesses vary from `dz1` a the surface to `dz2` at
`z --> -infinity` such that the sum of the layer thicknesses
hits the target bottom depth in the given number of layers.

Docstrings and formatting have been improved.
---
 .../define_vertical_grid_64_layers.py         |  13 +-
 .../vertical_grid/make_vertical_grid.py       | 261 +++++++++++-------
 2 files changed, 171 insertions(+), 103 deletions(-)

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
index f3d1906933..7392d14f55 100755
--- 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
@@ -1,14 +1,7 @@
 #!/usr/bin/env python
 
-import make_vertical_grid
+from make_vertical_grid import create_vertical_grid
 
 # Standard 64 layer vertical grid
-make_vertical_grid.create_vertical_grid(
-     5500,            # bottom_depth
-     64,              # num_vert_levels
-     2,               # layer1_thickness
-     200,             # maxLayer_thickness
-     True,            # plot_vertical_grid
-     maxit=1000,
-     epsilon=1.0e-2,
-     outFile= 'vertical_grid.nc')# output_file_name  
+create_vertical_grid(bottom_depth=5500, nz=64, dz1_in=2., dz2_in=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
index 832690a59c..6f1db75529 100755
--- 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
@@ -1,12 +1,17 @@
 #!/usr/bin/env python
-'''
-This script creates the vertical grid for MPAS-Ocean and writes it to a netcdf file.
-'''
+"""
+This script creates the vertical grid for MPAS-Ocean and writes it to a netcdf
+file.
+"""
 # import modules
 # {{{
 from netCDF4 import Dataset
 import numpy as np
 import argparse
+from scipy.optimize import root_scalar
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
 # }}}
 
 
@@ -26,7 +31,7 @@ def main():
         action='store_true')
     parser.add_argument(
         '-bd', '--bottom_depth', dest='bottom_depth',
-        default=5000,
+        default=5000.0,
         help='bottom depth for the chosen vertical coordinate [m]',
         type=float)
     parser.add_argument(
@@ -44,107 +49,81 @@ def main():
         default=250.0,
         help='Target maximum thickness in column [m]',
         type=float)
-    parser.add_argument(
-        '-eps', '--error_tolerance', dest='epsilon',
-        default=1e-2,
-        help='Threshold for iterations', type=float)
-    parser.add_argument(
-        '-maxit',
-        '--max_iterations',
-        dest='maxit',
-        default=1000,
-        help='maximum number of iterations for grid convergences',
-        type=int)
     args = parser.parse_args()
 
     create_vertical_grid(args.bottom_depth, args.nz, args.dz1,
                          args.dz2, args.plot_vertical_grid,
-                         maxit=args.maxit, epsilon=args.epsilon,
-                         outFile=args.output_filename_name)
+                         outfile=args.output_filename_name)
 # }}}
 
 
-def create_vertical_grid(
-        bottom_depth,
-        nz,
-        dz1_in,
-        dz2_in,
-        plot_vertical_grid,
-        maxit,
-        epsilon,
-        outFile):
+def create_vertical_grid(bottom_depth=5000.0, nz=64, dz1_in=2.0, dz2_in=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]
+
+    nz : int, optional
+        Number of vertical levels for the grid
+
+    dz1_in : float, optional
+        Target thickness of the first layer [m]
+
+    dz2_in : 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 ', nz, ' layers...')
     dz1 = dz1_in
     dz2 = dz2_in
     # open a new netCDF file for writing.
-    ncfile = Dataset(outFile, 'w')
+    ncfile = Dataset(outfile, 'w')
     # create the depth_t dimension.
     ncfile.createDimension('nVertLevels', nz)
 
     refBottomDepth = ncfile.createVariable(
-        'refBottomDepth', np.dtype('float64').char, ('nVertLevels'))
+        'refBottomDepth', np.dtype('float64').char, ('nVertLevels',))
     refMidDepth = ncfile.createVariable(
-        'refMidDepth', np.dtype('float64').char, ('nVertLevels'))
+        'refMidDepth', np.dtype('float64').char, ('nVertLevels',))
     refLayerThickness = ncfile.createVariable(
-        'refLayerThickness', np.dtype('float64').char, ('nVertLevels'))
-
-    Hmax = bottom_depth
-    nLayers = 0
-
-    layerThickness = np.zeros(nz)
-    dz = [epsilon]
-    z = [0]
-    count = 0
-
-    while nLayers != nz and count < maxit:
-        zval = -epsilon
-        dz = [epsilon]
-        z = [0]
-        nLayers = 0
-        while zval > -Hmax:
-            difference = dz_z(zval, Hmax, epsilon, dz2) - zval
-            while abs(difference) > 0.3:
-                zval -= epsilon
-                difference = dz_z(zval, Hmax, epsilon, dz2) - z[nLayers] + zval
-            z.append(zval)
-            dz.append(dz_z(zval, Hmax, epsilon, dz2))
-            nLayers += 1
-            zval -= epsilon
-
-        dz_arr = np.asarray(dz)
-        ind = abs(dz_arr - dz1).argmin()
-
-        dztemp = dz_arr[ind:]
-        nLayers = len(dztemp)
-        change = nz - nLayers
-        dz2 -= float(change)
-        count += 1
-
-    layerThickness[:nLayers] = dztemp
+        '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 = np.zeros(nVertLevels)
-    midDepth = np.zeros(nVertLevels)
-    botDepth[0] = layerThickness[0]
-    midDepth[0] = 0.5 * layerThickness[0]
-
-    for i in range(1, nVertLevels):
-        botDepth[i] = botDepth[i - 1] + layerThickness[i]
-        midDepth[i] = midDepth[i - 1] + 0.5 * \
-            (layerThickness[i] + layerThickness[i - 1])
-
-    if count >= maxit:
-        print('Error: grid did not converge, adjust parameters')
-    else:
-        refBottomDepth[:] = botDepth
-        refMidDepth[:] = midDepth
-        refLayerThickness[:] = layerThickness[:nVertLevels]
-        ncfile.close()
+    botDepth = -z[1:]
+    midDepth = -0.5*(z[0:-1] + z[1:])
+
+    refBottomDepth[:] = botDepth
+    refMidDepth[:] = midDepth
+    refLayerThickness[:] = layerThickness[:nVertLevels]
+    ncfile.close()
 
     if plot_vertical_grid:
-        import matplotlib
-        import matplotlib.pyplot as plt
-        matplotlib.use('Agg')
         fig = plt.figure()
         fig.set_size_inches(16.0, 8.0)
         zInd = np.arange(1, nVertLevels + 1)
@@ -171,13 +150,13 @@ def create_vertical_grid(
         plt.grid()
 
         txt = \
-            'number layers: ' + str(nz) + '\n' + \
-            'bottom depth requested:  ' + '{:8.2f}'.format(bottom_depth) + '\n' +  \
-            'bottom depth actual:     ' + '{:8.2f}'.format(np.amax(botDepth)) + '\n' +  \
-            'min thickness reqeusted: ' + '{:8.2f}'.format(dz1_in) + '\n' + \
-            'min thickness actual:    ' + '{:8.2f}'.format(np.amin(layerThickness[:])) + '\n' + \
-            'max thickness reqeusted: ' + '{:8.2f}'.format(dz2_in) + '\n' + \
-            'max thickness actual:    ' + '{:8.2f}'.format(np.amax(layerThickness[:]))
+            '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(dz1_in) + \
+            'min thickness actual:    {:8.2f}\n'.format(np.amin(layerThickness[:])) + \
+            'max thickness reqeusted: {:8.2f}\n'.format(dz2_in) + \
+            'max thickness actual:    {:8.2f}'.format(np.amax(layerThickness[:]))
         print(txt)
         plt.subplot(2, 2, 4)
         plt.text(0, 0, txt, fontsize=12)
@@ -187,8 +166,104 @@ def create_vertical_grid(
 # }}}
 
 
-def dz_z(z, Hmax, epsilon, dzmax):
-    return dzmax * np.tanh(-z * np.pi / Hmax) + epsilon
+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__':

From 1b7023110913133a2025c037b09bf47981a5126d Mon Sep 17 00:00:00 2001
From: Mark Petersen <mpetersen@lanl.gov>
Date: Wed, 1 Apr 2020 13:32:40 -0600
Subject: [PATCH 4/9] Add new vertical grid to all cases

remove link for config_initial_state.xml
copy ocean/develop version of config_initial_state.xml in
---
 .../init/config_initial_state_60_layers.xml   |  1 +
 .../init/config_initial_state_60_layers.xml   |  1 +
 .../config_initial_state_60_layers.xml        |  2 +-
 .../config_files_ISC/config_initial_state.xml | 10 +--
 .../define_vertical_grid_64_layers.py         |  9 ++-
 .../define_vertical_grid_80_layers.py         | 12 ++++
 .../vertical_grid/make_vertical_grid.py       | 63 ++++++++++---------
 7 files changed, 59 insertions(+), 39 deletions(-)
 create mode 120000 testing_and_setup/compass/ocean/global_ocean/EC60to30/init/config_initial_state_60_layers.xml
 create mode 120000 testing_and_setup/compass/ocean/global_ocean/QU240/init/config_initial_state_60_layers.xml
 create mode 100755 testing_and_setup/compass/ocean/scripts/vertical_grid/define_vertical_grid_80_layers.py

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/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/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 @@
 <?xml version="1.0"?>
-<config case="initial_state">
+<config case="initial_state_60_layer">
 	<get_file dest_path="initial_condition_database" file_name="Salinity.01.filled.60levels.PHC.151106.nc">
 		<mirror protocol="wget" url="https://web.lcrc.anl.gov/public/e3sm/mpas_standalonedata/mpas-ocean/initial_condition_database/"/>
 	</get_file>
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 f410a9920f..9fbc9cb23d 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 @@
 <?xml version="1.0"?>
 <config case="initial_state">
-	<get_file dest_path="initial_condition_database" file_name="layer_depth.80Layer.180619.nc">
-		<mirror protocol="wget" url="https://web.lcrc.anl.gov/public/e3sm/mpas_standalonedata/mpas-ocean/initial_condition_database/"/>
-	</get_file>
-
 	<get_file dest_path="initial_condition_database" file_name="PTemp.Jan_p3.filled.mpas100levs.160127.nc">
 		<mirror protocol="wget" url="https://web.lcrc.anl.gov/public/e3sm/mpas_standalonedata/mpas-ocean/initial_condition_database/"/>
 	</get_file>
@@ -28,25 +24,25 @@
 	<add_link source="../culled_mesh/culled_mesh.nc" dest="mesh.nc"/>
 	<add_link source="../culled_mesh/culled_graph.info" dest="graph.info"/>
 	<add_link source="../culled_mesh/critical_passages_mask_final.nc" dest="critical_passages.nc"/>
-	<copy_file source_path="script_test_dir" source="define_vertical_grid.py" dest="define_vertical_grid.py"/>
+	<add_link source="../culled_mesh/land_ice_mask.nc" dest="land_ice_mask.nc"/>
 	<add_link source_path="script_core_dir" source="scripts/vertical_grid/make_vertical_grid.py" dest="make_vertical_grid.py"/>
+	<copy_file source_path="script_test_dir" source="define_vertical_grid.py" dest="define_vertical_grid.py"/>
 	<add_link source_path="initial_condition_database" source="PTemp.Jan_p3.filled.mpas100levs.160127.nc" dest="temperature.nc"/>
 	<add_link source_path="initial_condition_database" source="Salt.Jan_p3.noBlackCaspian.filled.mpas100levs.160127.nc" dest="salinity.nc"/>
 	<add_link source_path="initial_condition_database" source="windStress.ncep_1958-2000avg.interp3600x2431.151106.nc" dest="wind_stress.nc"/>
 	<add_link source_path="bathymetry_database" source="BedMachineAntarctica_and_GEBCO_2019_0.05_degree.200128.nc" dest="topography.nc"/>
 	<add_link source_path="initial_condition_database" source="chlorophyllA_monthly_averages_1deg.151201.nc" dest="swData.nc"/>
-	<add_link source="../culled_mesh/land_ice_mask.nc" dest="land_ice_mask.nc"/>
 
 	<namelist name="namelist.ocean" mode="init">
 		<template file="template_initial_state.xml" path_base="script_configuration_dir"/>
 		<template file="template_critical_passages.xml" path_base="script_core_dir" path="global_ocean"/>
 		<template file="template_init_with_land_ice.xml" path_base="script_configuration_dir"/>
+		<option name="config_iterative_init_variable">'landIcePressure_from_top_density'</option>
 		<option name="config_global_ocean_depth_file">'vertical_grid.nc'</option>
 		<option name="config_global_ocean_depth_dimname">'nVertLevels'</option>
 		<option name="config_global_ocean_depth_varname">'refMidDepth'</option>
 		<option name="config_global_ocean_minimum_depth">10</option>
 		<option name="config_global_ocean_deepen_critical_passages">.false.</option>
-		<option name="config_iterative_init_variable">'landIcePressure_from_top_density'</option>
 	</namelist>
 
 	<streams name="streams.ocean" keep="immutable" mode="init">
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
index 7392d14f55..51a00d73f1 100755
--- 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
@@ -3,5 +3,10 @@
 from make_vertical_grid import create_vertical_grid
 
 # Standard 64 layer vertical grid
-create_vertical_grid(bottom_depth=5500, nz=64, dz1_in=2., dz2_in=200.,
-                     plot_vertical_grid=True, outfile= 'vertical_grid.nc')
+create_vertical_grid(
+    num_vert_levels=64,
+    bottom_depth=6000,
+    min_layer_thickness=2.,
+    max_layer_thickness=200.,
+    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..ff9b204a22
--- /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=6500,
+    min_layer_thickness=1.0,
+    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
index 6f1db75529..a663f5613c 100755
--- 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
@@ -5,13 +5,13 @@
 """
 # 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')
-import matplotlib.pyplot as plt
 # }}}
 
 
@@ -22,25 +22,17 @@ def main():
         description=__doc__,
         formatter_class=argparse.RawTextHelpFormatter)
     parser.add_argument(
-        '-o', '--output_file_name', dest='output_filename_name',
-        default='MPAS-Ocean_vertical_grid.nc',
-        help='MPAS file name for output of vertical grid.',
-        metavar='NAME')
-    parser.add_argument(
-        '-p', '--plot_vertical_grid', dest='plot_vertical_grid',
-        action='store_true')
+        '-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(
-        '-nz', '--num_vert_levels', dest='nz',
-        default=64,
-        help='Number of vertical levels for the grid',
-        type=int)
-    parser.add_argument(
-        '-dz1', '--layer1_thickness', dest='dz1',
+        '-dz1', '--min_layer_thickness', dest='dz1',
         default=2.0,
         help='Target thickness of the first layer [m]',
         type=float)
@@ -49,6 +41,14 @@ def main():
         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_filename_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,
@@ -57,9 +57,13 @@ def main():
 # }}}
 
 
-def create_vertical_grid(bottom_depth=5000.0, nz=64, dz1_in=2.0, dz2_in=250.0,
-                         plot_vertical_grid=False,
-                         outfile='MPAS-Ocean_vertical_grid.nc'):
+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
@@ -71,13 +75,13 @@ def create_vertical_grid(bottom_depth=5000.0, nz=64, dz1_in=2.0, dz2_in=250.0,
     bottom_depth : float, optional
         bottom depth for the chosen vertical coordinate [m]
 
-    nz : int, optional
+    num_vert_levels : int, optional
         Number of vertical levels for the grid
 
-    dz1_in : float, optional
+    min_layer_thickness : float, optional
         Target thickness of the first layer [m]
 
-    dz2_in : float, optional
+    max_layer_thickness : float, optional
         Target maximum thickness in column [m]
 
     plot_vertical_grid : bool, optional
@@ -87,9 +91,10 @@ def create_vertical_grid(bottom_depth=5000.0, nz=64, dz1_in=2.0, dz2_in=250.0,
         MPAS file name for output of vertical grid
     """
 
-    print('Creating mesh with ', nz, ' layers...')
-    dz1 = dz1_in
-    dz2 = dz2_in
+    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.
@@ -109,14 +114,14 @@ def create_vertical_grid(bottom_depth=5000.0, nz=64, dz1_in=2.0, dz2_in=250.0,
     # 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],
+                      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:])
+    midDepth = -0.5 * (z[0:-1] + z[1:])
 
     refBottomDepth[:] = botDepth
     refMidDepth[:] = midDepth
@@ -153,9 +158,9 @@ def create_vertical_grid(bottom_depth=5000.0, nz=64, dz1_in=2.0, dz2_in=250.0,
             '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(dz1_in) + \
+            '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(dz2_in) + \
+            '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)
@@ -232,10 +237,10 @@ def cumsum_z(delta, nz, dz1, dz2):
         elements)
     """
     dz = np.zeros(nz)
-    z = np.zeros(nz+1)
+    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]
+        z[zindex + 1] = z[zindex] - dz[zindex]
     return dz, z
 
 

From e6a6aad0df824a5f1945c78ebb23ed6ee4c305da Mon Sep 17 00:00:00 2001
From: Mark Petersen <mpetersen@lanl.gov>
Date: Wed, 1 Apr 2020 13:51:48 -0600
Subject: [PATCH 5/9] Remove EC60to30wISC/init/config_initial_state.xml

---
 .../init/config_initial_state.xml             | 60 -------------------
 1 file changed, 60 deletions(-)
 delete mode 100644 testing_and_setup/compass/ocean/global_ocean/EC60to30wISC/init/config_initial_state.xml

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 @@
-<?xml version="1.0"?>
-<config case="initial_state">
-	<get_file dest_path="initial_condition_database" file_name="PotentialTemperature.01.filled.60levels.PHC.151106.nc">
-		<mirror protocol="wget" url="https://web.lcrc.anl.gov/public/e3sm/mpas_standalonedata/mpas-ocean/initial_condition_database/"/>
-	</get_file>
-
-	<get_file dest_path="initial_condition_database" file_name="Salinity.01.filled.60levels.PHC.151106.nc">
-		<mirror protocol="wget" url="https://web.lcrc.anl.gov/public/e3sm/mpas_standalonedata/mpas-ocean/initial_condition_database/"/>
-	</get_file>
-
-	<get_file dest_path="initial_condition_database" file_name="windStress.ncep_1958-2000avg.interp3600x2431.151106.nc">
-		<mirror protocol="wget" url="https://web.lcrc.anl.gov/public/e3sm/mpas_standalonedata/mpas-ocean/initial_condition_database/"/>
-	</get_file>
-
-	<get_file dest_path="bathymetry_database" file_name="BedMachineAntarctica_and_GEBCO_2019_0.05_degree.200128.nc">
-		<mirror protocol="wget" url="https://web.lcrc.anl.gov/public/e3sm/mpas_standalonedata/mpas-ocean/bathymetry_database/"/>
-	</get_file>
-
-	<get_file dest_path="initial_condition_database" file_name="chlorophyllA_monthly_averages_1deg.151201.nc">
-		<mirror protocol="wget" url="https://web.lcrc.anl.gov/public/e3sm/mpas_standalonedata/mpas-ocean/initial_condition_database/"/>
-	</get_file>
-
-
-	<add_link source="../culled_mesh/culled_mesh.nc" dest="mesh.nc"/>
-	<add_link source="../culled_mesh/culled_graph.info" dest="graph.info"/>
-	<add_link source="../culled_mesh/critical_passages_mask_final.nc" dest="critical_passages.nc"/>
-	<add_link source="../culled_mesh/land_ice_mask.nc" dest="land_ice_mask.nc"/>
-    <add_link source_path="initial_condition_database" source="Salinity.01.filled.60levels.PHC.151106.nc" dest="layer_depth.nc"/>
-	<add_link source_path="initial_condition_database" source="PotentialTemperature.01.filled.60levels.PHC.151106.nc" dest="temperature.nc"/>
-	<add_link source_path="initial_condition_database" source="Salinity.01.filled.60levels.PHC.151106.nc" dest="salinity.nc"/>
-	<add_link source_path="initial_condition_database" source="windStress.ncep_1958-2000avg.interp3600x2431.151106.nc" dest="wind_stress.nc"/>
-	<add_link source_path="bathymetry_database" source="BedMachineAntarctica_and_GEBCO_2019_0.05_degree.200128.nc" dest="topography.nc"/>
-	<add_link source_path="initial_condition_database" source="chlorophyllA_monthly_averages_1deg.151201.nc" dest="swData.nc"/>
-
-	<namelist name="namelist.ocean" mode="init">
-		<template file="template_initial_state.xml" path_base="script_configuration_dir"/>
-		<template file="template_critical_passages.xml" path_base="script_core_dir" path="global_ocean"/>
-		<template file="template_init_with_land_ice.xml" path_base="script_configuration_dir"/>
-		<option name="config_iterative_init_variable">'landIcePressure_from_top_density'</option>
-		<option name="config_global_ocean_depth_file">'layer_depth.nc'</option>
-		<option name="config_global_ocean_minimum_depth">10</option>
-		<option name="config_global_ocean_deepen_critical_passages">.false.</option>
-        	<option name="config_global_ocean_depth_conversion_factor">0.01</option>
-        	<option name="config_global_ocean_tracer_depth_conversion_factor">0.01</option>
-	</namelist>
-
-	<streams name="streams.ocean" keep="immutable" mode="init">
-		<template file="template_initial_state.xml" path_base="script_configuration_dir"/>
-		<template file="template_critical_passages.xml" path_base="script_core_dir" path="global_ocean"/>
-		<template file="template_init_with_land_ice.xml" path_base="script_configuration_dir"/>
-	</streams>
-
-	<run_script name="run.py">
-		<step executable="gpmetis">
-			<argument flag="graph.info">4</argument>
-		</step>
-		<model_run procs="4" threads="1" namelist="namelist.ocean" streams="streams.ocean"/>
-    </run_script>
-
-</config>

From d4202b83710c206344d16ae0f3bc81cc77641ad9 Mon Sep 17 00:00:00 2001
From: Mark Petersen <mpetersen@lanl.gov>
Date: Wed, 1 Apr 2020 13:52:43 -0600
Subject: [PATCH 6/9] Add link for EC60to30wISC/init/config_initial_state.xml

---
 .../global_ocean/EC60to30wISC/init/config_initial_state.xml      | 1 +
 1 file changed, 1 insertion(+)
 create mode 120000 testing_and_setup/compass/ocean/global_ocean/EC60to30wISC/init/config_initial_state.xml

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

From e378b8ed569a4fc0bcf9c269d37cd008df80cabe Mon Sep 17 00:00:00 2001
From: Mark Petersen <mpetersen@lanl.gov>
Date: Wed, 1 Apr 2020 14:05:29 -0600
Subject: [PATCH 7/9] Add old 60 layer to other cases

---
 .../ARM60to10/init/config_initial_state_60_layers.xml            | 1 +
 .../ocean/global_ocean/ARM60to10/init/define_vertical_grid.py    | 1 +
 .../ARM60to6/init/config_initial_state_60_layers.xml             | 1 +
 .../ocean/global_ocean/ARM60to6/init/define_vertical_grid.py     | 1 +
 .../global_ocean/CUSP12/init/config_initial_state_60_layers.xml  | 1 +
 .../ocean/global_ocean/CUSP12/init/define_vertical_grid.py       | 1 +
 .../global_ocean/CUSP8/init/config_initial_state_60_layers.xml   | 1 +
 .../ocean/global_ocean/CUSP8/init/define_vertical_grid.py        | 1 +
 .../ocean/global_ocean/SO60to10wISC/init/define_vertical_grid.py | 1 +
 9 files changed, 9 insertions(+)
 create mode 120000 testing_and_setup/compass/ocean/global_ocean/ARM60to10/init/config_initial_state_60_layers.xml
 create mode 120000 testing_and_setup/compass/ocean/global_ocean/ARM60to10/init/define_vertical_grid.py
 create mode 120000 testing_and_setup/compass/ocean/global_ocean/ARM60to6/init/config_initial_state_60_layers.xml
 create mode 120000 testing_and_setup/compass/ocean/global_ocean/ARM60to6/init/define_vertical_grid.py
 create mode 120000 testing_and_setup/compass/ocean/global_ocean/CUSP12/init/config_initial_state_60_layers.xml
 create mode 120000 testing_and_setup/compass/ocean/global_ocean/CUSP12/init/define_vertical_grid.py
 create mode 120000 testing_and_setup/compass/ocean/global_ocean/CUSP8/init/config_initial_state_60_layers.xml
 create mode 120000 testing_and_setup/compass/ocean/global_ocean/CUSP8/init/define_vertical_grid.py
 create mode 120000 testing_and_setup/compass/ocean/global_ocean/SO60to10wISC/init/define_vertical_grid.py

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_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/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

From 80630e9aa597848ec27ea3d24100826affc591f8 Mon Sep 17 00:00:00 2001
From: Mark Petersen <mpetersen@lanl.gov>
Date: Wed, 1 Apr 2020 17:07:53 -0600
Subject: [PATCH 8/9] Set 64 layer to 6000m max depth, varying from 2m to 210m

---
 .../global_ocean/QU240/init/define_vertical_grid.py  |  2 +-
 .../config_files/config_initial_state.xml            |  4 ++--
 .../config_files_ISC/config_initial_state.xml        |  4 ++--
 .../vertical_grid/define_vertical_grid_16_layers.py  | 12 ++++++++++++
 .../vertical_grid/define_vertical_grid_64_layers.py  |  2 +-
 .../vertical_grid/define_vertical_grid_80_layers.py  |  4 ++--
 6 files changed, 20 insertions(+), 8 deletions(-)
 create mode 100755 testing_and_setup/compass/ocean/scripts/vertical_grid/define_vertical_grid_16_layers.py

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
index 4e1cac8e2b..e92ea74842 120000
--- 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
@@ -1 +1 @@
-../../../scripts/vertical_grid/define_vertical_grid_64_layers.py
\ No newline at end of file
+../../../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/config_files/config_initial_state.xml b/testing_and_setup/compass/ocean/global_ocean/config_files/config_initial_state.xml
index 20a429653e..680ac155d7 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
@@ -53,11 +53,11 @@
 
 	<run_script name="run.py">
 		<step executable="gpmetis">
-			<argument flag="graph.info">8</argument>
+			<argument flag="graph.info">4</argument>
 		</step>
 		<step executable="./define_vertical_grid.py">
 		</step>
-		<model_run procs="8" threads="1" namelist="namelist.ocean" streams="streams.ocean"/>
+		<model_run procs="4" threads="1" namelist="namelist.ocean" streams="streams.ocean"/>
 	</run_script>
 
 </config>
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 9fbc9cb23d..6a9b88c5b2 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
@@ -53,11 +53,11 @@
 
 	<run_script name="run.py">
 		<step executable="gpmetis">
-			<argument flag="graph.info">8</argument>
+			<argument flag="graph.info">4</argument>
 		</step>
 		<step executable="./define_vertical_grid.py">
 		</step>
-		<model_run procs="8" threads="1" namelist="namelist.ocean" streams="streams.ocean"/>
+		<model_run procs="4" threads="1" namelist="namelist.ocean" streams="streams.ocean"/>
 	</run_script>
 
 </config>
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
index 51a00d73f1..36f9997a91 100755
--- 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
@@ -7,6 +7,6 @@
     num_vert_levels=64,
     bottom_depth=6000,
     min_layer_thickness=2.,
-    max_layer_thickness=200.,
+    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
index ff9b204a22..ebf3bc8ef7 100755
--- 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
@@ -5,8 +5,8 @@
 # Standard 80 layer vertical grid
 create_vertical_grid(
     num_vert_levels=80,
-    bottom_depth=6500,
-    min_layer_thickness=1.0,
+    bottom_depth=6000,
+    min_layer_thickness=1.,
     max_layer_thickness=200.,
     plot_vertical_grid=True,
     outfile='vertical_grid.nc')

From d44bb16714e066769756f1e28932b921e8a3517a Mon Sep 17 00:00:00 2001
From: Mark Petersen <mpetersen@lanl.gov>
Date: Wed, 1 Apr 2020 22:37:44 -0600
Subject: [PATCH 9/9] Add histogram plots of initial state

---
 .../config_files/config_initial_state.xml     |  3 +
 .../config_files_ISC/config_initial_state.xml |  5 +-
 .../config_ssh_adjustment.xml                 |  2 +-
 .../2km/Ocean0/config_adjust_ssh.xml          |  2 +-
 .../2km/Ocean1/config_adjust_ssh.xml          |  2 +-
 .../2km/Ocean2/config_adjust_ssh.xml          |  2 +-
 .../time_varying_Ocean0/config_adjust_ssh.xml |  2 +-
 .../5km/Ocean0/config_adjust_ssh.xml          |  2 +-
 .../scripts/{ => plots}/plot_globalStats.py   |  0
 .../ocean/scripts/plots/plot_initial_state.py | 80 +++++++++++++++++++
 .../vertical_grid/make_vertical_grid.py       |  4 +-
 11 files changed, 95 insertions(+), 9 deletions(-)
 rename testing_and_setup/compass/ocean/scripts/{ => plots}/plot_globalStats.py (100%)
 create mode 100755 testing_and_setup/compass/ocean/scripts/plots/plot_initial_state.py

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 680ac155d7..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
@@ -30,6 +30,7 @@
 	<add_link source="../culled_mesh/critical_passages_mask_final.nc" dest="critical_passages.nc"/>
 	<copy_file source_path="script_test_dir" source="define_vertical_grid.py" dest="define_vertical_grid.py"/>
 	<add_link source_path="script_core_dir" source="scripts/vertical_grid/make_vertical_grid.py" dest="make_vertical_grid.py"/>
+	<add_link source_path="script_core_dir" source="scripts/plots/plot_initial_state.py" dest="plot_initial_state.py"/>
 	<add_link source_path="initial_condition_database" source="PTemp.Jan_p3.filled.mpas100levs.160127.nc" dest="temperature.nc"/>
 	<add_link source_path="initial_condition_database" source="Salt.Jan_p3.noBlackCaspian.filled.mpas100levs.160127.nc" dest="salinity.nc"/>
 	<add_link source_path="initial_condition_database" source="windStress.ncep_1958-2000avg.interp3600x2431.151106.nc" dest="wind_stress.nc"/>
@@ -58,6 +59,8 @@
 		<step executable="./define_vertical_grid.py">
 		</step>
 		<model_run procs="4" threads="1" namelist="namelist.ocean" streams="streams.ocean"/>
+		<step executable="./plot_initial_state.py">
+		</step>
 	</run_script>
 
 </config>
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 6a9b88c5b2..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
@@ -24,9 +24,10 @@
 	<add_link source="../culled_mesh/culled_mesh.nc" dest="mesh.nc"/>
 	<add_link source="../culled_mesh/culled_graph.info" dest="graph.info"/>
 	<add_link source="../culled_mesh/critical_passages_mask_final.nc" dest="critical_passages.nc"/>
+	<copy_file source_path="script_test_dir" source="define_vertical_grid.py" dest="define_vertical_grid.py"/>
 	<add_link source="../culled_mesh/land_ice_mask.nc" dest="land_ice_mask.nc"/>
 	<add_link source_path="script_core_dir" source="scripts/vertical_grid/make_vertical_grid.py" dest="make_vertical_grid.py"/>
-	<copy_file source_path="script_test_dir" source="define_vertical_grid.py" dest="define_vertical_grid.py"/>
+	<add_link source_path="script_core_dir" source="scripts/plots/plot_initial_state.py" dest="plot_initial_state.py"/>
 	<add_link source_path="initial_condition_database" source="PTemp.Jan_p3.filled.mpas100levs.160127.nc" dest="temperature.nc"/>
 	<add_link source_path="initial_condition_database" source="Salt.Jan_p3.noBlackCaspian.filled.mpas100levs.160127.nc" dest="salinity.nc"/>
 	<add_link source_path="initial_condition_database" source="windStress.ncep_1958-2000avg.interp3600x2431.151106.nc" dest="wind_stress.nc"/>
@@ -58,6 +59,8 @@
 		<step executable="./define_vertical_grid.py">
 		</step>
 		<model_run procs="4" threads="1" namelist="namelist.ocean" streams="streams.ocean"/>
+		<step executable="./plot_initial_state.py">
+		</step>
 	</run_script>
 
 </config>
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 @@
 	<add_link source="../initial_state/initial_state.nc" dest="init0.nc"/>
 
 	<add_link source_path="script_core_dir" source="iterative_ssh_landIcePressure_scripts/iterate_init.py" dest="iterate_init.py"/>
-	<add_link source_path="script_core_dir" source="scripts/plot_globalStats.py" dest="plot_globalStats.py"/>
+	<add_link source_path="script_core_dir" source="scripts/plots/plot_globalStats.py" dest="plot_globalStats.py"/>
 
 	<namelist name="namelist.ocean" mode="forward">
 		<template file="template_forward.xml" path_base="script_configuration_dir"/>
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 @@
 
 	<add_link source_path="script_core_dir" source="iterative_ssh_landIcePressure_scripts/iterate_init.py" dest="iterate_init.py"/>
 	<add_link source_path="script_core_dir" source="iterative_ssh_landIcePressure_scripts/plot_cart_ssh_landIcePressure.py" dest="plot_cart_ssh_landIcePressure.py"/>
-	<add_link source_path="script_core_dir" source="scripts/plot_globalStats.py" dest="plot_globalStats.py"/>
+	<add_link source_path="script_core_dir" source="scripts/plots/plot_globalStats.py" dest="plot_globalStats.py"/>
 
 	<namelist name="namelist.ocean" mode="forward">
 		<template file="template_forward.xml" path_base="script_configuration_dir"/>
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 @@
 
 	<add_link source_path="script_core_dir" source="iterative_ssh_landIcePressure_scripts/iterate_init.py" dest="iterate_init.py"/>
 	<add_link source_path="script_core_dir" source="iterative_ssh_landIcePressure_scripts/plot_cart_ssh_landIcePressure.py" dest="plot_cart_ssh_landIcePressure.py"/>
-	<add_link source_path="script_core_dir" source="scripts/plot_globalStats.py" dest="plot_globalStats.py"/>
+	<add_link source_path="script_core_dir" source="scripts/plots/plot_globalStats.py" dest="plot_globalStats.py"/>
 
 	<namelist name="namelist.ocean" mode="forward">
 		<template file="template_forward.xml" path_base="script_configuration_dir"/>
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 @@
 
 	<add_link source_path="script_core_dir" source="iterative_ssh_landIcePressure_scripts/iterate_init.py" dest="iterate_init.py"/>
 	<add_link source_path="script_core_dir" source="iterative_ssh_landIcePressure_scripts/plot_cart_ssh_landIcePressure.py" dest="plot_cart_ssh_landIcePressure.py"/>
-	<add_link source_path="script_core_dir" source="scripts/plot_globalStats.py" dest="plot_globalStats.py"/>
+	<add_link source_path="script_core_dir" source="scripts/plots/plot_globalStats.py" dest="plot_globalStats.py"/>
 
 	<namelist name="namelist.ocean" mode="forward">
 		<template file="template_forward.xml" path_base="script_configuration_dir"/>
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 @@
 
 	<add_link source_path="script_core_dir" source="iterative_ssh_landIcePressure_scripts/iterate_init.py" dest="iterate_init.py"/>
 	<add_link source_path="script_core_dir" source="iterative_ssh_landIcePressure_scripts/plot_cart_ssh_landIcePressure.py" dest="plot_cart_ssh_landIcePressure.py"/>
-	<add_link source_path="script_core_dir" source="scripts/plot_globalStats.py" dest="plot_globalStats.py"/>
+	<add_link source_path="script_core_dir" source="scripts/plots/plot_globalStats.py" dest="plot_globalStats.py"/>
 
 	<namelist name="namelist.ocean" mode="forward">
 		<template file="template_forward.xml" path_base="script_configuration_dir"/>
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 @@
 
 	<add_link source_path="script_core_dir" source="iterative_ssh_landIcePressure_scripts/iterate_init.py" dest="iterate_init.py"/>
 	<add_link source_path="script_core_dir" source="iterative_ssh_landIcePressure_scripts/plot_cart_ssh_landIcePressure.py" dest="plot_cart_ssh_landIcePressure.py"/>
-	<add_link source_path="script_core_dir" source="scripts/plot_globalStats.py" dest="plot_globalStats.py"/>
+	<add_link source_path="script_core_dir" source="scripts/plots/plot_globalStats.py" dest="plot_globalStats.py"/>
 
 	<namelist name="namelist.ocean" mode="forward">
 		<template file="template_forward.xml" path_base="script_configuration_dir"/>
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/make_vertical_grid.py b/testing_and_setup/compass/ocean/scripts/vertical_grid/make_vertical_grid.py
index a663f5613c..e9a59d0295 100755
--- 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
@@ -45,7 +45,7 @@ def main():
         '-p', '--plot_vertical_grid', dest='plot_vertical_grid',
         action='store_true')
     parser.add_argument(
-        '-o', '--output_file_name', dest='output_filename_name',
+        '-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')
@@ -53,7 +53,7 @@ def main():
 
     create_vertical_grid(args.bottom_depth, args.nz, args.dz1,
                          args.dz2, args.plot_vertical_grid,
-                         outfile=args.output_filename_name)
+                         outfile=args.output_file_name)
 # }}}