/[MITgcm]/manual/s_algorithm/text/mom_fluxform.tex
ViewVC logotype

Contents of /manual/s_algorithm/text/mom_fluxform.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.5 - (show annotations) (download) (as text)
Thu Oct 25 12:06:56 2001 UTC (23 years, 8 months ago) by cnh
Branch: MAIN
Changes since 1.4: +2 -1 lines
File MIME type: application/x-tex
More editing

1 % $Header: /u/u0/gcmpack/mitgcmdoc/part2/mom_fluxform.tex,v 1.4 2001/10/25 00:55:28 cnh Exp $
2 % $Name: $
3
4 \section{Flux-form momentum equations}
5 \label{sec:flux-form_momentum_eqautions}
6
7 The original finite volume model was based on the Eulerian flux form
8 momentum equations. This is the default though the vector invariant
9 form is optionally available (and recommended in some cases).
10
11 The ``G's'' (our colloquial name for all terms on rhs!) are broken
12 into the various advective, Coriolis, horizontal dissipation, vertical
13 dissipation and metric forces:
14 \marginpar{$G_u$: {\bf Gu} }
15 \marginpar{$G_v$: {\bf Gv} }
16 \marginpar{$G_w$: {\bf Gw} }
17 \begin{eqnarray}
18 G_u & = & G_u^{adv} + G_u^{cor} + G_u^{h-diss} + G_u^{v-diss} +
19 G_u^{metric} + G_u^{nh-metric} \label{eq:gsplit_momu} \\
20 G_v & = & G_v^{adv} + G_v^{cor} + G_v^{h-diss} + G_v^{v-diss} +
21 G_v^{metric} + G_v^{nh-metric} \label{eq:gsplit_momv} \\
22 G_w & = & G_w^{adv} + G_w^{cor} + G_w^{h-diss} + G_w^{v-diss} +
23 G_w^{metric} + G_w^{nh-metric} \label{eq:gsplit_momw}
24 \end{eqnarray}
25 In the hydrostatic limit, $G_w=0$ and $\epsilon_{nh}=0$, reducing the
26 vertical momentum to hydrostatic balance.
27
28 These terms are calculated in routines called from subroutine {\em
29 CALC\_MOM\_RHS} a collected into the global arrays {\bf Gu}, {\bf Gv},
30 and {\bf Gw}.
31
32 \fbox{ \begin{minipage}{4.75in}
33 {\em S/R CALC\_MOM\_RHS} ({\em pkg/mom\_fluxform/calc\_mom\_rhs.F})
34
35 $G_u$: {\bf Gu} ({\em DYNVARS.h})
36
37 $G_v$: {\bf Gv} ({\em DYNVARS.h})
38
39 $G_w$: {\bf Gw} ({\em DYNVARS.h})
40 \end{minipage} }
41
42
43 \subsection{Advection of momentum}
44
45 The advective operator is second order accurate in space:
46 \begin{eqnarray}
47 {\cal A}_w \Delta r_f h_w G_u^{adv} & = &
48 \delta_i \overline{ U }^i \overline{ u }^i
49 + \delta_j \overline{ V }^i \overline{ u }^j
50 + \delta_k \overline{ W }^i \overline{ u }^k \label{eq:discrete-momadvu} \\
51 {\cal A}_s \Delta r_f h_s G_v^{adv} & = &
52 \delta_i \overline{ U }^j \overline{ v }^i
53 + \delta_j \overline{ V }^j \overline{ v }^j
54 + \delta_k \overline{ W }^j \overline{ v }^k \label{eq:discrete-momadvv} \\
55 {\cal A}_c \Delta r_c G_w^{adv} & = &
56 \delta_i \overline{ U }^k \overline{ w }^i
57 + \delta_j \overline{ V }^k \overline{ w }^j
58 + \delta_k \overline{ W }^k \overline{ w }^k \label{eq:discrete-momadvw}
59 \end{eqnarray}
60 and because of the flux form does not contribute to the global budget
61 of linear momentum. The quantities $U$, $V$ and $W$ are volume fluxes
62 defined:
63 \marginpar{$U$: {\bf uTrans} }
64 \marginpar{$V$: {\bf vTrans} }
65 \marginpar{$W$: {\bf rTrans} }
66 \begin{eqnarray}
67 U & = & \Delta y_g \Delta r_f h_w u \label{eq:utrans} \\
68 V & = & \Delta x_g \Delta r_f h_s v \label{eq:vtrans} \\
69 W & = & {\cal A}_c w \label{eq:rtrans}
70 \end{eqnarray}
71 The advection of momentum takes the same form as the advection of
72 tracers but by a translated advective flow. Consequently, the
73 conservation of second moments, derived for tracers later, applies to
74 $u^2$ and $v^2$ and $w^2$ so that advection of momentum correctly
75 conserves kinetic energy.
76
77 \fbox{ \begin{minipage}{4.75in}
78 {\em S/R MOM\_U\_ADV\_UU} ({\em mom\_u\_adv\_uu.F})
79
80 {\em S/R MOM\_U\_ADV\_VU} ({\em mom\_u\_adv\_vu.F})
81
82 {\em S/R MOM\_U\_ADV\_WU} ({\em mom\_u\_adv\_wu.F})
83
84 {\em S/R MOM\_U\_ADV\_UV} ({\em mom\_u\_adv\_uv.F})
85
86 {\em S/R MOM\_U\_ADV\_VV} ({\em mom\_u\_adv\_vv.F})
87
88 {\em S/R MOM\_U\_ADV\_WV} ({\em mom\_u\_adv\_wv.F})
89
90 $uu$, $uv$, $vu$, $vv$: {\bf aF} (local to {\em calc\_mom\_rhs.F})
91 \end{minipage} }
92
93
94
95 \subsection{Coriolis terms}
96
97 The ``pure C grid'' Coriolis terms (i.e. in absence of C-D scheme) are
98 discretized:
99 \begin{eqnarray}
100 {\cal A}_w \Delta r_f h_w G_u^{Cor} & = &
101 \overline{ f {\cal A}_c \Delta r_f h_c \overline{ v }^j }^i
102 - \epsilon_{nh} \overline{ f' {\cal A}_c \Delta r_f h_c \overline{ w }^k }^i \\
103 {\cal A}_s \Delta r_f h_s G_v^{Cor} & = &
104 - \overline{ f {\cal A}_c \Delta r_f h_c \overline{ u }^i }^j \\
105 {\cal A}_c \Delta r_c G_w^{Cor} & = &
106 \epsilon_{nh} \overline{ f' {\cal A}_c \Delta r_f h_c \overline{ u }^i }^k
107 \end{eqnarray}
108 where the Coriolis parameters $f$ and $f'$ are defined:
109 \begin{eqnarray}
110 f & = & 2 \Omega \sin{\varphi} \\
111 f' & = & 2 \Omega \cos{\varphi}
112 \end{eqnarray}
113 where $\varphi$ is geographic latitude when using spherical geometry,
114 otherwise the $\beta$-plane definition is used:
115 \begin{eqnarray}
116 f & = & f_o + \beta y \\
117 f' & = & 0
118 \end{eqnarray}
119
120 This discretization globally conserves kinetic energy. It should be
121 noted that despite the use of this discretization in former
122 publications, all calculations to date have used the following
123 different discretization:
124 \begin{eqnarray}
125 G_u^{Cor} & = &
126 f_u \overline{ v }^{ji}
127 - \epsilon_{nh} f_u' \overline{ w }^{ik} \\
128 G_v^{Cor} & = &
129 - f_v \overline{ u }^{ij} \\
130 G_w^{Cor} & = &
131 \epsilon_{nh} f_w' \overline{ u }^{ik}
132 \end{eqnarray}
133 \marginpar{Need to change the default in code to match this}
134 where the subscripts on $f$ and $f'$ indicate evaluation of the
135 Coriolis parameters at the appropriate points in space. The above
136 discretization does {\em not} conserve anything, especially energy and
137 for historical reasons is the default for the code. A
138 flag controls this discretization: set run-time logical {\bf
139 useEnergyConservingCoriolis} to {\em true} which otherwise defaults to
140 {\em false}.
141
142 \fbox{ \begin{minipage}{4.75in}
143 {\em S/R MOM\_CDSCHEME} ({\em mom\_cdscheme.F})
144
145 {\em S/R MOM\_U\_CORIOLIS} ({\em mom\_u\_coriolis.F})
146
147 {\em S/R MOM\_V\_CORIOLIS} ({\em mom\_v\_coriolis.F})
148
149 $G_u^{Cor}$, $G_v^{Cor}$: {\bf cF} (local to {\em calc\_mom\_rhs.F})
150 \end{minipage} }
151
152
153 \subsection{Curvature metric terms}
154
155 The most commonly used coordinate system on the sphere is the
156 geographic system $(\lambda,\varphi)$. The curvilinear nature of these
157 coordinates on the sphere lead to some ``metric'' terms in the
158 component momentum equations. Under the thin-atmosphere and
159 hydrostatic approximations these terms are discretized:
160 \begin{eqnarray}
161 {\cal A}_w \Delta r_f h_w G_u^{metric} & = &
162 \overline{ \frac{ \overline{u}^i }{a} \tan{\varphi} {\cal A}_c \Delta r_f h_c \overline{ v }^j }^i \\
163 {\cal A}_s \Delta r_f h_s G_v^{metric} & = &
164 - \overline{ \frac{ \overline{u}^i }{a} \tan{\varphi} {\cal A}_c \Delta r_f h_c \overline{ u }^i }^j \\
165 G_w^{metric} & = & 0
166 \end{eqnarray}
167 where $a$ is the radius of the planet (sphericity is assumed) or the
168 radial distance of the particle (i.e. a function of height). It is
169 easy to see that this discretization satisfies all the properties of
170 the discrete Coriolis terms since the metric factor $\frac{u}{a}
171 \tan{\varphi}$ can be viewed as a modification of the vertical Coriolis
172 parameter: $f \rightarrow f+\frac{u}{a} \tan{\varphi}$.
173
174 However, as for the Coriolis terms, a non-energy conserving form has
175 exclusively been used to date:
176 \begin{eqnarray}
177 G_u^{metric} & = & \frac{u \overline{v}^{ij} }{a} \tan{\varphi} \\
178 G_v^{metric} & = & \frac{ \overline{u}^{ij} \overline{u}^{ij}}{a} \tan{\varphi}
179 \end{eqnarray}
180 where $\tan{\varphi}$ is evaluated at the $u$ and $v$ points
181 respectively.
182
183 \fbox{ \begin{minipage}{4.75in}
184 {\em S/R MOM\_U\_METRIC\_SPHERE} ({\em mom\_u\_metric\_sphere.F})
185
186 {\em S/R MOM\_V\_METRIC\_SPHERE} ({\em mom\_v\_metric\_sphere.F})
187
188 $G_u^{metric}$, $G_v^{metric}$: {\bf mT} (local to {\em calc\_mom\_rhs.F})
189 \end{minipage} }
190
191
192
193 \subsection{Non-hydrostatic metric terms}
194
195 For the non-hydrostatic equations, dropping the thin-atmosphere
196 approximation re-introduces metric terms involving $w$ and are
197 required to conserve anglular momentum:
198 \begin{eqnarray}
199 {\cal A}_w \Delta r_f h_w G_u^{metric} & = &
200 - \overline{ \frac{ \overline{u}^i \overline{w}^k }{a} {\cal A}_c \Delta r_f h_c }^i \\
201 {\cal A}_s \Delta r_f h_s G_v^{metric} & = &
202 - \overline{ \frac{ \overline{v}^j \overline{w}^k }{a} {\cal A}_c \Delta r_f h_c}^j \\
203 {\cal A}_c \Delta r_c G_w^{metric} & = &
204 \overline{ \frac{ {\overline{u}^i}^2 + {\overline{v}^j}^2}{a} {\cal A}_c \Delta r_f h_c }^k
205 \end{eqnarray}
206
207 Because we are always consistent, even if consistently wrong, we have,
208 in the past, used a different discretization in the model which is:
209 \begin{eqnarray}
210 G_u^{metric} & = &
211 - \frac{u}{a} \overline{w}^{ik} \\
212 G_v^{metric} & = &
213 - \frac{v}{a} \overline{w}^{jk} \\
214 G_w^{metric} & = &
215 \frac{1}{a} ( {\overline{u}^{ik}}^2 + {\overline{v}^{jk}}^2 )
216 \end{eqnarray}
217
218 \fbox{ \begin{minipage}{4.75in}
219 {\em S/R MOM\_U\_METRIC\_NH} ({\em mom\_u\_metric\_nh.F})
220
221 {\em S/R MOM\_V\_METRIC\_NH} ({\em mom\_v\_metric\_nh.F})
222
223 $G_u^{metric}$, $G_v^{metric}$: {\bf mT} (local to {\em calc\_mom\_rhs.F})
224 \end{minipage} }
225
226
227 \subsection{Lateral dissipation}
228
229 Historically, we have represented the SGS Reynolds stresses as simply
230 down gradient momentum fluxes, ignoring constraints on the stress
231 tensor such as symmetry.
232 \begin{eqnarray}
233 {\cal A}_w \Delta r_f h_w G_u^{h-diss} & = &
234 \delta_i \Delta y_f \Delta r_f h_c \tau_{11}
235 + \delta_j \Delta x_v \Delta r_f h_\zeta \tau_{12} \\
236 {\cal A}_s \Delta r_f h_s G_v^{h-diss} & = &
237 \delta_i \Delta y_u \Delta r_f h_\zeta \tau_{21}
238 + \delta_j \Delta x_f \Delta r_f h_c \tau_{22}
239 \end{eqnarray}
240 \marginpar{Check signs of stress definitions}
241
242 The lateral viscous stresses are discretized:
243 \begin{eqnarray}
244 \tau_{11} & = & A_h c_{11\Delta}(\varphi) \frac{1}{\Delta x_f} \delta_i u
245 -A_4 c_{11\Delta^2}(\varphi) \frac{1}{\Delta x_f} \delta_i \nabla^2 u \\
246 \tau_{12} & = & A_h c_{12\Delta}(\varphi) \frac{1}{\Delta y_u} \delta_j u
247 -A_4 c_{12\Delta^2}(\varphi)\frac{1}{\Delta y_u} \delta_j \nabla^2 u \\
248 \tau_{21} & = & A_h c_{21\Delta}(\varphi) \frac{1}{\Delta x_v} \delta_i v
249 -A_4 c_{21\Delta^2}(\varphi) \frac{1}{\Delta x_v} \delta_i \nabla^2 v \\
250 \tau_{22} & = & A_h c_{22\Delta}(\varphi) \frac{1}{\Delta y_f} \delta_j v
251 -A_4 c_{22\Delta^2}(\varphi) \frac{1}{\Delta y_f} \delta_j \nabla^2 v
252 \end{eqnarray}
253 where the non-dimensional factors $c_{lm\Delta^n}(\varphi), \{l,m,n\} \in
254 \{1,2\}$ define the ``cosine'' scaling with latitude which can be
255 applied in various ad-hoc ways. For instance, $c_{11\Delta} =
256 c_{21\Delta} = (\cos{\varphi})^{3/2}$, $c_{12\Delta}=c_{22\Delta}=0$ would
257 represent the an-isotropic cosine scaling typically used on the
258 ``lat-lon'' grid for Laplacian viscosity.
259 \marginpar{Need to tidy up method for controlling this in code}
260
261 It should be noted that dispite the ad-hoc nature of the scaling, some
262 scaling must be done since on a lat-lon grid the converging meridians
263 make it very unlikely that a stable viscosity parameter exists across
264 the entire model domain.
265
266 The Laplacian viscosity coefficient, $A_h$ ({\bf viscAh}), has units
267 of $m^2 s^{-1}$. The bi-harmonic viscosity coefficient, $A_4$ ({\bf
268 viscA4}), has units of $m^4 s^{-1}$.
269
270 \fbox{ \begin{minipage}{4.75in}
271 {\em S/R MOM\_U\_XVISCFLUX} ({\em mom\_u\_xviscflux.F})
272
273 {\em S/R MOM\_U\_YVISCFLUX} ({\em mom\_u\_yviscflux.F})
274
275 {\em S/R MOM\_V\_XVISCFLUX} ({\em mom\_v\_xviscflux.F})
276
277 {\em S/R MOM\_V\_YVISCFLUX} ({\em mom\_v\_yviscflux.F})
278
279 $\tau_{11}$, $\tau_{12}$, $\tau_{22}$, $\tau_{22}$: {\bf vF}, {\bf
280 v4F} (local to {\em calc\_mom\_rhs.F})
281 \end{minipage} }
282
283 Two types of lateral boundary condition exist for the lateral viscous
284 terms, no-slip and free-slip.
285
286 The free-slip condition is most convenient to code since it is
287 equivalent to zero-stress on boundaries. Simple masking of the stress
288 components sets them to zero. The fractional open stress is properly
289 handled using the lopped cells.
290
291 The no-slip condition defines the normal gradient of a tangential flow
292 such that the flow is zero on the boundary. Rather than modify the
293 stresses by using complicated functions of the masks and ``ghost''
294 points (see \cite{Adcroft+Marshall98}) we add the boundary stresses as
295 an additional source term in cells next to solid boundaries. This has
296 the advantage of being able to cope with ``thin walls'' and also makes
297 the interior stress calculation (code) independent of the boundary
298 conditions. The ``body'' force takes the form:
299 \begin{eqnarray}
300 G_u^{side-drag} & = &
301 \frac{4}{\Delta z_f} \overline{ (1-h_\zeta) \frac{\Delta x_v}{\Delta y_u} }^j
302 \left( A_h c_{12\Delta}(\varphi) u - A_4 c_{12\Delta^2}(\varphi) \nabla^2 u \right)
303 \\
304 G_v^{side-drag} & = &
305 \frac{4}{\Delta z_f} \overline{ (1-h_\zeta) \frac{\Delta y_u}{\Delta x_v} }^i
306 \left( A_h c_{21\Delta}(\varphi) v - A_4 c_{21\Delta^2}(\varphi) \nabla^2 v \right)
307 \end{eqnarray}
308
309 In fact, the above discretization is not quite complete because it
310 assumes that the bathymetry at velocity points is deeper than at
311 neighbouring vorticity points, e.g. $1-h_w < 1-h_\zeta$
312
313 \fbox{ \begin{minipage}{4.75in}
314 {\em S/R MOM\_U\_SIDEDRAG} ({\em mom\_u\_sidedrag.F})
315
316 {\em S/R MOM\_V\_SIDEDRAG} ({\em mom\_v\_sidedrag.F})
317
318 $G_u^{side-drag}$, $G_v^{side-drag}$: {\bf vF} (local to {\em calc\_mom\_rhs.F})
319 \end{minipage} }
320
321
322 \subsection{Vertical dissipation}
323
324 Vertical viscosity terms are discretized with only partial adherence
325 to the variable grid lengths introduced by the finite volume
326 formulation. This reduces the formal accuracy of these terms to just
327 first order but only next to boundaries; exactly where other terms
328 appear such as linar and quadratic bottom drag.
329 \begin{eqnarray}
330 G_u^{v-diss} & = &
331 \frac{1}{\Delta r_f h_w} \delta_k \tau_{13} \\
332 G_v^{v-diss} & = &
333 \frac{1}{\Delta r_f h_s} \delta_k \tau_{23} \\
334 G_w^{v-diss} & = & \epsilon_{nh}
335 \frac{1}{\Delta r_f h_d} \delta_k \tau_{33}
336 \end{eqnarray}
337 represents the general discrete form of the vertical dissipation terms.
338
339 In the interior the vertical stresses are discretized:
340 \begin{eqnarray}
341 \tau_{13} & = & A_v \frac{1}{\Delta r_c} \delta_k u \\
342 \tau_{23} & = & A_v \frac{1}{\Delta r_c} \delta_k v \\
343 \tau_{33} & = & A_v \frac{1}{\Delta r_f} \delta_k w
344 \end{eqnarray}
345 It should be noted that in the non-hydrostatic form, the stress tensor
346 is even less consistent than for the hydrostatic (see Wazjowicz
347 \cite{Waojz}). It is well known how to do this properly (see Griffies
348 \cite{Griffies}) and is on the list of to-do's.
349
350 \fbox{ \begin{minipage}{4.75in}
351 {\em S/R MOM\_U\_RVISCLFUX} ({\em mom\_u\_rviscflux.F})
352
353 {\em S/R MOM\_V\_RVISCLFUX} ({\em mom\_v\_rviscflux.F})
354
355 $\tau_{13}$: {\bf urf} (local to {\em calc\_mom\_rhs.F})
356
357 $\tau_{23}$: {\bf vrf} (local to {\em calc\_mom\_rhs.F})
358 \end{minipage} }
359
360
361 As for the lateral viscous terms, the free-slip condition is
362 equivalent to simply setting the stress to zero on boundaries. The
363 no-slip condition is implemented as an additional term acting on top
364 of the interior and free-slip stresses. Bottom drag represents
365 additional friction, in addition to that imposed by the no-slip
366 condition at the bottom. The drag is cast as a stress expressed as a
367 linear or quadratic function of the mean flow in the layer above the
368 topography:
369 \begin{eqnarray}
370 \tau_{13}^{bottom-drag} & = &
371 \left(
372 2 A_v \frac{1}{\Delta r_c}
373 + r_b
374 + C_d \sqrt{ \overline{2 KE}^i }
375 \right) u \\
376 \tau_{23}^{bottom-drag} & = &
377 \left(
378 2 A_v \frac{1}{\Delta r_c}
379 + r_b
380 + C_d \sqrt{ \overline{2 KE}^j }
381 \right) v
382 \end{eqnarray}
383 where these terms are only evaluated immediately above topography.
384 $r_b$ ({\bf bottomDragLinear}) has units of $m s^{-1}$ and a typical value
385 of the order 0.0002 $m s^{-1}$. $C_d$ ({\bf bottomDragQuadratic}) is
386 dimensionless with typical values in the range 0.001--0.003.
387
388 \fbox{ \begin{minipage}{4.75in}
389 {\em S/R MOM\_U\_BOTTOMDRAG} ({\em mom\_u\_bottomdrag.F})
390
391 {\em S/R MOM\_V\_BOTTOMDRAG} ({\em mom\_v\_bottomdrag.F})
392
393 $\tau_{13}^{bottom-drag}$, $\tau_{23}^{bottom-drag}$: {\bf vf} (local to {\em calc\_mom\_rhs.F})
394 \end{minipage} }
395
396 \subsection{Derivation of discrete energy conservation}
397
398 These discrete equations conserve kinetic plus potential energy using the
399 following definitions:
400 \begin{equation}
401 KE = \frac{1}{2} \left( \overline{ u^2 }^i + \overline{ v^2 }^j +
402 \epsilon_{nh} \overline{ w^2 }^k \right)
403 \end{equation}
404

  ViewVC Help
Powered by ViewVC 1.1.22