% Governing Equations.

\chapter{Model Descriptions}
\label{chap:equation}
 
The numerical detail and formulations of TIMCOM with Boussinsq approximation 
are explained in this chapter. The hydrostatic, rigid-lid approximation version 
documented herein is freely available. Hydrostatic free surface and 
nonhydrostatic members are available by license agreement. The free surface 
version uses an efficient `truly leapfrog' (interlaced in time) approach for the 
time-split barotropic `fast mode'. The fast mode is governed by the shallow 
water equations, including forcing terms from the baroclinic mode. This chapter 
provides the governing equation and some numerical procedures. 

\section{Basic TIMCOM Ocean Model Description}

TIMCOM is a $Z$-level model with geopotential is its vertical coordinate. 
Geopotential, longitude and latitude iso-surfaces define a set of `control 
volume boxes' or `cells' in a spherical coordinate system. The conservation laws 
are integrated in each cell to get `control volume' (cell averaged) equations 
governing cell mean values. Reynolds stress turbulence closure for subgrid-scale 
transports across the cell faces is parameterized by user-defined eddy viscosity 
and diffusivities, which are defined as: 

a) for vertical cell faces, a modified~\citet{Pacanowski81} approach 
specifies vertical viscosity and diffusivities based on the Richardson number. 

b) for horizontal cell faces, the default is constant horizontal viscosity and 
diffusivities; however, the user may specify values in a three dimensional 
scratch array, which is then interpolated to the cell faces. 

The same set of cells is used for all conservation quantities; momentum density 
(velocity) component cells are NOT staggered relative to pressure, temperature, 
salinity and other predicted scalar quantities. This avoids interpolation errors 
in evaluating the large Coriolis term. 

\section{Computational Grid}
\label{sec:computational_grid}

The resolution parameters, or the number of control volumes in each direction in 
TIMCOM are defined as $I0$, $J0$, and $K0$, 
as shown in Figure~\ref{fig:TIMCOM_Grid}. 

\begin{figure}[htbp]
\centerline{\scalebox{.4}{\includegraphics{FIGS/TIMCOM_Grid.eps}}}
\caption{The grid resolution parameters in TIMCOM are $I0$, $J0$, $K0$ in $X$, 
$Y$ and $Z$ direction, 
respectively, which define the number of control volume along
the three directions. Thus, the number of control volume interfaces, or the 
`layers', along the three directions,
could be derived as: $I1=I0-1$, $J1=J0-1$ and $K1=K0-1$. The colored-grids
represent the stair-approximation of the bathymetric map.}
\label{fig:TIMCOM_Grid}
\end{figure}

This includes `ghost zones' at $I=1$, $I=I0$, $J=1$ and $J=J0$, which are used 
for convenient specification of inflows and outflows at open boundary points. No 
ghost zones are used in the vertical direction. 

The conservation equations are applied only to the interior: control volumes 
that actually contain water. In the modified `A' grid TIMCOM model these control 
volumes are the same for all fields. The `C' grid SOMS model, from which TIMCOM 
evolved~\citep{Dietrich87}, is similar except the control volumes for the
horizontal momentum conservation are laterally staggered $1/2$ grid interval 
from 
the control volumes for scalar quantities pressure, temperature and salinity. 
The modeled fluid, to which the collocated-control-volume-integrated 
conservation laws are applied, is contained in the between $I=2$ and $I=I1$
in the $X$ direction, between $J=2$ and 
$J1$ in the $Y$ direction, and between $K=1$ and $K1$ in the $Z$ direction. 
Thus, the number of interior control volumes is: 

\begin{itemize}
\item $I0-2$ ($I2$ in model variables) in the $X$ direction 

\item $J0-2$ ($J2$ in model variables) in the $Y$ direction 

\item $K0-2$ ($K2$ in model variables) in the vertical ($Z$) direction 
\end{itemize}

The definition of the 3-D `A' grid model is shown in 
Figure ~\ref{fig:diecast_cell}.
The temperature, salinity, and pressure are given as the average values over the 
volume of each box or model cell, and a control volume approach is used in which 
the average values of quantities within each cell are updated in accordance with 
the fluxes across the six cell faces. 
$U$ is determined for the cell faces normal to the $X$-direction, $V$ for cell 
faces normal to the $Y$-direction and $W$ for cell faces normal to the $Z$-
direction. 
The horizontal velocity components $U1$, $U2$, $V1$, $V$2 are specified at the 
same locations as $T$, $S$, and $P$, representing averages over the cell volume.

\begin{figure}[htbp]
\centerline{\scalebox{.4}{\includegraphics{FIGS/diecast_cell}}}
\caption{The collocated arrangement of `A' grid model.
$Q(i,j,k)$ may be any cell-averaged quantity located in the center of the cell,
i.e. $\rho$, $T$, $p$, $u$, $v$, $w$, while $U$, $V$, $W$ are face-averaged velocities, located in the center
of each cell face.}
\label{fig:diecast_cell}
\end{figure}

\subsection{Model Mask Array}

The following is critical to the understanding and application of TIMCOM. A 
master mask array called `IN' defines all boundaries. Together with the metrics 
arrays, this defines the region modeled. For scalar cell $I$, $J$, $K$ filled 
with water, $IN(I,J,K)=1$, while $IN$ is zero for cell on land. Thus, land and 
topography are specified by 
a stairstep approximation. The other mask arrays $IU$, $IV$ and $IW$ are derived 
from IN. All mask arrays are INTEGERS*2 and 3-dimensional. 

A row (column) of scalar ghost zones is included around the logically 
rectangular perimeter of the modeled region. These zones are used to store 
temperature, salinity, and momentum data that can be advected into the modeled 
region through inflow cell faces. The data is often specified from climatology, 
but other specifications can be used depending on application. 

At interior points adjacent to boundaries, normal open-boundary fluxes are 
specified by an upwind approximation at the boundary. For inward flow, the 
aforementioned ghost zone data is advected inward. For outward flow, the 
interior data is advected outward. The normal velocity on the boundary is 
updated separately, by: 

a) combining psuedo-Neuman (for outward flow) and climatological restoring 
effects or 

b) by near-boundary data from an overlap region of a separate model grid 
(as in nesting). 

In the former case, if the first interior point has outward directed normal 
velocity component, the boundary point may be nudged toward the nearest interior 
value on a short (synoptic) time scale rather than setting it directly to Neuman 
(which would imply infinite signal propagation speed). 

This simple approach works well, although a more accurate scheme might be 
better. The normal flow on the logically rectangular boundary may be 
continuously relaxed toward a climatologically specified normal velocity on a 
longer time scale. The climatologically specified baroclinic component of the 
normal velocity should be consistent with a geostrophic thermal wind based on 
the climatological salinity and temperature data except very close to the 
equator. The barotropic component is independent of this data, but the total 
velocity must be such that there is no net volume inflow from lateral boundaries 
because of the rigid-lid approximation. 

\section{Governing Equations}

The governing equations of the TIMCOM model are based on
the three-dimensional (3D) primitive equations for an incompressible, stratified fluid.
Both non-hydrostatic and hydrostatic versions are supported. Here, we only describe
the hydrostatic model using the rigid-lid and Boussinesq approximations.
The governing equations of TIMCOM in the spherical coordinate 
are listed as follows. $\lambda$ is the longitudinal variable, 
$\phi$ is the latitudinal variable,
and $z$ is the vertical variable, respectively.
\begin{itemize}
\item Conservation of mass:
\begin{eqnarray}
\frac{1}{R\cos \phi}\left(
\frac{\pd u}{\pd \lambda}
+\frac{\pd (v\cos \phi)}{\pd \phi}
\right)
+ \frac{\partial w}{\partial z}
= 0, \label{CONTI}
\end{eqnarray}
\item Horizontal momentun equations:
\begin{eqnarray}
\frac{\partial u}{\partial t}=
-\mathcal{L} u + \left(f + \frac{u \tan \phi}{R} \right) v 
- \frac{1}{\rho_{0}R \cos \phi}\frac{\partial p}{\partial \lambda}
+ \mathcal{D}_{m} u +  
\frac{\pd }{\pd z} \left(
A_{u} \frac{\pd u}{\pd z}
\right) \label{UME}
\notag\\\\
\frac{\partial v}{\partial t}=
-\mathcal{L}v - \left(
f + \frac{u \tan \phi}{R} 
\right) u
-\frac{1}{\rho_{0}R}\frac{\partial p}{\partial \phi}
+ \mathcal{D}_{m} v 
+ \frac{\pd}{\pd z}
\left(
A_{v}\frac{\pd v}{\pd z}
\right) , \label{VME}
\notag\\
\end{eqnarray}
\item Conservation of potential temperature 
\begin{eqnarray}
\frac{\pd T}{\pd t}=
-\mathcal{L} T + \mathcal{D}_{h} T 
+ \frac{\pd}{\pd z} \left(
K_{T} \frac{\pd T}{\pd z}
\right),
\label{RTMEQ}
\end{eqnarray}
\item Conservation of Salinity 
\begin{eqnarray}
\frac{\pd S}{\pd t}=
-\mathcal{L} S + \mathcal{D}_{h} S
+ \frac{\pd}{\pd z} \left(
K_{S} \frac{\pd S}{\pd z}
\right),
\label{RSMEQ}
\end{eqnarray}
\item Hydrostatic equation:
\begin{eqnarray}
\frac{\partial p}{\partial z}=-(\rho - \bar \rho)g, \label{HYDRO}
\end{eqnarray}
\item Equation of state:
\begin{eqnarray}
\rho=\rho(S,T,p),
\end{eqnarray}
\end{itemize}
where 
$u$ and $v$ are the velocity components in $\lambda$ and $\phi$ directions, respectively,
and $w$ is vertical velocity.
$S$ is the salinity, and $T$ is the potential temperature.
$f$ is Coriolis parameter; $\rho_{0}$ is the mean density;
$\bar \rho$ is the horizontal average of density
at depth $z$; $p$ is the pressure. 
$A_{u}$ and $A_{v}$ are vertical
eddy viscosity coefficients for horizontal momentum equations.
$K_{S}$ and $K_{T}$ are the vertical eddy diffusivity for
salinity and temperature equations, respectively.
The convection and horizontal diffusion operators $\mathcal{L}, \mathcal{D}_{m(h)}$ are defined as:
\begin{eqnarray}
\mathcal{L} &=& \frac{u}{R\cos \phi}
\frac{\pd}{\pd \lambda} + \frac{v}{R}
\frac{\pd}{\pd \phi} + w\frac{\pd}{\pd z}, \\
D_{m(h)} &=&
\frac{A_{m(h)}}{R^{2}}
\left(
\frac{1}{\cos^{2} \phi} \frac{\pd^{2}}{\pd \lambda^{2}}
-\tan \phi \frac{\pd}{\pd \phi}
+\frac{\pd^{2}}{\pd \phi^{2}}
\right),  \label{DIFF}
\end{eqnarray}
where $A_{m(h)}$ is the horizontal eddy viscosity(diffusivity) coefficient 
for momentum(or scalar) equations.

Assuming that the hydrostatic approximation is valid, the pressure field can
be computing by integrating (\ref{HYDRO})
 from rigid-lid surface to any depth $z$:
\begin{eqnarray}
p=p_{s}+p_{b}, \hspace{0.5cm}  
p_{b}=g \int_{z}^{0} \rho dz,
\end{eqnarray}
where $p_{s}$ is the surface pressure due to atmosphere 
or other physical phenomena.

\section{Temporal Discretization}

The three-dimensional, primitive equations are solved on a combined
Arakawa A and C grid.
A variant of projection method (pressure correction) is used here, which splits the numerical operators
and enforces continuity by solving a pressure equation.
Thus, the explicit terms are discretized with either Euler forward (diffusion terms) scheme
or Leap-Frog time stepping (advection and pressure variations) scheme.
The implicit term (Coriolis) is discretized with  Crank-Nicolson scheme.
The time discretized equations (\ref{UME}), (\ref{VME}), (\ref{RTMEQ})
and (\ref{RSMEQ}) can be written as below:
\begin{eqnarray}
\frac{u^{n+1}-u^{n-1}}{2\Delta t} &=&  
\mathcal{D}_{m} u^{n} - \mathcal{L}u^{n} 
+ (\frac{u^{n} \tan \phi}{R} 
+ f)v^{n+1} \newline 
-\frac{1}{\rho_{0}R\cos \phi}\frac{\pd p^{n}}{\pd \lambda} 
+ \frac{\pd}{\pd z} \left( 
A_{u}\frac{\pd u^{n+1}}{\pd z} \right), \notag \label{DMEQU} \\
\frac{v^{n+1}-v^{n-1}}{2\Delta t} &=& 
\mathcal{D}_{m} v^{n} - \mathcal{L}v^{n}
- (\frac{u^{n} \tan \phi}{R} 
+ f)u^{n+1} 
-\frac{1}{\rho_{0}R}\frac{\pd p^{n}}{\pd \phi}
+ \frac{\pd}{\pd z} \left(
A_{v}\frac{\pd v^{n+1}}{\pd z} \right),  \label{DMEQV} \\
\frac{T^{n+1}-T^{n-1}}{2\Delta t} &=& 
\mathcal{D}_{h} T^{n} -\mathcal{L}T^{n}+
\frac{\pd}{\pd z} \left( K_{T}
\frac{\pd T^{n+1}}{\pd z}
\right), \label{DMEQT}\\
\frac{S^{n+1}-S^{n-1}}{2\Delta t} &=&
\mathcal{D}_{h} S^{n} -\mathcal{L}S^{n}+
\frac{\pd}{\pd z} \left( K_{S}
\frac{\pd S^{n+1}}{\pd z}
\right), \label{DMEQS}
\end{eqnarray}

%Keep the following for future use
%When the bottom Ekman layer is resolved, an extremely short time step would be required
%if an explicit scheme were used for the vertical diffusion terms.
%To avoid this restriction, the vertical mixing should be treated implicitly.
Since the Coriolis term and vertical diffusion term 
are updated implicitly, the state variables cannot be updated
in the straightworward manner. 
The other explicit terms are updated initially, including
horizontal diffusion, 
advection and pressure terms for momentum equations (and scalar equaitons).
This leads to the intermediate variables
$u^{*},v^{*}$, $T^{*}$, and $S^{*}$ as follows,
\begin{eqnarray}
\frac{u^{*}-u^{n-1}}{2\Delta t} &=&
\mathcal{D}_{m} u^{n} -\mathcal{L}u^{n}
-\frac{1}{\rho_{0}R\cos \phi}
\frac{\pd p^{n}}{\pd \lambda}, \label{EQMU}\\
\frac{v^{*}-v^{n-1}}{2\Delta t} &=&
\mathcal{D}_{m} v^{n} -\mathcal{L}v^{n}
-\frac{1}{\rho_{0}R}
\frac{\pd p^{n}}{\pd \phi} , \label{EQMV}\\
\frac{T^{*}-T^{n-1}}{2\Delta t} &=&
\mathcal{D}_{h}T^{n}-\mathcal{L}T^{n}, \label{EQMT}\\
\frac{S^{*}-S^{n-1}}{2\Delta t} &=&
\mathcal{D}_{h}S^{n}-\mathcal{L}S^{n}.\label{EQMS}
\end{eqnarray}
Since the right hand side of equations (\ref{EQMU}) to (\ref{EQMS})
are expressed in terms of previously determined state variables,
$*$ intermediate state variables can be updated.

The implicit terms are then used to update 
horizontal vertical diffusion and Coriolis terms.
The variables
$u^{**}$, $v^{**}$, $T^{**}$, and $S^{**}$ are asked to satisfy,
\begin{eqnarray} 
\frac{u^{**}-u^{*}}{2\Delta t}
&=&
\frac{\pd }{\pd z} \left(
A_{u}\frac{\pd u^{**}}{\pd z}
\right), \label{VDU}\\
\frac{v^{**}-v^{*}}{2\Delta t}
&=&
\frac{\pd }{\pd z} \left(
A_{v}\frac{\pd v^{**}}{\pd z}
\right), \\
\frac{T^{**}-T^{*}}{2\Delta t}
&=&
\frac{\pd }{\pd z} \left(
K_{T}\frac{\pd T^{**}}{\pd z}
\right), \\
\frac{S^{**}-S^{*}}{2\Delta t}
&=&
\frac{\pd }{\pd z} \left(
K_{S}\frac{\pd T^{**}}{\pd z}
\right).
\end{eqnarray}
Details of updations for these $**$ variables
are illuminated in the next section.
After updating vertical diffusion terms, Coriolis term is going to update. 
The detailed trapezoidal procedure is provided in the next section.
At this point, the scalar equations (temperature and salinity) are completely
updated, that is $T^{n+1}=T^{**}$ and $S^{n+1}=S^{**}$,
but the velocity fields are based on the use of previous pressure and
generally do not satisfy the depth integrated continuity equations, 
i.e. the divergence may not be free.
We have equations that is for updation for Coriolis terms,
\begin{eqnarray}
\frac{\tilde{u}^{n+1}-u^{**}}{2\Delta t} = 
\left(
\frac{u^{n}\tan \phi}{R} + f
\right) \tilde{v}^{n+1}, \\
\frac{\tilde{v}^{n+1}-v^{**}}{2\Delta t} =
-\left(
\frac{u^{n}\tan \phi}{R} + f
\right) \tilde{u}^{n+1}
\end{eqnarray}
The last step for the pressure correction method is then to determine the
correction terms for the pressure and horizontal velocity fields.
This ensures the incompressibility is unconditionally hold
so that the appropriate (more accurate) velocity and pressure fields are achieved at time step $n+1$. 
Assuming the final pressure and velocity fields at time step $n+1$ are
$p^{n}=\tilde{p}^{n}+\Delta p$,
$u^{n+1}=\tilde{u}^{n+1}+\Delta u$, and
$v^{n+1}=\tilde{v}^{n+1}+\Delta v$,
we can substitute them into equations 
(\ref{DMEQU}) and (\ref{DMEQV}), respectively, This leads to the following corrections.
\begin{eqnarray}
\Delta u &=& -2\Delta t \times \frac{1}{R \cos \phi}
\frac{\pd \Delta p}{\pd \lambda}, \label{CORRECU}\\
\Delta v&=& -2\Delta t \times
\frac{1}{R}
\frac{\pd \Delta p}{\pd \phi}. \label{CORRECV}
\end{eqnarray}
Note that the correction $\Delta u$ and $\Delta v$ are depth independent.
Note that the surface rigid-lid approximation and no normal flow condition
at the bottom ensures the vertical velocity $w$ as follows,
\begin{eqnarray}
w_{0}=0, \hspace{0.5cm} w_{-h}=0
\end{eqnarray}
where the subscripts $0$ and $-h$ denote 
the sea surface and bottom, respectively.
Integating (\ref{CONTI}) with the boundary conditions, we get
\begin{eqnarray}
\int^{0}_{-h} \frac{1}{R \cos \phi}
\frac{\pd u}{\pd \lambda}
+ \frac{\pd (v\cos \phi)}{\partial \phi} dz
=-w_{0}+w_{-h}=0 \label{INTDIV}
\end{eqnarray}
Using (\ref{INTDIV})  and the relationship
between velocity and pressure corrections in (\ref{CORRECU}) and (\ref{CORRECV}),
we get:
\begin{eqnarray}
0&=&\int^{0}_{-h} \frac{1}{R \cos \phi}
\left(\frac{\pd \tilde{u}^{n+1}}{\pd \lambda}
+\frac{\pd (\tilde{v}^{n+1}\cos \phi)}{\pd \phi}
\right) dz
=-\int^{0}_{-h} \frac{1}{R \cos \phi}
(\frac{\pd \Delta u}{\pd \lambda}
+\frac{\pd (\Delta v\cos \phi)}{\pd \phi}) dz \notag \\
&=&\frac{2\Delta t}{\rho_{0}R^{2} \cos \phi}
\int^{0}_{-h}\left(
\frac{\pd}{\pd \lambda} \frac{1}{\cos \phi}
\frac{\pd \Delta p}{\pd \lambda}
+\frac{\pd}{\pd \phi} \left(
\frac{\pd \Delta p}{\pd \phi}
\right) \cos \phi
\right) dz.
\end{eqnarray}
Rearranging the terms in the above equation results in
\begin{eqnarray}
\frac{2\Delta t}{\rho_{0}R^{2} \cos \phi}
\int^{0}_{-h}\left(
\frac{\pd}{\pd \lambda} \frac{1}{\cos \phi}
\frac{\pd \Delta p}{\pd \lambda}
+\frac{\pd}{\pd \phi} \left(
\frac{\pd \Delta p}{\pd \phi}
\right) \cos \phi
\right) dz = 0.
\end{eqnarray}
Solving $\Delta p$ is equivalent to solving a 
2-D pressure Poisson equation for the surface pressure correction, where the forcing is
determined by the net divergence implied by the final updated velocity fields
over each vertical water column.
An efficient EVP solver is applied for .
The $\Delta p$ represents the required surface pressure change to adjust
the vertical velocity at the bottom to be zero, as defined by the non-slip bottom boundary.
Zero normal velocity is also specified at the bottom boundary condition.
After the $\Delta p$ is obtained,  the pressure and horizontal velocity fields
at the next time step $n$ are updated. 
Finally, the hydrostatic vertical velocity $\hat{w}^{n+1}$  can be calculated diagnostically
by equation (\ref{CONTI}), indicating the divergence-free property is always hold.
The whole state variables are now fully updated.

Let's summarize the updates of these velocity variables:
\begin{itemize}  
\item $u^{*}, v^{*}$ 
are the variables which update explicitly without Coriolis terms and curvature effects,
as well as vertical diffusion terms. 
They are determined by equations (\ref{EQMU}) and (\ref{EQMV}).
We might call this updation explicit updation.
\item $u^{**}, v^{**}$ and $\tilde{u}^{n+1}, \tilde{v}^{n+1}$  
are the updated variables that including the implicit contribution 
of vertical diffusion.
They are determined in the following section.
We can call these two updation steps implicit updation
\item $\Delta{u}, \Delta{v}$ are the velocity correction terms,
determined by equations (\ref{CORRECU}) and (\ref{CORRECV}) 
after computing the pressure correction term $\Delta{p}$ by an efficient EVP solver.
\end{itemize}
After these updation, we can write 
\begin{eqnarray}
u^{n+1} = \tilde{u}^{n+1} + \Delta u, \\
v^{n+1} = \tilde{v}^{n+1} + \Delta v,
\end{eqnarray}
which are final velocity fields.

\subsection{Treatment of the Coriolis Term}

Coriolis and curvature terms 
requires a careful treatment since numerical instability may result.
We adopted trapezoidal treatment for Coriolis term, 
which discusses detailly in the following section.
%There are several methods for treating Coriolis term such as backward or
%forward treatment.
Here, we use trapezoidal treatment for Coriolis term. 
For convinience, let $f+\frac{u^{n}\tan \phi}{R}\equiv \tilde{F}$.
Consider the below momentum equations:
\begin{eqnarray}
\tilde{u}^{n+1} &=& u^{**}+
2\Delta t \tilde{F}
(\frac{\tilde{v}^{n+1}+v^{n-1}}{2}) 
= u^{**} + 
\tilde{F}\Delta t  v^{n-1} + 
\tilde{F}\Delta t \tilde{v}^{n+1} \label{FU}\\
\tilde{v}^{n+1} &=& v^{**}-
2\Delta t \tilde{F}
(\frac{\tilde{u}^{n+1}+u^{n-1}}{2}) 
= v^{**} - 
\tilde{F} \Delta t u^{n-1} - 
\tilde{F} \Delta t  \tilde{u}^{n+1} \label{FV},
\end{eqnarray}
where $f$ is Coriolis parameter.\\
We treat $n-1$ terms by using previus step values, and treat $n+1$ terms
by using values that is the newest updating values. So, we can separate
them into two updation steps.\\
In first step, consider equations associated to $**$ intermediate step terms,
\begin{eqnarray}
\tilde{u}^{n+1}_{int} = u^{**} + 
\tilde{F} \Delta t \times v^{n-1}, \\
\tilde{v}^{n+1}_{int} = v^{**} - 
\tilde{F} \Delta t \times u^{n-1},
\end{eqnarray}
where $int$ denotes intermediate variable. For simlicity, we define
$QU \equiv \tilde{u}^{n+1}_{int}$ and $QV \equiv \tilde{v}^{n+1}_{int}$.\\
After compute $QU$ and $QV$, we can proceed second updation step.
We want to update second step by using newest velocity variables, so 
we can write down equations:
\begin{eqnarray}
\tilde{u}^{n+1} = \tilde{u}_{int}^{n+1}+ 
\tilde{F} \Delta t \times \tilde{v}^{n+1}, \label{FU2}\\
\tilde{v}^{n+1} = \tilde{v}_{int}^{n+1}- 
\tilde{F} \Delta t \times \tilde{u}^{n+1}, \label{FV2}
\end{eqnarray}
Substitute (\ref{FV2}) into (\ref{FU2}) we can get
\begin{eqnarray}
\tilde{u}^{n+1}=\frac{1}{1+(\tilde{F}\Delta t)^{2}}
[\tilde{u}^{n+1}_{int}+\tilde{F} 
\Delta t \times \tilde{v}_{int}^{n+1}]
=\frac{1}{1+(\tilde{F}\Delta t)^{2}}(QU + 
\tilde{F}\Delta t \times QV).
\end{eqnarray}
Similarly, substitute (\ref{FU2}) into (\ref{FV2}), we can get
\begin{eqnarray}
\tilde{v}^{n+1}=\frac{1}{1+(\tilde{F}\Delta t)^{2}}
[\tilde{v}^{n+1}_{int}-\tilde{F} \Delta t \times \tilde{u}_{int}^{n+1}]
=\frac{1}{1+(\tilde{F}\Delta t)^{2}}(QV - \tilde{F}\Delta t \times QU)
\end{eqnarray}
Right hand side of these equations are all known variables, thus we can 
update velocity variables including Coriolis terms.\\
Implicit treatment method for Coriolis terms gives 
such as trapezoindal method
better performance of velocity fields, it might reduce oscillation results, 
but numerical experiments need to be operated in order to 
get appropriate choice for Coriolis terms.


\subsection{Implicit Treatment for Vertical Diffusion}

Concerning about stability, 
the most sensitive term in momentum equations is
diffusion term, especially vertical diffusion term,
so we might treat vertical diffusion terms implicitly 
in order to get stability numerical scheme.
In this section, we consider only momentum equation of $u$, 
that is (\ref{UME}), 
other governing equations can be done by the similar method.
Consider equation (\ref{VDU}) again,
\begin{eqnarray}
\frac{u^{**}-u^{*}}{2\Delta t}
=
\frac{\pd }{\pd z}
\left(
A_{u}\frac{\pd u^{**}}{\pd z}
\right)
\end{eqnarray}
Write discretized form of equation (\ref{VDU}) in space
and after some arrangement, we can get equation
\begin{eqnarray}
a_{i,j,k}^{n+1}u^{**}_{i,j,k+1} +
b_{i,j,k}^{n+1}u^{**}_{i,j,k} +
c_{i,j,k}^{n+1}u^{**}_{i,j,k-1} = 
 u^{*}_{i,j,k-1}.  \label{IMPU}
\end{eqnarray}
where 
\begin{eqnarray}
a_{i,j,k} &\equiv &
-\frac{EV_{i,j,k}}{\Delta z_{k}\Delta z_{k-1}} 
\times 2\Delta t,  \notag \\
b_{i,j,k} &\equiv &
2\Delta t \times \left(
\frac{EV_{i,j,k}}{\Delta z_{k}\Delta z_{k-1}}
+\frac{EV_{i,j,k-1}}{\Delta z_{k-1}^{2}} \right) \notag \\
c_{i,j,k} &\equiv &
1-\frac{EV_{i,j,k-1}}{\Delta z_{k-1}^{2}}\times 2\Delta t.
\end{eqnarray}
Note that $k$ is form $1$ to $K1$, 
$1$ means sea surface layer,
$K1$ means bottom layer.
$k=0$ and $k=K0$ variables are ghost zone variables\\
Insert $k$ from $1$ to $K1$ into equation (\ref{IMPU}), 
we can formulate a system of equations of number $K1$,
write it as matrix form $Mx=G$:
\begin{eqnarray}
\left[
\begin{array}{cccc}
b_{1} & a_{1} &  \cdots & 0 \\
c_{2} & b_{2} & a_{2} & \cdots \\
\multicolumn{4}{c}\dotfill\\
0 & \cdots & c_{K1} & b_{K1}
\end{array}
\right]
\left[
\begin{array}{c}
u_{1}^{**} \\
u_{2}^{**} \\
\vdots      \\
u_{K1}^{**}
\end{array}
\right]
=
\left[
\begin{array}{c}
u_{0}^{*} \\
u_{1}^{*} \\
\vdots      \\
u_{K1-1}^{*}
\end{array}
\right]
+
\left[
\begin{array}{c}
-c_{1}{u}_{0}^{**} \\
0           \\
\vdots      \\
-a_{K1}u^{**}_{K0}
\end{array}
\right],
\end{eqnarray}
where $M$ is a tridiagonal square matrix.
For simplicity we use subscript number to denote 
which vertical layer it is.
Subscript $0$ is the ghost zone variable,
which is used for implementation of boundary condition.
Since matrix $M$ is a tridiagonal matrix with three stencil,
we can use Thomas algorithm to solve this system,
which is an efficient solver to solve tridiagonal
matrix system.

After solving $u^{**}$, we can proceed to do Coriolis updation
as previous section illustrated.

\subsection{Time Filter}
 
Time filter is used for stability control 
when we solve a time-dependent problem.
Since we are using Leap-frog time advancement structure as a basis,
filter used here needs to concern Leap-frog scheme.
So, after updating variables in the previous section,
we operate filter scheme for time advacement.
Consider now we are in the $n$ step, variables at time step $n-1$ are
what we already have from the previous updating.
The goal is to gain variables after filering at $n$ step.
By classical Leap-frog sheme, we have:
\begin{eqnarray}
u^{n+1} - u^{n-1} = 2\Delta t \times F(u^{n},t_{n}),
\end{eqnarray}
where $F$ is some known source term.
Insert filter into the scheme, we have:
\begin{eqnarray}
u^{n}_{*} = OFLTW \times u^{n} + FLTW \times (u^{n-1}+u^{n+1}),
\end{eqnarray}
where $*$ means variables that have been filtered 
at time step $n$, and $OFLTW=0.9$, $FLTW=0.05$ are set in the model.
In DieCast ocean model, we accomplish this by
assigning variables $U1$ and $ULF$ cyclingly,
where $U1$, $ULF$, $U2$ denote old, filtered (in the central time step),
and new variables respectively.
Convergent oreder of this scheme is first order.
Convergent order analysis by Taylar expansion is shown below:
\begin{eqnarray*}
\frac{\Delta t}{20} \times u^{''n}.
\end{eqnarray*}

The Leap-Frog could also be replaced by a third-order Adam-Bashford method soon. 
The general formula of leap-frog is shown as follows,
\begin{eqnarray}
\frac{dF}{dt}=S(F(x,t),x,t),
\end{eqnarray}
where right-hand side $S$ is a specified source term (or partial differential equations),
$x$($F$) is the variable (or function) that we are interested in. 
The discretized equation from time step $t_{n}$ to the next step $t_{n+1}$
(superscripts of the variables) is written as:
\begin{eqnarray}
\frac{F^{n+1}_{i} - F^{n-1}_{i}}{2\Delta t}
= S(F(x_{i},t_{n}),x_{i},t_{n}),
\end{eqnarray}
where the superscripts $n-1$ and $n+1$ are the actual time step and
the subscript $i$ denotes the spatial discretization.


\section{Spatial Discretization}

In this section, we take look at how to solve every terms in momentum
equations in space, that is how to slove $\frac{\pd p}{\pd x}$,
$L$ operator, dissipation terms etc.
We use notations such as $\Delta x_{j}$, $\Delta y_{j}$ for
increments in $\lambda$ and $\phi$ direction, subscript $j$
denotes increments that includes curvilture effect.
\subsection{Pressure Gradient}
After integrating hydrostatic equation (\ref{HYDRO}), we get pressure
field $p$ in the whole domain. Then we compute pressure gradient term
by using fourth-order interpolation scheme:
\begin{eqnarray}
\left(\frac{\pd p^{n}}{\pd x}\right)_{i,j,k}
=\frac{7}{12} 
\left(\frac{p_{i+1,j,k}^{n}-p_{i,j,k}^{n}}{\Delta x_{j}}
+\frac{p_{i,j,k}^{n}-p_{i-1,j,k}^{n}}{\Delta x_{j}}   
\right)-
\frac{1}{12} \left(
\frac{p_{i+2,j,k}^{n}-p_{i+1,j,k}^{n}}{\Delta x_{j}}
+\frac{p_{i-1,j,k}^{n}-p_{i-2,j,k}^{n}}{\Delta x_{j}}
\right),
\end{eqnarray}
with truncation error. \\

\subsection{Fluxes}
$\mathcal{L}$ operator as before described can be discretized as the form:
\begin{eqnarray}
\mathcal{L} q^{n}_{i,j,k}&=&(
\frac{U^{n}_{i,j,k}q^{n}_{i+1/2,j,k}-
U^{n}_{i-1,j,k}q^{n}_{i-1/2,j,k}}{\Delta x_{j}}
+
\frac{V^{n}_{i,j,k}q^{n}_{i,j+1/2,k}-
V^{n}_{i,j-1,k}q^{n}_{i,j-1/2,k}}{\Delta y_{j}} \notag \\
&+&
\frac{W^{n}_{i,j,k+1}q^{n}_{i,j,k+1/2}-
W^{n}_{i,j,k-1}q^{n}_{i,j,k-1/2}}{\Delta z_{k}}
),
\end{eqnarray}
where $q$ denotes some variable of $u,v$ or $T$, subscripts $+1/2$ and $-1/2$
means face quantities of the variable. So, in order to get
face quantities, we use the fourth-order interpolation scheme
as used for deriving pressure gradient in the previous section:
\begin{eqnarray}
q_{i+1/2}=\frac{7}{12}(Q_{i,j,k}+Q_{i+1,j,k})-
\frac{1}{12}(Q_{i-1,j,k}+Q_{i+2,j,k}),
\end{eqnarray}
where uppercase $Q$ means cell-quantity, while lowercase $q$ means
face-quantity.\\
For determineing vertical flux, first we interpolate cell-quantities
in order to get face-quantities.
\begin{eqnarray}
u_{i,j,k+1/2}^{n}=\frac{7}{12}(U_{i,j,k}^{n}+U_{i,j,k+1}^{n})-
\frac{1}{12}(U_{i,j,k-1}^{n}+U_{i,j,k+2}^{n}), \label{FDU}
\end{eqnarray}
where $n$ demotes at $n$ time step.
For stability and better oscillation control, experiences gave
modified equation:
\begin{eqnarray}
u_{i,j,k+1/2}^{n}=\frac{7}{12}(U_{i,j,k}^{n}+U_{i,j,k+1}^{n})-
\frac{1}{12}(U_{i,j,k-1}^{*}+U_{i,j,k+2}^{n}),
\end{eqnarray}
$*$ denotes intermediate step quantity between time step $n$ and $n+1$.\\
In DieCast ocean model, we use the original equation (\ref{FDU}).
Then multiply face-quantity $u_{i,j,k+1/2}^{n}$ with
cell-quantity $U_{i,j,k+1}^{n}$
($u_{i,j,k-1/2}^{n}$ with
cell-quantity $U_{i,j,k}^{n}$),
we can get flux quantitities at verticl direction $Wu_{i,j,k}^{n}$
before computing derivative.\\
Use simple finite difference shceme, we can get vertical flux:
\begin{eqnarray}
F_{Z}q_{i,j,k}^{n}=\frac{W^{n}_{i,j,k+1}q^{n}_{i,j,k+1/2}-
W^{n}_{i,j,k-1}q^{n}_{i,j,k-1/2}}{\Delta z_{k}},
\end{eqnarray}
where $q$ can be horizontal velocity variables $u,v$ or temperature $T$.
$F_{Z}q^{n}_{i,j,k}$ is the approximation to vertical flux 
$\left(w\frac{\pd q}{\pd z}\right)^{n}_{i,j,k}$.\\
Similarly, `Longitudinal Flux' and `Latitudinal Flux' are determined by
\begin{eqnarray}
F_{X}q_{i,j,k}^{n}=\frac{U^{n}_{i,j,k+1}q^{n}_{i,j,k+1/2}-
U^{n}_{i,j,k-1}q^{n}_{i,j,k-1/2}}{\Delta x_{j}}, \\
F_{Y}q_{i,j,k}^{n}=\frac{V^{n}_{i,j,k+1}q^{n}_{i,j,k+1/2}-
V^{n}_{i,j,k-1}q^{n}_{i,j,k-1/2}}{\Delta y_{j}}.
\end{eqnarray}
So, $\mathcal{L}$ operator terms in equations 
(\ref{DMEQU}), (\ref{DMEQV}), and (\ref{DMEQT})
are to be determined as:
\begin{eqnarray}
\mathcal{L} u^{n}_{i,j,k}= F_{X}u^{n}_{i,j,k}+F_{Y}u^{n}_{i,j,k}+F_{Z}u^{n}_{i,j,k}, \\
\mathcal{L} v^{n}_{i,j,k}= F_{X}v^{n}_{i,j,k}+F_{Y}v^{n}_{i,j,k}+F_{Z}v^{n}_{i,j,k}, \\
\mathcal{L} T^{n}_{i,j,k}= F_{X}T^{n}_{i,j,k}+F_{Y}T^{n}_{i,j,k}+F_{Z}T^{n}_{i,j,k}, \\
\mathcal{L} S^{n}_{i,j,k}= F_{X}S^{n}_{i,j,k}+F_{Y}S^{n}_{i,j,k}+F_{Z}S^{n}_{i,j,k},
\end{eqnarray}


\subsection{Dissipation and Diffusion Terms}

\subsubsection{Vertical Diffusion}
We separate diffusion terms into vertical diffusion contribution
and horizontal diffusion contribution.
Parameters $EV$ and $HV$ are vertical diffusion
coefficiets for $Km$ and $K_{h}$
in equations (\ref{DMEQU}), (\ref{DMEQV}) and (\ref{DMEQT}) respectively.
We compute them by using previous time step $n-1$ variables.\\
First compute 
\begin{eqnarray}
\left(A_{u}\frac{\pd u^{n-1}}{\pd z}\right) 
\approx EV_{i,j,k}\times 
\left(\frac{u^{n-1}_{i,j,k+1}-u^{n-1}_{i,j,k}}{\Delta z_{k}}\right)
\equiv D_{Z}u^{n}_{i,j,k}, \\
\left(A_{v}\frac{\pd v^{n-1}}{\pd z}\right)
\approx EV_{i,j,k}\times 
\left(\frac{v^{n-1}_{i,j,k+1}-v^{n-1}_{i,j,k}}{\Delta z_{k}}\right)
\equiv D_{Z}v^{n}_{i,j,k}, \\
\left(K_{T}\frac{\pd T^{n-1}}{\pd z}\right)
\approx HV_{i,j,k}\times 
\left(\frac{T^{n-1}_{i,j,k+1}-T^{n-1}_{i,j,k}}{\Delta z_{k}}\right)
\equiv D_{Z}T^{n}_{i,j,k}, \\
\left(K_{S}\frac{\pd S^{n-1}}{\pd z}\right)
\approx HV_{i,j,k}\times
\left(\frac{S^{n-1}_{i,j,k+1}-S^{n-1}_{i,j,k}}{\Delta z_{k}}\right)
\equiv D_{Z}S^{n}_{i,j,k}. 
\end{eqnarray}
Then compute
\begin{eqnarray}
\frac{\pd}{\pd z} \left(
A_{u}\frac{\pd u^{n-1}}{\pd z} \right) \approx
\left( \frac{D_{Z}u_{i,j,k+1}^{n}-D_{Z}u_{i,j,k}^{n}}{\Delta z_{k}}
\right), \\
\frac{\pd }{\pd z} \left(
A_{v}\frac{\pd v^{n-1}}{\pd z} \right) \approx
\left( \frac{D_{Z}v_{i,j,k+1}^{n}-D_{Z}v_{i,j,k}^{n}}{\Delta z_{k}}
\right), \\
\frac{\pd }{\pd z} \left(
K_{T}\frac{\pd T^{n-1}}{\pd z} \right) \approx
 \left( \frac{D_{Z}T_{i,j,k+1}^{n}-D_{Z}T_{i,j,k}^{n}}{\Delta z_{k}}
\right), \\
\frac{\pd }{\pd z} \left(
K_{S}\frac{\pd S^{n-1}}{\pd z} \right) \approx
 \left( \frac{D_{Z}S_{i,j,k+1}^{n}-D_{Z}S_{i,j,k}^{n}}{\Delta z_{k}}
\right).
\end{eqnarray}

\subsubsection{Horizontal Diffusion}
Suppose $A_{m}$ and $A_{h}$ are constants in equations
(\ref{DIFF}).
As similar as previous section, we separate diffusion
%\mathcal{D}_{m}% and %\mathcal{D}_{h}% into several parts.
First we compute x-component quantities
\begin{eqnarray}
\left(\mathcal{D}_{m}\frac{\pd u^{n-1}}{\pd x}\right)
\approx DMX\times
\left(\frac{u^{n-1}_{i+1,j,k}-u^{n-1}_{i,j,k}}{\Delta x_{j}}\right)
\equiv D_{X}u^{n}_{i,j,k}, \\
\left(\mathcal{D}_{m}\frac{\pd v^{n-1}}{\pd x}\right)
\approx DMX\times
\left(\frac{v^{n-1}_{i+1,j,k}-v^{n-1}_{i,j,k}}{\Delta x_{j}}\right)
\equiv D_{X}v^{n}_{i,j,k}, \\
\left(\mathcal{D}_{h}\frac{\pd T^{n-1}}{\pd x}\right)
\approx DHX\times
\left(\frac{T^{n-1}_{i+1,j,k}-T^{n-1}_{i,j,k}}{\Delta x_{j}}\right)
\equiv D_{X}T^{n}_{i,j,k}, \\
\left(\mathcal{D}_{h}\frac{\pd S^{n-1}}{\pd x}\right)
\approx DHX\times
\left(\frac{S^{n-1}_{i+1,j,k}-S^{n-1}_{i,j,k}}{\Delta x_{j}}\right)
\equiv D_{X}S^{n}_{i,j,k}.
\end{eqnarray}

Then we compute y-component quantities
\begin{eqnarray}
\left(\mathcal{D}_{m}\frac{\pd u^{n-1}}{\pd y}\right)
\approx DMY\times
\left(\frac{u^{n-1}_{i,j+1,k}-u^{n-1}_{i,j,k}}{\Delta y_{j}}\right)
\equiv D_{Y}u^{n}_{i,j,k}, \\
\left(\mathcal{D}_{m}\frac{\pd v^{n-1}}{\pd y}\right)
\approx DMY\times
\left(\frac{v^{n-1}_{i,j+1,k}-v^{n-1}_{i,j,k}}{\Delta y_{j}}\right)
\equiv D_{Y}v^{n}_{i,j,k}, \\
\left(D_{h}\frac{\pd T^{n-1}}{\pd y}\right)
\approx DHY\times
\left(\frac{T^{n-1}_{i,j+1,k}-T^{n-1}_{i,j,k}}{\Delta y_{j}}\right)
\equiv D_{Y}T^{n}_{i,j,k}, \\
\left(D_{h}\frac{\pd S^{n-1}}{\pd y}\right)
\approx DHY\times
\left(\frac{S^{n-1}_{i,j+1,k}-S^{n-1}_{i,j,k}}{\Delta y_{j}}\right)
\equiv D_{Y}S^{n}_{i,j,k}. 
\end{eqnarray}

And then
\begin{eqnarray}
D_{m} u^{n-1}_{i,j,k}
=\frac{D_{X}u^{n}_{i,j,k}-D_{X}u^{n}_{i-1,j,k}}{\Delta x_{j}}
+\frac{D_{Y}u^{n}_{i,j,k}-D_{Y}u^{n}_{i,j-1,k}}{\Delta y_{j}}, \\
D_{m} v^{n-1}_{i,j,k}
=\frac{D_{X}v^{n}_{i,j,k}-D_{X}v^{n}_{i-1,j,k}}{\Delta x_{j}}
+\frac{D_{Y}v^{n}_{i,j,k}-D_{Y}v^{n}_{i,j-1,k}}{\Delta y_{j}}, \\
D_{h} T^{n-1}_{i,j,k}
=\frac{D_{X}T^{n}_{i,j,k}-D_{X}T^{n}_{i-1,j,k}}{\Delta x_{j}}
+\frac{D_{Y}T^{n}_{i,j,k}-D_{Y}T^{n}_{i,j-1,k}}{\Delta y_{j}}, \\
D_{h} S^{n-1}_{i,j,k}
=\frac{D_{X}S^{n}_{i,j,k}-D_{X}S^{n}_{i-1,j,k}}{\Delta x_{j}}
+\frac{D_{Y}S^{n}_{i,j,k}-D_{Y}S^{n}_{i,j-1,k}}{\Delta y_{j}}. 
\end{eqnarray}

RHS of equations (\ref{DMEQU}), (\ref{DMEQV}), and (\ref{DMEQT})
are all be determined at time step $n$.


\section{EVP Poisson Solver}
This might be corrected when we use CCD?

\section{Boundary Conditions}

Boundary conditions are required to close the governing equations.
There are several boundary options for the lateral boundary conditions in TIMCOM.
For a close boundary condition, no normal flow across solid boundary is imposed.
Either non-slip or free-slip is satisfied. In the case of non-slip boundary condition,
the tangential velocity at the boundary is set to zero, while the tangential stress
is set to zero in the case of free-slip boundary condition.

At present, one of the greatest limitations in regional modeling is the
open boundary conditions (OBCs). 
The open boundary should allow perturbations generated inside the
computational domain to leave it without deterioration of the inner model
solution and allow physically important external information to advect inward.
The rigid-lid approximation simplifies the treatment of open boundaries
by greatly reducing the range of frequencies that must be addressed.
The Orlanski-like radiative approximation is excellent when dealing with
a single dominant phase speed~\cite{Orlanski76}.
The Orlanski radiation condition solves a one-dimensional wave
equation at the boundary, thus allowing waves to propagate out of the domain,
but it is not accurate when a wide range of phase
speeds must be addressed.
A simple and effective upwind open boundary option is provided in TIMCOM,
which works nicely in several applications~\cite{Tseng05,Tseng06,Tseng09}.

All open boundary conditions are based on boundary fluxes.
A pure upwind advective scheme is used at the three lateral open boundaries
(north, south, and west) for all variables. The normal boundary advective
flux is expressed as
\begin{equation}
F_{n}=U_{n}\frac{\partial \phi}{\partial n}
\end{equation}
where
\begin{equation}
\frac{\partial \phi}{\partial n}=
\begin{cases}
(\phi-\phi_{o})/\Delta x_{n} & U_{n} \ge 0 \\
(\phi_{i}-\phi)/\Delta x_{n} & U_{n} < 0
\end{cases}
\end{equation}
and $U_{n}$ is the normal velocity on the open boundary.
$\phi$ represents any of the three velocity components, temperature or salinity
at the boundary. $\phi_{o}$ is the variable on the open boundary
and $\phi_{i}$ is the variable at one grid point inside the open boundary,
$\Delta x_{n}$ is the grid spacing in the direction normal to the boundary.
Thus, large scale data are advected inward at an inflow boundary
and the interior data is advected outward at an outflow boundary.
The net lateral boundary flow through the open boundary is adjusted
at every time step at the open boundary.
This insures there is no net inflow into the modeled region
and is required by incompressibility.
It has been argued that the hydrostatic, primitive equations are
ill-posed when an inappropriate open boundary condition is used.
The ill-posedness implies no single set of boundary conditions
can be defined to guarantee existence of a stable, unique solution.
If the proper number of boundary conditions is not specified,
the solution of the primitive equations will lead to the
exponential growth of energy and numerical instability~\cite{Oliger78}.
Note that no other sponge layers are used, which
dissipate the errors generated by the over-specification of boundary
conditions in most ocean models.

The default bottom topography in TIMCOM is the stair-step representation.
Vertical velocity $w$ is based on rigid-lid approximation at surface and the condition 
of no normal flow at the bottom by setting:
\begin{eqnarray}
w_{0}=0, \hspace{0.5cm} w_{-h}=0
\end{eqnarray}
where the subscripts $0$ and $-h$ denote 
the sea surface and bottom, respectively.
Since the real bathymetry is always complex, the option of immesred boundary methods
developed by~\cite{Tseng03}can be used to better represent the bottom topography. 

\newpage
