/[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.16 by edhill, Thu Aug 7 18:27:52 2003 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 \\  \filelink{FORWARD\_STEP}{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 algorithm
95      (\filelink{FORWARD\_STEP}{model-src-forward_step.F})}
96  \label{fig:call-tree-pressure-method}  \label{fig:call-tree-pressure-method}
97  \end{figure}  \end{figure}
98    
# Line 102  flow boundary conditions applied, become Line 114  flow boundary conditions applied, become
114  \label{eq:rigid-lid-continuity}  \label{eq:rigid-lid-continuity}
115  \end{equation}  \end{equation}
116  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$,
117  similarly for $H\widehat{v}$. The rigid-lid approzimation sets $w=0$  similarly for $H\widehat{v}$. The rigid-lid approximation sets $w=0$
118  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
119  exerted on the fluid by the lid. The horizontal momentum equations and  exerted on the fluid by the lid. The horizontal momentum equations and
120  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 142  pressure gradient is explicit and so can
142  $G$.  $G$.
143    
144  Substituting the two momentum equations into the depth integrated  Substituting the two momentum equations into the depth integrated
145  continuity equation eliminates $u^{n+1}$ and $v^{n+1}$ yeilding an  continuity equation eliminates $u^{n+1}$ and $v^{n+1}$ yielding an
146  elliptic equation for $\eta^{n+1}$. Equations  elliptic equation for $\eta^{n+1}$. Equations
147  \ref{eq:discrete-time-u}, \ref{eq:discrete-time-v} and  \ref{eq:discrete-time-u}, \ref{eq:discrete-time-v} and
148  \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 169  Fig.~\ref{fig:pressure-method-rigid-lid}
169  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
170  have been drawn as staggered in time with the flow.  have been drawn as staggered in time with the flow.
171    
172  The correspondance to the code is as follows:  The correspondence to the code is as follows:
173  \begin{itemize}  \begin{itemize}
174  \item  \item
175  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},
176  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
177  {\em TIMESTEP.F}  \filelink{TIMESTEP()}{model-src-timestep.F}
178  \item  \item
179  the vertical integration, $H \widehat{u^*}$ and $H  the vertical integration, $H \widehat{u^*}$ and $H
180  \widehat{v^*}$, divergence and invertion of the elliptic operator in  \widehat{v^*}$, divergence and inversion of the elliptic operator in
181  equation \ref{eq:elliptic} is coded in {\em  equation \ref{eq:elliptic} is coded in
182  SOLVE\_FOR\_PRESSURE.F}  \filelink{SOLVE\_FOR\_PRESSURE()}{model-src-solve_for_pressure.F}
183  \item  \item
184  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
185  \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
186    \filelink{CORRECTION\_STEP()}{model-src-correction_step.F}.
187  \end{itemize}  \end{itemize}
188  The calling tree for these routines is given in  The calling tree for these routines is given in
189  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 263  to~\ref{eq:vn+1-backward-free-surface},
263  the pressure method algorithm with a backward implicit, linearized  the pressure method algorithm with a backward implicit, linearized
264  free surface. The method is still formerly a pressure method because  free surface. The method is still formerly a pressure method because
265  in the limit of large $\Delta t$ the rigid-lid method is  in the limit of large $\Delta t$ the rigid-lid method is
266  reovered. However, the implicit treatment of the free-surface allows  recovered. However, the implicit treatment of the free-surface allows
267  the flow to be divergent and for the surface pressure/elevation to  the flow to be divergent and for the surface pressure/elevation to
268  respond on a finite time-scale (as opposed to instantly). To recovere  respond on a finite time-scale (as opposed to instantly). To recover
269  the rigid-lid formulation, we introduced a switch-like parameter,  the rigid-lid formulation, we introduced a switch-like parameter,
270  $\epsilon_{fs}$, which selects between the free-surface and rigid-lid;  $\epsilon_{fs}$, which selects between the free-surface and rigid-lid;
271  $\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 344  A stability analysis for a relaxation eq
344    
345  Vertical diffusion and viscosity can be treated implicitly in time  Vertical diffusion and viscosity can be treated implicitly in time
346  using the backward method which is an intrinsic scheme. For tracers,  using the backward method which is an intrinsic scheme. For tracers,
347  the time discrretized equation is:  the time discretized equation is:
348  \begin{equation}  \begin{equation}
349  \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} =
350  \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 390  implicit and are thus cast as a an expli
390  \caption{  \caption{
391  A schematic of the explicit Adams-Bashforth and implicit time-stepping  A schematic of the explicit Adams-Bashforth and implicit time-stepping
392  phases of the algorithm. All prognostic variables are co-located in  phases of the algorithm. All prognostic variables are co-located in
393  time. Explicit tendancies are evaluated at time level $n$ as a  time. Explicit tendencies are evaluated at time level $n$ as a
394  function of the state at that time level (dotted arrow). The explicit  function of the state at that time level (dotted arrow). The explicit
395  tendancy from the previous time level, $n-1$, is used to extrapolate  tendency from the previous time level, $n-1$, is used to extrapolate
396  tendancies to $n+1/2$ (dashed arrow). This extrapolated tendancy  tendencies to $n+1/2$ (dashed arrow). This extrapolated tendency
397  allows variables to be stably integrated forward-in-time to render an  allows variables to be stably integrated forward-in-time to render an
398  estimate ($*$-variables) at the $n+1$ time level (solid  estimate ($*$-variables) at the $n+1$ time level (solid
399  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 430  Adams-Bashforth time-stepping.}
430  \label{fig:call-tree-adams-bashforth-sync}  \label{fig:call-tree-adams-bashforth-sync}
431  \end{figure}  \end{figure}
432    
433  The Adams-Bashforth extrapolation of explicit tendancies fits neatly  The Adams-Bashforth extrapolation of explicit tendencies fits neatly
434  into the pressure method algorithm when all state variables are  into the pressure method algorithm when all state variables are
435  co-locacted in time. Fig.~\ref{fig:adams-bashforth-sync} illustrates  co-located in time. Fig.~\ref{fig:adams-bashforth-sync} illustrates
436  the location of variables in time and the evolution of the algorithm  the location of variables in time and the evolution of the algorithm
437  with time. The algorithm can be represented by the sequential solution  with time. The algorithm can be represented by the sequential solution
438  of the follow equations:  of the follow equations:
# Line 453  G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil Line 466  G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil
466  \end{eqnarray}  \end{eqnarray}
467  Fig.~\ref{fig:adams-bashforth-sync} illustrates the location of  Fig.~\ref{fig:adams-bashforth-sync} illustrates the location of
468  variables in time and evolution of the algorithm with time. The  variables in time and evolution of the algorithm with time. The
469  Adams-Bashforth extrapolation of the tracer tendancies is illustrated  Adams-Bashforth extrapolation of the tracer tendencies is illustrated
470  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
471  solid arc. Inversion of the implicit terms, ${\cal  solid arc. Inversion of the implicit terms, ${\cal
472  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
473  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 493  time-step. The corresponding calling tre
493  \caption{  \caption{
494  A schematic of the explicit Adams-Bashforth and implicit time-stepping  A schematic of the explicit Adams-Bashforth and implicit time-stepping
495  phases of the algorithm but with staggering in time of thermodynamic  phases of the algorithm but with staggering in time of thermodynamic
496  variables with the flow. Explicit thermodynamics tendancies are  variables with the flow. Explicit thermodynamics tendencies are
497  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
498  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
499  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
500  extrapolate tendancies to $n$ (dashed arrow). This extrapolated  extrapolate tendencies to $n$ (dashed arrow). This extrapolated
501  tendancy allows thermo-dynamics variables to be stably integrated  tendency allows thermo-dynamics variables to be stably integrated
502  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$
503  time level (solid arc-arrow). The implicit-in-time operator ${\cal  time level (solid arc-arrow). The implicit-in-time operator ${\cal
504  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 515  limiting process for determining a stabl
515  circumstance, it is more efficient to stagger in time the  circumstance, it is more efficient to stagger in time the
516  thermodynamic variables with the flow  thermodynamic variables with the flow
517  variables. Fig.~\ref{fig:adams-bashforth-staggered} illustrates the  variables. Fig.~\ref{fig:adams-bashforth-staggered} illustrates the
518  staggering and algorith. The key difference between this and  staggering and algorithm. The key difference between this and
519  Fig.~\ref{fig:adams-bashforth-sync} is that the new thermodynamics  Fig.~\ref{fig:adams-bashforth-sync} is that the new thermodynamics
520  fields are used to compute the hydrostatic pressure at time level  fields are used to compute the hydrostatic pressure at time level
521  $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 528  algorithm involves replacing equation \r
528  \phi_{hyd}^n = \int b(\theta^{n+1},S^{n+1}) dr  \phi_{hyd}^n = \int b(\theta^{n+1},S^{n+1}) dr
529  \end{displaymath}  \end{displaymath}
530  but the pressure gradient must also be taken out of the  but the pressure gradient must also be taken out of the
531  Adams-Bashforth extrapoltion. Also, retaining the integer time-levels,  Adams-Bashforth extrapolation. Also, retaining the integer time-levels,
532  $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
533  located in time.  Instead, we re-write the entire algorithm,  located in time.  Instead, we re-write the entire algorithm,
534  \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 567  is activated with the run-time flag {\bf
567  {\em PARM01} of {\em data}.  {\em PARM01} of {\em data}.
568    
569  The only difficulty with this approach is apparent in equation  The only difficulty with this approach is apparent in equation
570  $\ref{eq:Gt-n-staggered}$ and illustrated by the dotted arrow  \ref{eq:Gt-n-staggered} and illustrated by the dotted arrow
571  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
572  tracers around is not naturally located in time. This could be avoided  tracers around is not naturally located in time. This could be avoided
573  by applying the Adams-Bashforth extrapolation to the tracer field  by applying the Adams-Bashforth extrapolation to the tracer field
574  itself and advection that around but this is not yet available. We're  itself and advecting that around but this approach is not yet
575  not aware of any detrimental effect of this feature. The difficulty  available. We're not aware of any detrimental effect of this
576  lies mainly in interpretation of what time-level variables and terms  feature. The difficulty lies mainly in interpretation of what
577  correspond to.  time-level variables and terms correspond to.
578    
579    
580  \section{Non-hydrostatic formulation}  \section{Non-hydrostatic formulation}
581  \label{sect:non-hydrostatic}  \label{sect:non-hydrostatic}
582    
583  [to be written...]  The non-hydrostatic formulation re-introduces the full vertical
584    momentum equation and requires the solution of a 3-D elliptic
585    equations for non-hydrostatic pressure perturbation. We still
586    intergrate vertically for the hydrostatic pressure and solve a 2-D
587    elliptic equation for the surface pressure/elevation for this reduces
588    the amount of work needed to solve for the non-hydrostatic pressure.
589    
590    The momentum equations are discretized in time as follows:
591    \begin{eqnarray}
592    \frac{1}{\Delta t} u^{n+1} + g \partial_x \eta^{n+1} + \partial_x \phi_{nh}^{n+1}
593    & = & \frac{1}{\Delta t} u^{n} + G_u^{(n+1/2)} \label{eq:discrete-time-u-nh} \\
594    \frac{1}{\Delta t} v^{n+1} + g \partial_y \eta^{n+1} + \partial_y \phi_{nh}^{n+1}
595    & = & \frac{1}{\Delta t} v^{n} + G_v^{(n+1/2)} \label{eq:discrete-time-v-nh} \\
596    \frac{1}{\Delta t} w^{n+1} + \partial_r \phi_{nh}^{n+1}
597    & = & \frac{1}{\Delta t} w^{n} + G_w^{(n+1/2)} \label{eq:discrete-time-w-nh} \\
598    \end{eqnarray}
599    which must satisfy the discrete-in-time depth integrated continuity,
600    equation~\ref{eq:discrete-time-backward-free-surface} and the local continuity equation
601    \begin{equation}
602    \partial_x u^{n+1} + \partial_y v^{n+1} + \partial_r w^{n+1} = 0
603    \label{eq:non-divergence-nh}
604    \end{equation}
605    As before, the explicit predictions for momentum are consolidated as:
606    \begin{eqnarray*}
607    u^* & = & u^n + \Delta t G_u^{(n+1/2)} \\
608    v^* & = & v^n + \Delta t G_v^{(n+1/2)} \\
609    w^* & = & w^n + \Delta t G_w^{(n+1/2)}
610    \end{eqnarray*}
611    but this time we introduce an intermediate step by splitting the
612    tendancy of the flow as follows:
613    \begin{eqnarray}
614    u^{n+1} = u^{**} - \Delta t \partial_x \phi_{nh}^{n+1}
615    & &
616    u^{**} = u^{*} - \Delta t g \partial_x \eta^{n+1} \\
617    v^{n+1} = v^{**} - \Delta t \partial_y \phi_{nh}^{n+1}
618    & &
619    v^{**} = v^{*} - \Delta t g \partial_y \eta^{n+1}
620    \end{eqnarray}
621    Substituting into the depth integrated continuity
622    (equation~\ref{eq:discrete-time-backward-free-surface}) gives
623    \begin{equation}
624    \partial_x H \partial_x \left( g \eta^{n+1} + \widehat{\phi}_{nh}^{n+1} \right)
625    +
626    \partial_y H \partial_y \left( g \eta^{n+1} + \widehat{\phi}_{nh}^{n+1} \right)
627     - \frac{\epsilon_{fs}\eta^*}{\Delta t^2}
628    = - \frac{\eta^*}{\Delta t^2}
629    \end{equation}
630    which is approximated by equation
631    \ref{eq:elliptic-backward-free-surface} on the basis that i)
632    $\phi_{nh}^{n+1}$ is not yet known and ii) $\nabla \widehat{\phi}_{nh}
633    << g \nabla \eta$. If \ref{eq:elliptic-backward-free-surface} is
634    solved accurately then the implication is that $\widehat{\phi}_{nh}
635    \approx 0$ so that thet non-hydrostatic pressure field does not drive
636    barotropic motion.
637    
638    The flow must satisfy non-divergence
639    (equation~\ref{eq:non-divergence-nh}) locally, as well as depth
640    integrated, and this constraint is used to form a 3-D elliptic
641    equations for $\phi_{nh}^{n+1}$:
642    \begin{equation}
643    \partial_{xx} \phi_{nh}^{n+1} + \partial_{yy} \phi_{nh}^{n+1} +
644    \partial_{rr} \phi_{nh}^{n+1} =
645    \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*}
646    \end{equation}
647    
648    The entire algorithm can be summarized as the sequential solution of
649    the following equations:
650    \begin{eqnarray}
651    u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2)} \label{eq:ustar-nh} \\
652    v^{*} & = & v^{n} + \Delta t G_v^{(n+1/2)} \label{eq:vstar-nh} \\
653    w^{*} & = & w^{n} + \Delta t G_w^{(n+1/2)} \label{eq:wstar-nh} \\
654    \eta^* & = & \epsilon_{fs} \left( \eta^{n} +P-E+R \right)- \Delta t
655      \partial_x H \widehat{u^{*}}
656    + \partial_y H \widehat{v^{*}}
657    \\
658      \partial_x g H \partial_x \eta^{n+1}
659    + \partial_y g H \partial_y \eta^{n+1}
660    - \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2}
661    & = &
662    - \frac{\eta^*}{\Delta t^2}
663    \label{eq:elliptic-nh}
664    \\
665    u^{**} & = & u^{*} - \Delta t g \partial_x \eta^{n+1} \label{eq:unx-nh}\\
666    v^{**} & = & v^{*} - \Delta t g \partial_y \eta^{n+1} \label{eq:vnx-nh}\\
667    \partial_{xx} \phi_{nh}^{n+1} + \partial_{yy} \phi_{nh}^{n+1} +
668    \partial_{rr} \phi_{nh}^{n+1} & = &
669    \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*} \\
670    u^{n+1} & = & u^{**} - \Delta t \partial_x \phi_{nh}^{n+1} \label{eq:un+1-nh}\\
671    v^{n+1} & = & v^{**} - \Delta t \partial_y \phi_{nh}^{n+1} \label{eq:vn+1-nh}\\
672    \partial_r w^{n+1} & = & - \partial_x u^{n+1} - \partial_y v^{n+1}
673    \end{eqnarray}
674    where the last equation is solved by vertically integrating for
675    $w^{n+1}$.
676    
677    
678    
679  \section{Variants on the Free Surface}  \section{Variants on the Free Surface}
680    
681  We now descibe the various formulations of the free-surface that  We now describe the various formulations of the free-surface that
682  include non-linear forms, implicit in time using Crank-Nicholson,  include non-linear forms, implicit in time using Crank-Nicholson,
683  explicit and [one day] split-explicit. First, we'll reiterate the  explicit and [one day] split-explicit. First, we'll reiterate the
684  underlying algorithm but this time using the notation consistent with  underlying algorithm but this time using the notation consistent with
685  the more general vertical coordinate $r$. The elliptic equation for  the more general vertical coordinate $r$. The elliptic equation for
686  free-surface coordinate (units of $r$), correpsonding to  free-surface coordinate (units of $r$), corresponding to
687  \ref{eq:discrete-time-backward-free-surface}, and  \ref{eq:discrete-time-backward-free-surface}, and
688  assuming no non-hydrostatic effects ($\epsilon_{nh} = 0$) is:  assuming no non-hydrostatic effects ($\epsilon_{nh} = 0$) is:
689  \begin{eqnarray}  \begin{eqnarray}
# Line 611  $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h) Line 715  $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)
715    
716    
717  Once ${\eta}^{n+1}$ has been found, substituting into  Once ${\eta}^{n+1}$ has been found, substituting into
718  \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
719  hydrostatic ($\epsilon_{nh}=0$):  hydrostatic ($\epsilon_{nh}=0$):
720  $$  $$
721  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}
# Line 621  $$ Line 725  $$
725  This is known as the correction step. However, when the model is  This is known as the correction step. However, when the model is
726  non-hydrostatic ($\epsilon_{nh}=1$) we need an additional step and an  non-hydrostatic ($\epsilon_{nh}=1$) we need an additional step and an
727  additional equation for $\phi'_{nh}$. This is obtained by substituting  additional equation for $\phi'_{nh}$. This is obtained by substituting
728  \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}
729  \ref{eq-tDsC-cont}:  into continuity:
730  \begin{equation}  \begin{equation}
731  \left[ {\bf \nabla}_h^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1}  \left[ {\bf \nabla}_h^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1}
732  = \frac{1}{\Delta t} \left(  = \frac{1}{\Delta t} \left(
# Line 702  In the code, $\beta,\gamma$ are defined Line 806  In the code, $\beta,\gamma$ are defined
806  {\it implicSurfPress}, {\it implicDiv2DFlow}. They are read from  {\it implicSurfPress}, {\it implicDiv2DFlow}. They are read from
807  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.
808    
809  Equations \ref{eq-tDsC-Hmom} and \ref{eq-tDsC-eta} are modified as follows:  Equations \ref{eq:ustar-backward-free-surface} --
810    \ref{eq:vn+1-backward-free-surface} are modified as follows:
811  $$  $$
812  \frac{ \vec{\bf v}^{n+1} }{ \Delta t }  \frac{ \vec{\bf v}^{n+1} }{ \Delta t }
813  + {\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.16
  ViewVC Help
Powered by ViewVC 1.1.22