/[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.3 by jmc,
Thu Aug  9 20:06:18 2001 UTC
+revision 1.5 by jmc,
Mon Sep 24 19:30:40 2001 UTC
 Line 21 
 The continuous form of the model equatio
  \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}_h b_s \eta
  + \epsilon_{nh} {\bf \nabla}_h \phi'_{nh}
  & = & \vec{\bf G}_{\vec{\bf v}}
  - {\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}_h \cdot \vec{\bf v} + \partial_r \dot{r}
  & = & 0
- \label{eq-r-cont}
+ \label{eq-tCsC-cont}
  \end{eqnarray}
  where
  \begin{eqnarray*}
-Line 71 
 all the RHS terms except the pressure /
+Line 73 
 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 equation \ref{eq-r-split-rdot} vanishes.
+ and equation \ref{eq-tCsC-Vmom} vanishes.
  The equation for $\eta$ is obtained by integrating the
  continuity equation over the entire depth of the fluid,
- from $R_{min}(x,y)$ up to $R_o(x,y)$
+ from $R_{fixed}(x,y)$ up to $R_o(x,y)$
  (Linear free surface):
  \begin{eqnarray}
  \epsilon_{fs} \partial_t \eta =
  \left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) =
- - {\bf \nabla} \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)
- \label{eq-cont-2D}
+ \label{eq-tCsC-eta}
  \end{eqnarray}
  Where $\epsilon_{fs}$,$\epsilon_{fw}$ are two flags to
-Line 93 
 at the ocean surface (natural BC) ($\eps
+Line 95 
 at the ocean surface (natural BC) ($\eps
  or not ($\epsilon_{fw} = 0$).
  The hydrostatic potential is found by
- integrating \ref{eq-r-split-hyd} with the boundary condition that
+ integrating \ref{eq-tCsC-hyd} with the boundary condition that
  $\phi'_{hyd}(r=R_o) = 0$:
  \begin{eqnarray*}
  & &
-Line 107 
 $\phi'_{hyd}(r=R_o) = 0$:
+Line 109 
 $\phi'_{hyd}(r=R_o) = 0$:
  \subsection{General method}
+ An overview of the general method is presented hereafter,
+ with explicit references to the Fortran code. This part
+ often refers to the discretized equations of the model
+ that are detailed in the following sections.
  The general algorithm consist in  a "predictor step" that computes
  the forward tendencies ("G" terms") and all
- the "first guess" values star notation):
+ the "first guess" values (star notation):
- $\vec{\bf v}^*, \theta^*, S^*$ (and $\dot{r}^*$
+ $\theta^*, S^*, \vec{\bf v}^*$ (and $\dot{r}^*$
- in non-hydrostatic mode). This is done in the routine {\it DYNAMICS}.
+ in non-hydrostatic mode). This is done in the two routines
+ {\it THERMODYNAMICS} and {\it DYNAMICS}.
- Then the implicit terms that appear here on the left hand side (LHS),
+ Then the implicit terms that appear on the left hand side (LHS)
+ of equations \ref{eq-tDsC-theta} - \ref{eq-tDsC-cont},
  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.
+ and incorporated at the end of the {\it THERMODYNAMICS} and
+ {\it DYNAMICS} routines.
  Then the surface pressure and non hydrostatic pressure
  are evaluated ({\it SOLVE\_FOR\_PRESSURE});
-Line 140 
 the "target residual" that could be lowe
+Line 150 
 the "target residual" that could be lowe
  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}).
+ is performed, according to \ref{eq-tDsC-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
+ that 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
+ the additional 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.
+ [more about C\_grid, A\_grid conversion \& drag term]\\
  \subsection{Standard synchronous time stepping}
  In the standard formulation, the surface pressure is
-Line 169 
 as follows:
+Line 181 
 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}_h b_s {\eta}^{n+1}
-Line 185 
 S^{n+1} & = & S^*
+Line 199 
 S^{n+1} & = & S^*
  - \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}_h \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}_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 228 
 written separately here.
+Line 242 
 written separately here.
  To ensure a second order time discretization for both
  momentum and tracer,
- The "G" terms are "extrapolated" forward in time
+ The "{\it G}" terms are "extrapolated" forward in time
  (Adams Bashforth time stepping)
  from the values computed at time step $n$ and $n-1$
  to the time $n+1/2$:
-Line 286 
 And
+Line 300 
 And
  %\\
  %\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}
+ %\label{eq-tDsC-hyd}
  \\
  \vec{\bf v} ^{n+1}
  + \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
-Line 294 
 And
+Line 308 
 And
  - \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}_h \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}
  \\
-Line 305 
 And
+Line 319 
 And
  + \Delta t \partial_r {\phi'_{nh}} ^{n+1}
  \right)
  & = & \epsilon_{nh} \dot{r}^*
- %\label{eq-rtd-rdot}
+ %\label{eq-tDsC-Vmom}
  \\
  {\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*}
  with
  \begin{eqnarray*}
-Line 325 
 with
+Line 339 
 with
  \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}_h \cdot \Delta t^2 (R_o-R_{min})
+ {\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{fixed})
  {\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}_h \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{solve_2d_rhs}
+ \label{eq-solve2D_rhs}
  \end{eqnarray}
- Once ${\eta}^{n+1}$ has been found substituting into \ref{eq-rtd-hmom}
+ Once ${\eta}^{n+1}$ has been found substituting into \ref{eq-tDsC-Hmom}
  would yield $\vec{\bf v}^{n+1}$ if the model is hydrostatic
  ($\epsilon_{nh}=0$):
  $$
-Line 353 
 $$
+Line 367 
 $$
  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
- substituting \ref{eq-rtd-hmom} and \ref{eq-rtd-rdot} into
+ substituting \ref{eq-tDsC-Hmom} and \ref{eq-tDsC-Vmom} into
- \ref{eq-rtd-cont}:
+ \ref{eq-tDsC-cont}:
  \begin{equation}
  \left[ {\bf \nabla}_h^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1}
  = \frac{1}{\Delta t} \left(
-Line 384 
 without any consequence on the solution.
+Line 398 
 without any consequence on the solution.
  Regarding the implementation, all those computation are done
  within the routine {\it SOLVE\_FOR\_PRESSURE} and its dependent
  {\it CALL}s.
- The standard method to solve the 2D elliptic problem (\ref{solve_2d})
+ The standard method to solve the 2D elliptic problem (\ref{eq-solve2D})
  uses the conjugate gradient method (routine {\it CG2D}); The
  the solver matrix and conjugate gradient operator are only function
  of the discretized domain and are therefore evaluated separately,
-Line 417 
 In the code, $\beta,\gamma$ are defined
+Line 431 
 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}_h b_s [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ]
-Line 426 
 $$
+Line 440 
 $$
  $$
  $$
  \epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t}
- + {\bf \nabla}_h \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 439 
 where:
+Line 453 
 where:
  \\
  {\eta}^* & = &
  \epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta t (P-E)
- - \Delta t {\bf \nabla}_h \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*}
  \\
-Line 447 
 In the hydrostatic case ($\epsilon_{nh}=
+Line 461 
 In the hydrostatic case ($\epsilon_{nh}=
  this allow to find ${\eta}^{n+1}$, according to:
  $$
  \epsilon_{fs} {\eta}^{n+1} -
- {\bf \nabla}_h \cdot \beta\gamma \Delta t^2 b_s (R_o - R_{min})
+ {\bf \nabla}_h \cdot \beta\gamma \Delta t^2 b_s (R_o - R_{fixed})
  {\bf \nabla}_h {\eta}^{n+1}
  = {\eta}^*
  $$

 Legend:



Removed from v.1.3
 


changed lines


 
Added in v.1.5
 Legend:



Removed from v.1.3
 


changed lines


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

	ViewVC Help
Powered by ViewVC 1.1.22