Numerical Solution Technique

Numerical Solution Technique

Horizontal Discretization

<wikitex> In the horizontal, the ROMS governing equations are discretized over a boundary-fitted, orthogonal curvilinear coordinates ($\xi$, $\eta$) grid. The general formulation of the curvilinear coordinates system allows Cartesian, polar and spherical coordinates applications. The transformation of any of these coordinates to ROMS ($\xi$, $\eta$) 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.

ROMS staggered horizontal grid

In ROMS 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$-type	1:Lm(ng),1:Mm(ng)	0:L(ng),0:M(ng)
$\psi$-type	2:Lm(ng),2:Mm(ng)	1:L(ng),1:M(ng)
u-type	2:Lm(ng),1:Mm(ng)	1:L(ng),0:M(ng)
v-type	1:Lm(ng),2:Mm(ng)	0:L(ng),1:M(ng)

</wikitex>

Coastal boundaries can also be specified as a finite-discretized grid via land/sea masking arrays (rmask, umask, and vmask).

Vertical Discretization

<wikitex> The ROMS 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 cell. See diagram below.

Vieste-Dubrovnik Transect

ROMS staggered vertical grid

The total thickness of the water column is $\zeta(i,j)+h(i,j)$. The bathymetry (h) is usually time invariant whereas the free-surface (zeta) evolves in time. However, in sediment applications h changes with time when SED_MORPH is activated. 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. </wikitex>

Masking of Land Areas

<wikitex> ROMS has the ability to work with interior land areas, although the computations occur over the entire model domain. The grid structure as a whole is shown above while two land cells are shown here:

Land-Sea Masking Strategy

The process of defining which areas are to be masked happens during the preparation of the grid file, while this section describes how the masking affects the computation of the various terms in the equations of motion. </wikitex>

Velocity points

<wikitex> At the end of every timestep, the values of many variables within the masked region are set to zero by multiplying by the mask for either the $u$, $v$ or $\rho$ points. This is appropriate for the $v$ points E and L in the figure, since the flow in and out of the land should be zero. It is likewise appropriate for the $u$ point at I, but is not necessarily correct for point G. The only term in the $u$ equation that requires the $u$ value at point G is the horizontal viscosity, which has a term of the form ${\partial \over \partial \eta} \nu {\partial u \over \partial \eta}$. Since point G is used in this term by both points A and M, it is not sufficient to replace its value with that of the image point for A. Instead, the term ${\partial u \over \partial \eta}$ is computed and the values at points D and K are replaced with the values appropriate for either free-slip or no-slip boundary conditions. Likewise, the term ${\partial \over \partial \xi} \nu {\partial v \over \partial \xi}$ in the $v$ equation must be corrected at the mask boundaries.

This is accomplished by having a fourth mask array defined at the $\psi$ points, in which the values are set to be no-slip in metrics.F. For no-slip boundaries, we count on the values inside the land (point G) having been zeroed out. For point D, the image point at G should contain minus the value of $u$ at point A. The desired value of ${\partial u \over \partial \eta}$ is therefore $2 u_{\bf A}$ while instead we have simply $u_{\bf A}$. In order to achieve the correct result, we multiply by a mask which contains the value 2 at point D. It also contains a 2 at point K so that ${\partial u \over \partial \eta}$ there will acquire the desired value of $-2 u_{\bf M}$. The corner point F is set to have a value of 1. </wikitex>

Temperature, Salinity and Surface Elevation

<wikitex> The handling of masks by the temperature, salinity and surface elevation equations is similar to that in the momentum equations, and is in fact simpler. Values of $T$, $S$ and $\zeta$ inside the land masks, such as point H in the figure above, are set to zero after every timestep. This point would be used by the horizontal diffusion term for points B, J, and N. This is handled by setting the first derivative terms at points E, I, and L to zero, to be consistent with a no-flux boundary condition. Note that the equation of state must be able to handle $T = S = 0$ since this is the value inside masked regions.

Wetting and Drying

There is now an option to have wetting and drying in the model, in which a cell can switch between being wet or being dry as the tides come in and go out, for instance. Cells which are masked out as above are never allowed to be wet, however.

• In the case of wetting and drying, a critical depth, $D_{crit}$, is supplied by the user.
• The total water depth ($D=h+\zeta$) is compared to $D_{crit}$. If the water level is less than this depth, no flux is allowed out of that cell. Water can always flow in and resubmerge the cell.
• The wetting and drying only happens during the 2-D computations; the 3-D computations see a depth of $D_{crit}$ in the "dry" areas.

</wikitex>

Nesting Grids

<wikitex> ROMS has a nested design structure for the temporal coupling between a hierarchy of parent and child grids. All the state model variables are dynamically allocated and passed as arguments to the computational routines via Fortran-90 dereferenced pointer structures. All private arrays are automatic; their size is determined when the procedure is entered. This code structure facilitates computations over nested grids of different topologies. There are three types of nesting capabilities in ROMS:

• Refinement Grids which provide increased resolution (3:1 or 5:1) in a specific region,
• Mosaics Grids which connect several grids along their edges, and
• Composite Grids which allow overlap regions of aligned and non-aligned grids.

</wikitex>

Refinement Grids

<wikitex>

Blue and Red Grids

Red Grid Contact Points

Blue Grid Contact Points

</wikitex>

Mosaic Grids

<wikitex>

Blue and Red Grids

Red Grid Contact Points

Blue Grid Contact Points

</wikitex>

Composite Grids

<wikitex>

Red and Blue Grids

Red Grid Contact Points

Blue Grid Contact Points

</wikitex>

Time-Stepping Overview

<wikitex> 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:

1. Take a predictor step for at least the 3-D tracers to time $n+1/2$.
2. Compute $\overline{\rho}$ and $\rho^*$ for use in the depth-integrated time-steps, from the density either at time $n$ or time $n+1/2$.
3. Depth integrate the 3-D momentum right-hand side terms at time $n+1/2$ for use in the depth-integrated time-steps (or extrapolate to obtain an estimate of those terms).
4. Take all the depth-integrated steps. Store weighted time-means of the $\overline{u}$, $\overline{v}$ fields centered at both time $n+1/2$ and time $n+1$ (plus $\zeta$ at time $n+1$). The latter requires this time-stepping to extend past time $n+1$, using $M^*$ steps rather than just $M$.
5. Use the weighted time-means from depth-integrated fields to complete the corrector step for the 3-D fields to time $n+1$.

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 McWilliams (2005) and Shchepetkin and McWilliams (2008b); the relevant ones are described in Time-stepping Schemes.

The following figures show the time-stepping and mode coupling used in the various ROMS versions. In all, 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:

ROMS AGRIF:

UCLA ROMS:

Non-hydrostatic ROMS:

Timestepping table from Shchepetkin and McWilliams (2008a):

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 number, the first one is for momenta, the second for tracers. </wikitex>

Conservation Properties

<wikitex> From Shchepetkin and McWilliams (2005), we have a tracer concentration equation in advective form: $$\frac{\partial C}{\partial t} + (u \cdot \nabla) C = 0 \eqno{(1)}$$ and also a tracer concentration equation in conservation form: $$\frac{\partial C}{\partial t} + \nabla \cdot (u C) = 0. \eqno{(2)}$$ The continuity equation: $$( \nabla \cdot u) = 0 \eqno{(3)}$$ 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: $$ \eqalign{ {\partial \over \partial t} \left( {H_z C \over m n} \right) &+ \delta_{\xi} \left( {u \overline{H_z}^\xi \overline{C}^\xi \over \overline{n}^\xi} \right) + \delta_{\eta} \left( {v \overline{H_z}^\eta \overline{C}^\eta \over \overline{m}^\eta} \right) + \delta_\sigma \left( \overline{C}^\sigma {H_z \Omega \over m n} \right) =\cr &{ 1 \over mn} {\partial \over \partial \sigma} \left( {K_m \over \Delta z} {\partial C \over \partial \sigma} \right) + {\cal D}_C + {\cal F}_C} \eqno{(4)}$$

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{ ( \qquad )}^{\xi}$, $\overline{ ( \qquad )}^{\eta}$ and $\overline{ ( \qquad)}^\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(\frac{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$, with 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}{m n} \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}{m n} \right) = 0 \eqno{(5)}$$

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 ROMS, the discrete continuity equation is needed to compute vertical velocity rather than grid-box volume $\frac{H_z}{m n}$ (the latter is controlled by changes in $\zeta$ in the barotropic mode computations). Here, $\frac{H_z \Omega}{m n}$ 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 \eqno{(6)}$$

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.

</wikitex>

Depth-Integrated Equations

<wikitex>The depth average of a quantity $A$ is given by: $$\overline{A} = \frac{1}{D} \int_{-1}^0 H_z A d\sigma \eqno{(7)} $$ where the overbar indicates a vertically averaged quantity and $$ D \equiv \zeta(\xi, \eta, t) + h(\xi, \eta) \eqno{(8)} $$ is the total depth of the water column. The vertical integral of the momentum equations are: $$ \eqalign{ \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} \cr
 - \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) \cr
   + \frac{ D}{mn}
  \left( \overline{\cal F}_u + \overline{\cal D}_{h_u} \right) 
  &+ \frac{1}{mn} \left( \tau^{\xi}_s - \tau^{\xi}_b \right) } \eqno{(9)} $$

and $$ \eqalign{

   \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} \cr
     + \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) \cr
 + \frac{ D}{mn}
  \left( \overline{\cal F}_v + \overline{\cal D}_{h_v} \right) 
  &+ \frac{1}{mn} \left( \tau^{\eta}_s - \tau^{\eta}_b \right)} \eqno{(10)} $$

where $\phi_2$ includes the $\frac{\partial z}{\partial \xi}$ term, $\overline{\cal 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_{\rm slow}}$ and $R_{v_{\rm slow}}$, equations (9) and (10) become: $$ \eqalign{

  \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} \cr
  - \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_{\rm slow}} -
  \frac{gD}{n} \frac{\partial \zeta}{\partial \xi}
   + \frac{D}{mn} {\cal D}_{\overline{u}}
  - \frac{1}{mn} \tau^{\xi}_b} \eqno{(11)} $$

and $$ \eqalign{

  \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} \cr
  + \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_{\rm slow}} -
  \frac{gD}{m} \frac{\partial \zeta}{\partial \eta}
   + \frac{D}{mn} {\cal D}_{\overline{v}}
  - \frac{1}{mn} \tau^{\eta}_b .} \eqno{(12)} $$

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_{\rm slow}}$ and $R_{v_{\rm 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 $(n+$ dtfast $\times M^\star)$ 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 vertical 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 \sum_{m=1}^{M^\star} 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 \sum_{m=1}^{M^\star} b_m \overline{u}^m$. In order to maintain constancy preservation, this relation must hold: $$ \eqalign{ \langle \zeta \rangle_{i,j}^{n+1} = \langle \zeta \rangle_{i,j}^n &- \cr

 (mn)_{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]} \eqno{(13)} $$

Shchepetkin and McWilliams (2005) introduce a range of possible weights, but the ones used here have a shape function: $$ A(\tau) = 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\} \eqno{(14)} $$

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^\star} \tau^n A(\tau) d \tau = 1, \quad n=0,1,2 \eqno{(15)} $$ using Newton iterations. $\tau^\star$ is the upper limit of $\tau$ with $A(\tau) \geq 0$. In practice we initially set $$ A_0 = 1, r = 0, \quad \hbox{and} \quad \tau = \frac{(p+2)(p+q+2)}{(p+1)(p+q+1)} $$ compute $A(\tau)$ using eq.~(14), normalize using: $$ \sum_{m=1}^{M^\star} a_m \equiv 1, \quad

  \sum_{m=1}^{M^\star} a_m\frac{m}{M} \equiv 1,  \eqno{(16)} $$

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^\star$ to be less than $1.5M$.

ROMS also supports an older cosine weighting option, which isn't recommended since it is only first-order accurate. </wikitex>

Pressure Gradient Terms in Mode Coupling

<wikitex>Equation (11) contains the term $R_{u_{\rm 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{g D}{n} \frac{\partial \zeta}{\partial \xi} +

  \left[\frac{g D}{n} \frac{\partial \zeta}{\partial \xi} + {\cal F}
  \right] \eqno{(17)} $$

where the term in square brackets is the mode coupling term and is held fixed over all the barotropic steps and $$ {\cal F} = - \frac{1}{\rho_0 n} \int_{-h}^\zeta \frac{\partial

 P}{\partial \xi} dz \eqno{(18)} $$

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 , \quad

 \rho^\star = \frac{1}{\frac{1}{2} D^2} \int_{-h}^\zeta
 \left\{ \int_z^\zeta \rho dz^{\prime} \right\} dz \eqno{(19)} $$

Changing the vertical coordinate to $\sigma$ yields: $$ \overline{\rho} = \int_{-1}^0 \rho d\sigma , \quad

 \rho^\star = 2 \int_{-1}^0
 \left\{ \int_\sigma^0 \rho d\sigma^{\prime} \right\} d\sigma \eqno{(20)} $$

which implies that $\overline{\rho}$ and $\rho^\star$ are actually independent of $\zeta$ as long as the density profile $\rho = \rho(\sigma)$ does not change. The vertically integrated pressure gradient becomes: $$ -\frac{1}{\rho_0} \frac{g}{n} \left\{ \frac{\partial}{\partial \xi}

  \left( \frac{\rho^\star D^2}{2} \right) - \overline{\rho} D
  \frac{\partial h}{\partial \xi} \right\} = -\frac{1}{\rho_0}
  \frac{g}{n} D \left\{ \rho^\star \frac{\partial \zeta}{\partial \xi} +
  \frac{D}{2} \frac{\partial \rho^\star}{\partial \xi} + (\rho^\star -
  \overline{\rho}) \frac{\partial h}{\partial \xi} \right\} \eqno{(21)} $$

In the case of uniform density $\rho_0$, we obtain $\rho^\star \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). </wikitex>

Time Stepping: Internal Velocity Modes and Tracers

<wikitex>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 (see Time-stepping Schemes). These stored up r.h.s. values can be used to extrapolate to a value at time $n+\frac{1}{2}$ for use in the barotropic steps as shown in #Time-Stepping Overview.

The tracer concentration equation is advanced in a predictor-corrector leapfrog-trapezoidal step, with great care taken to optimize both the conservation and constancy-preserving properties of the continuous equations. The corrector step can maintain both, as long as it uses velocities and column depths which satisfy eq. (13). This also requires tracer values centered 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. </wikitex>

Horizontal and Vertical Advection

<wikitex>The advection of a tracer $C$ has an equation of the form $$

 \frac{\partial}{\partial t} \frac{ H_z C}{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: $$ \eqalign{

  F^\xi &= \frac{H_z u C}{n} \cr
  F^\eta &= \frac{H_z v C}{m} \cr
  F^\sigma &= \frac{H_z \Omega C}{mn} .}

$$ </wikitex>

Second-order Centered

<wikitex>The simplest form of the advective fluxes is the centered second-order: $$ \eqalign{

  F^\xi &= \frac{\overline{H_z}^\xi u
              \overline{C}^\xi}{\overline{n}^\xi} \cr
  F^\eta &= \frac{\overline{H_z}^\eta v
              \overline{C}^\eta}{\overline{m}^\eta} \cr
  F^\sigma &= \frac{\overline{H_z}^\sigma \Omega
              \overline{C}^\sigma}{mn} .}

$$ This scheme is known to have some unfortunate properties in the presence of strong gradients, such as large over- and under-shoots of tracers, leading to the need for large amounts of horizontal smoothing. ROMS provides alternative advection schemes with better behavior in many situations, but retains this one for comparison purposes.</wikitex>

Fourth-order Centered

<wikitex>The barotropic advection is centered fourth-order unless you specifically pick centered second-order as your horizontal advection scheme. To get fourth-order, create gradient terms: $$ \eqalign{

    G^{\xi} &= \overline{\left(\frac{\partial C}{\partial
           \xi}\right)}^\xi \cr
    G^{\eta} &= \overline{\left(\frac{\partial C}{\partial
           \eta}\right)}^\eta \cr
    G^{\sigma} &= \overline{\left(\frac{\partial C}{\partial
           \sigma}\right)}^\sigma . }

$$ The fluxes now become: $$ \eqalign{

  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)  \cr
  F^\eta &= \frac{\overline{H_z}}{\overline{m}^\eta} ^\eta v \left(
              \overline{C}^\eta-
          \frac{1}{3} \frac{\partial
              G^{\eta}}{\partial \eta} \right) \cr
  F^\sigma &= \frac{\overline{H_z}^\sigma} {mn} \Omega \left(
              \overline{C}^\sigma -
          \frac{1}{3} \frac{\partial
              G^{\sigma}}{\partial \sigma} \right). } \eqno{(22)} 

$$ </wikitex>

Fourth-order Akima

<wikitex>An alternate fourth-order algorithm is that by Akima: $$ \eqalign{

    G_{\xi} &= 2 {\frac{\partial C}{\partial \xi}_i
        \frac{\partial C}{\partial \xi}_{i+1}} \, \bigg/ \left(
        \frac{\partial C}{\partial \xi}_i +
         \frac{\partial C}{\partial \xi}_{i+1} \right) \cr
    G_{\eta} &= 2 {\frac{\partial C}{\partial \eta}_j
        \frac{\partial C}{\partial \eta}_{j+1}} \, \bigg/ \left(
        \frac{\partial C}{\partial \eta}_j +
         \frac{\partial C}{\partial \eta}_{j+1} \right) \cr
    G_{\sigma} &= 2 {\frac{\partial C}{\partial \sigma}_k
        \frac{\partial C}{\partial \sigma}_{k-1}} \, \bigg/ \left(
        \frac{\partial C}{\partial \sigma}_k +
         \frac{\partial C}{\partial \sigma}_{k-1} \right) }

$$ With the fluxes as in eq. (22). </wikitex>

Third-order Upwind

<wikitex>There is a class of third-order upwind advection schemes, both one-dimensional (Leonard \cite{Leonard79}) and two-dimensional (Rasch \cite{Rasch94} and Shchepetkin and McWilliams \cite{SS98}). This scheme is known as UTOPIA (Uniformly Third-Order Polynomial Interpolation Algorithm). Applying flux limiters to UTOPIA is explored in Thuburn \cite{Thuburn96}, although it is not implemented in ROMS. 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 ROMS, leaving one extra upwind term in the computation of the advective fluxes: $$ \eqalign{

  F^\xi &= \frac{H_z u}{n} \left( C - \gamma \frac{\partial^2
  C}{\partial \xi^2} \right) \cr
  F^\eta &= \frac{H_z v}{m} \left( C - \gamma \frac{\partial^2
  C}{\partial \eta^2} \right)}

$$ 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}^\xi}{\overline{n}^\xi}
  \left[ \max(0,u_{i,j,k}) C_{i-1,j,k} +
  \min(0,u_{i,j,k}) C_{i,j,k} \right] .  \eqno{(23)}

$$ The value of $\gamma$ in the model is $\frac{1}{8}$ while that in Rasch \cite{Rasch94} 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_z w}{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: $$ \eqalign{

  F^\xi &= \left(u - \gamma \frac{\partial^2 u}{\partial \xi^2} \right)
  \left[ \frac{H_z u}{n} - \gamma \frac{\partial^2}{\partial \xi^2}
  \left( \frac{H_z u}{n} \right) \right] \cr
  F^\eta &= \left(u - \gamma \frac{\partial^2 u}{\partial \eta^2}
    \right)
  \left[ \frac{H_z v}{m} - \gamma \frac{\partial^2}{\partial \xi^2}
  \left( \frac{H_z v}{m} \right) \right] \cr
  F^\sigma &= \frac{H_z w}{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] }

$$ while for the $v$-velocity we have: $$ \eqalign{

  F^\xi &= \left(v - \gamma \frac{\partial^2 v}{\partial \xi^2} \right)
  \left[ \frac{H_z u}{n} - \gamma \frac{\partial^2}{\partial \eta^2}
  \left( \frac{H_z u}{n} \right) \right] \cr
  F^\eta &= \left(v - \gamma \frac{\partial^2 v}{\partial \eta^2}
    \right)
  \left[ \frac{H_z v}{m} - \gamma \frac{\partial^2}{\partial \eta^2}
  \left( \frac{H_z v}{m} \right) \right] \cr
  F^\sigma &= \frac{H_z w}{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]  }

$$ In all these terms, the second derivatives are evaluated at an upstream location. </wikitex>

MPDATA

<wikitex></wikitex>

Vertical Velocity

<wikitex>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 combining equations (\ref{st17}) and (\ref{zeta1}) to obtain: $$

  \frac{\partial}{\partial \xi} \left( \frac{H_z u}{n} \right) +
  \frac{\partial}{\partial \eta} \left( \frac{H_z v}{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?] </wikitex>

Vertical Friction and Diffusion

<wikitex></wikitex>

	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)