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

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

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


Revision 1.3 - (hide annotations) (download) (as text)
Fri Sep 28 03:28:32 2001 UTC (22 years, 7 months ago) by adcroft
Branch: MAIN
Changes since 1.2: +6 -6 lines
File MIME type: application/x-tex
Updated coriolis flag to match code.

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

  ViewVC Help
Powered by ViewVC 1.1.22