| 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, |
| 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 \\ |
| 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 algorihtm} |
\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 |
|
|
| 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 inversion 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}. |
| 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 |
Equation for $w^{n+1}$ will be here as will 3-D elliptic equations. |
The entire algorithm can be summarized as the sequential solution of |
| 649 |
\label{eq:discrete-time-w} |
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 |
|
|
| 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:discrete-time-u}, \ref{eq:discrete-time-v} and \ref{eq:discrete-time-w} |
\ref{eq:discrete-time-u-nh}, \ref{eq:discrete-time-v-nh} and \ref{eq:discrete-time-w-nh} |
| 729 |
into continuity: |
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} |
Parent Directory
| 
Revision Log
| 
Revision Graph
| 
Patch