/[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.2 by jmc,
Thu Aug  9 16:22:09 2001 UTC
+revision 1.3 by jmc,
Thu Aug  9 20:06:18 2001 UTC
 Line 8 
 at time $t = n \Delta t$ . We used the 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
 Line 30 
 b' & = & b'(\theta,S,r)
  \label{eq-r-split-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}
  \\
  \epsilon_{nh} \frac {\partial{\dot{r}}}{\partial{t}}
 Line 41 
 b' & = & b'(\theta,S,r)
  & = & \epsilon_{nh} G_{\dot{r}}
  \label{eq-r-split-rdot}
  \\
- {\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}
  \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
 Line 70 
 all the RHS terms except the pressure /
  The switch $\epsilon_{nh}$ allows to activate the non hydrostatic
  mode ($\epsilon_{nh}=1$) for the ocean model. Otherwise,
- in the hydrostatic limit $\epsilon_{nh} = 0$ and the 3rd equation vanishes.
+ in the hydrostatic limit $\epsilon_{nh} = 0$
+ and equation \ref{eq-r-split-rdot} vanishes.
  The equation for $\eta$ is obtained by integrating the
  continuity equation over the entire depth of the fluid,
-Line 79 
 from $R_{min}(x,y)$ up to $R_o(x,y)$
+Line 80 
 from $R_{min}(x,y)$ up to $R_o(x,y)$
  \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_{min}}^{R_o} \vec{\bf v} dr
  + \epsilon_{fw} (P-E)
  \label{eq-cont-2D}
  \end{eqnarray}
-Line 104 
 $\phi'_{hyd}(r=R_o) = 0$:
+Line 105 
 $\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}
+ The general algorithm consist in  a "predictor step" that computes
+ the forward tendencies ("G" terms") and all
+ the "first guess" values star notation):
+ $\vec{\bf v}^*, \theta^*, S^*$ (and $\dot{r}^*$
+ 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.
+ [more details needed]\\
+ The atmospheric physics follows this general scheme.
+ \subsection{Standard synchronous time stepping}
  In the standard formulation, the surface pressure is
  evaluated at time step n+1 (implicit method).
-Line 124 
 S^{n+1} & = & S^*
+Line 180 
 S^{n+1} & = & S^*
  \label{eq-rtd-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}
  \\
  \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_{min}}^{R_o} \vec{\bf v}^{n+1} dr
  & = &
      \epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta_t (P-E)^{n}
  \nonumber
-Line 146 
 S^{n+1} & = & S^*
+Line 202 
 S^{n+1} & = & S^*
  & = & \epsilon_{nh} \dot{r}^*
  \label{eq-rtd-rdot}
  \\
- {\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}
  \end{eqnarray}
-Line 160 
 S ^{n} + \Delta t G_{S} ^{(n+1/2)}
+Line 216 
 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)}
-Line 184 
 In the standard non-stagger formulation,
+Line 240 
 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}).
+ $\vec{\bf G}_{\vec{\bf v}}$ arrays ({\bf gU,gV}).
- \subsection{general method}
- The general algorithm consist in  a "predictor step" that computes
- the forward tendencies ("G" terms") and all
- the "first guess" values star notation):
- $\vec{\bf v}^*, \theta^*, S^*$ (and $\dot{r}^*$
- 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]
- The atmospheric physics follows this general scheme.
- \subsection{stagger baroclinic time stepping}
+ \subsection{Stagger baroclinic time stepping}
  An alternative is to evaluate $\phi'_{hyd}$ with the
  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}").
+ {\bf 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,
-Line 288 
 And
+Line 289 
 And
  %\label{eq-rtd-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}_h {\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
+ {\bf \nabla}_h \cdot \int_{R_{min}}^{R_o} \vec{\bf v}^{n+1} dr
  & = &
  \epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta_t (P-E)^{n}
  \\
-Line 306 
 And
+Line 307 
 And
  & = & \epsilon_{nh} \dot{r}^*
  %\label{eq-rtd-rdot}
  \\
- {\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}
  \end{eqnarray*}
-Line 314 
 with
+Line 315 
 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}
+ + \Delta t  {\bf \nabla}_h {\phi'_{hyd}}^{n+1/2}
  \\
  \dot{r}^* & = &
  \dot{r} ^{n} + \Delta t G_{\dot{r}} ^{(n+1/2)}
-Line 322 
 with
+Line 323 
 with
  %---------------------------------------------------------------------
- \subsection{surface pressure}
+ \subsection{Surface pressure}
  Substituting \ref{eq-rtd-hmom} into \ref{eq-rtd-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_{min})
- {\bf \nabla}_r B_o {\eta}^{n+1}
+ {\bf \nabla}_h b_s {\eta}^{n+1}
  = {\eta}^*
  \label{solve_2d}
  \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_{min}}^{R_o} \vec{\bf v}^* dr
  \: + \: \epsilon_{fw} \Delta_t (P-E)^{n}
  \label{solve_2d_rhs}
  \end{eqnarray}
-Line 346 
 would yield $\vec{\bf v}^{n+1}$ if the m
+Line 347 
 would yield $\vec{\bf v}^{n+1}$ if the m
  ($\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
-Line 355 
 additional equation for $\phi'_{nh}$. Th
+Line 356 
 additional equation for $\phi'_{nh}$. Th
  substituting \ref{eq-rtd-hmom} and \ref{eq-rtd-rdot} into
  \ref{eq-rtd-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 371 
 is evaluated as $(\eta^{n+1} - \eta^n) /
+Line 372 
 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.
-Line 419 
 the main data file "{\it data}" and are
+Line 420 
 the main data file "{\it data}" and are
  Equations \ref{eq-rtd-hmom} and \ref{eq-rtd-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_{min}}^{R_o}
  [ \gamma \vec{\bf v}^{n+1} + (1-\gamma) \vec{\bf v}^{n}] dr
  = \epsilon_{fw} (P-E)
  $$
-Line 433 
 where:
+Line 434 
 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_{min}}^{R_o}
  [ \gamma \vec{\bf v}^* + (1-\gamma) \vec{\bf v}^{n}] dr
  \end{eqnarray*}
  \\
-Line 446 
 In the hydrostatic case ($\epsilon_{nh}=
+Line 447 
 In the hydrostatic case ($\epsilon_{nh}=
  this allow to find ${\eta}^{n+1}$, according to:
  $$
  \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_{min})
- {\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.

 Legend:



Removed from v.1.2
 


changed lines


 
Added in v.1.3
 Legend:



Removed from v.1.2
 


changed lines


 
Added in v.1.3
-Removed from v.1.2
+Added in v.1.3

	ViewVC Help
Powered by ViewVC 1.1.22