/[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.27 by jmc,
Fri Aug 27 13:08:18 2010 UTC
+revision 1.34 by jmc,
Mon Sep 25 18:11:06 2017 UTC
 Line 13 
 passive and dynamically active tracers a
  \input{s_algorithm/text/notation}
  \section{Time-stepping}
+ \label{sec:time_stepping}
  \begin{rawhtml}
  <!-- CMIREDIR:time-stepping: -->
  \end{rawhtml}
-Line 50 
 evaluated explicitly in time. Since the
+Line 51 
 evaluated explicitly in time. Since the
  independent of the particular time-stepping scheme chosen we will
  describe first the over-arching algorithm, known as the pressure
  method, with a rigid-lid model in section
- \ref{sect:pressure-method-rigid-lid}. This algorithm is essentially
+ \ref{sec: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 flavors
+ section \ref{sec: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
+ in section \ref{sec:adams-bashforth} and put into the context of the
- overall algorithm in sections \ref{sect:adams-bashforth-sync} and
+ overall algorithm in sections \ref{sec:adams-bashforth-sync} and
- \ref{sect:adams-bashforth-staggered}. Inclusion of non-hydrostatic
+ \ref{sec:adams-bashforth-staggered}. Inclusion of non-hydrostatic
  terms requires applying the pressure method in three dimensions
  instead of two and this algorithm modification is described in section
- \ref{sect:non-hydrostatic}. Finally, the free-surface equation may be
+ \ref{sec:non-hydrostatic}. Finally, the free-surface equation may be
  treated more exactly, including non-linear terms, and this is
- described in section \ref{sect:nonlinear-freesurface}.
+ described in section \ref{sec:nonlinear-freesurface}.
  \section{Pressure method with rigid-lid}
- \label{sect:pressure-method-rigid-lid}
+ \label{sec:pressure-method-rigid-lid}
  \begin{rawhtml}
  <!-- CMIREDIR:pressure_method_rigid_lid: -->
  \end{rawhtml}
-Line 89 
 temporary.}
+Line 90 
 temporary.}
  \begin{figure}
  \begin{center} \fbox{ \begin{minipage}{4.5in} \begin{tabbing}
  aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill
  \filelink{FORWARD\_STEP}{model-src-forward_step.F} \\
  \> DYNAMICS \\
  \>\> TIMESTEP \` $u^*$,$v^*$ (\ref{eq:ustar-rigid-lid},\ref{eq:vstar-rigid-lid}) \\
-Line 100 
 aaa \= aaa \= aaa \= aaa \= aaa \= aaa \
+Line 101 
 aaa \= aaa \= aaa \= aaa \= aaa \= aaa \
  \>\> 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 algorithm
    (\filelink{FORWARD\_STEP}{model-src-forward_step.F})}
  \label{fig:call-tree-pressure-method}
  \end{figure}
-Line 187 
 stepping forward $u^n$ and $v^n$ to $u^{
+Line 188 
 stepping forward $u^n$ and $v^n$ to $u^{
  \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
  \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
  \filelink{CORRECTION\_STEP()}{model-src-correction_step.F}.
  \end{itemize}
  The calling tree for these routines is given in
-Line 199 
 Fig.~\ref{fig:call-tree-pressure-method}
+Line 200 
 Fig.~\ref{fig:call-tree-pressure-method}
  %\paragraph{Need to discuss implicit viscosity somewhere:}
  In general, the horizontal momentum time-stepping can contain some terms
  that are treated implicitly in time,
  such as the vertical viscosity when using the backward time-stepping scheme
  (\varlink{implicitViscosity}{implicitViscosity} {\it =.TRUE.}).
  The method used to solve those implicit terms is provided in
- section \ref{sect:implicit-backward-stepping}, and modifies
+ section \ref{sec:implicit-backward-stepping}, and modifies
  equations \ref{eq:discrete-time-u} and \ref{eq:discrete-time-v} to
  give:
  \begin{eqnarray}
-Line 217 
 v^{n+1} - \Delta t \partial_z A_v \parti
+Line 218 
 v^{n+1} - \Delta t \partial_z A_v \parti
  \section{Pressure method with implicit linear free-surface}
- \label{sect:pressure-method-linear-backward}
+ \label{sec:pressure-method-linear-backward}
  \begin{rawhtml}
  <!-- CMIREDIR:pressure_method_linear_backward: -->
  \end{rawhtml}
-Line 253 
 pressure method is then replaced by the
+Line 254 
 pressure method is then replaced by the
  where the use of flow at time level $n+1$ makes the method implicit
  and backward in time. This is the preferred scheme since it still
  filters the fast unresolved wave motions by damping them. A centered
- scheme, such as Crank-Nicholson (see section \ref{sect:freesurf-CrankNick}),
+ scheme, such as Crank-Nicholson (see section \ref{sec:freesurf-CrankNick}),
  would alias the energy of the fast modes onto slower modes of motion.
  As for the rigid-lid pressure method, equations
-Line 262 
 As for the rigid-lid pressure method, eq
+Line 263 
 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} + \Delta t (P-E) \right)- \Delta t
+ \eta^* & = & \epsilon_{fs} \left( \eta^{n} + \Delta t (P-E) \right)
-   \partial_x H \widehat{u^{*}}
+          - \Delta t \left( \partial_x H \widehat{u^{*}}
- + \partial_y H \widehat{v^{*}}
+                          + \partial_y H \widehat{v^{*}} \right)
  \\
    \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-backward-free-surface}
  \\
-Line 285 
 recovered. However, the implicit treatme
+Line 286 
 recovered. However, the implicit treatme
  the flow to be divergent and for the surface pressure/elevation to
  respond on a finite time-scale (as opposed to instantly). To recover
  the rigid-lid formulation, we introduced a switch-like parameter,
  $\epsilon_{fs}$ (\varlink{freesurfFac}{freesurfFac}),
  which selects between the free-surface and rigid-lid;
  $\epsilon_{fs}=1$ allows the free-surface to evolve; $\epsilon_{fs}=0$
  imposes the rigid-lid. The evolution in time and location of variables
-Line 297 
 pressure-method.
+Line 298 
 pressure-method.
  \section{Explicit time-stepping: Adams-Bashforth}
- \label{sect:adams-bashforth}
+ \label{sec:adams-bashforth}
  \begin{rawhtml}
  <!-- CMIREDIR:adams_bashforth: -->
  \end{rawhtml}
-Line 307 
 time discretization of the explicit term
+Line 308 
 time discretization of the explicit term
  the quasi-second order Adams-Bashforth method for all explicit terms
  in both the momentum and tracer equations. This is still the default
  mode of operation but it is now possible to use alternate schemes for
- tracers (see section \ref{sect:tracer-advection}).
+ tracers (see section \ref{sec:tracer-advection}).
  \begin{figure}
  \begin{center} \fbox{ \begin{minipage}{4.5in} \begin{tabbing}
  aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill
  FORWARD\_STEP \\
  \> THERMODYNAMICS \\
  \>\> CALC\_GT \\
-Line 357 
 simpler to include these terms and this
+Line 358 
 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 tracer and momentum
  forcing terms and the dissipation terms outside the
  Adams-Bashforth extrapolation, by turning off the logical flags
  \varlink{forcing\_In\_AB}{forcing_In_AB}
  (parameter file {\em data}, namelist {\em PARM01}, default value = True).
  (\varlink{tracForcingOutAB}{tracForcingOutAB}, default=0,
-Line 396 
 The right column represents the damping
+Line 397 
 The right column represents the damping
  \section{Implicit time-stepping: backward method}
- \label{sect:implicit-backward-stepping}
+ \label{sec:implicit-backward-stepping}
  \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.
  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:
-Line 424 
 using the Adams-Bashforth method as desc
+Line 425 
 using the Adams-Bashforth method as desc
  \end{eqnarray}
  where ${\cal L}_\tau^{-1}$ is the inverse of the operator
  \begin{equation}
- {\cal L} = \left[ 1 + \Delta t \partial_r \kappa_v \partial_r \right]
+ {\cal L}_\tau = \left[ 1 + \Delta t \partial_r \kappa_v \partial_r \right]
  \end{equation}
  Equation \ref{eq:taustar-implicit} looks exactly as \ref{eq:taustar}
  while \ref{eq:tau-n+1-implicit} involves an operator or matrix
-Line 444 
 unconditionally stable, no-slip boundary
+Line 445 
 unconditionally stable, no-slip boundary
  implicit and are thus cast as a an explicit drag term.
  \section{Synchronous time-stepping: variables co-located in time}
- \label{sect:adams-bashforth-sync}
+ \label{sec:adams-bashforth-sync}
  \begin{rawhtml}
  <!-- CMIREDIR:adams_bashforth_sync: -->
  \end{rawhtml}
-Line 469 
 is solved to yield the state variables a
+Line 470 
 is solved to yield the state variables a
  \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 \\
-Line 500 
 FORWARD\_STEP \\
+Line 501 
 FORWARD\_STEP \\
  \end{tabbing} \end{minipage} } \end{center}
  \caption{
  Calling tree for the overall synchronous 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.
  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 535 
 G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil
+Line 536 
 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} + \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}
  ~ = ~ - \frac{\eta^*}{\Delta t^2}
  \label{eq:elliptic-sync} \\
- \vec{\bf v}^{n+1} & = & \vec{\bf v}^{*} - \Delta t g \nabla \eta^{n+1}
+ \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
-Line 564 
 time-step. The corresponding calling tre
+Line 565 
 time-step. The corresponding calling tre
  \ref{fig:call-tree-adams-bashforth-sync}.
  \section{Staggered baroclinic time-stepping}
- \label{sect:adams-bashforth-staggered}
+ \label{sec:adams-bashforth-staggered}
  \begin{rawhtml}
  <!-- CMIREDIR:adams_bashforth_staggered: -->
  \end{rawhtml}
-Line 576 
 time-step. The corresponding calling tre
+Line 577 
 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 momentum tendencies are evaluated at time level $n-1/2$ as a
  function of the flow field at that time level $n-1/2$.
  The explicit tendency from the previous time level, $n-3/2$, is used to
  extrapolate tendencies to $n$ (dashed arrow).
  The hydrostatic pressure/geo-potential $\phi_{hyd}$ is evaluated directly
  at time level $n$ (vertical arrows) and used with the extrapolated tendencies
  to step forward the flow variables from $n-1/2$ to $n+1/2$ (solid arc-arrow).
  The implicit-in-time operator ${\cal L}_{\bf u,v}$ (vertical arrows) is
  then applied to the previous estimation of the the flow field ($*$-variables)
  and yields to the two velocity components $u,v$ at time level $n+1/2$.
  These are then used to calculate the advection term (dashed arc-arrow)
  of the thermo-dynamics tendencies at time step $n$.
  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}
-Line 602 
 circumstance, it is more efficient to st
+Line 603 
 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 thermodynamic variables
  are solved after the dynamics, using the recently updated flow field.
  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 that the
  thermodynamics solver is delayed from half a time step,
  allowing the use of the most recent velocities to compute
  the advection terms. Once the thermodynamics fields are
  updated, the hydrostatic pressure is computed
- to step forwrad the dynamics.
+ to step forward the dynamics.
  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
-Line 630 
 position in time of variables appropriat
+Line 631 
 position in time of variables appropriat
  \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-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/2} & - & \frac{\epsilon_{fs} \eta^{n+1/2}}{\Delta t^2}
  ~ = ~ - \frac{\eta^*}{\Delta t^2}
  \label{eq:elliptic-staggered} \\
- \vec{\bf v}^{n+1/2} & = & \vec{\bf v}^{*} - \Delta t g \nabla \eta^{n+1/2}
+ \vec{\bf v}^{n+1/2} & = & \vec{\bf v}^{**} - \Delta t g \nabla \eta^{n+1/2}
  \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} \\
-Line 649 
 G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil
+Line 650 
 G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil
  \end{eqnarray}
  The corresponding calling tree is given in
  \ref{fig:call-tree-adams-bashforth-staggered}.
  The staggered algorithm is activated with the run-time flag
  {\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)\\
-Line 676 
 FORWARD\_STEP \\
+Line 677 
 FORWARD\_STEP \\
  \>\> 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}) \\
-Line 689 
 FORWARD\_STEP \\
+Line 690 
 FORWARD\_STEP \\
  \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}
-Line 709 
 time-level variables and terms correspon
+Line 710 
 time-level variables and terms correspon
  \section{Non-hydrostatic formulation}
- \label{sect:non-hydrostatic}
+ \label{sec:non-hydrostatic}
  \begin{rawhtml}
  <!-- CMIREDIR:non-hydrostatic_formulation: -->
  \end{rawhtml}
-Line 717 
 time-level variables and terms correspon
+Line 718 
 time-level variables and terms correspon
  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
+ integrate 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.
-Line 758 
 Substituting into the depth integrated c
+Line 759 
 Substituting into the depth integrated c
  \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{\epsilon_{fs}\eta^{n+1}}{\Delta t^2}
  = - \frac{\eta^*}{\Delta t^2}
  \end{equation}
  which is approximated by equation
-Line 766 
 which is approximated by equation
+Line 767 
 which is approximated by equation
  $\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
+ \approx 0$ so that the non-hydrostatic pressure field does not drive
  barotropic motion.
  The flow must satisfy non-divergence
-Line 774 
 The flow must satisfy non-divergence
+Line 775 
 The flow must satisfy non-divergence
  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_{rr} \phi_{nh}^{n+1} = \left(
  \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*}
+ \right) / \Delta t
  \end{equation}
  The entire algorithm can be summarized as the sequential solution of
-Line 786 
 u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2
+Line 788 
 u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2
  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
+ & - & \Delta t \left( \partial_x H \widehat{u^{*}}
-   \partial_x H \widehat{u^{*}}
+                     + \partial_y H \widehat{v^{*}} \right)
- + \partial_y H \widehat{v^{*}}
  \\
    \partial_x g H \partial_x \eta^{n+1}
  + \partial_y g H \partial_y \eta^{n+1}
-Line 799 
 w^{*} & = & w^{n} + \Delta t G_w^{(n+1/2
+Line 800 
 w^{*} & = & w^{n} + \Delta t G_w^{(n+1/2
  \\
  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_{rr} \phi_{nh}^{n+1} & = & \left(
- \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*} \\
+ \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*}
+ \right) / \Delta t  \label{eq:phi-nh}\\
  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}
+ \label{sec:free-surface}
  We now describe the various formulations of the free-surface that
  include non-linear forms, implicit in time using Crank-Nicholson,
 Line 832 
 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}
  \label{eq-solve2D_rhs}
  \end{eqnarray}
 Line 851 
 $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)
  Once ${\eta}^{n+1}$ has been found, substituting into
  \ref{eq:discrete-time-u}, \ref{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 874 
 where
  \end{displaymath}
  Note that $\eta^{n+1}$ is also used to update the second RHS term
  $\partial_r \dot{r}^* $ since
  the vertical velocity at the surface ($\dot{r}_{surf}$)
  is evaluated as $(\eta^{n+1} - \eta^n) / \Delta t$.
  Finally, the horizontal velocities at the new time level are found by:
 Line 924 
 at the same point in the code.
- \subsection{Crank-Nickelson barotropic time stepping}
+ \subsection{Crank-Nicolson barotropic time stepping}
- \label{sect:freesurf-CrankNick}
+ \label{sec:freesurf-CrankNick}
  The full implicit time stepping described previously is
  unconditionally stable but damps the fast gravity waves, resulting in
 Line 936 
 for the barotropic flow divergence ($\ga
  \\
  For instance, $\beta=\gamma=1$ is the previous fully implicit scheme;
  $\beta=\gamma=1/2$ is the non damping (energy conserving), unconditionally
- stable, Crank-Nickelson scheme; $(\beta,\gamma)=(1,0)$ or $=(0,1)$
+ stable, Crank-Nicolson scheme; $(\beta,\gamma)=(1,0)$ or $=(0,1)$
  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
  {\bf implicSurfPress}, {\bf implicDiv2DFlow}. They are read from
- the main parameter file "{\em data}" and are set by default to 1,1.
+ the main parameter file "{\em data}" (namelist {\em PARM01})
+ and are set by default to 1,1.
  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 }
+  = \frac{ \vec{\bf v}^{n} }{ \Delta t }
+  + \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
+  + {\bf \nabla}_h {\phi'_{hyd}}^{(n+1/2)}
  \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:
+ We set
  \begin{eqnarray*}
  \vec{\bf v}^* & = &
  \vec{\bf v} ^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
-Line 966 
 where:
+Line 969 
 where:
  + \Delta t {\bf \nabla}_h {\phi'_{hyd}}^{(n+1/2)}
  \\
  {\eta}^* & = &
  \epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta t (P-E)
  - \Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o}
  [ \gamma \vec{\bf v}^* + (1-\gamma) \vec{\bf v}^{n}] dr
  \end{eqnarray*}
  \\
-Line 978 
 $$
+Line 981 
 $$
  {\bf \nabla}_h \cdot \beta\gamma \Delta t^2 b_s (R_o - R_{fixed})
  {\bf \nabla}_h {\eta}^{n+1}
  = {\eta}^*
  $$
  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.
  \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:
+ affected by the Crank-Nicolson 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
+ %\item The non-hydrostatic part of the code has not yet been
- updated, and therefore 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
+ \item The stability criteria with Crank-Nicolson time stepping
  for the pure linear gravity wave problem in cartesian coordinates is:
  \begin{itemize}
  \item $\beta + \gamma < 1$ : unstable
  \item $\beta \geq 1/2$ and $ \gamma \geq 1/2$ : stable
  \item $\beta + \gamma \geq 1$ : stable if
  $$
  c_{max}^2 (\beta - 1/2)(\gamma - 1/2) + 1 \geq 0
  $$
  $$
  \mbox{with }~
  %c^2 = 2 g H {\Delta t}^2
  %(\frac{1-cos 2 \pi / k}{\Delta x^2}
  %+\frac{1-cos 2 \pi / l}{\Delta y^2})
  %$$
  %Practically, the most stringent condition is obtained with $k=l=2$ :
  %$$
  c_{max} =  2 \Delta t \: \sqrt{g H} \:
  \sqrt{ \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} }
  $$
  \end{itemize}
+ \item A similar mixed forward/backward time-stepping is also available for
+ the non-hydrostatic algorithm,
+ with a fraction $\beta_{nh}$ ($ 0 < \beta_{nh} \leq 1$)
+ of the non-hydrostatic pressure gradient being evaluated at time step $n+1$
+ (backward in time) and the remaining part ($1 - \beta_{nh}$) being evaluated
+ at time step $n$ (forward in time).
+ The run-time parameter {\bf implicitNHPress} corresponding to the implicit
+ fraction $\beta_{nh}$ of the non-hydrostatic pressure is set by default to
+ the implicit fraction $\beta$ of surface pressure ({\bf implicSurfPress}),
+ but can also be specified independently (in main parameter file {\em data},
+ namelist {\em PARM01}).
  \end{enumerate}

 Legend:



Removed from v.1.27
 


changed lines


 
Added in v.1.34
 Legend:



Removed from v.1.27
 


changed lines


 
Added in v.1.34
-Removed from v.1.27
+Added in v.1.34

	ViewVC Help
Powered by ViewVC 1.1.22