/[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.17 by jmc, Wed Oct 13 18:50:54 2004 UTC revision 1.26 by jmc, Thu Jun 29 01:45:32 2006 UTC
# Line 10  describe the spatial discretization. The Line 10  describe the spatial discretization. The
10  terms are described first, afterwards the schemes that apply to  terms are described first, afterwards the schemes that apply to
11  passive and dynamically active tracers are described.  passive and dynamically active tracers are described.
12    
13    \input{part2/notation}
14    
15  \section{Time-stepping}  \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 61  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 189  The calling tree for these routines is g Line 198  The calling tree for these routines is g
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 214  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 222  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 244  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 267  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 279  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 302  FORWARD\_STEP \\ Line 324  FORWARD\_STEP \\
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 333  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 forcing term outside the  The model offers the possibility to leave the tracer and momentum
361  Adams-Bashforth extrapolation, by turning off the logical flag  forcing terms and the dissipation terms outside the
362  {\bf forcing\_In\_AB } (parameter file {\em data}, namelist {\em PARM01},  Adams-Bashforth extrapolation, by turning off the logical flags
363  default value = True).  \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 344  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.  using the backward method which is an intrinsic scheme.
# Line 391  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 422  FORWARD\_STEP \\ Line 479  FORWARD\_STEP \\
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})\\
# Line 445  FORWARD\_STEP \\ Line 502  FORWARD\_STEP \\
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  The place where the model geometry
505  ({\em hFac} factors) is updated is added here but is only relevant  ({\bf hFac} factors) is updated is added here but is only relevant
506  for the non-linear free-surface algorithm.  for the non-linear free-surface algorithm.
507  For completeness, the external forcing,  For completeness, the external forcing,
508  ocean and atmospheric physics have been added, although they are mainly  ocean and atmospheric physics have been added, although they are mainly
# Line 508  time-step. The corresponding calling tre Line 565  time-step. The corresponding calling tre
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 516  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$ 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+1/2$ (dotted arrow). The  function of the flow field at that time level $n-1/2$.
582  explicit tendency from the previous time level, $n-1$, is used to  The explicit tendency from the previous time level, $n-3/2$, is used to
583  extrapolate tendencies to $n+1/2$ (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$  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$. 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 at time level  These are then used to calculate the advection term (dashed arc-arrow)
591  $n+1$ in stepping forward the flow variables from $n+1/2$ to $n+3/2$  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 549  thermodynamics solver is delayed from ha Line 612  thermodynamics solver is delayed from ha
612  allowing the use of the most recent velocities to compute  allowing the use of the most recent velocities to compute
613  the advection terms. Once the thermodynamics fields are  the advection terms. Once the thermodynamics fields are
614  updated, the hydrostatic pressure is computed  updated, the hydrostatic pressure is computed
615  to step frowrad the dynamics  to step forwrad the dynamics.
616  Note that 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
# Line 557  located in time.  Instead, we re-write t Line 620  located in time.  Instead, we re-write t
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    \phi^{n}_{hyd} & = & \int b(\theta^{n},S^{n}) dr
624    \label{eq:phi-hyd-staggered} \\
625  \vec{\bf G}_{\vec{\bf v}}^{n-1/2} & = & \vec{\bf G}_{\vec{\bf v}} ( \vec{\bf v}^{n-1/2} )  \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)} & = & (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}  \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} \\
 \phi^{n}_{hyd} & = & \int b(\theta^{n},S^{n}) dr  
 \label{eq:phi-hyd-staggered} \\  
629  \vec{\bf v}^{*} & = & \vec{\bf v}^{n-1/2} + \Delta t \left( \vec{\bf G}_{\vec{\bf v}}^{(n)} - \nabla \phi_{hyd}^{n} \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}^* )
# Line 582  G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil Line 645  G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil
645  (\theta^*,S^*) & = & (\theta^{n},S^{n}) + \Delta t G_{\theta,S}^{(n+1/2)}  (\theta^*,S^*) & = & (\theta^{n},S^{n}) + \Delta t G_{\theta,S}^{(n+1/2)}
646  \label{eq:tstar-staggered} \\  \label{eq:tstar-staggered} \\
647  (\theta^{n+1},S^{n+1}) & = & {\cal L}^{-1}_{\theta,S} (\theta^*,S^*)  (\theta^{n+1},S^{n+1}) & = & {\cal L}^{-1}_{\theta,S} (\theta^*,S^*)
648  \label{eq:t-n+1-staggered} \\  \label{eq:t-n+1-staggered}
649  \end{eqnarray}  \end{eqnarray}
650  The corresponding calling tree is given in  The corresponding calling tree is given in
651  \ref{fig:call-tree-adams-bashforth-staggered}.  \ref{fig:call-tree-adams-bashforth-staggered}.
652  The staggered algorithm is activated with the run-time flag  The staggered algorithm is activated with the run-time flag
653  {\bf staggerTimeStep=.TRUE.} in parameter file {\em data},  {\bf staggerTimeStep}{\em=.TRUE.} in parameter file {\em data},
654  namelist {\em PARM01}.  namelist {\em PARM01}.
655    
656  \begin{figure}  \begin{figure}
# Line 617  FORWARD\_STEP \\ Line 680  FORWARD\_STEP \\
680       (\ref{eq:Gt-n-staggered})\\       (\ref{eq:Gt-n-staggered})\\
681  \>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\  \>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\
682  \>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:Gt-n+5-staggered}) \\  \>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:Gt-n+5-staggered}) \\
683  \>\> TIMESTEP\_TRACER \` $\tau^*$ (\ref{eq:tstar-staggered}) \\  \>\> TIMESTEP\_TRACER \` $\theta^*$ (\ref{eq:tstar-staggered}) \\
684  \>\> IMPLDIFF \` $\tau^{(n+1)}$ (\ref{eq:t-n+1-staggered}) \\  \>\> IMPLDIFF \` $\theta^{(n+1)}$ (\ref{eq:t-n+1-staggered}) \\
685  \> TRACERS\_CORRECTION\_STEP  \\  \> TRACERS\_CORRECTION\_STEP  \\
686  \>\> CYCLE\_TRACER \` $\theta^{n+1}$ \\  \>\> CYCLE\_TRACER \` $\theta^{n+1}$ \\
687  \>\> FILTER \` Shapiro Filter, Zonal Filter (FFT) \\  \>\> FILTER \` Shapiro Filter, Zonal Filter (FFT) \\
# Line 628  FORWARD\_STEP \\ Line 691  FORWARD\_STEP \\
691  Calling tree for the overall staggered algorithm using  Calling tree for the overall staggered algorithm using
692  Adams-Bashforth time-stepping.  Adams-Bashforth time-stepping.
693  The place where the model geometry  The place where the model geometry
694  ({\em hFac} factors) is updated is added here but is only relevant  ({\bf hFac} factors) is updated is added here but is only relevant
695  for the non-linear free-surface algorithm.  for the non-linear free-surface algorithm.
696  }  }
697  \label{fig:call-tree-adams-bashforth-staggered}  \label{fig:call-tree-adams-bashforth-staggered}
# Line 636  for the non-linear free-surface algorith Line 699  for the non-linear free-surface algorith
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 advecting that around but this approach is not yet  itself and advecting that around but this approach is not yet
# Line 647  time-level variables and terms correspon Line 710  time-level variables and terms correspon
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  The non-hydrostatic formulation re-introduces the full vertical
718  momentum equation and requires the solution of a 3-D elliptic  momentum equation and requires the solution of a 3-D elliptic
# Line 662  The momentum equations are discretized i Line 728  The momentum equations are discretized i
728  \frac{1}{\Delta t} v^{n+1} + g \partial_y \eta^{n+1} + \partial_y \phi_{nh}^{n+1}  \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} \\  & = & \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}  \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} \\  & = & \frac{1}{\Delta t} w^{n} + G_w^{(n+1/2)} \label{eq:discrete-time-w-nh}
732  \end{eqnarray}  \end{eqnarray}
733  which must satisfy the discrete-in-time depth integrated continuity,  which must satisfy the discrete-in-time depth integrated continuity,
734  equation~\ref{eq:discrete-time-backward-free-surface} and the local continuity equation  equation~\ref{eq:discrete-time-backward-free-surface} and the local continuity equation
# Line 746  $w^{n+1}$. Line 812  $w^{n+1}$.
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 765  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 784  $\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 826  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 858  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 872  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 910  $$ 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.17  
changed lines
  Added in v.1.26

  ViewVC Help
Powered by ViewVC 1.1.22