added seiche test case

config.ocean.nonhydro

@@ -0,0 +1,50 @@
+# This file is the ocean core's configuration file. It is specific to the ocean
+# core, and a specific machine. Each machine will configure this file
+# differently, but it can be used to point on version of the testing
+# infrastructure at a different version of the model.
+
+
+# The namelists section defines paths to template namelists that will be used
+# to generate specific namelists. Typically these will point to the forward and
+# init namelists in the default_inputs directory after a successful build of
+# the ocean model.
+[namelists]
+forward = /usr/projects/climate/scalandrini/repos/nonhydrostatic/namelist.ocean.forward
+init = /usr/projects/climate/scalandrini/repos/nonhydrostatic/namelist.ocean.init
+
+# The streams section defines paths to template streams files that will be used
+# to generate specific streams files. Typically these will point to the forward and
+# init streams files in the default_inputs directory after a successful build of
+# the ocean model.
+[streams]
+forward = /usr/projects/climate/scalandrini/repos/nonhydrostatic/streams.ocean.forward
+init = /usr/projects/climate/scalandrini/repos/nonhydrostatic/streams.ocean.init
+
+# The executables section defines paths to required executables. These
+# executables are provided for use by specific test cases.
+# Full paths should be provided in order to access the executables from
+# anywhere on the machine.
+[executables]
+model = /usr/projects/climate/scalandrini/repos/nonhydrostatic/ocean_model
+
+# The paths section describes paths that are used within the ocean core test
+# cases.
+[paths]
+
+# The mesh_database and the initial_condition_database are locations where
+# meshes / initial conditions might be found on a specific machine. They can be
+# the same directory, or different directory. Additionally, if they are empty
+# some test cases might download data into them, which will then be reused if
+# the test case is run again later.
+mpas_model = /usr/projects/climate/scalandrini/repos/nonhydrostatic
+mesh_database = /usr/projects/regionalclimate/COMMON_MPAS/ocean/grids/mesh_database
+initial_condition_database = /usr/projects/regionalclimate/COMMON_MPAS/ocean/grids/initial_condition_database
+bathymetry_database = /lustre/scratch4/turquoise/mpeterse/bathymetry_database
+
+# The options in this section relate to configuring and running with a compass
+# conda environment
+[conda]
+
+# Whether a link to load_compass_env.sh should be included with each test case
+# (can also be specified with the --link_load_compass flag)
+link_load_compass = False

config.ocean.training

@@ -0,0 +1,50 @@
+# This file is the ocean core's configuration file. It is specific to the ocean
+# core, and a specific machine. Each machine will configure this file
+# differently, but it can be used to point on version of the testing
+# infrastructure at a different version of the model.
+
+
+# The namelists section defines paths to template namelists that will be used
+# to generate specific namelists. Typically these will point to the forward and
+# init namelists in the default_inputs directory after a successful build of
+# the ocean model.
+[namelists]
+forward = /usr/projects/climate/scalandrini/repos/HCm_Ocean/namelist.ocean.forward
+init = /usr/projects/climate/scalandrini/repos/HCm_Ocean/namelist.ocean.init
+
+# The streams section defines paths to template streams files that will be used
+# to generate specific streams files. Typically these will point to the forward and
+# init streams files in the default_inputs directory after a successful build of
+# the ocean model.
+[streams]
+forward = /usr/projects/climate/scalandrini/repos/HCm_Ocean/streams.ocean.forward
+init = /usr/projects/climate/scalandrini/repos/HCm_Ocean/streams.ocean.init
+
+# The executables section defines paths to required executables. These
+# executables are provided for use by specific test cases.
+# Full paths should be provided in order to access the executables from
+# anywhere on the machine.
+[executables]
+model = /usr/projects/climate/scalandrini/repos/HCm_Ocean/ocean_model
+
+# The paths section describes paths that are used within the ocean core test
+# cases.
+[paths]
+
+# The mesh_database and the initial_condition_database are locations where
+# meshes / initial conditions might be found on a specific machine. They can be
+# the same directory, or different directory. Additionally, if they are empty
+# some test cases might download data into them, which will then be reused if
+# the test case is run again later.
+mpas_model = /usr/projects/climate/scalandrini/repos/HCm_Ocean
+mesh_database = /usr/projects/regionalclimate/COMMON_MPAS/ocean/grids/mesh_database
+initial_condition_database = /usr/projects/regionalclimate/COMMON_MPAS/ocean/grids/initial_condition_database
+bathymetry_database = /lustre/scratch4/turquoise/mpeterse/bathymetry_database
+
+# The options in this section relate to configuring and running with a compass
+# conda environment
+[conda]
+
+# Whether a link to load_compass_env.sh should be included with each test case
+# (can also be specified with the --link_load_compass flag)
+link_load_compass = False

ocean/seiche/stratified/default/add_initial_state.py

@@ -0,0 +1,243 @@
+#!/usr/bin/env python
+
+'''
+This script creates an initial condition file for MPAS-Ocean.
+
+Steady state ETD test case
+see:
+Demange, J., Debreu, L., Marchesiello, P., Lemarié, F., Blayo, E., Eldred, C., 2019. Stability analysis of split-explicit free surface ocean models: Implication of the depth-independent barotropic mode approximation. 
+Journal of Computational Physics 398, 108875. https://doi.org/10.1016/j.jcp.2019.108875
+'''
+
+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=40,
+                        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 = 16
+    vertical_coordinate = 'uniform'  # sigma or z or uniform
+
+    # 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)
+    print(xCellMax)
+    print(xCellMin)   
+
+    # 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', #'bottomDepthObserved',
+        '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']
+    #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]
+
+    # Marsaleix et al 2008 page 81
+    # Gaussian function in depth for deep sea ridge 
+    #bottomDepthObserved[:] = maxDepth - 1000.0 * np.exp(-((xCell[:] - xMid) / 15e3)**2)
+    # SSH varies from -1.0 to 1.0m across the domain
+    #ssh[:] = (xCell[:] - xMid)/ xMid
+    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, 0:maxLevelCell[iCell] ] = refLayerThickness[0:maxLevelCell[iCell]]
+            layerThickness[0, iCell, :] = refLayerThickness[:]
+            layerThickness[0, iCell, 0] += ssh[iCell]
+        restingThickness[:, :] = layerThickness[0, :, :]
+    #Compute maxLevelCell and layerThickness for uniform
+    elif (vertical_coordinate=='uniform'):  
+        vertCoordMovementWeights[:] = 1.0
+        vertCoordMovementWeights[0] = 1.0
+        for iCell in range(0, nCells): 
+            maxLevelCell[iCell] = nVertLevels - 1
+            bottomDepth[iCell] = refBottomDepth[nVertLevels - 1]
+            layerThickness[0, iCell, :] = refLayerThickness[:] + ssh[iCell]/nVertLevels 
+        #restingThickness[:, :] = layerThickness[0, :, :]
+        restingThickness[:, :] = refLayerThickness[:]
+
+    # 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) #TODO: not sure it is correct
+    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]
+        density[0, activeCells, k] = rho0 + (10.0/2)*(1.0 - np.tanh((2/1)*np.arctanh(0.99)*(zMid[0, :, k] + 0.5*maxDepth - 0.1*np.cos((np.pi/10)*xCell[:]))))
+        # 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): #TODO: what values should we give to the temperature?    
+    #    if (xCell[iCell] < xCellMid): 
+    #        temperature[0, iCell, :] = 20
+    #    else: 
+    #        temperature[0, iCell, :] = 20  
+    #    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):  
+        normalVelocity[0, iEdge, :] = u[0, iEdge, :] * math.cos(angleEdge[iEdge])
+
+    # 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()

ocean/seiche/stratified/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', '64', '--ny', '4', '--dc',
+                       '0.15625', '--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'])

ocean/seiche/stratified/default/config_base_mesh.xml

@@ -0,0 +1,10 @@
+<?xml version="1.0"?>
+<config case="base_mesh">
+	<copy_file source_path="script_test_dir" source="build_base_mesh.py" dest="build_base_mesh.py"/>
+
+	<run_script name="run.py">
+                <step executable="python">
+                        <argument flag="build_base_mesh.py"></argument>  
+                </step>
+	</run_script>
+</config>

ocean/seiche/stratified/default/config_driver.xml

@@ -0,0 +1,14 @@
+<driver_script name="run.py">
+	<case name="base_mesh">
+		<step executable="./run.py" quiet="true" pre_message=" * Creating 2D global base mesh with jigsaw..." post_message="     complete!  Created file:  base_mesh/base_mesh.nc"/>
+	</case>
+	<case name="initial_state">
+		<step executable="./run.py" quiet="true" pre_message=" * Initializing ocean state with bathymetry and tracers..." post_message="     complete!  Created file:  initial_state/initial_state.nc"/>
+	</case>
+	<case name="forward">
+		<step executable="./run.py"/>
+	</case>
+        <case name="visualize">
+		<step executable="./run.py" quiet="true" pre_message=" * Visualize velocity and tracer fields..." post_message="     complete!  Created file:  visualize/section.png"/>
+	</case>
+</driver_script>