1 |
% $Header: /u/gcmpack/mitgcmdoc/part2/mom_vecinv.tex,v 1.2 2001/10/25 18:36:53 cnh Exp $ |
2 |
% $Name: $ |
3 |
|
4 |
\section{Vector invariant momentum equations} |
5 |
|
6 |
The finite volume method lends itself to describing the continuity and |
7 |
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 |
|
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 |
numerical stability; alternative definitions break the conservation |
61 |
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 |
distinction between using absolute vorticity or relative vorticity is |
136 |
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} } |