/[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.9 by cnh,
Thu Oct 25 18:36:53 2001 UTC
+revision 1.16 by edhill,
Thu Aug  7 18:27:52 2003 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 \\
-Line 80 
 FORWARD\_STEP \\
+Line 91 
 FORWARD\_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 554 
 is activated with the run-time flag {\bf
+Line 567 
 is activated with the run-time flag {\bf
  {\em PARM01} of {\em data}.
  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
  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.
+ 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} +P-E+R \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 611 
 $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)
+Line 715 
 $\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,eq:discrete-time-v} yields $\vec{\bf v}^{n+1}$ if the model is
  hydrostatic ($\epsilon_{nh}=0$):
  $$
  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}
-Line 621 
 $$
+Line 725 
 $$
  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 702 
 In the code, $\beta,\gamma$ are defined
+Line 806 
 In the code, $\beta,\gamma$ are defined
  {\it implicSurfPress}, {\it implicDiv2DFlow}. They are read from
  the main data file "{\it 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:
  $$
  \frac{ \vec{\bf v}^{n+1} }{ \Delta t }
  + {\bf \nabla}_h b_s [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ]

 Legend:



Removed from v.1.9
 


changed lines


 
Added in v.1.16
 Legend:



Removed from v.1.9
 


changed lines


 
Added in v.1.16
-Removed from v.1.9
+Added in v.1.16

	ViewVC Help
Powered by ViewVC 1.1.22