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

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

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

-revision 1.8 by adcroft,
Fri Oct  5 19:41:08 2001 UTC
+revision 1.25 by jmc,
Wed Jun 28 16:55:53 2006 UTC
 Line 1
  % $Header$
  % $Name$
+ This chapter lays out the numerical schemes that are
+ employed in the core MITgcm algorithm. Whenever possible
+ links are made to actual program code in the MITgcm implementation.
+ The chapter begins with a discussion of the temporal discretization
+ used in MITgcm. This discussion is followed by sections that
+ describe the spatial discretization. The schemes employed for momentum
+ terms are described first, afterwards the schemes that apply to
+ passive and dynamically active tracers are described.
+ \input{part2/notation}
+ \section{Time-stepping}
+ \begin{rawhtml}
+ <!-- CMIREDIR:time-stepping: -->
+ \end{rawhtml}
  The equations of motion integrated by the model involve four
  prognostic equations for flow, $u$ and $v$, temperature, $\theta$, and
  salt/moisture, $S$, and three diagnostic equations for vertical flow,
-Line 13 
 and diagnostic equations requires a mode
+Line 29 
 and diagnostic equations requires a mode
  forward prognostic variables while satisfying constraints imposed by
  diagnostic equations.
- Since the model comes in several flavours and formulation, it would be
+ Since the model comes in several flavors and formulation, it would be
  confusing to present the model algorithm exactly as written into code
  along with all the switches and optional terms. Instead, we present
  the algorithm for each of the basic formulations which are:
-Line 37 
 method, with a rigid-lid model in sectio
+Line 53 
 method, with a rigid-lid model in sectio
  \ref{sect:pressure-method-rigid-lid}. This algorithm is essentially
  unchanged, apart for some coefficients, when the rigid lid assumption
  is replaced with a linearized implicit free-surface, described in
- section \ref{sect:pressure-method-linear-backward}. These two flavours
+ section \ref{sect:pressure-method-linear-backward}. These two flavors
  of the pressure-method encompass all formulations of the model as it
  exists today. The integration of explicit in time terms is out-lined
  in section \ref{sect:adams-bashforth} and put into the context of the
-Line 50 
 treated more exactly, including non-line
+Line 66 
 treated more exactly, including non-line
  described in section \ref{sect:nonlinear-freesurface}.
- \section{Pressure method with rigid-lid} \label{sect:pressure-method-rigid-lid}
+ \section{Pressure method with rigid-lid}
+ \label{sect:pressure-method-rigid-lid}
+ \begin{rawhtml}
+ <!-- CMIREDIR:pressure_method_rigid_lid: -->
+ \end{rawhtml}
  \begin{figure}
  \begin{center}
-Line 61 
 A schematic of the evolution in time of
+Line 81 
 A schematic of the evolution in time of
  algorithm. A prediction for the flow variables at time level $n+1$ is
  made based only on the explicit terms, $G^{(n+^1/_2)}$, and denoted
  $u^*$, $v^*$. Next, a pressure field is found such that $u^{n+1}$,
- $v^{n+1}$ will be non-divergent. Conceptually, the $*$ quantitites
+ $v^{n+1}$ will be non-divergent. Conceptually, the $*$ quantities
  exist at time level $n+1$ but they are intermediate and only
  temporary.}
  \label{fig:pressure-method-rigid-lid}
-Line 70 
 temporary.}
+Line 90 
 temporary.}
  \begin{figure}
  \begin{center} \fbox{ \begin{minipage}{4.5in} \begin{tabbing}
  aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill
- FORWARD\_STEP \\
+ \filelink{FORWARD\_STEP}{model-src-forward_step.F} \\
  \> DYNAMICS \\
  \>\> TIMESTEP \` $u^*$,$v^*$ (\ref{eq:ustar-rigid-lid},\ref{eq:vstar-rigid-lid}) \\
  \> SOLVE\_FOR\_PRESSURE \\
  \>\> CALC\_DIV\_GHAT \` $H\widehat{u^*}$,$H\widehat{v^*}$ (\ref{eq:elliptic}) \\
  \>\> CG2D \` $\eta^{n+1}$ (\ref{eq:elliptic}) \\
- \> THE\_CORRECTION\_STEP  \\
+ \> MOMENTUM\_CORRECTION\_STEP  \\
  \>\> CALC\_GRAD\_PHI\_SURF \` $\nabla \eta^{n+1}$ \\
  \>\> CORRECTION\_STEP \` $u^{n+1}$,$v^{n+1}$ (\ref{eq:un+1-rigid-lid},\ref{eq:vn+1-rigid-lid})
  \end{tabbing} \end{minipage} } \end{center}
- \caption{Calling tree for the pressure method alogtihm}
+ \caption{Calling tree for the pressure method algorithm
+   (\filelink{FORWARD\_STEP}{model-src-forward_step.F})}
  \label{fig:call-tree-pressure-method}
  \end{figure}
-Line 102 
 flow boundary conditions applied, become
+Line 123 
 flow boundary conditions applied, become
  \label{eq:rigid-lid-continuity}
  \end{equation}
  Here, $H\widehat{u} = \int_H u dz$ is the depth integral of $u$,
- similarly for $H\widehat{v}$. The rigid-lid approzimation sets $w=0$
+ similarly for $H\widehat{v}$. The rigid-lid approximation sets $w=0$
  at the lid so that it does not move but allows a pressure to be
  exerted on the fluid by the lid. The horizontal momentum equations and
  vertically integrated continuity equation are be discretized in time
-Line 130 
 pressure gradient is explicit and so can
+Line 151 
 pressure gradient is explicit and so can
  $G$.
  Substituting the two momentum equations into the depth integrated
- continuity equation eliminates $u^{n+1}$ and $v^{n+1}$ yeilding an
+ continuity equation eliminates $u^{n+1}$ and $v^{n+1}$ yielding an
  elliptic equation for $\eta^{n+1}$. Equations
  \ref{eq:discrete-time-u}, \ref{eq:discrete-time-v} and
  \ref{eq:discrete-time-cont-rigid-lid} can then be re-arranged as follows:
-Line 157 
 Fig.~\ref{fig:pressure-method-rigid-lid}
+Line 178 
 Fig.~\ref{fig:pressure-method-rigid-lid}
  with the future flow field (time level $n+1$) but it could equally
  have been drawn as staggered in time with the flow.
- The correspondance to the code is as follows:
+ The correspondence to the code is as follows:
  \begin{itemize}
  \item
  the prognostic phase, equations \ref{eq:ustar-rigid-lid} and \ref{eq:vstar-rigid-lid},
  stepping forward $u^n$ and $v^n$ to $u^{*}$ and $v^{*}$ is coded in
- {\em TIMESTEP.F}
+ \filelink{TIMESTEP()}{model-src-timestep.F}
  \item
  the vertical integration, $H \widehat{u^*}$ and $H
- \widehat{v^*}$, divergence and invertion of the elliptic operator in
+ \widehat{v^*}$, divergence and inversion of the elliptic operator in
- equation \ref{eq:elliptic} is coded in {\em
+ equation \ref{eq:elliptic} is coded in
- SOLVE\_FOR\_PRESSURE.F}
+ \filelink{SOLVE\_FOR\_PRESSURE()}{model-src-solve_for_pressure.F}
  \item
  finally, the new flow field at time level $n+1$ given by equations
- \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
+ \filelink{CORRECTION\_STEP()}{model-src-correction_step.F}.
  \end{itemize}
  The calling tree for these routines is given in
  Fig.~\ref{fig:call-tree-pressure-method}.
-Line 190 
 G_u^{(n+1/2)}
+Line 212 
 G_u^{(n+1/2)}
  \section{Pressure method with implicit linear free-surface}
  \label{sect:pressure-method-linear-backward}
+ \begin{rawhtml}
+ <!-- CMIREDIR:pressure_method_linear_backward: -->
+ \end{rawhtml}
  The rigid-lid approximation filters out external gravity waves
  subsequently modifying the dispersion relation of barotropic Rossby
-Line 201 
 The rigid-lid approximation can be easil
+Line 226 
 The rigid-lid approximation can be easil
  of the free-surface equation which can be written:
  \begin{equation}
  \partial_t \eta + \partial_x H \widehat{u} + \partial_y H \widehat{v} = P-E+R
- \label{eq:linear-free-surface=P-E+R}
+ \label{eq:linear-free-surface=P-E}
  \end{equation}
  which differs from the depth integrated continuity equation with
  rigid-lid (\ref{eq:rigid-lid-continuity}) by the time-dependent term
-Line 209 
 and fresh-water source term.
+Line 234 
 and fresh-water source term.
  Equation \ref{eq:discrete-time-cont-rigid-lid} in the rigid-lid
  pressure method is then replaced by the time discretization of
- \ref{eq:linear-free-surface=P-E+R} which is:
+ \ref{eq:linear-free-surface=P-E} which is:
  \begin{equation}
  \eta^{n+1}
  + \Delta t \partial_x H \widehat{u^{n+1}}
  + \Delta t \partial_y H \widehat{v^{n+1}}
  =
  \eta^{n}
- + \Delta t ( P - E + R )
+ + \Delta t ( P - E )
  \label{eq:discrete-time-backward-free-surface}
  \end{equation}
  where the use of flow at time level $n+1$ makes the method implicit
-Line 231 
 As for the rigid-lid pressure method, eq
+Line 256 
 As for the rigid-lid pressure method, eq
  \begin{eqnarray}
  u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2)} \label{eq:ustar-backward-free-surface} \\
  v^{*} & = & v^{n} + \Delta t G_v^{(n+1/2)} \label{eq:vstar-backward-free-surface} \\
- \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
    \partial_x H \widehat{u^{*}}
  + \partial_y H \widehat{v^{*}}
  \\
    \partial_x g H \partial_x \eta^{n+1}
- + \partial_y g H \partial_y \eta^{n+1}
+ & + & \partial_y g H \partial_y \eta^{n+1}
- - \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2}
+  - \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2}
- & = &
+  =
  - \frac{\eta^*}{\Delta t^2}
  \label{eq:elliptic-backward-free-surface}
  \\
-Line 250 
 to~\ref{eq:vn+1-backward-free-surface},
+Line 275 
 to~\ref{eq:vn+1-backward-free-surface},
  the pressure method algorithm with a backward implicit, linearized
  free surface. The method is still formerly a pressure method because
  in the limit of large $\Delta t$ the rigid-lid method is
- reovered. However, the implicit treatment of the free-surface allows
+ recovered. However, the implicit treatment of the free-surface allows
  the flow to be divergent and for the surface pressure/elevation to
- respond on a finite time-scale (as opposed to instantly). To recovere
+ respond on a finite time-scale (as opposed to instantly). To recover
  the rigid-lid formulation, we introduced a switch-like parameter,
  $\epsilon_{fs}$, which selects between the free-surface and rigid-lid;
  $\epsilon_{fs}=1$ allows the free-surface to evolve; $\epsilon_{fs}=0$
-Line 266 
 pressure-method.
+Line 291 
 pressure-method.
  \section{Explicit time-stepping: Adams-Bashforth}
  \label{sect:adams-bashforth}
+ \begin{rawhtml}
+ <!-- CMIREDIR:adams_bashforth: -->
+ \end{rawhtml}
  In describing the the pressure method above we deferred describing the
  time discretization of the explicit terms. We have historically used
-Line 281 
 FORWARD\_STEP \\
+Line 309 
 FORWARD\_STEP \\
  \> THERMODYNAMICS \\
  \>\> CALC\_GT \\
  \>\>\> GAD\_CALC\_RHS \` $G_\theta^n = G_\theta( u, \theta^n )$ \\
- \>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\
+ \>either\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\
  \>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:adams-bashforth2}) \\
+ \>or\>\> EXTERNAL\_FORCING \` $G_\theta^{(n+1/2)} = G_\theta^{(n+1/2)} + {\cal Q}$ \\
  \>\> TIMESTEP\_TRACER \` $\tau^*$ (\ref{eq:taustar}) \\
  \>\> IMPLDIFF \` $\tau^{(n+1)}$ (\ref{eq:tau-n+1-implicit})
  \end{tabbing} \end{minipage} } \end{center}
-Line 319 
 simpler to include these terms and this
+Line 348 
 simpler to include these terms and this
  and forcing evolves smoothly. Problems can, and do, arise when forcing
  or motions are high frequency and this corresponds to a reduced
  stability compared to a simple forward time-stepping of such terms.
+ The model offers the possibility to leave the forcing term outside the
+ Adams-Bashforth extrapolation, by turning off the logical flag
+ {\bf forcing\_In\_AB } (parameter file {\em data}, namelist {\em PARM01},
+ default value = True).
  A stability analysis for an oscillation equation should be given at this point.
  \marginpar{AJA needs to find his notes on this...}
-Line 328 
 A stability analysis for a relaxation eq
+Line 361 
 A stability analysis for a relaxation eq
  \section{Implicit time-stepping: backward method}
+ \begin{rawhtml}
+ <!-- CMIREDIR:implicit_time-stepping_backward: -->
+ \end{rawhtml}
  Vertical diffusion and viscosity can be treated implicitly in time
- using the backward method which is an intrinsic scheme. For tracers,
+ using the backward method which is an intrinsic scheme.
- the time discrretized equation is:
+ Recently, the option to treat the vertical advection
+ implicitly has been added, but not yet tested; therefore,
+ the description hereafter is limited to diffusion and viscosity.
+ For tracers,
+ the time discretized equation is:
  \begin{equation}
  \tau^{n+1} - \Delta t \partial_r \kappa_v \partial_r \tau^{n+1} =
  \tau^{n} + \Delta t G_\tau^{(n+1/2)}
-Line 361 
 Fig.~\ref{fig:call-tree-adams-bashforth}
+Line 401 
 Fig.~\ref{fig:call-tree-adams-bashforth}
  stepping forward a tracer variable such as temperature.
  In order to fit within the pressure method, the implicit viscosity
- 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
  redistribute momentum in the vertical. The upshot of this is that
  although vertical viscosity may be backward implicit and
  unconditionally stable, no-slip boundary conditions may not be made
-Line 369 
 implicit and are thus cast as a an expli
+Line 409 
 implicit and are thus cast as a an expli
  \section{Synchronous time-stepping: variables co-located in time}
  \label{sect:adams-bashforth-sync}
+ \begin{rawhtml}
+ <!-- CMIREDIR:adams_bashforth_sync: -->
+ \end{rawhtml}
  \begin{figure}
  \begin{center}
-Line 377 
 implicit and are thus cast as a an expli
+Line 420 
 implicit and are thus cast as a an expli
  \caption{
  A schematic of the explicit Adams-Bashforth and implicit time-stepping
  phases of the algorithm. All prognostic variables are co-located in
- time. Explicit tendancies are evaluated at time level $n$ as a
+ time. Explicit tendencies are evaluated at time level $n$ as a
  function of the state at that time level (dotted arrow). The explicit
- tendancy from the previous time level, $n-1$, is used to extrapolate
+ tendency from the previous time level, $n-1$, is used to extrapolate
- tendancies to $n+1/2$ (dashed arrow). This extrapolated tendancy
+ tendencies to $n+1/2$ (dashed arrow). This extrapolated tendency
  allows variables to be stably integrated forward-in-time to render an
  estimate ($*$-variables) at the $n+1$ time level (solid
  arc-arrow). The operator ${\cal L}$ formed from implicit-in-time terms
-Line 389 
 is solved to yield the state variables a
+Line 432 
 is solved to yield the state variables a
  \end{figure}
  \begin{figure}
- \begin{center} \fbox{ \begin{minipage}{4.5in} \begin{tabbing}
+ \begin{center} \fbox{ \begin{minipage}{4.7in} \begin{tabbing}
  aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill
  FORWARD\_STEP \\
+ \>\> EXTERNAL\_FIELDS\_LOAD\\
+ \>\> DO\_ATMOSPHERIC\_PHYS \\
+ \>\> DO\_OCEANIC\_PHYS \\
  \> THERMODYNAMICS \\
  \>\> CALC\_GT \\
  \>\>\> GAD\_CALC\_RHS \` $G_\theta^n = G_\theta( u, \theta^n )$ (\ref{eq:Gt-n-sync})\\
  \>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\
  \>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:Gt-n+5-sync}) \\
- \>\> TIMESTEP\_TRACER \` $\tau^*$ (\ref{eq:tstar-sync}) \\
+ \>\> TIMESTEP\_TRACER \` $\theta^*$ (\ref{eq:tstar-sync}) \\
- \>\> IMPLDIFF \` $\tau^{(n+1)}$ (\ref{eq:t-n+1-sync}) \\
+ \>\> IMPLDIFF \` $\theta^{(n+1)}$ (\ref{eq:t-n+1-sync}) \\
  \> DYNAMICS \\
  \>\> CALC\_PHI\_HYD \` $\phi_{hyd}^n$ (\ref{eq:phi-hyd-sync}) \\
  \>\> MOM\_FLUXFORM or MOM\_VECINV \` $G_{\vec{\bf v}}^n$ (\ref{eq:Gv-n-sync})\\
  \>\> TIMESTEP \` $\vec{\bf v}^*$ (\ref{eq:Gv-n+5-sync}, \ref{eq:vstar-sync}) \\
  \>\> IMPLDIFF \` $\vec{\bf v}^{**}$ (\ref{eq:vstarstar-sync}) \\
+ \> UPDATE\_R\_STAR or UPDATE\_SURF\_DR \` (NonLin-FS only)\\
  \> SOLVE\_FOR\_PRESSURE \\
  \>\> CALC\_DIV\_GHAT \` $\eta^*$ (\ref{eq:nstar-sync}) \\
  \>\> CG2D \` $\eta^{n+1}$ (\ref{eq:elliptic-sync}) \\
- \> THE\_CORRECTION\_STEP  \\
+ \> MOMENTUM\_CORRECTION\_STEP  \\
  \>\> CALC\_GRAD\_PHI\_SURF \` $\nabla \eta^{n+1}$ \\
- \>\> 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})\\
+ \> TRACERS\_CORRECTION\_STEP  \\
+ \>\> CYCLE\_TRACER \` $\theta^{n+1}$ \\
+ \>\> FILTER \` Shapiro Filter, Zonal Filter (FFT) \\
+ \>\> CONVECTIVE\_ADJUSTMENT \` \\
  \end{tabbing} \end{minipage} } \end{center}
  \caption{
  Calling tree for the overall synchronous algorithm using
- Adams-Bashforth time-stepping.}
+ Adams-Bashforth time-stepping.
+ The place where the model geometry
+ ({\bf hFac} factors) is updated is added here but is only relevant
+ for the non-linear free-surface algorithm.
+ For completeness, the external forcing,
+ ocean and atmospheric physics have been added, although they are mainly
+ optional}
  \label{fig:call-tree-adams-bashforth-sync}
  \end{figure}
- The Adams-Bashforth extrapolation of explicit tendancies fits neatly
+ The Adams-Bashforth extrapolation of explicit tendencies fits neatly
  into the pressure method algorithm when all state variables are
- co-locacted in time. Fig.~\ref{fig:adams-bashforth-sync} illustrates
+ co-located in time. Fig.~\ref{fig:adams-bashforth-sync} illustrates
  the location of variables in time and the evolution of the algorithm
  with time. The algorithm can be represented by the sequential solution
  of the follow equations:
-Line 442 
 G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil
+Line 499 
 G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil
  \label{eq:vstar-sync} \\
  \vec{\bf v}^{**} & = & {\cal L}_{\vec{\bf v}}^{-1} ( \vec{\bf v}^* )
  \label{eq:vstarstar-sync} \\
- \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
    \nabla \cdot H \widehat{ \vec{\bf v}^{**} }
  \label{eq:nstar-sync} \\
- \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}
- & = & - \frac{\eta^*}{\Delta t^2}
+ ~ = ~ - \frac{\eta^*}{\Delta t^2}
  \label{eq:elliptic-sync} \\
  \vec{\bf v}^{n+1} & = & \vec{\bf v}^{*} - \Delta t g \nabla \eta^{n+1}
  \label{eq:v-n+1-sync}
  \end{eqnarray}
  Fig.~\ref{fig:adams-bashforth-sync} illustrates the location of
  variables in time and evolution of the algorithm with time. The
- Adams-Bashforth extrapolation of the tracer tendancies is illustrated
+ Adams-Bashforth extrapolation of the tracer tendencies is illustrated
- byt the dashed arrow, the prediction at $n+1$ is indicated by the
+ by the dashed arrow, the prediction at $n+1$ is indicated by the
  solid arc. Inversion of the implicit terms, ${\cal
  L}^{-1}_{\theta,S}$, then yields the new tracer fields at $n+1$. All
  these operations are carried out in subroutine {\em THERMODYNAMICS} an
-Line 465 
 accelerations, stepping forward and solv
+Line 522 
 accelerations, stepping forward and solv
  surface pressure gradient terms, corresponding to equations
  \ref{eq:Gv-n-sync} to \ref{eq:v-n+1-sync}.
  These operations are carried out in subroutines {\em DYNAMCIS}, {\em
- SOLVE\_FOR\_PRESSURE} and {\em THE\_CORRECTION\_STEP}. This, then,
+ SOLVE\_FOR\_PRESSURE} and {\em MOMENTUM\_CORRECTION\_STEP}. This, then,
  represents an entire algorithm for stepping forward the model one
  time-step. The corresponding calling tree is given in
  \ref{fig:call-tree-adams-bashforth-sync}.
  \section{Staggered baroclinic time-stepping}
  \label{sect:adams-bashforth-staggered}
+ \begin{rawhtml}
+ <!-- CMIREDIR:adams_bashforth_staggered: -->
+ \end{rawhtml}
  \begin{figure}
  \begin{center}
-Line 480 
 time-step. The corresponding calling tre
+Line 540 
 time-step. The corresponding calling tre
  \caption{
  A schematic of the explicit Adams-Bashforth and implicit time-stepping
  phases of the algorithm but with staggering in time of thermodynamic
- variables with the flow. Explicit thermodynamics tendancies are
+ variables with the flow.
- 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
- 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$.
- explicit tendancy 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
- extrapolate tendancies to $n$ (dashed arrow). This extrapolated
+ extrapolate tendencies to $n$ (dashed arrow).
- tendancy allows thermo-dynamics variables to be stably integrated
+ The hydrostatic pressure/geo-potential $\phi_{hyd}$ is evaluated directly
- 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
- 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).
- 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
- level $n+1/2$. These are then used to calculate the hydrostatic
+ then applied to the previous estimation of the the flow field ($*$-variables)
- pressure/geo-potential, $\phi_{hyd}$ (vertical arrows). The
+ and yields to the two velocity components $u,v$ at time level $n+1/2$.
- hydrostatic pressure gradient is evaluated directly an time level
+ These are then used to calculate the advection term (dashed arc-arrow)
- $n+1/2$ in stepping forward the flow variables from $n$ to $n+1$
+ of the thermo-dynamics tendencies at time step $n$.
- (solid arc-arrow). }
+ The extrapolated thermodynamics tendency, from time level $n-1$ and $n$
+ to $n+1/2$, allows thermodynamic variables to be stably integrated
+ forward-in-time (solid arc-arrow) up to time level $n+1$.
+ }
  \label{fig:adams-bashforth-staggered}
  \end{figure}
-Line 502 
 limiting process for determining a stabl
+Line 565 
 limiting process for determining a stabl
  circumstance, it is more efficient to stagger in time the
  thermodynamic variables with the flow
  variables. Fig.~\ref{fig:adams-bashforth-staggered} illustrates the
- staggering and algorith. The key difference between this and
+ staggering and algorithm. The key difference between this and
- Fig.~\ref{fig:adams-bashforth-sync} is that the new thermodynamics
+ Fig.~\ref{fig:adams-bashforth-sync} is that the thermodynamic variables
- fields are used to compute the hydrostatic pressure at time level
+ are solved after the dynamics, using the recently updated flow field.
- $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
  time giving second order accuracy and more stability.
- The essential change in the staggered algorithm is the calculation of
+ The essential change in the staggered algorithm is that the
- hydrostatic pressure which, in the context of the synchronous
+ thermodynamics solver is delayed from half a time step,
- algorithm involves replacing equation \ref{eq:phi-hyd-sync} with
+ allowing the use of the most recent velocities to compute
- \begin{displaymath}
+ the advection terms. Once the thermodynamics fields are
- \phi_{hyd}^n = \int b(\theta^{n+1},S^{n+1}) dr
+ updated, the hydrostatic pressure is computed
- \end{displaymath}
+ to step forwrad the dynamics.
- but the pressure gradient must also be taken out of the
+ Note that the pressure gradient must also be taken out of the
- Adams-Bashforth extrapoltion. Also, retaining the integer time-levels,
+ Adams-Bashforth extrapolation. Also, retaining the integer time-levels,
  $n$ and $n+1$, does not give a user the sense of where variables are
  located in time.  Instead, we re-write the entire algorithm,
  \ref{eq:Gt-n-sync} to \ref{eq:v-n+1-sync}, annotating the
  position in time of variables appropriately:
  \begin{eqnarray}
- 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
  \label{eq:phi-hyd-staggered} \\
- \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} )
  \label{eq:Gv-n-staggered} \\
- \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}
  \label{eq:Gv-n+5-staggered} \\
- \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)
  \label{eq:vstar-staggered} \\
  \vec{\bf v}^{**} & = & {\cal L}_{\vec{\bf v}}^{-1} ( \vec{\bf v}^* )
  \label{eq:vstarstar-staggered} \\
- \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
    \nabla \cdot H \widehat{ \vec{\bf v}^{**} }
  \label{eq:nstar-staggered} \\
- \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}
- & = & - \frac{\eta^*}{\Delta t^2}
+ ~ = ~ - \frac{\eta^*}{\Delta t^2}
  \label{eq:elliptic-staggered} \\
- \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}
- \label{eq:v-n+1-staggered}
+ \label{eq:v-n+1-staggered} \\
+ G_{\theta,S}^{n} & = & G_{\theta,S} ( u^{n+1/2}, \theta^{n}, S^{n} )
+ \label{eq:Gt-n-staggered} \\
+ G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsilon_{AB}) G_{\theta,S}^{n}-(1/2+\epsilon_{AB}) G_{\theta,S}^{n-1}
+ \label{eq:Gt-n+5-staggered} \\
+ (\theta^*,S^*) & = & (\theta^{n},S^{n}) + \Delta t G_{\theta,S}^{(n+1/2)}
+ \label{eq:tstar-staggered} \\
+ (\theta^{n+1},S^{n+1}) & = & {\cal L}^{-1}_{\theta,S} (\theta^*,S^*)
+ \label{eq:t-n+1-staggered}
  \end{eqnarray}
- The calling sequence is unchanged from
+ The corresponding calling tree is given in
- Fig.~\ref{fig:call-tree-adams-bashforth-sync}. The staggered algorithm
+ \ref{fig:call-tree-adams-bashforth-staggered}.
- is activated with the run-time flag {\bf staggerTimeStep=.TRUE.} in
+ The staggered algorithm is activated with the run-time flag
- {\em PARM01} of {\em data}.
+ {\bf staggerTimeStep}{\em=.TRUE.} in parameter file {\em data},
+ namelist {\em PARM01}.
+ \begin{figure}
+ \begin{center} \fbox{ \begin{minipage}{4.7in} \begin{tabbing}
+ aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill
+ FORWARD\_STEP \\
+ \>\> EXTERNAL\_FIELDS\_LOAD\\
+ \>\> DO\_ATMOSPHERIC\_PHYS \\
+ \>\> DO\_OCEANIC\_PHYS \\
+ \> DYNAMICS \\
+ \>\> CALC\_PHI\_HYD \` $\phi_{hyd}^n$ (\ref{eq:phi-hyd-staggered}) \\
+ \>\> MOM\_FLUXFORM or MOM\_VECINV \` $G_{\vec{\bf v}}^{n-1/2}$
+     (\ref{eq:Gv-n-staggered})\\
+ \>\> TIMESTEP \` $\vec{\bf v}^*$ (\ref{eq:Gv-n+5-staggered},
+                                   \ref{eq:vstar-staggered}) \\
+ \>\> IMPLDIFF \` $\vec{\bf v}^{**}$ (\ref{eq:vstarstar-staggered}) \\
+ \> UPDATE\_R\_STAR or UPDATE\_SURF\_DR \` (NonLin-FS only)\\
+ \> SOLVE\_FOR\_PRESSURE \\
+ \>\> CALC\_DIV\_GHAT \` $\eta^*$ (\ref{eq:nstar-staggered}) \\
+ \>\> CG2D \` $\eta^{n+1/2}$ (\ref{eq:elliptic-staggered}) \\
+ \> MOMENTUM\_CORRECTION\_STEP  \\
+ \>\> CALC\_GRAD\_PHI\_SURF \` $\nabla \eta^{n+1/2}$ \\
+ \>\> CORRECTION\_STEP \` $u^{n+1/2}$,$v^{n+1/2}$ (\ref{eq:v-n+1-staggered})\\
+ \> THERMODYNAMICS \\
+ \>\> CALC\_GT \\
+ \>\>\> GAD\_CALC\_RHS \` $G_\theta^n = G_\theta( u, \theta^n )$
+      (\ref{eq:Gt-n-staggered})\\
+ \>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\
+ \>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:Gt-n+5-staggered}) \\
+ \>\> TIMESTEP\_TRACER \` $\theta^*$ (\ref{eq:tstar-staggered}) \\
+ \>\> IMPLDIFF \` $\theta^{(n+1)}$ (\ref{eq:t-n+1-staggered}) \\
+ \> TRACERS\_CORRECTION\_STEP  \\
+ \>\> CYCLE\_TRACER \` $\theta^{n+1}$ \\
+ \>\> FILTER \` Shapiro Filter, Zonal Filter (FFT) \\
+ \>\> CONVECTIVE\_ADJUSTMENT \` \\
+ \end{tabbing} \end{minipage} } \end{center}
+ \caption{
+ Calling tree for the overall staggered algorithm using
+ Adams-Bashforth time-stepping.
+ The place where the model geometry
+ ({\bf hFac} factors) is updated is added here but is only relevant
+ for the non-linear free-surface algorithm.
+ }
+ \label{fig:call-tree-adams-bashforth-staggered}
+ \end{figure}
  The only difficulty with this approach is apparent in equation
- $\ref{eq:Gt-n-staggered}$ and illustrated by the dotted arrow
+ \ref{eq:Gt-n-staggered} and illustrated by the dotted arrow
- 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
  tracers around is not naturally located in time. This could be avoided
  by applying the Adams-Bashforth extrapolation to the tracer field
- itself and advection that around but this is not yet available. We're
+ itself and advecting that around but this approach is not yet
- not aware of any detrimental effect of this feature. The difficulty
+ available. We're not aware of any detrimental effect of this
- lies mainly in interpretation of what time-level variables and terms
+ feature. The difficulty lies mainly in interpretation of what
- correspond to.
+ time-level variables and terms correspond to.
  \section{Non-hydrostatic formulation}
  \label{sect:non-hydrostatic}
+ \begin{rawhtml}
+ <!-- CMIREDIR:non-hydrostatic_formulation: -->
+ \end{rawhtml}
+ The non-hydrostatic formulation re-introduces the full vertical
+ momentum equation and requires the solution of a 3-D elliptic
+ equations for non-hydrostatic pressure perturbation. We still
+ intergrate vertically for the hydrostatic pressure and solve a 2-D
+ elliptic equation for the surface pressure/elevation for this reduces
+ the amount of work needed to solve for the non-hydrostatic pressure.
- [to be written...]
+ The momentum equations are discretized in time as follows:
+ \begin{eqnarray}
+ \frac{1}{\Delta t} u^{n+1} + g \partial_x \eta^{n+1} + \partial_x \phi_{nh}^{n+1}
+ & = & \frac{1}{\Delta t} u^{n} + G_u^{(n+1/2)} \label{eq:discrete-time-u-nh} \\
+ \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} + G_v^{(n+1/2)} \label{eq:discrete-time-v-nh} \\
+ \frac{1}{\Delta t} w^{n+1} + \partial_r \phi_{nh}^{n+1}
+ & = & \frac{1}{\Delta t} w^{n} + G_w^{(n+1/2)} \label{eq:discrete-time-w-nh}
+ \end{eqnarray}
+ which must satisfy the discrete-in-time depth integrated continuity,
+ equation~\ref{eq:discrete-time-backward-free-surface} and the local continuity equation
+ \begin{equation}
+ \partial_x u^{n+1} + \partial_y v^{n+1} + \partial_r w^{n+1} = 0
+ \label{eq:non-divergence-nh}
+ \end{equation}
+ As before, the explicit predictions for momentum are consolidated as:
+ \begin{eqnarray*}
+ u^* & = & u^n + \Delta t G_u^{(n+1/2)} \\
+ v^* & = & v^n + \Delta t G_v^{(n+1/2)} \\
+ w^* & = & w^n + \Delta t G_w^{(n+1/2)}
+ \end{eqnarray*}
+ but this time we introduce an intermediate step by splitting the
+ tendancy of the flow as follows:
+ \begin{eqnarray}
+ u^{n+1} = u^{**} - \Delta t \partial_x \phi_{nh}^{n+1}
+ & &
+ u^{**} = u^{*} - \Delta t g \partial_x \eta^{n+1} \\
+ v^{n+1} = v^{**} - \Delta t \partial_y \phi_{nh}^{n+1}
+ & &
+ v^{**} = v^{*} - \Delta t g \partial_y \eta^{n+1}
+ \end{eqnarray}
+ Substituting into the depth integrated continuity
+ (equation~\ref{eq:discrete-time-backward-free-surface}) gives
+ \begin{equation}
+ \partial_x H \partial_x \left( g \eta^{n+1} + \widehat{\phi}_{nh}^{n+1} \right)
+ +
+ \partial_y H \partial_y \left( g \eta^{n+1} + \widehat{\phi}_{nh}^{n+1} \right)
+  - \frac{\epsilon_{fs}\eta^*}{\Delta t^2}
+ = - \frac{\eta^*}{\Delta t^2}
+ \end{equation}
+ which is approximated by equation
+ \ref{eq:elliptic-backward-free-surface} on the basis that i)
+ $\phi_{nh}^{n+1}$ is not yet known and ii) $\nabla \widehat{\phi}_{nh}
+ << g \nabla \eta$. If \ref{eq:elliptic-backward-free-surface} is
+ solved accurately then the implication is that $\widehat{\phi}_{nh}
+ \approx 0$ so that thet non-hydrostatic pressure field does not drive
+ barotropic motion.
+ The flow must satisfy non-divergence
+ (equation~\ref{eq:non-divergence-nh}) locally, as well as depth
+ integrated, and this constraint is used to form a 3-D elliptic
+ equations for $\phi_{nh}^{n+1}$:
+ \begin{equation}
+ \partial_{xx} \phi_{nh}^{n+1} + \partial_{yy} \phi_{nh}^{n+1} +
+ \partial_{rr} \phi_{nh}^{n+1} =
+ \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*}
+ \end{equation}
+ The entire algorithm can be summarized as the sequential solution of
+ the following equations:
+ \begin{eqnarray}
+ u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2)} \label{eq:ustar-nh} \\
+ v^{*} & = & v^{n} + \Delta t G_v^{(n+1/2)} \label{eq:vstar-nh} \\
+ w^{*} & = & w^{n} + \Delta t G_w^{(n+1/2)} \label{eq:wstar-nh} \\
+ \eta^* ~ = ~ \epsilon_{fs} \left( \eta^{n} + \Delta t (P-E) \right)
+ & - & \Delta t
+   \partial_x H \widehat{u^{*}}
+ + \partial_y H \widehat{v^{*}}
+ \\
+   \partial_x g H \partial_x \eta^{n+1}
+ + \partial_y g H \partial_y \eta^{n+1}
+ & - & \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2}
+ ~ = ~
+ - \frac{\eta^*}{\Delta t^2}
+ \label{eq:elliptic-nh}
+ \\
+ u^{**} & = & u^{*} - \Delta t g \partial_x \eta^{n+1} \label{eq:unx-nh}\\
+ v^{**} & = & v^{*} - \Delta t g \partial_y \eta^{n+1} \label{eq:vnx-nh}\\
+ \partial_{xx} \phi_{nh}^{n+1} + \partial_{yy} \phi_{nh}^{n+1} +
+ \partial_{rr} \phi_{nh}^{n+1} & = &
+ \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*} \\
+ u^{n+1} & = & u^{**} - \Delta t \partial_x \phi_{nh}^{n+1} \label{eq:un+1-nh}\\
+ v^{n+1} & = & v^{**} - \Delta t \partial_y \phi_{nh}^{n+1} \label{eq:vn+1-nh}\\
+ \partial_r w^{n+1} & = & - \partial_x u^{n+1} - \partial_y v^{n+1}
+ \end{eqnarray}
+ where the last equation is solved by vertically integrating for
+ $w^{n+1}$.
  \section{Variants on the Free Surface}
+ \label{sect:free-surface}
- We now descibe the various formulations of the free-surface that
+ We now describe the various formulations of the free-surface that
  include non-linear forms, implicit in time using Crank-Nicholson,
  explicit and [one day] split-explicit. First, we'll reiterate the
  underlying algorithm but this time using the notation consistent with
  the more general vertical coordinate $r$. The elliptic equation for
- free-surface coordinate (units of $r$), correpsonding to
+ free-surface coordinate (units of $r$), corresponding to
  \ref{eq:discrete-time-backward-free-surface}, and
  assuming no non-hydrostatic effects ($\epsilon_{nh} = 0$) is:
  \begin{eqnarray}
-Line 592 
 where
+Line 796 
 where
  \begin{eqnarray}
  {\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
  \Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o} \vec{\bf v}^* dr
- \: + \: \epsilon_{fw} \Delta_t (P-E)^{n}
+ \: + \: \epsilon_{fw} \Delta t (P-E)^{n}
  \label{eq-solve2D_rhs}
  \end{eqnarray}
  \fbox{ \begin{minipage}{4.75in}
  {\em S/R SOLVE\_FOR\_PRESSURE} ({\em solve\_for\_pressure.F})
- $u^*$: {\bf GuNm1} ({\em DYNVARS.h})
+ $u^*$: {\bf gU} ({\em DYNVARS.h})
- $v^*$: {\bf GvNm1} ({\em DYNVARS.h})
+ $v^*$: {\bf gV} ({\em DYNVARS.h})
  $\eta^*$: {\bf cg2d\_b} (\em SOLVE\_FOR\_PRESSURE.h)
-Line 611 
 $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)
+Line 815 
 $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)
  Once ${\eta}^{n+1}$ has been found, substituting into
- \ref{eq-tDsC-Hmom} 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}$
- hydrostatic ($\epsilon_{nh}=0$):
+ if the model is hydrostatic ($\epsilon_{nh}=0$):
  $$
  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}
  - \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
-Line 621 
 $$
+Line 825 
 $$
  This is known as the correction step. However, when the model is
  non-hydrostatic ($\epsilon_{nh}=1$) we need an additional step and an
  additional equation for $\phi'_{nh}$. This is obtained by substituting
- \ref{eq-tDsC-Hmom} and \ref{eq-tDsC-Vmom} into
+ \ref{eq:discrete-time-u-nh}, \ref{eq:discrete-time-v-nh} and \ref{eq:discrete-time-w-nh}
- \ref{eq-tDsC-cont}:
+ into continuity:
  \begin{equation}
  \left[ {\bf \nabla}_h^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1}
  = \frac{1}{\Delta t} \left(
-Line 653 
 without any consequence on the solution.
+Line 857 
 without any consequence on the solution.
  $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)
- $\phi_{nh}^{n+1}$: {\bf phi\_nh} (\em DYNVARS.h)
+ $\phi_{nh}^{n+1}$: {\bf phi\_nh} (\em NH\_VARS.h)
- $u^*$: {\bf GuNm1} ({\em DYNVARS.h})
+ $u^*$: {\bf gU} ({\em DYNVARS.h})
- $v^*$: {\bf GvNm1} ({\em DYNVARS.h})
+ $v^*$: {\bf gV} ({\em DYNVARS.h})
  $u^{n+1}$: {\bf uVel} ({\em DYNVARS.h})
-Line 685 
 at the same point in the code.
+Line 889 
 at the same point in the code.
  \subsection{Crank-Nickelson barotropic time stepping}
+ \label{sect:freesurf-CrankNick}
  The full implicit time stepping described previously is
  unconditionally stable but damps the fast gravity waves, resulting in
-Line 699 
 stable, Crank-Nickelson scheme; $(\beta,
+Line 904 
 stable, Crank-Nickelson scheme; $(\beta,
  corresponds to the forward - backward scheme that conserves energy but is
  only stable for small time steps.\\
  In the code, $\beta,\gamma$ are defined as parameters, respectively
- {\it implicSurfPress}, {\it implicDiv2DFlow}. They are read from
+ {\bf implicSurfPress}, {\bf implicDiv2DFlow}. They are read from
- 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.
- Equations \ref{eq-tDsC-Hmom} and \ref{eq-tDsC-eta} are modified as follows:
+ Equations \ref{eq:ustar-backward-free-surface} --
- $$
+ \ref{eq:vn+1-backward-free-surface} are modified as follows:
+ \begin{eqnarray*}
  \frac{ \vec{\bf v}^{n+1} }{ \Delta t }
  + {\bf \nabla}_h b_s [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ]
  + \epsilon_{nh} {\bf \nabla}_h {\phi'_{nh}}^{n+1}
   = \frac{ \vec{\bf v}^* }{ \Delta t }
- $$
+ \end{eqnarray*}
- $$
+ \begin{eqnarray}
  \epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t}
  + {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o}
  [ \gamma \vec{\bf v}^{n+1} + (1-\gamma) \vec{\bf v}^{n}] dr
  = \epsilon_{fw} (P-E)
- $$
+ \label{eq:eta-n+1-CrankNick}
+ \end{eqnarray}
  where:
  \begin{eqnarray*}
  \vec{\bf v}^* & = &
-Line 736 
 $$
+Line 943 
 $$
  {\bf \nabla}_h {\eta}^{n+1}
  = {\eta}^*
  $$
- and then to compute (correction step):
+ and then to compute ({\em CORRECTION\_STEP}):
  $$
  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}
  - \beta \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
  $$
- The non-hydrostatic part is solved as described previously.
+ %The non-hydrostatic part is solved as described previously.
- Note that:
+ \noindent
+ Notes:
  \begin{enumerate}
+ \item The RHS term of equation \ref{eq:eta-n+1-CrankNick}
+ corresponds the contribution of fresh water flux (P-E)
+ to the free-surface variations ($\epsilon_{fw}=1$,
+ {\bf useRealFreshWater}{\em=TRUE} in parameter file {\em data}).
+ In order to remain consistent with the tracer equation, specially in
+ the non-linear free-surface formulation, this term is also
+ affected by the Crank-Nickelson time stepping. The RHS reads:
+ $\epsilon_{fw} ( \gamma (P-E)^{n+1/2} + (1-\gamma) (P-E)^{n-1/2} )$
  \item The non-hydrostatic part of the code has not yet been
- 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)$.
  \item The stability criteria with Crank-Nickelson time stepping
  for the pure linear gravity wave problem in cartesian coordinates is:
  \begin{itemize}

 Legend:



Removed from v.1.8
 


changed lines


 
Added in v.1.25
 Legend:



Removed from v.1.8
 


changed lines


 
Added in v.1.25
-Removed from v.1.8
+Added in v.1.25

	ViewVC Help
Powered by ViewVC 1.1.22