%\section{Numerical procedures of Diecast model}
\chapter{Numerical Procedures}
\label{app:numerics}

We present the details of DieCAST hydrostatic model
in this appendix. The governing equations have

Conservation of mass:
\be
\nabla \cdot \mathbf{V}=0
\label{eq:3Dmass}
\ee

Horizontal momentum equations:
\begin{align}
\fc{\pd u}{\pd t}&= -\nabla \cdot u\mathbf{V}+fv
-\fc{1}{\rho_0}\fc{\pd p}{\pd x}+\nabla_{h} \cdot A_h\nabla_{h}u
+\fc{\pd}{\pd z}(A_v\fc{\pd u}{\pd z}) \label{eq:3DNS_U}\\
\fc{\pd v}{\pd t}&= -\nabla \cdot v\mathbf{V}-fu
-\fc{1}{\rho_0}\fc{\pd p}{\pd y}+\nabla_{h} \cdot A_h\nabla_{h}v
+\fc{\pd}{\pd z}(A_v\fc{\pd v}{\pd z}) \label{eq:3DNS_V}
\end{align}

Conservation of scalar (salt or potential temperature):
\be
\fc{\pd S}{\pd t}= -\nabla \cdot S\mathbf{V}
+\nabla_{h}\cdot K_h\nabla_{h} S
+\fc{\pd}{\pd z}(K_v\fc{\pd S}{\pd z}) \label{eq:3DNS_S}
\ee

Hydrostatic equation:
\be
\fc{\pd p}{\pd z}=-(\rho-\overline{\rho})g
\label{eq:Hydrostatic}
\ee

Equation of state:
\be
\rho=\rho(S,T)
\ee
where $u$ and $v$ are the $x$ and $y$ components of velocity vector
$\mathbf{V}=(u,v,w)$ respectively; $f$ is the Coriolis parameter;
$\rho_0$ is the mean density; $\overline{\rho}$ is the horizontally
averaged density at depth $z$; $p$ is the pressure; $A_h$, $A_v$ are the
horizontal and vertical eddy viscosity respectively; $S$ is the salinity;
$K_h$, $K_v$ are the horizontal and vertical eddy diffusivity respectively;
$T$ is the potential temperature.

%Rewrite the equation~(\ref{eq:3DNS_U}) in numerical discretization form
%and the staggered grid arrangement is shown in Figure~\ref{fig:diecast_
From the staggered grid arrangement shown in
Figure~\ref{fig:diecast_cell},
the control volume and time integrated form of equation (\ref{eq:3DNS_U}) may
be written,
\be
\begin{split}
\fc{\tilde{u}_{i,j,k}^{n+1}-u_{i,j,k}^{n-1}}{\Delta t}=
&-(\fc{U_{i,j,k}^{n}u_{i+1/2,j,k}^{n}-U_{i-1,j,k}^{n}u_{i-1/2,j,k}^{n}}{\Delta x_{j}}\\
&+\fc{V_{i,j,k}^{n}u_{i,j+1/2,k}^{n}-V_{i,j-1,k}^{n}u_{i,j-1/2,k}^{n}}{\Delta y_{j}}\\
&+\fc{W_{i,j,k+1}^{n}u_{i,j,k+1/2}^{n}-W_{i,j,k}^{n}u_{i,j,k-1/2}^{n}}{\Delta z_{k}})\\
&+[2\Omega_{e}sin(\phi_{j})+u_{i,j,k}^{n}\tan(\fc{\phi_{j}}{r_{e}})]v_{i,j,k}^{n} \\
&-(\fc{\pd \tilde{p}}{\pd x})_{i,j,k}^{n}+ dissipation
\end{split}
\label{eq:numerical_discritization}
\ee
where the pressure derivative term is
\be
\begin{split}
\fc{\pd \tilde{p}}{\pd x}_{i,j,k}^{n}&=\fc{7}{12}
(\fc{\tilde{p}_{i+1,j,k}-\tilde{p}_{i,j,k}}{\Delta x_{j}}
+\fc{\tilde{p}_{i,j,k}-\tilde{p}_{i-1,j,k}}{\Delta x_{j}}) \\
&-\fc{1}{12}
(\fc{\tilde{p}_{i+2,j,k}-\tilde{p}_{i+1,j,k}}{\Delta x_{j}}
+\fc{\tilde{p}_{i-1,j,k}-\tilde{p}_{i-2,j,k}}{\Delta x_{j}})\\
\end{split}
\ee
and the dissipation term is the viscous and quadratic bottom drag
(bottom layer only) terms.
In addition, a wind stress is applied to the top layer.
Although standard Earth spherical coordinates are assumed, we do not
include the cosine factors here (horizontal grid increments vary
with latitude only and vertical grid increments with depth only).
The quantity $\tilde{q}$ means the trial value of field $q$ derived from
the trial pressure at the rigid-lid. The rigid-lid pressure
is hydrostatically related to the slow mode surface height.

In order to evaluate a cell face average quantity,
fourth order accurate interpolation is used
(see Appendix~\ref{app:fourth_order}):
Then $\overline{q}_{i+1/2,j,k}$ represents the cell face quantity derived 
from the four nearest control volumes, for example:
\be
\overline{q}_{i+1/2,j,k}=\fc{7}{12}(Q_{i,j,k}+Q_{i+1,j,k})
-\fc{1}{12}(Q_{i-1,j,k}+Q_{i+2,j,k})
\label{eq:interp}
\ee
where $\overline{q}_{i+1/2,j,k}$ represents the cell face average quantity
derived from its four nearest control volume averages on the right hand
side of equation~(\ref{eq:interp}). 
After determining the pressure field $\tilde{p}$
from the hydrostatic equation (\ref{eq:Hydrostatic})
using a trial (or most recent) value of rigid-lid pressure,
the control volume conservation equations are used to get
$\tilde{u},\tilde{v},S,T$ at the advanced time level $(n+1)$.
All terms except diffusion are evaluated from fields known at leapfrog time
level $n$; diffusion terms are at time level $n-1$.

The numerical procedure for the incompressible correction is summarized
in the following:
\bn

\item Guess the trial surface pressure $\tilde{p}_{s}^{n}$ from the previous
time step ($\tilde{p}_{s}^{n} =p_{s}^{n-1}$) and integrate the hydrostatic 
equation~(\ref{eq:Hydrostatic}) to get the intermediate pressure field 
$\tilde{p}^{n}$ over the whole domain.

\item Update the trial integral average $\tilde{u}^{n+1}$,$\tilde{v}^{n+1}$
in the control volume using equation~(\ref{eq:numerical_discritization}).

\item Interpolate $\tilde{u}^{n+1},\tilde{u}^{n+1}$ to 
$\tilde{U}_{i,j,k}^{n+1},\tilde{V}_{i,j,k}^{n+1}$
at the cell face using $4^{th}$ order interpolation from
equation~(\ref{eq:interp}).

\item Adjust the rigid-lid pressure such that the divergence
of the vertically integrated horizontal velocity matches the specified value
(near zero, but small precipitation/evaporation and river sources are included
to give non-zero $W$ at the rigid lid which is slightly porous to fresh water).
This procedure uses mass conservation and the detailed procedure is listed
as follows.
\bi
\item Set the final pressure $p^{n}=\tilde{p}^{n}+\Delta \tilde{p}$
where $\Delta \tilde{p}$ is due to the change of rigid-lid pressure
and thus independent of depth. $\tilde{p}^{n}$ is derived hydrostatically
from trial rigid-lid pressure in the first step. Then the final velocities 
can be written as
\begin{align}
U^{n+1}&=\tilde{U}^{n+1}+\Delta \tilde{U} \label{eq:U_final} \\
V^{n+1}&=\tilde{V}^{n+1}+\Delta \tilde{V},
\label{eq:V_final} 
\end{align}
where
$\Delta \tilde{U}=-\Delta t \fc{\pd \Delta \tilde{p}}{\pd x}$ and
$\Delta \tilde{V}=-\Delta t \fc{\pd \Delta \tilde{p}}{\pd y}$. 
\item The final velocity field satisfies mass conservation from equation
(\ref{eq:3Dmass}) so that the vertically-integrated horizontal velocity
can be determined by vertical velocity at the top and bottom:
\be
\int_{0}^{D}(\fc{\pd U^{n+1}}{\pd x}+\fc{\pd V^{n+1}}{\pd y})dz=
W^{n+1}(0)-W^{n+1}(D)=0
\ee

\item Substitute the final velocity integration in terms of trial velocity and
the correction in equations (\ref{eq:U_final}) and (\ref{eq:V_final}),
thus deriving a 2-D Poisson equation for rigid-lid pressure adjustment. Solve
the Poisson equation (\ref{eq:2d_poisson}) using an efficient EVP elliptic solver~\cite{Roache95}.
\be
\int_{0}^{D}(\fc{\pd \tilde{U}^{n+1}}{\pd x}+\fc{\pd \tilde{V}^{n+1}}{\pd y})dz=
\int_{0}^{D}-(\fc{\pd \Delta \tilde{U}}{\pd x}+\fc{\pd \Delta \tilde{V}}{\pd y})dz
\label{eq:2d_poisson}
\ee
\ei

\item The velocity $U^{n+1},V^{n+1}$ is corrected from the pressure adjustment
$\Delta \tilde{p}$ in the last step.

\item Finally, interpolate cell average changes
$\Delta \tilde{u},\Delta \tilde{v}$
from $\Delta \tilde{U}$,$\Delta \tilde{V}$ using $4^{th}$ order
interpolation, then
$u^{n+1}=\tilde{u}^{n+1}+\Delta \tilde{u}$ and
$v^{n+1}=\tilde{v}^{n+1}+\Delta \tilde{v}$.
\en

\begin{figure}[htbp]
\centerline{
\scalebox{.4}{\includegraphics{FIGS/diecast_cell.eps}}
}
\vskip 0.1in
\caption{The staggered and collocated arrangement of grids.
$Q(i,j,k)$ is any cell-averaged quantity, $u,v,s,t,p$ or $\rho$.
$U,V,W$ are the face-averaged quantities.}
\label{fig:diecast_cell}
\end{figure}
