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	jmc	1.7	% $Header: /u/gcmpack/manual/s_algorithm/text/mom_vecinv.tex,v 1.6 2006/06/26 01:03:47 jmc Exp $
2	cnh	1.2	% $Name: $
3	adcroft	1.1
4			\section{Vector invariant momentum equations}
5	jmc	1.7	\label{sec:vect-inv_momentum_equations}
6	edhill	1.4	\begin{rawhtml}
7			<!-- CMIREDIR:vector_invariant_momentum_eqautions: -->
8			\end{rawhtml}
9	adcroft	1.1
10			The finite volume method lends itself to describing the continuity and
11	adcroft	1.3	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	adcroft	1.1
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	jmc	1.6	{\em S/R MOM\_VECINV} ({\em pkg/mom\_vecinv/mom\_vecinv.F})
52	adcroft	1.1
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	cnh	1.2	numerical stability; alternative definitions break the conservation
65	adcroft	1.1	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	jmc	1.6	$\zeta_3$: {\bf vort3} (local to {\em mom\_vecinv.F})
79	adcroft	1.1	\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	jmc	1.6	$KE$: {\bf KE} (local to {\em mom\_vecinv.F})
94	adcroft	1.1	\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	cnh	1.2	distinction between using absolute vorticity or relative vorticity is
140	adcroft	1.1	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	jmc	1.6	$G_u^{fv}$, $G_u^{\zeta_3 v}$: {\bf uCf} (local to {\em mom\_vecinv.F})
154	adcroft	1.1
155	jmc	1.6	$G_v^{fu}$, $G_v^{\zeta_3 u}$: {\bf vCf} (local to {\em mom\_vecinv.F})
156	adcroft	1.1	\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	jmc	1.6	$G_u^{\zeta_2 w}$: {\bf uCf} (local to {\em mom\_vecinv.F})
183	adcroft	1.1
184	jmc	1.6	$G_v^{\zeta_1 w}$: {\bf vCf} (local to {\em mom\_vecinv.F})
185	adcroft	1.1	\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	jmc	1.6	$G_u^{\partial_x KE}$: {\bf uCf} (local to {\em mom\_vecinv.F})
206	adcroft	1.1
207	jmc	1.6	$G_v^{\partial_y KE}$: {\bf vCf} (local to {\em mom\_vecinv.F})
208	adcroft	1.1	\end{minipage} }
209
210
211
212	jmc	1.6	\subsection{Horizontal divergence}
213	adcroft	1.1
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	jmc	1.6	$D$: {\bf hDiv} (local to {\em mom\_vecinv.F})
227	adcroft	1.1	\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	jmc	1.6	$G_u^{h-dissip}$: {\bf uDiss} (local to {\em mom\_vecinv.F})
258	adcroft	1.1
259	jmc	1.6	$G_v^{h-dissip}$: {\bf vDiss} (local to {\em mom\_vecinv.F})
260	adcroft	1.1	\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	jmc	1.6	$\tau_{13}$: {\bf urf} (local to {\em mom\_vecinv.F})
286	adcroft	1.1
287	jmc	1.6	$\tau_{23}$: {\bf vrf} (local to {\em mom\_vecinv.F})
288	adcroft	1.1	\end{minipage} }