s_algorithm/text/mom_vecinv.tex

% $Header: /u/gcmpack/mitgcmdoc/part2/mom_vecinv.tex,v 1.2 2001/10/25 18:36:53 cnh Exp $
% $Name:  $

\section{Vector invariant momentum equations}

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