1 |
% $Header$ |
% $Header$ |
2 |
% $Name$ |
% $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 |
|
\input{part2/notation} |
14 |
|
|
15 |
|
\section{Time-stepping} |
16 |
|
\begin{rawhtml} |
17 |
|
<!-- CMIREDIR:time-stepping: --> |
18 |
|
\end{rawhtml} |
19 |
|
|
20 |
The equations of motion integrated by the model involve four |
The equations of motion integrated by the model involve four |
21 |
prognostic equations for flow, $u$ and $v$, temperature, $\theta$, and |
prognostic equations for flow, $u$ and $v$, temperature, $\theta$, and |
22 |
salt/moisture, $S$, and three diagnostic equations for vertical flow, |
salt/moisture, $S$, and three diagnostic equations for vertical flow, |
66 |
described in section \ref{sect:nonlinear-freesurface}. |
described in section \ref{sect:nonlinear-freesurface}. |
67 |
|
|
68 |
|
|
69 |
\section{Pressure method with rigid-lid} \label{sect:pressure-method-rigid-lid} |
\section{Pressure method with rigid-lid} |
70 |
|
\label{sect:pressure-method-rigid-lid} |
71 |
|
\begin{rawhtml} |
72 |
|
<!-- CMIREDIR:pressure_method_rigid_lid: --> |
73 |
|
\end{rawhtml} |
74 |
|
|
75 |
\begin{figure} |
\begin{figure} |
76 |
\begin{center} |
\begin{center} |
90 |
\begin{figure} |
\begin{figure} |
91 |
\begin{center} \fbox{ \begin{minipage}{4.5in} \begin{tabbing} |
\begin{center} \fbox{ \begin{minipage}{4.5in} \begin{tabbing} |
92 |
aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill |
aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill |
93 |
FORWARD\_STEP \\ |
\filelink{FORWARD\_STEP}{model-src-forward_step.F} \\ |
94 |
\> DYNAMICS \\ |
\> DYNAMICS \\ |
95 |
\>\> TIMESTEP \` $u^*$,$v^*$ (\ref{eq:ustar-rigid-lid},\ref{eq:vstar-rigid-lid}) \\ |
\>\> TIMESTEP \` $u^*$,$v^*$ (\ref{eq:ustar-rigid-lid},\ref{eq:vstar-rigid-lid}) \\ |
96 |
\> SOLVE\_FOR\_PRESSURE \\ |
\> SOLVE\_FOR\_PRESSURE \\ |
97 |
\>\> CALC\_DIV\_GHAT \` $H\widehat{u^*}$,$H\widehat{v^*}$ (\ref{eq:elliptic}) \\ |
\>\> CALC\_DIV\_GHAT \` $H\widehat{u^*}$,$H\widehat{v^*}$ (\ref{eq:elliptic}) \\ |
98 |
\>\> CG2D \` $\eta^{n+1}$ (\ref{eq:elliptic}) \\ |
\>\> CG2D \` $\eta^{n+1}$ (\ref{eq:elliptic}) \\ |
99 |
\> THE\_CORRECTION\_STEP \\ |
\> MOMENTUM\_CORRECTION\_STEP \\ |
100 |
\>\> CALC\_GRAD\_PHI\_SURF \` $\nabla \eta^{n+1}$ \\ |
\>\> CALC\_GRAD\_PHI\_SURF \` $\nabla \eta^{n+1}$ \\ |
101 |
\>\> CORRECTION\_STEP \` $u^{n+1}$,$v^{n+1}$ (\ref{eq:un+1-rigid-lid},\ref{eq:vn+1-rigid-lid}) |
\>\> CORRECTION\_STEP \` $u^{n+1}$,$v^{n+1}$ (\ref{eq:un+1-rigid-lid},\ref{eq:vn+1-rigid-lid}) |
102 |
\end{tabbing} \end{minipage} } \end{center} |
\end{tabbing} \end{minipage} } \end{center} |
103 |
\caption{Calling tree for the pressure method algorihtm} |
\caption{Calling tree for the pressure method algorithm |
104 |
|
(\filelink{FORWARD\_STEP}{model-src-forward_step.F})} |
105 |
\label{fig:call-tree-pressure-method} |
\label{fig:call-tree-pressure-method} |
106 |
\end{figure} |
\end{figure} |
107 |
|
|
183 |
\item |
\item |
184 |
the prognostic phase, equations \ref{eq:ustar-rigid-lid} and \ref{eq:vstar-rigid-lid}, |
the prognostic phase, equations \ref{eq:ustar-rigid-lid} and \ref{eq:vstar-rigid-lid}, |
185 |
stepping forward $u^n$ and $v^n$ to $u^{*}$ and $v^{*}$ is coded in |
stepping forward $u^n$ and $v^n$ to $u^{*}$ and $v^{*}$ is coded in |
186 |
{\em TIMESTEP.F} |
\filelink{TIMESTEP()}{model-src-timestep.F} |
187 |
\item |
\item |
188 |
the vertical integration, $H \widehat{u^*}$ and $H |
the vertical integration, $H \widehat{u^*}$ and $H |
189 |
\widehat{v^*}$, divergence and inversion of the elliptic operator in |
\widehat{v^*}$, divergence and inversion of the elliptic operator in |
190 |
equation \ref{eq:elliptic} is coded in {\em |
equation \ref{eq:elliptic} is coded in |
191 |
SOLVE\_FOR\_PRESSURE.F} |
\filelink{SOLVE\_FOR\_PRESSURE()}{model-src-solve_for_pressure.F} |
192 |
\item |
\item |
193 |
finally, the new flow field at time level $n+1$ given by equations |
finally, the new flow field at time level $n+1$ given by equations |
194 |
\ref{eq:un+1-rigid-lid} and \ref{eq:vn+1-rigid-lid} is calculated in {\em CORRECTION\_STEP.F}. |
\ref{eq:un+1-rigid-lid} and \ref{eq:vn+1-rigid-lid} is calculated in |
195 |
|
\filelink{CORRECTION\_STEP()}{model-src-correction_step.F}. |
196 |
\end{itemize} |
\end{itemize} |
197 |
The calling tree for these routines is given in |
The calling tree for these routines is given in |
198 |
Fig.~\ref{fig:call-tree-pressure-method}. |
Fig.~\ref{fig:call-tree-pressure-method}. |
199 |
|
|
200 |
|
|
201 |
|
%\paragraph{Need to discuss implicit viscosity somewhere:} |
202 |
\paragraph{Need to discuss implicit viscosity somewhere:} |
In general, the horizontal momentum time-stepping can contain some terms |
203 |
|
that are treated implicitly in time, |
204 |
|
such as the vertical viscosity when using the backward time-stepping scheme |
205 |
|
(\varlink{implicitViscosity}{implicitViscosity} {\it =.TRUE.}). |
206 |
|
The method used to solve those implicit terms is provided in |
207 |
|
section \ref{sect:implicit-backward-stepping}, and modifies |
208 |
|
equations \ref{eq:discrete-time-u} and \ref{eq:discrete-time-v} to |
209 |
|
give: |
210 |
\begin{eqnarray} |
\begin{eqnarray} |
211 |
\frac{1}{\Delta t} u^{n+1} - \partial_z A_v \partial_z u^{n+1} |
u^{n+1} - \Delta t \partial_z A_v \partial_z u^{n+1} |
212 |
+ g \partial_x \eta^{n+1} & = & \frac{1}{\Delta t} u^{n} + |
+ \Delta t g \partial_x \eta^{n+1} & = & u^{n} + \Delta t G_u^{(n+1/2)} |
|
G_u^{(n+1/2)} |
|
213 |
\\ |
\\ |
214 |
\frac{1}{\Delta t} v^{n+1} - \partial_z A_v \partial_z v^{n+1} |
v^{n+1} - \Delta t \partial_z A_v \partial_z v^{n+1} |
215 |
+ g \partial_y \eta^{n+1} & = & \frac{1}{\Delta t} v^{n} + G_v^{(n+1/2)} |
+ \Delta t g \partial_y \eta^{n+1} & = & v^{n} + \Delta t G_v^{(n+1/2)} |
216 |
\end{eqnarray} |
\end{eqnarray} |
217 |
|
|
218 |
|
|
219 |
\section{Pressure method with implicit linear free-surface} |
\section{Pressure method with implicit linear free-surface} |
220 |
\label{sect:pressure-method-linear-backward} |
\label{sect:pressure-method-linear-backward} |
221 |
|
\begin{rawhtml} |
222 |
|
<!-- CMIREDIR:pressure_method_linear_backward: --> |
223 |
|
\end{rawhtml} |
224 |
|
|
225 |
The rigid-lid approximation filters out external gravity waves |
The rigid-lid approximation filters out external gravity waves |
226 |
subsequently modifying the dispersion relation of barotropic Rossby |
subsequently modifying the dispersion relation of barotropic Rossby |
232 |
of the free-surface equation which can be written: |
of the free-surface equation which can be written: |
233 |
\begin{equation} |
\begin{equation} |
234 |
\partial_t \eta + \partial_x H \widehat{u} + \partial_y H \widehat{v} = P-E+R |
\partial_t \eta + \partial_x H \widehat{u} + \partial_y H \widehat{v} = P-E+R |
235 |
\label{eq:linear-free-surface=P-E+R} |
\label{eq:linear-free-surface=P-E} |
236 |
\end{equation} |
\end{equation} |
237 |
which differs from the depth integrated continuity equation with |
which differs from the depth integrated continuity equation with |
238 |
rigid-lid (\ref{eq:rigid-lid-continuity}) by the time-dependent term |
rigid-lid (\ref{eq:rigid-lid-continuity}) by the time-dependent term |
240 |
|
|
241 |
Equation \ref{eq:discrete-time-cont-rigid-lid} in the rigid-lid |
Equation \ref{eq:discrete-time-cont-rigid-lid} in the rigid-lid |
242 |
pressure method is then replaced by the time discretization of |
pressure method is then replaced by the time discretization of |
243 |
\ref{eq:linear-free-surface=P-E+R} which is: |
\ref{eq:linear-free-surface=P-E} which is: |
244 |
\begin{equation} |
\begin{equation} |
245 |
\eta^{n+1} |
\eta^{n+1} |
246 |
+ \Delta t \partial_x H \widehat{u^{n+1}} |
+ \Delta t \partial_x H \widehat{u^{n+1}} |
247 |
+ \Delta t \partial_y H \widehat{v^{n+1}} |
+ \Delta t \partial_y H \widehat{v^{n+1}} |
248 |
= |
= |
249 |
\eta^{n} |
\eta^{n} |
250 |
+ \Delta t ( P - E + R ) |
+ \Delta t ( P - E ) |
251 |
\label{eq:discrete-time-backward-free-surface} |
\label{eq:discrete-time-backward-free-surface} |
252 |
\end{equation} |
\end{equation} |
253 |
where the use of flow at time level $n+1$ makes the method implicit |
where the use of flow at time level $n+1$ makes the method implicit |
254 |
and backward in time. The is the preferred scheme since it still |
and backward in time. This is the preferred scheme since it still |
255 |
filters the fast unresolved wave motions by damping them. A centered |
filters the fast unresolved wave motions by damping them. A centered |
256 |
scheme, such as Crank-Nicholson, would alias the energy of the fast |
scheme, such as Crank-Nicholson (see section \ref{sect:freesurf-CrankNick}), |
257 |
modes onto slower modes of motion. |
would alias the energy of the fast modes onto slower modes of motion. |
258 |
|
|
259 |
As for the rigid-lid pressure method, equations |
As for the rigid-lid pressure method, equations |
260 |
\ref{eq:discrete-time-u}, \ref{eq:discrete-time-v} and |
\ref{eq:discrete-time-u}, \ref{eq:discrete-time-v} and |
262 |
\begin{eqnarray} |
\begin{eqnarray} |
263 |
u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2)} \label{eq:ustar-backward-free-surface} \\ |
u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2)} \label{eq:ustar-backward-free-surface} \\ |
264 |
v^{*} & = & v^{n} + \Delta t G_v^{(n+1/2)} \label{eq:vstar-backward-free-surface} \\ |
v^{*} & = & v^{n} + \Delta t G_v^{(n+1/2)} \label{eq:vstar-backward-free-surface} \\ |
265 |
\eta^* & = & \epsilon_{fs} \left( \eta^{n} +P-E+R \right)- \Delta t |
\eta^* & = & \epsilon_{fs} \left( \eta^{n} + \Delta t (P-E) \right)- \Delta t |
266 |
\partial_x H \widehat{u^{*}} |
\partial_x H \widehat{u^{*}} |
267 |
+ \partial_y H \widehat{v^{*}} |
+ \partial_y H \widehat{v^{*}} |
268 |
\\ |
\\ |
269 |
\partial_x g H \partial_x \eta^{n+1} |
\partial_x g H \partial_x \eta^{n+1} |
270 |
+ \partial_y g H \partial_y \eta^{n+1} |
& + & \partial_y g H \partial_y \eta^{n+1} |
271 |
- \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2} |
- \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2} |
272 |
& = & |
= |
273 |
- \frac{\eta^*}{\Delta t^2} |
- \frac{\eta^*}{\Delta t^2} |
274 |
\label{eq:elliptic-backward-free-surface} |
\label{eq:elliptic-backward-free-surface} |
275 |
\\ |
\\ |
285 |
the flow to be divergent and for the surface pressure/elevation to |
the flow to be divergent and for the surface pressure/elevation to |
286 |
respond on a finite time-scale (as opposed to instantly). To recover |
respond on a finite time-scale (as opposed to instantly). To recover |
287 |
the rigid-lid formulation, we introduced a switch-like parameter, |
the rigid-lid formulation, we introduced a switch-like parameter, |
288 |
$\epsilon_{fs}$, which selects between the free-surface and rigid-lid; |
$\epsilon_{fs}$ (\varlink{freesurfFac}{freesurfFac}), |
289 |
|
which selects between the free-surface and rigid-lid; |
290 |
$\epsilon_{fs}=1$ allows the free-surface to evolve; $\epsilon_{fs}=0$ |
$\epsilon_{fs}=1$ allows the free-surface to evolve; $\epsilon_{fs}=0$ |
291 |
imposes the rigid-lid. The evolution in time and location of variables |
imposes the rigid-lid. The evolution in time and location of variables |
292 |
is exactly as it was for the rigid-lid model so that |
is exactly as it was for the rigid-lid model so that |
298 |
|
|
299 |
\section{Explicit time-stepping: Adams-Bashforth} |
\section{Explicit time-stepping: Adams-Bashforth} |
300 |
\label{sect:adams-bashforth} |
\label{sect:adams-bashforth} |
301 |
|
\begin{rawhtml} |
302 |
|
<!-- CMIREDIR:adams_bashforth: --> |
303 |
|
\end{rawhtml} |
304 |
|
|
305 |
In describing the the pressure method above we deferred describing the |
In describing the the pressure method above we deferred describing the |
306 |
time discretization of the explicit terms. We have historically used |
time discretization of the explicit terms. We have historically used |
316 |
\> THERMODYNAMICS \\ |
\> THERMODYNAMICS \\ |
317 |
\>\> CALC\_GT \\ |
\>\> CALC\_GT \\ |
318 |
\>\>\> GAD\_CALC\_RHS \` $G_\theta^n = G_\theta( u, \theta^n )$ \\ |
\>\>\> GAD\_CALC\_RHS \` $G_\theta^n = G_\theta( u, \theta^n )$ \\ |
319 |
\>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\ |
\>either\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\ |
320 |
\>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:adams-bashforth2}) \\ |
\>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:adams-bashforth2}) \\ |
321 |
|
\>or\>\> EXTERNAL\_FORCING \` $G_\theta^{(n+1/2)} = G_\theta^{(n+1/2)} + {\cal Q}$ \\ |
322 |
\>\> TIMESTEP\_TRACER \` $\tau^*$ (\ref{eq:taustar}) \\ |
\>\> TIMESTEP\_TRACER \` $\tau^*$ (\ref{eq:taustar}) \\ |
323 |
\>\> IMPLDIFF \` $\tau^{(n+1)}$ (\ref{eq:tau-n+1-implicit}) |
\>\> IMPLDIFF \` $\tau^{(n+1)}$ (\ref{eq:tau-n+1-implicit}) |
324 |
\end{tabbing} \end{minipage} } \end{center} |
\end{tabbing} \end{minipage} } \end{center} |
325 |
\caption{ |
\caption{ |
326 |
Calling tree for the Adams-Bashforth time-stepping of temperature with |
Calling tree for the Adams-Bashforth time-stepping of temperature with |
327 |
implicit diffusion.} |
implicit diffusion. |
328 |
|
(\filelink{THERMODYNAMICS}{model-src-thermodynamics.F}, |
329 |
|
\filelink{ADAMS\_BASHFORTH2}{model-src-adams_bashforth2.F})} |
330 |
\label{fig:call-tree-adams-bashforth} |
\label{fig:call-tree-adams-bashforth} |
331 |
\end{figure} |
\end{figure} |
332 |
|
|
357 |
and forcing evolves smoothly. Problems can, and do, arise when forcing |
and forcing evolves smoothly. Problems can, and do, arise when forcing |
358 |
or motions are high frequency and this corresponds to a reduced |
or motions are high frequency and this corresponds to a reduced |
359 |
stability compared to a simple forward time-stepping of such terms. |
stability compared to a simple forward time-stepping of such terms. |
360 |
|
The model offers the possibility to leave the tracer and momentum |
361 |
|
forcing terms and the dissipation terms outside the |
362 |
|
Adams-Bashforth extrapolation, by turning off the logical flags |
363 |
|
\varlink{forcing\_In\_AB}{forcing_In_AB} |
364 |
|
(parameter file {\em data}, namelist {\em PARM01}, default value = True). |
365 |
|
(\varlink{tracForcingOutAB}{tracForcingOutAB}, default=0, |
366 |
|
\varlink{momForcingOutAB}{momForcingOutAB}, default=0) |
367 |
|
and \varlink{momDissip\_In\_AB}{momDissip_In_AB} |
368 |
|
(parameter file {\em data}, namelist {\em PARM01}, default value = True). |
369 |
|
respectively. |
370 |
|
|
371 |
A stability analysis for an oscillation equation should be given at this point. |
A stability analysis for an oscillation equation should be given at this point. |
372 |
\marginpar{AJA needs to find his notes on this...} |
\marginpar{AJA needs to find his notes on this...} |
374 |
A stability analysis for a relaxation equation should be given at this point. |
A stability analysis for a relaxation equation should be given at this point. |
375 |
\marginpar{...and for this too.} |
\marginpar{...and for this too.} |
376 |
|
|
377 |
|
\begin{figure} |
378 |
|
\begin{center} |
379 |
|
\resizebox{5.5in}{!}{\includegraphics{part2/oscil+damp_AB2.eps}} |
380 |
|
\end{center} |
381 |
|
\caption{ |
382 |
|
Oscillatory and damping response of |
383 |
|
quasi-second order Adams-Bashforth scheme for different values |
384 |
|
of the $\epsilon_{AB}$ parameter (0., 0.1, 0.25, from top to bottom) |
385 |
|
The analytical solution (in black), the physical mode (in blue) |
386 |
|
and the numerical mode (in red) are represented with a CFL |
387 |
|
step of 0.1. |
388 |
|
The left column represents the oscillatory response |
389 |
|
on the complex plane for CFL ranging from 0.1 up to 0.9. |
390 |
|
The right column represents the damping response amplitude |
391 |
|
(y-axis) function of the CFL (x-axis). |
392 |
|
} |
393 |
|
\label{fig:adams-bashforth-respons} |
394 |
|
\end{figure} |
395 |
|
|
396 |
|
|
397 |
|
|
398 |
\section{Implicit time-stepping: backward method} |
\section{Implicit time-stepping: backward method} |
399 |
|
\label{sect:implicit-backward-stepping} |
400 |
|
\begin{rawhtml} |
401 |
|
<!-- CMIREDIR:implicit_time-stepping_backward: --> |
402 |
|
\end{rawhtml} |
403 |
|
|
404 |
Vertical diffusion and viscosity can be treated implicitly in time |
Vertical diffusion and viscosity can be treated implicitly in time |
405 |
using the backward method which is an intrinsic scheme. For tracers, |
using the backward method which is an intrinsic scheme. |
406 |
|
Recently, the option to treat the vertical advection |
407 |
|
implicitly has been added, but not yet tested; therefore, |
408 |
|
the description hereafter is limited to diffusion and viscosity. |
409 |
|
For tracers, |
410 |
the time discretized equation is: |
the time discretized equation is: |
411 |
\begin{equation} |
\begin{equation} |
412 |
\tau^{n+1} - \Delta t \partial_r \kappa_v \partial_r \tau^{n+1} = |
\tau^{n+1} - \Delta t \partial_r \kappa_v \partial_r \tau^{n+1} = |
437 |
stepping forward a tracer variable such as temperature. |
stepping forward a tracer variable such as temperature. |
438 |
|
|
439 |
In order to fit within the pressure method, the implicit viscosity |
In order to fit within the pressure method, the implicit viscosity |
440 |
must not alter the barotropic flow. In other words, it can on ly |
must not alter the barotropic flow. In other words, it can only |
441 |
redistribute momentum in the vertical. The upshot of this is that |
redistribute momentum in the vertical. The upshot of this is that |
442 |
although vertical viscosity may be backward implicit and |
although vertical viscosity may be backward implicit and |
443 |
unconditionally stable, no-slip boundary conditions may not be made |
unconditionally stable, no-slip boundary conditions may not be made |
445 |
|
|
446 |
\section{Synchronous time-stepping: variables co-located in time} |
\section{Synchronous time-stepping: variables co-located in time} |
447 |
\label{sect:adams-bashforth-sync} |
\label{sect:adams-bashforth-sync} |
448 |
|
\begin{rawhtml} |
449 |
|
<!-- CMIREDIR:adams_bashforth_sync: --> |
450 |
|
\end{rawhtml} |
451 |
|
|
452 |
\begin{figure} |
\begin{figure} |
453 |
\begin{center} |
\begin{center} |
468 |
\end{figure} |
\end{figure} |
469 |
|
|
470 |
\begin{figure} |
\begin{figure} |
471 |
\begin{center} \fbox{ \begin{minipage}{4.5in} \begin{tabbing} |
\begin{center} \fbox{ \begin{minipage}{4.7in} \begin{tabbing} |
472 |
aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill |
aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill |
473 |
FORWARD\_STEP \\ |
FORWARD\_STEP \\ |
474 |
|
\>\> EXTERNAL\_FIELDS\_LOAD\\ |
475 |
|
\>\> DO\_ATMOSPHERIC\_PHYS \\ |
476 |
|
\>\> DO\_OCEANIC\_PHYS \\ |
477 |
\> THERMODYNAMICS \\ |
\> THERMODYNAMICS \\ |
478 |
\>\> CALC\_GT \\ |
\>\> CALC\_GT \\ |
479 |
\>\>\> GAD\_CALC\_RHS \` $G_\theta^n = G_\theta( u, \theta^n )$ (\ref{eq:Gt-n-sync})\\ |
\>\>\> GAD\_CALC\_RHS \` $G_\theta^n = G_\theta( u, \theta^n )$ (\ref{eq:Gt-n-sync})\\ |
480 |
\>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\ |
\>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\ |
481 |
\>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:Gt-n+5-sync}) \\ |
\>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:Gt-n+5-sync}) \\ |
482 |
\>\> TIMESTEP\_TRACER \` $\tau^*$ (\ref{eq:tstar-sync}) \\ |
\>\> TIMESTEP\_TRACER \` $\theta^*$ (\ref{eq:tstar-sync}) \\ |
483 |
\>\> IMPLDIFF \` $\tau^{(n+1)}$ (\ref{eq:t-n+1-sync}) \\ |
\>\> IMPLDIFF \` $\theta^{(n+1)}$ (\ref{eq:t-n+1-sync}) \\ |
484 |
\> DYNAMICS \\ |
\> DYNAMICS \\ |
485 |
\>\> CALC\_PHI\_HYD \` $\phi_{hyd}^n$ (\ref{eq:phi-hyd-sync}) \\ |
\>\> CALC\_PHI\_HYD \` $\phi_{hyd}^n$ (\ref{eq:phi-hyd-sync}) \\ |
486 |
\>\> MOM\_FLUXFORM or MOM\_VECINV \` $G_{\vec{\bf v}}^n$ (\ref{eq:Gv-n-sync})\\ |
\>\> MOM\_FLUXFORM or MOM\_VECINV \` $G_{\vec{\bf v}}^n$ (\ref{eq:Gv-n-sync})\\ |
487 |
\>\> TIMESTEP \` $\vec{\bf v}^*$ (\ref{eq:Gv-n+5-sync}, \ref{eq:vstar-sync}) \\ |
\>\> TIMESTEP \` $\vec{\bf v}^*$ (\ref{eq:Gv-n+5-sync}, \ref{eq:vstar-sync}) \\ |
488 |
\>\> IMPLDIFF \` $\vec{\bf v}^{**}$ (\ref{eq:vstarstar-sync}) \\ |
\>\> IMPLDIFF \` $\vec{\bf v}^{**}$ (\ref{eq:vstarstar-sync}) \\ |
489 |
|
\> UPDATE\_R\_STAR or UPDATE\_SURF\_DR \` (NonLin-FS only)\\ |
490 |
\> SOLVE\_FOR\_PRESSURE \\ |
\> SOLVE\_FOR\_PRESSURE \\ |
491 |
\>\> CALC\_DIV\_GHAT \` $\eta^*$ (\ref{eq:nstar-sync}) \\ |
\>\> CALC\_DIV\_GHAT \` $\eta^*$ (\ref{eq:nstar-sync}) \\ |
492 |
\>\> CG2D \` $\eta^{n+1}$ (\ref{eq:elliptic-sync}) \\ |
\>\> CG2D \` $\eta^{n+1}$ (\ref{eq:elliptic-sync}) \\ |
493 |
\> THE\_CORRECTION\_STEP \\ |
\> MOMENTUM\_CORRECTION\_STEP \\ |
494 |
\>\> CALC\_GRAD\_PHI\_SURF \` $\nabla \eta^{n+1}$ \\ |
\>\> CALC\_GRAD\_PHI\_SURF \` $\nabla \eta^{n+1}$ \\ |
495 |
\>\> CORRECTION\_STEP \` $u^{n+1}$,$v^{n+1}$ (\ref{eq:v-n+1-sync}) |
\>\> CORRECTION\_STEP \` $u^{n+1}$,$v^{n+1}$ (\ref{eq:v-n+1-sync})\\ |
496 |
|
\> TRACERS\_CORRECTION\_STEP \\ |
497 |
|
\>\> CYCLE\_TRACER \` $\theta^{n+1}$ \\ |
498 |
|
\>\> FILTER \` Shapiro Filter, Zonal Filter (FFT) \\ |
499 |
|
\>\> CONVECTIVE\_ADJUSTMENT \` \\ |
500 |
\end{tabbing} \end{minipage} } \end{center} |
\end{tabbing} \end{minipage} } \end{center} |
501 |
\caption{ |
\caption{ |
502 |
Calling tree for the overall synchronous algorithm using |
Calling tree for the overall synchronous algorithm using |
503 |
Adams-Bashforth time-stepping.} |
Adams-Bashforth time-stepping. |
504 |
|
The place where the model geometry |
505 |
|
({\bf hFac} factors) is updated is added here but is only relevant |
506 |
|
for the non-linear free-surface algorithm. |
507 |
|
For completeness, the external forcing, |
508 |
|
ocean and atmospheric physics have been added, although they are mainly |
509 |
|
optional} |
510 |
\label{fig:call-tree-adams-bashforth-sync} |
\label{fig:call-tree-adams-bashforth-sync} |
511 |
\end{figure} |
\end{figure} |
512 |
|
|
535 |
\label{eq:vstar-sync} \\ |
\label{eq:vstar-sync} \\ |
536 |
\vec{\bf v}^{**} & = & {\cal L}_{\vec{\bf v}}^{-1} ( \vec{\bf v}^* ) |
\vec{\bf v}^{**} & = & {\cal L}_{\vec{\bf v}}^{-1} ( \vec{\bf v}^* ) |
537 |
\label{eq:vstarstar-sync} \\ |
\label{eq:vstarstar-sync} \\ |
538 |
\eta^* & = & \epsilon_{fs} \left( \eta^{n} +P-E+R \right)- \Delta t |
\eta^* & = & \epsilon_{fs} \left( \eta^{n} + \Delta t (P-E) \right)- \Delta t |
539 |
\nabla \cdot H \widehat{ \vec{\bf v}^{**} } |
\nabla \cdot H \widehat{ \vec{\bf v}^{**} } |
540 |
\label{eq:nstar-sync} \\ |
\label{eq:nstar-sync} \\ |
541 |
\nabla \cdot g H \nabla \eta^{n+1} - \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2} |
\nabla \cdot g H \nabla \eta^{n+1} & - & \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2} |
542 |
& = & - \frac{\eta^*}{\Delta t^2} |
~ = ~ - \frac{\eta^*}{\Delta t^2} |
543 |
\label{eq:elliptic-sync} \\ |
\label{eq:elliptic-sync} \\ |
544 |
\vec{\bf v}^{n+1} & = & \vec{\bf v}^{*} - \Delta t g \nabla \eta^{n+1} |
\vec{\bf v}^{n+1} & = & \vec{\bf v}^{*} - \Delta t g \nabla \eta^{n+1} |
545 |
\label{eq:v-n+1-sync} |
\label{eq:v-n+1-sync} |
558 |
surface pressure gradient terms, corresponding to equations |
surface pressure gradient terms, corresponding to equations |
559 |
\ref{eq:Gv-n-sync} to \ref{eq:v-n+1-sync}. |
\ref{eq:Gv-n-sync} to \ref{eq:v-n+1-sync}. |
560 |
These operations are carried out in subroutines {\em DYNAMCIS}, {\em |
These operations are carried out in subroutines {\em DYNAMCIS}, {\em |
561 |
SOLVE\_FOR\_PRESSURE} and {\em THE\_CORRECTION\_STEP}. This, then, |
SOLVE\_FOR\_PRESSURE} and {\em MOMENTUM\_CORRECTION\_STEP}. This, then, |
562 |
represents an entire algorithm for stepping forward the model one |
represents an entire algorithm for stepping forward the model one |
563 |
time-step. The corresponding calling tree is given in |
time-step. The corresponding calling tree is given in |
564 |
\ref{fig:call-tree-adams-bashforth-sync}. |
\ref{fig:call-tree-adams-bashforth-sync}. |
565 |
|
|
566 |
\section{Staggered baroclinic time-stepping} |
\section{Staggered baroclinic time-stepping} |
567 |
\label{sect:adams-bashforth-staggered} |
\label{sect:adams-bashforth-staggered} |
568 |
|
\begin{rawhtml} |
569 |
|
<!-- CMIREDIR:adams_bashforth_staggered: --> |
570 |
|
\end{rawhtml} |
571 |
|
|
572 |
\begin{figure} |
\begin{figure} |
573 |
\begin{center} |
\begin{center} |
576 |
\caption{ |
\caption{ |
577 |
A schematic of the explicit Adams-Bashforth and implicit time-stepping |
A schematic of the explicit Adams-Bashforth and implicit time-stepping |
578 |
phases of the algorithm but with staggering in time of thermodynamic |
phases of the algorithm but with staggering in time of thermodynamic |
579 |
variables with the flow. Explicit thermodynamics tendencies are |
variables with the flow. |
580 |
evaluated at time level $n-1/2$ as a function of the thermodynamics |
Explicit momentum tendencies are evaluated at time level $n-1/2$ as a |
581 |
state at that time level $n$ and flow at time $n$ (dotted arrow). The |
function of the flow field at that time level $n-1/2$. |
582 |
explicit tendency from the previous time level, $n-3/2$, is used to |
The explicit tendency from the previous time level, $n-3/2$, is used to |
583 |
extrapolate tendencies to $n$ (dashed arrow). This extrapolated |
extrapolate tendencies to $n$ (dashed arrow). |
584 |
tendency allows thermo-dynamics variables to be stably integrated |
The hydrostatic pressure/geo-potential $\phi_{hyd}$ is evaluated directly |
585 |
forward-in-time to render an estimate ($*$-variables) at the $n+1/2$ |
at time level $n$ (vertical arrows) and used with the extrapolated tendencies |
586 |
time level (solid arc-arrow). The implicit-in-time operator ${\cal |
to step forward the flow variables from $n-1/2$ to $n+1/2$ (solid arc-arrow). |
587 |
L_{\theta,S}}$ is solved to yield the thermodynamic variables at time |
The implicit-in-time operator ${\cal L}_{\bf u,v}$ (vertical arrows) is |
588 |
level $n+1/2$. These are then used to calculate the hydrostatic |
then applied to the previous estimation of the the flow field ($*$-variables) |
589 |
pressure/geo-potential, $\phi_{hyd}$ (vertical arrows). The |
and yields to the two velocity components $u,v$ at time level $n+1/2$. |
590 |
hydrostatic pressure gradient is evaluated directly an time level |
These are then used to calculate the advection term (dashed arc-arrow) |
591 |
$n+1/2$ in stepping forward the flow variables from $n$ to $n+1$ |
of the thermo-dynamics tendencies at time step $n$. |
592 |
(solid arc-arrow). } |
The extrapolated thermodynamics tendency, from time level $n-1$ and $n$ |
593 |
|
to $n+1/2$, allows thermodynamic variables to be stably integrated |
594 |
|
forward-in-time (solid arc-arrow) up to time level $n+1$. |
595 |
|
} |
596 |
\label{fig:adams-bashforth-staggered} |
\label{fig:adams-bashforth-staggered} |
597 |
\end{figure} |
\end{figure} |
598 |
|
|
602 |
thermodynamic variables with the flow |
thermodynamic variables with the flow |
603 |
variables. Fig.~\ref{fig:adams-bashforth-staggered} illustrates the |
variables. Fig.~\ref{fig:adams-bashforth-staggered} illustrates the |
604 |
staggering and algorithm. The key difference between this and |
staggering and algorithm. The key difference between this and |
605 |
Fig.~\ref{fig:adams-bashforth-sync} is that the new thermodynamics |
Fig.~\ref{fig:adams-bashforth-sync} is that the thermodynamic variables |
606 |
fields are used to compute the hydrostatic pressure at time level |
are solved after the dynamics, using the recently updated flow field. |
607 |
$n+1/2$. The essentially allows the gravity wave terms to leap-frog in |
This essentially allows the gravity wave terms to leap-frog in |
608 |
time giving second order accuracy and more stability. |
time giving second order accuracy and more stability. |
609 |
|
|
610 |
The essential change in the staggered algorithm is the calculation of |
The essential change in the staggered algorithm is that the |
611 |
hydrostatic pressure which, in the context of the synchronous |
thermodynamics solver is delayed from half a time step, |
612 |
algorithm involves replacing equation \ref{eq:phi-hyd-sync} with |
allowing the use of the most recent velocities to compute |
613 |
\begin{displaymath} |
the advection terms. Once the thermodynamics fields are |
614 |
\phi_{hyd}^n = \int b(\theta^{n+1},S^{n+1}) dr |
updated, the hydrostatic pressure is computed |
615 |
\end{displaymath} |
to step forwrad the dynamics. |
616 |
but the pressure gradient must also be taken out of the |
Note that the pressure gradient must also be taken out of the |
617 |
Adams-Bashforth extrapolation. Also, retaining the integer time-levels, |
Adams-Bashforth extrapolation. Also, retaining the integer time-levels, |
618 |
$n$ and $n+1$, does not give a user the sense of where variables are |
$n$ and $n+1$, does not give a user the sense of where variables are |
619 |
located in time. Instead, we re-write the entire algorithm, |
located in time. Instead, we re-write the entire algorithm, |
620 |
\ref{eq:Gt-n-sync} to \ref{eq:v-n+1-sync}, annotating the |
\ref{eq:Gt-n-sync} to \ref{eq:v-n+1-sync}, annotating the |
621 |
position in time of variables appropriately: |
position in time of variables appropriately: |
622 |
\begin{eqnarray} |
\begin{eqnarray} |
623 |
G_{\theta,S}^{n-1/2} & = & G_{\theta,S} ( u^n, \theta^{n-1/2}, S^{n-1/2} ) |
\phi^{n}_{hyd} & = & \int b(\theta^{n},S^{n}) dr |
|
\label{eq:Gt-n-staggered} \\ |
|
|
G_{\theta,S}^{(n)} & = & (3/2+\epsilon_{AB}) G_{\theta,S}^{n-1/2}-(1/2+\epsilon_{AB}) G_{\theta,S}^{n-3/2} |
|
|
\label{eq:Gt-n+5-staggered} \\ |
|
|
(\theta^*,S^*) & = & (\theta^{n},S^{n}) + \Delta t G_{\theta,S}^{(n)} |
|
|
\label{eq:tstar-staggered} \\ |
|
|
(\theta^{n+1/2},S^{n+1/2}) & = & {\cal L}^{-1}_{\theta,S} (\theta^*,S^*) |
|
|
\label{eq:t-n+1-staggered} \\ |
|
|
\phi^{n+1/2}_{hyd} & = & \int b(\theta^{n+1/2},S^{n+1/2}) dr |
|
624 |
\label{eq:phi-hyd-staggered} \\ |
\label{eq:phi-hyd-staggered} \\ |
625 |
\vec{\bf G}_{\vec{\bf v}}^{n} & = & \vec{\bf G}_{\vec{\bf v}} ( \vec{\bf v}^n ) |
\vec{\bf G}_{\vec{\bf v}}^{n-1/2} & = & \vec{\bf G}_{\vec{\bf v}} ( \vec{\bf v}^{n-1/2} ) |
626 |
\label{eq:Gv-n-staggered} \\ |
\label{eq:Gv-n-staggered} \\ |
627 |
\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} |
\vec{\bf G}_{\vec{\bf v}}^{(n)} & = & (3/2 + \epsilon_{AB} ) \vec{\bf G}_{\vec{\bf v}}^{n-1/2} - (1/2 + \epsilon_{AB} ) \vec{\bf G}_{\vec{\bf v}}^{n-3/2} |
628 |
\label{eq:Gv-n+5-staggered} \\ |
\label{eq:Gv-n+5-staggered} \\ |
629 |
\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) |
\vec{\bf v}^{*} & = & \vec{\bf v}^{n-1/2} + \Delta t \left( \vec{\bf G}_{\vec{\bf v}}^{(n)} - \nabla \phi_{hyd}^{n} \right) |
630 |
\label{eq:vstar-staggered} \\ |
\label{eq:vstar-staggered} \\ |
631 |
\vec{\bf v}^{**} & = & {\cal L}_{\vec{\bf v}}^{-1} ( \vec{\bf v}^* ) |
\vec{\bf v}^{**} & = & {\cal L}_{\vec{\bf v}}^{-1} ( \vec{\bf v}^* ) |
632 |
\label{eq:vstarstar-staggered} \\ |
\label{eq:vstarstar-staggered} \\ |
633 |
\eta^* & = & \epsilon_{fs} \left( \eta^{n} +P-E+R \right)- \Delta t |
\eta^* & = & \epsilon_{fs} \left( \eta^{n-1/2} + \Delta t (P-E)^n \right)- \Delta t |
634 |
\nabla \cdot H \widehat{ \vec{\bf v}^{**} } |
\nabla \cdot H \widehat{ \vec{\bf v}^{**} } |
635 |
\label{eq:nstar-staggered} \\ |
\label{eq:nstar-staggered} \\ |
636 |
\nabla \cdot g H \nabla \eta^{n+1} - \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2} |
\nabla \cdot g H \nabla \eta^{n+1/2} & - & \frac{\epsilon_{fs} \eta^{n+1/2}}{\Delta t^2} |
637 |
& = & - \frac{\eta^*}{\Delta t^2} |
~ = ~ - \frac{\eta^*}{\Delta t^2} |
638 |
\label{eq:elliptic-staggered} \\ |
\label{eq:elliptic-staggered} \\ |
639 |
\vec{\bf v}^{n+1} & = & \vec{\bf v}^{*} - \Delta t g \nabla \eta^{n+1} |
\vec{\bf v}^{n+1/2} & = & \vec{\bf v}^{*} - \Delta t g \nabla \eta^{n+1/2} |
640 |
\label{eq:v-n+1-staggered} |
\label{eq:v-n+1-staggered} \\ |
641 |
|
G_{\theta,S}^{n} & = & G_{\theta,S} ( u^{n+1/2}, \theta^{n}, S^{n} ) |
642 |
|
\label{eq:Gt-n-staggered} \\ |
643 |
|
G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsilon_{AB}) G_{\theta,S}^{n}-(1/2+\epsilon_{AB}) G_{\theta,S}^{n-1} |
644 |
|
\label{eq:Gt-n+5-staggered} \\ |
645 |
|
(\theta^*,S^*) & = & (\theta^{n},S^{n}) + \Delta t G_{\theta,S}^{(n+1/2)} |
646 |
|
\label{eq:tstar-staggered} \\ |
647 |
|
(\theta^{n+1},S^{n+1}) & = & {\cal L}^{-1}_{\theta,S} (\theta^*,S^*) |
648 |
|
\label{eq:t-n+1-staggered} |
649 |
\end{eqnarray} |
\end{eqnarray} |
650 |
The calling sequence is unchanged from |
The corresponding calling tree is given in |
651 |
Fig.~\ref{fig:call-tree-adams-bashforth-sync}. The staggered algorithm |
\ref{fig:call-tree-adams-bashforth-staggered}. |
652 |
is activated with the run-time flag {\bf staggerTimeStep=.TRUE.} in |
The staggered algorithm is activated with the run-time flag |
653 |
{\em PARM01} of {\em data}. |
{\bf staggerTimeStep}{\em=.TRUE.} in parameter file {\em data}, |
654 |
|
namelist {\em PARM01}. |
655 |
|
|
656 |
|
\begin{figure} |
657 |
|
\begin{center} \fbox{ \begin{minipage}{4.7in} \begin{tabbing} |
658 |
|
aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill |
659 |
|
FORWARD\_STEP \\ |
660 |
|
\>\> EXTERNAL\_FIELDS\_LOAD\\ |
661 |
|
\>\> DO\_ATMOSPHERIC\_PHYS \\ |
662 |
|
\>\> DO\_OCEANIC\_PHYS \\ |
663 |
|
\> DYNAMICS \\ |
664 |
|
\>\> CALC\_PHI\_HYD \` $\phi_{hyd}^n$ (\ref{eq:phi-hyd-staggered}) \\ |
665 |
|
\>\> MOM\_FLUXFORM or MOM\_VECINV \` $G_{\vec{\bf v}}^{n-1/2}$ |
666 |
|
(\ref{eq:Gv-n-staggered})\\ |
667 |
|
\>\> TIMESTEP \` $\vec{\bf v}^*$ (\ref{eq:Gv-n+5-staggered}, |
668 |
|
\ref{eq:vstar-staggered}) \\ |
669 |
|
\>\> IMPLDIFF \` $\vec{\bf v}^{**}$ (\ref{eq:vstarstar-staggered}) \\ |
670 |
|
\> UPDATE\_R\_STAR or UPDATE\_SURF\_DR \` (NonLin-FS only)\\ |
671 |
|
\> SOLVE\_FOR\_PRESSURE \\ |
672 |
|
\>\> CALC\_DIV\_GHAT \` $\eta^*$ (\ref{eq:nstar-staggered}) \\ |
673 |
|
\>\> CG2D \` $\eta^{n+1/2}$ (\ref{eq:elliptic-staggered}) \\ |
674 |
|
\> MOMENTUM\_CORRECTION\_STEP \\ |
675 |
|
\>\> CALC\_GRAD\_PHI\_SURF \` $\nabla \eta^{n+1/2}$ \\ |
676 |
|
\>\> CORRECTION\_STEP \` $u^{n+1/2}$,$v^{n+1/2}$ (\ref{eq:v-n+1-staggered})\\ |
677 |
|
\> THERMODYNAMICS \\ |
678 |
|
\>\> CALC\_GT \\ |
679 |
|
\>\>\> GAD\_CALC\_RHS \` $G_\theta^n = G_\theta( u, \theta^n )$ |
680 |
|
(\ref{eq:Gt-n-staggered})\\ |
681 |
|
\>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\ |
682 |
|
\>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:Gt-n+5-staggered}) \\ |
683 |
|
\>\> TIMESTEP\_TRACER \` $\theta^*$ (\ref{eq:tstar-staggered}) \\ |
684 |
|
\>\> IMPLDIFF \` $\theta^{(n+1)}$ (\ref{eq:t-n+1-staggered}) \\ |
685 |
|
\> TRACERS\_CORRECTION\_STEP \\ |
686 |
|
\>\> CYCLE\_TRACER \` $\theta^{n+1}$ \\ |
687 |
|
\>\> FILTER \` Shapiro Filter, Zonal Filter (FFT) \\ |
688 |
|
\>\> CONVECTIVE\_ADJUSTMENT \` \\ |
689 |
|
\end{tabbing} \end{minipage} } \end{center} |
690 |
|
\caption{ |
691 |
|
Calling tree for the overall staggered algorithm using |
692 |
|
Adams-Bashforth time-stepping. |
693 |
|
The place where the model geometry |
694 |
|
({\bf hFac} factors) is updated is added here but is only relevant |
695 |
|
for the non-linear free-surface algorithm. |
696 |
|
} |
697 |
|
\label{fig:call-tree-adams-bashforth-staggered} |
698 |
|
\end{figure} |
699 |
|
|
700 |
The only difficulty with this approach is apparent in equation |
The only difficulty with this approach is apparent in equation |
701 |
\ref{eq:Gt-n-staggered} and illustrated by the dotted arrow |
\ref{eq:Gt-n-staggered} and illustrated by the dotted arrow |
702 |
connecting $u,v^n$ with $G_\theta^{n-1/2}$. The flow used to advect |
connecting $u,v^{n+1/2}$ with $G_\theta^{n}$. The flow used to advect |
703 |
tracers around is not naturally located in time. This could be avoided |
tracers around is not naturally located in time. This could be avoided |
704 |
by applying the Adams-Bashforth extrapolation to the tracer field |
by applying the Adams-Bashforth extrapolation to the tracer field |
705 |
itself and advection that around but this is not yet available. We're |
itself and advecting that around but this approach is not yet |
706 |
not aware of any detrimental effect of this feature. The difficulty |
available. We're not aware of any detrimental effect of this |
707 |
lies mainly in interpretation of what time-level variables and terms |
feature. The difficulty lies mainly in interpretation of what |
708 |
correspond to. |
time-level variables and terms correspond to. |
709 |
|
|
710 |
|
|
711 |
\section{Non-hydrostatic formulation} |
\section{Non-hydrostatic formulation} |
712 |
\label{sect:non-hydrostatic} |
\label{sect:non-hydrostatic} |
713 |
|
\begin{rawhtml} |
714 |
|
<!-- CMIREDIR:non-hydrostatic_formulation: --> |
715 |
|
\end{rawhtml} |
716 |
|
|
717 |
|
The non-hydrostatic formulation re-introduces the full vertical |
718 |
|
momentum equation and requires the solution of a 3-D elliptic |
719 |
|
equations for non-hydrostatic pressure perturbation. We still |
720 |
|
intergrate vertically for the hydrostatic pressure and solve a 2-D |
721 |
|
elliptic equation for the surface pressure/elevation for this reduces |
722 |
|
the amount of work needed to solve for the non-hydrostatic pressure. |
723 |
|
|
724 |
[to be written...] |
The momentum equations are discretized in time as follows: |
725 |
|
\begin{eqnarray} |
726 |
|
\frac{1}{\Delta t} u^{n+1} + g \partial_x \eta^{n+1} + \partial_x \phi_{nh}^{n+1} |
727 |
|
& = & \frac{1}{\Delta t} u^{n} + G_u^{(n+1/2)} \label{eq:discrete-time-u-nh} \\ |
728 |
|
\frac{1}{\Delta t} v^{n+1} + g \partial_y \eta^{n+1} + \partial_y \phi_{nh}^{n+1} |
729 |
|
& = & \frac{1}{\Delta t} v^{n} + G_v^{(n+1/2)} \label{eq:discrete-time-v-nh} \\ |
730 |
|
\frac{1}{\Delta t} w^{n+1} + \partial_r \phi_{nh}^{n+1} |
731 |
|
& = & \frac{1}{\Delta t} w^{n} + G_w^{(n+1/2)} \label{eq:discrete-time-w-nh} |
732 |
|
\end{eqnarray} |
733 |
|
which must satisfy the discrete-in-time depth integrated continuity, |
734 |
|
equation~\ref{eq:discrete-time-backward-free-surface} and the local continuity equation |
735 |
|
\begin{equation} |
736 |
|
\partial_x u^{n+1} + \partial_y v^{n+1} + \partial_r w^{n+1} = 0 |
737 |
|
\label{eq:non-divergence-nh} |
738 |
|
\end{equation} |
739 |
|
As before, the explicit predictions for momentum are consolidated as: |
740 |
|
\begin{eqnarray*} |
741 |
|
u^* & = & u^n + \Delta t G_u^{(n+1/2)} \\ |
742 |
|
v^* & = & v^n + \Delta t G_v^{(n+1/2)} \\ |
743 |
|
w^* & = & w^n + \Delta t G_w^{(n+1/2)} |
744 |
|
\end{eqnarray*} |
745 |
|
but this time we introduce an intermediate step by splitting the |
746 |
|
tendancy of the flow as follows: |
747 |
|
\begin{eqnarray} |
748 |
|
u^{n+1} = u^{**} - \Delta t \partial_x \phi_{nh}^{n+1} |
749 |
|
& & |
750 |
|
u^{**} = u^{*} - \Delta t g \partial_x \eta^{n+1} \\ |
751 |
|
v^{n+1} = v^{**} - \Delta t \partial_y \phi_{nh}^{n+1} |
752 |
|
& & |
753 |
|
v^{**} = v^{*} - \Delta t g \partial_y \eta^{n+1} |
754 |
|
\end{eqnarray} |
755 |
|
Substituting into the depth integrated continuity |
756 |
|
(equation~\ref{eq:discrete-time-backward-free-surface}) gives |
757 |
|
\begin{equation} |
758 |
|
\partial_x H \partial_x \left( g \eta^{n+1} + \widehat{\phi}_{nh}^{n+1} \right) |
759 |
|
+ |
760 |
|
\partial_y H \partial_y \left( g \eta^{n+1} + \widehat{\phi}_{nh}^{n+1} \right) |
761 |
|
- \frac{\epsilon_{fs}\eta^*}{\Delta t^2} |
762 |
|
= - \frac{\eta^*}{\Delta t^2} |
763 |
|
\end{equation} |
764 |
|
which is approximated by equation |
765 |
|
\ref{eq:elliptic-backward-free-surface} on the basis that i) |
766 |
|
$\phi_{nh}^{n+1}$ is not yet known and ii) $\nabla \widehat{\phi}_{nh} |
767 |
|
<< g \nabla \eta$. If \ref{eq:elliptic-backward-free-surface} is |
768 |
|
solved accurately then the implication is that $\widehat{\phi}_{nh} |
769 |
|
\approx 0$ so that thet non-hydrostatic pressure field does not drive |
770 |
|
barotropic motion. |
771 |
|
|
772 |
|
The flow must satisfy non-divergence |
773 |
|
(equation~\ref{eq:non-divergence-nh}) locally, as well as depth |
774 |
|
integrated, and this constraint is used to form a 3-D elliptic |
775 |
|
equations for $\phi_{nh}^{n+1}$: |
776 |
|
\begin{equation} |
777 |
|
\partial_{xx} \phi_{nh}^{n+1} + \partial_{yy} \phi_{nh}^{n+1} + |
778 |
|
\partial_{rr} \phi_{nh}^{n+1} = |
779 |
|
\partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*} |
780 |
|
\end{equation} |
781 |
|
|
782 |
|
The entire algorithm can be summarized as the sequential solution of |
783 |
|
the following equations: |
784 |
|
\begin{eqnarray} |
785 |
|
u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2)} \label{eq:ustar-nh} \\ |
786 |
|
v^{*} & = & v^{n} + \Delta t G_v^{(n+1/2)} \label{eq:vstar-nh} \\ |
787 |
|
w^{*} & = & w^{n} + \Delta t G_w^{(n+1/2)} \label{eq:wstar-nh} \\ |
788 |
|
\eta^* ~ = ~ \epsilon_{fs} \left( \eta^{n} + \Delta t (P-E) \right) |
789 |
|
& - & \Delta t |
790 |
|
\partial_x H \widehat{u^{*}} |
791 |
|
+ \partial_y H \widehat{v^{*}} |
792 |
|
\\ |
793 |
|
\partial_x g H \partial_x \eta^{n+1} |
794 |
|
+ \partial_y g H \partial_y \eta^{n+1} |
795 |
|
& - & \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2} |
796 |
|
~ = ~ |
797 |
|
- \frac{\eta^*}{\Delta t^2} |
798 |
|
\label{eq:elliptic-nh} |
799 |
|
\\ |
800 |
|
u^{**} & = & u^{*} - \Delta t g \partial_x \eta^{n+1} \label{eq:unx-nh}\\ |
801 |
|
v^{**} & = & v^{*} - \Delta t g \partial_y \eta^{n+1} \label{eq:vnx-nh}\\ |
802 |
|
\partial_{xx} \phi_{nh}^{n+1} + \partial_{yy} \phi_{nh}^{n+1} + |
803 |
|
\partial_{rr} \phi_{nh}^{n+1} & = & |
804 |
|
\partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*} \\ |
805 |
|
u^{n+1} & = & u^{**} - \Delta t \partial_x \phi_{nh}^{n+1} \label{eq:un+1-nh}\\ |
806 |
|
v^{n+1} & = & v^{**} - \Delta t \partial_y \phi_{nh}^{n+1} \label{eq:vn+1-nh}\\ |
807 |
|
\partial_r w^{n+1} & = & - \partial_x u^{n+1} - \partial_y v^{n+1} |
808 |
|
\end{eqnarray} |
809 |
|
where the last equation is solved by vertically integrating for |
810 |
|
$w^{n+1}$. |
811 |
|
|
812 |
|
|
813 |
|
|
814 |
\section{Variants on the Free Surface} |
\section{Variants on the Free Surface} |
815 |
|
\label{sect:free-surface} |
816 |
|
|
817 |
We now describe the various formulations of the free-surface that |
We now describe the various formulations of the free-surface that |
818 |
include non-linear forms, implicit in time using Crank-Nicholson, |
include non-linear forms, implicit in time using Crank-Nicholson, |
832 |
\begin{eqnarray} |
\begin{eqnarray} |
833 |
{\eta}^* = \epsilon_{fs} \: {\eta}^{n} - |
{\eta}^* = \epsilon_{fs} \: {\eta}^{n} - |
834 |
\Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o} \vec{\bf v}^* dr |
\Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o} \vec{\bf v}^* dr |
835 |
\: + \: \epsilon_{fw} \Delta_t (P-E)^{n} |
\: + \: \epsilon_{fw} \Delta t (P-E)^{n} |
836 |
\label{eq-solve2D_rhs} |
\label{eq-solve2D_rhs} |
837 |
\end{eqnarray} |
\end{eqnarray} |
838 |
|
|
839 |
\fbox{ \begin{minipage}{4.75in} |
\fbox{ \begin{minipage}{4.75in} |
840 |
{\em S/R SOLVE\_FOR\_PRESSURE} ({\em solve\_for\_pressure.F}) |
{\em S/R SOLVE\_FOR\_PRESSURE} ({\em solve\_for\_pressure.F}) |
841 |
|
|
842 |
$u^*$: {\bf GuNm1} ({\em DYNVARS.h}) |
$u^*$: {\bf gU} ({\em DYNVARS.h}) |
843 |
|
|
844 |
$v^*$: {\bf GvNm1} ({\em DYNVARS.h}) |
$v^*$: {\bf gV} ({\em DYNVARS.h}) |
845 |
|
|
846 |
$\eta^*$: {\bf cg2d\_b} (\em SOLVE\_FOR\_PRESSURE.h) |
$\eta^*$: {\bf cg2d\_b} (\em SOLVE\_FOR\_PRESSURE.h) |
847 |
|
|
851 |
|
|
852 |
|
|
853 |
Once ${\eta}^{n+1}$ has been found, substituting into |
Once ${\eta}^{n+1}$ has been found, substituting into |
854 |
\ref{eq:discrete-time-u,eq:discrete-time-v} yields $\vec{\bf v}^{n+1}$ if the model is |
\ref{eq:discrete-time-u}, \ref{eq:discrete-time-v} yields $\vec{\bf v}^{n+1}$ |
855 |
hydrostatic ($\epsilon_{nh}=0$): |
if the model is hydrostatic ($\epsilon_{nh}=0$): |
856 |
$$ |
$$ |
857 |
\vec{\bf v}^{n+1} = \vec{\bf v}^{*} |
\vec{\bf v}^{n+1} = \vec{\bf v}^{*} |
858 |
- \Delta t {\bf \nabla}_h b_s {\eta}^{n+1} |
- \Delta t {\bf \nabla}_h b_s {\eta}^{n+1} |
861 |
This is known as the correction step. However, when the model is |
This is known as the correction step. However, when the model is |
862 |
non-hydrostatic ($\epsilon_{nh}=1$) we need an additional step and an |
non-hydrostatic ($\epsilon_{nh}=1$) we need an additional step and an |
863 |
additional equation for $\phi'_{nh}$. This is obtained by substituting |
additional equation for $\phi'_{nh}$. This is obtained by substituting |
864 |
\ref{eq:discrete-time-u}, \ref{eq:discrete-time-v} and \ref{eq:discrete-time-w} |
\ref{eq:discrete-time-u-nh}, \ref{eq:discrete-time-v-nh} and \ref{eq:discrete-time-w-nh} |
865 |
into continuity: |
into continuity: |
866 |
\begin{equation} |
\begin{equation} |
867 |
\left[ {\bf \nabla}_h^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1} |
\left[ {\bf \nabla}_h^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1} |
893 |
|
|
894 |
$\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h) |
$\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h) |
895 |
|
|
896 |
$\phi_{nh}^{n+1}$: {\bf phi\_nh} (\em DYNVARS.h) |
$\phi_{nh}^{n+1}$: {\bf phi\_nh} (\em NH\_VARS.h) |
897 |
|
|
898 |
$u^*$: {\bf GuNm1} ({\em DYNVARS.h}) |
$u^*$: {\bf gU} ({\em DYNVARS.h}) |
899 |
|
|
900 |
$v^*$: {\bf GvNm1} ({\em DYNVARS.h}) |
$v^*$: {\bf gV} ({\em DYNVARS.h}) |
901 |
|
|
902 |
$u^{n+1}$: {\bf uVel} ({\em DYNVARS.h}) |
$u^{n+1}$: {\bf uVel} ({\em DYNVARS.h}) |
903 |
|
|
925 |
|
|
926 |
|
|
927 |
\subsection{Crank-Nickelson barotropic time stepping} |
\subsection{Crank-Nickelson barotropic time stepping} |
928 |
|
\label{sect:freesurf-CrankNick} |
929 |
|
|
930 |
The full implicit time stepping described previously is |
The full implicit time stepping described previously is |
931 |
unconditionally stable but damps the fast gravity waves, resulting in |
unconditionally stable but damps the fast gravity waves, resulting in |
940 |
corresponds to the forward - backward scheme that conserves energy but is |
corresponds to the forward - backward scheme that conserves energy but is |
941 |
only stable for small time steps.\\ |
only stable for small time steps.\\ |
942 |
In the code, $\beta,\gamma$ are defined as parameters, respectively |
In the code, $\beta,\gamma$ are defined as parameters, respectively |
943 |
{\it implicSurfPress}, {\it implicDiv2DFlow}. They are read from |
{\bf implicSurfPress}, {\bf implicDiv2DFlow}. They are read from |
944 |
the main data file "{\it data}" and are set by default to 1,1. |
the main parameter file "{\em data}" and are set by default to 1,1. |
945 |
|
|
946 |
Equations \ref{eq:ustar-backward-free-surface} -- |
Equations \ref{eq:ustar-backward-free-surface} -- |
947 |
\ref{eq:vn+1-backward-free-surface} are modified as follows: |
\ref{eq:vn+1-backward-free-surface} are modified as follows: |
948 |
$$ |
\begin{eqnarray*} |
949 |
\frac{ \vec{\bf v}^{n+1} }{ \Delta t } |
\frac{ \vec{\bf v}^{n+1} }{ \Delta t } |
950 |
+ {\bf \nabla}_h b_s [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ] |
+ {\bf \nabla}_h b_s [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ] |
951 |
+ \epsilon_{nh} {\bf \nabla}_h {\phi'_{nh}}^{n+1} |
+ \epsilon_{nh} {\bf \nabla}_h {\phi'_{nh}}^{n+1} |
952 |
= \frac{ \vec{\bf v}^* }{ \Delta t } |
= \frac{ \vec{\bf v}^* }{ \Delta t } |
953 |
$$ |
\end{eqnarray*} |
954 |
$$ |
\begin{eqnarray} |
955 |
\epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t} |
\epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t} |
956 |
+ {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o} |
+ {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o} |
957 |
[ \gamma \vec{\bf v}^{n+1} + (1-\gamma) \vec{\bf v}^{n}] dr |
[ \gamma \vec{\bf v}^{n+1} + (1-\gamma) \vec{\bf v}^{n}] dr |
958 |
= \epsilon_{fw} (P-E) |
= \epsilon_{fw} (P-E) |
959 |
$$ |
\label{eq:eta-n+1-CrankNick} |
960 |
|
\end{eqnarray} |
961 |
where: |
where: |
962 |
\begin{eqnarray*} |
\begin{eqnarray*} |
963 |
\vec{\bf v}^* & = & |
\vec{\bf v}^* & = & |
979 |
{\bf \nabla}_h {\eta}^{n+1} |
{\bf \nabla}_h {\eta}^{n+1} |
980 |
= {\eta}^* |
= {\eta}^* |
981 |
$$ |
$$ |
982 |
and then to compute (correction step): |
and then to compute ({\em CORRECTION\_STEP}): |
983 |
$$ |
$$ |
984 |
\vec{\bf v}^{n+1} = \vec{\bf v}^{*} |
\vec{\bf v}^{n+1} = \vec{\bf v}^{*} |
985 |
- \beta \Delta t {\bf \nabla}_h b_s {\eta}^{n+1} |
- \beta \Delta t {\bf \nabla}_h b_s {\eta}^{n+1} |
986 |
$$ |
$$ |
987 |
|
|
988 |
The non-hydrostatic part is solved as described previously. |
%The non-hydrostatic part is solved as described previously. |
989 |
|
|
990 |
Note that: |
\noindent |
991 |
|
Notes: |
992 |
\begin{enumerate} |
\begin{enumerate} |
993 |
|
\item The RHS term of equation \ref{eq:eta-n+1-CrankNick} |
994 |
|
corresponds the contribution of fresh water flux (P-E) |
995 |
|
to the free-surface variations ($\epsilon_{fw}=1$, |
996 |
|
{\bf useRealFreshWater}{\em=TRUE} in parameter file {\em data}). |
997 |
|
In order to remain consistent with the tracer equation, specially in |
998 |
|
the non-linear free-surface formulation, this term is also |
999 |
|
affected by the Crank-Nickelson time stepping. The RHS reads: |
1000 |
|
$\epsilon_{fw} ( \gamma (P-E)^{n+1/2} + (1-\gamma) (P-E)^{n-1/2} )$ |
1001 |
\item The non-hydrostatic part of the code has not yet been |
\item The non-hydrostatic part of the code has not yet been |
1002 |
updated, so that this option cannot be used with $(\beta,\gamma) \neq (1,1)$. |
updated, and therefore cannot be used with $(\beta,\gamma) \neq (1,1)$. |
1003 |
\item The stability criteria with Crank-Nickelson time stepping |
\item The stability criteria with Crank-Nickelson time stepping |
1004 |
for the pure linear gravity wave problem in cartesian coordinates is: |
for the pure linear gravity wave problem in cartesian coordinates is: |
1005 |
\begin{itemize} |
\begin{itemize} |