# Numerical Solution Technique

Numerical Solution Technique

## Horizontal Discretization

In the horizontal, the ROMS governing equations are discretized over a boundary-fitted, orthogonal curvilinear coordinates (, ) 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 (, ) 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.

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

Variable Interior Range Full Range
-type 1:Lm(ng),1:Mm(ng) 0:L(ng),0:M(ng)
-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)

## Vertical Discretization

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.

The total thickness of the water column is . 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.

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:

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.

### Velocity points

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 , or points. This is appropriate for the 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 point at I, but is not necessarily correct for point G. The only term in the equation that requires the value at point G is the horizontal viscosity, which has a term of the form . 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 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 in the equation must be corrected at the mask boundaries.

This is accomplished by having a fourth mask array defined at the 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 at point A. The desired value of is therefore while instead we have simply . 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 there will acquire the desired value of . The corner point F is set to have a value of 1.

### Temperature, Salinity and Surface Elevation

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 , and 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 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, , is supplied by the user.
• The total water depth () is compared to . 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 in the "dry" areas.

## Time-Stepping Overview

While time-stepping the model, we have a stored history of the model fields at time , an estimate of the fields at the current time , and we need to come up with an estimate for time . 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 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 .
2. Compute and for use in the depth-integrated time-steps, from the density either at time or time .
3. Depth integrate the 3-D momentum right-hand side terms at time 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 , fields centered at both time and time (plus at time ). The latter requires this time-stepping to extend past time , using steps rather than just .
5. Use the weighted time-means from depth-integrated fields to complete the corrector step for the 3-D fields to time .

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 (2009); the relevant ones are described in Time-stepping Schemes Review.

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 , while the light red function has weights to produce an estimate at time .

Non-hydrostatic 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

Gen. FB

(AB3-AM4)

Gen. FB

(AB3-AM4)

2-D , iter. , (2) 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)

, advect. 0.72 0.72 1.587 1.587 0.78
, Cor. 0.72 0.72 1.587 1.587 0.78
, int. w. 0.72, (1) 1.14, (1,2) 1.85, (2) 1.85, (2) 1.78, (1) Note: The generalized FB barotropic mode was ported into the newest AGRIF code at the end of 2007.

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.

## Conservation Properties

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

 $StartFraction normal partial-differential upper C Over normal partial-differential t EndFraction plus left-parenthesis u dot normal nabla right-parenthesis upper C equals 0$ (1)

and also a tracer concentration equation in conservation form:

 $StartFraction normal partial-differential upper C Over normal partial-differential t EndFraction plus normal nabla dot left-parenthesis u upper C right-parenthesis equals 0 period$ (2)

The continuity equation:

 $left-parenthesis normal nabla dot u right-parenthesis equals 0$ (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:

 $StartLayout 1st Row 1st Column StartFraction normal partial-differential Over normal partial-differential t EndFraction left-parenthesis StartFraction upper H Subscript z Baseline upper C Over m n EndFraction right-parenthesis 2nd Column plus delta Subscript xi Baseline left-parenthesis StartFraction u upper H Subscript z Baseline overbar Superscript xi Baseline upper C overbar Superscript xi Baseline Over n overbar Superscript xi Baseline EndFraction right-parenthesis plus delta Subscript eta Baseline left-parenthesis StartFraction v upper H Subscript z Baseline overbar Superscript eta Baseline upper C overbar Superscript eta Baseline Over m overbar Superscript eta Baseline EndFraction right-parenthesis plus delta Subscript sigma Baseline left-parenthesis upper C overbar Superscript sigma Baseline StartFraction upper H Subscript z Baseline normal upper Omega Over m n EndFraction right-parenthesis equals 2nd Row 1st Column Blank 2nd Column StartFraction 1 Over m n EndFraction StartFraction normal partial-differential Over normal partial-differential sigma EndFraction left-parenthesis StartFraction upper K Subscript m Baseline Over normal upper Delta z EndFraction StartFraction normal partial-differential upper C Over normal partial-differential sigma EndFraction right-parenthesis plus script upper D Subscript upper C Baseline plus script upper F Subscript upper C EndLayout$ (4)

Here , and denote simple centered finite-difference approximations to , and with the differences taken over the distances , and , respectively. is the vertical distance from one point to another. , and represent averages taken over the distances , and .

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

The quantity 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 to time , with the other terms being evaluated at time for second-order accuracy. Setting to 1 everywhere reduces eq. (4) to:

 $StartFraction normal partial-differential Over normal partial-differential t EndFraction left-parenthesis StartFraction upper H Subscript z Baseline Over m n EndFraction right-parenthesis plus delta Subscript xi Baseline left-parenthesis StartFraction u upper H Subscript z Baseline overbar Superscript xi Baseline Over n overbar Superscript xi Baseline EndFraction right-parenthesis plus delta Subscript eta Baseline left-parenthesis StartFraction v upper H Subscript z Baseline overbar Superscript eta Baseline Over m overbar Superscript eta Baseline EndFraction right-parenthesis plus delta Subscript sigma Baseline left-parenthesis StartFraction upper H Subscript z Baseline normal upper Omega Over m n EndFraction right-parenthesis equals 0$ (5)

If this equation holds true for the step from time to time , 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 (the latter is controlled by changes in in the barotropic mode computations). Here, is the finite-volume flux across the moving grid-box interface, vertically on the grid.

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

 $StartFraction normal partial-differential Over normal partial-differential t EndFraction left-parenthesis StartFraction zeta Over m n EndFraction right-parenthesis plus delta Subscript xi Baseline left-parenthesis StartFraction u overbar upper D overbar Superscript xi Baseline Over n overbar Superscript xi Baseline EndFraction right-parenthesis plus delta Subscript eta Baseline left-parenthesis StartFraction v overbar upper D overbar Superscript eta Baseline Over m overbar Superscript eta Baseline EndFraction right-parenthesis equals 0$ (6)

where is the surface elevation, is the total depth, and 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 is given by:

 $upper A overbar equals StartFraction 1 Over upper D EndFraction integral Subscript negative 1 Superscript 0 Baseline upper H Subscript z Baseline upper A d sigma$ (7)

where the overbar indicates a vertically averaged quantity and

 $upper D identical-to zeta left-parenthesis xi comma eta comma t right-parenthesis plus h left-parenthesis xi comma eta right-parenthesis$ (8)

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

 $StartLayout 1st Row 1st Column StartFraction normal partial-differential Over normal partial-differential t EndFraction left-parenthesis StartFraction upper D u overbar Over m n EndFraction right-parenthesis plus StartFraction normal partial-differential Over normal partial-differential xi EndFraction left-parenthesis StartFraction upper D ModifyingAbove u u With bar Over n EndFraction right-parenthesis plus StartFraction normal partial-differential Over normal partial-differential eta EndFraction left-parenthesis StartFraction upper D ModifyingAbove u v With bar Over m EndFraction right-parenthesis minus 2nd Column StartFraction upper D f v overbar Over m n EndFraction 2nd Row 1st Column minus left-bracket ModifyingAbove v v With bar StartFraction normal partial-differential Over normal partial-differential xi EndFraction left-parenthesis StartFraction 1 Over n EndFraction right-parenthesis minus ModifyingAbove u v With bar StartFraction normal partial-differential Over normal partial-differential eta EndFraction left-parenthesis StartFraction 1 Over m EndFraction right-parenthesis right-bracket upper D equals minus StartFraction upper D Over n EndFraction 2nd Column left-parenthesis StartFraction normal partial-differential ModifyingAbove phi 2 With bar Over normal partial-differential xi EndFraction plus g StartFraction normal partial-differential zeta Over normal partial-differential xi EndFraction right-parenthesis 3rd Row 1st Column plus StartFraction upper D Over m n EndFraction left-parenthesis script upper F overbar Subscript u Baseline plus script upper D overbar Subscript h Sub Subscript u Subscript Baseline right-parenthesis 2nd Column plus StartFraction 1 Over m n EndFraction left-parenthesis tau Subscript s Superscript xi Baseline minus tau Subscript b Superscript xi Baseline right-parenthesis EndLayout$ (9)

and

 $StartLayout 1st Row 1st Column StartFraction normal partial-differential Over normal partial-differential t EndFraction left-parenthesis StartFraction upper D v overbar Over m n EndFraction right-parenthesis plus StartFraction normal partial-differential Over normal partial-differential xi EndFraction left-parenthesis StartFraction upper D ModifyingAbove u v With bar Over n EndFraction right-parenthesis plus StartFraction normal partial-differential Over normal partial-differential eta EndFraction left-parenthesis StartFraction upper D ModifyingAbove v v With bar Over m EndFraction right-parenthesis plus 2nd Column StartFraction upper D f u overbar Over m n EndFraction 2nd Row 1st Column plus left-bracket ModifyingAbove u v With bar StartFraction normal partial-differential Over normal partial-differential xi EndFraction left-parenthesis StartFraction 1 Over n EndFraction right-parenthesis minus ModifyingAbove u u With bar StartFraction normal partial-differential Over normal partial-differential eta EndFraction left-parenthesis StartFraction 1 Over m EndFraction right-parenthesis right-bracket upper D equals minus StartFraction upper D Over m EndFraction 2nd Column left-parenthesis StartFraction normal partial-differential ModifyingAbove phi 2 With bar Over normal partial-differential eta EndFraction plus g StartFraction normal partial-differential zeta Over normal partial-differential eta EndFraction right-parenthesis 3rd Row 1st Column plus StartFraction upper D Over m n EndFraction left-parenthesis script upper F overbar Subscript v Baseline plus script upper D overbar Subscript h Sub Subscript v Subscript Baseline right-parenthesis 2nd Column plus StartFraction 1 Over m n EndFraction left-parenthesis tau Subscript s Superscript eta Baseline minus tau Subscript b Superscript eta Baseline right-parenthesis EndLayout$ (10)

where includes the term, 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 . 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 and 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 and , equations (9) and (10) become:

 $StartLayout 1st Row 1st Column StartFraction normal partial-differential Over normal partial-differential t EndFraction left-parenthesis StartFraction upper D u overbar Over m n EndFraction right-parenthesis plus StartFraction normal partial-differential Over normal partial-differential xi EndFraction left-parenthesis StartFraction upper D u overbar u overbar Over n EndFraction right-parenthesis plus StartFraction normal partial-differential Over normal partial-differential eta EndFraction left-parenthesis StartFraction upper D u overbar v overbar Over m EndFraction right-parenthesis 2nd Column minus StartFraction upper D f v overbar Over m n EndFraction 2nd Row 1st Column minus left-bracket v overbar v overbar StartFraction normal partial-differential Over normal partial-differential xi EndFraction left-parenthesis StartFraction 1 Over n EndFraction right-parenthesis minus u overbar v overbar StartFraction normal partial-differential Over normal partial-differential eta EndFraction left-parenthesis StartFraction 1 Over m EndFraction right-parenthesis right-bracket upper D 2nd Column equals upper R Subscript u Sub Subscript normal s normal l normal o normal w Baseline minus StartFraction g upper D Over n EndFraction StartFraction normal partial-differential zeta Over normal partial-differential xi EndFraction plus StartFraction upper D Over m n EndFraction script upper D Subscript u overbar Baseline minus StartFraction 1 Over m n EndFraction tau Subscript b Superscript xi EndLayout$ (11)

and

 $StartLayout 1st Row 1st Column StartFraction normal partial-differential Over normal partial-differential t EndFraction left-parenthesis StartFraction upper D v overbar Over m n EndFraction right-parenthesis plus StartFraction normal partial-differential Over normal partial-differential xi EndFraction left-parenthesis StartFraction upper D u overbar v overbar Over n EndFraction right-parenthesis plus StartFraction normal partial-differential Over normal partial-differential eta EndFraction left-parenthesis StartFraction upper D v overbar v overbar Over m EndFraction right-parenthesis 2nd Column plus StartFraction upper D f u overbar Over m n EndFraction 2nd Row 1st Column plus left-bracket u overbar v overbar StartFraction normal partial-differential Over normal partial-differential xi EndFraction left-parenthesis StartFraction 1 Over n EndFraction right-parenthesis minus u overbar u overbar StartFraction normal partial-differential Over normal partial-differential eta EndFraction left-parenthesis StartFraction 1 Over m EndFraction right-parenthesis right-bracket upper D 2nd Column equals upper R Subscript v Sub Subscript normal s normal l normal o normal w Subscript Baseline minus StartFraction g upper D Over m EndFraction StartFraction normal partial-differential zeta Over normal partial-differential eta EndFraction plus StartFraction upper D Over m n EndFraction script upper D Subscript v overbar Baseline minus StartFraction 1 Over m n EndFraction tau Subscript b Superscript eta Baseline period EndLayout$ (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 and . This was found to be the easiest way to retain the baroclinic contributions of the non-linear terms such as .

The model is time stepped from time to time by using short time steps on equations (11), (12) and (6). Equation (6) is time stepped first, so that an estimate of the new 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 dtfast and while maintaining weighted averages of the values of , and . The averages are used to replace the values at time in both the baroclinic and barotropic modes, and for recomputing the vertical grid spacing . The following figure shows one option for how these weights might look: The primary weights, , are used to compute . There is a related set of secondary weights , used as . In order to maintain constancy preservation, this relation must hold:

 $StartLayout 1st Row 1st Column mathematical left-angle zeta mathematical right-angle Subscript i comma j Superscript n plus 1 Baseline equals mathematical left-angle zeta mathematical right-angle Subscript i comma j Superscript n Baseline 2nd Column minus 2nd Row 1st Column left-parenthesis m n right-parenthesis Subscript i comma j Baseline normal upper Delta t 2nd Column left-bracket mathematical left-angle mathematical left-angle StartFraction upper D u overbar Over n EndFraction mathematical right-angle mathematical right-angle Subscript i plus one-half comma j Superscript n plus one-half Baseline minus mathematical left-angle mathematical left-angle StartFraction upper D u overbar Over n EndFraction mathematical right-angle mathematical right-angle Subscript i minus one-half comma j Superscript n plus one-half Baseline plus mathematical left-angle mathematical left-angle StartFraction upper D v overbar Over m EndFraction mathematical right-angle mathematical right-angle Subscript i comma j plus one-half Superscript n plus one-half Baseline minus mathematical left-angle mathematical left-angle StartFraction upper D v overbar Over m EndFraction mathematical right-angle mathematical right-angle Subscript i comma j minus one-half Superscript n plus one-half Baseline right-bracket EndLayout$ (13)

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

 $upper A left-parenthesis tau right-parenthesis equals upper A 0 StartSet left-parenthesis StartFraction tau Over tau 0 EndFraction right-parenthesis Superscript p Baseline left-bracket 1 minus left-parenthesis StartFraction tau Over tau 0 EndFraction right-parenthesis Superscript q Baseline right-bracket minus r StartFraction tau Over tau 0 EndFraction EndSet$ (14)

where are parameters and , and are chosen to satisfy normalization, consistency, and second-order accuracy conditions,

 $upper I Subscript n Baseline equals integral Subscript 0 Superscript tau Superscript star Baseline Baseline tau Superscript n Baseline upper A left-parenthesis tau right-parenthesis d tau equals 1 comma n equals 0 comma 1 comma 2$ (15)

using Newton iterations. is the upper limit of with . In practice we initially set

compute using eq.~(14), normalize using:

 $sigma-summation Underscript m equals 1 Overscript upper M Superscript star Baseline Endscripts a Subscript m Baseline identical-to 1 comma sigma-summation Underscript m equals 1 Overscript upper M Superscript star Baseline Endscripts a Subscript m Baseline StartFraction m Over upper M EndFraction identical-to 1 comma$ (16)

and adjust iteratively to satisfy the condition of (15). We are using values of , , and . This form allows some negative weights for small , allowing to be less than .

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

## Pressure Gradient Terms in Mode Coupling

Equation (11) contains the term , 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:

 $minus StartFraction g upper D Over n EndFraction StartFraction normal partial-differential zeta Over normal partial-differential xi EndFraction plus left-bracket StartFraction g upper D Over n EndFraction StartFraction normal partial-differential zeta Over normal partial-differential xi EndFraction plus script upper F right-bracket$ (17)

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

 $script upper F equals minus StartFraction 1 Over rho 0 n EndFraction integral Subscript negative h Superscript zeta Baseline StartFraction normal partial-differential upper P Over normal partial-differential xi EndFraction d z$ (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:

 $rho overbar equals StartFraction 1 Over upper D EndFraction integral Subscript negative h Superscript zeta Baseline rho d z comma rho Superscript star Baseline equals StartFraction 1 Over one-half upper D squared EndFraction integral Subscript negative h Superscript zeta Baseline StartSet integral Subscript z Superscript zeta Baseline rho d z Superscript prime Baseline EndSet d z$ (19)

Changing the vertical coordinate to yields:

 $rho overbar equals integral Subscript negative 1 Superscript 0 Baseline rho d sigma comma rho Superscript star Baseline equals 2 integral Subscript negative 1 Superscript 0 Baseline StartSet integral Subscript sigma Superscript 0 Baseline rho d sigma Superscript prime Baseline EndSet d sigma$ (20)

which implies that and are actually independent of as long as the density profile does not change. The vertically integrated pressure gradient becomes:

 $minus StartFraction 1 Over rho 0 EndFraction StartFraction g Over n EndFraction StartSet StartFraction normal partial-differential Over normal partial-differential xi EndFraction left-parenthesis StartFraction rho Superscript star Baseline upper D squared Over 2 EndFraction right-parenthesis minus rho overbar upper D StartFraction normal partial-differential h Over normal partial-differential xi EndFraction EndSet equals minus StartFraction 1 Over rho 0 EndFraction StartFraction g Over n EndFraction upper D StartSet rho Superscript star Baseline StartFraction normal partial-differential zeta Over normal partial-differential xi EndFraction plus StartFraction upper D Over 2 EndFraction StartFraction normal partial-differential rho Superscript star Baseline Over normal partial-differential xi EndFraction plus left-parenthesis rho Superscript star Baseline minus rho overbar right-parenthesis StartFraction normal partial-differential h Over normal partial-differential xi EndFraction EndSet$ (21)

In the case of uniform density , we obtain , 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).

## 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 and . The depth-averaged component is then removed and replaced by the and 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 Review). These stored up r.h.s. values can be used to extrapolate to a value at time 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 , 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.

The advection of a tracer has an equation of the form

where we have introduced the advective fluxes:

### Second-order Centered

The simplest form of the advective fluxes is the centered second-order:

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.

### Fourth-order Centered

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:

The fluxes now become:

 $StartLayout 1st Row 1st Column upper F Superscript xi 2nd Column equals StartFraction upper H Subscript z Baseline overbar Superscript xi Baseline Over n overbar Superscript xi Baseline EndFraction u left-parenthesis upper C overbar Superscript xi Baseline minus one-third StartFraction normal partial-differential upper G Superscript xi Baseline Over normal partial-differential xi EndFraction right-parenthesis 2nd Row 1st Column upper F Superscript eta 2nd Column equals StartFraction upper H Subscript z Baseline overbar Over m overbar Superscript eta Baseline EndFraction Superscript eta Baseline v left-parenthesis upper C overbar Superscript eta Baseline minus one-third StartFraction normal partial-differential upper G Superscript eta Baseline Over normal partial-differential eta EndFraction right-parenthesis 3rd Row 1st Column upper F Superscript sigma 2nd Column equals StartFraction upper H Subscript z Baseline overbar Superscript sigma Baseline Over m n EndFraction normal upper Omega left-parenthesis upper C overbar Superscript sigma Baseline minus one-third StartFraction normal partial-differential upper G Superscript sigma Baseline Over normal partial-differential sigma EndFraction right-parenthesis period EndLayout$ (22)

### Fourth-order Akima

An alternate fourth-order algorithm is that by Akima:

With the fluxes as in eq. (22).

### Third-order Upwind

There is a class of third-order upwind advection schemes, both one-dimensional (Leonard, 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 ROMS. The two-dimensional formulation in Rasch contains terms of order and , including cross terms (). The terms which are nonlinear in velocity have been dropped in ROMS, leaving one extra upwind term in the computation of the advective fluxes:

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

 $upper F Subscript i comma j comma k Superscript xi Baseline equals StartFraction upper H Subscript z Baseline overbar Superscript xi Baseline Over n overbar Superscript xi Baseline EndFraction left-bracket max left-parenthesis 0 comma u Subscript i comma j comma k Baseline right-parenthesis upper C Subscript i minus 1 comma j comma k Baseline plus min left-parenthesis 0 comma u Subscript i comma j comma k Baseline right-parenthesis upper C Subscript i comma j comma k Baseline right-bracket period$ (23)

The value of in the model is while that in Rasch (1994) is .

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:

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 -velocity, we have:

while for the -velocity we have:

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

## Vertical Velocity

Having obtained a complete specification of the and 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:

and combining it with equation (6) to obtain:

Solving for and using the semi-discrete notation we obtain:

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. ROMS provides a choice of a nonlinear equation of state or a linear equation of state . 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 () 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.