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	edhill	1.4	% $Header: /u/gcmpack/manual/part2/mom_vecinv.tex,v 1.3 2001/11/13 15:20:12 adcroft Exp $
2	cnh	1.2	% $Name: $
3	adcroft	1.1
4			\section{Vector invariant momentum equations}
5	edhill	1.4	\begin{rawhtml}
6			<!-- CMIREDIR:vector_invariant_momentum_eqautions: -->
7			\end{rawhtml}
8	adcroft	1.1
9			The finite volume method lends itself to describing the continuity and
10	adcroft	1.3	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	adcroft	1.1
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	cnh	1.2	numerical stability; alternative definitions break the conservation
64	adcroft	1.1	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	cnh	1.2	distinction between using absolute vorticity or relative vorticity is
139	adcroft	1.1	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} }