s_algorithm/text/mom_fluxform.tex

% $Header: /u/gcmpack/manual/part2/mom_fluxform.tex,v 1.7 2001/11/13 18:15:26 adcroft Exp $
% $Name:  $

\section{Flux-form momentum equations}
\label{sec:flux-form_momentum_eqautions}
\begin{rawhtml}
<!-- CMIREDIR:flux-form_momentum_eqautions: -->
\end{rawhtml}

The original finite volume model was based on the Eulerian flux form
momentum equations. This is the default though the vector invariant
form is optionally available (and recommended in some cases).

The ``G's'' (our colloquial name for all terms on rhs!) are broken
into the various advective, Coriolis, horizontal dissipation, vertical
dissipation and metric forces:
\marginpar{$G_u$: {\bf Gu} }
\marginpar{$G_v$: {\bf Gv} }
\marginpar{$G_w$: {\bf Gw} }
\begin{eqnarray}
G_u & = & G_u^{adv} + G_u^{cor} + G_u^{h-diss} + G_u^{v-diss} +
G_u^{metric} + G_u^{nh-metric} \label{eq:gsplit_momu} \\
G_v & = & G_v^{adv} + G_v^{cor} + G_v^{h-diss} + G_v^{v-diss} +
G_v^{metric} + G_v^{nh-metric} \label{eq:gsplit_momv} \\
G_w & = & G_w^{adv} + G_w^{cor} + G_w^{h-diss} + G_w^{v-diss} +
G_w^{metric} + G_w^{nh-metric} \label{eq:gsplit_momw}
\end{eqnarray}
In the hydrostatic limit, $G_w=0$ and $\epsilon_{nh}=0$, reducing the
vertical momentum to hydrostatic balance.

These terms are calculated in routines called from subroutine {\em
CALC\_MOM\_RHS} a collected into the global arrays {\bf Gu}, {\bf Gv},
and {\bf Gw}.

\fbox{ \begin{minipage}{4.75in}
{\em S/R CALC\_MOM\_RHS} ({\em pkg/mom\_fluxform/calc\_mom\_rhs.F})

$G_u$: {\bf Gu} ({\em DYNVARS.h})

$G_v$: {\bf Gv} ({\em DYNVARS.h})

$G_w$: {\bf Gw} ({\em DYNVARS.h})
\end{minipage} }


\subsection{Advection of momentum}

The advective operator is second order accurate in space:
\begin{eqnarray}
{\cal A}_w \Delta r_f h_w G_u^{adv} & = &
  \delta_i \overline{ U }^i \overline{ u }^i
+ \delta_j \overline{ V }^i \overline{ u }^j
+ \delta_k \overline{ W }^i \overline{ u }^k \label{eq:discrete-momadvu} \\
{\cal A}_s \Delta r_f h_s G_v^{adv} & = &
  \delta_i \overline{ U }^j \overline{ v }^i
+ \delta_j \overline{ V }^j \overline{ v }^j
+ \delta_k \overline{ W }^j \overline{ v }^k \label{eq:discrete-momadvv} \\
{\cal A}_c \Delta r_c G_w^{adv} & = &
  \delta_i \overline{ U }^k \overline{ w }^i
+ \delta_j \overline{ V }^k \overline{ w }^j
+ \delta_k \overline{ W }^k \overline{ w }^k \label{eq:discrete-momadvw}
\end{eqnarray}
and because of the flux form does not contribute to the global budget
of linear momentum. The quantities $U$, $V$ and $W$ are volume fluxes
defined:
\marginpar{$U$: {\bf uTrans} }
\marginpar{$V$: {\bf vTrans} }
\marginpar{$W$: {\bf rTrans} }
\begin{eqnarray}
U & = & \Delta y_g \Delta r_f h_w u \label{eq:utrans} \\
V & = & \Delta x_g \Delta r_f h_s v \label{eq:vtrans} \\
W & = & {\cal A}_c w \label{eq:rtrans}
\end{eqnarray}
The advection of momentum takes the same form as the advection of
tracers but by a translated advective flow. Consequently, the
conservation of second moments, derived for tracers later, applies to
$u^2$ and $v^2$ and $w^2$ so that advection of momentum correctly
conserves kinetic energy.

\fbox{ \begin{minipage}{4.75in}
{\em S/R MOM\_U\_ADV\_UU} ({\em mom\_u\_adv\_uu.F})

{\em S/R MOM\_U\_ADV\_VU} ({\em mom\_u\_adv\_vu.F})

{\em S/R MOM\_U\_ADV\_WU} ({\em mom\_u\_adv\_wu.F})

{\em S/R MOM\_U\_ADV\_UV} ({\em mom\_u\_adv\_uv.F})

{\em S/R MOM\_U\_ADV\_VV} ({\em mom\_u\_adv\_vv.F})

{\em S/R MOM\_U\_ADV\_WV} ({\em mom\_u\_adv\_wv.F})

$uu$, $uv$, $vu$, $vv$: {\bf aF} (local to {\em calc\_mom\_rhs.F})
\end{minipage} }


\subsection{Coriolis terms}

The ``pure C grid'' Coriolis terms (i.e. in absence of C-D scheme) are
discretized:
\begin{eqnarray}
{\cal A}_w \Delta r_f h_w G_u^{Cor} & = &
  \overline{ f {\cal A}_c \Delta r_f h_c \overline{ v }^j }^i
- \epsilon_{nh} \overline{ f' {\cal A}_c \Delta r_f h_c \overline{ w }^k }^i \\
{\cal A}_s \Delta r_f h_s G_v^{Cor} & = &
- \overline{ f {\cal A}_c \Delta r_f h_c \overline{ u }^i }^j \\
{\cal A}_c \Delta r_c G_w^{Cor} & = &
 \epsilon_{nh} \overline{ f' {\cal A}_c \Delta r_f h_c \overline{ u }^i }^k
\end{eqnarray}
where the Coriolis parameters $f$ and $f'$ are defined:
\begin{eqnarray}
f & = & 2 \Omega \sin{\varphi} \\
f' & = & 2 \Omega \cos{\varphi}
\end{eqnarray}
where $\varphi$ is geographic latitude when using spherical geometry,
otherwise the $\beta$-plane definition is used:
\begin{eqnarray}
f & = & f_o + \beta y \\
f' & = & 0
\end{eqnarray}

This discretization globally conserves kinetic energy. It should be
noted that despite the use of this discretization in former
publications, all calculations to date have used the following
different discretization:
\begin{eqnarray}
G_u^{Cor} & = &
  f_u \overline{ v }^{ji}
- \epsilon_{nh} f_u' \overline{ w }^{ik} \\
G_v^{Cor} & = &
- f_v \overline{ u }^{ij} \\
G_w^{Cor} & = &
 \epsilon_{nh} f_w' \overline{ u }^{ik}
\end{eqnarray}
\marginpar{Need to change the default in code to match this}
where the subscripts on $f$ and $f'$ indicate evaluation of the
Coriolis parameters at the appropriate points in space. The above
discretization does {\em not} conserve anything, especially energy and
for historical reasons is the default for the code. A
flag controls this discretization: set run-time logical {\bf
useEnergyConservingCoriolis} to {\em true} which otherwise defaults to
{\em false}.

\fbox{ \begin{minipage}{4.75in}
{\em S/R MOM\_CDSCHEME} ({\em mom\_cdscheme.F})

{\em S/R MOM\_U\_CORIOLIS} ({\em mom\_u\_coriolis.F})

{\em S/R MOM\_V\_CORIOLIS} ({\em mom\_v\_coriolis.F})

$G_u^{Cor}$, $G_v^{Cor}$: {\bf cF} (local to {\em calc\_mom\_rhs.F})
\end{minipage} }


\subsection{Curvature metric terms}

The most commonly used coordinate system on the sphere is the
geographic system $(\lambda,\varphi)$. The curvilinear nature of these
coordinates on the sphere lead to some ``metric'' terms in the
component momentum equations. Under the thin-atmosphere and
hydrostatic approximations these terms are discretized:
\begin{eqnarray}
{\cal A}_w \Delta r_f h_w G_u^{metric} & = &
  \overline{ \frac{ \overline{u}^i }{a} \tan{\varphi} {\cal A}_c \Delta r_f h_c \overline{ v }^j }^i \\
{\cal A}_s \Delta r_f h_s G_v^{metric} & = &
- \overline{ \frac{ \overline{u}^i }{a} \tan{\varphi} {\cal A}_c \Delta r_f h_c \overline{ u }^i }^j \\
G_w^{metric} & = & 0
\end{eqnarray}
where $a$ is the radius of the planet (sphericity is assumed) or the
radial distance of the particle (i.e. a function of height).  It is
easy to see that this discretization satisfies all the properties of
the discrete Coriolis terms since the metric factor $\frac{u}{a}
\tan{\varphi}$ can be viewed as a modification of the vertical Coriolis
parameter: $f \rightarrow f+\frac{u}{a} \tan{\varphi}$.

However, as for the Coriolis terms, a non-energy conserving form has
exclusively been used to date:
\begin{eqnarray}
G_u^{metric} & = & \frac{u \overline{v}^{ij} }{a} \tan{\varphi} \\
G_v^{metric} & = & \frac{ \overline{u}^{ij} \overline{u}^{ij}}{a} \tan{\varphi}
\end{eqnarray}
where $\tan{\varphi}$ is evaluated at the $u$ and $v$ points
respectively.

\fbox{ \begin{minipage}{4.75in}
{\em S/R MOM\_U\_METRIC\_SPHERE} ({\em mom\_u\_metric\_sphere.F})

{\em S/R MOM\_V\_METRIC\_SPHERE} ({\em mom\_v\_metric\_sphere.F})

$G_u^{metric}$, $G_v^{metric}$: {\bf mT} (local to {\em calc\_mom\_rhs.F})
\end{minipage} }


\subsection{Non-hydrostatic metric terms}

For the non-hydrostatic equations, dropping the thin-atmosphere
approximation re-introduces metric terms involving $w$ and are
required to conserve angular momentum:
\begin{eqnarray}
{\cal A}_w \Delta r_f h_w G_u^{metric} & = &
- \overline{ \frac{ \overline{u}^i \overline{w}^k }{a} {\cal A}_c \Delta r_f h_c }^i \\
{\cal A}_s \Delta r_f h_s G_v^{metric} & = &
- \overline{ \frac{ \overline{v}^j \overline{w}^k }{a} {\cal A}_c \Delta r_f h_c}^j \\
{\cal A}_c \Delta r_c G_w^{metric} & = &
  \overline{ \frac{ {\overline{u}^i}^2 + {\overline{v}^j}^2}{a} {\cal A}_c \Delta r_f h_c }^k
\end{eqnarray}

Because we are always consistent, even if consistently wrong, we have,
in the past, used a different discretization in the model which is:
\begin{eqnarray}
G_u^{metric} & = &
- \frac{u}{a} \overline{w}^{ik} \\
G_v^{metric} & = &
- \frac{v}{a} \overline{w}^{jk} \\
G_w^{metric} & = &
  \frac{1}{a} ( {\overline{u}^{ik}}^2 + {\overline{v}^{jk}}^2 )
\end{eqnarray}

\fbox{ \begin{minipage}{4.75in}
{\em S/R MOM\_U\_METRIC\_NH} ({\em mom\_u\_metric\_nh.F})

{\em S/R MOM\_V\_METRIC\_NH} ({\em mom\_v\_metric\_nh.F})

$G_u^{metric}$, $G_v^{metric}$: {\bf mT} (local to {\em calc\_mom\_rhs.F})
\end{minipage} }


\subsection{Lateral dissipation}

Historically, we have represented the SGS Reynolds stresses as simply
down gradient momentum fluxes, ignoring constraints on the stress
tensor such as symmetry.
\begin{eqnarray}
{\cal A}_w \Delta r_f h_w G_u^{h-diss} & = &
  \delta_i  \Delta y_f \Delta r_f h_c \tau_{11}
+ \delta_j  \Delta x_v \Delta r_f h_\zeta \tau_{12} \\
{\cal A}_s \Delta r_f h_s G_v^{h-diss} & = &
  \delta_i  \Delta y_u \Delta r_f h_\zeta \tau_{21}
+ \delta_j  \Delta x_f \Delta r_f h_c \tau_{22}
\end{eqnarray}
\marginpar{Check signs of stress definitions}

The lateral viscous stresses are discretized:
\begin{eqnarray}
\tau_{11} & = & A_h c_{11\Delta}(\varphi) \frac{1}{\Delta x_f} \delta_i u
               -A_4 c_{11\Delta^2}(\varphi) \frac{1}{\Delta x_f} \delta_i \nabla^2 u \\
\tau_{12} & = & A_h c_{12\Delta}(\varphi) \frac{1}{\Delta y_u} \delta_j u
               -A_4 c_{12\Delta^2}(\varphi)\frac{1}{\Delta y_u} \delta_j \nabla^2 u \\
\tau_{21} & = & A_h c_{21\Delta}(\varphi) \frac{1}{\Delta x_v} \delta_i v
               -A_4 c_{21\Delta^2}(\varphi) \frac{1}{\Delta x_v} \delta_i \nabla^2 v \\
\tau_{22} & = & A_h c_{22\Delta}(\varphi) \frac{1}{\Delta y_f} \delta_j v
               -A_4 c_{22\Delta^2}(\varphi) \frac{1}{\Delta y_f} \delta_j \nabla^2 v
\end{eqnarray}
where the non-dimensional factors $c_{lm\Delta^n}(\varphi), \{l,m,n\} \in
\{1,2\}$ define the ``cosine'' scaling with latitude which can be
applied in various ad-hoc ways. For instance, $c_{11\Delta} =
c_{21\Delta} = (\cos{\varphi})^{3/2}$, $c_{12\Delta}=c_{22\Delta}=0$ would
represent the an-isotropic cosine scaling typically used on the
``lat-lon'' grid for Laplacian viscosity.
\marginpar{Need to tidy up method for controlling this in code}

It should be noted that despite the ad-hoc nature of the scaling, some
scaling must be done since on a lat-lon grid the converging meridians
make it very unlikely that a stable viscosity parameter exists across
the entire model domain.

The Laplacian viscosity coefficient, $A_h$ ({\bf viscAh}), has units
of $m^2 s^{-1}$. The bi-harmonic viscosity coefficient, $A_4$ ({\bf
viscA4}), has units of $m^4 s^{-1}$.

\fbox{ \begin{minipage}{4.75in}
{\em S/R MOM\_U\_XVISCFLUX} ({\em mom\_u\_xviscflux.F})

{\em S/R MOM\_U\_YVISCFLUX} ({\em mom\_u\_yviscflux.F})

{\em S/R MOM\_V\_XVISCFLUX} ({\em mom\_v\_xviscflux.F})

{\em S/R MOM\_V\_YVISCFLUX} ({\em mom\_v\_yviscflux.F})

$\tau_{11}$, $\tau_{12}$, $\tau_{22}$, $\tau_{22}$: {\bf vF}, {\bf
v4F} (local to {\em calc\_mom\_rhs.F})
\end{minipage} }

Two types of lateral boundary condition exist for the lateral viscous
terms, no-slip and free-slip.

The free-slip condition is most convenient to code since it is
equivalent to zero-stress on boundaries. Simple masking of the stress
components sets them to zero. The fractional open stress is properly
handled using the lopped cells.

The no-slip condition defines the normal gradient of a tangential flow
such that the flow is zero on the boundary. Rather than modify the
stresses by using complicated functions of the masks and ``ghost''
points (see \cite{adcroft:98}) we add the boundary stresses as
an additional source term in cells next to solid boundaries. This has
the advantage of being able to cope with ``thin walls'' and also makes
the interior stress calculation (code) independent of the boundary
conditions. The ``body'' force takes the form:
\begin{eqnarray}
G_u^{side-drag} & = &
\frac{4}{\Delta z_f} \overline{ (1-h_\zeta) \frac{\Delta x_v}{\Delta y_u} }^j
\left( A_h c_{12\Delta}(\varphi) u - A_4 c_{12\Delta^2}(\varphi) \nabla^2 u \right)
\\
G_v^{side-drag} & = &
\frac{4}{\Delta z_f} \overline{ (1-h_\zeta) \frac{\Delta y_u}{\Delta x_v} }^i
\left( A_h c_{21\Delta}(\varphi) v - A_4 c_{21\Delta^2}(\varphi) \nabla^2 v \right)
\end{eqnarray}

In fact, the above discretization is not quite complete because it
assumes that the bathymetry at velocity points is deeper than at
neighboring vorticity points, e.g. $1-h_w < 1-h_\zeta$

\fbox{ \begin{minipage}{4.75in}
{\em S/R MOM\_U\_SIDEDRAG} ({\em mom\_u\_sidedrag.F})

{\em S/R MOM\_V\_SIDEDRAG} ({\em mom\_v\_sidedrag.F})

$G_u^{side-drag}$, $G_v^{side-drag}$: {\bf vF} (local to {\em calc\_mom\_rhs.F})
\end{minipage} }


\subsection{Vertical dissipation}

Vertical viscosity terms are discretized with only partial adherence
to the variable grid lengths introduced by the finite volume
formulation. This reduces the formal accuracy of these terms to just
first order but only next to boundaries; exactly where other terms
appear such as linear and quadratic bottom drag.
\begin{eqnarray}
G_u^{v-diss} & = &
\frac{1}{\Delta r_f h_w} \delta_k \tau_{13} \\
G_v^{v-diss} & = &
\frac{1}{\Delta r_f h_s} \delta_k \tau_{23} \\
G_w^{v-diss} & = & \epsilon_{nh}
\frac{1}{\Delta r_f h_d} \delta_k \tau_{33}
\end{eqnarray}
represents the general discrete form of the vertical dissipation terms.

In the interior the vertical stresses are discretized:
\begin{eqnarray}
\tau_{13} & = & A_v \frac{1}{\Delta r_c} \delta_k u \\
\tau_{23} & = & A_v \frac{1}{\Delta r_c} \delta_k v \\
\tau_{33} & = & A_v \frac{1}{\Delta r_f} \delta_k w
\end{eqnarray}
It should be noted that in the non-hydrostatic form, the stress tensor
is even less consistent than for the hydrostatic (see
\cite{wajsowicz:93}). It is well known how to do this properly (see 
\cite{griffies:00}) and is on the list of to-do's.

\fbox{ \begin{minipage}{4.75in}
{\em S/R MOM\_U\_RVISCLFUX} ({\em mom\_u\_rviscflux.F})

{\em S/R MOM\_V\_RVISCLFUX} ({\em mom\_v\_rviscflux.F})

$\tau_{13}$: {\bf urf} (local to {\em calc\_mom\_rhs.F})

$\tau_{23}$: {\bf vrf} (local to {\em calc\_mom\_rhs.F})
\end{minipage} }


As for the lateral viscous terms, the free-slip condition is
equivalent to simply setting the stress to zero on boundaries.  The
no-slip condition is implemented as an additional term acting on top
of the interior and free-slip stresses. Bottom drag represents
additional friction, in addition to that imposed by the no-slip
condition at the bottom. The drag is cast as a stress expressed as a
linear or quadratic function of the mean flow in the layer above the
topography:
\begin{eqnarray}
\tau_{13}^{bottom-drag} & = &
\left(
2 A_v \frac{1}{\Delta r_c}
+ r_b
+ C_d \sqrt{ \overline{2 KE}^i }
\right) u \\
\tau_{23}^{bottom-drag} & = &
\left(
2 A_v \frac{1}{\Delta r_c}
+ r_b
+ C_d \sqrt{ \overline{2 KE}^j }
\right) v
\end{eqnarray}
where these terms are only evaluated immediately above topography.
$r_b$ ({\bf bottomDragLinear}) has units of $m s^{-1}$ and a typical value
of the order 0.0002 $m s^{-1}$. $C_d$ ({\bf bottomDragQuadratic}) is
dimensionless with typical values in the range 0.001--0.003.

\fbox{ \begin{minipage}{4.75in}
{\em S/R MOM\_U\_BOTTOMDRAG} ({\em mom\_u\_bottomdrag.F})

{\em S/R MOM\_V\_BOTTOMDRAG} ({\em mom\_v\_bottomdrag.F})

$\tau_{13}^{bottom-drag}$, $\tau_{23}^{bottom-drag}$: {\bf vf} (local to {\em calc\_mom\_rhs.F})
\end{minipage} }

\subsection{Derivation of discrete energy conservation}

These discrete equations conserve kinetic plus potential energy using the
following definitions:
\begin{equation}
KE = \frac{1}{2} \left( \overline{ u^2 }^i + \overline{ v^2 }^j +
\epsilon_{nh} \overline{ w^2 }^k \right)
\end{equation}

1	edhill	1.8	% $Header: /u/gcmpack/manual/part2/mom_fluxform.tex,v 1.7 2001/11/13 18:15:26 adcroft Exp $
2	adcroft	1.2	% $Name: $
3	adcroft	1.1
4			\section{Flux-form momentum equations}
5	cnh	1.5	\label{sec:flux-form_momentum_eqautions}
6	edhill	1.8	\begin{rawhtml}
7			<!-- CMIREDIR:flux-form_momentum_eqautions: -->
8			\end{rawhtml}
9	adcroft	1.1
10			The original finite volume model was based on the Eulerian flux form
11			momentum equations. This is the default though the vector invariant
12			form is optionally available (and recommended in some cases).
13
14			The ``G's'' (our colloquial name for all terms on rhs!) are broken
15			into the various advective, Coriolis, horizontal dissipation, vertical
16			dissipation and metric forces:
17			\marginpar{$G_u$: {\bf Gu} }
18			\marginpar{$G_v$: {\bf Gv} }
19			\marginpar{$G_w$: {\bf Gw} }
20			\begin{eqnarray}
21			G_u & = & G_u^{adv} + G_u^{cor} + G_u^{h-diss} + G_u^{v-diss} +
22	adcroft	1.2	G_u^{metric} + G_u^{nh-metric} \label{eq:gsplit_momu} \\
23	adcroft	1.1	G_v & = & G_v^{adv} + G_v^{cor} + G_v^{h-diss} + G_v^{v-diss} +
24	adcroft	1.2	G_v^{metric} + G_v^{nh-metric} \label{eq:gsplit_momv} \\
25	adcroft	1.1	G_w & = & G_w^{adv} + G_w^{cor} + G_w^{h-diss} + G_w^{v-diss} +
26	adcroft	1.2	G_w^{metric} + G_w^{nh-metric} \label{eq:gsplit_momw}
27	adcroft	1.1	\end{eqnarray}
28			In the hydrostatic limit, $G_w=0$ and $\epsilon_{nh}=0$, reducing the
29			vertical momentum to hydrostatic balance.
30
31			These terms are calculated in routines called from subroutine {\em
32			CALC\_MOM\_RHS} a collected into the global arrays {\bf Gu}, {\bf Gv},
33			and {\bf Gw}.
34
35			\fbox{ \begin{minipage}{4.75in}
36			{\em S/R CALC\_MOM\_RHS} ({\em pkg/mom\_fluxform/calc\_mom\_rhs.F})
37
38			$G_u$: {\bf Gu} ({\em DYNVARS.h})
39
40			$G_v$: {\bf Gv} ({\em DYNVARS.h})
41
42			$G_w$: {\bf Gw} ({\em DYNVARS.h})
43			\end{minipage} }
44
45
46			\subsection{Advection of momentum}
47
48			The advective operator is second order accurate in space:
49			\begin{eqnarray}
50			{\cal A}_w \Delta r_f h_w G_u^{adv} & = &
51			\delta_i \overline{ U }^i \overline{ u }^i
52			+ \delta_j \overline{ V }^i \overline{ u }^j
53	adcroft	1.2	+ \delta_k \overline{ W }^i \overline{ u }^k \label{eq:discrete-momadvu} \\
54	adcroft	1.1	{\cal A}_s \Delta r_f h_s G_v^{adv} & = &
55			\delta_i \overline{ U }^j \overline{ v }^i
56			+ \delta_j \overline{ V }^j \overline{ v }^j
57	adcroft	1.2	+ \delta_k \overline{ W }^j \overline{ v }^k \label{eq:discrete-momadvv} \\
58	adcroft	1.1	{\cal A}_c \Delta r_c G_w^{adv} & = &
59			\delta_i \overline{ U }^k \overline{ w }^i
60			+ \delta_j \overline{ V }^k \overline{ w }^j
61	adcroft	1.2	+ \delta_k \overline{ W }^k \overline{ w }^k \label{eq:discrete-momadvw}
62	adcroft	1.1	\end{eqnarray}
63			and because of the flux form does not contribute to the global budget
64			of linear momentum. The quantities $U$, $V$ and $W$ are volume fluxes
65			defined:
66			\marginpar{$U$: {\bf uTrans} }
67			\marginpar{$V$: {\bf vTrans} }
68			\marginpar{$W$: {\bf rTrans} }
69			\begin{eqnarray}
70	adcroft	1.2	U & = & \Delta y_g \Delta r_f h_w u \label{eq:utrans} \\
71			V & = & \Delta x_g \Delta r_f h_s v \label{eq:vtrans} \\
72			W & = & {\cal A}_c w \label{eq:rtrans}
73	adcroft	1.1	\end{eqnarray}
74			The advection of momentum takes the same form as the advection of
75			tracers but by a translated advective flow. Consequently, the
76			conservation of second moments, derived for tracers later, applies to
77			$u^2$ and $v^2$ and $w^2$ so that advection of momentum correctly
78			conserves kinetic energy.
79
80			\fbox{ \begin{minipage}{4.75in}
81			{\em S/R MOM\_U\_ADV\_UU} ({\em mom\_u\_adv\_uu.F})
82
83			{\em S/R MOM\_U\_ADV\_VU} ({\em mom\_u\_adv\_vu.F})
84
85			{\em S/R MOM\_U\_ADV\_WU} ({\em mom\_u\_adv\_wu.F})
86
87			{\em S/R MOM\_U\_ADV\_UV} ({\em mom\_u\_adv\_uv.F})
88
89			{\em S/R MOM\_U\_ADV\_VV} ({\em mom\_u\_adv\_vv.F})
90
91			{\em S/R MOM\_U\_ADV\_WV} ({\em mom\_u\_adv\_wv.F})
92
93			$uu$, $uv$, $vu$, $vv$: {\bf aF} (local to {\em calc\_mom\_rhs.F})
94			\end{minipage} }
95
96
97
98			\subsection{Coriolis terms}
99
100			The ``pure C grid'' Coriolis terms (i.e. in absence of C-D scheme) are
101			discretized:
102			\begin{eqnarray}
103			{\cal A}_w \Delta r_f h_w G_u^{Cor} & = &
104			\overline{ f {\cal A}_c \Delta r_f h_c \overline{ v }^j }^i
105			- \epsilon_{nh} \overline{ f' {\cal A}_c \Delta r_f h_c \overline{ w }^k }^i \\
106			{\cal A}_s \Delta r_f h_s G_v^{Cor} & = &
107			- \overline{ f {\cal A}_c \Delta r_f h_c \overline{ u }^i }^j \\
108			{\cal A}_c \Delta r_c G_w^{Cor} & = &
109			\epsilon_{nh} \overline{ f' {\cal A}_c \Delta r_f h_c \overline{ u }^i }^k
110			\end{eqnarray}
111			where the Coriolis parameters $f$ and $f'$ are defined:
112			\begin{eqnarray}
113	cnh	1.4	f & = & 2 \Omega \sin{\varphi} \\
114			f' & = & 2 \Omega \cos{\varphi}
115	adcroft	1.1	\end{eqnarray}
116	cnh	1.4	where $\varphi$ is geographic latitude when using spherical geometry,
117	adcroft	1.2	otherwise the $\beta$-plane definition is used:
118	adcroft	1.1	\begin{eqnarray}
119			f & = & f_o + \beta y \\
120			f' & = & 0
121			\end{eqnarray}
122
123			This discretization globally conserves kinetic energy. It should be
124			noted that despite the use of this discretization in former
125			publications, all calculations to date have used the following
126			different discretization:
127			\begin{eqnarray}
128			G_u^{Cor} & = &
129			f_u \overline{ v }^{ji}
130			- \epsilon_{nh} f_u' \overline{ w }^{ik} \\
131			G_v^{Cor} & = &
132			- f_v \overline{ u }^{ij} \\
133			G_w^{Cor} & = &
134			\epsilon_{nh} f_w' \overline{ u }^{ik}
135			\end{eqnarray}
136			\marginpar{Need to change the default in code to match this}
137			where the subscripts on $f$ and $f'$ indicate evaluation of the
138			Coriolis parameters at the appropriate points in space. The above
139	adcroft	1.3	discretization does {\em not} conserve anything, especially energy and
140			for historical reasons is the default for the code. A
141			flag controls this discretization: set run-time logical {\bf
142			useEnergyConservingCoriolis} to {\em true} which otherwise defaults to
143			{\em false}.
144	adcroft	1.1
145			\fbox{ \begin{minipage}{4.75in}
146			{\em S/R MOM\_CDSCHEME} ({\em mom\_cdscheme.F})
147
148			{\em S/R MOM\_U\_CORIOLIS} ({\em mom\_u\_coriolis.F})
149
150			{\em S/R MOM\_V\_CORIOLIS} ({\em mom\_v\_coriolis.F})
151
152			$G_u^{Cor}$, $G_v^{Cor}$: {\bf cF} (local to {\em calc\_mom\_rhs.F})
153			\end{minipage} }
154
155
156			\subsection{Curvature metric terms}
157
158			The most commonly used coordinate system on the sphere is the
159	cnh	1.4	geographic system $(\lambda,\varphi)$. The curvilinear nature of these
160	adcroft	1.1	coordinates on the sphere lead to some ``metric'' terms in the
161			component momentum equations. Under the thin-atmosphere and
162			hydrostatic approximations these terms are discretized:
163			\begin{eqnarray}
164			{\cal A}_w \Delta r_f h_w G_u^{metric} & = &
165	cnh	1.4	\overline{ \frac{ \overline{u}^i }{a} \tan{\varphi} {\cal A}_c \Delta r_f h_c \overline{ v }^j }^i \\
166	adcroft	1.1	{\cal A}_s \Delta r_f h_s G_v^{metric} & = &
167	cnh	1.4	- \overline{ \frac{ \overline{u}^i }{a} \tan{\varphi} {\cal A}_c \Delta r_f h_c \overline{ u }^i }^j \\
168	adcroft	1.1	G_w^{metric} & = & 0
169			\end{eqnarray}
170			where $a$ is the radius of the planet (sphericity is assumed) or the
171			radial distance of the particle (i.e. a function of height). It is
172			easy to see that this discretization satisfies all the properties of
173			the discrete Coriolis terms since the metric factor $\frac{u}{a}
174	cnh	1.4	\tan{\varphi}$ can be viewed as a modification of the vertical Coriolis
175			parameter: $f \rightarrow f+\frac{u}{a} \tan{\varphi}$.
176	adcroft	1.1
177			However, as for the Coriolis terms, a non-energy conserving form has
178			exclusively been used to date:
179			\begin{eqnarray}
180	cnh	1.4	G_u^{metric} & = & \frac{u \overline{v}^{ij} }{a} \tan{\varphi} \\
181			G_v^{metric} & = & \frac{ \overline{u}^{ij} \overline{u}^{ij}}{a} \tan{\varphi}
182	adcroft	1.1	\end{eqnarray}
183	cnh	1.4	where $\tan{\varphi}$ is evaluated at the $u$ and $v$ points
184	adcroft	1.1	respectively.
185
186			\fbox{ \begin{minipage}{4.75in}
187			{\em S/R MOM\_U\_METRIC\_SPHERE} ({\em mom\_u\_metric\_sphere.F})
188
189			{\em S/R MOM\_V\_METRIC\_SPHERE} ({\em mom\_v\_metric\_sphere.F})
190
191			$G_u^{metric}$, $G_v^{metric}$: {\bf mT} (local to {\em calc\_mom\_rhs.F})
192			\end{minipage} }
193
194
195
196			\subsection{Non-hydrostatic metric terms}
197
198			For the non-hydrostatic equations, dropping the thin-atmosphere
199			approximation re-introduces metric terms involving $w$ and are
200	cnh	1.6	required to conserve angular momentum:
201	adcroft	1.1	\begin{eqnarray}
202			{\cal A}_w \Delta r_f h_w G_u^{metric} & = &
203			- \overline{ \frac{ \overline{u}^i \overline{w}^k }{a} {\cal A}_c \Delta r_f h_c }^i \\
204			{\cal A}_s \Delta r_f h_s G_v^{metric} & = &
205			- \overline{ \frac{ \overline{v}^j \overline{w}^k }{a} {\cal A}_c \Delta r_f h_c}^j \\
206			{\cal A}_c \Delta r_c G_w^{metric} & = &
207			\overline{ \frac{ {\overline{u}^i}^2 + {\overline{v}^j}^2}{a} {\cal A}_c \Delta r_f h_c }^k
208			\end{eqnarray}
209
210			Because we are always consistent, even if consistently wrong, we have,
211			in the past, used a different discretization in the model which is:
212			\begin{eqnarray}
213			G_u^{metric} & = &
214			- \frac{u}{a} \overline{w}^{ik} \\
215			G_v^{metric} & = &
216			- \frac{v}{a} \overline{w}^{jk} \\
217			G_w^{metric} & = &
218			\frac{1}{a} ( {\overline{u}^{ik}}^2 + {\overline{v}^{jk}}^2 )
219			\end{eqnarray}
220
221			\fbox{ \begin{minipage}{4.75in}
222			{\em S/R MOM\_U\_METRIC\_NH} ({\em mom\_u\_metric\_nh.F})
223
224			{\em S/R MOM\_V\_METRIC\_NH} ({\em mom\_v\_metric\_nh.F})
225
226			$G_u^{metric}$, $G_v^{metric}$: {\bf mT} (local to {\em calc\_mom\_rhs.F})
227			\end{minipage} }
228
229
230			\subsection{Lateral dissipation}
231
232			Historically, we have represented the SGS Reynolds stresses as simply
233			down gradient momentum fluxes, ignoring constraints on the stress
234			tensor such as symmetry.
235			\begin{eqnarray}
236			{\cal A}_w \Delta r_f h_w G_u^{h-diss} & = &
237			\delta_i \Delta y_f \Delta r_f h_c \tau_{11}
238			+ \delta_j \Delta x_v \Delta r_f h_\zeta \tau_{12} \\
239			{\cal A}_s \Delta r_f h_s G_v^{h-diss} & = &
240			\delta_i \Delta y_u \Delta r_f h_\zeta \tau_{21}
241			+ \delta_j \Delta x_f \Delta r_f h_c \tau_{22}
242			\end{eqnarray}
243			\marginpar{Check signs of stress definitions}
244
245			The lateral viscous stresses are discretized:
246			\begin{eqnarray}
247	cnh	1.4	\tau_{11} & = & A_h c_{11\Delta}(\varphi) \frac{1}{\Delta x_f} \delta_i u
248			-A_4 c_{11\Delta^2}(\varphi) \frac{1}{\Delta x_f} \delta_i \nabla^2 u \\
249			\tau_{12} & = & A_h c_{12\Delta}(\varphi) \frac{1}{\Delta y_u} \delta_j u
250			-A_4 c_{12\Delta^2}(\varphi)\frac{1}{\Delta y_u} \delta_j \nabla^2 u \\
251			\tau_{21} & = & A_h c_{21\Delta}(\varphi) \frac{1}{\Delta x_v} \delta_i v
252			-A_4 c_{21\Delta^2}(\varphi) \frac{1}{\Delta x_v} \delta_i \nabla^2 v \\
253			\tau_{22} & = & A_h c_{22\Delta}(\varphi) \frac{1}{\Delta y_f} \delta_j v
254			-A_4 c_{22\Delta^2}(\varphi) \frac{1}{\Delta y_f} \delta_j \nabla^2 v
255	adcroft	1.1	\end{eqnarray}
256	cnh	1.4	where the non-dimensional factors $c_{lm\Delta^n}(\varphi), \{l,m,n\} \in
257	adcroft	1.1	\{1,2\}$ define the ``cosine'' scaling with latitude which can be
258			applied in various ad-hoc ways. For instance, $c_{11\Delta} =
259	cnh	1.4	c_{21\Delta} = (\cos{\varphi})^{3/2}$, $c_{12\Delta}=c_{22\Delta}=0$ would
260	adcroft	1.1	represent the an-isotropic cosine scaling typically used on the
261			``lat-lon'' grid for Laplacian viscosity.
262			\marginpar{Need to tidy up method for controlling this in code}
263
264	cnh	1.6	It should be noted that despite the ad-hoc nature of the scaling, some
265	adcroft	1.1	scaling must be done since on a lat-lon grid the converging meridians
266			make it very unlikely that a stable viscosity parameter exists across
267			the entire model domain.
268
269			The Laplacian viscosity coefficient, $A_h$ ({\bf viscAh}), has units
270			of $m^2 s^{-1}$. The bi-harmonic viscosity coefficient, $A_4$ ({\bf
271			viscA4}), has units of $m^4 s^{-1}$.
272
273			\fbox{ \begin{minipage}{4.75in}
274			{\em S/R MOM\_U\_XVISCFLUX} ({\em mom\_u\_xviscflux.F})
275
276			{\em S/R MOM\_U\_YVISCFLUX} ({\em mom\_u\_yviscflux.F})
277
278			{\em S/R MOM\_V\_XVISCFLUX} ({\em mom\_v\_xviscflux.F})
279
280			{\em S/R MOM\_V\_YVISCFLUX} ({\em mom\_v\_yviscflux.F})
281
282			$\tau_{11}$, $\tau_{12}$, $\tau_{22}$, $\tau_{22}$: {\bf vF}, {\bf
283			v4F} (local to {\em calc\_mom\_rhs.F})
284			\end{minipage} }
285
286			Two types of lateral boundary condition exist for the lateral viscous
287			terms, no-slip and free-slip.
288
289			The free-slip condition is most convenient to code since it is
290			equivalent to zero-stress on boundaries. Simple masking of the stress
291			components sets them to zero. The fractional open stress is properly
292			handled using the lopped cells.
293
294			The no-slip condition defines the normal gradient of a tangential flow
295			such that the flow is zero on the boundary. Rather than modify the
296			stresses by using complicated functions of the masks and ``ghost''
297	adcroft	1.7	points (see \cite{adcroft:98}) we add the boundary stresses as
298	adcroft	1.1	an additional source term in cells next to solid boundaries. This has
299			the advantage of being able to cope with ``thin walls'' and also makes
300			the interior stress calculation (code) independent of the boundary
301			conditions. The ``body'' force takes the form:
302			\begin{eqnarray}
303			G_u^{side-drag} & = &
304			\frac{4}{\Delta z_f} \overline{ (1-h_\zeta) \frac{\Delta x_v}{\Delta y_u} }^j
305	cnh	1.4	\left( A_h c_{12\Delta}(\varphi) u - A_4 c_{12\Delta^2}(\varphi) \nabla^2 u \right)
306	adcroft	1.1	\\
307			G_v^{side-drag} & = &
308			\frac{4}{\Delta z_f} \overline{ (1-h_\zeta) \frac{\Delta y_u}{\Delta x_v} }^i
309	cnh	1.4	\left( A_h c_{21\Delta}(\varphi) v - A_4 c_{21\Delta^2}(\varphi) \nabla^2 v \right)
310	adcroft	1.1	\end{eqnarray}
311
312			In fact, the above discretization is not quite complete because it
313			assumes that the bathymetry at velocity points is deeper than at
314	cnh	1.6	neighboring vorticity points, e.g. $1-h_w < 1-h_\zeta$
315	adcroft	1.1
316			\fbox{ \begin{minipage}{4.75in}
317			{\em S/R MOM\_U\_SIDEDRAG} ({\em mom\_u\_sidedrag.F})
318
319			{\em S/R MOM\_V\_SIDEDRAG} ({\em mom\_v\_sidedrag.F})
320
321			$G_u^{side-drag}$, $G_v^{side-drag}$: {\bf vF} (local to {\em calc\_mom\_rhs.F})
322			\end{minipage} }
323
324
325			\subsection{Vertical dissipation}
326
327			Vertical viscosity terms are discretized with only partial adherence
328			to the variable grid lengths introduced by the finite volume
329			formulation. This reduces the formal accuracy of these terms to just
330			first order but only next to boundaries; exactly where other terms
331	cnh	1.6	appear such as linear and quadratic bottom drag.
332	adcroft	1.1	\begin{eqnarray}
333			G_u^{v-diss} & = &
334			\frac{1}{\Delta r_f h_w} \delta_k \tau_{13} \\
335			G_v^{v-diss} & = &
336			\frac{1}{\Delta r_f h_s} \delta_k \tau_{23} \\
337			G_w^{v-diss} & = & \epsilon_{nh}
338			\frac{1}{\Delta r_f h_d} \delta_k \tau_{33}
339			\end{eqnarray}
340			represents the general discrete form of the vertical dissipation terms.
341
342			In the interior the vertical stresses are discretized:
343			\begin{eqnarray}
344			\tau_{13} & = & A_v \frac{1}{\Delta r_c} \delta_k u \\
345			\tau_{23} & = & A_v \frac{1}{\Delta r_c} \delta_k v \\
346			\tau_{33} & = & A_v \frac{1}{\Delta r_f} \delta_k w
347			\end{eqnarray}
348			It should be noted that in the non-hydrostatic form, the stress tensor
349	adcroft	1.7	is even less consistent than for the hydrostatic (see
350			\cite{wajsowicz:93}). It is well known how to do this properly (see
351			\cite{griffies:00}) and is on the list of to-do's.
352	adcroft	1.1
353			\fbox{ \begin{minipage}{4.75in}
354			{\em S/R MOM\_U\_RVISCLFUX} ({\em mom\_u\_rviscflux.F})
355
356			{\em S/R MOM\_V\_RVISCLFUX} ({\em mom\_v\_rviscflux.F})
357
358			$\tau_{13}$: {\bf urf} (local to {\em calc\_mom\_rhs.F})
359
360			$\tau_{23}$: {\bf vrf} (local to {\em calc\_mom\_rhs.F})
361			\end{minipage} }
362
363
364			As for the lateral viscous terms, the free-slip condition is
365			equivalent to simply setting the stress to zero on boundaries. The
366			no-slip condition is implemented as an additional term acting on top
367			of the interior and free-slip stresses. Bottom drag represents
368			additional friction, in addition to that imposed by the no-slip
369			condition at the bottom. The drag is cast as a stress expressed as a
370			linear or quadratic function of the mean flow in the layer above the
371			topography:
372			\begin{eqnarray}
373			\tau_{13}^{bottom-drag} & = &
374			\left(
375			2 A_v \frac{1}{\Delta r_c}
376			+ r_b
377			+ C_d \sqrt{ \overline{2 KE}^i }
378			\right) u \\
379			\tau_{23}^{bottom-drag} & = &
380			\left(
381			2 A_v \frac{1}{\Delta r_c}
382			+ r_b
383			+ C_d \sqrt{ \overline{2 KE}^j }
384			\right) v
385			\end{eqnarray}
386			where these terms are only evaluated immediately above topography.
387			$r_b$ ({\bf bottomDragLinear}) has units of $m s^{-1}$ and a typical value
388			of the order 0.0002 $m s^{-1}$. $C_d$ ({\bf bottomDragQuadratic}) is
389			dimensionless with typical values in the range 0.001--0.003.
390
391			\fbox{ \begin{minipage}{4.75in}
392			{\em S/R MOM\_U\_BOTTOMDRAG} ({\em mom\_u\_bottomdrag.F})
393
394			{\em S/R MOM\_V\_BOTTOMDRAG} ({\em mom\_v\_bottomdrag.F})
395
396			$\tau_{13}^{bottom-drag}$, $\tau_{23}^{bottom-drag}$: {\bf vf} (local to {\em calc\_mom\_rhs.F})
397			\end{minipage} }
398
399			\subsection{Derivation of discrete energy conservation}
400
401			These discrete equations conserve kinetic plus potential energy using the
402			following definitions:
403			\begin{equation}
404			KE = \frac{1}{2} \left( \overline{ u^2 }^i + \overline{ v^2 }^j +
405			\epsilon_{nh} \overline{ w^2 }^k \right)
406			\end{equation}
407