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

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

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

revision 1.12 by adcroft, Tue Nov 13 20:27:54 2001 UTC revision 1.26 by jmc, Thu Jun 29 01:45:32 2006 UTC
# Line 1  Line 1 
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,
# Line 50  treated more exactly, including non-line Line 66  treated more exactly, including non-line
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}
# Line 70  temporary.} Line 90  temporary.}
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    
# Line 162  The correspondence to the code is as fol Line 183  The correspondence to the code is as fol
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
# Line 201  The rigid-lid approximation can be easil Line 232  The rigid-lid approximation can be easil
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
# Line 209  and fresh-water source term. Line 240  and fresh-water source 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
# Line 231  As for the rigid-lid pressure method, eq Line 262  As for the rigid-lid pressure method, eq
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  \\  \\
# Line 254  recovered. However, the implicit treatme Line 285  recovered. However, the implicit treatme
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
# Line 266  pressure-method. Line 298  pressure-method.
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
# Line 281  FORWARD\_STEP \\ Line 316  FORWARD\_STEP \\
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    
# Line 319  simpler to include these terms and this Line 357  simpler to include these terms and this
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...}
# Line 326  A stability analysis for an oscillation Line 374  A stability analysis for an oscillation
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} =
# Line 361  Fig.~\ref{fig:call-tree-adams-bashforth} Line 437  Fig.~\ref{fig:call-tree-adams-bashforth}
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
# Line 369  implicit and are thus cast as a an expli Line 445  implicit and are thus cast as a an expli
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}
# Line 389  is solved to yield the state variables a Line 468  is solved to yield the state variables a
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    
# Line 442  G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil Line 535  G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil
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}
# Line 465  accelerations, stepping forward and solv Line 558  accelerations, stepping forward and solv
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}
# Line 480  time-step. The corresponding calling tre Line 576  time-step. The corresponding calling tre
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    
# Line 503  circumstance, it is more efficient to st Line 602  circumstance, it is more efficient to st
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,
# Line 592  where Line 832  where
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    
# Line 611  $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h) Line 851  $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)
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}
# Line 621  $$ Line 861  $$
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}
# Line 653  without any consequence on the solution. Line 893  without any consequence on the solution.
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    
# Line 685  at the same point in the code. Line 925  at the same point in the code.
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
# Line 699  stable, Crank-Nickelson scheme; $(\beta, Line 940  stable, Crank-Nickelson scheme; $(\beta,
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}^* & = &
# Line 737  $$ Line 979  $$
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}

Legend:
Removed from v.1.12  
changed lines
  Added in v.1.26

  ViewVC Help
Powered by ViewVC 1.1.22