s_algorithm/text/mom_vecinv.tex

% $Header: /u/gcmpack/manual/part2/mom_vecinv.tex,v 1.3 2001/11/13 15:20:12 adcroft Exp $
% $Name:  $

\section{Vector invariant 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 CALC\_MOM\_RHS} ({\em pkg/mom\_vecinv/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{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 calc\_mom\_rhs.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 calc\_mom\_rhs.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 calc\_mom\_rhs.F})

$G_v^{fu}$, $G_v^{\zeta_3 u}$: {\bf vCf} (local to {\em calc\_mom\_rhs.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 calc\_mom\_rhs.F})

$G_v^{\zeta_1 w}$: {\bf vCf} (local to {\em calc\_mom\_rhs.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 calc\_mom\_rhs.F})

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


\subsection{Horizontal dissipation}

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 calc\_mom\_rhs.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 calc\_mom\_rhs.F})

$G_v^{h-dissip}$: {\bf vDiss} (local to {\em calc\_mom\_rhs.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 calc\_mom\_rhs.F})

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