s_algorithm/text/mom_vecinv.tex

% $Header: /u/gcmpack/manual/s_algorithm/text/mom_vecinv.tex,v 1.6 2006/06/26 01:03:47 jmc Exp $
% $Name:  $

\section{Vector invariant momentum equations}
\label{sec:vect-inv_momentum_equations}
\begin{rawhtml}
<!-- CMIREDIR:vector_invariant_momentum_eqautions: -->
\end{rawhtml}

The finite volume method lends itself to describing the continuity and
tracer equations in curvilinear coordinate systems. However, in
curvilinear coordinates many new metric terms appear in the momentum
equations (written in Lagrangian or flux-form) making generalization
far from elegant. Fortunately, an alternative form of the equations,
the vector invariant equations are exactly that; invariant under
coordinate transformations so that they can be applied uniformly in
any orthogonal curvilinear coordinate system such as spherical
coordinates, boundary following or the conformal spherical cube
system.

The non-hydrostatic vector invariant equations read:
\begin{equation}
\partial_t \vec{v} + ( 2\vec{\Omega} + \vec{\zeta}) \wedge \vec{v}
- b \hat{r}
+ \vec{\nabla} B = \vec{\nabla} \cdot \vec{\bf \tau}
\end{equation}
which describe motions in any orthogonal curvilinear coordinate
system. Here, $B$ is the Bernoulli function and $\vec{\zeta}=\nabla
\wedge \vec{v}$ is the vorticity vector. We can take advantage of the
elegance of these equations when discretizing them and use the
discrete definitions of the grad, curl and divergence operators to
satisfy constraints. We can also consider the analogy to forming
derived equations, such as the vorticity equation, and examine how the
discretization can be adjusted to give suitable vorticity advection
among other things.

The underlying algorithm is the same as for the flux form
equations. All that has changed is the contents of the ``G's''. For
the time-being, only the hydrostatic terms have been coded but we will
indicate the points where non-hydrostatic contributions will enter:
\begin{eqnarray}
G_u & = & G_u^{fv} + G_u^{\zeta_3 v} + G_u^{\zeta_2 w} + G_u^{\partial_x B}
+ G_u^{\partial_z \tau^x} + G_u^{h-dissip} + G_u^{v-dissip} \\
G_v & = & G_v^{fu} + G_v^{\zeta_3 u} + G_v^{\zeta_1 w} + G_v^{\partial_y B}
+ G_v^{\partial_z \tau^y} + G_v^{h-dissip} + G_v^{v-dissip} \\
G_w & = & G_w^{fu} + G_w^{\zeta_1 v} + G_w^{\zeta_2 u} + G_w^{\partial_z B}
+ G_w^{h-dissip} + G_w^{v-dissip}
\end{eqnarray}

\fbox{ \begin{minipage}{4.75in}
{\em S/R MOM\_VECINV} ({\em pkg/mom\_vecinv/mom\_vecinv.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{Relative vorticity}

The vertical component of relative vorticity is explicitly calculated
and use in the discretization. The particular form is crucial for
numerical stability; alternative definitions break the conservation
properties of the discrete equations.

Relative vorticity is defined:
\begin{equation}
\zeta_3 = \frac{\Gamma}{A_\zeta}
= \frac{1}{{\cal A}_\zeta} ( \delta_i \Delta y_c v - \delta_j \Delta x_c u )
\end{equation}
where ${\cal A}_\zeta$ is the area of the vorticity cell presented in
the vertical and $\Gamma$ is the circulation about that cell.

\fbox{ \begin{minipage}{4.75in}
{\em S/R MOM\_VI\_CALC\_RELVORT3} ({\em mom\_vi\_calc\_relvort3.F})

$\zeta_3$: {\bf vort3} (local to {\em mom\_vecinv.F})
\end{minipage} }


\subsection{Kinetic energy}

The kinetic energy, denoted $KE$, is defined:
\begin{equation}
KE = \frac{1}{2} ( \overline{ u^2 }^i + \overline{ v^2 }^j 
+ \epsilon_{nh} \overline{ w^2 }^k )
\end{equation}

\fbox{ \begin{minipage}{4.75in}
{\em S/R MOM\_VI\_CALC\_KE} ({\em mom\_vi\_calc\_ke.F})

$KE$: {\bf KE} (local to {\em mom\_vecinv.F})
\end{minipage} }


\subsection{Coriolis terms}

The potential enstrophy conserving form of the linear Coriolis terms
are written:
\begin{eqnarray}
G_u^{fv} & = &
\frac{1}{\Delta x_c}
\overline{ \frac{f}{h_\zeta} }^j \overline{ \overline{ \Delta x_g h_s v }^j }^i \\
G_v^{fu} & = & -
\frac{1}{\Delta y_c}
\overline{ \frac{f}{h_\zeta} }^i \overline{ \overline{ \Delta y_g h_w u }^i }^j
\end{eqnarray}
Here, the Coriolis parameter $f$ is defined at vorticity (corner)
points.
\marginpar{$f$: {\bf fCoriG}}
\marginpar{$h_\zeta$: {\bf hFacZ}}

The potential enstrophy conserving form of the non-linear Coriolis
terms are written:
\begin{eqnarray}
G_u^{\zeta_3 v} & = &
\frac{1}{\Delta x_c}
\overline{ \frac{\zeta_3}{h_\zeta} }^j \overline{ \overline{ \Delta x_g h_s v }^j }^i \\
G_v^{\zeta_3 u} & = & -
\frac{1}{\Delta y_c}
\overline{ \frac{\zeta_3}{h_\zeta} }^i \overline{ \overline{ \Delta y_g h_w u }^i }^j
\end{eqnarray}
\marginpar{$\zeta_3$: {\bf vort3}}

The Coriolis terms can also be evaluated together and expressed in
terms of absolute vorticity $f+\zeta_3$. The potential enstrophy
conserving form using the absolute vorticity is written:
\begin{eqnarray}
G_u^{fv} + G_u^{\zeta_3 v} & = &
\frac{1}{\Delta x_c}
\overline{ \frac{f + \zeta_3}{h_\zeta} }^j \overline{ \overline{ \Delta x_g h_s v }^j }^i \\
G_v^{fu} + G_v^{\zeta_3 u} & = & -
\frac{1}{\Delta y_c}
\overline{ \frac{f + \zeta_3}{h_\zeta} }^i \overline{ \overline{ \Delta y_g h_w u }^i }^j
\end{eqnarray}

\marginpar{Run-time control needs to be added for these options} The
distinction between using absolute vorticity or relative vorticity is
useful when constructing higher order advection schemes; monotone
advection of relative vorticity behaves differently to monotone
advection of absolute vorticity. Currently the choice of
relative/absolute vorticity, centered/upwind/high order advection is
available only through commented subroutine calls.

\fbox{ \begin{minipage}{4.75in}
{\em S/R MOM\_VI\_CORIOLIS} ({\em mom\_vi\_coriolis.F})

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

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

$G_u^{fv}$, $G_u^{\zeta_3 v}$: {\bf uCf} (local to {\em mom\_vecinv.F})

$G_v^{fu}$, $G_v^{\zeta_3 u}$: {\bf vCf} (local to {\em mom\_vecinv.F})
\end{minipage} }


\subsection{Shear terms}

The shear terms ($\zeta_2w$ and $\zeta_1w$) are are discretized to
guarantee that no spurious generation of kinetic energy is possible;
the horizontal gradient of Bernoulli function has to be consistent
with the vertical advection of shear:
\marginpar{N-H terms have not been tried!}
\begin{eqnarray}
G_u^{\zeta_2 w} & = &
\frac{1}{ {\cal A}_w \Delta r_f h_w } \overline{
\overline{ {\cal A}_c w }^i ( \delta_k u - \epsilon_{nh} \delta_j w )
}^k \\
G_v^{\zeta_1 w} & = &
\frac{1}{ {\cal A}_s \Delta r_f h_s } \overline{
\overline{ {\cal A}_c w }^i ( \delta_k u - \epsilon_{nh} \delta_j w )
}^k
\end{eqnarray}

\fbox{ \begin{minipage}{4.75in}
{\em S/R MOM\_VI\_U\_VERTSHEAR} ({\em mom\_vi\_u\_vertshear.F})

{\em S/R MOM\_VI\_V\_VERTSHEAR} ({\em mom\_vi\_v\_vertshear.F})

$G_u^{\zeta_2 w}$: {\bf uCf} (local to {\em mom\_vecinv.F})

$G_v^{\zeta_1 w}$: {\bf vCf} (local to {\em mom\_vecinv.F})
\end{minipage} }


\subsection{Gradient of Bernoulli function}

\begin{eqnarray}
G_u^{\partial_x B} & = &
\frac{1}{\Delta x_c} \delta_i ( \phi' + KE ) \\
G_v^{\partial_y B} & = &
\frac{1}{\Delta x_y} \delta_j ( \phi' + KE )
%G_w^{\partial_z B} & = &
%\frac{1}{\Delta r_c} h_c \delta_k ( \phi' + KE )
\end{eqnarray}

\fbox{ \begin{minipage}{4.75in}
{\em S/R MOM\_VI\_U\_GRAD\_KE} ({\em mom\_vi\_u\_grad\_ke.F})

{\em S/R MOM\_VI\_V\_GRAD\_KE} ({\em mom\_vi\_v\_grad\_ke.F})

$G_u^{\partial_x KE}$: {\bf uCf} (local to {\em mom\_vecinv.F})

$G_v^{\partial_y KE}$: {\bf vCf} (local to {\em mom\_vecinv.F})
\end{minipage} }


\subsection{Horizontal divergence}

The horizontal divergence, a complimentary quantity to relative
vorticity, is used in parameterizing the Reynolds stresses and is
discretized:
\begin{equation}
D = \frac{1}{{\cal A}_c h_c} (
  \delta_i \Delta y_g h_w u
+ \delta_j \Delta x_g h_s v )
\end{equation}

\fbox{ \begin{minipage}{4.75in}
{\em S/R MOM\_VI\_CALC\_HDIV} ({\em mom\_vi\_calc\_hdiv.F})

$D$: {\bf hDiv} (local to {\em mom\_vecinv.F})
\end{minipage} }


\subsection{Horizontal dissipation}

The following discretization of horizontal dissipation conserves
potential vorticity (thickness weighted relative vorticity) and
divergence and dissipates energy, enstrophy and divergence squared:
\begin{eqnarray}
G_u^{h-dissip} & = &
  \frac{1}{\Delta x_c} \delta_i ( A_D D - A_{D4} D^*)
- \frac{1}{\Delta y_u h_w} \delta_j h_\zeta ( A_\zeta \zeta - A_{\zeta4} \zeta^* )
\\
G_v^{h-dissip} & = &
  \frac{1}{\Delta x_v h_s} \delta_i h_\zeta ( A_\zeta \zeta - A_\zeta \zeta^* )
+ \frac{1}{\Delta y_c} \delta_j ( A_D D - A_{D4} D^* )
\end{eqnarray}
where
\begin{eqnarray}
D^* & = & \frac{1}{{\cal A}_c h_c} (
  \delta_i \Delta y_g h_w \nabla^2 u
+ \delta_j \Delta x_g h_s \nabla^2 v ) \\
\zeta^* & = & \frac{1}{{\cal A}_\zeta} (
  \delta_i \Delta y_c \nabla^2 v
- \delta_j \Delta x_c \nabla^2 u )
\end{eqnarray}

\fbox{ \begin{minipage}{4.75in}
{\em S/R MOM\_VI\_HDISSIP} ({\em mom\_vi\_hdissip.F})

$G_u^{h-dissip}$: {\bf uDiss} (local to {\em mom\_vecinv.F})

$G_v^{h-dissip}$: {\bf vDiss} (local to {\em mom\_vecinv.F})
\end{minipage} }


\subsection{Vertical dissipation}

Currently, this is exactly the same code as the flux form equations.
\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}
\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
\end{eqnarray}

\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 mom\_vecinv.F})

$\tau_{23}$: {\bf vrf} (local to {\em mom\_vecinv.F})
\end{minipage} }
1	% $Header: /u/gcmpack/manual/s_algorithm/text/mom_vecinv.tex,v 1.6 2006/06/26 01:03:47 jmc Exp $
2	% $Name: $
3
4	\section{Vector invariant momentum equations}
5	\label{sec:vect-inv_momentum_equations}
6	\begin{rawhtml}
7	<!-- CMIREDIR:vector_invariant_momentum_eqautions: -->
8	\end{rawhtml}
9
10	The finite volume method lends itself to describing the continuity and
11	tracer equations in curvilinear coordinate systems. However, in
12	curvilinear coordinates many new metric terms appear in the momentum
13	equations (written in Lagrangian or flux-form) making generalization
14	far from elegant. Fortunately, an alternative form of the equations,
15	the vector invariant equations are exactly that; invariant under
16	coordinate transformations so that they can be applied uniformly in
17	any orthogonal curvilinear coordinate system such as spherical
18	coordinates, boundary following or the conformal spherical cube
19	system.
20
21	The non-hydrostatic vector invariant equations read:
22	\begin{equation}
23	\partial_t \vec{v} + ( 2\vec{\Omega} + \vec{\zeta}) \wedge \vec{v}
24	- b \hat{r}
25	+ \vec{\nabla} B = \vec{\nabla} \cdot \vec{\bf \tau}
26	\end{equation}
27	which describe motions in any orthogonal curvilinear coordinate
28	system. Here, $B$ is the Bernoulli function and $\vec{\zeta}=\nabla
29	\wedge \vec{v}$ is the vorticity vector. We can take advantage of the
30	elegance of these equations when discretizing them and use the
31	discrete definitions of the grad, curl and divergence operators to
32	satisfy constraints. We can also consider the analogy to forming
33	derived equations, such as the vorticity equation, and examine how the
34	discretization can be adjusted to give suitable vorticity advection
35	among other things.
36
37	The underlying algorithm is the same as for the flux form
38	equations. All that has changed is the contents of the ``G's''. For
39	the time-being, only the hydrostatic terms have been coded but we will
40	indicate the points where non-hydrostatic contributions will enter:
41	\begin{eqnarray}
42	G_u & = & G_u^{fv} + G_u^{\zeta_3 v} + G_u^{\zeta_2 w} + G_u^{\partial_x B}
43	+ G_u^{\partial_z \tau^x} + G_u^{h-dissip} + G_u^{v-dissip} \\
44	G_v & = & G_v^{fu} + G_v^{\zeta_3 u} + G_v^{\zeta_1 w} + G_v^{\partial_y B}
45	+ G_v^{\partial_z \tau^y} + G_v^{h-dissip} + G_v^{v-dissip} \\
46	G_w & = & G_w^{fu} + G_w^{\zeta_1 v} + G_w^{\zeta_2 u} + G_w^{\partial_z B}
47	+ G_w^{h-dissip} + G_w^{v-dissip}
48	\end{eqnarray}
49
50	\fbox{ \begin{minipage}{4.75in}
51	{\em S/R MOM\_VECINV} ({\em pkg/mom\_vecinv/mom\_vecinv.F})
52
53	$G_u$: {\bf Gu} ({\em DYNVARS.h})
54
55	$G_v$: {\bf Gv} ({\em DYNVARS.h})
56
57	$G_w$: {\bf Gw} ({\em DYNVARS.h})
58	\end{minipage} }
59
60	\subsection{Relative vorticity}
61
62	The vertical component of relative vorticity is explicitly calculated
63	and use in the discretization. The particular form is crucial for
64	numerical stability; alternative definitions break the conservation
65	properties of the discrete equations.
66
67	Relative vorticity is defined:
68	\begin{equation}
69	\zeta_3 = \frac{\Gamma}{A_\zeta}
70	= \frac{1}{{\cal A}_\zeta} ( \delta_i \Delta y_c v - \delta_j \Delta x_c u )
71	\end{equation}
72	where ${\cal A}_\zeta$ is the area of the vorticity cell presented in
73	the vertical and $\Gamma$ is the circulation about that cell.
74
75	\fbox{ \begin{minipage}{4.75in}
76	{\em S/R MOM\_VI\_CALC\_RELVORT3} ({\em mom\_vi\_calc\_relvort3.F})
77
78	$\zeta_3$: {\bf vort3} (local to {\em mom\_vecinv.F})
79	\end{minipage} }
80
81
82	\subsection{Kinetic energy}
83
84	The kinetic energy, denoted $KE$, is defined:
85	\begin{equation}
86	KE = \frac{1}{2} ( \overline{ u^2 }^i + \overline{ v^2 }^j
87	+ \epsilon_{nh} \overline{ w^2 }^k )
88	\end{equation}
89
90	\fbox{ \begin{minipage}{4.75in}
91	{\em S/R MOM\_VI\_CALC\_KE} ({\em mom\_vi\_calc\_ke.F})
92
93	$KE$: {\bf KE} (local to {\em mom\_vecinv.F})
94	\end{minipage} }
95
96
97	\subsection{Coriolis terms}
98
99	The potential enstrophy conserving form of the linear Coriolis terms
100	are written:
101	\begin{eqnarray}
102	G_u^{fv} & = &
103	\frac{1}{\Delta x_c}
104	\overline{ \frac{f}{h_\zeta} }^j \overline{ \overline{ \Delta x_g h_s v }^j }^i \\
105	G_v^{fu} & = & -
106	\frac{1}{\Delta y_c}
107	\overline{ \frac{f}{h_\zeta} }^i \overline{ \overline{ \Delta y_g h_w u }^i }^j
108	\end{eqnarray}
109	Here, the Coriolis parameter $f$ is defined at vorticity (corner)
110	points.
111	\marginpar{$f$: {\bf fCoriG}}
112	\marginpar{$h_\zeta$: {\bf hFacZ}}
113
114	The potential enstrophy conserving form of the non-linear Coriolis
115	terms are written:
116	\begin{eqnarray}
117	G_u^{\zeta_3 v} & = &
118	\frac{1}{\Delta x_c}
119	\overline{ \frac{\zeta_3}{h_\zeta} }^j \overline{ \overline{ \Delta x_g h_s v }^j }^i \\
120	G_v^{\zeta_3 u} & = & -
121	\frac{1}{\Delta y_c}
122	\overline{ \frac{\zeta_3}{h_\zeta} }^i \overline{ \overline{ \Delta y_g h_w u }^i }^j
123	\end{eqnarray}
124	\marginpar{$\zeta_3$: {\bf vort3}}
125
126	The Coriolis terms can also be evaluated together and expressed in
127	terms of absolute vorticity $f+\zeta_3$. The potential enstrophy
128	conserving form using the absolute vorticity is written:
129	\begin{eqnarray}
130	G_u^{fv} + G_u^{\zeta_3 v} & = &
131	\frac{1}{\Delta x_c}
132	\overline{ \frac{f + \zeta_3}{h_\zeta} }^j \overline{ \overline{ \Delta x_g h_s v }^j }^i \\
133	G_v^{fu} + G_v^{\zeta_3 u} & = & -
134	\frac{1}{\Delta y_c}
135	\overline{ \frac{f + \zeta_3}{h_\zeta} }^i \overline{ \overline{ \Delta y_g h_w u }^i }^j
136	\end{eqnarray}
137
138	\marginpar{Run-time control needs to be added for these options} The
139	distinction between using absolute vorticity or relative vorticity is
140	useful when constructing higher order advection schemes; monotone
141	advection of relative vorticity behaves differently to monotone
142	advection of absolute vorticity. Currently the choice of
143	relative/absolute vorticity, centered/upwind/high order advection is
144	available only through commented subroutine calls.
145
146	\fbox{ \begin{minipage}{4.75in}
147	{\em S/R MOM\_VI\_CORIOLIS} ({\em mom\_vi\_coriolis.F})
148
149	{\em S/R MOM\_VI\_U\_CORIOLIS} ({\em mom\_vi\_u\_coriolis.F})
150
151	{\em S/R MOM\_VI\_V\_CORIOLIS} ({\em mom\_vi\_v\_coriolis.F})
152
153	$G_u^{fv}$, $G_u^{\zeta_3 v}$: {\bf uCf} (local to {\em mom\_vecinv.F})
154
155	$G_v^{fu}$, $G_v^{\zeta_3 u}$: {\bf vCf} (local to {\em mom\_vecinv.F})
156	\end{minipage} }
157
158
159	\subsection{Shear terms}
160
161	The shear terms ($\zeta_2w$ and $\zeta_1w$) are are discretized to
162	guarantee that no spurious generation of kinetic energy is possible;
163	the horizontal gradient of Bernoulli function has to be consistent
164	with the vertical advection of shear:
165	\marginpar{N-H terms have not been tried!}
166	\begin{eqnarray}
167	G_u^{\zeta_2 w} & = &
168	\frac{1}{ {\cal A}_w \Delta r_f h_w } \overline{
169	\overline{ {\cal A}_c w }^i ( \delta_k u - \epsilon_{nh} \delta_j w )
170	}^k \\
171	G_v^{\zeta_1 w} & = &
172	\frac{1}{ {\cal A}_s \Delta r_f h_s } \overline{
173	\overline{ {\cal A}_c w }^i ( \delta_k u - \epsilon_{nh} \delta_j w )
174	}^k
175	\end{eqnarray}
176
177	\fbox{ \begin{minipage}{4.75in}
178	{\em S/R MOM\_VI\_U\_VERTSHEAR} ({\em mom\_vi\_u\_vertshear.F})
179
180	{\em S/R MOM\_VI\_V\_VERTSHEAR} ({\em mom\_vi\_v\_vertshear.F})
181
182	$G_u^{\zeta_2 w}$: {\bf uCf} (local to {\em mom\_vecinv.F})
183
184	$G_v^{\zeta_1 w}$: {\bf vCf} (local to {\em mom\_vecinv.F})
185	\end{minipage} }
186
187
188
189	\subsection{Gradient of Bernoulli function}
190
191	\begin{eqnarray}
192	G_u^{\partial_x B} & = &
193	\frac{1}{\Delta x_c} \delta_i ( \phi' + KE ) \\
194	G_v^{\partial_y B} & = &
195	\frac{1}{\Delta x_y} \delta_j ( \phi' + KE )
196	%G_w^{\partial_z B} & = &
197	%\frac{1}{\Delta r_c} h_c \delta_k ( \phi' + KE )
198	\end{eqnarray}
199
200	\fbox{ \begin{minipage}{4.75in}
201	{\em S/R MOM\_VI\_U\_GRAD\_KE} ({\em mom\_vi\_u\_grad\_ke.F})
202
203	{\em S/R MOM\_VI\_V\_GRAD\_KE} ({\em mom\_vi\_v\_grad\_ke.F})
204
205	$G_u^{\partial_x KE}$: {\bf uCf} (local to {\em mom\_vecinv.F})
206
207	$G_v^{\partial_y KE}$: {\bf vCf} (local to {\em mom\_vecinv.F})
208	\end{minipage} }
209
210
211
212	\subsection{Horizontal divergence}
213
214	The horizontal divergence, a complimentary quantity to relative
215	vorticity, is used in parameterizing the Reynolds stresses and is
216	discretized:
217	\begin{equation}
218	D = \frac{1}{{\cal A}_c h_c} (
219	\delta_i \Delta y_g h_w u
220	+ \delta_j \Delta x_g h_s v )
221	\end{equation}
222
223	\fbox{ \begin{minipage}{4.75in}
224	{\em S/R MOM\_VI\_CALC\_HDIV} ({\em mom\_vi\_calc\_hdiv.F})
225
226	$D$: {\bf hDiv} (local to {\em mom\_vecinv.F})
227	\end{minipage} }
228
229
230	\subsection{Horizontal dissipation}
231
232	The following discretization of horizontal dissipation conserves
233	potential vorticity (thickness weighted relative vorticity) and
234	divergence and dissipates energy, enstrophy and divergence squared:
235	\begin{eqnarray}
236	G_u^{h-dissip} & = &
237	\frac{1}{\Delta x_c} \delta_i ( A_D D - A_{D4} D^*)
238	- \frac{1}{\Delta y_u h_w} \delta_j h_\zeta ( A_\zeta \zeta - A_{\zeta4} \zeta^* )
239	\\
240	G_v^{h-dissip} & = &
241	\frac{1}{\Delta x_v h_s} \delta_i h_\zeta ( A_\zeta \zeta - A_\zeta \zeta^* )
242	+ \frac{1}{\Delta y_c} \delta_j ( A_D D - A_{D4} D^* )
243	\end{eqnarray}
244	where
245	\begin{eqnarray}
246	D^* & = & \frac{1}{{\cal A}_c h_c} (
247	\delta_i \Delta y_g h_w \nabla^2 u
248	+ \delta_j \Delta x_g h_s \nabla^2 v ) \\
249	\zeta^* & = & \frac{1}{{\cal A}_\zeta} (
250	\delta_i \Delta y_c \nabla^2 v
251	- \delta_j \Delta x_c \nabla^2 u )
252	\end{eqnarray}
253
254	\fbox{ \begin{minipage}{4.75in}
255	{\em S/R MOM\_VI\_HDISSIP} ({\em mom\_vi\_hdissip.F})
256
257	$G_u^{h-dissip}$: {\bf uDiss} (local to {\em mom\_vecinv.F})
258
259	$G_v^{h-dissip}$: {\bf vDiss} (local to {\em mom\_vecinv.F})
260	\end{minipage} }
261
262
263	\subsection{Vertical dissipation}
264
265	Currently, this is exactly the same code as the flux form equations.
266	\begin{eqnarray}
267	G_u^{v-diss} & = &
268	\frac{1}{\Delta r_f h_w} \delta_k \tau_{13} \\
269	G_v^{v-diss} & = &
270	\frac{1}{\Delta r_f h_s} \delta_k \tau_{23}
271	\end{eqnarray}
272	represents the general discrete form of the vertical dissipation terms.
273
274	In the interior the vertical stresses are discretized:
275	\begin{eqnarray}
276	\tau_{13} & = & A_v \frac{1}{\Delta r_c} \delta_k u \\
277	\tau_{23} & = & A_v \frac{1}{\Delta r_c} \delta_k v
278	\end{eqnarray}
279
280	\fbox{ \begin{minipage}{4.75in}
281	{\em S/R MOM\_U\_RVISCLFUX} ({\em mom\_u\_rviscflux.F})
282
283	{\em S/R MOM\_V\_RVISCLFUX} ({\em mom\_v\_rviscflux.F})
284
285	$\tau_{13}$: {\bf urf} (local to {\em mom\_vecinv.F})
286
287	$\tau_{23}$: {\bf vrf} (local to {\em mom\_vecinv.F})
288	\end{minipage} }