/[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.1 - (show annotations) (download) (as text)
Thu Aug 9 19:48:39 2001 UTC (22 years, 9 months ago) by adcroft
Branch: MAIN
File MIME type: application/x-tex
Moved things around for more sections fewer subsections
Split spatial-discrete.tex into several smaller files.

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

  ViewVC Help
Powered by ViewVC 1.1.22