/[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.8 by adcroft, Fri Oct 5 19:41:08 2001 UTC revision 1.15 by cnh, Thu Feb 28 19:32:19 2002 UTC
# Line 1  Line 1 
1  % $Header$  % $Header$
2  % $Name$  % $Name$
3    
4    This chapter lays out the numerical schemes that are
5    employed in the core MITgcm algorithm. Whenever possible
6    links are made to actual program code in the MITgcm implementation.
7    The chapter begins with a discussion of the temporal discretization
8    used in MITgcm. This discussion is followed by sections that
9    describe the spatial discretization. The schemes employed for momentum
10    terms are described first, afterwards the schemes that apply to
11    passive and dynamically active tracers are described.
12    
13    
14    \section{Time-stepping}
15  The equations of motion integrated by the model involve four  The equations of motion integrated by the model involve four
16  prognostic equations for flow, $u$ and $v$, temperature, $\theta$, and  prognostic equations for flow, $u$ and $v$, temperature, $\theta$, and
17  salt/moisture, $S$, and three diagnostic equations for vertical flow,  salt/moisture, $S$, and three diagnostic equations for vertical flow,
# Line 13  and diagnostic equations requires a mode Line 24  and diagnostic equations requires a mode
24  forward prognostic variables while satisfying constraints imposed by  forward prognostic variables while satisfying constraints imposed by
25  diagnostic equations.  diagnostic equations.
26    
27  Since the model comes in several flavours and formulation, it would be  Since the model comes in several flavors and formulation, it would be
28  confusing to present the model algorithm exactly as written into code  confusing to present the model algorithm exactly as written into code
29  along with all the switches and optional terms. Instead, we present  along with all the switches and optional terms. Instead, we present
30  the algorithm for each of the basic formulations which are:  the algorithm for each of the basic formulations which are:
# Line 37  method, with a rigid-lid model in sectio Line 48  method, with a rigid-lid model in sectio
48  \ref{sect:pressure-method-rigid-lid}. This algorithm is essentially  \ref{sect:pressure-method-rigid-lid}. This algorithm is essentially
49  unchanged, apart for some coefficients, when the rigid lid assumption  unchanged, apart for some coefficients, when the rigid lid assumption
50  is replaced with a linearized implicit free-surface, described in  is replaced with a linearized implicit free-surface, described in
51  section \ref{sect:pressure-method-linear-backward}. These two flavours  section \ref{sect:pressure-method-linear-backward}. These two flavors
52  of the pressure-method encompass all formulations of the model as it  of the pressure-method encompass all formulations of the model as it
53  exists today. The integration of explicit in time terms is out-lined  exists today. The integration of explicit in time terms is out-lined
54  in section \ref{sect:adams-bashforth} and put into the context of the  in section \ref{sect:adams-bashforth} and put into the context of the
# Line 61  A schematic of the evolution in time of Line 72  A schematic of the evolution in time of
72  algorithm. A prediction for the flow variables at time level $n+1$ is  algorithm. A prediction for the flow variables at time level $n+1$ is
73  made based only on the explicit terms, $G^{(n+^1/_2)}$, and denoted  made based only on the explicit terms, $G^{(n+^1/_2)}$, and denoted
74  $u^*$, $v^*$. Next, a pressure field is found such that $u^{n+1}$,  $u^*$, $v^*$. Next, a pressure field is found such that $u^{n+1}$,
75  $v^{n+1}$ will be non-divergent. Conceptually, the $*$ quantitites  $v^{n+1}$ will be non-divergent. Conceptually, the $*$ quantities
76  exist at time level $n+1$ but they are intermediate and only  exist at time level $n+1$ but they are intermediate and only
77  temporary.}  temporary.}
78  \label{fig:pressure-method-rigid-lid}  \label{fig:pressure-method-rigid-lid}
# Line 70  temporary.} Line 81  temporary.}
81  \begin{figure}  \begin{figure}
82  \begin{center} \fbox{ \begin{minipage}{4.5in} \begin{tabbing}  \begin{center} \fbox{ \begin{minipage}{4.5in} \begin{tabbing}
83  aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill  aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill
84  FORWARD\_STEP \\  \proclink{FORWARD\_STEP}{../code/._model_src_forward_step.F} \\
85  \> DYNAMICS \\  \> DYNAMICS \\
86  \>\> TIMESTEP \` $u^*$,$v^*$ (\ref{eq:ustar-rigid-lid},\ref{eq:vstar-rigid-lid}) \\  \>\> TIMESTEP \` $u^*$,$v^*$ (\ref{eq:ustar-rigid-lid},\ref{eq:vstar-rigid-lid}) \\
87  \> SOLVE\_FOR\_PRESSURE \\  \> SOLVE\_FOR\_PRESSURE \\
# Line 80  FORWARD\_STEP \\ Line 91  FORWARD\_STEP \\
91  \>\> CALC\_GRAD\_PHI\_SURF \` $\nabla \eta^{n+1}$ \\  \>\> CALC\_GRAD\_PHI\_SURF \` $\nabla \eta^{n+1}$ \\
92  \>\> CORRECTION\_STEP \` $u^{n+1}$,$v^{n+1}$ (\ref{eq:un+1-rigid-lid},\ref{eq:vn+1-rigid-lid})  \>\> CORRECTION\_STEP \` $u^{n+1}$,$v^{n+1}$ (\ref{eq:un+1-rigid-lid},\ref{eq:vn+1-rigid-lid})
93  \end{tabbing} \end{minipage} } \end{center}  \end{tabbing} \end{minipage} } \end{center}
94  \caption{Calling tree for the pressure method alogtihm}  \caption{Calling tree for the pressure method algorihtm}
95  \label{fig:call-tree-pressure-method}  \label{fig:call-tree-pressure-method}
96  \end{figure}  \end{figure}
97    
# Line 102  flow boundary conditions applied, become Line 113  flow boundary conditions applied, become
113  \label{eq:rigid-lid-continuity}  \label{eq:rigid-lid-continuity}
114  \end{equation}  \end{equation}
115  Here, $H\widehat{u} = \int_H u dz$ is the depth integral of $u$,  Here, $H\widehat{u} = \int_H u dz$ is the depth integral of $u$,
116  similarly for $H\widehat{v}$. The rigid-lid approzimation sets $w=0$  similarly for $H\widehat{v}$. The rigid-lid approximation sets $w=0$
117  at the lid so that it does not move but allows a pressure to be  at the lid so that it does not move but allows a pressure to be
118  exerted on the fluid by the lid. The horizontal momentum equations and  exerted on the fluid by the lid. The horizontal momentum equations and
119  vertically integrated continuity equation are be discretized in time  vertically integrated continuity equation are be discretized in time
# Line 130  pressure gradient is explicit and so can Line 141  pressure gradient is explicit and so can
141  $G$.  $G$.
142    
143  Substituting the two momentum equations into the depth integrated  Substituting the two momentum equations into the depth integrated
144  continuity equation eliminates $u^{n+1}$ and $v^{n+1}$ yeilding an  continuity equation eliminates $u^{n+1}$ and $v^{n+1}$ yielding an
145  elliptic equation for $\eta^{n+1}$. Equations  elliptic equation for $\eta^{n+1}$. Equations
146  \ref{eq:discrete-time-u}, \ref{eq:discrete-time-v} and  \ref{eq:discrete-time-u}, \ref{eq:discrete-time-v} and
147  \ref{eq:discrete-time-cont-rigid-lid} can then be re-arranged as follows:  \ref{eq:discrete-time-cont-rigid-lid} can then be re-arranged as follows:
# Line 157  Fig.~\ref{fig:pressure-method-rigid-lid} Line 168  Fig.~\ref{fig:pressure-method-rigid-lid}
168  with the future flow field (time level $n+1$) but it could equally  with the future flow field (time level $n+1$) but it could equally
169  have been drawn as staggered in time with the flow.  have been drawn as staggered in time with the flow.
170    
171  The correspondance to the code is as follows:  The correspondence to the code is as follows:
172  \begin{itemize}  \begin{itemize}
173  \item  \item
174  the prognostic phase, equations \ref{eq:ustar-rigid-lid} and \ref{eq:vstar-rigid-lid},  the prognostic phase, equations \ref{eq:ustar-rigid-lid} and \ref{eq:vstar-rigid-lid},
175  stepping forward $u^n$ and $v^n$ to $u^{*}$ and $v^{*}$ is coded in  stepping forward $u^n$ and $v^n$ to $u^{*}$ and $v^{*}$ is coded in
176  {\em TIMESTEP.F}  \proclink{TIMESTEP}{../code/._model_src_timestep.F}
177  \item  \item
178  the vertical integration, $H \widehat{u^*}$ and $H  the vertical integration, $H \widehat{u^*}$ and $H
179  \widehat{v^*}$, divergence and invertion of the elliptic operator in  \widehat{v^*}$, divergence and inversion of the elliptic operator in
180  equation \ref{eq:elliptic} is coded in {\em  equation \ref{eq:elliptic} is coded in
181  SOLVE\_FOR\_PRESSURE.F}  \proclink{SOLVE\_FOR\_PRESSURE}{../code/._model_src_solve_for_pressure.F}
182  \item  \item
183  finally, the new flow field at time level $n+1$ given by equations  finally, the new flow field at time level $n+1$ given by equations
184  \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
185    \proclink{CORRECTION\_STEP}{../code/._model_src_correction_step.F}.
186  \end{itemize}  \end{itemize}
187  The calling tree for these routines is given in  The calling tree for these routines is given in
188  Fig.~\ref{fig:call-tree-pressure-method}.  Fig.~\ref{fig:call-tree-pressure-method}.
# Line 250  to~\ref{eq:vn+1-backward-free-surface}, Line 262  to~\ref{eq:vn+1-backward-free-surface},
262  the pressure method algorithm with a backward implicit, linearized  the pressure method algorithm with a backward implicit, linearized
263  free surface. The method is still formerly a pressure method because  free surface. The method is still formerly a pressure method because
264  in the limit of large $\Delta t$ the rigid-lid method is  in the limit of large $\Delta t$ the rigid-lid method is
265  reovered. However, the implicit treatment of the free-surface allows  recovered. However, the implicit treatment of the free-surface allows
266  the flow to be divergent and for the surface pressure/elevation to  the flow to be divergent and for the surface pressure/elevation to
267  respond on a finite time-scale (as opposed to instantly). To recovere  respond on a finite time-scale (as opposed to instantly). To recover
268  the rigid-lid formulation, we introduced a switch-like parameter,  the rigid-lid formulation, we introduced a switch-like parameter,
269  $\epsilon_{fs}$, which selects between the free-surface and rigid-lid;  $\epsilon_{fs}$, which selects between the free-surface and rigid-lid;
270  $\epsilon_{fs}=1$ allows the free-surface to evolve; $\epsilon_{fs}=0$  $\epsilon_{fs}=1$ allows the free-surface to evolve; $\epsilon_{fs}=0$
# Line 331  A stability analysis for a relaxation eq Line 343  A stability analysis for a relaxation eq
343    
344  Vertical diffusion and viscosity can be treated implicitly in time  Vertical diffusion and viscosity can be treated implicitly in time
345  using the backward method which is an intrinsic scheme. For tracers,  using the backward method which is an intrinsic scheme. For tracers,
346  the time discrretized equation is:  the time discretized equation is:
347  \begin{equation}  \begin{equation}
348  \tau^{n+1} - \Delta t \partial_r \kappa_v \partial_r \tau^{n+1} =  \tau^{n+1} - \Delta t \partial_r \kappa_v \partial_r \tau^{n+1} =
349  \tau^{n} + \Delta t G_\tau^{(n+1/2)}  \tau^{n} + \Delta t G_\tau^{(n+1/2)}
# Line 377  implicit and are thus cast as a an expli Line 389  implicit and are thus cast as a an expli
389  \caption{  \caption{
390  A schematic of the explicit Adams-Bashforth and implicit time-stepping  A schematic of the explicit Adams-Bashforth and implicit time-stepping
391  phases of the algorithm. All prognostic variables are co-located in  phases of the algorithm. All prognostic variables are co-located in
392  time. Explicit tendancies are evaluated at time level $n$ as a  time. Explicit tendencies are evaluated at time level $n$ as a
393  function of the state at that time level (dotted arrow). The explicit  function of the state at that time level (dotted arrow). The explicit
394  tendancy from the previous time level, $n-1$, is used to extrapolate  tendency from the previous time level, $n-1$, is used to extrapolate
395  tendancies to $n+1/2$ (dashed arrow). This extrapolated tendancy  tendencies to $n+1/2$ (dashed arrow). This extrapolated tendency
396  allows variables to be stably integrated forward-in-time to render an  allows variables to be stably integrated forward-in-time to render an
397  estimate ($*$-variables) at the $n+1$ time level (solid  estimate ($*$-variables) at the $n+1$ time level (solid
398  arc-arrow). The operator ${\cal L}$ formed from implicit-in-time terms  arc-arrow). The operator ${\cal L}$ formed from implicit-in-time terms
# Line 417  Adams-Bashforth time-stepping.} Line 429  Adams-Bashforth time-stepping.}
429  \label{fig:call-tree-adams-bashforth-sync}  \label{fig:call-tree-adams-bashforth-sync}
430  \end{figure}  \end{figure}
431    
432  The Adams-Bashforth extrapolation of explicit tendancies fits neatly  The Adams-Bashforth extrapolation of explicit tendencies fits neatly
433  into the pressure method algorithm when all state variables are  into the pressure method algorithm when all state variables are
434  co-locacted in time. Fig.~\ref{fig:adams-bashforth-sync} illustrates  co-located in time. Fig.~\ref{fig:adams-bashforth-sync} illustrates
435  the location of variables in time and the evolution of the algorithm  the location of variables in time and the evolution of the algorithm
436  with time. The algorithm can be represented by the sequential solution  with time. The algorithm can be represented by the sequential solution
437  of the follow equations:  of the follow equations:
# Line 453  G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil Line 465  G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil
465  \end{eqnarray}  \end{eqnarray}
466  Fig.~\ref{fig:adams-bashforth-sync} illustrates the location of  Fig.~\ref{fig:adams-bashforth-sync} illustrates the location of
467  variables in time and evolution of the algorithm with time. The  variables in time and evolution of the algorithm with time. The
468  Adams-Bashforth extrapolation of the tracer tendancies is illustrated  Adams-Bashforth extrapolation of the tracer tendencies is illustrated
469  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
470  solid arc. Inversion of the implicit terms, ${\cal  solid arc. Inversion of the implicit terms, ${\cal
471  L}^{-1}_{\theta,S}$, then yields the new tracer fields at $n+1$. All  L}^{-1}_{\theta,S}$, then yields the new tracer fields at $n+1$. All
472  these operations are carried out in subroutine {\em THERMODYNAMICS} an  these operations are carried out in subroutine {\em THERMODYNAMICS} an
# Line 480  time-step. The corresponding calling tre Line 492  time-step. The corresponding calling tre
492  \caption{  \caption{
493  A schematic of the explicit Adams-Bashforth and implicit time-stepping  A schematic of the explicit Adams-Bashforth and implicit time-stepping
494  phases of the algorithm but with staggering in time of thermodynamic  phases of the algorithm but with staggering in time of thermodynamic
495  variables with the flow. Explicit thermodynamics tendancies are  variables with the flow. Explicit thermodynamics tendencies are
496  evaluated at time level $n-1/2$ as a function of the thermodynamics  evaluated at time level $n-1/2$ as a function of the thermodynamics
497  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$ (dotted arrow). The
498  explicit tendancy from the previous time level, $n-3/2$, is used to  explicit tendency from the previous time level, $n-3/2$, is used to
499  extrapolate tendancies to $n$ (dashed arrow). This extrapolated  extrapolate tendencies to $n$ (dashed arrow). This extrapolated
500  tendancy allows thermo-dynamics variables to be stably integrated  tendency allows thermo-dynamics variables to be stably integrated
501  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/2$
502  time level (solid arc-arrow). The implicit-in-time operator ${\cal  time level (solid arc-arrow). The implicit-in-time operator ${\cal
503  L_{\theta,S}}$ is solved to yield the thermodynamic variables at time  L_{\theta,S}}$ is solved to yield the thermodynamic variables at time
# Line 502  limiting process for determining a stabl Line 514  limiting process for determining a stabl
514  circumstance, it is more efficient to stagger in time the  circumstance, it is more efficient to stagger in time the
515  thermodynamic variables with the flow  thermodynamic variables with the flow
516  variables. Fig.~\ref{fig:adams-bashforth-staggered} illustrates the  variables. Fig.~\ref{fig:adams-bashforth-staggered} illustrates the
517  staggering and algorith. The key difference between this and  staggering and algorithm. The key difference between this and
518  Fig.~\ref{fig:adams-bashforth-sync} is that the new thermodynamics  Fig.~\ref{fig:adams-bashforth-sync} is that the new thermodynamics
519  fields are used to compute the hydrostatic pressure at time level  fields are used to compute the hydrostatic pressure at time level
520  $n+1/2$. The essentially allows the gravity wave terms to leap-frog in  $n+1/2$. The essentially allows the gravity wave terms to leap-frog in
# Line 515  algorithm involves replacing equation \r Line 527  algorithm involves replacing equation \r
527  \phi_{hyd}^n = \int b(\theta^{n+1},S^{n+1}) dr  \phi_{hyd}^n = \int b(\theta^{n+1},S^{n+1}) dr
528  \end{displaymath}  \end{displaymath}
529  but the pressure gradient must also be taken out of the  but the pressure gradient must also be taken out of the
530  Adams-Bashforth extrapoltion. Also, retaining the integer time-levels,  Adams-Bashforth extrapolation. Also, retaining the integer time-levels,
531  $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
532  located in time.  Instead, we re-write the entire algorithm,  located in time.  Instead, we re-write the entire algorithm,
533  \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
# Line 554  is activated with the run-time flag {\bf Line 566  is activated with the run-time flag {\bf
566  {\em PARM01} of {\em data}.  {\em PARM01} of {\em data}.
567    
568  The only difficulty with this approach is apparent in equation  The only difficulty with this approach is apparent in equation
569  $\ref{eq:Gt-n-staggered}$ and illustrated by the dotted arrow  \ref{eq:Gt-n-staggered} and illustrated by the dotted arrow
570  connecting $u,v^n$ with $G_\theta^{n-1/2}$. The flow used to advect  connecting $u,v^n$ with $G_\theta^{n-1/2}$. The flow used to advect
571  tracers around is not naturally located in time. This could be avoided  tracers around is not naturally located in time. This could be avoided
572  by applying the Adams-Bashforth extrapolation to the tracer field  by applying the Adams-Bashforth extrapolation to the tracer field
573  itself and advection that around but this is not yet available. We're  itself and advecting that around but this approach is not yet
574  not aware of any detrimental effect of this feature. The difficulty  available. We're not aware of any detrimental effect of this
575  lies mainly in interpretation of what time-level variables and terms  feature. The difficulty lies mainly in interpretation of what
576  correspond to.  time-level variables and terms correspond to.
577    
578    
579  \section{Non-hydrostatic formulation}  \section{Non-hydrostatic formulation}
580  \label{sect:non-hydrostatic}  \label{sect:non-hydrostatic}
581    
582  [to be written...]  The non-hydrostatic formulation re-introduces the full vertical
583    momentum equation and requires the solution of a 3-D elliptic
584    equations for non-hydrostatic pressure perturbation. We still
585    intergrate vertically for the hydrostatic pressure and solve a 2-D
586    elliptic equation for the surface pressure/elevation for this reduces
587    the amount of work needed to solve for the non-hydrostatic pressure.
588    
589    The momentum equations are discretized in time as follows:
590    \begin{eqnarray}
591    \frac{1}{\Delta t} u^{n+1} + g \partial_x \eta^{n+1} + \partial_x \phi_{nh}^{n+1}
592    & = & \frac{1}{\Delta t} u^{n} + G_u^{(n+1/2)} \label{eq:discrete-time-u-nh} \\
593    \frac{1}{\Delta t} v^{n+1} + g \partial_y \eta^{n+1} + \partial_y \phi_{nh}^{n+1}
594    & = & \frac{1}{\Delta t} v^{n} + G_v^{(n+1/2)} \label{eq:discrete-time-v-nh} \\
595    \frac{1}{\Delta t} w^{n+1} + \partial_r \phi_{nh}^{n+1}
596    & = & \frac{1}{\Delta t} w^{n} + G_w^{(n+1/2)} \label{eq:discrete-time-w-nh} \\
597    \end{eqnarray}
598    which must satisfy the discrete-in-time depth integrated continuity,
599    equation~\ref{eq:discrete-time-backward-free-surface} and the local continuity equation
600    \begin{equation}
601    \partial_x u^{n+1} + \partial_y v^{n+1} + \partial_r w^{n+1} = 0
602    \label{eq:non-divergence-nh}
603    \end{equation}
604    As before, the explicit predictions for momentum are consolidated as:
605    \begin{eqnarray*}
606    u^* & = & u^n + \Delta t G_u^{(n+1/2)} \\
607    v^* & = & v^n + \Delta t G_v^{(n+1/2)} \\
608    w^* & = & w^n + \Delta t G_w^{(n+1/2)}
609    \end{eqnarray*}
610    but this time we introduce an intermediate step by splitting the
611    tendancy of the flow as follows:
612    \begin{eqnarray}
613    u^{n+1} = u^{**} - \Delta t \partial_x \phi_{nh}^{n+1}
614    & &
615    u^{**} = u^{*} - \Delta t g \partial_x \eta^{n+1} \\
616    v^{n+1} = v^{**} - \Delta t \partial_y \phi_{nh}^{n+1}
617    & &
618    v^{**} = v^{*} - \Delta t g \partial_y \eta^{n+1}
619    \end{eqnarray}
620    Substituting into the depth integrated continuity
621    (equation~\ref{eq:discrete-time-backward-free-surface}) gives
622    \begin{equation}
623    \partial_x H \partial_x \left( g \eta^{n+1} + \widehat{\phi}_{nh}^{n+1} \right)
624    +
625    \partial_y H \partial_y \left( g \eta^{n+1} + \widehat{\phi}_{nh}^{n+1} \right)
626     - \frac{\epsilon_{fs}\eta^*}{\Delta t^2}
627    = - \frac{\eta^*}{\Delta t^2}
628    \end{equation}
629    which is approximated by equation
630    \ref{eq:elliptic-backward-free-surface} on the basis that i)
631    $\phi_{nh}^{n+1}$ is not yet known and ii) $\nabla \widehat{\phi}_{nh}
632    << g \nabla \eta$. If \ref{eq:elliptic-backward-free-surface} is
633    solved accurately then the implication is that $\widehat{\phi}_{nh}
634    \approx 0$ so that thet non-hydrostatic pressure field does not drive
635    barotropic motion.
636    
637    The flow must satisfy non-divergence
638    (equation~\ref{eq:non-divergence-nh}) locally, as well as depth
639    integrated, and this constraint is used to form a 3-D elliptic
640    equations for $\phi_{nh}^{n+1}$:
641    \begin{equation}
642    \partial_{xx} \phi_{nh}^{n+1} + \partial_{yy} \phi_{nh}^{n+1} +
643    \partial_{rr} \phi_{nh}^{n+1} =
644    \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*}
645    \end{equation}
646    
647    The entire algorithm can be summarized as the sequential solution of
648    the following equations:
649    \begin{eqnarray}
650    u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2)} \label{eq:ustar-nh} \\
651    v^{*} & = & v^{n} + \Delta t G_v^{(n+1/2)} \label{eq:vstar-nh} \\
652    w^{*} & = & w^{n} + \Delta t G_w^{(n+1/2)} \label{eq:wstar-nh} \\
653    \eta^* & = & \epsilon_{fs} \left( \eta^{n} +P-E+R \right)- \Delta t
654      \partial_x H \widehat{u^{*}}
655    + \partial_y H \widehat{v^{*}}
656    \\
657      \partial_x g H \partial_x \eta^{n+1}
658    + \partial_y g H \partial_y \eta^{n+1}
659    - \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2}
660    & = &
661    - \frac{\eta^*}{\Delta t^2}
662    \label{eq:elliptic-nh}
663    \\
664    u^{**} & = & u^{*} - \Delta t g \partial_x \eta^{n+1} \label{eq:unx-nh}\\
665    v^{**} & = & v^{*} - \Delta t g \partial_y \eta^{n+1} \label{eq:vnx-nh}\\
666    \partial_{xx} \phi_{nh}^{n+1} + \partial_{yy} \phi_{nh}^{n+1} +
667    \partial_{rr} \phi_{nh}^{n+1} & = &
668    \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*} \\
669    u^{n+1} & = & u^{**} - \Delta t \partial_x \phi_{nh}^{n+1} \label{eq:un+1-nh}\\
670    v^{n+1} & = & v^{**} - \Delta t \partial_y \phi_{nh}^{n+1} \label{eq:vn+1-nh}\\
671    \partial_r w^{n+1} & = & - \partial_x u^{n+1} - \partial_y v^{n+1}
672    \end{eqnarray}
673    where the last equation is solved by vertically integrating for
674    $w^{n+1}$.
675    
676    
677    
678  \section{Variants on the Free Surface}  \section{Variants on the Free Surface}
679    
680  We now descibe the various formulations of the free-surface that  We now describe the various formulations of the free-surface that
681  include non-linear forms, implicit in time using Crank-Nicholson,  include non-linear forms, implicit in time using Crank-Nicholson,
682  explicit and [one day] split-explicit. First, we'll reiterate the  explicit and [one day] split-explicit. First, we'll reiterate the
683  underlying algorithm but this time using the notation consistent with  underlying algorithm but this time using the notation consistent with
684  the more general vertical coordinate $r$. The elliptic equation for  the more general vertical coordinate $r$. The elliptic equation for
685  free-surface coordinate (units of $r$), correpsonding to  free-surface coordinate (units of $r$), corresponding to
686  \ref{eq:discrete-time-backward-free-surface}, and  \ref{eq:discrete-time-backward-free-surface}, and
687  assuming no non-hydrostatic effects ($\epsilon_{nh} = 0$) is:  assuming no non-hydrostatic effects ($\epsilon_{nh} = 0$) is:
688  \begin{eqnarray}  \begin{eqnarray}
# Line 611  $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h) Line 714  $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)
714    
715    
716  Once ${\eta}^{n+1}$ has been found, substituting into  Once ${\eta}^{n+1}$ has been found, substituting into
717  \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
718  hydrostatic ($\epsilon_{nh}=0$):  hydrostatic ($\epsilon_{nh}=0$):
719  $$  $$
720  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}
# Line 621  $$ Line 724  $$
724  This is known as the correction step. However, when the model is  This is known as the correction step. However, when the model is
725  non-hydrostatic ($\epsilon_{nh}=1$) we need an additional step and an  non-hydrostatic ($\epsilon_{nh}=1$) we need an additional step and an
726  additional equation for $\phi'_{nh}$. This is obtained by substituting  additional equation for $\phi'_{nh}$. This is obtained by substituting
727  \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}
728  \ref{eq-tDsC-cont}:  into continuity:
729  \begin{equation}  \begin{equation}
730  \left[ {\bf \nabla}_h^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1}  \left[ {\bf \nabla}_h^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1}
731  = \frac{1}{\Delta t} \left(  = \frac{1}{\Delta t} \left(
# Line 702  In the code, $\beta,\gamma$ are defined Line 805  In the code, $\beta,\gamma$ are defined
805  {\it implicSurfPress}, {\it implicDiv2DFlow}. They are read from  {\it implicSurfPress}, {\it implicDiv2DFlow}. They are read from
806  the main data file "{\it data}" and are set by default to 1,1.  the main data file "{\it data}" and are set by default to 1,1.
807    
808  Equations \ref{eq-tDsC-Hmom} and \ref{eq-tDsC-eta} are modified as follows:  Equations \ref{eq:ustar-backward-free-surface} --
809    \ref{eq:vn+1-backward-free-surface} are modified as follows:
810  $$  $$
811  \frac{ \vec{\bf v}^{n+1} }{ \Delta t }  \frac{ \vec{\bf v}^{n+1} }{ \Delta t }
812  + {\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} ]
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.22