% Governing Equations.

\chapter{Governing Equations and Numerical 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. 

\subsection{The TIMCOM 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{Governing Equations}

Since TIMCOM divides the ocean domain into a large number of cubical finite control 
volumes, or cells, most of them contain water, while others may cover land or 
locate below the bottom of the ocean. To reduce wasted computation in masked land areas, one may use 
domain decomposition with domain coupling via the robust TIMCOM grid coupling 
approach. 

The conservation equations are applied in integral form to each individual 
control volume; this means that time changes are calculated from a sum of flux 
approximations across the six sides of the cubical control volumes plus special 
terms such as pressure gradient and Coriolis. The side faces are exactly 
vertical. The top and bottom faces are exactly horizontal. Bottom and lateral 
boundaries are approximated by a stairstep approximation. The sides are assumed 
to meet at right angles, as in Cartesian and spherical coordinates. We present 
the details of TIMCOM hydrostatic model in this section. The governing equations 
are: \newline 

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}

Vertical momentum equation:
\begin{align}
\fc{\pd w}{\pd t}&= -\nabla \cdot w\mathbf{V}
-\fc{1}{\rho_0}\fc{\pd p}{\pd z}+\nabla_{h} \cdot A_h\nabla_{h}w
+\fc{\pd}{\pd z}(A_v\fc{\pd w}{\pd z})-(\fc{\rho-\overline{\rho}}{\rho_0})g 
\label{eq:3DNS_W}\\
\end{align}

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

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

Equation of state:
\be
\rho=\rho(S,T)
\ee

\begin{table}
\caption{Parameters in control equations}
\centering
\begin{tabular}{|c|l|}
\hline
Parameter	& Description \\
\hline
\hline
$u$ & $x$ components of velocity vector  $V= (u, v, w)$ \\
$v$	& $y$ components of $V= (u, v, w)$ \\
$w$	& $z$ components of $V= (u, v, w)$ \\
$f$	& the Coriolis parameter \\ 
$\rho_0$ & is the background density \\ 
$\rho$ & the horizontally averaged density at depth $z$ \\ 
$p$	& pressure \\
$A_h$ & horizontal eddy viscosity \\ 
$A_v$	& vertical eddy viscosity \\ 
$S$	& salinity \\ 
$T$	& potential temperature \\
$K_h$	& horizontal eddy diffusivity \\
$K_v$ &	vertical eddy diffusivity \\ 
\hline
\end{tabular}
\label{parameters}
\end{table}

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.

\subsection{Discretization}

From the collocated grid arrangement, 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

\noindent 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 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}).   

\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. 

\subsection{Leapfrog Time Stepping}
 
Leap-frog time-stepping uses three time-levels. Consider a differential equation      			 
\begin{equation}
du/dt=F(u,t)
\end{equation}
which can be written using second-order accurate centered time-differencing 
\begin{equation}
U3-U1=2\times \Delta t \times F(U2,t2)
\end{equation} 
so that $U3=U1+2\times \Delta t \times F2$ . To achieve second-order accuracy, we 
clearly need to obtain $U2$ at the central time level. This is done by cycling 
the variables so that $U1=U2$ and $U2=U3$ at the end of the time step. Note that 
the update from the old time level $U1$ leaps over the central level $U2$ to the 
leap-frog level $U3$. Generally, time-splitting of the solution will result and 
this is corrected by adding an Asselin filter~\citep{Asselin72} with weight FLTW 
before 
cycling the variables so $U2=(1-FLTW)\times U2+FLTW\times (U1+U3)$ which mixes 
some of the 
information in from other time levels. This Aselin filter formally reduces the 
time-stepping to first-order accuracy~\citep{Roache88} -- albeit a 
generally accurate first-order, because very small FLTW is adequate to eliminate 
time splitting. 

\subsection{Open Boundary Treatment} 

TIMCOM uses the following open boundary treatment. All boundaries are defined by 
cell faces that divide water-containing cells from land or external fluid. If 
there is external fluid outside the model domain, special open boundary 
treatment is required. 

All model predicted field variables are internal cell averaged quantities by 
definition. In a given time step these model variables are first updated 
excluding effects of open boundary fluxes. This is automatically done by mask 
array factors that effectively give zero normal boundary fluxes on all open 
boundary cell faces, including islands and seamounts as well as the logically 
rectangular domain boundaries. The mask arrays are defined by the model 
preprocessor. 

Full control of (and responsibility for) open boundary conditions is left to the 
user. Such control makes TIMCOM general, efficient and relocatable but requires 
some expertise and care by the user as with any general ocean model. 

The TIMCOM model code includes open boundary effects along all lateral 
boundaries that fall on its logically rectangular grid boundaries. The open 
boundary coding is designed for efficient use of vector and massively parallel 
computers, as is the rest of TIMCOM. 

TIMCOM open boundary conditions are specified or derived fluxes of all 
transported quantities into the adjacent-to-boundary control volumes, as 
required at open boundaries~\citep{Dietrich87}. This approach is simple and 
robust. For example, in nesting a one minute (1/60$^\circ$) resolution TIMCOM 
Santa 
Barbara Channel application inside a coarser five minute (1/12$^\circ$) 
California 
Current `BIG picture' application (developed in collaboration with Drs. Bob 
Haney and Bob Hale), using no sponge layer, very little boundary noise is 
generated and the fine grid solution is consistent with the coarse grid but 
richer in detail. We have extended that basic approach to a two-way-coupled duo 
grid North Atlantic Ocean adaptation of DieCAST~\citep{Dietrich04a}. 

TIMCOM open boundary conditions and nested grid coupling are based on upwind 
approximations. All open boundary conditions are based on boundary fluxes. Users 
can refer to~\citet{Tseng05a} for more detailed one or two way couplings. 
Point sources, such as from rivers smaller than one grid point, are treated by 
specifying salinity (zero), heat, momentum and volume inflows at specified 
source locations; this is done through vertical velocity and fluxes specified 
through a `semi-porous' rigid-lid but may be generalized for deeper river 
inflows. `e-p' (evaporation - precipitation) is similarly specified by vertical 
velocity relative to the modeled ocean surface (semi-porous rigid-lid). 

\section{Numerical Methods} 

The basic computation procedure during a model time step is as follows. 
The first step is to calculate the local cell averaged density, from potential 
temperature, salinity and pressure-depth. This uses an equation-of-state 
provided by Dan Wright. Then, starting with a trial `surface' pressure (or 
hydrostatically equivalent surface height) at the top of the surface layer of 
cells, the hydrostatic relation is integrated analytically to get face-averaged 
pressures at all other top and bottom cell faces. These are then used to get 
cell-averaged pressures within the all cells. Next, the latter are used to get 
face-averaged pressures at the four lateral faces of all cells; this gives the 
pressure force in the cell-averaged momentum equations. These steps are all 
calculated with fourth-order-accuracy. 

After partially updating horizontal velocity by adding the divergence of 
momentum flux plus diffusion fluxes to the pressure force effects, a top layer 
pressure adjustment is calculated such that, when added to the guessed pressure, 
the divergence of the vertically integrated horizontal velocity is zero (or, in 
the case of net specified e-p, a small specified value). This (semi-porous) 
rigid-lid pressure (hydrostatically equivalent to slow mode surface height) 
correction is governed by a Poisson equation for the top layer pressure 
adjustment; one must be careful that the boundary conditions do not violate 
global incompressibility -- there must be no net inflow in order for an 
incompressible velocity field to exist and for the Poisson equation for the top 
layer pressure adjustment to be well posed. After the face averaged velocities 
are adjusted according to the rigid-lid pressure adjustment, their CHANGES are 
interpolated back to give CHANGES of the cell averages, again using rigorously 
fourth-order-accurate interpolations. The algebraic details of this procedure 
are given in Appendices. 

Vertical velocity is derived simply by vertical integration of the divergence of 
the horizontal velocity, as per the incompressibility approximation to the mass 
conservation equation. 

The numerical details of this straightforward, rigorously fourth-order- accurate 
procedure are also documented by extensive comments among the lines of model 
code, and by Dietrich (1997). This formal accuracy gives an accurate, robust 
model. Only the barotropic mode pressure adjustment is second-order-accurate; 
yet, it is quite accurate due to the staggered grid arrangement between 
barotropic mode pressure and velocity. 

The above procedure includes the divergence of the barotropic mode of advection 
velocity, e.g., as implied by non-zero e-p in the framework of surface boundary 
conditions that do not require transient-damping nudging toward climatological 
data~\citep{Dietrich04b}. This new approach is implemented in recent TIMCOM 
model adaptations including the North Atlantic Ocean~\citep{Dietrich04a}. 

\section{Model Features} 

This documented TIMCOM version uses standard Boussinesq, hydrostatic and 
rigid-lid approximations. It includes: 

\begin{itemize}

\item Full spherical metrics with variable grid increments; 

\item A standard nonlinear equation of state relating in-situ density to 
potential temperature, salinity and pressure-depth; 

\item 3-D turbulent vertical diffusivity arrays; in the SOMS model from which 
TIMCOM evolved, these were determined by computation- and storage- intensive 
Mellor-Yamada level 2.5 turbulence closure~\citep{Dietrich87}; the present 
default is the~\citet{Pacanowski81} approach; the model user may apply 
any alternate preferred scheme. 

\item 3-D turbulent horizontal diffusivity arrays. Although no diffusive 
turbulence 
parameterization is deemed effective for horizontal Reynolds stress 
representation in general, these arrays are available for numerical and 
weakly physical based usages. 

\item Model program written in a modular structure, allowing fully customized by 
the user.

\end{itemize}

For efficiency in general circulation (slow mode) applications, TIMCOM uses a 
rigid-lid approximation. 

\section{Advanced Model features}
 
\subsection{Free Surface}
 
An explicit time-split free surface version exists (available under license) 
that integrates interlaced leap-frog shallow water equations using multiple `
fast mode' time steps each `slow mode' time step; this version addresses fast 
free surface modes (tides, wind surge events, etc) and is efficient, accurate 
and robust, including vertically averaged baroclinic pressure gradient and 
momentum flux forcings. The motivation for this is the first derivative coupling 
of different fields, similar to the motivation for spatial grids. 

Incidentally, the shallow water equations are a one-layer application of sigma 
coordinates. They are useful for the `fast modes' of ocean dynamics. Generally, 
these fast modes have only a very minor effect on the ocean `slow modes' that 
dominate ocean material trajectories and general circulation, other than perhaps 
turbulent dissipation from breaking wave effects in localized coastal regions. 
The extent to which such localized dissipation affects the general circulation, 
and thus climate, is not well understood, but short term local effects are quite 
small. 

Thus, a general circulation model may run independently, creating vertically 
averaged baroclinic pressure gradients and fluctuation correlations that may be 
used by a more specialized free surface `fast mode' model, while ignoring at 
least short term feedbacks from the barotropic mode (shallow water equations 
solution) to the baroclinic mode. Thus, models using the rigid-lid approximation 
can be quite useful, as proven by their many demonstrated successes in general 
circulation studies. 

\subsection{Nesting and Grid Coupling} 

TIMCOM supports Multiple gridded one- and two-way grid coupling. 
See~\citet{Dietrich04a} for basic methodology and results from application 
to North Atlantic Ocean. Model details available under licensing agreement.

\subsection{Hydrostatic and Non-hydrostatic}
 
TIMCOM Supports hydrostatic and non-hydrostatic approximations. For basic 
methodology, see~\citet{Dietrich02} and~\citet{Tseng03} for the numerical 
procedure. Model details available under licensing agreement.

\subsection{Other Advanced Features under Developing}

\begin{itemize}

\item Turbulence parameterization options 

\item Bio-geochemistry dynamics

\item Immersed boundary methods for topography( to be supported by 2010)

\item Lagrangian tracers releases capacity

\item Hindcast wind forcing capability

\item New data and code structures, NetCDF and customized I/O interface

\item Ensemble simulation capability

\end{itemize} 

\newpage
