/[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.1.1.1 by adcroft,
Wed Aug  8 16:15:16 2001 UTC
+revision 1.7 by adcroft,
Fri Sep 28 14:09:56 2001 UTC
 Line 2
  % $Name$
  The convention used in this section is as follows:
- Time is "discretize" using a time step $\Delta t$
+ Time is ``discretized'' using a time step $\Delta t$
  and $\Phi^n$ refers to the variable $\Phi$
- at time $t = n \Delta t$ . We used the notation $\Phi^{(n)}$
+ at time $t = n \Delta t$ . We use the notation $\Phi^{(n)}$
  when time interpolation is required to estimate the value of $\phi$
  at the time $n \Delta t$.
- \section{Time Integration}
+ \section{Time integration}
  The discretization in time of the model equations (cf section I )
  does not depend of the discretization in space of each
- term, so that this section can be read independently.
+ term and so  this section can be read independently.
- Also for this purpose, we will refers to the continuous
- space-derivative form of model equations, and focus on
- the time discretization.
  The continuous form of the model equations is:
  \begin{eqnarray}
  \partial_t \theta & = & G_\theta
+ \label{eq-tCsC-theta}
  \\
  \partial_t S & = & G_s
+ \label{eq-tCsC-salt}
  \\
  b' & = & b'(\theta,S,r)
  \\
  \partial_r \phi'_{hyd} & = & -b'
- \label{eq-r-split-hyd}
+ \label{eq-tCsC-hyd}
  \\
  \partial_t \vec{\bf v}
- + {\bf \nabla}_r B_o \eta
+ + {\bf \nabla}_h b_s \eta
- + \epsilon_{nh} {\bf \nabla}_r \phi'_{nh}
+ + \epsilon_{nh} {\bf \nabla}_h \phi'_{nh}
  & = & \vec{\bf G}_{\vec{\bf v}}
- - {\bf \nabla}_r \phi'_{hyd}
+ - {\bf \nabla}_h \phi'_{hyd}
- \label{eq-r-split-hmom}
+ \label{eq-tCsC-Hmom}
  \\
  \epsilon_{nh} \frac {\partial{\dot{r}}}{\partial{t}}
  + \epsilon_{nh} \partial_r \phi'_{nh}
  & = & \epsilon_{nh} G_{\dot{r}}
- \label{eq-r-split-rdot}
+ \label{eq-tCsC-Vmom}
  \\
- {\bf \nabla}_r \cdot \vec{\bf v} + \partial_r \dot{r}
+ {\bf \nabla}_h \cdot \vec{\bf v} + \partial_r \dot{r}
  & = & 0
- \label{eq-r-cont}
+ \label{eq-tCsC-cont}
  \end{eqnarray}
  where
  \begin{eqnarray*}
  G_\theta & = &
- - \vec{\bf v} \cdot {\bf \nabla}_r \theta + {\cal Q}_\theta
+ - \vec{\bf v} \cdot {\bf \nabla} \theta + {\cal Q}_\theta
  \\
  G_S & = &
- - \vec{\bf v} \cdot {\bf \nabla}_r S + {\cal Q}_S
+ - \vec{\bf v} \cdot {\bf \nabla} S + {\cal Q}_S
  \\
  \vec{\bf G}_{\vec{\bf v}}
  & = &
- - \vec{\bf v} \cdot {\bf \nabla}_r \vec{\bf v}
+ - \vec{\bf v} \cdot {\bf \nabla} \vec{\bf v}
  - f \hat{\bf k} \wedge \vec{\bf v}
  + \vec{\cal F}_{\vec{\bf v}}
  \\
  G_{\dot{r}}
  & = &
- - \vec{\bf v} \cdot {\bf \nabla}_r \dot{r}
+ - \vec{\bf v} \cdot {\bf \nabla} \dot{r}
  + {\cal F}_{\dot{r}}
  \end{eqnarray*}
- The exact form of all the "{\it G}"s terms is described in the next
+ The exact form of all the ``{\it G}''s terms is described in the next
- section (). Here its sufficient to mention that they contains
+ section \ref{sect:discrete}. Here its sufficient to mention that they contains
- all the RHS terms except the pressure / geo- potential terms.
+ all the RHS terms except the pressure/geo-potential terms.
- The switch $\epsilon_{nh}$ allows to activate the non hydrostatic
+ The switch $\epsilon_{nh}$ allows one to activate the non-hydrostatic
- mode ($\epsilon_{nh}=1$) for the ocean model. Otherwise,
+ mode ($\epsilon_{nh}=1$) for the ocean model. Otherwise, in the
- in the hydrostatic limit $\epsilon_{nh} = 0$ and the 3rd equation vanishes.
+ hydrostatic limit $\epsilon_{nh} = 0$ and equation \ref{eq-tCsC-Vmom}
+ is not used.
- The equation for $\eta$ is obtained by integrating the
- continuity equation over the entire depth of the fluid,
+ As discussed in section \ref{sect:1.3.6.2}, the equation for $\eta$ is
- from $R_{min}(x,y)$ up to $R_o(x,y)$
+ obtained by integrating the continuity equation over the entire depth
- (Linear free surface):
+ of the fluid, from $R_{fixed}(x,y)$ up to $R_o(x,y)$. The linear free
- \begin{displaymath}
+ surface evolves according to:
+ \begin{eqnarray}
  \epsilon_{fs} \partial_t \eta =
  \left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) =
- - {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o} \vec{\bf v} dr
+ - {\bf \nabla} \cdot \int_{R_{fixed}}^{R_o} \vec{\bf v} dr
  + \epsilon_{fw} (P-E)
- \end{displaymath}
+ \label{eq-tCsC-eta}
+ \end{eqnarray}
- Where $\epsilon_{fs}$,$\epsilon_{fw}$ are two flags to
+ Here, $\epsilon_{fs}$ is a flag to switch between the free-surface,
- distinguish between a free-surface equation ($\epsilon_{fs}=1$)
+ $\epsilon_{fs}=1$, and a rigid-lid, $\epsilon_{fs}=0$. The flag
- or the rigid-lid approximation ($\epsilon_{fs}=0$);
+ $\epsilon_{fw}$ determines whether an exchange of fresh water is
- and to distinguish when exchange of Fresh-Water is included
+ included at the ocean surface (natural BC) ($\epsilon_{fw} = 1$) or
- at the ocean surface (natural BC) ($\epsilon_{fw} = 1$)
+ not ($\epsilon_{fw} = 0$).
- or not ($\epsilon_{fw} = 0$).
- The hydrostatic potential is found by
+ The hydrostatic potential is found by integrating (equation
- integrating \ref{eq-r-split-hyd} with the boundary condition that
+ \ref{eq-tCsC-hyd}) with the boundary condition that
  $\phi'_{hyd}(r=R_o) = 0$:
  \begin{eqnarray*}
  & &
 Line 103 
 $\phi'_{hyd}(r=R_o) = 0$:
  \phi'_{hyd}(x,y,r') = \int_{r'}^{R_o} b' dr
  \end{eqnarray*}
- \subsection{standard synchronous time stepping}
+ \subsection{General method}
+ An overview of the general method is now presented with explicit
+ references to the Fortran code. We often refer to the discretized
+ equations of the model that are detailed in the following sections.
+ The general algorithm consist of a ``predictor step'' that computes
+ the forward tendencies ("G" terms") comprising of explicit-in-time
+ terms and the ``first guess'' values (star notation): $\theta^*, S^*,
+ \vec{\bf v}^*$ (and $\dot{r}^*$ in non-hydrostatic mode). This is done
+ in the two routines {\it THERMODYNAMICS} and {\it DYNAMICS}.
+ Terms that are integrated implicitly in time are handled at the end of
+ the {\it THERMODYNAMICS} and {\it DYNAMICS} routines. Then the
+ surface pressure and non hydrostatic pressure are solved for in ({\it
+ SOLVE\_FOR\_PRESSURE}).
+ Finally, in the ``corrector step'', (routine {\it
+ THE\_CORRECTION\_STEP}) the new values of $u,v,w,\theta,S$ are
+ determined (see details in \ref{sect:II.1.3}).
+ At this point, the regular time stepping process is complete. However,
+ after the correction step there are optional adjustments such as
+ convective adjustment or filters (zonal FFT filter, shapiro filter)
+ that can be applied on both momentum and tracer fields, just prior to
+ incrementing the time level to $n+1$.
+ Since the pressure solver precision is of the order of the ``target
+ residual'' and can be lower than the the computer truncation error,
+ and also because some filters might alter the divergence part of the
+ flow field, a final evaluation of the surface r anomaly $\eta^{n+1}$
+ is performed in {\it CALC\_EXACT\_ETA}. This ensures exact volume
+ conservation. Note that there is no need for an equivalent
+ non-hydrostatic ``exact conservation'' step, since by default $w$ is
+ already computed after the filters are applied.
+ Optional forcing terms (usually part of a physics ``package''), that
+ account for a specific source or sink process (e.g. condensation as a
+ sink of water vapor Q) are generally incorporated in the main
+ algorithm as follows: at the the beginning of the time step, the
+ additional tendencies are computed as a function of the present state
+ (time step $n$) and external forcing; then within the main part of
+ model, only those new tendencies are added to the model variables.
+ [more details needed]\\
- In the standard formulation, the surface pressure is
+ The atmospheric physics follows this general scheme.
- evaluated at time step n+1 (implicit method).
- The above set of equations is then discretized in time
+ [more about C\_grid, A\_grid conversion \& drag term]\\
- as follows:
+ \subsection{Standard synchronous time stepping}
+ In the standard formulation, the surface pressure is evaluated at time
+ step n+1 (an implicit method).  Equations \ref{eq-tCsC-theta} to
+ \ref{eq-tCsC-cont} are then discretized in time as follows:
  \begin{eqnarray}
  \left[ 1 - \partial_r \kappa_v^\theta \partial_r \right]
  \theta^{n+1} & = & \theta^*
+ \label{eq-tDsC-theta}
  \\
  \left[ 1 - \partial_r \kappa_v^S \partial_r \right]
  S^{n+1} & = & S^*
+ \label{eq-tDsC-salt}
  \\
  %{b'}^{n} & = & b'(\theta^{n},S^{n},r)
  %\partial_r {\phi'_{hyd}}^{n} & = & {-b'}^{n}
  %\\
  {\phi'_{hyd}}^{n} & = & \int_{r'}^{R_o} b'(\theta^{n},S^{n},r) dr
- \label{eq-rtd-hyd}
+ \label{eq-tDsC-hyd}
  \\
  \vec{\bf v} ^{n+1}
- + \Delta t {\bf \nabla}_r B_o {\eta}^{n+1}
+ + \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
- + \epsilon_{nh} \Delta t {\bf \nabla}_r {\phi'_{nh}}^{n+1}
+ + \epsilon_{nh} \Delta t {\bf \nabla} {\phi'_{nh}}^{n+1}
  - \partial_r A_v \partial_r \vec{\bf v}^{n+1}
  & = &
  \vec{\bf v}^*
- \label{eq-rtd-hmom}
+ \label{eq-tDsC-Hmom}
  \\
  \epsilon_{fs} {\eta}^{n+1} + \Delta t
- {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o} \vec{\bf v}^{n+1} dr
+ {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o} \vec{\bf v}^{n+1} dr
  & = &
      \epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta_t (P-E)^{n}
  \nonumber
  \\
  % = \epsilon_{fs} {\eta}^{n} & + & \epsilon_{fw} \Delta_t (P-E)^{n}
- \label{eq-rtd-eta}
+ \label{eq-tDsC-eta}
  \\
  \epsilon_{nh} \left( \dot{r} ^{n+1}
  + \Delta t \partial_r {\phi'_{nh}} ^{n+1}
  \right)
  & = & \epsilon_{nh} \dot{r}^*
- \label{eq-rtd-rdot}
+ \label{eq-tDsC-Vmom}
  \\
- {\bf \nabla}_r \cdot \vec{\bf v}^{n+1} + \partial_r \dot{r}^{n+1}
+ {\bf \nabla}_h \cdot \vec{\bf v}^{n+1} + \partial_r \dot{r}^{n+1}
  & = & 0
- \label{eq-rtd-cont}
+ \label{eq-tDsC-cont}
  \end{eqnarray}
  where
  \begin{eqnarray}
-Line 159 
 S ^{n} + \Delta t G_{S} ^{(n+1/2)}
+Line 212 
 S ^{n} + \Delta t G_{S} ^{(n+1/2)}
  \\
  \vec{\bf v}^* & = &
  \vec{\bf v}^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
- + \Delta t  {\bf \nabla}_r {\phi'_{hyd}}^{(n+1/2)}
+ + \Delta t  {\bf \nabla}_h {\phi'_{hyd}}^{(n+1/2)}
  \\
  \dot{r}^* & = &
  \dot{r} ^{n} + \Delta t G_{\dot{r}} ^{(n+1/2)}
  \end{eqnarray}
- Note that implicit vertical terms (viscosity and diffusivity) are
+ Note that implicit vertical viscosity and diffusivity terms are not
- not considered as part of the "{\it G}" terms, but are
+ considered as part of the ``{\it G}'' terms, but are written
- written separately here.
+ separately here.
- To ensure a second order time discretization for both
+ The default time-stepping method is the Adams-Bashforth quasi-second
- momentum and tracer,
+ order scheme in which the ``G'' terms are extrapolated forward in time
- The "G" terms are "extrapolated" forward in time
+ from time-levels $n-1$ and $n$ to time-level $n+1/2$ and provides a
- (Adams Bashforth time stepping)
+ simple but stable algorithm:
- from the values computed at time step $n$ and $n-1$
+ \begin{equation}
- to the time $n+1/2$:
+ G^{(n+1/2)} = G^n + (1/2+\epsilon_{AB}) (G^n - G^{n-1})
- $$G^{(n+1/2)} = G^n + (1/2+\epsilon_{AB}) (G^n - G^{n-1})$$
+ \end{equation}
- A small number for the parameter $\epsilon_{AB}$ is generally used
+ where $\epsilon_{AB}$ is a small number used to stabilize the time
- to stabilize this time stepping.
+ stepping.
- In the standard non-stagger formulation,
- the Adams-Bashforth time stepping is also applied to the
- hydrostatic (pressure / geo-) potential term $\nabla_h \Phi'_{hyd}$.
- Note that presently, this term is in fact incorporated to the
- $\vec{\bf G}_{\vec{\bf v}}$ arrays ({\it gU,gV}).
- \subsection{general method}
+ In the standard non-staggered formulation, the Adams-Bashforth time
+ stepping is also applied to the hydrostatic pressure/geo-potential
- The general algorithm consist in  a "predictor step" that computes
+ gradient, $\nabla_h \Phi'_{hyd}$.  Note that presently, this term is in
- the forward tendencies ("G" terms") and all
+ fact incorporated to the $\vec{\bf G}_{\vec{\bf v}}$ arrays ({\bf
- the "first guess" values star notation):
+ gU,gV}).
- $\vec{\bf v}^*, \theta^*, S^*$ (and $\dot{r}^*$
+ \marginpar{JMC: Clarify ``this term''?}
- in non-hydrostatic mode). This is done in the routine {\it DYNAMICS}.
- Then the implicit terms that appear here on the left hand side (LHS),
- are solved as follows:
- Since implicit vertical diffusion and viscosity terms
- are independent from the barotropic flow adjustment,
- they are computed first, solving a 3 diagonal Nr x Nr linear system,
- and incorporated at the end of the {\it DYNAMICS} routines.
- Then the surface pressure and non hydrostatic pressure
- are evaluated ({\it SOLVE\_FOR\_PRESSURE});
- Finally, within a "corrector step',
- (routine {\it THE\_CORRECTION\_STEP})
- the new values of $u,v,w,\theta,S$
- are derived according to the above equations
- (see details in II.1.3).
- At this point, the regular time step is over, but
- the correction step contains also other optional
- adjustments such as convective adjustment algorithm, or filters
- (zonal FFT filter, shapiro filter)
- that applied on both momentum and tracer fields.
- just before setting the $n+1$ new time step value.
- Since the pressure solver precision is of the order of
- the "target residual" that could be lower than the
- the computer truncation error, and also because some filters
- might alter the divergence part of the flow field,
- a final evaluation of the surface r anomaly $\eta^{n+1}$
- is performed, according to \ref{eq-rtd-eta} ({\it CALC\_EXACT\_ETA}).
- This ensures a perfect volume conservation.
- Note that there is no need for an equivalent Non-hydrostatic
- "exact conservation" step, since W is already computed after
- the filters applied.
- optionnal forcing terms (package):\\
- Regarding optional forcing terms (usually part of a "package"),
- that a account for a specific source or sink term (e.g.: condensation
- as a sink of water vapor Q), they are generally incorporated
- in the main algorithm as follows;
- At the the beginning of the time step,
- the additionnal tendencies are computed
- as function of the present state (time step $n$) and external forcing ;
- Then within the main part of model,
- only those new tendencies are added to the model variables.
- [mode details needed]
+ \fbox{ \begin{minipage}{4.75in}
- The atmospheric physics follows this general scheme.
+ {\em S/R TIMESTEP} ({\em timestep.F})
- \subsection{stagger baroclinic time stepping}
+ $G_u^n$: {\bf Gu} ({\em DYNVARS.h})
- An alternative is to evaluate $\phi'_{hyd}$ with the
+ $G_u^{n-1}, u^*$: {\bf GuNm1} ({\em DYNVARS.h})
- new tracer fields, and step forward the momentum after.
- This option is known as stagger baroclinic time stepping,
- since tracer and momentum are step forward in time one after the other.
- It can be activated turning on a running flag parameter
- {\it staggerTimeStep} in file "{\it data}").
- The main advantage of this time stepping compared to a synchronous one,
- is a better stability, specially regarding internal gravity waves,
- and a very natural implementation of a 2nd order in time
- hydrostatic pressure / geo- potential term.
- In the other hand, a synchronous time step might be  better
- for convection problems; Its also make simpler time dependent forcing
- and diagnostic implementation ; and allows a more efficient threading.
- Although the stagger time step does not affect deeply the
- structure of the code --- a switch allows to evaluate the
- hydrostatic pressure / geo- potential from new $\theta,S$
- instead of the Adams-Bashforth estimation ---
- this affect the way the time discretization is presented :
- \begin{eqnarray*}
+ $G_v^n$: {\bf Gv} ({\em DYNVARS.h})
- \left[ 1 - \partial_r \kappa_v^\theta \partial_r \right]
- \theta^{n+1/2} & = & \theta^*
+ $G_v^{n-1}, v^*$: {\bf GvNm1} ({\em DYNVARS.h})
- \\
- \left[ 1 - \partial_r \kappa_v^S \partial_r \right]
+ $G_u^{(n+1/2)}$: {\bf GuTmp} (local)
- S^{n+1/2} & = & S^*
- \end{eqnarray*}
+ $G_v^{(n+1/2)}$: {\bf GvTmp} (local)
- with
+ \end{minipage} }
+ \subsection{Stagger baroclinic time stepping}
+ An alternative to synchronous time-stepping is to stagger the momentum
+ and tracer fields in time. This allows the evaluation and gradient of
+ $\phi'_{hyd}$ to be centered in time with out needing to use the
+ Adams-Bashforth extrapoltion. This option is known as staggered
+ baroclinic time stepping because tracer and momentum are stepped
+ forward-in-time one after the other.  It can be activated by turning
+ on a run-time parameter {\bf staggerTimeStep} in namelist ``{\it
+ PARM01}''.
+ The main advantage of staggered time-stepping compared to synchronous,
+ is improved stability to internal gravity wave motions and a very
+ natural implementation of a 2nd order in time hydrostatic
+ pressure/geo-potential gradient term. However, synchronous
+ time-stepping might be better for convection problems, time dependent
+ forcing and diagnosing the model are also easier and it allows a more
+ efficient parallel decomposition.
+ The staggered baroclinic time-stepping scheme is equations
+ \ref{eq-tDsC-theta} to \ref{eq-tDsC-cont} except that \ref{eq-tDsC-hyd} is replaced with:
+ \begin{equation}
+ {\phi'_{hyd}}^{n+1/2} = \int_{r'}^{R_o} b'(\theta^{n+1/2},S^{n+1/2},r)
+ dr
+ \end{equation}
+ and the time-level for tracers has been shifted back by half:
  \begin{eqnarray*}
  \theta^* & = &
  \theta ^{(n-1/2)} + \Delta t G_{\theta} ^{(n)}
  \\
  S^* & = &
  S ^{(n-1/2)} + \Delta t G_{S} ^{(n)}
- \end{eqnarray*}
- And
- \begin{eqnarray*}
- %{b'}^{n+1/2} & = & b'(\theta^{n+1/2},S^{n+1/2},r)
- %\\
- %\partial_r {\phi'_{hyd}}^{n+1/2} & = & {-b'}^{n+1/2}
- {\phi'_{hyd}}^{n+1/2} & = & \int_{r'}^{R_o} b'(\theta^{n+1/2},S^{n+1/2},r) dr
- %\label{eq-rtd-hyd}
- \\
- \vec{\bf v} ^{n+1}
- + \Delta t {\bf \nabla}_r B_o {\eta}^{n+1}
- + \epsilon_{nh} \Delta t {\bf \nabla}_r {\phi'_{nh}}^{n+1}
- - \partial_r A_v \partial_r \vec{\bf v}^{n+1}
- & = &
- \vec{\bf v}^*
- %\label{eq-rtd-hmom}
- \\
- \epsilon_{fs} {\eta}^{n+1} + \Delta t
- {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o} \vec{\bf v}^{n+1} dr
- & = &
- \epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta_t (P-E)^{n}
  \\
- \epsilon_{nh} \left( \dot{r} ^{n+1}
+ \left[ 1 - \partial_r \kappa_v^\theta \partial_r \right]
- + \Delta t \partial_r {\phi'_{nh}} ^{n+1}
+ \theta^{n+1/2} & = & \theta^*
- \right)
- & = & \epsilon_{nh} \dot{r}^*
- %\label{eq-rtd-rdot}
- \\
- {\bf \nabla}_r \cdot \vec{\bf v}^{n+1} + \partial_r \dot{r}^{n+1}
- & = & 0
- %\label{eq-rtd-cont}
- \end{eqnarray*}
- with
- \begin{eqnarray*}
- \vec{\bf v}^* & = &
- \vec{\bf v} ^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
- + \Delta t  {\bf \nabla}_r {\phi'_{hyd}}^{n+1/2}
  \\
- \dot{r}^* & = &
+ \left[ 1 - \partial_r \kappa_v^S \partial_r \right]
- \dot{r} ^{n} + \Delta t G_{\dot{r}} ^{(n+1/2)}
+ S^{n+1/2} & = & S^*
  \end{eqnarray*}
- %---------------------------------------------------------------------
- \subsection{surface pressure}
+ \subsection{Surface pressure}
- Substituting \ref{eq-rtd-hmom} into \ref{eq-rtd-cont}, assuming
+ Substituting \ref{eq-tDsC-Hmom} into \ref{eq-tDsC-cont}, assuming
  $\epsilon_{nh} = 0$ yields a Helmholtz equation for ${\eta}^{n+1}$:
- $$
+ \begin{eqnarray}
  \epsilon_{fs} {\eta}^{n+1} -
- {\bf \nabla}_r \cdot \Delta t^2 (R_o-R_{min})
+ {\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{fixed})
- {\bf \nabla}_r B_o {\eta}^{n+1}
+ {\bf \nabla}_h b_s {\eta}^{n+1}
  = {\eta}^*
- \label{solve_2d}
+ \label{eq-solve2D}
- $$
+ \end{eqnarray}
  where
- $$
+ \begin{eqnarray}
  {\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
- \Delta t {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o} \vec{\bf v}^* dr
+ \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}
+ \fbox{ \begin{minipage}{4.75in}
+ {\em S/R SOLVE\_FOR\_PRESSURE} ({\em solve\_for\_pressure.F})
+ $u^*$: {\bf GuNm1} ({\em DYNVARS.h})
+ $v^*$: {\bf GvNm1} ({\em DYNVARS.h})
- Once ${\eta}^{n+1}$ has been found substituting into \ref{eq-rtd-hmom}
+ $\eta^*$: {\bf cg2d\_b} (\em SOLVE\_FOR\_PRESSURE.h)
- would yield $\vec{\bf v}^{n+1}$ if the model is hydrostatic
- ($\epsilon_{nh}=0$):
+ $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)
+ \end{minipage} }
+ Once ${\eta}^{n+1}$ has been found, substituting into
+ \ref{eq-tDsC-Hmom} yields $\vec{\bf v}^{n+1}$ if the model is
+ hydrostatic ($\epsilon_{nh}=0$):
  $$
  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}
- - \Delta t {\bf \nabla}_r B_o {\eta}^{n+1}
+ - \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
  $$
  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
+ additional equation for $\phi'_{nh}$. This is obtained by substituting
- substituting \ref{eq-rtd-hmom} and \ref{eq-rtd-rdot} into
+ \ref{eq-tDsC-Hmom} and \ref{eq-tDsC-Vmom} into
- \ref{eq-rtd-cont}:
+ \ref{eq-tDsC-cont}:
  \begin{equation}
- \left[ {\bf \nabla}_r^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1}
+ \left[ {\bf \nabla}_h^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1}
  = \frac{1}{\Delta t} \left(
- {\bf \nabla}_r \cdot \vec{\bf v}^{**} + \partial_r \dot{r}^* \right)
+ {\bf \nabla}_h \cdot \vec{\bf v}^{**} + \partial_r \dot{r}^* \right)
  \end{equation}
  where
  \begin{displaymath}
- \vec{\bf v}^{**} = \vec{\bf v}^* - \Delta t {\bf \nabla}_r B_o {\eta}^{n+1}
+ \vec{\bf v}^{**} = \vec{\bf v}^* - \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
  \end{displaymath}
  Note that $\eta^{n+1}$ is also used to update the second RHS term
  $\partial_r \dot{r}^* $ since
-Line 369 
 is evaluated as $(\eta^{n+1} - \eta^n) /
+Line 364 
 is evaluated as $(\eta^{n+1} - \eta^n) /
  Finally, the horizontal velocities at the new time level are found by:
  \begin{equation}
  \vec{\bf v}^{n+1} = \vec{\bf v}^{**}
- - \epsilon_{nh} \Delta t {\bf \nabla}_r {\phi'_{nh}}^{n+1}
+ - \epsilon_{nh} \Delta t {\bf \nabla}_h {\phi'_{nh}}^{n+1}
  \end{equation}
  and the vertical velocity is found by integrating the continuity
- equation vertically.
+ equation vertically.  Note that, for the convenience of the restart
- Note that, for convenience regarding the restart procedure,
+ procedure, the vertical integration of the continuity equation has
- the integration of the continuity equation has been
+ been moved to the beginning of the time step (instead of at the end),
- moved at the beginning of the time step (instead of at the end),
  without any consequence on the solution.
- Regarding the implementation, all those computation are done
+ \fbox{ \begin{minipage}{4.75in}
- within the routine {\it SOLVE\_FOR\_PRESSURE} and its dependent
+ {\em S/R CORRECTION\_STEP} ({\em correction\_step.F})
- {\it CALL}s.
- The standard method to solve the 2D elliptic problem (\ref{solve_2d})
+ $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)
- uses the conjugate gradient method (routine {\it CG2D}); The
- the solver matrix and conjugate gradient operator are only function
+ $\phi_{nh}^{n+1}$: {\bf phi\_nh} (\em DYNVARS.h)
- of the discretized domain and are therefore evaluated separately,
- before the time iteration loop, within {\it INI\_CG2D}.
+ $u^*$: {\bf GuNm1} ({\em DYNVARS.h})
- The computation of the RHS $\eta^*$ is partly
- done in {\it CALC\_DIV\_GHAT} and in {\it SOLVE\_FOR\_PRESSURE}.
+ $v^*$: {\bf GvNm1} ({\em DYNVARS.h})
- The same method is applied for the non hydrostatic part, using
+ $u^{n+1}$: {\bf uVel} ({\em DYNVARS.h})
- a conjugate gradient 3D solver ({\it CG3D}) that is initialized
- in {\it INI\_CG3D}. The RHS terms of 2D and 3D problems
+ $v^{n+1}$: {\bf vVel} ({\em DYNVARS.h})
- are computed together, within the same part of the code.
+ \end{minipage} }
+ Regarding the implementation of the surface pressure solver, all
+ computation are done within the routine {\it SOLVE\_FOR\_PRESSURE} and
+ its dependent calls.  The standard method to solve the 2D elliptic
+ problem (\ref{eq-solve2D}) uses the conjugate gradient method (routine
+ {\it CG2D}); the solver matrix and conjugate gradient operator are
+ only function of the discretized domain and are therefore evaluated
+ separately, before the time iteration loop, within {\it INI\_CG2D}.
+ The computation of the RHS $\eta^*$ is partly done in {\it
+ CALC\_DIV\_GHAT} and in {\it SOLVE\_FOR\_PRESSURE}.
+ The same method is applied for the non hydrostatic part, using a
+ conjugate gradient 3D solver ({\it CG3D}) that is initialized in {\it
+ INI\_CG3D}. The RHS terms of 2D and 3D problems are computed together
+ at the same point in the code.
- \newpage
- %-----------------------------------------------------------------------------------
  \subsection{Crank-Nickelson barotropic time stepping}
- The full implicit time stepping described previously is unconditionally stable
+ The full implicit time stepping described previously is
- but damps the fast gravity waves, resulting in a loss of
+ unconditionally stable but damps the fast gravity waves, resulting in
- gravity potential energy.
+ a loss of potential energy.  The modification presented now allows one
- The modification presented hereafter allows to combine an implicit part
+ to combine an implicit part ($\beta,\gamma$) and an explicit part
- ($\beta,\gamma$) and an explicit part ($1-\beta,1-\gamma$) for the surface
+ ($1-\beta,1-\gamma$) for the surface pressure gradient ($\beta$) and
- pressure gradient ($\beta$) and for the barotropic flow divergence ($\gamma$).
+ for the barotropic flow divergence ($\gamma$).
  \\
  For instance, $\beta=\gamma=1$ is the previous fully implicit scheme;
  $\beta=\gamma=1/2$ is the non damping (energy conserving), unconditionally
-Line 414 
 In the code, $\beta,\gamma$ are defined
+Line 426 
 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-rtd-hmom} and \ref{eq-rtd-eta} are modified as follows:
+ Equations \ref{eq-tDsC-Hmom} and \ref{eq-tDsC-eta} are modified as follows:
  $$
  \frac{ \vec{\bf v}^{n+1} }{ \Delta t }
- + {\bf \nabla}_r B_o [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ]
+ + {\bf \nabla}_h b_s [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ]
- + \epsilon_{nh} {\bf \nabla}_r {\phi'_{nh}}^{n+1}
+ + \epsilon_{nh} {\bf \nabla}_h {\phi'_{nh}}^{n+1}
   = \frac{ \vec{\bf v}^* }{ \Delta t }
  $$
  $$
  \epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t}
- + {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o}
+ + {\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)
  $$
-Line 431 
 where:
+Line 443 
 where:
  \begin{eqnarray*}
  \vec{\bf v}^* & = &
  \vec{\bf v} ^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
- + (\beta-1) \Delta t {\bf \nabla}_r B_o {\eta}^{n}
+ + (\beta-1) \Delta t {\bf \nabla}_h b_s {\eta}^{n}
- + \Delta t {\bf \nabla}_r {\phi'_{hyd}}^{(n+1/2)}
+ + \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}_r \cdot \int_{R_{min}}^{R_o}
+ - \Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o}
  [ \gamma \vec{\bf v}^* + (1-\gamma) \vec{\bf v}^{n}] dr
  \end{eqnarray*}
  \\
- In the hydrostatic case ($\epsilon_{nh}=0$),
+ In the hydrostatic case ($\epsilon_{nh}=0$), allowing us to find
- this allow to find ${\eta}^{n+1}$, according to:
+ ${\eta}^{n+1}$, thus:
  $$
  \epsilon_{fs} {\eta}^{n+1} -
- {\bf \nabla}_r \cdot \beta\gamma \Delta t^2 B_o (R_o - R_{min})
+ {\bf \nabla}_h \cdot \beta\gamma \Delta t^2 b_s (R_o - R_{fixed})
- {\bf \nabla}_r {\eta}^{n+1}
+ {\bf \nabla}_h {\eta}^{n+1}
  = {\eta}^*
  $$
  and then to compute (correction step):
  $$
  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}
- - \beta \Delta t {\bf \nabla}_r B_o {\eta}^{n+1}
+ - \beta \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
  $$
  The non-hydrostatic part is solved as described previously.
- \\ \\
- N.B:
+ Note that:
- \\
+ \begin{enumerate}
-  a) The non-hydrostatic part of the code has not yet been
+ \item The non-hydrostatic part of the code has not yet been
- updated, %since it falls out of the purpose of this test,
+ updated, so that this option cannot be used with $(\beta,\gamma) \neq (1,1)$.
- so that this option cannot be used with $(\beta,\gamma) \neq (1,1)$.
+ \item The stability criteria with Crank-Nickelson time stepping
- \\
+ for the pure linear gravity wave problem in cartesian coordinates is:
- b) To remind, the stability criteria with the Crank-Nickelson time stepping
+ \begin{itemize}
- for the pure linear gravity wave problem in cartesian coordinate is:
+ \item $\beta + \gamma < 1$ : unstable
- \\
+ \item $\beta \geq 1/2$ and $ \gamma \geq 1/2$ : stable
- $\star$~ $\beta + \gamma < 1$ : unstable
+ \item $\beta + \gamma \geq 1$ : stable if
- \\
- $\star$~ $\beta \geq 1/2$ and $ \gamma \geq 1/2$ : stable
- \\
- $\star$~ $\beta + \gamma \geq 1$ : stable if
- %, for all wave length $(k\Delta x,l\Delta y)$
  $$
  c_{max}^2 (\beta - 1/2)(\gamma - 1/2) + 1 \geq 0
  $$
-Line 485 
 $$
+Line 492 
 $$
  c_{max} =  2 \Delta t \: \sqrt{g H} \:
  \sqrt{ \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} }
  $$
+ \end{itemize}
+ \end{enumerate}

 Legend:



Removed from v.1.1.1.1
 


changed lines


 
Added in v.1.7
 Legend:



Removed from v.1.1.1.1
 


changed lines


 
Added in v.1.7
-Removed from v.1.1.1.1
+Added in v.1.7

	ViewVC Help
Powered by ViewVC 1.1.22