/[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.13 by adcroft,
Tue Nov 13 20:51:36 2001 UTC
+revision 1.17 by jmc,
Wed Oct 13 18:50:54 2004 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.
+ \section{Time-stepping}
  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 70 
 temporary.}
+Line 81 
 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 algorihtm}
+ \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 162 
 The correspondence to the code is as fol
+Line 174 
 The correspondence to the code is as fol
  \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 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 281 
 FORWARD\_STEP \\
+Line 294 
 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 333 
 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 330 
 A stability analysis for a relaxation eq
+Line 348 
 A stability analysis for a relaxation eq
  \section{Implicit time-stepping: backward method}
  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.
+ 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} =
-Line 361 
 Fig.~\ref{fig:call-tree-adams-bashforth}
+Line 383 
 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 389 
 is solved to yield the state variables a
+Line 411 
 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})\\
-Line 404 
 FORWARD\_STEP \\
+Line 429 
 FORWARD\_STEP \\
  \>\> 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
+ ({\em 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}
-Line 442 
 G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil
+Line 478 
 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}
-Line 465 
 accelerations, stepping forward and solv
+Line 501 
 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}.
-Line 481 
 time-step. The corresponding calling tre
+Line 517 
 time-step. The corresponding calling tre
  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 tendencies are
- evaluated at time level $n-1/2$ as a function of the thermodynamics
+ evaluated at time level $n$ as a function of the thermodynamics
- state at that time level $n$ and flow at time $n$ (dotted arrow). The
+ state at that time level $n$ and flow at time $n+1/2$ (dotted arrow). The
- explicit tendency from the previous time level, $n-3/2$, is used to
+ explicit tendency from the previous time level, $n-1$, is used to
- extrapolate tendencies to $n$ (dashed arrow). This extrapolated
+ extrapolate tendencies to $n+1/2$ (dashed arrow). This extrapolated
  tendency allows thermo-dynamics variables to be stably integrated
- forward-in-time to render an estimate ($*$-variables) at the $n+1/2$
+ forward-in-time to render an estimate ($*$-variables) at the $n+1$
  time level (solid arc-arrow). The implicit-in-time operator ${\cal
  L_{\theta,S}}$ is solved to yield the thermodynamic variables at time
- level $n+1/2$. These are then used to calculate the hydrostatic
+ level $n+1$. These are then used to calculate the hydrostatic
  pressure/geo-potential, $\phi_{hyd}$ (vertical arrows). The
- hydrostatic pressure gradient is evaluated directly an time level
+ hydrostatic pressure gradient is evaluated directly at time level
- $n+1/2$ in stepping forward the flow variables from $n$ to $n+1$
+ $n+1$ in stepping forward the flow variables from $n+1/2$ to $n+3/2$
  (solid arc-arrow). }
  \label{fig:adams-bashforth-staggered}
  \end{figure}
-Line 503 
 circumstance, it is more efficient to st
+Line 539 
 circumstance, it is more efficient to st
  thermodynamic variables with the flow
  variables. Fig.~\ref{fig:adams-bashforth-staggered} illustrates the
  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 frowrad 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 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} )
+ \vec{\bf G}_{\vec{\bf v}}^{n-1/2} & = & \vec{\bf G}_{\vec{\bf v}} ( \vec{\bf v}^{n-1/2} )
- \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 )
  \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)
+ \phi^{n}_{hyd} & = & \int b(\theta^{n},S^{n}) dr
+ \label{eq:phi-hyd-staggered} \\
+ \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=.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 \` $\tau^*$ (\ref{eq:tstar-staggered}) \\
+ \>\> IMPLDIFF \` $\tau^{(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
+ ({\em 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
  connecting $u,v^n$ with $G_\theta^{n-1/2}$. 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}
- [to be written...]
+ 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.
- Equation for $w^{n+1}$ will be here as will 3-D elliptic equations.
+ The momentum equations are discretized in time as follows:
- \label{eq:discrete-time-w}
+ \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}$.
-Line 623 
 $$
+Line 794 
 $$
  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: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}
  into continuity:
  \begin{equation}
  \left[ {\bf \nabla}_h^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1}

 Legend:



Removed from v.1.13
 


changed lines


 
Added in v.1.17
 Legend:



Removed from v.1.13
 


changed lines


 
Added in v.1.17
-Removed from v.1.13
+Added in v.1.17

	ViewVC Help
Powered by ViewVC 1.1.22