REMORA

REMORA is currently under development as a next-generation version of the Regional Ocean Modeling System (ROMS).

REMORA is built on AMReX, an adaptive mesh refinement software framework, which provides the underlying software infrastructure for block structured AMR operations. The full AMReX documentation can be found here, and the tutorials can be found here.

REMORA is designed to run on machines from laptops to multicore CPU and hybrid CPU/GPU systems.

This documentation is currently under development, there are detailed resources available on the ROMS Documentation Portal for the Regional Ocean Modeling System.

In addition to this documentation, there is API documentation for REMORA generated by Doxygen.

Code of Conduct

Our Pledge

In the interest of fostering an open and welcoming environment, we as contributors and maintainers pledge to making participation in our project and our community a harassment-free experience for everyone, regardless of age, body size, disability, ethnicity, sex characteristics, gender identity and expression, level of experience, education, socio-economic status, nationality, personal appearance, race, religion, or sexual identity and orientation.

Our Standards

Examples of behavior that contributes to creating a positive environment include:

Using welcoming and inclusive language
Being respectful of differing viewpoints and experiences
Gracefully accepting constructive criticism
Focusing on what is best for the community
Showing empathy towards other community members

Examples of unacceptable behavior by participants include:

The use of sexualized language or imagery and unwelcome sexual attention or advances
Trolling, insulting/derogatory comments, and personal or political attacks
Public or private harassment
Publishing others’ private information, such as a physical or electronic address, without explicit permission
Other conduct which could reasonably be considered inappropriate in a professional setting

Our Responsibilities

Project maintainers are responsible for clarifying the standards of acceptable behavior and are expected to take appropriate and fair corrective action in response to any instances of unacceptable behavior.

Project maintainers have the right and responsibility to remove, edit, or reject comments, commits, code, wiki edits, issues, and other contributions that are not aligned to this Code of Conduct, or to ban temporarily or permanently any contributor for other behaviors that they deem inappropriate, threatening, offensive, or harmful.

Scope

This Code of Conduct applies both within project spaces and in public spaces when an individual is representing the project or its community. Examples of representing a project or community include using an official project e-mail address, posting via an official social media account, or acting as an appointed representative at an online or offline event. Representation of a project may be further defined and clarified by project maintainers.

Enforcement

Instances of abusive, harassing, or otherwise unacceptable behavior may be reported by contacting the project team. All complaints will be reviewed and investigated and will result in a response that is deemed necessary and appropriate to the circumstances. The project team is obligated to maintain confidentiality with regard to the reporter of an incident. Further details of specific enforcement policies may be posted separately.

Project maintainers who do not follow or enforce the Code of Conduct in good faith may face temporary or permanent repercussions as determined by other members of the project’s leadership.

Attribution

This Code of Conduct is adapted from the Contributor Covenant, version 1.4, available at https://www.contributor-covenant.org/version/1/4/code-of-conduct.html

For answers to common questions about this code of conduct, see https://www.contributor-covenant.org/faq

Getting Started

Downloading the code

First, make sure that git is installed on your machine.

Then download the REMORA repository by typing:

git clone https://github.com/seahorce-scidac/REMORA.git

Or, to automatically include the AMReX submodule when downloading REMORA, type:

git clone --recursive https://github.com/seahorce-scidac/REMORA.git

Git Submodules

Ideally, when using the submodule to build it is properly updated to match what is in the repository. Depending on Git version, different commands and options to ensure these match are available. An example workflow is to run git pull to get the latest commit on your current branch, and then run git submodule update to explicitly update the submodule. This should work for all versions of git which support submodules.

The following example demonstrates a shorter form that combines both commands and requires Git 2.14 or newer:

# Replaces your git pull to use both the updated code and the updated submodule
git pull --recurse-submodules

The following example demonstrates setting defaults in the config file for the repository and requires Git 2.15 or newer:

# Set this once
git config submodule.recurse true
# When configured as true, this will use both the updated code and the updated submodule
git pull
# This option will override any configuration and not update the submodule
git pull --no-recurse-submodules

These example also apply to git checkout. For more details see Git Tools Submodules: https://git-scm.com/book/en/v2/Git-Tools-Submodules

Building

REMORA can be built using either GNU Make or CMake.

GNU Make

The GNU Make system is best for use on large computing facility machines and production runs. With the GNU Make implementation, the build system will inspect the machine and use known compiler optimizations explicit to that machine if possible. These explicit settings are kept up-to-date by the AMReX project.

Using the GNU Make build system involves first setting environment variables for the directories of the dependencies of REMORA which is the repository of AMReX. AMReX is provided as a git submodule in REMORA and can be populated by using git submodule init; git submodule update in the REMORA repo, or before cloning by using git clone --recursive <remora_repo>. Although submodules of these projects are provided, they can be placed externally as long as the <REPO_HOME> environment variables for each dependency is set correctly. An example of setting the <REPO_HOME> environment variables in the user’s .bashrc is shown below:

export REMORA_HOME=${HOME}/REMORA
export AMREX_HOME=${REMORA_HOME}/Submodules/AMReX

The GNU Make system is set up to use the path to AMReX submodule by default, so it is not necessary to set these paths explicitly, unless it is desired to do so. It is also possible to use an external version of AMReX, downloaded by running

git clone https://github.com/amrex-codes/amrex.git

in which case the AMREX_HOME environment variable must point to the location where AMReX has been downloaded, which will take precedence over the default path to the submodule. If using bash shell,

export AMREX_HOME=/path/to/external/amrex

or if using tcsh,

setenv AMREX_HOME /path/to/external/amrex

cd to the desired build directory, e.g. REMORA/Exec/Upwelling/

Edit the GNUmakefile; options include

Option name	Description	Possible values	Default value
COMP	Compiler (gnu or intel)	gnu / intel	None
USE_MPI	Whether to enable MPI	TRUE / FALSE	FALSE
USE_OMP	Whether to enable OpenMP	TRUE / FALSE	FALSE
USE_CUDA	Whether to enable CUDA	TRUE / FALSE	FALSE
DEBUG	Whether to use DEBUG mode	TRUE / FALSE	FALSE
USE_NETCDF	Whether to compile with NETCDF	TRUE / FALSE	FALSE
USE_PARTICLES	Whether to compile with particle functionality enabled	TRUE / FALSE	FALSE
PROFILE	Include profiling info	TRUE / FALSE	FALSE
TINY_PROFILE	Include tiny profiling info	TRUE / FALSE	FALSE
COMM_PROFILE	Include comm profiling info	TRUE / FALSE	FALSE
TRACE_PROFILE	Include trace profiling info	TRUE / FALSE	FALSE

Note

Do not set both USE_OMP and USE_CUDA to true.

Information on using other compilers can be found in the AMReX documentation at https://amrex-codes.github.io/amrex/docs_html/BuildingAMReX.html .

Make the executable by typing
```
make
```
The name of the resulting executable (generated by the GNUmake system) encodes several of the build characteristics, including dimensionality of the problem, compiler name, and whether MPI and/or OpenMP were linked with the executable. Thus, several different build configurations may coexist simultaneously in a problem folder. For example, the default build in REMORA/Exec/Upwelling will look like REMORA3d.gnu.MPI.ex, indicating that this is a 3-d version of the code, made with COMP=gnu, and USE_MPI=TRUE.

Job info

The build information can be accessed by typing

./REMORA*ex --describe

in the directory where the executable has been built.

CMake

CMake is often preferred by developers of REMORA; CMake allows for building as well as easy testing and verification of REMORA through the use of CTest which is included in CMake. CTest functionality requires additional options, described in Testing and Verification.

Using CMake involves an additional configure step before using the make command. It is also expected that the user has cloned the REMORA repo with the --recursive option or performed git submodule init; git submodule update in the REMORA repo to populate its submodules.

To build with CMake, a user typically creates a build directory in the project directory and in that directory the cmake <options> .. command is used to configure the project before building it. REMORA provides an example build directory called Build with example scripts for performing the CMake configure. Once the CMake configure step is done, then the make command will build the executable.

An example CMake configure command to build REMORA with MPI is listed below:

cmake -DCMAKE_BUILD_TYPE:STRING=Release \
      -DREMORA_ENABLE_MPI:BOOL=ON \
      -DCMAKE_CXX_COMPILER:STRING=mpicxx \
      -DCMAKE_C_COMPILER:STRING=mpicc \
      -DCMAKE_Fortran_COMPILER:STRING=mpifort \
      .. && make

An example CMake configure command to build REMORA with MPI and particles is listed below:

cmake -DCMAKE_BUILD_TYPE:STRING=Release \
      -DREMORA_ENABLE_MPI:BOOL=ON \
      -DCMAKE_CXX_COMPILER:STRING=mpicxx \
      -DCMAKE_C_COMPILER:STRING=mpicc \
      -DCMAKE_Fortran_COMPILER:STRING=mpifort \
      -DREMORA_ENABLE_PARTICLES:BOOL=ON \
      .. && make

Note that CMake is able to generate makefiles for the Ninja build system as well which will allow for faster building of the executable(s).

Perlmutter (NERSC)

Recall the GNU Make system is best for use on large computing facility machines and production runs. With the GNU Make implementation, the build system will inspect the machine and use known compiler optimizations explicit to that machine if possible. These explicit settings are kept up-to-date by the AMReX project.

For Perlmutter at NERSC, look at the general instructions for building REMORA using GNU Make, and then you can initialize your environment by sourcing or running the saul-env.sh script in the Build directory. GNU Make may complain that it cannot find NetCDF. This is fine.

Then build REMORA as, for example (specify your own path to the AMReX submodule in REMORA/Submodules/AMReX):

make -j 4 COMP=gnu USE_MPI=TRUE USE_OMP=FALSE USE_CUDA=TRUE AMREX_HOME=/global/u2/d/dwillcox/dev-remora.REMORA/Submodules/AMReX USE_SUNDIALS=FALSE

Finally, you can prepare your SLURM job script, using the following as a guide:

#!/bin/bash

## specify your allocation (with the _g) and that you want GPU nodes
#SBATCH -A mXXXX_g
#SBATCH -C gpu

## the job will be named "REMORA" in the queue and will save stdout to remora_[job ID].out
#SBATCH -J REMORA
#SBATCH -o remora_%j.out

## set the max walltime
#SBATCH -t 10

## specify the number of nodes you want
#SBATCH -N 2

## we use the same number of MPI ranks per node as GPUs per node
#SBATCH --ntasks-per-node=4
#SBATCH --gpus-per-node=4
#SBATCH --gpu-bind=none

# pin to closest NIC to GPU
export MPICH_OFI_NIC_POLICY=GPU

# use GPU-aware MPI
#GPU_AWARE_MPI=""
GPU_AWARE_MPI="amrex.use_gpu_aware_mpi=1"

# the -n argument is (--ntasks-per-node) * (-N) = (number of MPI ranks per node) * (number of nodes)
# set ordering of CUDA visible devices inverse to local task IDs for optimal GPU-aware MPI
srun -n 8 --cpus-per-task=32 --cpu-bind=cores bash -c "
  export CUDA_VISIBLE_DEVICES=\$((3-SLURM_LOCALID));
  ./REMORA3d.gnu.MPI.CUDA.ex inputs ${GPU_AWARE_MPI}" \
> test.out

To submit your job script, do sbatch [your job script] and you can check its status by doing squeue -u [your username].

Running

The input file specified on the command line is a free-format text file, one entry per row, that specifies input data processed by the AMReX ParmParse module.

This file needs to be specified along with the executable as an argv option, for example:

mpirun -np 64 ./REMORA3d.xxx.yyy.ex inputs

Also, any entry that can be specified in the inputs file can also be specified on the command line; values specified on the command line override values in the inputs file, e.g.:

mpirun -np 64 ./REMORA3d.gnu.DEBUG.MPI.ex inputs amr.restart=chk0030 remora.use_gravity=true

See Inputs for details on run-time options that can be specified.

Testing and Verification

Testing and verification of REMORA can be performed using CTest, which is included in the CMake build system. If one builds REMORA with CMake, the testing suite, and the verification suite, can be enabled during the CMake configure step. A nightly test is reflected on the dashboard here .

An example cmake configure command performed in the Build directory in REMORA is shown below with options relevant to the testing suite:

cmake -DCMAKE_INSTALL_PREFIX:PATH=./install \
      -DCMAKE_BUILD_TYPE:STRING=Release \
      -DREMORA_ENABLE_MPI:BOOL=ON \
      -DCMAKE_CXX_COMPILER:STRING=mpicxx \
      -DCMAKE_C_COMPILER:STRING=mpicc \
      -DCMAKE_Fortran_COMPILER:STRING=mpifort \
      -DREMORA_ENABLE_FCOMPARE:BOOL=ON \
      -DREMORA_ENABLE_TESTS:BOOL=ON \
      -DREMORA_USE_CPP:BOOL=ON \
      ..

While performing a cmake -LAH .. command will give descriptions of every option for the CMake project. Descriptions of particular options regarding the testing suite are listed below:

REMORA_ENABLE_FCOMPARE – builds the fcompare utility from AMReX as well as the executable(s), to allow for testing differences between plot files

REMORA_ENABLE_TESTS – enables the base level regression test suite that will check whether each test will run its executable to completion successfully

Building the Tests

Once the user has performed the CMake configure step, the make command will build every executable required for each test. In this step, it is highly beneficial for the user to use the -j option for make to build source files in parallel.

Running the Tests

Once the test executables are built, CTest also creates working directories for each test within the Build directory where plot files will be output, etc. This directory is analogous to the source location of the tests in Tests/test_files.

To run the test suite, run ctest in the Build directory. CTest will run the tests and report their exit status. Useful options for CTest are -VV which runs in a verbose mode where the output of each test can be seen. -R where a regex string can be used to run specific sets of tests. -j where CTest will bin pack and run tests in parallel based on how many processes each test is specified to use and fit them into the amount of cores available on the machine. -L where the subset of tests containing a particular label will be run. (We note that using -L 'regression' will run all the tests that do not use SUNDIALS.) Output for the last set of tests run is available in the Build directory in Testing/Temporary/LastTest.log.

Adding Tests

Developers are encouraged to add tests to REMORA and in this section we describe how the tests are organized in the CTest framework. The locations (relative to the REMORA code base) of the tests are in Tests. To add a test, first create a problem directory with a name in Exec/<prob_name>. This problem directory is meant for a production run where the simulation is run until convergence or a solution is developed. This problem setup could comprise of a more complex physics than the corresponding tests for regression at Tests/test_files/<test_name>. Prepare toned down versions of the input file(s) for each combination of physics that a regression test is desired. For example, Upwelling problem with input file Exec/Upwelling/inputs solves the double gyre problem. The corresponding regression test is driven by the input files Tests/test_files/Upwelling/Upwelling.i.

Any file in the test directory will be copied during CMake configure to the test’s working directory. The input files meant for regression test run only until a few time steps. The reference solution that the regression test will refer to should be placed in Tests/REMORA-Gold-Files/<test_name>. Next, edit the Exec/CMakeLists.txt and Tests/CTestList.cmake files, add the problem and the corresponding tests to the list. Note that there are different categories of tests and if your test falls outside of these categories, a new function to add the test will need to be created. After these steps, your test will be automatically added to the test suite database when doing the CMake configure with the testing suite enabled.

Inputs

The REMORA executable reads run-time information from an “inputs” file which you put on the command line. This section describes the inputs which can be specified either in the inputs file or on the command line. If a value is specified on the command line, that value will override a value specified in the inputs file.

Problem Geometry

List of Parameters

Parameter	Definition	Acceptable Values	Default
geometry.prob_lo	physical location of low corner of the domain	Real	must be set
geometry.prob_hi	physical location of high corner of the domain	Real	must be set
geometry.is_periodic	is the domain periodic in this direction	0 if false, 1 if true	0 0 0

Examples of Usage

geometry.prob_lo = 0 0 0 defines the low corner of the domain at (0,0,0) in physical space.
geometry.prob_hi = 1.e8 2.e8 2.e8 defines the high corner of the domain at (1.e8,2.e8,2.e8) in physical space.
geometry.is_periodic = 0 1 0 says the domain is periodic in the y-direction only.

Domain Boundary Conditions

List of Parameters

Parameter	Definition	Default
xlo.type	boundary type of xlo face	must be set if not periodic
xhi.type	boundary type of xhi face	must be set if not periodic
ylo.type	boundary type of ylo face	must be set if not periodic
yhi.type	boundary type of yhi face	must be set if not periodic
zlo.type	boundary type of zlo face	must be set if not periodic
zhi.type	boundary type of zhi face	must be set if not periodic

Resolution

List of Parameters

Parameter	Definition	Acceptable Values	Default
amr.n_cell	number of cells in each direction at the coarsest level	Integer > 0	must be set
amr.max_level	number of levels of refinement above the coarsest level	Integer >= 0	must be set
amr.ref_ratio	ratio of coarse to fine grid spacing between subsequent levels	2 / 3 / 4 (one per level)	2 for all levels
amr.ref_ratio_vect	ratio of coarse to fine grid spacing between subsequent levels	3 integers (one per dir) 2 / 3 / 4	2 for all directions
amr.regrid_int	how often to regrid	Integer > 0	must be set
amr.regrid_on_restart	should we regrid immediately after restarting	0 or 1	0

Note: if amr.max_level = 0 then you do not need to set amr.ref_ratio or amr.regrid_int.

Examples of Usage

amr.n_cell = 32 64 64

would define the domain to have 32 cells in the x-direction, 64 cells in the y-direction, and 64 cells in the z-direction at the coarsest level.
amr.max_level = 2

would allow a maximum of 2 refined levels in addition to the coarse level. Note that these additional levels will only be created only if the tagging criteria are such that cells are flagged as needing refinement. The number of refined levels in a calculation must be \(\leq\) amr.max_level, but can change in time and need not always be equal to amr.max_level.
amr.ref_ratio = 2 3

would set factor 2 refinement between levels 0 and 1, and factor 3 refinement between levels 1 and 2. Note that you must have at least amr.max_level values of amr.ref_ratio (Additional values may appear in that line and they will be ignored).
amr.ref_ratio_vect = 2 4 3

would set factor {2 in x-dir, 4 in y-dir, 3 in z-dir} refinement between all adjacent levels. Note that you must specify 3 values, one for each coordinate direction.
amr.regrid_int = 2 2

tells the code to regrid every 2 steps. Thus in this example, new level-1 grids will be created every 2 level-0 time steps, and new level-2 grids will be created every 2 level-1 time steps.

Regridding

Overview

The user defines how to tag individual cells at a given level for refinement. This list of tagged cells is sent to a grid generation routine, which uses the Berger–Rigoutsos algorithm to create rectangular grids that contain the tagged cells.

See Mesh Refinement for more details on how to specify regions for refinement.

List of Parameters

Parameter	Definition	Acceptable Values	Default
amr.regrid_file	name of file from which to read the grids	text	no file
amr.grid_eff	grid efficiency at coarse level at which grids are created	Real > 0, < 1	0.7
amr.n_error_buf	radius of additional tagging around already tagged cells	Integer >= 0	1
amr.max_grid_size	maximum size of a grid in any direction	Integer > 0	32
amr.max_grid_size	maximum size	Integer	32
amr.blocking_factor	grid size must be a multiple of this	Integer > 0	2
amr.refine_grid_layout	refine grids more if # of processors \(>\) # of grids	0 if false, 1 if true	1

Notes

amr.n_error_buf, amr.max_grid_size and amr.blocking_factor can be read in as a single value which is assigned to every level, or as multiple values, one for each level
amr.max_grid_size at every level must be even
amr.blocking_factor at every level must be a power of 2
the domain size amr.n_cell must be a multiple of amr.blocking_factor at level 0
amr.max_grid_size must be a multiple of amr.blocking_factor at every level

Examples of Usage

amr.regrid_file = fixed_grids

In this case the list of grids at each fine level are contained in the file fixed_grids, which will be read during the gridding procedure. These grids must not violate the amr.max_grid_size criterion. The rest of the gridding procedure described below will not occur if amr.regrid_file is set.
amr.grid_eff = 0.9

During the grid creation process, at least 90% of the cells in each grid at the level at which the grid creation occurs must be tagged cells. Note that this is applied at the coarsened level at which the grids are actually made, and before amr.max_grid_size is imposed.
amr.max_grid_size = 64

The final grids will be no longer than 64 cells on a side at every level.
amr.max_grid_size = 64 32 16

The final grids will be no longer than 64 cells on a side at level 0, 32 cells on a side at level 1, and 16 cells on a side at level 2.
amr.blocking_factor = 32

The dimensions of all the final grids will be multiples of 32 at all levels.
amr.blocking_factor = 32 16 8

The dimensions of all the final grids will be multiples of 32 at level 0, multiples of 16 at level 1, and multiples of 8 at level 2.

Gridding and Load Balancing

The REMORA gridding and load balancing strategy is based on that in AMReX. See the Gridding section of the AMReX documentation for details.

Simulation Time

List of Parameters

Parameter	Definition	Acceptable Values	Default
max_step	maximum number of level 0 time steps	Integer >= 0	-1
stop_time	final simulation time	Real >= 0	-1.0
remora.start_time	initial simulation time	Real >= 0	0.0

Notes

To control the number of time steps, you can limit by the maximum number of level-0 time steps (max_step), or the final simulation time (stop_time), or both. The code will stop at whichever criterion comes first. Note that if the code reaches stop_time then the final time step will be shortened so as to end exactly at stop_time, not pass it.

Examples of Usage

max_step = 1000
stop_time = 1.0

will end the calculation when either the simulation time reaches 1.0 or the number of level-0 steps taken equals 1000, whichever comes first.

Time Step

List of Parameters for Single-Rate

Parameter	Definition	Acceptable Values	Default
remora.cfl	CFL number for hydro	Real > 0 and <= 1	0.8
remora.fixed_dt	set level 0 dt as this value regardless of cfl or other settings	Real > 0	unused if not set
remora.fixed_fast_dt	set fast dt as this value	real > 0	inferred from fixed_dt and fixed_ndfast_ratio if not set
remora.fixed_ndfast_ratio	set fast dt as slow dt / this ratio	int	inferred from fixed_dt and fixed_fast_dt if not set
remora.change_max	factor by which dt can grow in subsequent steps	Real >= 1	1.1

Examples of Usage

remora.cfl = 0.9

defines the timestep as dt = cfl * dx / (u+c). Only relevant if fixed_dt not set
remora.change_max = 1.1

allows the time step to increase by no more than 10% in this case. Note that the time step can shrink by any factor; this only controls the extent to which it can grow.
remora.fixed_dt = 1.e-4

sets the level-0 time step to be 1.e-4 for the entire simulation, ignoring the other timestep controls. Note that if remora.init_shrink \(\neq 1\) then the first time step will in fact be remora.init_shrink * remora.fixed_dt.

Restart Capability

See Checkpoint / Restart for how to control the checkpoint/restart capability.

PlotFiles

See Plotfiles for how to control the types and frequency of plotfile generation.

Screen Output

List of Parameters

Parameter	Definition	Acceptable Values	Default
amr.v	verbosity of Amr.cpp	0 or 1	0
remora.v	verbosity of REMORA functions	0: none 1: integrated/max quantities 2: print boxes	0
remora.sum_interval	how often (in level-0 time steps) to compute integral quantities	Integer	-1

Examples of Usage

remora.sum_interval = 2

if remora.sum_interval \(> 0\) then the code computes and prints certain integral quantities, such as total mass, momentum and energy in the domain every remora.sum_interval level-0 steps. In this example the code will print these quantities every two coarse time steps. The print statements have the form

TIME= 1.91717746 MASS= 1.792410279e+34

for example. If this line is commented out or if remora.v \(<= 0\) then it will not compute and print these quanitities.

Included terms

List of Parameters

Parameter	Definition	Acceptable Values	Default
remora.use_coriolis	Include Coriolis terms.	true / false	false
remora.flat_bathymetry	Use flat bathymetry.	true / false	true
remora.use_prestep	Do prestep terms. Only for debugging purposes.	true / false	true
remora.use_uv3dmix	Include harmonic viscosity. Only for debugging purposes.	true / false	true
remora.use_barotropic	Include 2d barotropic step. Only for debugging purposes.	true / false	true

Physics Parameters

List of Parameters

Parameter	Definition	Acceptable Values	Default
remora.ggrav	Gravitational field strength [kg m/s^2]	Real number	9.81
remora.R0	Background density [kg/m^3] used in Linear Equation of State. May be used in setup of some problems.	Real number	1028
remora.S0	Background salinity (nondimensional) used in Linear Equation of State State. May be used in setup of some problems.	Real number	35
remora.T0	Background temperature (Celsius) used in Linear Equation of State State. May be used in setup of some problems.	Real number	5
remora.Tcoef	Linear EOS parameter (1/Celsius)	Real number	1.7e-4
remora.Scoef	Linear EOS parameter (nondimensional)	Real number	0.0
remora.rho0	Mean density (kg/m^3) used when Boussinesq approx is inferred	Real number	1025
remora.coriolis_type	Type of Coriolis forcing. `beta_plane` uses a linear approximation. `custom` is calculated from a function in `prob.cpp`, and `real` is read from the netcdf grid file	`beta_plane` / `custom` / `real`	`beta_plane`
remora.coriolis_f0	f-plane constant for Coriolis param \(f = f_0 + \beta y\) when using beta plane Coriolis type	Real number	0.0
remora.coriolis_beta	beta-plane constant for Coriolis param \(f = f_0 + \beta y\) when using beta plane Coriolis type	Real number	0.0

Numerical Algorithms

List of Parameters

Parameter	Definition	Acceptable Values	Default
remora.horizontal_advection_scheme	Scheme for horizontal advection	upstream3, centered4	upstream3

Vertical Stretch prameters

List of Parameters

Parameter	Definition	Acceptable Values	Default
remora.theta_s	Stretching parameter for surface refinement of vertical S-grid	\(0 \leq \theta_S \leq 10\)	3.0
remora.theta_b	Stretching parameter for bottom refinement of vertical S-grid	\(0 \leq \theta_B \leq 4\)	0.0
remora.tcline	Surface/bottom layer width (m) in vertical S-grid	Positive number	150

These parameters are used to calculate the vertical S-grid stretch/transform functions detailed in Vertical S-Coordinate.

Mesh Refinement

REMORA allows both static and dynamic mesh refinement, as well as the choice of one-way or two-way coupling.

Note that any tagged region will be covered by one or more boxes. The user may specify the refinement criteria and/or region to be covered, but not the decomposition of the region into individual grids.

See the Gridding section of the AMReX documentation for details of how individual grids are created.

Coupling Types

REMORA supports one-way and two-way coupling between levels; this is a run-time input

remora.coupling_type = "OneWay" or "TwoWay"

By one-way coupling, we mean that between each pair of refinement levels, the coarse level communicates data to the fine level to serve as boundary conditions for the time advance of the fine solution. For cell-centered quantities, and face-baced normal momenta on the coarse-fine interface, the coarse data is conservatively interpolated to the fine level.

The interpolated data is utilized to specify ghost cell data (outside of the valid fine region).

By two-way coupling, we mean that in additional to interpolating data from the coarser level to supply boundary conditions for the fine regions, the fine level also communicates data back to the coarse level in two ways:

The fine cell-centered data are conservatively averaged onto the coarse mesh covered by fine mesh.
The fine momenta are conservatively averaged onto the coarse faces covered by fine mesh.
A “reflux” operation is performed for all cell-centered data; this updates values on the coarser level outside of regions covered by the finer level.

Advected quantities which are advanced in conservation form will lose conservation with one-way coupling. Two-way coupling ensures conservation of the advective contribution to all scalar updates but does not account for loss of conservation due to diffusive or source terms.

Checkpoint / Restart

REMORA has a standard sort of checkpointing and restarting capability and uses the native AMReX format for reading and writing checkpoints. In the inputs file, the following options control the generation of checkpoint files (which are really directories):

Writing the Checkpoint “Files”

List of Parameters

Parameter	Definition	Acceptable Values	Default
remora.check_file	prefix for restart files	String	“chk”
remora.check_int	how often (by level-0 time steps) to write restart files	Integer \(> 0\)	-1

Restarting

Parameter	Definition	Acceptable Values	Default
remora.restart	name of the file (directory) from which to restart files	String	not used if not set

Examples of Usage

remora.check_file = chk_run
remora.check_int = 10

means that restart files (really directories) starting with the prefix “chk_run” will be generated every 10 level-0 time steps. The directory names will be chk_run00000, chk_run00010, chk_run00020, etc.

To restart from chk_run00061,for example, then set

amr.restart = chk_run00061

Plotfiles

Controlling PlotFile Generation

“Plotfiles” can be written very efficiently in parallel in a native AMReX format or in HDF5. They can also be written in NetCDF.

The following options in the inputs file control the generation of plotfiles.

List of Parameters

Parameter	Definition	Acceptable Values	Default
remora.plotfile_type	AMReX or NETCDF	“amrex” or “netcdf / “NetCDF”	“amrex”
remora.write_history_file	do we write netcdf files at each timestep or one file for all timesteps?	false or true	false
remora.plot_file	prefix for plotfiles	String	“plt_”
remora.plot_int	how often (by level-0 time steps) to write plot files	Integer \(> 0\)	-1
remora.plot_vars	name of variables to include in plotfiles	list of names	None

Notes

The NeTCDF option is only available if REMORA has been built with USE_NETCDF enabled.
The write_history_file option is only available if plotfile_type = netcdf
If plotfile_type = netcdf and write_history_file = false, the frequency will be determined by plot_int
Velocity components are defined on faces within the REMORA code, but are averaged onto cell centers when written in amrex/native plotfiles. They are not averaged when writing NetCDF files.
File prefixes can include directories.

Examples of Usage

remora.plotfile_type = amrex
remora.plot_file = out/plt_run
remora.plot_int = 10

means that native plot files (actually directories) starting with the prefix “plt_run” will be generated every 10 level-0 time steps in the directory out. If using amrex format, that directory names will be plt_run00000, plt_run00010, plt_run00020, etc. If using NetCDF format, the names will have “.nc” appended.

AMReX plotfiles will contain data at all of the refinement levels. NetCDF files will not be output if there is more than one level.

Visualization

By default, REMORA currently generates plotfile in the native AMReX format, but they can also be written in HDF5 or NetCDF.

There are several visualization tools that can be used for AMReX plotfiles, specifically ParaView, VisIt and yt.

ParaView

The open source visualization package ParaView v5.10 and later can be used to view REMORA plotfiles with and without terrain. You can download the paraview executable at https://www.paraview.org/.

To open a plotfile

Run ParaView v5.10, then select “File” \(\rightarrow\) “Open”.
Navigate to your run directory, and select either a single plotfile or a set of plotfiles. Open multiple plotfile at once by selecting plt.. Paraview will load the plotfiles as a time series. ParaView will ask you about the file type – choose “AMReX/BoxLib Grid Reader”.
If you have run the REMORA executable with terrain, then the mapped grid information will be stored as nodal data. Choose the “point data” called “nu”, then click on “Warp by Vector” which can be found via Filters–>Alphabetical. This wil then plot data onto the mapped grid locations.
Under the “Cell Arrays” field, select a variable (e.g., “x_velocity”) and click “Apply”. Note that the default number of refinement levels loaded and vizualized is 1. Change to the required number of AMR level before clicking “Apply”.
For “Representation” select “Surface”.
For “Coloring” select the variable you chose above.
To add planes, near the top left you will see a cube icon with a green plane slicing through it. If you hover your mouse over it, it will say “Slice”. Click that button.
You can play with the Plane Parameters to define a plane of data to view, as shown in 1.

: Plotfile image generated with ParaView

VisIt

AMReX data can also be visualized by VisIt, an open source visualization and analysis software. To follow along with this example, first build and run the first heat equation tutorial code.

Next, download and install VisIt from https://wci.llnl.gov/simulation/computer-codes/visit. To open a single plotfile, run VisIt, then select “File” \(\rightarrow\) “Open file …”, then select the Header file associated the the plotfile of interest (e.g., plt00000/Header). Assuming you ran the simulation in 2D, here are instructions for making a simple plot:

To view the data, select “Add” \(\rightarrow\) “Pseudocolor” \(\rightarrow\) “phi”, and then select “Draw”.
To view the grid structure (not particularly interesting yet, but when we add AMR it will be), select “Add” \(\rightarrow\) “Subset” \(\rightarrow\) “levels”. Then double-click the text “Subset - levels”, enable the “Wireframe” option, select “Apply”, select “Dismiss”, and then select “Draw”.
To save the image, select “File” \(\rightarrow\) “Set save options”, then customize the image format to your liking, then click “Save”.

Your image should look similar to the left side of 1.

: 2D (left) and 3D (right) images generated using VisIt.

In 3D, you must apply the “Operators” \(\rightarrow\) “Slicing” \(\rightarrow\) “ThreeSlice”, with the “ThreeSlice operator attribute” set to x=0.25, y=0.25, and z=0.25. You can left-click and drag over the image to rotate the image to generate something similar to right side of 1.

To make a movie, you must first create a text file named movie.visit with a list of the Header files for the individual frames. This can most easily be done using the command:

~/amrex/Tutorials/Basic/HeatEquation_EX1_C> ls -1 plt*/Header | tee movie.visit
plt00000/Header
plt01000/Header
plt02000/Header
plt03000/Header
plt04000/Header
plt05000/Header
plt06000/Header
plt07000/Header
plt08000/Header
plt09000/Header
plt10000/Header

The next step is to run VisIt, select “File” \(\rightarrow\) “Open file…”, then select movie.visit. Create an image to your liking and press the “play” button on the VCR-like control panel to preview all the frames. To save the movie, choose “File” \(\rightarrow\) “Save movie …”, and follow the on-screen instructions.

Caveat:

The Visit reader determines “Cycle” from the name of the plotfile (directory), specifically from the integer that follows the string “plt” in the plotfile name.

So … if you call it plt00100 or myplt00100 or this_is_my_plt00100 then it will correctly recognize and print Cycle: 100.

If you call it plt00100_old it will also correctly recognize and print Cycle: 100

But, if you do not have “plt” followed immediately by the number, e.g. you name it pltx00100, then VisIt will not be able to correctly recognize and print the value for “Cycle”. (It will still read and display the data itself.)

yt

yt, an open source Python package available at http://yt-project.org/, can be used for analyzing and visualizing mesh and particle data generated by AMReX codes. Some of the AMReX developers are also yt project members. Below we describe how to use on both a local workstation, as well as at the NERSC HPC facility for high-throughput visualization of large data sets.

Note - AMReX datasets require yt version 3.4 or greater.

Using on a local workstation

Running yt on a local system generally provides good interactivity, but limited performance. Consequently, this configuration is best when doing exploratory visualization (e.g., experimenting with camera angles, lighting, and color schemes) of small data sets.

To use yt on an AMReX plot file, first start a Jupyter notebook or an IPython kernel, and import the yt module:

In [1]: import yt

In [2]: print(yt.__version__)
3.4-dev

Next, load a plot file; in this example we use a plot file from the Nyx cosmology application:

In [3]: ds = yt.load("plt00401")
yt : [INFO     ] 2017-05-23 10:03:56,182 Parameters: current_time              = 0.00605694344696544
yt : [INFO     ] 2017-05-23 10:03:56,182 Parameters: domain_dimensions         = [128 128 128]
yt : [INFO     ] 2017-05-23 10:03:56,182 Parameters: domain_left_edge          = [ 0.  0.  0.]
yt : [INFO     ] 2017-05-23 10:03:56,183 Parameters: domain_right_edge         = [ 14.24501  14.24501  14.24501]

In [4]: ds.field_list
Out[4]:
[('DM', 'particle_mass'),
 ('DM', 'particle_position_x'),
 ('DM', 'particle_position_y'),
 ('DM', 'particle_position_z'),
 ('DM', 'particle_velocity_x'),
 ('DM', 'particle_velocity_y'),
 ('DM', 'particle_velocity_z'),
 ('all', 'particle_mass'),
 ('all', 'particle_position_x'),
 ('all', 'particle_position_y'),
 ('all', 'particle_position_z'),
 ('all', 'particle_velocity_x'),
 ('all', 'particle_velocity_y'),
 ('all', 'particle_velocity_z'),
 ('boxlib', 'density'),
 ('boxlib', 'particle_mass_density')]

From here one can make slice plots, 3-D volume renderings, etc. An example of the slice plot feature is shown below:

In [9]: slc = yt.SlicePlot(ds, "z", "density")
yt : [INFO     ] 2017-05-23 10:08:25,358 xlim = 0.000000 14.245010
yt : [INFO     ] 2017-05-23 10:08:25,358 ylim = 0.000000 14.245010
yt : [INFO     ] 2017-05-23 10:08:25,359 xlim = 0.000000 14.245010
yt : [INFO     ] 2017-05-23 10:08:25,359 ylim = 0.000000 14.245010

In [10]: slc.show()

In [11]: slc.save()
yt : [INFO     ] 2017-05-23 10:08:34,021 Saving plot plt00401_Slice_z_density.png
Out[11]: ['plt00401_Slice_z_density.png']

The resulting image is 2. One can also make volume renderings with ; an example is show below:

: Slice plot of \(128^3\) Nyx simulation using yt.

In [12]: sc = yt.create_scene(ds, field="density", lens_type="perspective")

In [13]: source = sc[0]

In [14]: source.tfh.set_bounds((1e8, 1e15))

In [15]: source.tfh.set_log(True)

In [16]: source.tfh.grey_opacity = True

In [17]: sc.show()
<Scene Object>:
Sources:
    source_00: <Volume Source>:YTRegion (plt00401): , center=[  1.09888770e+25   1.09888770e+25   1.09888770e+25] cm, left_edge=[ 0.  0.  0.] cm, right_edge=[  2.19777540e+25   2.19777540e+25   2.19777540e+25] cm transfer_function:None
Camera:
    <Camera Object>:
    position:[ 14.24501  14.24501  14.24501] code_length
    focus:[ 7.122505  7.122505  7.122505] code_length
    north_vector:[ 0.81649658 -0.40824829 -0.40824829]
    width:[ 21.367515  21.367515  21.367515] code_length
    light:None
    resolution:(512, 512)
Lens: <Lens Object>:
    lens_type:perspective
    viewpoint:[ 0.95423473  0.95423473  0.95423473] code_length

In [19]: sc.save()
yt : [INFO     ] 2017-05-23 10:15:07,825 Rendering scene (Can take a while).
yt : [INFO     ] 2017-05-23 10:15:07,825 Creating volume
yt : [INFO     ] 2017-05-23 10:15:07,996 Creating transfer function
yt : [INFO     ] 2017-05-23 10:15:07,997 Calculating data bounds. This may take a while.
Set the TransferFunctionHelper.bounds to avoid this.
yt : [INFO     ] 2017-05-23 10:15:16,471 Saving render plt00401_Render_density.png

The output of this is 3.

Volume rendering of \(128^3\) Nyx simulation using yt. This corresponds to the same plot file used to generate the slice plot in 2.

Problem Setup

Each problem setup with a different initial e.g. temperature profile and bathymetry has its own subdirectory within Exec. To create a new problem, create a new subdirectory and copy the following files from another subdirectory:

prob.cpp
prob.H
inputs
GNUmakefile (change the directory in REMORA_PROBLEM_DIR to match)
CMakeLists.txt (change the problem name to match)
Make.package
amrvis.defaults (for visualization with AMRVis)

The file prob.cpp contains a number of functions that set the initial temeprature profile, as well as surface stress, mixing coefficients, and bathymetry. New problem-specific input parameters can be defined by adding a variable to the ProbParm class in prob.H, and reading in the value in amrex_probinit in prob.cpp. See the AMReX documentation on ParmParse for how to add parameters.

Variable Definitions

These variable definitions are as given in the ROMS documentation

Space and Time Variables
\(t\)	time
\(x,y\)	horizontal coordinates
\(z\)	vertical coordinate

Prescribed Inputs
\(h(x,y)\)	bottom depth
\(g\)	magnitude of gravity
\(f(x,y)\)	Coriolis parameter
\(\nu, \nu_\theta\)	molecular viscosity and diffusivity

Prognostic Variables
\(\vec{v} = (u,v,w)\)	velocity vector and x-, y-, z-components
\(T(x,y,z,t)\)	potential temperature
\(S(x,y,z,t)\)	salinity
\(\zeta(x,y,t)\)	surface elevation

Diagnostic Variables
\(\rho_0+\rho(x,y,z,t)\)	total in situ density
\(P\)	total pressure \(\approx -\rho_0 g z\)
\(\phi(x,y,z,t)\)	dynamic pressure \(\phi = (P/\rho_0)\)
\(K_m, K_C\)	vertical eddy viscosity and diffusivity

Additional terms
\(\cal{F}_u, \cal{F}_v, \cal{F}_C\)	forcing terms
\(\cal{D}_u, \cal{D}_v, \cal{D}_C\)	diffusive terms

Prognostic Equations

These equations are as given in the ROMS documentation

The momentum balance in the \(x\)- and \(y\)-directions are:

\[ \begin{align}\begin{aligned}{\frac {\partial u}{\partial t}}+{\vec {v}}\cdot \nabla u-fv=-{\frac {\partial \phi }{\partial x}}-{\partial \over \partial z}\left({\overline{u'w'}}-\nu {\partial u \over \partial z}\right)+{\cal {F}}_{u}+{\cal {D}}_{u}\\{\frac {\partial v}{\partial t}}+{\vec {v}}\cdot \nabla v+fu=-{\frac {\partial \phi }{\partial y}}-{\partial \over \partial z}\left({\overline{v'w'}}-\nu {\partial v \over \partial z}\right)+{\cal {F}}_{v}+{\cal {D}}_{v}\end{aligned}\end{align} \]

The time evolution of a scalar concentration field, \(C(x,y,z,t)\), e.g., salinity, temperature, or nutrient species, is governed by the advective-diffusive equation:

\[{\frac {\partial C}{\partial t}}+{\vec {v}}\cdot \nabla C=-{\partial \over \partial z}\left({\overline{C'w'}}-\nu _{\theta }{\partial C \over \partial z}\right)+{\cal {F}}_{C}+{\cal {D}}_{C}\]

The equation of state is given by:

\[\rho =\rho (T,S,P)\]

In the Boussinesq approximation, density variations are neglected in the momentum equations except in their contribution to the buoyancy force in the vertical momentum equation. Under the hydrostatic approximation, it is further assumed that the vertical pressure gradient balances the buoyancy force:

\[{\frac {\partial \phi }{\partial z}}=-{\frac {\rho g}{\rho _{o}}}\]

The final equation expresses the continuity equation for an incompressible fluid:

\[{\frac {\partial u}{\partial x}}+{\frac {\partial v}{\partial y}}+{\frac {\partial w}{\partial z}} = 0\]

In the ocean, vertical mixing due to molecular viscosity is extremely weak compared to the turbulent mixing, so the terms involving \(\nu\) and \(\nu_\theta\) can be neglected.

To close these questions, parametrizations of the Reynolds stresses and turbulent tracer fluxes are introduced as functions of the other variables and vertical turbulent eddy viscosity and eddy diffusivity coefficients \(K_m\) and \(K_C\), respectively.

\[\overline{u^\prime w^\prime} = -K_M \frac{\partial u}{\partial z}; \hspace{0.5in} \overline{v^\prime w^\prime} = -K_M \frac{\partial v}{\partial z}; \hspace{0.5in} \overline{C^\prime w^\prime} = -K_C \frac{\partial C}{\partial z};\]

An overbar represents a time average and a prime represents a fluctuation about the mean.

See Vertical Mixing for how \(K_m\) and \(K_C\) are computed in REMORA.

Vertical Mixing

These equations are as given in the ROMS documentation

Particles

REMORA has the option to include Lagrangian particles in addition to the mesh-based solution. Currently the particle functionality is very simple – the particles are initialized one per mesh cell in a particular plane, and are advected by the velocity field.

However, the AMReX particle data structure is very general and particles may take on a number of different roles in future.

To enable the use of particles, set

USE_PARTICLES = TRUE

in the GNUmakefile if using gmake, or add

-DREMORA_ENABLE_PARTICLES:BOOL=ON \

to the cmake command if using cmake. (See, e.g., Build/cmake_with_particles.sh)

One must also set

remora.use_tracer_particles = true

in the inputs file or on the command line at runtime.

Currently, by default, the particles are initialized at cell centers, one per cell when the cell index is (3,j,k), with zero initial velocity. They are advanced in time every time step.

Caveat: the particle information is currently output when using the AMReX-native plotfile format, but not when using netcdf. Writing particles into the netcdf files is a WIP.

To see an example of using the particle functionality, build the executable using gmake in Exec/ParticlesOverSeaMount.

To visualize the number of particles per cell as a mesh-based variable, add tracer_particle_count to the line in the inputs file

remora.plot_vars_1 =

Arakawa C-Grid

Variables are located on Arakawa C-grid as pictured in the images below.

Numerical Solution Technique

Horizontal Discretization

In the horizontal, the REMORA governing equations are discretized over a boundary-fitted, orthogonal curvilinear coordinates \(\left(\xi,\eta\right)\) grid. The general formulation of the curvilinear coordinates system allows Cartesian, polar, and spherical coordinates applications. The transformation of any of these coordinates to REMORA \(\left(\xi,\eta\right)\) grid is specified in the metric terms (pm, pn).

The model state variables are staggered using an Arakawa C-grid. As illustrated below, the free-surface (zeta), density (rho), and active/passive tracers (t) are located at the center of the cell whereas the horizontal velocity (u and v) are located at the west/east and south/north edges of the cell, respectively. That is, the density is evaluated between points where the currents are evaluated.

Staggered Horizontal Grid

In REMORA all the state arrays are dimensioned the same size to facilitate parallelization. However, the computational ranges for all the state variables are:

Grid Cell

Variable	Interior Range	Full Range
\(\rho\text{-type}\)	1:`Lm(ng)`,1:`Mm(ng)`	0:`L(ng)`,0:`M(ng)`
\(\psi\text{-type}\)	2:`Lm(ng)`,2:`Mm(ng)`	1:`L(ng)`,1:`M(ng)`
\(\text{u-type}\)	2:`Lm(ng)`,1:`Mm(ng)`	1:`L(ng)`,0:`M(ng)`
\(\text{v-type}\)	1:`Lm(ng)`,2:`Mm(ng)`	0:`L(ng)`,1:`M(ng)`

Vertical Discretization

The REMORA governing equations are discretized over variable topography using a stretched, terrain-following, vertical coordinate. As a result, each grid cell may have different level thickness (Hz) and volume. The model state variables are vertically staggered so that horizontal momentum (u, v), (rho), and active/passive tracers (t) are located at the center of the grid cell. The vertical velocity (omega, w) and vertical mixing variables (Akt, Akv, etc) are located at the bottom and top faces of the cell. See diagram below.

Vieste-Dubrovnik Transect

Staggered Vertical Grid

In this diagram, indices are 1-indexed (as in ROMS), while the indices in REMORA are 0-indexed.

The total thickness of the water column is \(\zeta\left(i,j\right)+h\left(i,j\right)\). The bathymetry (h) is usually time invariant whereas the free-surface (zeta) evolves in time. At input and output, the bathymetry is always a positive quantity. However, the depths z_r(i,j,k) and z_w(i,j,k) are negative for all locations below the mean sea level.

Grid Variables

Variable	Variable in Code	Definition	Location	Origin
m		x-dir coordinate	corners
n		y-dir coordinate	corners
\(\xi\)		transformed x-dir orthogonal curvilinear coordinate	corners
\(\eta\)		transformed y-dir orthogonal curvilinear coordinate	corners
\(\zeta\)		free-surface	center
\(\rho\)	`rho`	density	center	equation of state
\(t\)		active/passive tracers	center
\(u\)	`vec_Huon`	x-dir horizontal velocity	west/east faces
\(v\)	`vec_Hvom`	y-dir horizontal velocity	north/south faces
\(\overline{u}\)	`vec_ubar`	x-dir vertically integrated momentum	west/east faces
\(\overline{v}\)	`vec_vbar`	y-dir vertically integrated momentum	west/east faces
\(\psi\)			corners
\(H_z\)	`vec_Hz`	level thickness	center
\(\omega\)		vertical velocity	bottom/top faces
\(w\)		vertical velocity	bottom/top faces
Akt	`vec_Akv`	vertical mixing	bottom/top faces
Akv		vertical mixing	bottom/top faces
\(h\)	`vec_hOfTheConfusingName`	bathymetry (always positive)
\(z_{r\left(i,j,k\right)}\)		depth (negative below sea level)	center
\(z_{w\left(i,j,k\right)}\)		depth (negative below sea level)	bottom/top faces
\(T\)	`vec_t3`	temperature
\(S\)	`vec_s3`	salinity
\(D_{crit}\)		critical depth
\(D\)		total water depth
\(C\)		concentration
\(\delta_{\xi}\)		centered finite-difference approximation of \(\Delta\xi\)
\(\delta_{\eta}\)		centered finite-difference approximation of \(\Delta\eta\)
\(\delta_{\sigma}\)		centered finite-difference approximation of \(\Delta\sigma\)
\(\Delta_{\xi}\)		transformed x-dir differential distance
\(\Delta_{\eta}\)		transformed y-dir differential distance
\(\Delta_{\sigma}\)	`vec_s_r`	transformed z-dir differential distance
\(\Delta_{z}\)		vertical distance from one \(\rho\) to another
\(\overline{\left(\qquad\right)}^{\xi}\)		average taken over \(\Delta\xi\)
\(\overline{\left(\qquad\right)}^{\eta}\)		average taken over \(\Delta\eta\)
\(\overline{\left(\qquad\right)}^{\sigma}\)		average taken over \(\Delta\sigma\)
\(\Delta V\)		grid box volume

Conservation Properties

From Shchepetkin and McWilliams (2005), we have a tracer concentration equation in advective form:

()\[\frac{\partial C}{\partial t}+\left(u\cdot\nabla\right)C=0\]

and also a tracer concentration equation in conservation form:

()\[\frac{\partial C}{\partial t}+\nabla\cdot\left(uC\right)=0.\]

The continuity equation:

()\[\left(\nabla\cdot u\right)=0\]

can be used to get from one tracer equation to the other. As a consequence of eq. (1), if the tracer is spatially uniform, it will remain so regardless of the velocity field (constancy preservation). On the other hand, as a consequence of (2), the volume integral of the tracer concentration is conserved in the absence of internal sources and fluxes through the boundary. Both properties are valuable and should be retained when constructing numerical ocean models.

The semi-discrete form of the tracer equation is:

()\[\begin{split}\frac{\partial}{\partial t}\left(\frac{H_zC}{mn}\right)+\delta_{\xi}\left(\frac{u\overline{H_z}^{\xi}\overline{C}^{\xi}}{\overline{n}^{\xi}}\right)+\delta_{\eta}\left(\frac{v\overline{H_z}^{\eta}\overline{C}^{\eta}}{\overline{m}^{\eta}}\right)+\delta_{\sigma}\left(\overline{C}^{\sigma}\frac{H_z\Omega}{mn}\right)= \\ \frac{1}{mn}\frac{\partial}{\partial\sigma}\left(\frac{K_m}{\Delta z}\frac{\partial C}{\partial\sigma}\right)+\mathcal{D}_C+\mathcal{F}_C\end{split}\]

Here \(\delta_{\xi},\delta_{\eta}\) and \(\delta_{\sigma}\) denote simple centered finite-difference approximations to \(\partial/\partial\xi,\partial/\partial\eta\) and \(\partial/\partial\sigma\) with the differences taken over the distances \(\Delta\xi,\Delta\eta\) and \(\Delta\sigma\), respectively. \(\Delta z\) is the vertical distance from one \(\rho\) point to another. \(\overline{\left(\qquad\right)}^{\xi}, \overline{\left(\qquad\right)}^{\eta}\) and \(\overline{\left(\qquad\right)}^{\sigma}\) represent averages taken over the distances \(\Delta\xi,\Delta\eta\) and \(\Delta\sigma\).

The finite volume version of the same equation is no different, except that a quantity \(C\) is defined as the volume-averaged concentration over the grid box \(\Delta V\):

\[C=\frac{mn}{H_z}\int_{\Delta V}\frac{H_z C}{mn}\delta\xi\ \delta\eta\ \delta\sigma\]

The quantity \(\left(\dfrac{u\overline{H_z}^{\xi}\overline{C}^{\xi}}{\overline{n}^{\xi}}\right)\) is the flux through an interface between adjacent grid boxes.

This method of averaging was chosen because it internally conserves first moments in the model domain, although it is still possible to exchange mass and energy through the open boundaries. The method is similar to that used in Arakawa and Lamb; though their scheme also conserves enstrophy. Instead, we will focus on (nearly) retaining constancy preservation while coupling the barotropic (depth-integrated) equations and the baroclinic equations.

The timestep in eq. (4) is assumed to be from time \(n\) to time \(n+1\), while the other terms being evaluated at time \(n+1/2\) for second-order accuracy. Setting \(C\) to \(1\) everywhere reduces eq. (4) to:

()\[\frac{\partial}{\partial t}\left(\frac{H_z}{mn}\right)+\delta_{\xi}\left(\frac{u\overline{H_z}^{\xi}}{\overline{n}^{\xi}}\right)+\delta_{\eta}\left(\frac{v\overline{H_z}^{\eta}}{\overline{m}^{\eta}}\right)+\delta_{\sigma}\left(\frac{H_z\Omega}{mn}\right)=0\]

If this equation holds true for the step from time \(n\) to time \(n+1\), then constancy preservation will hold.

In a hydrostatic model such as REMORA, the discrete continuity equation is needed to compute vertical velocity rather than grid-box volume \(\dfrac{H_z}{mn}\) (the latter is controlled by changes in \(\zeta\) in the barotropic mode computations). Here, \(\dfrac{H_z\Omega}{mn}\) is the finite-volume flux across the moving grid-box interface, vertically on the \(w\) grid.

The vertical integral of the continuity equation (5), using the vertical boundary conditions on \(\Omega\), is:

()\[\frac{\partial}{\partial t}\left(\frac{\zeta}{mn}\right)+\delta_{\xi}\left(\frac{\overline{u}\overline{D}^{\xi}}{\overline{n}^{\xi}}\right)+\delta_{\eta}\left(\frac{\overline{v}\overline{D}^{\eta}}{\overline{m}^{\eta}}\right)=0\]

where \(\zeta\) is the surface elevation, \(D=h+\zeta\) is the total depth, and \(\overline{u},\overline{v}\) are the depth-integrated horizontal velocities. This equation and the corresponding 2-D momentum equations are timestepped on a shorter timestep than eq.(4) and the other 3-D equations. Due to the details in the mode coupling, it is only possible to maintain constancy preservation to the accuracy of the barotropic timesteps.

Depth-Integrated Equations

The depth average of a quantity \(A\) is given by:

()\[\overline{A}=\frac{1}{D}\int_{-1}^0H_zA\ d\sigma\]

where the overbar indicates a vertically averaged quantity and

()\[D\equiv\zeta\left(\xi,\eta,t\right)+h\left(\xi,\eta\right)\]

is the total depth of the water column. The vertical integral of the momentum equations are:

()\[\begin{split}\frac{\partial}{\partial t}\left(\frac{D\overline{u}}{mn}\right)+\frac{\partial}{\partial\xi}\left(\frac{D\overline{uu}}{n}\right)+\frac{\partial}{\partial\eta}\left(\frac{D\overline{uv}}{m}\right)-&\frac{Df\overline{v}}{mn}\\ -\left[\overline{vv}\frac{\partial}{\partial\xi}\left(\frac{1}{n}\right)-\overline{uv}\frac{\partial}{\partial\eta}\left(\frac{1}{m}\right)\right]D=-\frac{D}{n}&\left(\frac{\partial\overline{\phi_2}}{\partial\xi}+g\frac{\partial\zeta}{\partial\xi}\right)\\ +\frac{D}{mn}\left(\overline{\mathcal{F}}_u+\overline{\mathcal{D}}_{h_u}\right)&+\frac{1}{mn}\left(\tau^{\xi}_s-\tau^{\xi}_b\right)\end{split}\]

and

()\[\begin{split}\frac{\partial}{\partial t}\left(\frac{D\overline{v}}{mn}\right)+\frac{\partial}{\partial\xi}\left(\frac{D\overline{uv}}{n}\right)+\frac{\partial}{\partial\eta}\left(\frac{D\overline{vv}}{m}\right)+&\frac{Df\overline{u}}{mn}\\ +\left[\overline{uv}\frac{\partial}{\partial\xi}\left(\frac{1}{n}\right)-\overline{uu}\frac{\partial}{\partial\eta}\left(\frac{1}{m}\right)\right]D=-\frac{D}{m}&\left(\frac{\partial\overline{\phi_2}}{\partial\eta}+g\frac{\partial\zeta}{\partial\eta}\right)\\ +\frac{D}{mn}\left(\overline{\mathcal{F}}_v+\overline{\mathcal{D}}_{h_v}\right)+&\frac{1}{mn}\left(\tau_s^{\eta}-\tau_b^{\eta}\right)\end{split}\]

where \(\phi_2\) includes the \(\frac{\partial z}{\partial\xi}\) term, \(\overline{\mathcal{D}}_{h_u}\) is the horizontal viscosity, and the vertical viscosity only contributes through the upper and lower boundary conditions. We also need the vertical integral of the continuity equation, shown above as eq. (6). The presence of a free surface introduces waves which propagate at a speed of \(\sqrt{gh}\). These waves usually impose a more severe time-step limit than any of the internal processes. We have therefore chosen to solve the full equations by means of a split time step. In other words, the depth integrated equations (9), (10), and (6) are integrated using a short time step and the values of \(\overline{u}\) and \(\overline{v}\) are used to replace those found by integrating the full equations on a longer time step. A diagram of the barotropic time stepping is shown here:

Some of the terms in equations (9) and (10) are updated on the short time step while others are not. The contributions from the slow terms are computed once per long time step and stored. If we call these terms \(R_{u_{\text{slow}}}\) and \(R_{v_{\text{slow}}}\), equations (9) and (10) become:

()\[\begin{split}\frac{\partial}{\partial t}\left(\frac{D\overline{u}}{mn}\right)+\frac{\partial}{\partial\xi}\left(\frac{D\overline{u}\,\overline{u}}{n}\right)+\frac{\partial}{\partial\eta}\left(\frac{D\overline{u}\,\overline{v}}{m}\right)&-\frac{Df\overline{v}}{mn}\\ -\left[\overline{v}\,\overline{v}\frac{\partial}{\partial\xi}\left(\frac{1}{n}\right)-\overline{u}\,\overline{v}\frac{\partial}{\partial\eta}\left(\frac{1}{m}\right)\right]D&=R_{u_{\text{slow}}}-\frac{gD}{n}\frac{\partial\zeta}{\partial\xi}+\frac{D}{mn}\mathcal{D}_{\overline{u}}-\frac{1}{mn}\tau_{b}^{\xi}\\\end{split}\]

and

()\[\begin{split}\frac{\partial}{\partial t}\left(\frac{D\overline{v}}{mn}\right)+\frac{\partial}{\partial\xi}\left(\frac{D\overline{u}\,\overline{v}}{n}\right)+\frac{\partial}{\partial\eta}\left(\frac{D\overline{v}\,\overline{v}}{m}\right)&+\frac{Df\overline{u}}{mn}\\ +\left[\overline{u}\,\overline{v}\frac{\partial}{\partial\xi}\left(\frac{1}{n}\right)-\overline{u}\,\overline{u}\frac{\partial}{\partial\eta}\left(\frac{1}{m}\right)\right]D&=R_{v_{\text{slow}}}-\frac{gD}{m}\frac{\partial\zeta}{\partial\eta}+\frac{D}{mn}\mathcal{D}_{\overline{v}}-\frac{1}{mn}\tau_{b}^{\eta}.\end{split}\]

When time stepping the model, we compute the right-hand-sides for the full 3-D momentum equations as well as the right-hand-sides for equations (11) and (12). The vertical integral of the 3-D right-hand-sides are obtained and then the 2-D right-hand-sides are subtracted. The resulting fields are the slow forcings \(R_{u_{\text{slow}}}\) and \(R_{v_{\text{slow}}}\). This was found to be the easiest way to retain the baroclinic contributions of the non-linear terms such as \(\overline{uu}-\overline{u}\,\overline{u}\). The model is time stepped from time \(n\) to time \(n+1\) by using short time steps on equations (11), (12) and (6). Equation (6) is time stepped first, so that an estimate of the new \(D\) is available for the time rate of change terms in equations (11) and (12). A third-order predictor-corrector time stepping is used. In practice, we actually time step all the way to time \(\left(n+\textbf{dtfast}\times M^*\right)\) and while maintaining weighted averages of the values of \(\overline{u},\overline{v}\) and \(\zeta\). The averages are used to replace the values at time \(n+1\) in both the baroclinic and barotropic modes, and for recomputing the vertial grid spacing \(H_z\). The following figure shows one option for how these weights might look:

The primary weights, \(a_m\), are used to compute \(\langle\zeta\rangle^{n+1}\equiv\overunderset{M^*}{m=1}{\sum}a_m\zeta^m\). There is a related set of secondary weights \(b_m\), used as \(\langle\!\langle\overline{u}\rangle\!\rangle^{n+\frac{1}{2}}\equiv\overunderset{M^*}{m=1}{\sum}b_m\overline{u}^m\). In order to maintain constancy preservation, this relation must hold:

()\[\begin{split}\langle\zeta\rangle_{i,j}^{n+1}=\langle\zeta\rangle_{i,j}^n&-\\ \left(mn\right)_{i,j}\Delta t&\left[\left\langle\!\!\left\langle\frac{D\overline{u}}{n}\right\rangle\!\!\right\rangle_{i+\frac{1}{2},j}^{n+\frac{1}{2}}-\left\langle\!\!\left\langle\frac{D\overline{u}}{n}\right\rangle\!\!\right\rangle_{i-\frac{1}{2},j}^{n+\frac{1}{2}}+\left\langle\!\!\left\langle\frac{D\overline{v}}{m}\right\rangle\!\!\right\rangle_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}-\left\langle\!\!\left\langle\frac{D\overline{v}}{m}\right\rangle\!\!\right\rangle_{i,j-\frac{1}{2}}^{n+\frac{1}{2}}\right]\end{split}\]

Shchepetkin and McWilliams (2005) introduce a range of possible weights, but the ones used here have a shape function:

()\[A\left(\tau\right)=A_0\left\{\left(\frac{\tau}{\tau_0}\right)^p\left[1-\left(\frac{\tau}{\tau_0}\right)^q\right]-r\frac{\tau}{\tau_0}\right\}\]

where \(p,q\) are parameters and \(A_0,\tau_0\), and \(r\) are chosen to satisfy normalization, consistency, and second-order accuracy conditions,

()\[I_n=\int_0^{\tau^*}\tau^nA\left(\tau\right)d\tau=1,\qquad n=0,1,2\]

using Newton iterations. \(\tau^*\) is the upper limit of \(\tau\) with \(A\left(\tau\right)\geq0\). In practice we initially set

\[A_0=1,r=0,\text{ and }\tau=\frac{\left(p+2\right)\left(p+q+2\right)}{\left(p+1\right)\left(p+q+1\right)}\]

compute \(A\left(\tau\right)\) using eq.(14), normalize using:

()\[\sum_{m=1}^{M^*}a_m\equiv1,\qquad\sum_{m=1}^{M^*}a_m\frac{m}{M}\equiv1,\]

and adjust \(r\) iteratively to satisfy the \(n=2\) condition of (15). We are using values of \(p=2,q=4\), and \(r=0.284\). This form allows some negative weights for small \(m\), allowing \(M^*\) to be less than \(1.5M\).

Pressure Gradient Terms in Mode Coupling

Equation (11) contains the term \(R_{u_{\text{slow}}}\), computed as the difference between the 3-D right-hand-side and the 2-D right-hand-side. The pressure gradient therefore has the form:

()\[-\frac{gD}{n}\frac{\partial\zeta}{\partial\xi}+\left[\frac{gD}{n}\frac{\partial\zeta}{\partial\xi}+\mathcal{F}\right]\]

where the term in square brackets is the mode coupling term and is held fixed over all the barotropic steps and

()\[\mathcal{F}=-\frac{1}{\rho_0n}\int_{-h}^{\zeta}\frac{\partial P}{\partial\xi}dz\]

is the vertically integrated pressure gradient. The latter is a function of the bathymetry, free surface gradient, and the free surface itself, as well as the vertical distribution of density.

The disadvantage of this approach is that after the barotropic time stepping is complete and the new free surface is substituted into the full baroclinic pressure gradient, its vertical integral will no longer be equal to the sum of the new surface slope term and the original coupling term based on the old free surface. This is one form of mode-splitting error which can lead to trouble because the vertically integrated pressure gradient is not in balance with the barotropic mass flux.

Instead, let us define the following:

()\[\overline{\rho}=\frac{1}{D}\int_{-h}^{\zeta}\rho\ dz,\qquad\rho^*=\frac{1}{\frac{1}{2}D^2}\int_{-h}^{\zeta}\left\{\int_z^{\zeta}\rho\ dz^{\prime}\right\}\ dz\]

Changing the vertical coordinate to \(\sigma\) yields:

()\[\overline{\rho}=\int_{-1}^0\rho\ d\sigma,\qquad\rho^*=2\int_{-1}^0\left\{\int_{\sigma}^0\rho\ d\sigma^{\prime}\right\}\ d\sigma\]

which implies that \(\overline{\rho}\) and \(\rho^*\) are actually independent of \(\zeta\) as long as the density profile \(\rho=\rho\left(\sigma\right)\) does not change. The vertically integrated pressure gradient becomes:

()\[-\frac{1}{\rho_0}\frac{g}{n}\left\{\frac{\partial}{\partial\xi}\left(\frac{\rho^*D^2}{2}\right)-\overline{\rho}D\frac{\partial h}{\partial\xi}\right\}=-\frac{1}{\rho_0}\frac{g}{n}D\left\{\rho^*\frac{\partial\zeta}{\partial\xi}+\frac{D}{2}\frac{\partial\rho^*}{\partial\xi}+\left(\rho^*-\overline{\rho}\right)\frac{\partial h}{\partial\xi}\right\}\]

In the case of uniform density \(\rho_0\), we obtain \(\rho^*\equiv\overline{\rho}\equiv\rho_0\), but we otherwise have two new terms. The accuracy of these terms depends on an accurate vertical integration of the density, as described in Shchepetkin and McWilliams (2005).

Horizontal and Vertical Advection

The advection of a tracer \(C\) has an equation of the form

\[\frac{\partial}{\partial t}\frac{H_zC}{mn}=-\frac{\partial}{\partial\xi}F^{\xi}-\frac{\partial}{\partial\eta}F^{\eta}-\frac{\partial}{\partial\sigma} F^{\sigma},\]

where we have introduced the advective fluxes:

\[\begin{split}F^{\xi}&=\frac{H_zuC}{n}\\ F^{\eta}&=\frac{H_zvC}{m}\\ F^{\sigma}&=\frac{H_z\Omega C}{mn}.\end{split}\]

Fourth-order Centered

The barotropic advection is centered fourth-order. Create gradient terms:

\[\begin{split}G^{\xi}&=\overline{\left(\frac{\partial C}{\partial\xi}\right)}^{\xi}\\ G^{\eta}&=\overline{\left(\frac{\partial C}{\partial\eta}\right)}^{\eta}\\ G^{\sigma}&=\overline{\left(\frac{\partial C}{\partial\sigma}\right)}^{\sigma}.\end{split}\]

The fluxes now become:

\[\begin{split}F^{\xi}&=\frac{\overline{H_z}^{\xi}}{\overline{n}^{\xi}}u\left(\overline{C}^{\xi}-\frac{1}{3}\frac{\partial G^{\xi}}{\partial\xi}\right)\\ F^{\eta}&=\frac{\overline{H_z}^{\eta}}{\overline{m}^{\eta}}v\left(\overline{C}^{\eta}-\frac{1}{3}\frac{\partial G^{\eta}}{\partial\eta}\right)\\ F^{\sigma}&=\frac{\overline{H_z}^{\sigma}}{mn}\Omega\left(\overline{C}^{\sigma}-\frac{1}{3}\frac{\partial G^{\sigma}}{\partial\sigma}\right).\end{split}\]

Third-order Upwind

There is a class of third-order upwind advection schemes, both one-dimensional (Leanord, 1979) and two-dimensional (Rasch, 1994 and Shchepetkin and McWilliams, 1998). This scheme is known as UTOPIA (Uniformly Third-Order Polynomial Interpolation Algorithm). Applying flux limiters to UTOPIA is explored in Thuburn (1995), although it is not implemented in REMORA. The two-dimensional formulation in Rasch contains terms of order \(u^2C\) and \(u^3C\), including cross terms (\(uvC\)). The terms which are nonlinear in velocity have been dropped in REMORA, leaving one extra upwind term in the computation of the advective fluxes:

\[\begin{split}F^{\xi}&=\frac{H_zu}{n}\left(C-\gamma\frac{\partial^2C}{\partial\xi^2}\right)\\ F^{\eta}&=\frac{H_zv}{m}\left(C-\gamma\frac{\partial^2C}{\partial\eta^2}\right)\end{split}\]

The second derivative terms are centered on a \(\rho\) point in the grid, but are needed at a \(u\) or \(v\) point in the flux. The upstream value is used:

\[F^{\xi}_{i,j,k}=\frac{\overline{H_z}^{\eta}}{\overline{n}^{\xi}}\left[\max\left(0,u_{i,j,k}\right)C_{i-1,j,k}+\min\left(0,u_{i,j,k}\right)C_{i,j,k}\right].\]

The value of \(\gamma\) in the model is \(\frac{1}{8}\) while that in Rasch(1994) is \(\frac{1}{6}\).

Because the third-order upwind scheme is designed to be two-dimensional, it is not used in the vertical (though one might argue that we are simply performing one-dimensional operations here). Instead we use a centered fourth-order scheme in the vertical when the third-order upwind option is turned on:

\[F^s=\frac{H_zw}{mn}\left[-\frac{1}{16}C_{i,j,k-1}+\frac{9}{16}C_{i,j,k}+\frac{9}{16}C_{i,j,k+1}-\frac{1}{16}C_{i,j,k+2}\right]\]

One advantage of UTOPIA over MPDATA is that it can be used on variables having both negative and positive values. Therefore, it can be used on velocity as well as scalars (is there a reference for this?). For the \(u\)-velocity, we have:

\[\begin{split}F^{\xi}&=\left(u-\gamma\frac{\partial^2u}{\partial\xi^2}\right)\left[\frac{H_zu}{n}-\gamma\frac{\partial^2}{\partial\xi^2}\left(\frac{H_zu}{n}\right)\right]\\ F^{\eta}&=\left(u-\gamma\frac{\partial^2u}{\partial\eta^2}\right)\left[\frac{H_zv}{m}-\gamma\frac{\partial^2}{\partial\xi^2}\left(\frac{H_zv}{m}\right)\right]\\ F^{\sigma}&=\frac{H_zw}{mn}\left[-\frac{1}{16}u_{i,j,k-1}+\frac{9}{16}u_{i,j,k}+\frac{9}{16}u_{i,j,k+1}-\frac{1}{16}u_{i,j,k+2}\right]\end{split}\]

while for the \(v\)-velocity we have:

\[\begin{split}F^{\xi}&=\left(v-\gamma\frac{\partial^2v}{\partial\xi^2}\right)\left[\frac{H_zu}{n}-\gamma\frac{\partial^2}{\partial\eta^2}\left(\frac{H_zu}{n}\right)\right]\\ F^{\eta}&=\left(v-\gamma\frac{\partial^2v}{\partial\eta ^2}\right)\left[\frac{H_zv}{m}-\gamma\frac{\partial^2}{\partial\eta^2}\left(\frac{H_zv}{m}\right)\right]\\ F^{\sigma}&=\frac{H_zw}{mn}\left[-\frac{1}{16}v_{i,j,k-1}+\frac{9}{16}v_{i,j,k}+\frac{9}{16}v_{i,j,k+1}-\frac{1}{16}v_{i,j,k+2}\right]\end{split}\]

In all these terms, the second derivatives are evaluated at an upstream location.

Vertical Velocity

Having obtained a complete specification of the \(u,v,T\), and \(S\) fields at the next time level by the methods outlined above, the vertical velocity and density fields can be calculated. the vertical velocity is obtained by rewriting equation (5) as:

\[\frac{\partial}{\partial t}\left(\frac{\zeta}{mn}\right)+\delta_{\xi}\left(\frac{u\overline{H_z}^{\xi}}{\overline{n}^{\xi}}\right)+\delta_{\eta}\left(\frac{v\overline{H_z}^{\eta}}{\overline{m}^{\eta}}\right)+\delta_{\sigma}\left(\frac{H_z\Omega}{mn}\right)=0\]

and combining it with equation (6) to obtain:

\[\frac{\partial}{\partial\xi}\left(\frac{H_zu}{n}\right)+\frac{\partial}{\partial\eta}\left(\frac{H_zv}{m}\right)+\frac{\partial}{\partial\sigma}\left(\frac{H_z\Omega}{mn}\right)-\frac{\partial}{\partial\xi}\left(\frac{D\overline{u}}{n}\right)-\frac{\partial}{\partial\eta}\left(\frac{D\overline{v}}{m}\right)=0.\]

Solving for \(H_z\Omega/mn\) and using the semi-discrete notation we obtain:

\[\frac{H_z\Omega}{mn}=\int\left[\delta_{\xi}\left(\frac{\overline{u}\,\overline{D}^{\xi}}{\overline{n}^{\xi}}\right)+\delta_{\eta}\left(\frac{\overline{v}\,\overline{D}^{\eta}}{\overline{m}^{\eta}}\right)-\delta_{\xi}\left(\frac{u\overline{H_z}^{\xi}}{\overline{n}^{\xi}}\right)-\delta_{\eta}\left(\frac{v\overline{H_z}^{\eta}}{\overline{m}^{\eta}}\right)\right]\ d\sigma.\]

The integral is actually computed as a sum from the bottom upwards and also as a sum from the top downwards. The value used is a linear combination of the two, weighted so that the surface down value is used near the surface while the other is used near the bottom. [Is this still done?]

Equation of State

The density is obtained from temperature and salinity via an equation of state. REMORA provides a choice of a nonlinear equation of state \(\rho=\rho\left(T,S,z\right)\) or a linear equation of state \(\rho=\rho\left(T\right)\). The nonlinear equation of state has been modified and now corresponds to the UNESCO equation of state as derived by Jackett and McDougall (1995). It computes in situ density as a function of potential temperature, salinity and pressure.

Warning: although we have used it quite extensively in the past, McDougall (personal communication) claims that the single-variable \(\left(\rho=\rho\left(T\right)\right)\) equation of state is not dynamically appropriate as is. He has worked out the extra source and sink terms required, arising from vertical motions and the compressibility of water. They are quite complicated and we have not implemented them to see if they alter the flow.

Time Stepping Overview

While time-stepping the model, we have a stored history of the model fields at time \(n-1\), an estimate of the fields at the current time \(n\), and we need to come up with an estimate for time \(n+1\). For reasons of efficiency, we choose to use a split-explicit time-step, integrating the depth-integrated equations with a shorter time-step than the full 3-D equations. There is an integer ratio \(M\) between the time-steps. The exact details of how the time-stepping is done vary from one version of ROMS to the next, with the east coast ROMS described here being older than other branches. Still, most versions have these steps:

Great care is taken to avoid the introduction of a mode-splitting instability due to the use of shorter time steps for the depth-integrated computations. The mode coupling has evolved through the various ROMS versions, as shown in the figures below [from Shchepetkin and McWilliams (2008a)]. The time-stepping schemes are also listed in the following table and described in detail in Shchepetkin and Williams (2005) and Shchepetkin and McWilliams (2009); the relevant ones are described in Time-stepping Schemes Review. REMORA uses the Rutgers Scheme.

The following figure shows the time-stepping and mode coupling used in REMORA. The curved arrows update the 3-D fields; those with “pillars” are leapfrog in nature with the pillar representing the r.h.s. terms. Straight arrows indicate exchange between the barotropic and baroclinic modes. The shape functions for the fast time steps show just one option out of many possibilities. The grey function has weights to produce an estimate at time \(n+1\), while the light red function has weights to produce an estimate at time \(n+1/2\).

Rutgers ROMS:

Timestepping table from Shchepetkin and McWilliams (2008a):

	SCRUM 3.0	Rutgers	AGRIF	UCLA	Non-hydrostatic
barotropic mode	LF-TR	LF-AM3 with FB feedback	LF-AM3 with FB feedback \(\dagger\)	Gen. FB (AB3-AM4)	Gen. FB (AB3-AM4)
2-D \(\alpha_{\max}\),iter.	\(\sqrt{2}, (2)\ddagger\)	\(1.85, (2)\)	\(1.85, (2)\)	\(1.78, (1)\)	\(1.78, (1)\)
3-D momenta	AB3	AB3	LF-AM3	LF-AM3	AB3 (mod)
Tracers	AB3	LF-TR	LF-AM3	LF-AM3	AB3 (mod)
internal waves	AB3	Gen. FB (AB3-TR)	LF-AM3, FB feedback	LF-AM3 FB feedback	Gen. FB (AB3-AM4)
\(\alpha_{\max}\),advect.	\(0.72\)	\(0.72\)	\(1.587\)	\(1.587\)	\(0.78\)
\(\alpha_{\max}\),Cor.	\(0.72\)	\(0.72\)	\(1.587\)	\(1.587\)	\(0.78\)
\(\alpha_{\max}\),int. w.	\(0.72, (1)\)	\(1.14, (1,2)\)	\(1.85, (2)\)	\(1.85, (2)\)	\(1.78, (1)\)

Note: \(\dagger\) The generalized FB barotropic mode was ported into the newest AGRIF code at the end of 2007. \(\ddagger\) The number in parentheses (e.g., 2) indicates the number of r.h.s. computations per time step. If there are two parenthesized numbers, the first one is for momenta, the second for tracers.

Time Stepping: Internal Velocity Modes and Tracers

The momentum equations are advanced before the tracer equation, by computing all the terms except the vertical viscosity and then using the implicit scheme described in #Vertical Friction and Diffusion to find the new values for \(u\) and \(v\). The depth-averaged component is then removed and replaced by the \(\langle \overline{u} \rangle\) and \(\langle \overline{v} \rangle\) computed as in #Depth-Integrated Equations. A third-order Adams-Bashforth (AB3) time step is used, requiring multiple right-hand-side time levels. These stored up r.h.s. values can be used to extrapolate to a value at time \(n+\frac{1}{2}\), obtained from the predictor step. The vertical diffusion is computed as in #Vertical Friction and Diffusion. The predictor step cannot be both constancy=preserving and conservative; it was therefore decided to make it constancy-preserving. Also, since it is only being used to compute the advection for the corrector step, the expensive diffusion operations are not carried out on the predictor step.

The preceeding notes on tracer advection refer to all but the MPDATA option. The MPDATA algorithm has its own predictor-corrector with emphasis on not allowing values to exceed their original range, and therefore gives up the constancy-preservation. This will be most noticeable in shallow areas with large tides.

Vertical S-coordinate

These equations are as given in the ROMS documentation

REMORA has a generalized vertical, terrain-following, coordinate system. Currently, the vertical transformation equation \(z=z \left( x,y, \sigma ,t \right)\), (Vtransform=2 from ROMS) is available which can support vertical stretching 1D-functions when several constraints are satisfied.

Transformation Equation

\[\begin{split}z\left( x,y,\sigma ,t\right) &=\zeta \left( x,y,t\right) +\left[ \zeta \left( x,y,t\right) +h\left( x,y\right) \right] S\left( x,y,\sigma \right) ,\\ S\left( x,y,\sigma \right) &=\frac{h_c\ \sigma +h \left( x,y\right) C\left( \sigma \right) }{h_c +h\left( x,y\right)}\end{split}\]

where \(S\left( x,y,\sigma \right)\) is a nonlinear vertical transformation functional, \(\zeta \left( x,y,t\right)\) is the time-varying free-surface, \(h\left( x,y\right)\) is the unperturbed water column thickness and \(z=-h\left( x,y\right)\) corresponds to the ocean bottom, \(\sigma\) is a fractional vertical stretching coordinate ranging from \(-1\leq \sigma \leq 0,\ C\left( \sigma \right)\) is a nondimensional, monotonic, vertical stretching function ranging from \(-1\leq C\left( \sigma \right) \leq 0\), and \(h_c\) is a positive thickness controlling the stretching.

We find it convenient to define:

\[H_z \equiv \frac{\partial z}{\partial \sigma }\]

where \(H_z = H_z \left( x,y,\sigma ,t\right)\) are the vertical grid thicknesses. In REMORA, \(H_z\) is computed discretely as \(\Delta z/\Delta \sigma\) since this leads to the vertical sum of \(H_z\) being exactly the total water column thickness \(D\).

Notice that,

\[\begin{split}S\left( x,y,\sigma \right) = \begin{cases} 0, & \text{if } \sigma = 0, & C\left( \sigma \right) = 0, & \text{at the free-surface;}\\ -1 , & \text{if } \sigma = -1, & C\left( \sigma \right) = -1, & \text{at the ocean bottom} \end{cases}\end{split}\]

This transformation offers several advantages:

Regardless of the design of \(C\left( \sigma \right)\), it behaves like equally-spaced sigma-coordinates in shallow regions, where \(h\left( x,y\right) \ll h_c\). This is advantageous because it avoids excessive resolution and associated CFL limitation in such areas.
Near-surface refinement behaves more or less like geopotential coordinates in deep regions (level thicknesses, \(H_z\), do not depend or weakly depend on bathymetry), while near-bottom like sigma coordinates (\(H_z\) is roughly proportional to depth). This reduces the extreme r-factors near the bottom and reduces pressure gradient errors.
The true sigma-coordinate system is recovered as \(h_c \rightarrow \infty\). This is useful when configuring applications with flat bathymetry and uniform level thickness. Practically, you can achieve this by setting tcline to 1.0d+16 in DataStruct.H. This will set \(h_c =1.0\times 10^{16}\). Although not necessary, we also recommend that you set \(\theta _S = 0\) and \(\theta _B =0\).

In an undisturbed ocean state, \(\zeta \equiv 0\), this transformation yields the following unperturbed depths, \(\hat{z}\),

\[\hat{z} \left( x,y,\sigma \right) \equiv h\left( x,y\right) S\left( x,y,\sigma \right) =h\left( x,y\right) \left[ \frac{h_c\ \sigma +h\left( x,y\right) C\left( \sigma \right)}{h_c +h\left( x,y\right) } \right]\]

and

\[d\hat{z} =d\sigma \ h\left( x,y\right) \left[ \frac{h_c}{h_c +h\left( x,y\right) } \right]\]

As a consequence, the uppermost grid box retains very little dependency from bathymetry in deep areas, where \(h_c \ll h\left( x,y\right)\). For example if \(h_c =250m\) and \(h\left( x,y\right)\) changes from \(2000\) to \(6000m\), the uppermost grid box changes only by a factor of 1.08 (less than 10%).

Vertical Stretching

The above generic vertical transformation design facilitates numerous vertical stretching functions \(C \left( \sigma \right)\) . This function is defined in terms of several parameters which are specified in the standard input file, inputs.

Stretching Function Properties

\(C\left( \sigma \right)\) is a dimensionless, nonlinear, monotonic function;
\(C\left( \sigma \right)\) is a continuous differentiable function, or a differentiable piecewise function with smooth transition;
\(C\left( \sigma \right)\) must be discretized in terms of the fractional stretched vertical coordinate \(\sigma\),

\[\begin{split}\sigma \left( k \right) = \begin{cases} \frac{k-N}{N}, & \text{at vertical }W\text{-points}, & k=0,\ldots ,N,\\ \frac{k-N-0.5}{N}, & \text{at vertical}\rho \text{-points}, & k=1,\ldots ,N \end{cases}\end{split}\]
\(C\left( \sigma \right)\) must be constrained by \(-1 \leq C\left( \sigma \right) \leq 0\), that is,

\[\begin{split}C\left( \sigma \right) = \begin{cases} 0, & \text{if } \sigma = 0, & C\left( \sigma \right) = 0, & \text{at the free-surface};\\ -1, & \text{if } \sigma = -1, & C\left( \sigma \right) = -1, & \text{at the ocean bottom}. \end{cases}\end{split}\]

Available Stretching Functions

(A. Shchepetkin (2010) Vstretching=4 from ROMS). \(C\left( \sigma \right)\) is defined as a continuous, double stretching function:

Surface refinement function:

\[C\left( \sigma \right) = \frac{1-\cosh \left( \theta _S \sigma \right) }{\cosh \left( \theta _S \right) -1}, \quad \text{for } \theta _S > 0, \quad C\left( \sigma \right) = - \sigma ^2, \quad \text{for }\theta _S \leq 0\]

Bottom refinement function:

\[C\left( \sigma \right) = \frac{\exp \left( \theta _B C\left( \sigma \right) \right) -1}{1-\exp \left( -\theta _B \right) }, \quad \text{for }\theta _B >0\]

Notice that the bottom function is the second stretching of an already stretched transform. The resulting stretching function is continuous with respect to \(\theta _S\) and \(\theta _B\) as their values approach zero. The range of meaningful values for \(\theta _S\) and \(\theta _B\) are:

\[0\leq \theta _S \leq 10 \quad \text{and} \quad 0\leq \theta _B \leq 4\]

However, users need to pay attention to extreme r-factor (rx1) values near the bottom.

REMORA

Code of Conduct

Our Pledge

Our Standards

Our Responsibilities

Scope

Enforcement

Attribution

Getting Started

Downloading the code

Git Submodules

Building

GNU Make

Job info

CMake

Perlmutter (NERSC)

Running

Testing and Verification

Building the Tests

Running the Tests

Adding Tests

Inputs

Problem Geometry

List of Parameters

Examples of Usage

Domain Boundary Conditions

List of Parameters

Resolution

List of Parameters

Examples of Usage

Regridding

Overview

List of Parameters

Notes

Examples of Usage

Gridding and Load Balancing

Simulation Time

List of Parameters

Notes

Examples of Usage

Time Step

List of Parameters for Single-Rate

Examples of Usage

Restart Capability

PlotFiles

Screen Output

List of Parameters

Examples of Usage

Included terms

List of Parameters

Physics Parameters

List of Parameters

Numerical Algorithms

List of Parameters

Vertical Stretch prameters

List of Parameters

Mesh Refinement

Static Mesh Refinement

Dynamic Mesh Refinement

Coupling Types

Checkpoint / Restart

Writing the Checkpoint “Files”

List of Parameters

Restarting

Examples of Usage

Plotfiles

Controlling PlotFile Generation

List of Parameters

Notes

Examples of Usage

Visualization

ParaView

VisIt

yt

Using on a local workstation

Problem Setup

Variable Definitions

Prognostic Equations

Vertical Mixing

Particles