diff --git a/ocean/merry_go_round/100Cells/default/add_initial_state.py b/ocean/merry_go_round/100Cells/default/add_initial_state.py
new file mode 100644
index 0000000000..2ba038c043
--- /dev/null
+++ b/ocean/merry_go_round/100Cells/default/add_initial_state.py
@@ -0,0 +1,238 @@
+#!/usr/bin/env python
+
+'''
+This script creates an initial condition file for MPAS-Ocean.
+
+Merry-go-round Test Case: steady state test case for horizontal
+and vertical advection.
+see:
+Calandrini, S., Pieper, K. and Gunzburger, M.D., 2020. Exponential time
+differencing for the tracer equations appearing in primitive equation
+ocean models. Computer Methods in Applied Mechanics and Engineering,
+365, p.113002.
+'''
+
+import os
+import shutil
+import numpy as np
+import xarray as xr
+from mpas_tools.io import write_netcdf
+import argparse
+import math
+import time
+verbose = True
+
+
+def main():
+ timeStart = time.time()
+ parser = argparse.ArgumentParser()
+ parser.add_argument('-i', '--input_file', dest='input_file',
+ default='base_mesh.nc',
+ help='Input file, containing base mesh'
+ )
+ parser.add_argument('-o', '--output_file', dest='output_file',
+ default='initial_state.nc',
+ help='Output file, containing initial variables'
+ )
+ parser.add_argument('-L', '--nVertLevels', dest='nVertLevels',
+ default=50,
+ help='Number of vertical levels'
+ )
+ nVertLevels = parser.parse_args().nVertLevels
+ ds = xr.open_dataset(parser.parse_args().input_file)
+
+# configuration settings
+ # reference vertical grid spacing
+ maxDepth = 500
+ vertical_coordinate = 'z' # sigma or z
+
+ # comment('obtain dimensions and mesh variables')
+ nCells = ds['nCells'].size
+ nEdges = ds['nEdges'].size
+ nVertices = ds['nVertices'].size
+
+ xCell = ds['xCell']
+ xEdge = ds['xEdge']
+ xVertex = ds['xVertex']
+ yCell = ds['yCell']
+ yEdge = ds['yEdge']
+ yVertex = ds['yVertex']
+ angleEdge = ds['angleEdge']
+ cellsOnEdge = ds['cellsOnEdge']
+ edgesOnCell = ds['edgesOnCell']
+
+ # Adjust coordinates so first edge is at zero in x and y
+ xOffset = xEdge.min()
+ xCell -= xOffset
+ xEdge -= xOffset
+ xVertex -= xOffset
+ yOffset = np.min(yEdge)
+ yCell -= yOffset
+ yEdge -= yOffset
+ yVertex -= yOffset
+
+ # x values for convenience
+ xCellMax = max(xCell)
+ xCellMin = min(xCell)
+ xCellMid = 0.5 * (xCellMax + xCellMin)
+ xEdgeMax = max(xEdge)
+
+ # initialize velocity field
+ u = np.zeros([1, nEdges, nVertLevels])
+
+ comment('create and initialize variables')
+ time1 = time.time()
+
+ varsZ = [ 'refLayerThickness', 'refBottomDepth', 'refZMid', 'vertCoordMovementWeights']
+ for var in varsZ:
+ globals()[var] = np.nan * np.ones(nVertLevels)
+
+ vars2D = ['ssh', 'bottomDepth', 'surfaceStress',
+ 'atmosphericPressure', 'boundaryLayerDepth']
+ for var in vars2D:
+ globals()[var] = np.nan * np.ones(nCells)
+ maxLevelCell = np.ones(nCells, dtype=np.int32)
+
+ vars3D = [ 'layerThickness','temperature', 'salinity',
+ 'zMid', 'density', 'tracer1', 'tracer2', 'tracer3']
+ for var in vars3D:
+ globals()[var] = np.nan * np.ones([1, nCells, nVertLevels])
+ restingThickness = np.nan * np.ones([nCells, nVertLevels])
+
+ # Note that this line shouldn't be required, but if layerThickness is
+ # initialized with nans, the simulation dies. It must multiply by a nan on
+ # a land cell on an edge, and then multiply by zero.
+ layerThickness[:] = -1e34
+
+ # equally spaced layers
+ refLayerThickness[:] = maxDepth / nVertLevels
+ refBottomDepth[0] = refLayerThickness[0]
+ refZMid[0] = -0.5 * refLayerThickness[0]
+ for k in range(1, nVertLevels):
+ refBottomDepth[k] = refBottomDepth[k - 1] + refLayerThickness[k]
+ refZMid[k] = -refBottomDepth[k - 1] - 0.5 * refLayerThickness[k]
+
+ # See surface height
+ ssh[:] = 0.0
+
+ # Compute maxLevelCell and layerThickness for z-level (variation only on top)
+ if (vertical_coordinate=='z'):
+ vertCoordMovementWeights[:] = 0.0
+ vertCoordMovementWeights[0] = 1.0
+ for iCell in range(0, nCells):
+ maxLevelCell[iCell] = nVertLevels - 1
+ bottomDepth[iCell] = refBottomDepth[nVertLevels - 1]
+ layerThickness[0, iCell, :] = refLayerThickness[:]
+ layerThickness[0, iCell, 0] += ssh[iCell]
+ restingThickness[:, :] = layerThickness[0, :, :]
+
+ # Compute maxLevelCell and layerThickness for sigma
+ #elif (vertical_coordinate=='sigma'):
+ # vertCoordMovementWeights[:] = 1.0
+ # maxLevelCell[:] = nVertLevels - 1 # maxLevelCell is zero-based within python code
+ # bottomDepth[:] = bottomDepthObserved[:]
+ # for iCell in range(0, nCells):
+ # restingThickness[iCell, :] = refLayerThickness[:]*bottomDepth[iCell]/maxDepth
+ # layerThickness[0, iCell, :] = refLayerThickness[:]*(ssh[iCell] + bottomDepth[iCell])/maxDepth
+
+ # Compute zMid (same, regardless of vertical coordinate)
+ for iCell in range(0, nCells):
+ k = maxLevelCell[iCell]
+ zMid[0, iCell, k] = -bottomDepth[iCell] + \
+ 0.5 * layerThickness[0, iCell, k]
+ for k in range(maxLevelCell[iCell] - 1, -1, -1):
+ zMid[0, iCell, k] = zMid[0, iCell, k + 1] + 0.5 * \
+ (layerThickness[0, iCell, k + 1] + layerThickness[0, iCell, k])
+
+ # initialize tracers
+ rho0 = 1000.0 # kg/m^3
+ rhoz = -2.0e-4 # kg/m^3/m in z
+ S0 = 35.0
+ r = 1
+
+ # linear equation of state
+ # rho = rho0 - alpha*(T-Tref) + beta*(S-Sref)
+ # set S=Sref
+ # T = Tref - (rho - rhoRef)/alpha
+ config_eos_linear_alpha = 0.2
+ config_eos_linear_beta = 0.8
+ config_eos_linear_Tref = 10.0
+ config_eos_linear_Sref = 35.0
+ config_eos_linear_densityref = 1000.0
+
+ for k in range(0, nVertLevels):
+ activeCells = k <= maxLevelCell
+ salinity[0, activeCells, k] = S0
+ density[0, activeCells, k] = rho0 + rhoz * zMid[0, activeCells, k]
+ # T = Tref - (rho - rhoRef)/alpha
+ #temperature[0, activeCells, k] = config_eos_linear_Tref \
+ # - (density[0, activeCells, k] - config_eos_linear_densityref) / \
+ # config_eos_linear_alpha
+ for iCell in range(0, nCells):
+ if (xCell[iCell] < xCellMid):
+ temperature[0, iCell, :] = 5
+ else:
+ temperature[0, iCell, :] = 30
+ for k in range(0, nVertLevels):
+ psi1 = 1.0 - (((xCell[iCell] - 0.5*xCellMax)**4)/((0.5*xCellMax)**4))
+ psi2 = 1.0 - (((zMid[0, iCell, k] + 0.5*maxDepth)**2)/((0.5*maxDepth)**2))
+ psi = psi1*psi2
+ tracer1[0,iCell,k] = 0.5*(1 + math.tanh( r*(2*psi - 1) ) )
+ tracer2[0,:,:] = 10
+ tracer3[0,:,:] = 20
+
+ # initial velocity on edges
+ ds['normalVelocity'] = (('Time', 'nEdges', 'nVertLevels',), np.zeros([1, nEdges, nVertLevels]))
+ normalVelocity = ds['normalVelocity']
+ for iEdge in range(0, nEdges):
+ cell1 = cellsOnEdge[iEdge,0] - 1
+ cell2 = cellsOnEdge[iEdge,1] - 1
+ for k in range(0, nVertLevels):
+ zMidEdge = 0.5*(zMid[0, cell1, k] + zMid[0, cell2, k])
+ dPsi = - (2.0*zMidEdge + maxDepth) / (0.5*maxDepth)**2
+ u[0, iEdge, k] = (1.0 - ((xEdge[iEdge] - 0.5*xCellMax)**4)/((0.5*(xCellMax-xCellMin))**4)) * dPsi
+ if ((cell1<0.5*nCells and cell1>0.25*nCells-1) and (cell2<0.5*nCells and cell2>0.25*nCells-1)):
+ u[0, iEdge, k] = (1.0 - (((xEdge[iEdge]-xCellMin) - 0.5*(xCellMax+xCellMin))**4)/((0.5*(xCellMax-xCellMin))**4)) * dPsi
+ if (cell1>0.75*nCells-1 and cell2>0.75*nCells-1):
+ u[0, iEdge, k] = (1.0 - (((xEdge[iEdge]-xCellMin) - 0.5*(xCellMax+xCellMin))**4)/((0.5*(xCellMax-xCellMin))**4)) * dPsi
+ normalVelocity[0, iEdge, :] = u[0, iEdge, :] * math.cos(angleEdge[iEdge])
+ if (xEdge[iEdge]<=xCellMin or xEdge[iEdge]>=xCellMax):
+ normalVelocity[0, iEdge, :] = 0.0
+
+ # Coriolis parameter
+ ds['fCell'] = (('nCells', 'nVertLevels',), np.zeros([nCells, nVertLevels]))
+ ds['fEdge'] = (('nEdges', 'nVertLevels',), np.zeros([nEdges, nVertLevels]))
+ ds['fVertex'] = (('nVertices', 'nVertLevels',), np.zeros([nVertices, nVertLevels]))
+
+ # surface fields
+ surfaceStress[:] = 0.0
+ atmosphericPressure[:] = 0.0
+ boundaryLayerDepth[:] = 0.0
+ print(' time: %f' % ((time.time() - time1)))
+
+ comment('finalize and write file')
+ time1 = time.time()
+ ds['maxLevelCell'] = (['nCells'], maxLevelCell + 1)
+ ds['restingThickness'] = (['nCells', 'nVertLevels'], restingThickness)
+ for var in varsZ:
+ ds[var] = (['nVertLevels'], globals()[var])
+ for var in vars2D:
+ ds[var] = (['nCells'], globals()[var])
+ for var in vars3D:
+ ds[var] = (['Time', 'nCells', 'nVertLevels'], globals()[var])
+ # If you prefer not to have NaN as the fill value, you should consider
+ # using mpas_tools.io.write_netcdf() instead
+ ds.to_netcdf('initial_state.nc', format='NETCDF3_64BIT_OFFSET')
+ # write_netcdf(ds,'initial_state.nc')
+ print(' time: %f' % ((time.time() - time1)))
+ print('Total time: %f' % ((time.time() - timeStart)))
+
+
+def comment(string):
+ if verbose:
+ print('*** ' + string)
+
+
+if __name__ == '__main__':
+ # If called as a primary module, run main
+ main()
diff --git a/ocean/merry_go_round/100Cells/default/build_base_mesh.py b/ocean/merry_go_round/100Cells/default/build_base_mesh.py
new file mode 100644
index 0000000000..ab6c2f60f8
--- /dev/null
+++ b/ocean/merry_go_round/100Cells/default/build_base_mesh.py
@@ -0,0 +1,26 @@
+#!/usr/bin/env python
+
+# This script was generated from setup_testcases.py as part of a config file
+
+import sys
+import os
+import shutil
+import glob
+import subprocess
+
+
+dev_null = open('/dev/null', 'w')
+
+# Run command is:
+# planar_hex --nx 1600 --ny 4 --dc 3e3 --nonperiodic_x -o planar_hex_mesh.nc
+subprocess.check_call(['planar_hex', '--nx', '100', '--ny', '4', '--dc',
+ '5', '--nonperiodic_x', '-o', 'planar_hex_mesh.nc'])
+
+# Run command is:
+# MpasCellCuller.x planar_hex_mesh.nc culled_mesh.nc
+subprocess.check_call(['MpasCellCuller.x', 'planar_hex_mesh.nc',
+ 'culled_mesh.nc'])
+
+# Run command is:
+# MpasMeshConverter.x culled_mesh.nc base_mesh.nc
+subprocess.check_call(['MpasMeshConverter.x', 'culled_mesh.nc', 'base_mesh.nc'])
diff --git a/ocean/merry_go_round/100Cells/default/config_base_mesh.xml b/ocean/merry_go_round/100Cells/default/config_base_mesh.xml
new file mode 100644
index 0000000000..3ac5df0612
--- /dev/null
+++ b/ocean/merry_go_round/100Cells/default/config_base_mesh.xml
@@ -0,0 +1,10 @@
+
+
+
+
+
+
+
+
+
+
diff --git a/ocean/merry_go_round/100Cells/default/config_driver.xml b/ocean/merry_go_round/100Cells/default/config_driver.xml
new file mode 100644
index 0000000000..457a7a033a
--- /dev/null
+++ b/ocean/merry_go_round/100Cells/default/config_driver.xml
@@ -0,0 +1,14 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/ocean/merry_go_round/100Cells/default/config_forward.xml b/ocean/merry_go_round/100Cells/default/config_forward.xml
new file mode 100644
index 0000000000..8b73cea440
--- /dev/null
+++ b/ocean/merry_go_round/100Cells/default/config_forward.xml
@@ -0,0 +1,23 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 4
+
+
+
+
diff --git a/ocean/merry_go_round/100Cells/default/config_initial_state.xml b/ocean/merry_go_round/100Cells/default/config_initial_state.xml
new file mode 100644
index 0000000000..3b9176f2c1
--- /dev/null
+++ b/ocean/merry_go_round/100Cells/default/config_initial_state.xml
@@ -0,0 +1,14 @@
+
+
+
+
+
+
+
+
+
+ base_mesh.nc
+ initial_state.nc
+
+
+
diff --git a/ocean/merry_go_round/100Cells/default/config_visualize.xml b/ocean/merry_go_round/100Cells/default/config_visualize.xml
new file mode 100644
index 0000000000..573eb219e8
--- /dev/null
+++ b/ocean/merry_go_round/100Cells/default/config_visualize.xml
@@ -0,0 +1,10 @@
+
+
+
+
+
+
+
+
+
+
diff --git a/ocean/merry_go_round/100Cells/default/fromSide.py b/ocean/merry_go_round/100Cells/default/fromSide.py
new file mode 100644
index 0000000000..eb5ee92dba
--- /dev/null
+++ b/ocean/merry_go_round/100Cells/default/fromSide.py
@@ -0,0 +1,62 @@
+import matplotlib
+import matplotlib.pyplot as plt
+from netCDF4 import Dataset
+import numpy as np
+matplotlib.use('Agg')
+
+fig = plt.gcf()
+fig.set_size_inches(8.0,10.0)
+
+ncfile = Dataset('../forward/output/output.0001-01-01_00.00.00.nc', 'r')
+var1 = ncfile.variables['velocityX']
+var2 = ncfile.variables['vertVelocityTop']
+var3 = ncfile.variables['tracer1']
+var11 = var1[:, 0:100, :]
+var22 = var2[:, 0:100, :]
+var33 = var3[:, 0:100, :]
+print(var11[1,:,0])
+for iCell in range(100):
+ sumVert = 0.0
+ for k in range(50):
+ sumVert = sumVert + var11[1, iCell, k]
+ print(sumVert)
+print(var22[1,0,:])
+for k in range(50):
+ sumVert = 0.0
+ for iCell in range(100):
+ sumVert = sumVert + var22[1, iCell, k]
+ print(sumVert)
+
+plt.title('merry-go-round')
+plt.subplot(2, 2, 1)
+ax = plt.imshow(var11[1, :, :].T)
+plt.jet()
+plt.xlabel('x, # of cells')
+#plt.ylabel('depth, m')
+plt.colorbar(ax, shrink=0.5)
+plt.title('horizontal velocity')
+
+plt.subplot(2, 2, 2)
+ax = plt.imshow(var22[1, :, :].T)
+plt.jet()
+plt.xlabel('x, # of cells')
+plt.colorbar(ax, shrink=0.5)
+plt.title('vertical velocity')
+
+plt.subplot(2, 2, 3)
+ax = plt.imshow(var33[0, :, :].T)
+plt.jet()
+plt.xlabel('x, # of cells')
+plt.colorbar(ax, shrink=0.5)
+plt.title('tracer1 at t=0')
+
+plt.subplot(2, 2, 4)
+ax = plt.imshow(var33[1, :, :].T)
+plt.jet()
+plt.xlabel('x, # of cells')
+plt.colorbar(ax, shrink=0.5)
+plt.title('tracer1 at 6h')
+
+ncfile.close()
+plt.savefig('section.png')
+
diff --git a/ocean/merry_go_round/200Cells/default/add_initial_state.py b/ocean/merry_go_round/200Cells/default/add_initial_state.py
new file mode 100644
index 0000000000..fa8866e796
--- /dev/null
+++ b/ocean/merry_go_round/200Cells/default/add_initial_state.py
@@ -0,0 +1,238 @@
+#!/usr/bin/env python
+
+'''
+This script creates an initial condition file for MPAS-Ocean.
+
+Merry-go-round Test Case: steady state test case for horizontal
+and vertical advection.
+see:
+Calandrini, S., Pieper, K. and Gunzburger, M.D., 2020. Exponential time
+differencing for the tracer equations appearing in primitive equation
+ocean models. Computer Methods in Applied Mechanics and Engineering,
+365, p.113002.
+'''
+
+import os
+import shutil
+import numpy as np
+import xarray as xr
+from mpas_tools.io import write_netcdf
+import argparse
+import math
+import time
+verbose = True
+
+
+def main():
+ timeStart = time.time()
+ parser = argparse.ArgumentParser()
+ parser.add_argument('-i', '--input_file', dest='input_file',
+ default='base_mesh.nc',
+ help='Input file, containing base mesh'
+ )
+ parser.add_argument('-o', '--output_file', dest='output_file',
+ default='initial_state.nc',
+ help='Output file, containing initial variables'
+ )
+ parser.add_argument('-L', '--nVertLevels', dest='nVertLevels',
+ default=100,
+ help='Number of vertical levels'
+ )
+ nVertLevels = parser.parse_args().nVertLevels
+ ds = xr.open_dataset(parser.parse_args().input_file)
+
+# configuration settings
+ # reference vertical grid spacing
+ maxDepth = 500
+ vertical_coordinate = 'z' # sigma or z
+
+ # comment('obtain dimensions and mesh variables')
+ nCells = ds['nCells'].size
+ nEdges = ds['nEdges'].size
+ nVertices = ds['nVertices'].size
+
+ xCell = ds['xCell']
+ xEdge = ds['xEdge']
+ xVertex = ds['xVertex']
+ yCell = ds['yCell']
+ yEdge = ds['yEdge']
+ yVertex = ds['yVertex']
+ angleEdge = ds['angleEdge']
+ cellsOnEdge = ds['cellsOnEdge']
+ edgesOnCell = ds['edgesOnCell']
+
+ # Adjust coordinates so first edge is at zero in x and y
+ xOffset = xEdge.min()
+ xCell -= xOffset
+ xEdge -= xOffset
+ xVertex -= xOffset
+ yOffset = np.min(yEdge)
+ yCell -= yOffset
+ yEdge -= yOffset
+ yVertex -= yOffset
+
+ # x values for convenience
+ xCellMax = max(xCell)
+ xCellMin = min(xCell)
+ xCellMid = 0.5 * (xCellMax + xCellMin)
+ xEdgeMax = max(xEdge)
+
+ # initialize velocity field
+ u = np.zeros([1, nEdges, nVertLevels])
+
+ comment('create and initialize variables')
+ time1 = time.time()
+
+ varsZ = [ 'refLayerThickness', 'refBottomDepth', 'refZMid', 'vertCoordMovementWeights']
+ for var in varsZ:
+ globals()[var] = np.nan * np.ones(nVertLevels)
+
+ vars2D = ['ssh', 'bottomDepth', 'surfaceStress',
+ 'atmosphericPressure', 'boundaryLayerDepth']
+ for var in vars2D:
+ globals()[var] = np.nan * np.ones(nCells)
+ maxLevelCell = np.ones(nCells, dtype=np.int32)
+
+ vars3D = [ 'layerThickness','temperature', 'salinity',
+ 'zMid', 'density', 'tracer1', 'tracer2', 'tracer3']
+ for var in vars3D:
+ globals()[var] = np.nan * np.ones([1, nCells, nVertLevels])
+ restingThickness = np.nan * np.ones([nCells, nVertLevels])
+
+ # Note that this line shouldn't be required, but if layerThickness is
+ # initialized with nans, the simulation dies. It must multiply by a nan on
+ # a land cell on an edge, and then multiply by zero.
+ layerThickness[:] = -1e34
+
+ # equally spaced layers
+ refLayerThickness[:] = maxDepth / nVertLevels
+ refBottomDepth[0] = refLayerThickness[0]
+ refZMid[0] = -0.5 * refLayerThickness[0]
+ for k in range(1, nVertLevels):
+ refBottomDepth[k] = refBottomDepth[k - 1] + refLayerThickness[k]
+ refZMid[k] = -refBottomDepth[k - 1] - 0.5 * refLayerThickness[k]
+
+ # See surface height
+ ssh[:] = 0.0
+
+ # Compute maxLevelCell and layerThickness for z-level (variation only on top)
+ if (vertical_coordinate=='z'):
+ vertCoordMovementWeights[:] = 0.0
+ vertCoordMovementWeights[0] = 1.0
+ for iCell in range(0, nCells):
+ maxLevelCell[iCell] = nVertLevels - 1
+ bottomDepth[iCell] = refBottomDepth[nVertLevels - 1]
+ layerThickness[0, iCell, :] = refLayerThickness[:]
+ layerThickness[0, iCell, 0] += ssh[iCell]
+ restingThickness[:, :] = layerThickness[0, :, :]
+
+ # Compute maxLevelCell and layerThickness for sigma
+ #elif (vertical_coordinate=='sigma'):
+ # vertCoordMovementWeights[:] = 1.0
+ # maxLevelCell[:] = nVertLevels - 1 # maxLevelCell is zero-based within python code
+ # bottomDepth[:] = bottomDepthObserved[:]
+ # for iCell in range(0, nCells):
+ # restingThickness[iCell, :] = refLayerThickness[:]*bottomDepth[iCell]/maxDepth
+ # layerThickness[0, iCell, :] = refLayerThickness[:]*(ssh[iCell] + bottomDepth[iCell])/maxDepth
+
+ # Compute zMid (same, regardless of vertical coordinate)
+ for iCell in range(0, nCells):
+ k = maxLevelCell[iCell]
+ zMid[0, iCell, k] = -bottomDepth[iCell] + \
+ 0.5 * layerThickness[0, iCell, k]
+ for k in range(maxLevelCell[iCell] - 1, -1, -1):
+ zMid[0, iCell, k] = zMid[0, iCell, k + 1] + 0.5 * \
+ (layerThickness[0, iCell, k + 1] + layerThickness[0, iCell, k])
+
+ # initialize tracers
+ rho0 = 1000.0 # kg/m^3
+ rhoz = -2.0e-4 # kg/m^3/m in z
+ S0 = 35.0
+ r = 1
+
+ # linear equation of state
+ # rho = rho0 - alpha*(T-Tref) + beta*(S-Sref)
+ # set S=Sref
+ # T = Tref - (rho - rhoRef)/alpha
+ config_eos_linear_alpha = 0.2
+ config_eos_linear_beta = 0.8
+ config_eos_linear_Tref = 10.0
+ config_eos_linear_Sref = 35.0
+ config_eos_linear_densityref = 1000.0
+
+ for k in range(0, nVertLevels):
+ activeCells = k <= maxLevelCell
+ salinity[0, activeCells, k] = S0
+ density[0, activeCells, k] = rho0 + rhoz * zMid[0, activeCells, k]
+ # T = Tref - (rho - rhoRef)/alpha
+ #temperature[0, activeCells, k] = config_eos_linear_Tref \
+ # - (density[0, activeCells, k] - config_eos_linear_densityref) / \
+ # config_eos_linear_alpha
+ for iCell in range(0, nCells):
+ if (xCell[iCell] < xCellMid):
+ temperature[0, iCell, :] = 5
+ else:
+ temperature[0, iCell, :] = 30
+ for k in range(0, nVertLevels):
+ psi1 = 1.0 - (((xCell[iCell] - 0.5*xCellMax)**4)/((0.5*xCellMax)**4))
+ psi2 = 1.0 - (((zMid[0, iCell, k] + 0.5*maxDepth)**2)/((0.5*maxDepth)**2))
+ psi = psi1*psi2
+ tracer1[0,iCell,k] = 0.5*(1 + math.tanh( r*(2*psi - 1) ) )
+ tracer2[0,:,:] = 10
+ tracer3[0,:,:] = 20
+
+ # initial velocity on edges
+ ds['normalVelocity'] = (('Time', 'nEdges', 'nVertLevels',), np.zeros([1, nEdges, nVertLevels]))
+ normalVelocity = ds['normalVelocity']
+ for iEdge in range(0, nEdges):
+ cell1 = cellsOnEdge[iEdge,0] - 1
+ cell2 = cellsOnEdge[iEdge,1] - 1
+ for k in range(0, nVertLevels):
+ zMidEdge = 0.5*(zMid[0, cell1, k] + zMid[0, cell2, k])
+ dPsi = - (2.0*zMidEdge + maxDepth) / (0.5*maxDepth)**2
+ u[0, iEdge, k] = (1.0 - ((xEdge[iEdge] - 0.5*xCellMax)**4)/((0.5*(xCellMax-xCellMin))**4)) * dPsi
+ if ((cell1<0.5*nCells and cell1>0.25*nCells-1) and (cell2<0.5*nCells and cell2>0.25*nCells-1)):
+ u[0, iEdge, k] = (1.0 - (((xEdge[iEdge]-xCellMin) - 0.5*(xCellMax+xCellMin))**4)/((0.5*(xCellMax-xCellMin))**4)) * dPsi
+ if (cell1>0.75*nCells-1 and cell2>0.75*nCells-1):
+ u[0, iEdge, k] = (1.0 - (((xEdge[iEdge]-xCellMin) - 0.5*(xCellMax+xCellMin))**4)/((0.5*(xCellMax-xCellMin))**4)) * dPsi
+ normalVelocity[0, iEdge, :] = u[0, iEdge, :] * math.cos(angleEdge[iEdge])
+ if (xEdge[iEdge]<=xCellMin or xEdge[iEdge]>=xCellMax):
+ normalVelocity[0, iEdge, :] = 0.0
+
+ # Coriolis parameter
+ ds['fCell'] = (('nCells', 'nVertLevels',), np.zeros([nCells, nVertLevels]))
+ ds['fEdge'] = (('nEdges', 'nVertLevels',), np.zeros([nEdges, nVertLevels]))
+ ds['fVertex'] = (('nVertices', 'nVertLevels',), np.zeros([nVertices, nVertLevels]))
+
+ # surface fields
+ surfaceStress[:] = 0.0
+ atmosphericPressure[:] = 0.0
+ boundaryLayerDepth[:] = 0.0
+ print(' time: %f' % ((time.time() - time1)))
+
+ comment('finalize and write file')
+ time1 = time.time()
+ ds['maxLevelCell'] = (['nCells'], maxLevelCell + 1)
+ ds['restingThickness'] = (['nCells', 'nVertLevels'], restingThickness)
+ for var in varsZ:
+ ds[var] = (['nVertLevels'], globals()[var])
+ for var in vars2D:
+ ds[var] = (['nCells'], globals()[var])
+ for var in vars3D:
+ ds[var] = (['Time', 'nCells', 'nVertLevels'], globals()[var])
+ # If you prefer not to have NaN as the fill value, you should consider
+ # using mpas_tools.io.write_netcdf() instead
+ ds.to_netcdf('initial_state.nc', format='NETCDF3_64BIT_OFFSET')
+ # write_netcdf(ds,'initial_state.nc')
+ print(' time: %f' % ((time.time() - time1)))
+ print('Total time: %f' % ((time.time() - timeStart)))
+
+
+def comment(string):
+ if verbose:
+ print('*** ' + string)
+
+
+if __name__ == '__main__':
+ # If called as a primary module, run main
+ main()
diff --git a/ocean/merry_go_round/200Cells/default/build_base_mesh.py b/ocean/merry_go_round/200Cells/default/build_base_mesh.py
new file mode 100644
index 0000000000..2d567b6a6c
--- /dev/null
+++ b/ocean/merry_go_round/200Cells/default/build_base_mesh.py
@@ -0,0 +1,26 @@
+#!/usr/bin/env python
+
+# This script was generated from setup_testcases.py as part of a config file
+
+import sys
+import os
+import shutil
+import glob
+import subprocess
+
+
+dev_null = open('/dev/null', 'w')
+
+# Run command is:
+# planar_hex --nx 1600 --ny 4 --dc 3e3 --nonperiodic_x -o planar_hex_mesh.nc
+subprocess.check_call(['planar_hex', '--nx', '200', '--ny', '4', '--dc',
+ '2.5', '--nonperiodic_x', '-o', 'planar_hex_mesh.nc'])
+
+# Run command is:
+# MpasCellCuller.x planar_hex_mesh.nc culled_mesh.nc
+subprocess.check_call(['MpasCellCuller.x', 'planar_hex_mesh.nc',
+ 'culled_mesh.nc'])
+
+# Run command is:
+# MpasMeshConverter.x culled_mesh.nc base_mesh.nc
+subprocess.check_call(['MpasMeshConverter.x', 'culled_mesh.nc', 'base_mesh.nc'])
diff --git a/ocean/merry_go_round/200Cells/default/config_base_mesh.xml b/ocean/merry_go_round/200Cells/default/config_base_mesh.xml
new file mode 100644
index 0000000000..3ac5df0612
--- /dev/null
+++ b/ocean/merry_go_round/200Cells/default/config_base_mesh.xml
@@ -0,0 +1,10 @@
+
+
+
+
+
+
+
+
+
+
diff --git a/ocean/merry_go_round/200Cells/default/config_driver.xml b/ocean/merry_go_round/200Cells/default/config_driver.xml
new file mode 100644
index 0000000000..457a7a033a
--- /dev/null
+++ b/ocean/merry_go_round/200Cells/default/config_driver.xml
@@ -0,0 +1,14 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/ocean/merry_go_round/200Cells/default/config_forward.xml b/ocean/merry_go_round/200Cells/default/config_forward.xml
new file mode 100644
index 0000000000..d587b3ae71
--- /dev/null
+++ b/ocean/merry_go_round/200Cells/default/config_forward.xml
@@ -0,0 +1,23 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 4
+
+
+
+
diff --git a/ocean/merry_go_round/200Cells/default/config_initial_state.xml b/ocean/merry_go_round/200Cells/default/config_initial_state.xml
new file mode 100644
index 0000000000..3b9176f2c1
--- /dev/null
+++ b/ocean/merry_go_round/200Cells/default/config_initial_state.xml
@@ -0,0 +1,14 @@
+
+
+
+
+
+
+
+
+
+ base_mesh.nc
+ initial_state.nc
+
+
+
diff --git a/ocean/merry_go_round/200Cells/default/config_visualize.xml b/ocean/merry_go_round/200Cells/default/config_visualize.xml
new file mode 100644
index 0000000000..573eb219e8
--- /dev/null
+++ b/ocean/merry_go_round/200Cells/default/config_visualize.xml
@@ -0,0 +1,10 @@
+
+
+
+
+
+
+
+
+
+
diff --git a/ocean/merry_go_round/200Cells/default/fromSide.py b/ocean/merry_go_round/200Cells/default/fromSide.py
new file mode 100644
index 0000000000..6e188aa6ff
--- /dev/null
+++ b/ocean/merry_go_round/200Cells/default/fromSide.py
@@ -0,0 +1,62 @@
+import matplotlib
+import matplotlib.pyplot as plt
+from netCDF4 import Dataset
+import numpy as np
+matplotlib.use('Agg')
+
+fig = plt.gcf()
+fig.set_size_inches(8.0,10.0)
+
+ncfile = Dataset('../forward/output/output.0001-01-01_00.00.00.nc', 'r')
+var1 = ncfile.variables['velocityX']
+var2 = ncfile.variables['vertVelocityTop']
+var3 = ncfile.variables['tracer1']
+var11 = var1[:, 0:200, :]
+var22 = var2[:, 0:200, :]
+var33 = var3[:, 0:200, :]
+print(var11[1,:,0])
+for iCell in range(200):
+ sumVert = 0.0
+ for k in range(100):
+ sumVert = sumVert + var11[1, iCell, k]
+ print(sumVert)
+print(var22[1,:,0])
+for k in range(100):
+ sumVert = 0.0
+ for iCell in range(200):
+ sumVert = sumVert + var22[1, iCell, k]
+ print(sumVert)
+
+plt.title('merry-go-round')
+plt.subplot(2, 2, 1)
+ax = plt.imshow(var11[1, :, :].T)
+plt.jet()
+plt.xlabel('x, # of cells')
+#plt.ylabel('depth, m')
+plt.colorbar(ax, shrink=0.5)
+plt.title('horizontal velocity')
+
+plt.subplot(2, 2, 2)
+ax = plt.imshow(var22[1, :, :].T)
+plt.jet()
+plt.xlabel('x, # of cells')
+plt.colorbar(ax, shrink=0.5)
+plt.title('vertical velocity')
+
+plt.subplot(2, 2, 3)
+ax = plt.imshow(var33[0, :, :].T)
+plt.jet()
+plt.xlabel('x, # of cells')
+plt.colorbar(ax, shrink=0.5)
+plt.title('tracer1 at t=0')
+
+plt.subplot(2, 2, 4)
+ax = plt.imshow(var33[1, :, :].T)
+plt.jet()
+plt.xlabel('x, # of cells')
+plt.colorbar(ax, shrink=0.5)
+plt.title('tracer1 at 6h')
+
+ncfile.close()
+plt.savefig('section.png')
+
diff --git a/ocean/merry_go_round/400Cells/default/add_initial_state.py b/ocean/merry_go_round/400Cells/default/add_initial_state.py
new file mode 100644
index 0000000000..29ef54d89a
--- /dev/null
+++ b/ocean/merry_go_round/400Cells/default/add_initial_state.py
@@ -0,0 +1,238 @@
+#!/usr/bin/env python
+
+'''
+This script creates an initial condition file for MPAS-Ocean.
+
+Merry-go-round Test Case: steady state test case for horizontal
+and vertical advection.
+see:
+Calandrini, S., Pieper, K. and Gunzburger, M.D., 2020. Exponential time
+differencing for the tracer equations appearing in primitive equation
+ocean models. Computer Methods in Applied Mechanics and Engineering,
+365, p.113002.
+'''
+
+import os
+import shutil
+import numpy as np
+import xarray as xr
+from mpas_tools.io import write_netcdf
+import argparse
+import math
+import time
+verbose = True
+
+
+def main():
+ timeStart = time.time()
+ parser = argparse.ArgumentParser()
+ parser.add_argument('-i', '--input_file', dest='input_file',
+ default='base_mesh.nc',
+ help='Input file, containing base mesh'
+ )
+ parser.add_argument('-o', '--output_file', dest='output_file',
+ default='initial_state.nc',
+ help='Output file, containing initial variables'
+ )
+ parser.add_argument('-L', '--nVertLevels', dest='nVertLevels',
+ default=200,
+ help='Number of vertical levels'
+ )
+ nVertLevels = parser.parse_args().nVertLevels
+ ds = xr.open_dataset(parser.parse_args().input_file)
+
+# configuration settings
+ # reference vertical grid spacing
+ maxDepth = 500
+ vertical_coordinate = 'z' # sigma or z
+
+ # comment('obtain dimensions and mesh variables')
+ nCells = ds['nCells'].size
+ nEdges = ds['nEdges'].size
+ nVertices = ds['nVertices'].size
+
+ xCell = ds['xCell']
+ xEdge = ds['xEdge']
+ xVertex = ds['xVertex']
+ yCell = ds['yCell']
+ yEdge = ds['yEdge']
+ yVertex = ds['yVertex']
+ angleEdge = ds['angleEdge']
+ cellsOnEdge = ds['cellsOnEdge']
+ edgesOnCell = ds['edgesOnCell']
+
+ # Adjust coordinates so first edge is at zero in x and y
+ xOffset = xEdge.min()
+ xCell -= xOffset
+ xEdge -= xOffset
+ xVertex -= xOffset
+ yOffset = np.min(yEdge)
+ yCell -= yOffset
+ yEdge -= yOffset
+ yVertex -= yOffset
+
+ # x values for convenience
+ xCellMax = max(xCell)
+ xCellMin = min(xCell)
+ xCellMid = 0.5 * (xCellMax + xCellMin)
+ xEdgeMax = max(xEdge)
+
+ # initialize velocity field
+ u = np.zeros([1, nEdges, nVertLevels])
+
+ comment('create and initialize variables')
+ time1 = time.time()
+
+ varsZ = [ 'refLayerThickness', 'refBottomDepth', 'refZMid', 'vertCoordMovementWeights']
+ for var in varsZ:
+ globals()[var] = np.nan * np.ones(nVertLevels)
+
+ vars2D = ['ssh', 'bottomDepth', 'surfaceStress',
+ 'atmosphericPressure', 'boundaryLayerDepth']
+ for var in vars2D:
+ globals()[var] = np.nan * np.ones(nCells)
+ maxLevelCell = np.ones(nCells, dtype=np.int32)
+
+ vars3D = [ 'layerThickness','temperature', 'salinity',
+ 'zMid', 'density', 'tracer1', 'tracer2', 'tracer3']
+ for var in vars3D:
+ globals()[var] = np.nan * np.ones([1, nCells, nVertLevels])
+ restingThickness = np.nan * np.ones([nCells, nVertLevels])
+
+ # Note that this line shouldn't be required, but if layerThickness is
+ # initialized with nans, the simulation dies. It must multiply by a nan on
+ # a land cell on an edge, and then multiply by zero.
+ layerThickness[:] = -1e34
+
+ # equally spaced layers
+ refLayerThickness[:] = maxDepth / nVertLevels
+ refBottomDepth[0] = refLayerThickness[0]
+ refZMid[0] = -0.5 * refLayerThickness[0]
+ for k in range(1, nVertLevels):
+ refBottomDepth[k] = refBottomDepth[k - 1] + refLayerThickness[k]
+ refZMid[k] = -refBottomDepth[k - 1] - 0.5 * refLayerThickness[k]
+
+ # See surface height
+ ssh[:] = 0.0
+
+ # Compute maxLevelCell and layerThickness for z-level (variation only on top)
+ if (vertical_coordinate=='z'):
+ vertCoordMovementWeights[:] = 0.0
+ vertCoordMovementWeights[0] = 1.0
+ for iCell in range(0, nCells):
+ maxLevelCell[iCell] = nVertLevels - 1
+ bottomDepth[iCell] = refBottomDepth[nVertLevels - 1]
+ layerThickness[0, iCell, :] = refLayerThickness[:]
+ layerThickness[0, iCell, 0] += ssh[iCell]
+ restingThickness[:, :] = layerThickness[0, :, :]
+
+ # Compute maxLevelCell and layerThickness for sigma
+ #elif (vertical_coordinate=='sigma'):
+ # vertCoordMovementWeights[:] = 1.0
+ # maxLevelCell[:] = nVertLevels - 1 # maxLevelCell is zero-based within python code
+ # bottomDepth[:] = bottomDepthObserved[:]
+ # for iCell in range(0, nCells):
+ # restingThickness[iCell, :] = refLayerThickness[:]*bottomDepth[iCell]/maxDepth
+ # layerThickness[0, iCell, :] = refLayerThickness[:]*(ssh[iCell] + bottomDepth[iCell])/maxDepth
+
+ # Compute zMid (same, regardless of vertical coordinate)
+ for iCell in range(0, nCells):
+ k = maxLevelCell[iCell]
+ zMid[0, iCell, k] = -bottomDepth[iCell] + \
+ 0.5 * layerThickness[0, iCell, k]
+ for k in range(maxLevelCell[iCell] - 1, -1, -1):
+ zMid[0, iCell, k] = zMid[0, iCell, k + 1] + 0.5 * \
+ (layerThickness[0, iCell, k + 1] + layerThickness[0, iCell, k])
+
+ # initialize tracers
+ rho0 = 1000.0 # kg/m^3
+ rhoz = -2.0e-4 # kg/m^3/m in z
+ S0 = 35.0
+ r = 1
+
+ # linear equation of state
+ # rho = rho0 - alpha*(T-Tref) + beta*(S-Sref)
+ # set S=Sref
+ # T = Tref - (rho - rhoRef)/alpha
+ config_eos_linear_alpha = 0.2
+ config_eos_linear_beta = 0.8
+ config_eos_linear_Tref = 10.0
+ config_eos_linear_Sref = 35.0
+ config_eos_linear_densityref = 1000.0
+
+ for k in range(0, nVertLevels):
+ activeCells = k <= maxLevelCell
+ salinity[0, activeCells, k] = S0
+ density[0, activeCells, k] = rho0 + rhoz * zMid[0, activeCells, k]
+ # T = Tref - (rho - rhoRef)/alpha
+ #temperature[0, activeCells, k] = config_eos_linear_Tref \
+ # - (density[0, activeCells, k] - config_eos_linear_densityref) / \
+ # config_eos_linear_alpha
+ for iCell in range(0, nCells):
+ if (xCell[iCell] < xCellMid):
+ temperature[0, iCell, :] = 5
+ else:
+ temperature[0, iCell, :] = 30
+ for k in range(0, nVertLevels):
+ psi1 = 1.0 - (((xCell[iCell] - 0.5*xCellMax)**4)/((0.5*xCellMax)**4))
+ psi2 = 1.0 - (((zMid[0, iCell, k] + 0.5*maxDepth)**2)/((0.5*maxDepth)**2))
+ psi = psi1*psi2
+ tracer1[0,iCell,k] = 0.5*(1 + math.tanh( r*(2*psi - 1) ) )
+ tracer2[0,:,:] = 10
+ tracer3[0,:,:] = 20
+
+ # initial velocity on edges
+ ds['normalVelocity'] = (('Time', 'nEdges', 'nVertLevels',), np.zeros([1, nEdges, nVertLevels]))
+ normalVelocity = ds['normalVelocity']
+ for iEdge in range(0, nEdges):
+ cell1 = cellsOnEdge[iEdge,0] - 1
+ cell2 = cellsOnEdge[iEdge,1] - 1
+ for k in range(0, nVertLevels):
+ zMidEdge = 0.5*(zMid[0, cell1, k] + zMid[0, cell2, k])
+ dPsi = - (2.0*zMidEdge + maxDepth) / (0.5*maxDepth)**2
+ u[0, iEdge, k] = (1.0 - ((xEdge[iEdge] - 0.5*xCellMax)**4)/((0.5*(xCellMax-xCellMin))**4)) * dPsi
+ if ((cell1<0.5*nCells and cell1>0.25*nCells-1) and (cell2<0.5*nCells and cell2>0.25*nCells-1)):
+ u[0, iEdge, k] = (1.0 - (((xEdge[iEdge]-xCellMin) - 0.5*(xCellMax+xCellMin))**4)/((0.5*(xCellMax-xCellMin))**4)) * dPsi
+ if (cell1>0.75*nCells-1 and cell2>0.75*nCells-1):
+ u[0, iEdge, k] = (1.0 - (((xEdge[iEdge]-xCellMin) - 0.5*(xCellMax+xCellMin))**4)/((0.5*(xCellMax-xCellMin))**4)) * dPsi
+ normalVelocity[0, iEdge, :] = u[0, iEdge, :] * math.cos(angleEdge[iEdge])
+ if (xEdge[iEdge]<=xCellMin or xEdge[iEdge]>=xCellMax):
+ normalVelocity[0, iEdge, :] = 0.0
+
+ # Coriolis parameter
+ ds['fCell'] = (('nCells', 'nVertLevels',), np.zeros([nCells, nVertLevels]))
+ ds['fEdge'] = (('nEdges', 'nVertLevels',), np.zeros([nEdges, nVertLevels]))
+ ds['fVertex'] = (('nVertices', 'nVertLevels',), np.zeros([nVertices, nVertLevels]))
+
+ # surface fields
+ surfaceStress[:] = 0.0
+ atmosphericPressure[:] = 0.0
+ boundaryLayerDepth[:] = 0.0
+ print(' time: %f' % ((time.time() - time1)))
+
+ comment('finalize and write file')
+ time1 = time.time()
+ ds['maxLevelCell'] = (['nCells'], maxLevelCell + 1)
+ ds['restingThickness'] = (['nCells', 'nVertLevels'], restingThickness)
+ for var in varsZ:
+ ds[var] = (['nVertLevels'], globals()[var])
+ for var in vars2D:
+ ds[var] = (['nCells'], globals()[var])
+ for var in vars3D:
+ ds[var] = (['Time', 'nCells', 'nVertLevels'], globals()[var])
+ # If you prefer not to have NaN as the fill value, you should consider
+ # using mpas_tools.io.write_netcdf() instead
+ ds.to_netcdf('initial_state.nc', format='NETCDF3_64BIT_OFFSET')
+ # write_netcdf(ds,'initial_state.nc')
+ print(' time: %f' % ((time.time() - time1)))
+ print('Total time: %f' % ((time.time() - timeStart)))
+
+
+def comment(string):
+ if verbose:
+ print('*** ' + string)
+
+
+if __name__ == '__main__':
+ # If called as a primary module, run main
+ main()
diff --git a/ocean/merry_go_round/400Cells/default/build_base_mesh.py b/ocean/merry_go_round/400Cells/default/build_base_mesh.py
new file mode 100644
index 0000000000..6709ae7255
--- /dev/null
+++ b/ocean/merry_go_round/400Cells/default/build_base_mesh.py
@@ -0,0 +1,26 @@
+#!/usr/bin/env python
+
+# This script was generated from setup_testcases.py as part of a config file
+
+import sys
+import os
+import shutil
+import glob
+import subprocess
+
+
+dev_null = open('/dev/null', 'w')
+
+# Run command is:
+# planar_hex --nx 1600 --ny 4 --dc 3e3 --nonperiodic_x -o planar_hex_mesh.nc
+subprocess.check_call(['planar_hex', '--nx', '400', '--ny', '4', '--dc',
+ '1.25', '--nonperiodic_x', '-o', 'planar_hex_mesh.nc'])
+
+# Run command is:
+# MpasCellCuller.x planar_hex_mesh.nc culled_mesh.nc
+subprocess.check_call(['MpasCellCuller.x', 'planar_hex_mesh.nc',
+ 'culled_mesh.nc'])
+
+# Run command is:
+# MpasMeshConverter.x culled_mesh.nc base_mesh.nc
+subprocess.check_call(['MpasMeshConverter.x', 'culled_mesh.nc', 'base_mesh.nc'])
diff --git a/ocean/merry_go_round/400Cells/default/config_base_mesh.xml b/ocean/merry_go_round/400Cells/default/config_base_mesh.xml
new file mode 100644
index 0000000000..3ac5df0612
--- /dev/null
+++ b/ocean/merry_go_round/400Cells/default/config_base_mesh.xml
@@ -0,0 +1,10 @@
+
+
+
+
+
+
+
+
+
+
diff --git a/ocean/merry_go_round/400Cells/default/config_driver.xml b/ocean/merry_go_round/400Cells/default/config_driver.xml
new file mode 100644
index 0000000000..457a7a033a
--- /dev/null
+++ b/ocean/merry_go_round/400Cells/default/config_driver.xml
@@ -0,0 +1,14 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/ocean/merry_go_round/400Cells/default/config_forward.xml b/ocean/merry_go_round/400Cells/default/config_forward.xml
new file mode 100644
index 0000000000..55c1bfc488
--- /dev/null
+++ b/ocean/merry_go_round/400Cells/default/config_forward.xml
@@ -0,0 +1,23 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 4
+
+
+
+
diff --git a/ocean/merry_go_round/400Cells/default/config_initial_state.xml b/ocean/merry_go_round/400Cells/default/config_initial_state.xml
new file mode 100644
index 0000000000..3b9176f2c1
--- /dev/null
+++ b/ocean/merry_go_round/400Cells/default/config_initial_state.xml
@@ -0,0 +1,14 @@
+
+
+
+
+
+
+
+
+
+ base_mesh.nc
+ initial_state.nc
+
+
+
diff --git a/ocean/merry_go_round/400Cells/default/config_visualize.xml b/ocean/merry_go_round/400Cells/default/config_visualize.xml
new file mode 100644
index 0000000000..573eb219e8
--- /dev/null
+++ b/ocean/merry_go_round/400Cells/default/config_visualize.xml
@@ -0,0 +1,10 @@
+
+
+
+
+
+
+
+
+
+
diff --git a/ocean/merry_go_round/400Cells/default/fromSide.py b/ocean/merry_go_round/400Cells/default/fromSide.py
new file mode 100644
index 0000000000..d37aa0ca40
--- /dev/null
+++ b/ocean/merry_go_round/400Cells/default/fromSide.py
@@ -0,0 +1,62 @@
+import matplotlib
+import matplotlib.pyplot as plt
+from netCDF4 import Dataset
+import numpy as np
+matplotlib.use('Agg')
+
+fig = plt.gcf()
+fig.set_size_inches(8.0,10.0)
+
+ncfile = Dataset('../forward/output/output.0001-01-01_00.00.00.nc', 'r')
+var1 = ncfile.variables['velocityX']
+var2 = ncfile.variables['vertVelocityTop']
+var3 = ncfile.variables['tracer1']
+var11 = var1[:, 0:400, :]
+var22 = var2[:, 0:400, :]
+var33 = var3[:, 0:400, :]
+print(var11[1,:,0])
+for iCell in range(400):
+ sumVert = 0.0
+ for k in range(200):
+ sumVert = sumVert + var11[1, iCell, k]
+ print(sumVert)
+print(var22[1,:,0])
+for k in range(200):
+ sumVert = 0.0
+ for iCell in range(400):
+ sumVert = sumVert + var22[1, iCell, k]
+ print(sumVert)
+
+plt.title('merry-go-round')
+plt.subplot(2, 2, 1)
+ax = plt.imshow(var11[1, :, :].T)
+plt.jet()
+plt.xlabel('x, # of cells')
+#plt.ylabel('depth, m')
+plt.colorbar(ax, shrink=0.5)
+plt.title('horizontal velocity')
+
+plt.subplot(2, 2, 2)
+ax = plt.imshow(var22[1, :, :].T)
+plt.jet()
+plt.xlabel('x, # of cells')
+plt.colorbar(ax, shrink=0.5)
+plt.title('vertical velocity')
+
+plt.subplot(2, 2, 3)
+ax = plt.imshow(var33[0, :, :].T)
+plt.jet()
+plt.xlabel('x, # of cells')
+plt.colorbar(ax, shrink=0.5)
+plt.title('tracer1 at t=0')
+
+plt.subplot(2, 2, 4)
+ax = plt.imshow(var33[1, :, :].T)
+plt.jet()
+plt.xlabel('x, # of cells')
+plt.colorbar(ax, shrink=0.5)
+plt.title('tracer1 at 6h')
+
+ncfile.close()
+plt.savefig('section.png')
+
diff --git a/ocean/merry_go_round/analysis/default/compute_rmse.py b/ocean/merry_go_round/analysis/default/compute_rmse.py
new file mode 100644
index 0000000000..3ef47035be
--- /dev/null
+++ b/ocean/merry_go_round/analysis/default/compute_rmse.py
@@ -0,0 +1,81 @@
+#!/usr/bin/env python
+'''
+This script plots results from MPAS-Ocean convergence test.
+'''
+import numpy as np
+import math
+from netCDF4 import Dataset
+import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.use('Agg')
+
+# mrp: read from file mesh after nx,ny attributes are added:
+#nx = 10 # ncfileMesh.getncattr('nx')
+#ny = 10 # ncfileMesh.getncattr('ny')
+#iz = 5
+
+dt = [12, 6, 3] # nRefinement
+order2 = [6.4, 1.6, 0.4]
+operators = ['tracer1']
+nOperators = len(operators)
+
+L2 = np.zeros([3])
+
+for k in range(nOperators):
+ ncfile12 = Dataset('../../../100Cells/default/forward/output' + '/output.0001-01-01_00.00.00.nc', 'r')
+ ncfile6 = Dataset('../../../200Cells/default/forward/output' + '/output.0001-01-01_00.00.00.nc', 'r')
+ ncfile3 = Dataset('../../../400Cells/default/forward/output' + '/output.0001-01-01_00.00.00.nc', 'r')
+ #
+ operator = operators[k]
+ areas12 = ncfile12.variables['areaCell'][:]
+ areas6 = ncfile6.variables['areaCell'][:]
+ areas3 = ncfile3.variables['areaCell'][:]
+ sol12 = ncfile12.variables[operator][1, :, 0]
+ sol6 = ncfile6.variables[operator][1, :, 0]
+ sol3 = ncfile3.variables[operator][1, :, 0]
+ ref12 = ncfile12.variables[operator][0, :, 0]
+ ref6 = ncfile6.variables[operator][0, :, 0]
+ ref3 = ncfile3.variables[operator][0, :, 0]
+ #
+ dif12 = abs(sol12 - ref12)
+ multDen12 = (ref12**2)*areas12
+ multNum12 = (dif12**2)*areas12
+ denL2_12 = np.sum(multDen12[:])/np.sum(areas12[:])
+ numL2_12 = np.sum(multNum12[:])/np.sum(areas12[:])
+ dif6 = abs(sol6 - ref6)
+ multDen6 = (ref6**2)*areas6
+ multNum6 = (dif6**2)*areas6
+ denL2_6 = np.sum(multDen6[:])/np.sum(areas6[:])
+ numL2_6 = np.sum(multNum6[:])/np.sum(areas6[:])
+ dif3 = abs(sol3 - ref3)
+ multDen3 = (ref3**2)*areas3
+ multNum3 = (dif3**2)*areas3
+ denL2_3 = np.sum(multDen3[:])/np.sum(areas3[:])
+ numL2_3 = np.sum(multNum3[:])/np.sum(areas3[:])
+ #
+ L2[0] = np.sqrt(numL2_12)/np.sqrt(denL2_12)
+ L2[1] = np.sqrt(numL2_6)/np.sqrt(denL2_6)
+ L2[2] = np.sqrt(numL2_3)/np.sqrt(denL2_3)
+ #
+ order = math.log2(L2[0]/L2[1])
+ print(order)
+ order = math.log2(L2[1]/L2[2])
+ print(order)
+ #
+ ncfile12.close()
+ ncfile6.close()
+ ncfile3.close()
+
+for k in range(len(operators)):
+ operator =operators[k]
+ plt.loglog(dt, L2[:], '-x', label='rate')
+plt.loglog(dt, order2, label='slope=-2')
+plt.title('Convergence to the exact solution')
+plt.ylabel('l_2 error norm')
+plt.legend()
+plt.grid()
+plt.xticks(dt,dt)
+plt.xlabel('time steps (in min)')
+
+plt.savefig('order.png')
+
diff --git a/ocean/merry_go_round/analysis/default/config_driver.xml b/ocean/merry_go_round/analysis/default/config_driver.xml
new file mode 100644
index 0000000000..a4c6ef926d
--- /dev/null
+++ b/ocean/merry_go_round/analysis/default/config_driver.xml
@@ -0,0 +1,8 @@
+
+
+
+
+
+
+
+
diff --git a/ocean/merry_go_round/analysis/default/config_rmse.xml b/ocean/merry_go_round/analysis/default/config_rmse.xml
new file mode 100644
index 0000000000..de2e9732e3
--- /dev/null
+++ b/ocean/merry_go_round/analysis/default/config_rmse.xml
@@ -0,0 +1,11 @@
+
+
+
+
+
+
+
+
+
+
+
diff --git a/ocean/merry_go_round/templates/template_forward_merry.xml b/ocean/merry_go_round/templates/template_forward_merry.xml
new file mode 100644
index 0000000000..d05728cbfb
--- /dev/null
+++ b/ocean/merry_go_round/templates/template_forward_merry.xml
@@ -0,0 +1,73 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ init.nc
+
+
+ init.nc
+
+
+ 0015_00:00:00
+
+
+ output
+ output/output.$Y-$M-$D_$h.$m.$s.nc
+ 0000_06:00:00
+ 01-00-00_00:00:00
+ truncate
+ 0001-01-01_00:00:00
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+