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

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

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


Revision 1.15 - (show annotations) (download) (as text)
Thu Feb 28 19:32:19 2002 UTC (23 years, 4 months ago) by cnh
Branch: MAIN
Changes since 1.14: +18 -6 lines
File MIME type: application/x-tex
Updates for special on-line version with
hyperlinked and animated figures

Separating tutorials and reference material

1 % $Header: /u/u0/gcmpack/manual/part2/time_stepping.tex,v 1.14 2001/11/14 21:07:13 adcroft Exp $
2 % $Name: $
3
4 This chapter lays out the numerical schemes that are
5 employed in the core MITgcm algorithm. Whenever possible
6 links are made to actual program code in the MITgcm implementation.
7 The chapter begins with a discussion of the temporal discretization
8 used in MITgcm. This discussion is followed by sections that
9 describe the spatial discretization. The schemes employed for momentum
10 terms are described first, afterwards the schemes that apply to
11 passive and dynamically active tracers are described.
12
13
14 \section{Time-stepping}
15 The equations of motion integrated by the model involve four
16 prognostic equations for flow, $u$ and $v$, temperature, $\theta$, and
17 salt/moisture, $S$, and three diagnostic equations for vertical flow,
18 $w$, density/buoyancy, $\rho$/$b$, and pressure/geo-potential,
19 $\phi_{hyd}$. In addition, the surface pressure or height may by
20 described by either a prognostic or diagnostic equation and if
21 non-hydrostatics terms are included then a diagnostic equation for
22 non-hydrostatic pressure is also solved. The combination of prognostic
23 and diagnostic equations requires a model algorithm that can march
24 forward prognostic variables while satisfying constraints imposed by
25 diagnostic equations.
26
27 Since the model comes in several flavors and formulation, it would be
28 confusing to present the model algorithm exactly as written into code
29 along with all the switches and optional terms. Instead, we present
30 the algorithm for each of the basic formulations which are:
31 \begin{enumerate}
32 \item the semi-implicit pressure method for hydrostatic equations
33 with a rigid-lid, variables co-located in time and with
34 Adams-Bashforth time-stepping, \label{it:a}
35 \item as \ref{it:a}. but with an implicit linear free-surface, \label{it:b}
36 \item as \ref{it:a}. or \ref{it:b}. but with variables staggered in time,
37 \label{it:c}
38 \item as \ref{it:a}. or \ref{it:b}. but with non-hydrostatic terms included,
39 \item as \ref{it:b}. or \ref{it:c}. but with non-linear free-surface.
40 \end{enumerate}
41
42 In all the above configurations it is also possible to substitute the
43 Adams-Bashforth with an alternative time-stepping scheme for terms
44 evaluated explicitly in time. Since the over-arching algorithm is
45 independent of the particular time-stepping scheme chosen we will
46 describe first the over-arching algorithm, known as the pressure
47 method, with a rigid-lid model in section
48 \ref{sect:pressure-method-rigid-lid}. This algorithm is essentially
49 unchanged, apart for some coefficients, when the rigid lid assumption
50 is replaced with a linearized implicit free-surface, described in
51 section \ref{sect:pressure-method-linear-backward}. These two flavors
52 of the pressure-method encompass all formulations of the model as it
53 exists today. The integration of explicit in time terms is out-lined
54 in section \ref{sect:adams-bashforth} and put into the context of the
55 overall algorithm in sections \ref{sect:adams-bashforth-sync} and
56 \ref{sect:adams-bashforth-staggered}. Inclusion of non-hydrostatic
57 terms requires applying the pressure method in three dimensions
58 instead of two and this algorithm modification is described in section
59 \ref{sect:non-hydrostatic}. Finally, the free-surface equation may be
60 treated more exactly, including non-linear terms, and this is
61 described in section \ref{sect:nonlinear-freesurface}.
62
63
64 \section{Pressure method with rigid-lid} \label{sect:pressure-method-rigid-lid}
65
66 \begin{figure}
67 \begin{center}
68 \resizebox{4.0in}{!}{\includegraphics{part2/pressure-method-rigid-lid.eps}}
69 \end{center}
70 \caption{
71 A schematic of the evolution in time of the pressure method
72 algorithm. A prediction for the flow variables at time level $n+1$ is
73 made based only on the explicit terms, $G^{(n+^1/_2)}$, and denoted
74 $u^*$, $v^*$. Next, a pressure field is found such that $u^{n+1}$,
75 $v^{n+1}$ will be non-divergent. Conceptually, the $*$ quantities
76 exist at time level $n+1$ but they are intermediate and only
77 temporary.}
78 \label{fig:pressure-method-rigid-lid}
79 \end{figure}
80
81 \begin{figure}
82 \begin{center} \fbox{ \begin{minipage}{4.5in} \begin{tabbing}
83 aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill
84 \proclink{FORWARD\_STEP}{../code/._model_src_forward_step.F} \\
85 \> DYNAMICS \\
86 \>\> TIMESTEP \` $u^*$,$v^*$ (\ref{eq:ustar-rigid-lid},\ref{eq:vstar-rigid-lid}) \\
87 \> SOLVE\_FOR\_PRESSURE \\
88 \>\> CALC\_DIV\_GHAT \` $H\widehat{u^*}$,$H\widehat{v^*}$ (\ref{eq:elliptic}) \\
89 \>\> CG2D \` $\eta^{n+1}$ (\ref{eq:elliptic}) \\
90 \> THE\_CORRECTION\_STEP \\
91 \>\> CALC\_GRAD\_PHI\_SURF \` $\nabla \eta^{n+1}$ \\
92 \>\> CORRECTION\_STEP \` $u^{n+1}$,$v^{n+1}$ (\ref{eq:un+1-rigid-lid},\ref{eq:vn+1-rigid-lid})
93 \end{tabbing} \end{minipage} } \end{center}
94 \caption{Calling tree for the pressure method algorihtm}
95 \label{fig:call-tree-pressure-method}
96 \end{figure}
97
98 The horizontal momentum and continuity equations for the ocean
99 (\ref{eq:ocean-mom} and \ref{eq:ocean-cont}), or for the atmosphere
100 (\ref{eq:atmos-mom} and \ref{eq:atmos-cont}), can be summarized by:
101 \begin{eqnarray}
102 \partial_t u + g \partial_x \eta & = & G_u \\
103 \partial_t v + g \partial_y \eta & = & G_v \\
104 \partial_x u + \partial_y v + \partial_z w & = & 0
105 \end{eqnarray}
106 where we are adopting the oceanic notation for brevity. All terms in
107 the momentum equations, except for surface pressure gradient, are
108 encapsulated in the $G$ vector. The continuity equation, when
109 integrated over the fluid depth, $H$, and with the rigid-lid/no normal
110 flow boundary conditions applied, becomes:
111 \begin{equation}
112 \partial_x H \widehat{u} + \partial_y H \widehat{v} = 0
113 \label{eq:rigid-lid-continuity}
114 \end{equation}
115 Here, $H\widehat{u} = \int_H u dz$ is the depth integral of $u$,
116 similarly for $H\widehat{v}$. The rigid-lid approximation sets $w=0$
117 at the lid so that it does not move but allows a pressure to be
118 exerted on the fluid by the lid. The horizontal momentum equations and
119 vertically integrated continuity equation are be discretized in time
120 and space as follows:
121 \begin{eqnarray}
122 u^{n+1} + \Delta t g \partial_x \eta^{n+1}
123 & = & u^{n} + \Delta t G_u^{(n+1/2)}
124 \label{eq:discrete-time-u}
125 \\
126 v^{n+1} + \Delta t g \partial_y \eta^{n+1}
127 & = & v^{n} + \Delta t G_v^{(n+1/2)}
128 \label{eq:discrete-time-v}
129 \\
130 \partial_x H \widehat{u^{n+1}}
131 + \partial_y H \widehat{v^{n+1}} & = & 0
132 \label{eq:discrete-time-cont-rigid-lid}
133 \end{eqnarray}
134 As written here, terms on the LHS all involve time level $n+1$ and are
135 referred to as implicit; the implicit backward time stepping scheme is
136 being used. All other terms in the RHS are explicit in time. The
137 thermodynamic quantities are integrated forward in time in parallel
138 with the flow and will be discussed later. For the purposes of
139 describing the pressure method it suffices to say that the hydrostatic
140 pressure gradient is explicit and so can be included in the vector
141 $G$.
142
143 Substituting the two momentum equations into the depth integrated
144 continuity equation eliminates $u^{n+1}$ and $v^{n+1}$ yielding an
145 elliptic equation for $\eta^{n+1}$. Equations
146 \ref{eq:discrete-time-u}, \ref{eq:discrete-time-v} and
147 \ref{eq:discrete-time-cont-rigid-lid} can then be re-arranged as follows:
148 \begin{eqnarray}
149 u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2)} \label{eq:ustar-rigid-lid} \\
150 v^{*} & = & v^{n} + \Delta t G_v^{(n+1/2)} \label{eq:vstar-rigid-lid} \\
151 \partial_x \Delta t g H \partial_x \eta^{n+1}
152 + \partial_y \Delta t g H \partial_y \eta^{n+1}
153 & = &
154 \partial_x H \widehat{u^{*}}
155 + \partial_y H \widehat{v^{*}} \label{eq:elliptic}
156 \\
157 u^{n+1} & = & u^{*} - \Delta t g \partial_x \eta^{n+1} \label{eq:un+1-rigid-lid}\\
158 v^{n+1} & = & v^{*} - \Delta t g \partial_y \eta^{n+1} \label{eq:vn+1-rigid-lid}
159 \end{eqnarray}
160 Equations \ref{eq:ustar-rigid-lid} to \ref{eq:vn+1-rigid-lid}, solved
161 sequentially, represent the pressure method algorithm used in the
162 model. The essence of the pressure method lies in the fact that any
163 explicit prediction for the flow would lead to a divergence flow field
164 so a pressure field must be found that keeps the flow non-divergent
165 over each step of the integration. The particular location in time of
166 the pressure field is somewhat ambiguous; in
167 Fig.~\ref{fig:pressure-method-rigid-lid} we depicted as co-located
168 with the future flow field (time level $n+1$) but it could equally
169 have been drawn as staggered in time with the flow.
170
171 The correspondence to the code is as follows:
172 \begin{itemize}
173 \item
174 the prognostic phase, equations \ref{eq:ustar-rigid-lid} and \ref{eq:vstar-rigid-lid},
175 stepping forward $u^n$ and $v^n$ to $u^{*}$ and $v^{*}$ is coded in
176 \proclink{TIMESTEP}{../code/._model_src_timestep.F}
177 \item
178 the vertical integration, $H \widehat{u^*}$ and $H
179 \widehat{v^*}$, divergence and inversion of the elliptic operator in
180 equation \ref{eq:elliptic} is coded in
181 \proclink{SOLVE\_FOR\_PRESSURE}{../code/._model_src_solve_for_pressure.F}
182 \item
183 finally, the new flow field at time level $n+1$ given by equations
184 \ref{eq:un+1-rigid-lid} and \ref{eq:vn+1-rigid-lid} is calculated in
185 \proclink{CORRECTION\_STEP}{../code/._model_src_correction_step.F}.
186 \end{itemize}
187 The calling tree for these routines is given in
188 Fig.~\ref{fig:call-tree-pressure-method}.
189
190
191
192 \paragraph{Need to discuss implicit viscosity somewhere:}
193 \begin{eqnarray}
194 \frac{1}{\Delta t} u^{n+1} - \partial_z A_v \partial_z u^{n+1}
195 + g \partial_x \eta^{n+1} & = & \frac{1}{\Delta t} u^{n} +
196 G_u^{(n+1/2)}
197 \\
198 \frac{1}{\Delta t} v^{n+1} - \partial_z A_v \partial_z v^{n+1}
199 + g \partial_y \eta^{n+1} & = & \frac{1}{\Delta t} v^{n} + G_v^{(n+1/2)}
200 \end{eqnarray}
201
202
203 \section{Pressure method with implicit linear free-surface}
204 \label{sect:pressure-method-linear-backward}
205
206 The rigid-lid approximation filters out external gravity waves
207 subsequently modifying the dispersion relation of barotropic Rossby
208 waves. The discrete form of the elliptic equation has some zero
209 eigen-values which makes it a potentially tricky or inefficient
210 problem to solve.
211
212 The rigid-lid approximation can be easily replaced by a linearization
213 of the free-surface equation which can be written:
214 \begin{equation}
215 \partial_t \eta + \partial_x H \widehat{u} + \partial_y H \widehat{v} = P-E+R
216 \label{eq:linear-free-surface=P-E+R}
217 \end{equation}
218 which differs from the depth integrated continuity equation with
219 rigid-lid (\ref{eq:rigid-lid-continuity}) by the time-dependent term
220 and fresh-water source term.
221
222 Equation \ref{eq:discrete-time-cont-rigid-lid} in the rigid-lid
223 pressure method is then replaced by the time discretization of
224 \ref{eq:linear-free-surface=P-E+R} which is:
225 \begin{equation}
226 \eta^{n+1}
227 + \Delta t \partial_x H \widehat{u^{n+1}}
228 + \Delta t \partial_y H \widehat{v^{n+1}}
229 =
230 \eta^{n}
231 + \Delta t ( P - E + R )
232 \label{eq:discrete-time-backward-free-surface}
233 \end{equation}
234 where the use of flow at time level $n+1$ makes the method implicit
235 and backward in time. The is the preferred scheme since it still
236 filters the fast unresolved wave motions by damping them. A centered
237 scheme, such as Crank-Nicholson, would alias the energy of the fast
238 modes onto slower modes of motion.
239
240 As for the rigid-lid pressure method, equations
241 \ref{eq:discrete-time-u}, \ref{eq:discrete-time-v} and
242 \ref{eq:discrete-time-backward-free-surface} can be re-arranged as follows:
243 \begin{eqnarray}
244 u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2)} \label{eq:ustar-backward-free-surface} \\
245 v^{*} & = & v^{n} + \Delta t G_v^{(n+1/2)} \label{eq:vstar-backward-free-surface} \\
246 \eta^* & = & \epsilon_{fs} \left( \eta^{n} +P-E+R \right)- \Delta t
247 \partial_x H \widehat{u^{*}}
248 + \partial_y H \widehat{v^{*}}
249 \\
250 \partial_x g H \partial_x \eta^{n+1}
251 + \partial_y g H \partial_y \eta^{n+1}
252 - \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2}
253 & = &
254 - \frac{\eta^*}{\Delta t^2}
255 \label{eq:elliptic-backward-free-surface}
256 \\
257 u^{n+1} & = & u^{*} - \Delta t g \partial_x \eta^{n+1} \label{eq:un+1-backward-free-surface}\\
258 v^{n+1} & = & v^{*} - \Delta t g \partial_y \eta^{n+1} \label{eq:vn+1-backward-free-surface}
259 \end{eqnarray}
260 Equations~\ref{eq:ustar-backward-free-surface}
261 to~\ref{eq:vn+1-backward-free-surface}, solved sequentially, represent
262 the pressure method algorithm with a backward implicit, linearized
263 free surface. The method is still formerly a pressure method because
264 in the limit of large $\Delta t$ the rigid-lid method is
265 recovered. However, the implicit treatment of the free-surface allows
266 the flow to be divergent and for the surface pressure/elevation to
267 respond on a finite time-scale (as opposed to instantly). To recover
268 the rigid-lid formulation, we introduced a switch-like parameter,
269 $\epsilon_{fs}$, which selects between the free-surface and rigid-lid;
270 $\epsilon_{fs}=1$ allows the free-surface to evolve; $\epsilon_{fs}=0$
271 imposes the rigid-lid. The evolution in time and location of variables
272 is exactly as it was for the rigid-lid model so that
273 Fig.~\ref{fig:pressure-method-rigid-lid} is still
274 applicable. Similarly, the calling sequence, given in
275 Fig.~\ref{fig:call-tree-pressure-method}, is as for the
276 pressure-method.
277
278
279 \section{Explicit time-stepping: Adams-Bashforth}
280 \label{sect:adams-bashforth}
281
282 In describing the the pressure method above we deferred describing the
283 time discretization of the explicit terms. We have historically used
284 the quasi-second order Adams-Bashforth method for all explicit terms
285 in both the momentum and tracer equations. This is still the default
286 mode of operation but it is now possible to use alternate schemes for
287 tracers (see section \ref{sect:tracer-advection}).
288
289 \begin{figure}
290 \begin{center} \fbox{ \begin{minipage}{4.5in} \begin{tabbing}
291 aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill
292 FORWARD\_STEP \\
293 \> THERMODYNAMICS \\
294 \>\> CALC\_GT \\
295 \>\>\> GAD\_CALC\_RHS \` $G_\theta^n = G_\theta( u, \theta^n )$ \\
296 \>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\
297 \>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:adams-bashforth2}) \\
298 \>\> TIMESTEP\_TRACER \` $\tau^*$ (\ref{eq:taustar}) \\
299 \>\> IMPLDIFF \` $\tau^{(n+1)}$ (\ref{eq:tau-n+1-implicit})
300 \end{tabbing} \end{minipage} } \end{center}
301 \caption{
302 Calling tree for the Adams-Bashforth time-stepping of temperature with
303 implicit diffusion.}
304 \label{fig:call-tree-adams-bashforth}
305 \end{figure}
306
307 In the previous sections, we summarized an explicit scheme as:
308 \begin{equation}
309 \tau^{*} = \tau^{n} + \Delta t G_\tau^{(n+1/2)}
310 \label{eq:taustar}
311 \end{equation}
312 where $\tau$ could be any prognostic variable ($u$, $v$, $\theta$ or
313 $S$) and $\tau^*$ is an explicit estimate of $\tau^{n+1}$ and would be
314 exact if not for implicit-in-time terms. The parenthesis about $n+1/2$
315 indicates that the term is explicit and extrapolated forward in time
316 and for this we use the quasi-second order Adams-Bashforth method:
317 \begin{equation}
318 G_\tau^{(n+1/2)} = ( 3/2 + \epsilon_{AB}) G_\tau^n
319 - ( 1/2 + \epsilon_{AB}) G_\tau^{n-1}
320 \label{eq:adams-bashforth2}
321 \end{equation}
322 This is a linear extrapolation, forward in time, to
323 $t=(n+1/2+{\epsilon_{AB}})\Delta t$. An extrapolation to the mid-point
324 in time, $t=(n+1/2)\Delta t$, corresponding to $\epsilon_{AB}=0$,
325 would be second order accurate but is weakly unstable for oscillatory
326 terms. A small but finite value for $\epsilon_{AB}$ stabilizes the
327 method. Strictly speaking, damping terms such as diffusion and
328 dissipation, and fixed terms (forcing), do not need to be inside the
329 Adams-Bashforth extrapolation. However, in the current code, it is
330 simpler to include these terms and this can be justified if the flow
331 and forcing evolves smoothly. Problems can, and do, arise when forcing
332 or motions are high frequency and this corresponds to a reduced
333 stability compared to a simple forward time-stepping of such terms.
334
335 A stability analysis for an oscillation equation should be given at this point.
336 \marginpar{AJA needs to find his notes on this...}
337
338 A stability analysis for a relaxation equation should be given at this point.
339 \marginpar{...and for this too.}
340
341
342 \section{Implicit time-stepping: backward method}
343
344 Vertical diffusion and viscosity can be treated implicitly in time
345 using the backward method which is an intrinsic scheme. For tracers,
346 the time discretized equation is:
347 \begin{equation}
348 \tau^{n+1} - \Delta t \partial_r \kappa_v \partial_r \tau^{n+1} =
349 \tau^{n} + \Delta t G_\tau^{(n+1/2)}
350 \label{eq:implicit-diffusion}
351 \end{equation}
352 where $G_\tau^{(n+1/2)}$ is the remaining explicit terms extrapolated
353 using the Adams-Bashforth method as described above. Equation
354 \ref{eq:implicit-diffusion} can be split split into:
355 \begin{eqnarray}
356 \tau^* & = & \tau^{n} + \Delta t G_\tau^{(n+1/2)}
357 \label{eq:taustar-implicit} \\
358 \tau^{n+1} & = & {\cal L}_\tau^{-1} ( \tau^* )
359 \label{eq:tau-n+1-implicit}
360 \end{eqnarray}
361 where ${\cal L}_\tau^{-1}$ is the inverse of the operator
362 \begin{equation}
363 {\cal L} = \left[ 1 + \Delta t \partial_r \kappa_v \partial_r \right]
364 \end{equation}
365 Equation \ref{eq:taustar-implicit} looks exactly as \ref{eq:taustar}
366 while \ref{eq:tau-n+1-implicit} involves an operator or matrix
367 inversion. By re-arranging \ref{eq:implicit-diffusion} in this way we
368 have cast the method as an explicit prediction step and an implicit
369 step allowing the latter to be inserted into the over all algorithm
370 with minimal interference.
371
372 Fig.~\ref{fig:call-tree-adams-bashforth} shows the calling sequence for
373 stepping forward a tracer variable such as temperature.
374
375 In order to fit within the pressure method, the implicit viscosity
376 must not alter the barotropic flow. In other words, it can on ly
377 redistribute momentum in the vertical. The upshot of this is that
378 although vertical viscosity may be backward implicit and
379 unconditionally stable, no-slip boundary conditions may not be made
380 implicit and are thus cast as a an explicit drag term.
381
382 \section{Synchronous time-stepping: variables co-located in time}
383 \label{sect:adams-bashforth-sync}
384
385 \begin{figure}
386 \begin{center}
387 \resizebox{5.0in}{!}{\includegraphics{part2/adams-bashforth-sync.eps}}
388 \end{center}
389 \caption{
390 A schematic of the explicit Adams-Bashforth and implicit time-stepping
391 phases of the algorithm. All prognostic variables are co-located in
392 time. Explicit tendencies are evaluated at time level $n$ as a
393 function of the state at that time level (dotted arrow). The explicit
394 tendency from the previous time level, $n-1$, is used to extrapolate
395 tendencies to $n+1/2$ (dashed arrow). This extrapolated tendency
396 allows variables to be stably integrated forward-in-time to render an
397 estimate ($*$-variables) at the $n+1$ time level (solid
398 arc-arrow). The operator ${\cal L}$ formed from implicit-in-time terms
399 is solved to yield the state variables at time level $n+1$. }
400 \label{fig:adams-bashforth-sync}
401 \end{figure}
402
403 \begin{figure}
404 \begin{center} \fbox{ \begin{minipage}{4.5in} \begin{tabbing}
405 aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill
406 FORWARD\_STEP \\
407 \> THERMODYNAMICS \\
408 \>\> CALC\_GT \\
409 \>\>\> GAD\_CALC\_RHS \` $G_\theta^n = G_\theta( u, \theta^n )$ (\ref{eq:Gt-n-sync})\\
410 \>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\
411 \>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:Gt-n+5-sync}) \\
412 \>\> TIMESTEP\_TRACER \` $\tau^*$ (\ref{eq:tstar-sync}) \\
413 \>\> IMPLDIFF \` $\tau^{(n+1)}$ (\ref{eq:t-n+1-sync}) \\
414 \> DYNAMICS \\
415 \>\> CALC\_PHI\_HYD \` $\phi_{hyd}^n$ (\ref{eq:phi-hyd-sync}) \\
416 \>\> MOM\_FLUXFORM or MOM\_VECINV \` $G_{\vec{\bf v}}^n$ (\ref{eq:Gv-n-sync})\\
417 \>\> TIMESTEP \` $\vec{\bf v}^*$ (\ref{eq:Gv-n+5-sync}, \ref{eq:vstar-sync}) \\
418 \>\> IMPLDIFF \` $\vec{\bf v}^{**}$ (\ref{eq:vstarstar-sync}) \\
419 \> SOLVE\_FOR\_PRESSURE \\
420 \>\> CALC\_DIV\_GHAT \` $\eta^*$ (\ref{eq:nstar-sync}) \\
421 \>\> CG2D \` $\eta^{n+1}$ (\ref{eq:elliptic-sync}) \\
422 \> THE\_CORRECTION\_STEP \\
423 \>\> CALC\_GRAD\_PHI\_SURF \` $\nabla \eta^{n+1}$ \\
424 \>\> CORRECTION\_STEP \` $u^{n+1}$,$v^{n+1}$ (\ref{eq:v-n+1-sync})
425 \end{tabbing} \end{minipage} } \end{center}
426 \caption{
427 Calling tree for the overall synchronous algorithm using
428 Adams-Bashforth time-stepping.}
429 \label{fig:call-tree-adams-bashforth-sync}
430 \end{figure}
431
432 The Adams-Bashforth extrapolation of explicit tendencies fits neatly
433 into the pressure method algorithm when all state variables are
434 co-located in time. Fig.~\ref{fig:adams-bashforth-sync} illustrates
435 the location of variables in time and the evolution of the algorithm
436 with time. The algorithm can be represented by the sequential solution
437 of the follow equations:
438 \begin{eqnarray}
439 G_{\theta,S}^{n} & = & G_{\theta,S} ( u^n, \theta^n, S^n )
440 \label{eq:Gt-n-sync} \\
441 G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsilon_{AB}) G_{\theta,S}^{n}-(1/2+\epsilon_{AB}) G_{\theta,S}^{n-1}
442 \label{eq:Gt-n+5-sync} \\
443 (\theta^*,S^*) & = & (\theta^{n},S^{n}) + \Delta t G_{\theta,S}^{(n+1/2)}
444 \label{eq:tstar-sync} \\
445 (\theta^{n+1},S^{n+1}) & = & {\cal L}^{-1}_{\theta,S} (\theta^*,S^*)
446 \label{eq:t-n+1-sync} \\
447 \phi^n_{hyd} & = & \int b(\theta^n,S^n) dr
448 \label{eq:phi-hyd-sync} \\
449 \vec{\bf G}_{\vec{\bf v}}^{n} & = & \vec{\bf G}_{\vec{\bf v}} ( \vec{\bf v}^n, \phi^n_{hyd} )
450 \label{eq:Gv-n-sync} \\
451 \vec{\bf G}_{\vec{\bf v}}^{(n+1/2)} & = & (3/2 + \epsilon_{AB} ) \vec{\bf G}_{\vec{\bf v}}^{n} - (1/2 + \epsilon_{AB} ) \vec{\bf G}_{\vec{\bf v}}^{n-1}
452 \label{eq:Gv-n+5-sync} \\
453 \vec{\bf v}^{*} & = & \vec{\bf v}^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}}^{(n+1/2)}
454 \label{eq:vstar-sync} \\
455 \vec{\bf v}^{**} & = & {\cal L}_{\vec{\bf v}}^{-1} ( \vec{\bf v}^* )
456 \label{eq:vstarstar-sync} \\
457 \eta^* & = & \epsilon_{fs} \left( \eta^{n} +P-E+R \right)- \Delta t
458 \nabla \cdot H \widehat{ \vec{\bf v}^{**} }
459 \label{eq:nstar-sync} \\
460 \nabla \cdot g H \nabla \eta^{n+1} - \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2}
461 & = & - \frac{\eta^*}{\Delta t^2}
462 \label{eq:elliptic-sync} \\
463 \vec{\bf v}^{n+1} & = & \vec{\bf v}^{*} - \Delta t g \nabla \eta^{n+1}
464 \label{eq:v-n+1-sync}
465 \end{eqnarray}
466 Fig.~\ref{fig:adams-bashforth-sync} illustrates the location of
467 variables in time and evolution of the algorithm with time. The
468 Adams-Bashforth extrapolation of the tracer tendencies is illustrated
469 by the dashed arrow, the prediction at $n+1$ is indicated by the
470 solid arc. Inversion of the implicit terms, ${\cal
471 L}^{-1}_{\theta,S}$, then yields the new tracer fields at $n+1$. All
472 these operations are carried out in subroutine {\em THERMODYNAMICS} an
473 subsidiaries, which correspond to equations \ref{eq:Gt-n-sync} to
474 \ref{eq:t-n+1-sync}.
475 Similarly illustrated is the Adams-Bashforth extrapolation of
476 accelerations, stepping forward and solving of implicit viscosity and
477 surface pressure gradient terms, corresponding to equations
478 \ref{eq:Gv-n-sync} to \ref{eq:v-n+1-sync}.
479 These operations are carried out in subroutines {\em DYNAMCIS}, {\em
480 SOLVE\_FOR\_PRESSURE} and {\em THE\_CORRECTION\_STEP}. This, then,
481 represents an entire algorithm for stepping forward the model one
482 time-step. The corresponding calling tree is given in
483 \ref{fig:call-tree-adams-bashforth-sync}.
484
485 \section{Staggered baroclinic time-stepping}
486 \label{sect:adams-bashforth-staggered}
487
488 \begin{figure}
489 \begin{center}
490 \resizebox{5.5in}{!}{\includegraphics{part2/adams-bashforth-staggered.eps}}
491 \end{center}
492 \caption{
493 A schematic of the explicit Adams-Bashforth and implicit time-stepping
494 phases of the algorithm but with staggering in time of thermodynamic
495 variables with the flow. Explicit thermodynamics tendencies are
496 evaluated at time level $n-1/2$ as a function of the thermodynamics
497 state at that time level $n$ and flow at time $n$ (dotted arrow). The
498 explicit tendency from the previous time level, $n-3/2$, is used to
499 extrapolate tendencies to $n$ (dashed arrow). This extrapolated
500 tendency allows thermo-dynamics variables to be stably integrated
501 forward-in-time to render an estimate ($*$-variables) at the $n+1/2$
502 time level (solid arc-arrow). The implicit-in-time operator ${\cal
503 L_{\theta,S}}$ is solved to yield the thermodynamic variables at time
504 level $n+1/2$. These are then used to calculate the hydrostatic
505 pressure/geo-potential, $\phi_{hyd}$ (vertical arrows). The
506 hydrostatic pressure gradient is evaluated directly an time level
507 $n+1/2$ in stepping forward the flow variables from $n$ to $n+1$
508 (solid arc-arrow). }
509 \label{fig:adams-bashforth-staggered}
510 \end{figure}
511
512 For well stratified problems, internal gravity waves may be the
513 limiting process for determining a stable time-step. In the
514 circumstance, it is more efficient to stagger in time the
515 thermodynamic variables with the flow
516 variables. Fig.~\ref{fig:adams-bashforth-staggered} illustrates the
517 staggering and algorithm. The key difference between this and
518 Fig.~\ref{fig:adams-bashforth-sync} is that the new thermodynamics
519 fields are used to compute the hydrostatic pressure at time level
520 $n+1/2$. The essentially allows the gravity wave terms to leap-frog in
521 time giving second order accuracy and more stability.
522
523 The essential change in the staggered algorithm is the calculation of
524 hydrostatic pressure which, in the context of the synchronous
525 algorithm involves replacing equation \ref{eq:phi-hyd-sync} with
526 \begin{displaymath}
527 \phi_{hyd}^n = \int b(\theta^{n+1},S^{n+1}) dr
528 \end{displaymath}
529 but the pressure gradient must also be taken out of the
530 Adams-Bashforth extrapolation. Also, retaining the integer time-levels,
531 $n$ and $n+1$, does not give a user the sense of where variables are
532 located in time. Instead, we re-write the entire algorithm,
533 \ref{eq:Gt-n-sync} to \ref{eq:v-n+1-sync}, annotating the
534 position in time of variables appropriately:
535 \begin{eqnarray}
536 G_{\theta,S}^{n-1/2} & = & G_{\theta,S} ( u^n, \theta^{n-1/2}, S^{n-1/2} )
537 \label{eq:Gt-n-staggered} \\
538 G_{\theta,S}^{(n)} & = & (3/2+\epsilon_{AB}) G_{\theta,S}^{n-1/2}-(1/2+\epsilon_{AB}) G_{\theta,S}^{n-3/2}
539 \label{eq:Gt-n+5-staggered} \\
540 (\theta^*,S^*) & = & (\theta^{n},S^{n}) + \Delta t G_{\theta,S}^{(n)}
541 \label{eq:tstar-staggered} \\
542 (\theta^{n+1/2},S^{n+1/2}) & = & {\cal L}^{-1}_{\theta,S} (\theta^*,S^*)
543 \label{eq:t-n+1-staggered} \\
544 \phi^{n+1/2}_{hyd} & = & \int b(\theta^{n+1/2},S^{n+1/2}) dr
545 \label{eq:phi-hyd-staggered} \\
546 \vec{\bf G}_{\vec{\bf v}}^{n} & = & \vec{\bf G}_{\vec{\bf v}} ( \vec{\bf v}^n )
547 \label{eq:Gv-n-staggered} \\
548 \vec{\bf G}_{\vec{\bf v}}^{(n+1/2)} & = & (3/2 + \epsilon_{AB} ) \vec{\bf G}_{\vec{\bf v}}^{n} - (1/2 + \epsilon_{AB} ) \vec{\bf G}_{\vec{\bf v}}^{n-1}
549 \label{eq:Gv-n+5-staggered} \\
550 \vec{\bf v}^{*} & = & \vec{\bf v}^{n} + \Delta t \left( \vec{\bf G}_{\vec{\bf v}}^{(n+1/2)} - \nabla \phi_{hyd}^{n+1/2} \right)
551 \label{eq:vstar-staggered} \\
552 \vec{\bf v}^{**} & = & {\cal L}_{\vec{\bf v}}^{-1} ( \vec{\bf v}^* )
553 \label{eq:vstarstar-staggered} \\
554 \eta^* & = & \epsilon_{fs} \left( \eta^{n} +P-E+R \right)- \Delta t
555 \nabla \cdot H \widehat{ \vec{\bf v}^{**} }
556 \label{eq:nstar-staggered} \\
557 \nabla \cdot g H \nabla \eta^{n+1} - \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2}
558 & = & - \frac{\eta^*}{\Delta t^2}
559 \label{eq:elliptic-staggered} \\
560 \vec{\bf v}^{n+1} & = & \vec{\bf v}^{*} - \Delta t g \nabla \eta^{n+1}
561 \label{eq:v-n+1-staggered}
562 \end{eqnarray}
563 The calling sequence is unchanged from
564 Fig.~\ref{fig:call-tree-adams-bashforth-sync}. The staggered algorithm
565 is activated with the run-time flag {\bf staggerTimeStep=.TRUE.} in
566 {\em PARM01} of {\em data}.
567
568 The only difficulty with this approach is apparent in equation
569 \ref{eq:Gt-n-staggered} and illustrated by the dotted arrow
570 connecting $u,v^n$ with $G_\theta^{n-1/2}$. The flow used to advect
571 tracers around is not naturally located in time. This could be avoided
572 by applying the Adams-Bashforth extrapolation to the tracer field
573 itself and advecting that around but this approach is not yet
574 available. We're not aware of any detrimental effect of this
575 feature. The difficulty lies mainly in interpretation of what
576 time-level variables and terms correspond to.
577
578
579 \section{Non-hydrostatic formulation}
580 \label{sect:non-hydrostatic}
581
582 The non-hydrostatic formulation re-introduces the full vertical
583 momentum equation and requires the solution of a 3-D elliptic
584 equations for non-hydrostatic pressure perturbation. We still
585 intergrate vertically for the hydrostatic pressure and solve a 2-D
586 elliptic equation for the surface pressure/elevation for this reduces
587 the amount of work needed to solve for the non-hydrostatic pressure.
588
589 The momentum equations are discretized in time as follows:
590 \begin{eqnarray}
591 \frac{1}{\Delta t} u^{n+1} + g \partial_x \eta^{n+1} + \partial_x \phi_{nh}^{n+1}
592 & = & \frac{1}{\Delta t} u^{n} + G_u^{(n+1/2)} \label{eq:discrete-time-u-nh} \\
593 \frac{1}{\Delta t} v^{n+1} + g \partial_y \eta^{n+1} + \partial_y \phi_{nh}^{n+1}
594 & = & \frac{1}{\Delta t} v^{n} + G_v^{(n+1/2)} \label{eq:discrete-time-v-nh} \\
595 \frac{1}{\Delta t} w^{n+1} + \partial_r \phi_{nh}^{n+1}
596 & = & \frac{1}{\Delta t} w^{n} + G_w^{(n+1/2)} \label{eq:discrete-time-w-nh} \\
597 \end{eqnarray}
598 which must satisfy the discrete-in-time depth integrated continuity,
599 equation~\ref{eq:discrete-time-backward-free-surface} and the local continuity equation
600 \begin{equation}
601 \partial_x u^{n+1} + \partial_y v^{n+1} + \partial_r w^{n+1} = 0
602 \label{eq:non-divergence-nh}
603 \end{equation}
604 As before, the explicit predictions for momentum are consolidated as:
605 \begin{eqnarray*}
606 u^* & = & u^n + \Delta t G_u^{(n+1/2)} \\
607 v^* & = & v^n + \Delta t G_v^{(n+1/2)} \\
608 w^* & = & w^n + \Delta t G_w^{(n+1/2)}
609 \end{eqnarray*}
610 but this time we introduce an intermediate step by splitting the
611 tendancy of the flow as follows:
612 \begin{eqnarray}
613 u^{n+1} = u^{**} - \Delta t \partial_x \phi_{nh}^{n+1}
614 & &
615 u^{**} = u^{*} - \Delta t g \partial_x \eta^{n+1} \\
616 v^{n+1} = v^{**} - \Delta t \partial_y \phi_{nh}^{n+1}
617 & &
618 v^{**} = v^{*} - \Delta t g \partial_y \eta^{n+1}
619 \end{eqnarray}
620 Substituting into the depth integrated continuity
621 (equation~\ref{eq:discrete-time-backward-free-surface}) gives
622 \begin{equation}
623 \partial_x H \partial_x \left( g \eta^{n+1} + \widehat{\phi}_{nh}^{n+1} \right)
624 +
625 \partial_y H \partial_y \left( g \eta^{n+1} + \widehat{\phi}_{nh}^{n+1} \right)
626 - \frac{\epsilon_{fs}\eta^*}{\Delta t^2}
627 = - \frac{\eta^*}{\Delta t^2}
628 \end{equation}
629 which is approximated by equation
630 \ref{eq:elliptic-backward-free-surface} on the basis that i)
631 $\phi_{nh}^{n+1}$ is not yet known and ii) $\nabla \widehat{\phi}_{nh}
632 << g \nabla \eta$. If \ref{eq:elliptic-backward-free-surface} is
633 solved accurately then the implication is that $\widehat{\phi}_{nh}
634 \approx 0$ so that thet non-hydrostatic pressure field does not drive
635 barotropic motion.
636
637 The flow must satisfy non-divergence
638 (equation~\ref{eq:non-divergence-nh}) locally, as well as depth
639 integrated, and this constraint is used to form a 3-D elliptic
640 equations for $\phi_{nh}^{n+1}$:
641 \begin{equation}
642 \partial_{xx} \phi_{nh}^{n+1} + \partial_{yy} \phi_{nh}^{n+1} +
643 \partial_{rr} \phi_{nh}^{n+1} =
644 \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*}
645 \end{equation}
646
647 The entire algorithm can be summarized as the sequential solution of
648 the following equations:
649 \begin{eqnarray}
650 u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2)} \label{eq:ustar-nh} \\
651 v^{*} & = & v^{n} + \Delta t G_v^{(n+1/2)} \label{eq:vstar-nh} \\
652 w^{*} & = & w^{n} + \Delta t G_w^{(n+1/2)} \label{eq:wstar-nh} \\
653 \eta^* & = & \epsilon_{fs} \left( \eta^{n} +P-E+R \right)- \Delta t
654 \partial_x H \widehat{u^{*}}
655 + \partial_y H \widehat{v^{*}}
656 \\
657 \partial_x g H \partial_x \eta^{n+1}
658 + \partial_y g H \partial_y \eta^{n+1}
659 - \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2}
660 & = &
661 - \frac{\eta^*}{\Delta t^2}
662 \label{eq:elliptic-nh}
663 \\
664 u^{**} & = & u^{*} - \Delta t g \partial_x \eta^{n+1} \label{eq:unx-nh}\\
665 v^{**} & = & v^{*} - \Delta t g \partial_y \eta^{n+1} \label{eq:vnx-nh}\\
666 \partial_{xx} \phi_{nh}^{n+1} + \partial_{yy} \phi_{nh}^{n+1} +
667 \partial_{rr} \phi_{nh}^{n+1} & = &
668 \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*} \\
669 u^{n+1} & = & u^{**} - \Delta t \partial_x \phi_{nh}^{n+1} \label{eq:un+1-nh}\\
670 v^{n+1} & = & v^{**} - \Delta t \partial_y \phi_{nh}^{n+1} \label{eq:vn+1-nh}\\
671 \partial_r w^{n+1} & = & - \partial_x u^{n+1} - \partial_y v^{n+1}
672 \end{eqnarray}
673 where the last equation is solved by vertically integrating for
674 $w^{n+1}$.
675
676
677
678 \section{Variants on the Free Surface}
679
680 We now describe the various formulations of the free-surface that
681 include non-linear forms, implicit in time using Crank-Nicholson,
682 explicit and [one day] split-explicit. First, we'll reiterate the
683 underlying algorithm but this time using the notation consistent with
684 the more general vertical coordinate $r$. The elliptic equation for
685 free-surface coordinate (units of $r$), corresponding to
686 \ref{eq:discrete-time-backward-free-surface}, and
687 assuming no non-hydrostatic effects ($\epsilon_{nh} = 0$) is:
688 \begin{eqnarray}
689 \epsilon_{fs} {\eta}^{n+1} -
690 {\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{fixed}) {\bf \nabla}_h b_s
691 {\eta}^{n+1} = {\eta}^*
692 \label{eq-solve2D}
693 \end{eqnarray}
694 where
695 \begin{eqnarray}
696 {\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
697 \Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o} \vec{\bf v}^* dr
698 \: + \: \epsilon_{fw} \Delta_t (P-E)^{n}
699 \label{eq-solve2D_rhs}
700 \end{eqnarray}
701
702 \fbox{ \begin{minipage}{4.75in}
703 {\em S/R SOLVE\_FOR\_PRESSURE} ({\em solve\_for\_pressure.F})
704
705 $u^*$: {\bf GuNm1} ({\em DYNVARS.h})
706
707 $v^*$: {\bf GvNm1} ({\em DYNVARS.h})
708
709 $\eta^*$: {\bf cg2d\_b} (\em SOLVE\_FOR\_PRESSURE.h)
710
711 $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)
712
713 \end{minipage} }
714
715
716 Once ${\eta}^{n+1}$ has been found, substituting into
717 \ref{eq:discrete-time-u,eq:discrete-time-v} yields $\vec{\bf v}^{n+1}$ if the model is
718 hydrostatic ($\epsilon_{nh}=0$):
719 $$
720 \vec{\bf v}^{n+1} = \vec{\bf v}^{*}
721 - \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
722 $$
723
724 This is known as the correction step. However, when the model is
725 non-hydrostatic ($\epsilon_{nh}=1$) we need an additional step and an
726 additional equation for $\phi'_{nh}$. This is obtained by substituting
727 \ref{eq:discrete-time-u-nh}, \ref{eq:discrete-time-v-nh} and \ref{eq:discrete-time-w-nh}
728 into continuity:
729 \begin{equation}
730 \left[ {\bf \nabla}_h^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1}
731 = \frac{1}{\Delta t} \left(
732 {\bf \nabla}_h \cdot \vec{\bf v}^{**} + \partial_r \dot{r}^* \right)
733 \end{equation}
734 where
735 \begin{displaymath}
736 \vec{\bf v}^{**} = \vec{\bf v}^* - \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
737 \end{displaymath}
738 Note that $\eta^{n+1}$ is also used to update the second RHS term
739 $\partial_r \dot{r}^* $ since
740 the vertical velocity at the surface ($\dot{r}_{surf}$)
741 is evaluated as $(\eta^{n+1} - \eta^n) / \Delta t$.
742
743 Finally, the horizontal velocities at the new time level are found by:
744 \begin{equation}
745 \vec{\bf v}^{n+1} = \vec{\bf v}^{**}
746 - \epsilon_{nh} \Delta t {\bf \nabla}_h {\phi'_{nh}}^{n+1}
747 \end{equation}
748 and the vertical velocity is found by integrating the continuity
749 equation vertically. Note that, for the convenience of the restart
750 procedure, the vertical integration of the continuity equation has
751 been moved to the beginning of the time step (instead of at the end),
752 without any consequence on the solution.
753
754 \fbox{ \begin{minipage}{4.75in}
755 {\em S/R CORRECTION\_STEP} ({\em correction\_step.F})
756
757 $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)
758
759 $\phi_{nh}^{n+1}$: {\bf phi\_nh} (\em DYNVARS.h)
760
761 $u^*$: {\bf GuNm1} ({\em DYNVARS.h})
762
763 $v^*$: {\bf GvNm1} ({\em DYNVARS.h})
764
765 $u^{n+1}$: {\bf uVel} ({\em DYNVARS.h})
766
767 $v^{n+1}$: {\bf vVel} ({\em DYNVARS.h})
768
769 \end{minipage} }
770
771
772
773 Regarding the implementation of the surface pressure solver, all
774 computation are done within the routine {\it SOLVE\_FOR\_PRESSURE} and
775 its dependent calls. The standard method to solve the 2D elliptic
776 problem (\ref{eq-solve2D}) uses the conjugate gradient method (routine
777 {\it CG2D}); the solver matrix and conjugate gradient operator are
778 only function of the discretized domain and are therefore evaluated
779 separately, before the time iteration loop, within {\it INI\_CG2D}.
780 The computation of the RHS $\eta^*$ is partly done in {\it
781 CALC\_DIV\_GHAT} and in {\it SOLVE\_FOR\_PRESSURE}.
782
783 The same method is applied for the non hydrostatic part, using a
784 conjugate gradient 3D solver ({\it CG3D}) that is initialized in {\it
785 INI\_CG3D}. The RHS terms of 2D and 3D problems are computed together
786 at the same point in the code.
787
788
789
790 \subsection{Crank-Nickelson barotropic time stepping}
791
792 The full implicit time stepping described previously is
793 unconditionally stable but damps the fast gravity waves, resulting in
794 a loss of potential energy. The modification presented now allows one
795 to combine an implicit part ($\beta,\gamma$) and an explicit part
796 ($1-\beta,1-\gamma$) for the surface pressure gradient ($\beta$) and
797 for the barotropic flow divergence ($\gamma$).
798 \\
799 For instance, $\beta=\gamma=1$ is the previous fully implicit scheme;
800 $\beta=\gamma=1/2$ is the non damping (energy conserving), unconditionally
801 stable, Crank-Nickelson scheme; $(\beta,\gamma)=(1,0)$ or $=(0,1)$
802 corresponds to the forward - backward scheme that conserves energy but is
803 only stable for small time steps.\\
804 In the code, $\beta,\gamma$ are defined as parameters, respectively
805 {\it implicSurfPress}, {\it implicDiv2DFlow}. They are read from
806 the main data file "{\it data}" and are set by default to 1,1.
807
808 Equations \ref{eq:ustar-backward-free-surface} --
809 \ref{eq:vn+1-backward-free-surface} are modified as follows:
810 $$
811 \frac{ \vec{\bf v}^{n+1} }{ \Delta t }
812 + {\bf \nabla}_h b_s [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ]
813 + \epsilon_{nh} {\bf \nabla}_h {\phi'_{nh}}^{n+1}
814 = \frac{ \vec{\bf v}^* }{ \Delta t }
815 $$
816 $$
817 \epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t}
818 + {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o}
819 [ \gamma \vec{\bf v}^{n+1} + (1-\gamma) \vec{\bf v}^{n}] dr
820 = \epsilon_{fw} (P-E)
821 $$
822 where:
823 \begin{eqnarray*}
824 \vec{\bf v}^* & = &
825 \vec{\bf v} ^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
826 + (\beta-1) \Delta t {\bf \nabla}_h b_s {\eta}^{n}
827 + \Delta t {\bf \nabla}_h {\phi'_{hyd}}^{(n+1/2)}
828 \\
829 {\eta}^* & = &
830 \epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta t (P-E)
831 - \Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o}
832 [ \gamma \vec{\bf v}^* + (1-\gamma) \vec{\bf v}^{n}] dr
833 \end{eqnarray*}
834 \\
835 In the hydrostatic case ($\epsilon_{nh}=0$), allowing us to find
836 ${\eta}^{n+1}$, thus:
837 $$
838 \epsilon_{fs} {\eta}^{n+1} -
839 {\bf \nabla}_h \cdot \beta\gamma \Delta t^2 b_s (R_o - R_{fixed})
840 {\bf \nabla}_h {\eta}^{n+1}
841 = {\eta}^*
842 $$
843 and then to compute (correction step):
844 $$
845 \vec{\bf v}^{n+1} = \vec{\bf v}^{*}
846 - \beta \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
847 $$
848
849 The non-hydrostatic part is solved as described previously.
850
851 Note that:
852 \begin{enumerate}
853 \item The non-hydrostatic part of the code has not yet been
854 updated, so that this option cannot be used with $(\beta,\gamma) \neq (1,1)$.
855 \item The stability criteria with Crank-Nickelson time stepping
856 for the pure linear gravity wave problem in cartesian coordinates is:
857 \begin{itemize}
858 \item $\beta + \gamma < 1$ : unstable
859 \item $\beta \geq 1/2$ and $ \gamma \geq 1/2$ : stable
860 \item $\beta + \gamma \geq 1$ : stable if
861 $$
862 c_{max}^2 (\beta - 1/2)(\gamma - 1/2) + 1 \geq 0
863 $$
864 $$
865 \mbox{with }~
866 %c^2 = 2 g H {\Delta t}^2
867 %(\frac{1-cos 2 \pi / k}{\Delta x^2}
868 %+\frac{1-cos 2 \pi / l}{\Delta y^2})
869 %$$
870 %Practically, the most stringent condition is obtained with $k=l=2$ :
871 %$$
872 c_{max} = 2 \Delta t \: \sqrt{g H} \:
873 \sqrt{ \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} }
874 $$
875 \end{itemize}
876 \end{enumerate}
877

  ViewVC Help
Powered by ViewVC 1.1.22