Difference between revisions of "Numerical Solution Technique"

From WikiROMS
Jump to navigationJump to search
(→‎Equation of State: fixed reference)   (change visibility)
 
(5 intermediate revisions by 3 users not shown)
Line 2: Line 2:


==Horizontal Discretization==
==Horizontal Discretization==
<wikitex>
In the horizontal, the ROMS governing equations are discretized over a boundary-fitted, orthogonal curvilinear coordinates (<math>\xi</math>, <math>\eta</math>) 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 (<math>\xi</math>, <math>\eta</math>) grid is specified in the metric terms ([[Variables#pm |pm]], [[Variables#pn |pn]]).
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 ([[Variables#pm |pm]], [[Variables#pn |pn]]).


The model state variables are staggered using an [[Bibliography#ArakawaA_1977a |Arakawa C-grid]]. As illustrated below, the free-surface ([[Variables#zeta |zeta]]), density ([[Variables#rho |rho]]), and active/passive tracers ([[Variables#t |t]]) are located at the center of the cell whereas the horizontal velocity ([[Variables#u |u]] and [[Variables#v |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.
The model state variables are staggered using an [[Bibliography#ArakawaA_1977a |Arakawa C-grid]]. As illustrated below, the free-surface ([[Variables#zeta |zeta]]), density ([[Variables#rho |rho]]), and active/passive tracers ([[Variables#t |t]]) are located at the center of the cell whereas the horizontal velocity ([[Variables#u |u]] and [[Variables#v |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.
Line 20: Line 19:
! Full Range
! Full Range
|-
|-
| $\rho$-type
| <math>\rho</math>-type
| 1:[[Variables#Lm|Lm(ng)]],1:[[Variables#Mm|Mm(ng)]]
| 1:[[Variables#Lm|Lm(ng)]],1:[[Variables#Mm|Mm(ng)]]
| 0:[[Variables#L|L(ng)]],0:[[Variables#M|M(ng)]]
| 0:[[Variables#L|L(ng)]],0:[[Variables#M|M(ng)]]
|-
|-
| $\psi$-type
| <math>\psi</math>-type
| 2:[[Variables#Lm|Lm(ng)]],2:[[Variables#Mm|Mm(ng)]]
| 2:[[Variables#Lm|Lm(ng)]],2:[[Variables#Mm|Mm(ng)]]
| 1:[[Variables#L|L(ng)]],1:[[Variables#M|M(ng)]]
| 1:[[Variables#L|L(ng)]],1:[[Variables#M|M(ng)]]
Line 39: Line 38:


<div style="clear:both;"></div>
<div style="clear:both;"></div>
</wikitex>
 


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


==Vertical Discretization==
==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 ([[Variables#Hz |Hz]]) and volume. The model state variables are vertically staggered so that horizontal momentum ([[Variables#u |u]], [[Variables#v |v]]), ([[Variables#rho |rho]]), and active/passive tracers ([[Variables#t |t]]) are located at the center of the grid cell. The vertical velocity ([[Variables#W |omega]], [[Variables#wvel |w]]) and vertical mixing variables ([[Variables#Akt |Akt]], [[Variables#Akv |Akv]], etc) are located at the bottom and top faces of cell. See diagram below.
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 ([[Variables#Hz |Hz]]) and volume. The model state variables are vertically staggered so that horizontal momentum ([[Variables#u |u]], [[Variables#v |v]]), ([[Variables#rho |rho]]), and active/passive tracers ([[Variables#t |t]]) are located at the center of the grid cell. The vertical velocity ([[Variables#W |omega]], [[Variables#wvel |w]]) and vertical mixing variables ([[Variables#Akt |Akt]], [[Variables#Akv |Akv]], etc) are located at the bottom and top faces of cell. See diagram below.


Line 53: Line 53:




The total thickness of the water column is $\zeta(i,j)+h(i,j)$. The bathymetry ([[Variables#h |h]]) is usually time invariant whereas the free-surface ([[Variables#zeta |zeta]]) evolves in time. However, in sediment applications [[Variables#h |h]] changes with time when [[Options#SED_MORPH | SED_MORPH]] is activated. At input and output, the bathymetry is always a positive quantity. However, the depths [[Variables#z_r |z_r(i,j,k)]] and [[Variables#z_w |z_w(i,j,k)]] are negative for all locations below the mean sea level.
The total thickness of the water column is <math>\zeta(i,j)+h(i,j)</math>. The bathymetry ([[Variables#h |h]]) is usually time invariant whereas the free-surface ([[Variables#zeta |zeta]]) evolves in time. However, in sediment applications [[Variables#h |h]] changes with time when [[Options#SED_MORPH | SED_MORPH]] is activated. At input and output, the bathymetry is always a positive quantity. However, the depths [[Variables#z_r |z_r(i,j,k)]] and [[Variables#z_w |z_w(i,j,k)]] are negative for all locations below the mean sea level.
</wikitex>
 
==Masking of Land Areas==
==Masking of Land Areas==
<wikitex>
 
ROMS has the ability to work with interior land areas, although the
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 [[#Horizontal Discretization| above]] while two land cells are shown here:
computations occur over the entire model domain. The grid structure as a whole is shown [[#Horizontal Discretization| above]] while two land cells are shown here:
Line 63: Line 63:


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.
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===
===Velocity points===
<wikitex>
 
At the end of every timestep, the values of many variables within the
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
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
<math>u</math>, <math>v</math> or <math>\rho</math> points.  This is appropriate for the <math>v</math> 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
land should be zero.  It is likewise appropriate for the <math>u</math> point at
'''I''', but is not necessarily correct for point '''G'''.  The only
'''I''', but is not necessarily correct for point '''G'''.  The only
term in the $u$ equation that requires the $u$ value at point '''G'''
term in the <math>u</math> equation that requires the <math>u</math> value at point '''G'''
is the horizontal viscosity, which has a term of the form ${\partial
is the horizontal viscosity, which has a term of the form <math>{\partial
\over \partial \eta} \nu {\partial u \over \partial \eta}$.  Since
\over \partial \eta} \nu {\partial u \over \partial \eta}</math>.  Since
point '''G''' is used in this term by both points '''A''' and '''M''',
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
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
for '''A'''.  Instead, the term <math>{\partial u \over \partial \eta}</math> is
computed and the values at points '''D''' and '''K''' are replaced with
computed and the values at points '''D''' and '''K''' are replaced with
the values appropriate for either free-slip or no-slip boundary
the values appropriate for either free-slip or no-slip boundary
conditions.  Likewise, the term ${\partial \over \partial \xi} \nu
conditions.  Likewise, the term <math>{\partial \over \partial \xi} \nu
{\partial v \over \partial \xi}$ in the $v$ equation must be corrected
{\partial v \over \partial \xi}</math> in the <math>v</math> equation must be corrected
at the mask boundaries.
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.
This is accomplished by having a fourth mask array defined at the <math>\psi</math> 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 <math>u</math> at point '''A'''.  The desired value of <math>{\partial u \over \partial \eta}</math> is therefore <math>2 u_{\bf A}</math> while instead we have simply <math>u_{\bf A}</math>. 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 <math>{\partial u \over \partial \eta}</math> there will acquire the desired value of <math>-2 u_{\bf M}</math>. The corner point '''F''' is set to have a value of 1.
</wikitex>
 
===Temperature, Salinity and Surface Elevation===
===Temperature, Salinity and Surface Elevation===
<wikitex>
 
The handling of masks by the temperature, salinity and surface
The handling of masks by the temperature, salinity and surface
elevation equations is similar to that in the momentum equations, and
elevation equations is similar to that in the momentum equations, and
is in fact simpler.  Values of $T$, $S$ and $\zeta$ inside the land
is in fact simpler.  Values of <math>T</math>, <math>S</math> and <math>\zeta</math> inside the land
masks, such as point '''H''' in the figure above, are set to zero
masks, such as point '''H''' in the figure above, are set to zero
after every timestep.  This point would be used by the horizontal
after every timestep.  This point would be used by the horizontal
Line 95: Line 95:
handled by setting the first derivative terms at points '''E''',  
handled by setting the first derivative terms at points '''E''',  
'''I''', and '''L''' to zero, to be consistent with a no-flux boundary
'''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$
condition. Note that the equation of state must be able to handle <math>T = S = 0</math>
since this is the value inside masked regions.
since this is the value inside masked regions.


Line 103: Line 103:
come in and go out, for instance. Cells which are masked out as above
come in and go out, for instance. Cells which are masked out as above
are never allowed to be wet, however.
are never allowed to be wet, however.
*In the case of wetting and drying, a critical depth, $D_{crit}$, is supplied by the user.
*In the case of wetting and drying, a critical depth, <math>D_{crit}</math>, 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 total water depth (<math>D=h+\zeta</math>) is compared to <math>D_{crit}</math>. 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.
*The wetting and drying only happens during the 2-D computations; the 3-D computations see a depth of <math>D_{crit}</math> 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,
==Time-Stepping Overview==
* '''Mosaics Grids''' which connect several grids along their edges, and
* '''Composite Grids''' which allow overlap regions of aligned and non-aligned grids.


</wikitex>
While time-stepping the model, we have a stored history of the model fields at time <math>n-1</math>, an estimate of the fields at the current time <math>n</math>, and we need to come up with an estimate for time <math>n+1</math>. 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
===Refinement Grids===
3-D equations. There is an integer ratio <math>M</math> 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:
<wikitex>
{|
|-
|[[Image:refinement.png|thumb|330px|<center>Blue and Red Grids</center>]]
|[[Image:refinement_grid1.png|thumb|330px|<center>Red Grid Contact Points</center>]]
|[[Image:refinement_grid2.png|thumb|330px|<center>Blue Grid Contact Points</center>]]
|}
</wikitex>
===Mosaic Grids===
<wikitex>
{|
|-
|[[Image:mosaic.png|thumb|330px|<center>Blue and Red Grids</center>]]
|[[Image:mosaic_grid2.png|thumb|330px|<center>Red Grid Contact Points</center>]]
|[[Image:mosaic_grid1.png|thumb|330px|<center>Blue Grid Contact Points</center>]]
|}
</wikitex>
===Composite Grids===
<wikitex>
{|
|-
|[[Image:composite.png|thumb|330px|<center>Red and Blue Grids</center>]]
|[[Image:composite_grid1.png|thumb|330px|<center>Red Grid Contact Points</center>]]
|[[Image:composite_grid2.png|thumb|330px|<center>Blue Grid Contact Points</center>]]
|}
</wikitex>


==Time-Stepping Overview==
#Take a predictor step for at least the 3-D tracers to time <math>n+1/2</math>.
<wikitex>
#Compute <math>\overline{\rho}</math> and <math>\rho^*</math> for use in the depth-integrated time-steps, from the density either at time <math>n</math> or time <math>n+1/2</math>.
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
#Depth integrate the 3-D momentum right-hand side terms at time <math>n+1/2</math> for use in the depth-integrated time-steps (or extrapolate to obtain an estimate of those terms).
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:
#Take all the depth-integrated steps. Store weighted time-means of the <math>\overline{u}</math>, <math>\overline{v}</math> fields centered at both time <math>n+1/2</math> and time <math>n+1</math> (plus <math>\zeta</math> at time <math>n+1</math>). The latter requires this time-stepping to extend past time <math>n+1</math>, using <math>M^*</math> steps rather than just <math>M</math>.
 
#Use the weighted time-means from depth-integrated fields to complete the corrector step for the 3-D fields to time <math>n+1</math>.
#Take a predictor step for at least the 3-D tracers to time $n+1/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$.
#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).
#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$.
#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
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
instability due to the use of shorter time steps for the depth-integrated
computations.
computations.


The mode coupling has evolved through the various ROMS versions, as shown in the figures below [from [[Bibliography#ShchepetkinAF_2008a | Shchepetkin and McWilliams (2008a)]]]. The time-stepping schemes are also listed in the following table and described in detail in [[Bibliography#ShchepetkinAF_2005a | Shchepetkin and McWilliams (2005)]] and [[Bibliography#ShchepetkinAF_2008b | Shchepetkin and McWilliams (2008b)]]; the relevant ones are described in [[Time-stepping Schemes Review]].
The mode coupling has evolved through the various ROMS versions, as shown in the figures below [from [[Bibliography#ShchepetkinAF_2008a | Shchepetkin and McWilliams (2008a)]]]. The time-stepping schemes are also listed in the following table and described in detail in [[Bibliography#ShchepetkinAF_2005a | Shchepetkin and McWilliams (2005)]] and [[Bibliography#ShchepetkinAF_2008b | 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 $n+1$, while the light red function has weights to produce an estimate at time $n+1/2$.
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 <math>n+1</math>, while the light red function has weights to produce an estimate at time <math>n+1/2</math>.


Rutgers ROMS:
Rutgers ROMS:
Line 191: Line 154:
feedback
feedback
| LF-AM3 with FB
| LF-AM3 with FB
feedback$\dagger$
feedback<math>\dagger</math>
| Gen. FB
| Gen. FB
(AB3-AM4)
(AB3-AM4)
Line 197: Line 160:
(AB3-AM4)
(AB3-AM4)
|-
|-
| 2-D $\alpha_{max}$, iter.
| 2-D <math>\alpha_{max}</math>, iter.
| $\sqrt{2}$, (2)$\ddagger$
| <math>\sqrt{2}</math>, (2)<math>\ddagger</math>
| 1.85, (2)
| 1.85, (2)
| 1.85, (2)
| 1.85, (2)
Line 230: Line 193:
(AB3-AM4)
(AB3-AM4)
|-
|-
| $\alpha_{max}$, advect.
| <math>\alpha_{max}</math>, advect.
| 0.72
| 0.72
| 0.72
| 0.72
Line 237: Line 200:
| 0.78
| 0.78
|-
|-
| $\alpha_{max}$, Cor.
| <math>\alpha_{max}</math>, Cor.
| 0.72
| 0.72
| 0.72
| 0.72
Line 244: Line 207:
| 0.78
| 0.78
|-
|-
| $\alpha_{max}$, int. w.
| <math>\alpha_{max}</math>, int. w.
| 0.72, (1)
| 0.72, (1)
| 1.14, (1,2)
| 1.14, (1,2)
Line 250: Line 213:
| 1.85, (2)
| 1.85, (2)
| 1.78, (1)
| 1.78, (1)
|-}
|}




{{note}} '''Note:'''
{{note}} '''Note:'''
$\dagger$The generalized FB barotropic mode was ported into the newest AGRIF code at the end of 2007.
<math>\dagger</math>The generalized FB barotropic mode was ported into the newest AGRIF code at the end of 2007.
 
<math>\ddagger</math>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.


$\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==
==Conservation Properties==
<wikitex>
 
From [[Bibliography#ShchepetkinAF_2005a | Shchepetkin and McWilliams (2005)]], we have a tracer concentration equation in advective form:
From [[Bibliography#ShchepetkinAF_2005a | 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)}$$
 
{| class="eqno"
|<math display="block">\frac{\partial C}{\partial t} + (u \cdot \nabla) C = 0 </math><!--\eqno{(1)}-->
|(1)
|}
 
and also a tracer concentration equation in conservation form:
and also a tracer concentration equation in conservation form:
$$\frac{\partial C}{\partial t} + \nabla \cdot (u C) = 0. \eqno{(2)}$$
 
{| class="eqno"
|<math display="block">\frac{\partial C}{\partial t} + \nabla \cdot (u C) = 0. </math><!--\eqno{(2)}-->
|(2)
|}
 
The continuity equation:
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
{| class="eqno"
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.
|<math display="block">( \nabla \cdot u) = 0 </math><!--\eqno{(3)}-->
|(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:
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
{| class="eqno"
centered finite-difference approximations to $\partial / \partial \xi$,
|<math display="block"> \begin{align} {\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) =\\
$\partial / \partial \eta$ and $\partial / \partial \sigma$ with the
&{ 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 \end{align}</math><!--\eqno{(4)}-->
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}$,
|(4)
$\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$:
Here <math>\delta_{\xi}</math>, <math>\delta_{\eta}</math> and <math>\delta_\sigma</math> denote simple centered finite-difference approximations to <math>\partial / \partial \xi</math>, <math>\partial / \partial \eta</math> and <math>\partial / \partial \sigma</math> with the differences taken over the distances <math>\Delta\xi</math>, <math>\Delta\eta</math> and <math>\Delta \sigma</math>, respectively.  <math>\Delta z</math> is the vertical distance from one <math>\rho</math> point to another. <math>\overline{ ( \qquad )}^{\xi}</math>, <math>\overline{ ( \qquad )}^{\eta}</math> and <math>\overline{ ( \qquad)}^\sigma</math> represent averages taken over the distances <math>\Delta\xi</math>, <math>\Delta \eta</math> and <math>\Delta \sigma</math>.
$$ C = \frac{mn}{H_z} \int_{\Delta V} \frac{H_z C}{mn} \delta \xi
 
   \, \delta \eta \, \delta \sigma$$
The finite volume version of the same equation is no different, except that a quantity <math>C</math> is defined as the volume-averaged concentration over the grid box <math>\Delta V</math>:
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.
 
<math display="block"> C = \frac{mn}{H_z} \int_{\Delta V} \frac{H_z C}{mn} \delta \xi
   \, \delta \eta \, \delta \sigma</math>
 
The quantity  <math>\left(\frac{u \overline{H_z}^\xi \overline{C}^\xi}{\overline{n}^\xi} \right)</math> 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 [[Bibliography#ArakawaA_1977a | 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.
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 [[Bibliography#ArakawaA_1977a | 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:
The timestep in eq. (4) is assumed to be from time <math>n</math> to time <math>n+1</math>, with the other terms being evaluated at time <math>n+1/2</math> for second-order accuracy. Setting <math>C</math> to 1 everywhere reduces eq. (4) to:
$$\frac{\partial}{\partial t} \left( \frac{H_z}{m n} \right)
 
{| class="eqno"
|<math display="block">\frac{\partial}{\partial t} \left( \frac{H_z}{m n} \right)
   + \delta_{\xi} \left(
   + \delta_{\xi} \left(
   \frac{u \overline{H_z}^\xi }{\overline{n}^\xi} \right)
   \frac{u \overline{H_z}^\xi }{\overline{n}^\xi} \right)
Line 296: Line 276:
   \frac{v \overline{H_z}^\eta}{\overline{m}^\eta} \right)
   \frac{v \overline{H_z}^\eta}{\overline{m}^\eta} \right)
   + \delta_\sigma \left(  
   + \delta_\sigma \left(  
   \frac{H_z \Omega}{m n} \right) = 0 \eqno{(5)}$$
   \frac{H_z \Omega}{m n} \right) = 0 </math><!--\eqno{(5)}-->
|(5)
|}
 
If this equation holds true for the step from time <math>n</math> to time <math>n+1</math>, then constancy preservation will hold.


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 <math>\frac{H_z}{m n}</math> (the latter is controlled by changes in <math>\zeta</math> in the barotropic mode computations). Here, <math>\frac{H_z \Omega}{m n}</math> is the finite-volume flux across the ''moving'' grid-box interface, vertically on the <math>w</math> grid.


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 <math>\Omega</math>, is:


The vertical integral of the continuity equation (5), using the vertical boundary conditions on $\Omega$, is:
{| class="eqno"
$$ \frac{\partial}{\partial t} \left( \frac{\zeta}{mn} \right) +
|<math display="block"> \frac{\partial}{\partial t} \left( \frac{\zeta}{mn} \right) +
   \delta_{\xi} \left(
   \delta_{\xi} \left(
   \frac{\overline{u} \overline{D}^\xi }{\overline{n}^\xi} \right)
   \frac{\overline{u} \overline{D}^\xi }{\overline{n}^\xi} \right)
   + \delta_{\eta} \left(
   + \delta_{\eta} \left(
   \frac{\overline{v} \overline{D}^\eta}{\overline{m}^\eta} \right)
   \frac{\overline{v} \overline{D}^\eta}{\overline{m}^\eta} \right)
   = 0 \eqno{(6)}$$
   = 0 </math><!--\eqno{(6)}-->
|(6)
|}
 
where <math>\zeta</math> is the surface elevation, <math>D= h+\zeta</math> is the total depth, and <math>\overline{u},\overline{v}</math> 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.


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==
==Depth-Integrated Equations==
<wikitex>The depth average of a quantity $A$ is given by:
The depth average of a quantity <math>A</math> is given by:
$$\overline{A} = \frac{1}{D} \int_{-1}^0 H_z A d\sigma \eqno{(7)} $$
 
{| class="eqno"
|<math display="block">\overline{A} = \frac{1}{D} \int_{-1}^0 H_z A d\sigma </math><!--\eqno{(7)}-->
|(7)
|}
 
where the overbar indicates a vertically averaged quantity and
where the overbar indicates a vertically averaged quantity and
$$ D \equiv \zeta(\xi, \eta, t) + h(\xi, \eta) \eqno{(8)} $$
 
{| class="eqno"
|<math display="block"> D \equiv \zeta(\xi, \eta, t) + h(\xi, \eta) </math><!--\eqno{(8)}-->
|(8)
|}
is the total depth of the water column.  The vertical integral of the momentum equations are:
is the total depth of the water column.  The vertical integral of the momentum equations are:
$$ \eqalign{  
 
{| class="eqno"
|<math display="block"> \begin{align}
\frac{\partial}{\partial t} \left( \frac{D \overline{u}}{mn} \right)
\frac{\partial}{\partial t} \left( \frac{D \overline{u}}{mn} \right)
   + \frac  
   + \frac  
Line 325: Line 322:
   + \frac  
   + \frac  
   {\partial}{\partial \eta} \left( \frac{D \overline{uv}}{m} \right)
   {\partial}{\partial \eta} \left( \frac{D \overline{uv}}{m} \right)
   -& \frac{Df\overline{v}}{mn} \cr
   -& \frac{Df\overline{v}}{mn} \\
   - \left[ \overline{vv}
   - \left[ \overline{vv}
   \frac{\partial}{\partial \xi}
   \frac{\partial}{\partial \xi}
Line 333: Line 330:
   - \frac{D}{n}
   - \frac{D}{n}
   &\left( \frac{\partial \overline{\phi_2}}{\partial \xi} +
   &\left( \frac{\partial \overline{\phi_2}}{\partial \xi} +
   g \frac{\partial \zeta}{\partial \xi} \right) \cr
   g \frac{\partial \zeta}{\partial \xi} \right) \\
     + \frac{ D}{mn}
     + \frac{ D}{mn}
   \left( \overline{\cal F}_u + \overline{\cal D}_{h_u} \right)  
   \left( \overline{\cal F}_u + \overline{\cal D}_{h_u} \right)  
   &+ \frac{1}{mn} \left( \tau^{\xi}_s - \tau^{\xi}_b \right) } \eqno{(9)} $$
   &+ \frac{1}{mn} \left( \tau^{\xi}_s - \tau^{\xi}_b \right) \end{align}
</math><!--\eqno{(9)}-->
|(9)
|} 
 
and
and
$$ \eqalign{
 
{| class="eqno"
|<math display="block"> \begin{align}
     \frac{\partial}{\partial t} \left( \frac{D \overline{v}}{mn} \right)
     \frac{\partial}{\partial t} \left( \frac{D \overline{v}}{mn} \right)
   + \frac  
   + \frac  
Line 344: Line 347:
   + \frac  
   + \frac  
   {\partial}{\partial \eta} \left( \frac{D \overline{vv}}{m} \right)
   {\partial}{\partial \eta} \left( \frac{D \overline{vv}}{m} \right)
   +& \frac{Df\overline{u}}{mn} \cr
   +& \frac{Df\overline{u}}{mn} \\
       + \left[ \overline{uv}
       + \left[ \overline{uv}
   \frac{\partial}{\partial \xi}
   \frac{\partial}{\partial \xi}
Line 352: Line 355:
   - \frac{D}{m}
   - \frac{D}{m}
   &\left( \frac{\partial \overline{\phi_2}}{\partial \eta} +
   &\left( \frac{\partial \overline{\phi_2}}{\partial \eta} +
   g \frac{\partial \zeta}{\partial \eta} \right) \cr
   g \frac{\partial \zeta}{\partial \eta} \right) \\
   + \frac{ D}{mn}
   + \frac{ D}{mn}
   \left( \overline{\cal F}_v + \overline{\cal D}_{h_v} \right)  
   \left( \overline{\cal F}_v + \overline{\cal D}_{h_v} \right)  
   &+ \frac{1}{mn} \left( \tau^{\eta}_s - \tau^{\eta}_b \right)} \eqno{(10)} $$
   &+ \frac{1}{mn} \left( \tau^{\eta}_s - \tau^{\eta}_b \right) \end{align}
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).
</math><!--\eqno{(10)}-->
|(10)
|}
 
where <math>\phi_2</math> includes the <math>\frac{\partial z}{\partial \xi}</math> term, <math>\overline{\cal D}_{h_u}</math> 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:
The presence of a free surface introduces waves which propagate at a speed of <math>\sqrt{gh}</math>.  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 <math>\overline{u}</math> and <math>\overline{v}</math> 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:
[[Image:Shortstep.png]]
[[Image:Shortstep.png]]


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:
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 <math>R_{u_{\rm slow}}</math> and <math>R_{v_{\rm slow}}</math>, equations (9) and (10) become:
$$ \eqalign{
 
{| class="eqno"
|<math display="block"> \begin{align}
   \frac{\partial}{\partial t} \left( \frac{D \overline{u}}{mn} \right)
   \frac{\partial}{\partial t} \left( \frac{D \overline{u}}{mn} \right)
   + \frac{\partial}{\partial \xi}
   + \frac{\partial}{\partial \xi}
Line 368: Line 377:
   + \frac{\partial}{\partial \eta}
   + \frac{\partial}{\partial \eta}
   \left( \frac{D \overline{u}\,\overline{v}}{m} \right)
   \left( \frac{D \overline{u}\,\overline{v}}{m} \right)
   &- \frac{Df\overline{v}}{mn} \cr
   &- \frac{Df\overline{v}}{mn} \\
   - \left[ \overline{v}\,\overline{v}
   - \left[ \overline{v}\,\overline{v}
   \frac{\partial}{\partial \xi}
   \frac{\partial}{\partial \xi}
Line 376: Line 385:
   \frac{gD}{n} \frac{\partial \zeta}{\partial \xi}
   \frac{gD}{n} \frac{\partial \zeta}{\partial \xi}
     + \frac{D}{mn} {\cal D}_{\overline{u}}
     + \frac{D}{mn} {\cal D}_{\overline{u}}
   - \frac{1}{mn} \tau^{\xi}_b} \eqno{(11)} $$
   - \frac{1}{mn} \tau^{\xi}_b \end{align}
</math><!--\eqno{(11)-->
|(11)
|}
 
and
and
$$ \eqalign{
 
{| class="eqno"
|<math display="block">\begin{align}
   \frac{\partial}{\partial t} \left( \frac{D \overline{v}}{mn} \right)
   \frac{\partial}{\partial t} \left( \frac{D \overline{v}}{mn} \right)
   + \frac{\partial}{\partial \xi}
   + \frac{\partial}{\partial \xi}
Line 384: Line 399:
   + \frac{\partial}{\partial \eta}
   + \frac{\partial}{\partial \eta}
   \left( \frac{D \overline{v}\,\overline{v}}{m} \right)
   \left( \frac{D \overline{v}\,\overline{v}}{m} \right)
   &+ \frac{Df\overline{u}}{mn} \cr
   &+ \frac{Df\overline{u}}{mn} \\
   + \left[ \overline{u}\,\overline{v}
   + \left[ \overline{u}\,\overline{v}
   \frac{\partial}{\partial \xi}
   \frac{\partial}{\partial \xi}
Line 392: Line 407:
   \frac{gD}{m} \frac{\partial \zeta}{\partial \eta}
   \frac{gD}{m} \frac{\partial \zeta}{\partial \eta}
     + \frac{D}{mn} {\cal D}_{\overline{v}}
     + \frac{D}{mn} {\cal D}_{\overline{v}}
   - \frac{1}{mn} \tau^{\eta}_b .} \eqno{(12)} $$
   - \frac{1}{mn} \tau^{\eta}_b . \end{align}
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}$.
</math><!--\eqno{(12)}-->
|(12)
|}


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:
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 <math>R_{u_{\rm slow}}</math> and <math>R_{v_{\rm slow}}</math>.  This was found to be the easiest way to retain the baroclinic contributions of the non-linear terms such as <math>\overline{uu} - \overline{u}\,\overline{u}</math>.
 
The model is time stepped from time <math>n</math> to time <math>n+1</math> by using short time steps on equations (11), (12) and (6). Equation (6) is time stepped first, so that an estimate of the new <math>D</math> 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 <math>(n+</math> '''dtfast''' <math>\times M^\star)</math> and while maintaining weighted averages of the values of <math>\overline{u}</math>, <math>\overline{v}</math> and <math>\zeta</math>.  The averages are used to replace the values at time <math>n+1</math> in both the baroclinic and barotropic modes, and for recomputing the vertical grid spacing <math>H_z</math>. The following figure shows one option for how these weights might look:
[[Image:Barostep.png]]
[[Image:Barostep.png]]


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:
The primary weights, <math>a_m</math>, are used to compute <math>\langle \zeta \rangle^{n+1} \equiv \sum_{m=1}^{M^\star} a_m \zeta^m</math>. There is a related set of secondary weights <math>b_m</math>, used as <math>\langle \! \langle \overline{u} \rangle \! \rangle^{n+\frac{1}{2}} \equiv \sum_{m=1}^{M^\star} b_m \overline{u}^m</math>. 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
 
{| class="eqno"
|<math display="block">\begin{align}
  \langle \zeta \rangle_{i,j}^{n+1} = \langle \zeta \rangle_{i,j}^n &- \\
   (mn)_{i,j} \Delta t &\left[ \left\langle \!\! \left\langle
   (mn)_{i,j} \Delta t &\left[ \left\langle \!\! \left\langle
   \frac{D\overline u}{n} \right\rangle \!\!
   \frac{D\overline u}{n} \right\rangle \!\!
Line 409: Line 431:
   \right\rangle_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}
   \right\rangle_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}
   - \left\langle \!\! \left\langle \frac{D\overline v}{m} \right\rangle
   - \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)} $$
   \!\! \right\rangle_{i,j-\frac{1}{2}}^{n+\frac{1}{2}} \right] \end{align}
</math><!--\eqno{(13)}-->
|(13)
|}
 
[[Bibliography#ShchepetkinAF_2005a | Shchepetkin and McWilliams (2005)]] introduce a range of possible weights, but the ones used here have a shape function:
[[Bibliography#ShchepetkinAF_2005a | 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-
 
{| class="eqno"
|<math display="block"> 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}
   \left(\frac{\tau}{\tau_0} \right)^q \right] - r \frac{\tau}{\tau_0}
   \right\} \eqno{(14)} $$
   \right\} </math><!--\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,
|(14)
$$ 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)} $$
where <math>p, q</math> are parameters and <math>A_0, \tau_0</math>, and <math>r</math> are chosen to satisfy normalization, consistency, and second-order accuracy conditions,
compute $A(\tau)$ using eq.~(14), normalize using:
 
$$ \sum_{m=1}^{M^\star} a_m \equiv 1, \quad
{| class="eqno"
   \sum_{m=1}^{M^\star} a_m\frac{m}{M} \equiv 1,  \eqno{(16)} $$
|<math display="block"> I_n = \int_0^{\tau^\star} \tau^n A(\tau) d \tau = 1, \quad n=0,1,2 </math><!--\eqno{(15)}-->
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$.
|(15)
|}
 
using Newton iterations. <math>\tau^\star</math> is the upper limit of <math>\tau</math> with <math>A(\tau) \geq 0</math>. In practice we initially set
 
<math display="block"> A_0 = 1, r = 0, \quad \hbox{and} \quad \tau = \frac{(p+2)(p+q+2)}{(p+1)(p+q+1)} </math>
 
compute <math>A(\tau)</math> using eq.~(14), normalize using:
 
{| class="eqno"
|<math display="block"> \sum_{m=1}^{M^\star} a_m \equiv 1, \quad
   \sum_{m=1}^{M^\star} a_m\frac{m}{M} \equiv 1,  </math><!--\eqno{(16)}-->
|(16)
|}
 
and adjust <math>r</math> iteratively to satisfy the <math>n=2</math> condition of (15). We are using values of <math>p=2</math>, <math>q=4</math>, and <math>r=0.284</math>. This form allows some negative weights for small <math>m</math>, allowing <math>M^\star</math> to be less than <math>1.5M</math>.


ROMS also supports an older cosine weighting option, which isn't recommended since it is only first-order accurate.
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 ==
== 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:
Equation (11) contains the term <math>R_{u_{\rm slow}}</math>, 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} +
 
{| class="eqno"
|<math display="block"> -\frac{g D}{n} \frac{\partial \zeta}{\partial \xi} +
   \left[\frac{g D}{n} \frac{\partial \zeta}{\partial \xi} + {\cal F}
   \left[\frac{g D}{n} \frac{\partial \zeta}{\partial \xi} + {\cal F}
   \right] \eqno{(17)} $$
   \right] </math><!--\eqno{(17)}-->
|(17)
|}
 
where the term in square brackets is the mode coupling term and is
where the term in square brackets is the mode coupling term and is
held fixed over all the barotropic steps and
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)} $$
{| class="eqno"
|<math display="block"> {\cal F} = - \frac{1}{\rho_0 n} \int_{-h}^\zeta \frac{\partial
   P}{\partial \xi} dz </math><!--\eqno{(18)}-->
|(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.
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.


Line 440: Line 493:


Instead, let us define the following:
Instead, let us define the following:
$$\overline{\rho} = \frac{1}{D} \int_{-h}^\zeta \rho dz , \quad
 
{| class="eqno"
|<math display="block">\overline{\rho} = \frac{1}{D} \int_{-h}^\zeta \rho dz , \quad
   \rho^\star = \frac{1}{\frac{1}{2} D^2} \int_{-h}^\zeta
   \rho^\star = \frac{1}{\frac{1}{2} D^2} \int_{-h}^\zeta
   \left\{ \int_z^\zeta \rho dz^{\prime} \right\} dz \eqno{(19)} $$
   \left\{ \int_z^\zeta \rho dz^{\prime} \right\} dz </math><!--\eqno{(19)}-->
Changing the vertical coordinate to $\sigma$ yields:
|(19)
$$ \overline{\rho} =  \int_{-1}^0 \rho d\sigma , \quad
|}
 
Changing the vertical coordinate to <math>\sigma</math> yields:
 
{| class="eqno"
|<math display="block"> \overline{\rho} =  \int_{-1}^0 \rho d\sigma , \quad
   \rho^\star = 2 \int_{-1}^0
   \rho^\star = 2 \int_{-1}^0
   \left\{ \int_\sigma^0 \rho d\sigma^{\prime} \right\} d\sigma \eqno{(20)} $$
   \left\{ \int_\sigma^0 \rho d\sigma^{\prime} \right\} d\sigma </math><!--\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:
|(20)
$$ -\frac{1}{\rho_0} \frac{g}{n} \left\{ \frac{\partial}{\partial \xi}
|}
 
which implies that <math>\overline{\rho}</math> and <math>\rho^\star</math> are actually independent of <math>\zeta</math> as long as the density profile <math>\rho = \rho(\sigma)</math> does not change. The vertically integrated pressure gradient becomes:
 
{| class="eqno"
|<math display="block"> -\frac{1}{\rho_0} \frac{g}{n} \left\{ \frac{\partial}{\partial \xi}
   \left( \frac{\rho^\star D^2}{2} \right) - \overline{\rho} D
   \left( \frac{\rho^\star D^2}{2} \right) - \overline{\rho} D
   \frac{\partial h}{\partial \xi} \right\} = -\frac{1}{\rho_0}
   \frac{\partial h}{\partial \xi} \right\} = -\frac{1}{\rho_0}
   \frac{g}{n} D \left\{ \rho^\star \frac{\partial \zeta}{\partial \xi} +
   \frac{g}{n} D \left\{ \rho^\star \frac{\partial \zeta}{\partial \xi} +
   \frac{D}{2} \frac{\partial \rho^\star}{\partial \xi} + (\rho^\star -
   \frac{D}{2} \frac{\partial \rho^\star}{\partial \xi} + (\rho^\star -
   \overline{\rho}) \frac{\partial h}{\partial \xi} \right\} \eqno{(21)} $$
   \overline{\rho}) \frac{\partial h}{\partial \xi} \right\} </math><!--\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 [[Bibliography#ShchepetkinAF_2005a | Shchepetkin and McWilliams (2005)]].
|(21)
</wikitex>
|}
 
In the case of uniform density <math>\rho_0</math>, we obtain <math>\rho^\star \equiv \overline{\rho} \equiv \rho_0</math>, but we otherwise have two new terms. The accuracy of these terms depends on an accurate vertical integration of the density, as described in [[Bibliography#ShchepetkinAF_2005a | Shchepetkin and McWilliams (2005)]].
 
== Time Stepping: Internal Velocity Modes and Tracers ==
== 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 Review]]). 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 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 <math>u</math> and <math>v</math>. The depth-averaged component is then removed and replaced by the <math>\langle \overline{u} \rangle</math> and <math>\langle \overline{v} \rangle</math> 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 <math>n+\frac{1}{2}</math> 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 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 <math>n+\frac{1}{2}</math>, 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 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.


{{warning}} 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.
{{warning}} 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==
==Horizontal and Vertical Advection==
<wikitex>The advection of a tracer $C$ has an equation of the form
The advection of a tracer <math>C</math> has an equation of the form
$$
 
<math display="block">
   \frac{\partial}{\partial t} \frac{ H_z C}{mn} =
   \frac{\partial}{\partial t} \frac{ H_z C}{mn} =
   - \frac{\partial}{\partial \xi} F^\xi
   - \frac{\partial}{\partial \xi} F^\xi
   - \frac{\partial}{\partial \eta} F^\eta
   - \frac{\partial}{\partial \eta} F^\eta
   - \frac{\partial}{\partial \sigma} F^\sigma ,
   - \frac{\partial}{\partial \sigma} F^\sigma ,
$$
</math>
 
where we have introduced the advective fluxes:
where we have introduced the advective fluxes:
$$ \eqalign{
 
   F^\xi &= \frac{H_z u C}{n} \cr
<math display="block"> \begin{align}
   F^\eta &= \frac{H_z v C}{m} \cr
   F^\xi &= \frac{H_z u C}{n} \\
   F^\sigma &= \frac{H_z \Omega C}{mn} .}
   F^\eta &= \frac{H_z v C}{m} \\
$$
   F^\sigma &= \frac{H_z \Omega C}{mn} .\end{align}
</wikitex>
</math>
 
===Second-order Centered===
===Second-order Centered===
<wikitex>The simplest form of the advective fluxes is the centered second-order:
The simplest form of the advective fluxes is the centered second-order:
$$ \eqalign{
 
<math display="block"> \begin{align}
   F^\xi &= \frac{\overline{H_z}^\xi u
   F^\xi &= \frac{\overline{H_z}^\xi u
               \overline{C}^\xi}{\overline{n}^\xi} \cr
               \overline{C}^\xi}{\overline{n}^\xi} \\
   F^\eta &= \frac{\overline{H_z}^\eta v
   F^\eta &= \frac{\overline{H_z}^\eta v
               \overline{C}^\eta}{\overline{m}^\eta} \cr
               \overline{C}^\eta}{\overline{m}^\eta} \\
   F^\sigma &= \frac{\overline{H_z}^\sigma \Omega
   F^\sigma &= \frac{\overline{H_z}^\sigma \Omega
               \overline{C}^\sigma}{mn} .}
               \overline{C}^\sigma}{mn} . \end{align}
$$
</math>
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>
 
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===
===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:
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{
 
<math display="block"> \begin{align}
     G^{\xi} &= \overline{\left(\frac{\partial C}{\partial
     G^{\xi} &= \overline{\left(\frac{\partial C}{\partial
             \xi}\right)}^\xi \cr
             \xi}\right)}^\xi \\
     G^{\eta} &= \overline{\left(\frac{\partial C}{\partial
     G^{\eta} &= \overline{\left(\frac{\partial C}{\partial
             \eta}\right)}^\eta \cr
             \eta}\right)}^\eta \\
     G^{\sigma} &= \overline{\left(\frac{\partial C}{\partial
     G^{\sigma} &= \overline{\left(\frac{\partial C}{\partial
             \sigma}\right)}^\sigma . }
             \sigma}\right)}^\sigma . \end{align}
$$
</math>
 
The fluxes now become:
The fluxes now become:
$$ \eqalign{
 
{| class="eqno"
|<math display="block"> \begin{align}
   F^\xi &= \frac{\overline{H_z}^\xi} {\overline{n}^\xi} u \left(
   F^\xi &= \frac{\overline{H_z}^\xi} {\overline{n}^\xi} u \left(
               \overline{C}^\xi -
               \overline{C}^\xi -
           \frac{1}{3} \frac{\partial
           \frac{1}{3} \frac{\partial
               G^{\xi}}{\partial \xi} \right)  \cr
               G^{\xi}}{\partial \xi} \right)  \\
   F^\eta &= \frac{\overline{H_z}}{\overline{m}^\eta} ^\eta v \left(
   F^\eta &= \frac{\overline{H_z}}{\overline{m}^\eta} ^\eta v \left(
               \overline{C}^\eta-
               \overline{C}^\eta-
           \frac{1}{3} \frac{\partial
           \frac{1}{3} \frac{\partial
               G^{\eta}}{\partial \eta} \right) \cr
               G^{\eta}}{\partial \eta} \right) \\
   F^\sigma &= \frac{\overline{H_z}^\sigma} {mn} \Omega \left(
   F^\sigma &= \frac{\overline{H_z}^\sigma} {mn} \Omega \left(
               \overline{C}^\sigma -
               \overline{C}^\sigma -
           \frac{1}{3} \frac{\partial
           \frac{1}{3} \frac{\partial
               G^{\sigma}}{\partial \sigma} \right). } \eqno{(22)}  
               G^{\sigma}}{\partial \sigma} \right). \end{align}
$$
</math><!--\eqno{(22)}-->
</wikitex>
|(22)
|}
 
===Fourth-order Akima===
===Fourth-order Akima===
<wikitex>An alternate fourth-order algorithm is that by Akima:
An alternate fourth-order algorithm is that by Akima:
$$ \eqalign{
 
<math display="block"> \begin{align}
     G_{\xi} &= 2 {\frac{\partial C}{\partial \xi}_i
     G_{\xi} &= 2 {\frac{\partial C}{\partial \xi}_i
         \frac{\partial C}{\partial \xi}_{i+1}} \, \bigg/ \left(
         \frac{\partial C}{\partial \xi}_{i+1}} \, \bigg/ \left(
         \frac{\partial C}{\partial \xi}_i +
         \frac{\partial C}{\partial \xi}_i +
           \frac{\partial C}{\partial \xi}_{i+1} \right) \cr
           \frac{\partial C}{\partial \xi}_{i+1} \right) \\
     G_{\eta} &= 2 {\frac{\partial C}{\partial \eta}_j
     G_{\eta} &= 2 {\frac{\partial C}{\partial \eta}_j
         \frac{\partial C}{\partial \eta}_{j+1}} \, \bigg/ \left(
         \frac{\partial C}{\partial \eta}_{j+1}} \, \bigg/ \left(
         \frac{\partial C}{\partial \eta}_j +
         \frac{\partial C}{\partial \eta}_j +
           \frac{\partial C}{\partial \eta}_{j+1} \right) \cr
           \frac{\partial C}{\partial \eta}_{j+1} \right) \\
     G_{\sigma} &= 2 {\frac{\partial C}{\partial \sigma}_k
     G_{\sigma} &= 2 {\frac{\partial C}{\partial \sigma}_k
         \frac{\partial C}{\partial \sigma}_{k-1}} \, \bigg/ \left(
         \frac{\partial C}{\partial \sigma}_{k-1}} \, \bigg/ \left(
         \frac{\partial C}{\partial \sigma}_k +
         \frac{\partial C}{\partial \sigma}_k +
           \frac{\partial C}{\partial \sigma}_{k-1} \right) }
           \frac{\partial C}{\partial \sigma}_{k-1} \right) \end{align}
$$
</math>
 
With the fluxes as in eq. (22).
With the fluxes as in eq. (22).
</wikitex>
 
===Third-order Upwind===
===Third-order Upwind===
<wikitex>There is a class of third-order upwind advection schemes, both one-dimensional ([[Bibliography#LeonardBP_1979 |Leonard, 1979]]) and two-dimensional ([[Bibliography#RaschP_1994 |Rasch, 1994]] and [[Bibliography#ShchepetkinAF_1998a | Shchepetkin and McWilliams, 1998]]). This scheme is known as UTOPIA (Uniformly Third-Order Polynomial Interpolation Algorithm). Applying flux limiters to UTOPIA is explored in [[Bibliography#ThuburnJ_1995 | Thuburn (1995)]], 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:
There is a class of third-order upwind advection schemes, both one-dimensional ([[Bibliography#LeonardBP_1979 |Leonard, 1979]]) and two-dimensional ([[Bibliography#RaschP_1994 |Rasch, 1994]] and [[Bibliography#ShchepetkinAF_1998a | Shchepetkin and McWilliams, 1998]]). This scheme is known as UTOPIA (Uniformly Third-Order Polynomial Interpolation Algorithm). Applying flux limiters to UTOPIA is explored in [[Bibliography#ThuburnJ_1995 | Thuburn (1995)]], although it is not implemented in ROMS. The two-dimensional formulation in Rasch contains terms of order <math>u^2C</math> and <math>u^3C</math>, including cross terms (<math>uvC</math>). 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{
 
<math display="block"> \begin{align}
   F^\xi &= \frac{H_z u}{n} \left( C - \gamma \frac{\partial^2
   F^\xi &= \frac{H_z u}{n} \left( C - \gamma \frac{\partial^2
   C}{\partial \xi^2} \right) \cr
   C}{\partial \xi^2} \right) \\
   F^\eta &= \frac{H_z v}{m} \left( C - \gamma \frac{\partial^2
   F^\eta &= \frac{H_z v}{m} \left( C - \gamma \frac{\partial^2
   C}{\partial \eta^2} \right)}
   C}{\partial \eta^2} \right) \end{align}
$$
</math>
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:
 
$$
The second derivative terms are centered on a <math>\rho</math> point in the grid, but are needed at a <math>u</math> or <math>v</math> point in the flux. The upstream value is used:
 
{| class="eqno"
|<math display="block">
   F^\xi_{i,j,k} = \frac{\overline{H_z}^\xi}{\overline{n}^\xi}
   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} +
   \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)}
   \min(0,u_{i,j,k}) C_{i,j,k} \right] .
$$
</math><!--\eqno{(23)}-->
The value of $\gamma$ in the model is
|(23)
$\frac{1}{8}$ while that in [[Bibliography#RaschP_1994 |Rasch (1994)]] is $\frac{1}{6}$.
|}
 
The value of <math>\gamma</math> in the model is
<math>\frac{1}{8}</math> while that in [[Bibliography#RaschP_1994 |Rasch (1994)]] is <math>\frac{1}{6}</math>.


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:
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:
$$
 
<math display="block">
   F^s = \frac{H_z w}{mn} \left[
   F^s = \frac{H_z w}{mn} \left[
     - \frac{1}{16} C_{i,j,k-1} + \frac{9}{16} C_{i,j,k} +
     - \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]
       \frac{9}{16} C_{i,j,k+1} - \frac{1}{16} C_{i,j,k+2} \right]
$$
</math>


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:
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 <math>u</math>-velocity, we have:
$$ \eqalign{
<math display="block"> \begin{align}
   F^\xi &= \left(u - \gamma \frac{\partial^2 u}{\partial \xi^2} \right)
   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} - \gamma \frac{\partial^2}{\partial \xi^2}
   \left( \frac{H_z u}{n} \right) \right] \cr
   \left( \frac{H_z u}{n} \right) \right] \\
   F^\eta &= \left(u - \gamma \frac{\partial^2 u}{\partial \eta^2}
   F^\eta &= \left(u - \gamma \frac{\partial^2 u}{\partial \eta^2}
     \right)
     \right)
   \left[ \frac{H_z v}{m} - \gamma \frac{\partial^2}{\partial \xi^2}
   \left[ \frac{H_z v}{m} - \gamma \frac{\partial^2}{\partial \xi^2}
   \left( \frac{H_z v}{m} \right) \right] \cr
   \left( \frac{H_z v}{m} \right) \right] \\
   F^\sigma &= \frac{H_z w}{mn} \left[
   F^\sigma &= \frac{H_z w}{mn} \left[
     - \frac{1}{16} u_{i,j,k-1} + \frac{9}{16} u_{i,j,k} +
     - \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] }
       \frac{9}{16} u_{i,j,k+1} - \frac{1}{16} u_{i,j,k+2} \right] \end{align}
$$
</math>
while for the $v$-velocity we have:
 
$$ \eqalign{
while for the <math>v</math>-velocity we have:
 
<math display="block"> \begin{align}
   F^\xi &= \left(v - \gamma \frac{\partial^2 v}{\partial \xi^2} \right)
   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} - \gamma \frac{\partial^2}{\partial \eta^2}
   \left( \frac{H_z u}{n} \right) \right] \cr
   \left( \frac{H_z u}{n} \right) \right] \\
   F^\eta &= \left(v - \gamma \frac{\partial^2 v}{\partial \eta^2}
   F^\eta &= \left(v - \gamma \frac{\partial^2 v}{\partial \eta^2}
     \right)
     \right)
   \left[ \frac{H_z v}{m} - \gamma \frac{\partial^2}{\partial \eta^2}
   \left[ \frac{H_z v}{m} - \gamma \frac{\partial^2}{\partial \eta^2}
   \left( \frac{H_z v}{m} \right) \right] \cr
   \left( \frac{H_z v}{m} \right) \right] \\
   F^\sigma &= \frac{H_z w}{mn} \left[
   F^\sigma &= \frac{H_z w}{mn} \left[
     - \frac{1}{16} v_{i,j,k-1} + \frac{9}{16} v_{i,j,k} +
     - \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] }
       \frac{9}{16} v_{i,j,k+1} - \frac{1}{16} v_{i,j,k+2} \right] \end{align}
$$
</math>
 
In all these terms, the second derivatives are evaluated at an upstream location.
In all these terms, the second derivatives are evaluated at an upstream location.
</wikitex>
 


===MPDATA===
===MPDATA===
<wikitex></wikitex>
 
==Vertical Velocity==
==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
Having obtained a complete specification of the <math>u,v,T,</math> and <math>S</math> 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:
is obtained by rewriting equation (5) as:
$$\frac{\partial}{\partial t} \left( \frac{\zeta}{m n} \right)
 
<math display="block">\frac{\partial}{\partial t} \left( \frac{\zeta}{m n} \right)
   + \delta_{\xi} \left(
   + \delta_{\xi} \left(
   \frac{u \overline{H_z}^\xi }{\overline{n}^\xi} \right)
   \frac{u \overline{H_z}^\xi }{\overline{n}^\xi} \right)
Line 600: Line 694:
   \frac{v \overline{H_z}^\eta}{\overline{m}^\eta} \right)
   \frac{v \overline{H_z}^\eta}{\overline{m}^\eta} \right)
   + \delta_\sigma \left(  
   + \delta_\sigma \left(  
   \frac{H_z \Omega}{m n} \right) = 0 $$
   \frac{H_z \Omega}{m n} \right) = 0 </math>
 
and combining it with equation (6) to obtain:
and combining it with equation (6) to obtain:
$$
 
<math display="block">
   \frac{\partial}{\partial \xi} \left( \frac{H_z u}{n} \right) +
   \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 \eta} \left( \frac{H_z v}{m} \right) +
Line 610: Line 706:
   \frac{\partial}{\partial \eta} \left( \frac{D \overline{v}}{m}
   \frac{\partial}{\partial \eta} \left( \frac{D \overline{v}}{m}
   \right)  = 0 .
   \right)  = 0 .
$$
</math>
Solving for $H_z \Omega / mn$ and using the semi-discrete notation we obtain:
 
$$
Solving for <math>H_z \Omega / mn</math> and using the semi-discrete notation we obtain:
 
<math display="block">
   \frac{H_z \Omega}{mn} =  \int  \left[
   \frac{H_z \Omega}{mn} =  \int  \left[
   \delta_{\xi} \left( \frac{\overline{u} \overline{D}^{\xi}}
   \delta_{\xi} \left( \frac{\overline{u} \overline{D}^{\xi}}
Line 622: Line 720:
   \delta_{\eta} \left( \frac{v \overline{H_z}^{\eta}}
   \delta_{\eta} \left( \frac{v \overline{H_z}^{\eta}}
   {\overline{m}^{\eta}} \right)  \right] d\sigma .
   {\overline{m}^{\eta}} \right)  \right] d\sigma .
$$
</math>
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?]
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>
 


==Equation of State==
==Equation of State==
<wikitex>The density is obtained from temperature and salinity via an equation of state. ROMS provides a choice of a nonlinear equation of state $\rho = \rho(T,S,z)$ or a linear equation of state $\rho = \rho(T)$. The nonlinear equation of state has been modified and now corresponds to the UNESCO equation of state as derived by [[Bibliography#JackettDR_1995a |Jackett and McDougall (1995)]]. It computes '' in situ'' density as a function of potential temperature, salinity and pressure.
The density is obtained from temperature and salinity via an equation of state. ROMS provides a choice of a nonlinear equation of state <math>\rho = \rho(T,S,z)</math> or a linear equation of state <math>\rho = \rho(T)</math>. The nonlinear equation of state has been modified and now corresponds to the UNESCO equation of state as derived by [[Bibliography#JackettDR_1995a |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 (<math>\rho = \rho(T)</math>) 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.


Warning: although we have used it quite extensively in the past, McDougall (personal communication) claims that the single-variable ($\rho = \rho(T)$) 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.
</wikitex>


==Vertical Friction and Diffusion==
==Vertical Friction and Diffusion==
<wikitex></wikitex>

Latest revision as of 13:38, 6 May 2016

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.


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
-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)


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


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.


Vieste-Dubrovnik Transect
ROMS staggered vertical grid


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.

Masking of Land Areas

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.

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 .

Rutgers ROMS: Timestep 1.png

ROMS AGRIF: Timestep 2.png

UCLA ROMS: Timestep 3.png

Non-hydrostatic ROMS: Timestep 4.png

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

(1)

and also a tracer concentration equation in conservation form:

(2)

The continuity equation:

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

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

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

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

(7)

where the overbar indicates a vertically averaged quantity and

(8)

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

(9)

and

(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: Shortstep.png

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:

(11)

and

(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: Barostep.png

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:

(13)

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

(14)

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

(15)

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

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

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

(17)

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

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

(19)

Changing the vertical coordinate to yields:

(20)

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

(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.

Warning 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.


Horizontal and Vertical Advection

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:

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

(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.


MPDATA

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.


Vertical Friction and Diffusion