% Code Structure

\chapter{Code Organization}
\label{chap:structure}

TIMCOM is written in Fortran 90 programming language and is executable on 
Windows or Unix-Like platform. The whole program contain: preprocessors,
main model program, and postprocessors. The preprocessors read 
bathymetric and climatological data required for the simulation, generate
grid coordinate according to the input data, initialize variables and set
boundary conditions. The main model program controls calculation, performs 
diagnostics and saves graphics data. Coupling block is adapted for the 
commuication between different coupled domains. Finally, third party 
postprocessors are used to view the output data.

To run a customized case, user-defined options can be applied to specify input 
data, computation parameters, processing, and output, by modifying the 
preconfigured option file. The current version of TIMCOM supports Intel and pgi 
compilers.

\begin{figure}[htbp]
\centerline{\scalebox{.35}{\includegraphics{FIGS/Code_Structure.eps}}}
\caption{Coding structure. ETOPO5, Levitus, Hellerman indicate data type of 
bathymetry, salinity/temperature and wind stress.}
\label{fig:Code_Structure}
\end{figure}

\begin{figure}[htbp]
\centerline{\scalebox{.35}{\includegraphics{FIGS/Process.eps}}}
\caption{Running procedure. Figure on the left illustrates the single domain 
case, and for the right hand one, domain TAI is coupled with domain NPB.}
\label{fig:Process}
\end{figure}

\section{Program Directories}

How the program is organized and distributed is explained in this section. 
Figure~\ref{fig:Dir_Structure_Single} and Figure~\ref{fig:Dir_Structure_Couple}
illustrate the file layout of TIMCOM of both single and coupled domain,
respectively.

\begin{figure}[htbp]
\centerline{\scalebox{.5}{\includegraphics{FIGS/Dir_Structure_Single.eps}}}
\caption{Directory tree for single domain model.}
\label{fig:Dir_Structure_Single}

\end{figure}

\begin{description}

\item[TIMCOMGLO]: Root directory, including all source, data, configuration and 
output files.

\item[BACKUP]: Backup of configuration files, {\bf ocn\_para.f90.bak}, 
{\bf precfg.txt.bak}, and {\bf run.sh.bak}. 
When users are about to run a new case, these files are 
copied and modified to satisfy specific requirements.

\item[CODE]: Where the source code files are stored and compiled.

\item[DATA]: The climatological data: bathymetry, salainity, 
temperature and wind stress reside within. 
Different types of data have their own subdirectory which are explained
in Chapter~\ref{sec:data_input}.

\item[EXAMPLES]: Upon running a new case, all configuration, temporary and 
output files are located inside a specific directory, which is named by the 
users, under the folder {\bf EXAMPLES}.

\item[EXAMPLES/CASE1]: Case name, e.g. home directory for a specified case,
which is customizable by the user.

\item[TEMP]: Stores temporary files during the simulation. 

\item[OUTPUT]: Where the output files generated by the model are stored,
within the following subdirecties:

\begin{itemize}
\item \textbf{OUTPUT/XY}: horizontal plots.
\item \textbf{OUTPUT/XZ}: vertical east-west cross-sections. 
\item \textbf{OUTPUT/YZ}: meridional cross-sections. 
\end{itemize}

In addition, by setting the optional flag \textbf{LMOVI} in the 
configuration file {\bf precfg.txt} from $0$ (off) to $1$ (on), 
computer-generated movies of TIMCOM output data is written to the \textbf{MOVIE} 
subdirectory: 

\begin{itemize}
\item \textbf{OUTPUT/MOVIE}: horizontal field animation (pressure contours and 
tracer streaklines) 
\end{itemize}

In order to minimize storage, output fields are scaled and written in binary 
\mbox{INTEGER*2} 
form and   locations are derived from the small basic metrics data sets 
rather than combined with the data which would be hugely wasteful of storage. 
Scale factors are used to recover original field values in graphics 
postprocessors.

\end{description}

\begin{figure}[htbp]
\centerline{\scalebox{.5}{\includegraphics{FIGS/Dir_Structure_Couple.eps}}}
\caption{Directory tree for coupled-domain model.}
\label{fig:Dir_Structure_Couple}
\end{figure}

In the case of coupled domain, the directory layout is very similar to the one
of single domain, except that the root directory is changed from {\bf TIMCOMGLO}
to {\bf TIMCOMNPNT} for distinguishing the case difference, and each domain has 
its own subdirectories under the case directory. In this case, 
they are {\bf D01} and {\bf D02}, which are customizable by the
user. More detail about this will be given later.

\section{Data Input}
\label{sec:data_input}

The following bathymetric and climatological data is required for TIMCOM to
carry out the simulation:

\begin{itemize}

\item Etopo5, the `5-Minute Gridded Global Relief Data Collection', maintained
by the National Geophysical Data Center of National Oceanic \& Atmospheric 
Administration. Etopo5 provides data of bathymetric map, or the land and 
sea-floor elevations, on a 5-minute latitude/longitude grid.

\item World Ocean Atlas, or the Levitus climatology, is the data product of the 
Ocean Climate Laboratory of the National Oceanographic Data Center, contains 
annual, seasonal, and monthly long term means for temperature and salinity at 
multiple depths.

\item Hellerman and Rosenstein Global Wind Stress Climatology, this data set 
contains monthly averaged wind stress over the globalocean for the years 1870 
through 1976.

\end{itemize}

% Appendix reference required!!

Once the user has prepared the data mentioned above, it is necessary to specifiy
the directories in where the data locates, and this is done in a shell script
file called {\bf run.sh}, which is already provided along with the source files. 
The content of {\bf run.sh} is self-explained, but an 
example of {\bf run.sh} with comment will be given in Appendix.

%Preprocessors read the data mentioned above and generate gridded input file
%for TIMCOM which include 
%source files \textbf{metgen.f90}, \textbf{inmet.f90}, \textbf{indata.f90} and 
%\textbf{\mbox{prep.f90}}. More detail of these files will be given in Section ~
%\ref{sec:preprocessor}.

%They are responsible for:

%\begin{description}

%\item[metgen.f90] generates grid coordinate for the computation of TIMCOM.

%\item[inmet.f90] read the Etopo5 data then interpolate it onto the grid 
%coordinate created
%by \textbf{metgen.f90}.

%\item[indata.f90] read the Levitus and Hellerman data then interpolate them onto 
%the grid coordinate created by \textbf{metgen.f90}. 
%Very similar to \textbf{inmet.f90}.

%\item[prep.f90] prepares and initialzes any other variables based on user's 
%condigurations, such as Coriolis parameter array, land/sea mask array, boundary 
%conditions, etc.

%\end{description}

\section{Program Structure}

Written in a modular approach, the program of TIMCOM is organized
in a well manner that each block of the sorce code performs their
own tasks without interference from, or reference to, the other blocks,
except essential data exchange. The modular approach makes the modification
and improvement of the source code easy as well. This block-structure is shown 
in Figure~\ref{fig:TIMCOM_Block}.

\begin{figure}[htbp]
\centerline{\scalebox{.4}{\includegraphics{FIGS/TIMCOM_Block.eps}}}
\caption{The structure of TIMCOM is simple. Data and variables are
encapsulated and protected in the corresponding block, while each block
is independent from others.}
\label{fig:TIMCOM_Block}
\end{figure}

\subsection{Variable Modules}
\label{sec:variable_modules}

Variables are distributed, declared and defined across source files,
encapsulated in different modules, 
depending on what purposes the variables are assigned.
This section explains the significant of each source file containing
variable modules. The suffix `{\bf xx}' in the source file name represents
in which domain the file is responsible for. For single domain simulation,
these suffixes are not required, but for the case of coupled-domain,
each domain may have its own corresponding variable, thus suffixes are
necessary.\newline

\noindent \textbf{ocn\_para.f90} \newline

In this file, the resolution parameters of TIMCOM are declared and
defined, such as 
$I0$, $J0$ and $K0$, as well as some derived constants: $I1=I0-1$,
$J1=J0-1$, $K1=K0-1$, as mentioned in~\ref{sec:computational_grid}. 
There're also dimension parameters for SEVP elliptic solver:
$NB0$ and $NB1$. \newline 

\noindent \textbf{grid\_var\_xx.f90}  \newline

Metric data describing the coordinate structure, such as control volume 
dimensions:$DX$ and $DY$; mask arrays: $IV$ and $IW$; cell-surface averaged 
velocities:$U$, $V$ and $W$; vertical turbulent viscosity and
diffusivity arrays: $EV$ and $HV$; density:$\rho$. \newline

\noindent \textbf{control\_var\_xx.f90} \newline

Parameters that control how the model runs reside in this file, including time 
step, feature on/off flag such as whether inlude the effect of wind forcing
and surface heating of air-sea exchange or not, 
and physical parameters such as bottom drag coefficient and 
direct turbulance viscocity. Most of these options and parameters are set 
through editing the preconfiguration file {\bf precfg.txt}. \newline

\noindent \textbf{common\_var.f90} \newline

This is the coupling block, only available when multi-domain model is used to 
simulate area with complex bathymetry. Each domain has similar code structure 
except independent variables inside their corrosponding domain. 
But some variables are still of global 
significance, such as case name, climatological data file, and the variables 
used to pass information between each domain. How model domains are 
coupled is illustrated in Figure~\ref{fig:coupled_domain_example}. \newline

\begin{figure}[htbp]
\centerline{\scalebox{.6}{\includegraphics{FIGS/NORPAC2_COUPLING.eps}}}
\caption{Domain NorPac and Tai are coupled in this case. 
As shown in the figure, a portion of 
grids of both domains are overlap and may even have different grid size.
Coupling block of the source code is adapted in this region to exchange 
information between different domains such that the numerical methods
are applied smoothly.}
\label{fig:coupled_domain_example}
\end{figure}

\noindent \textbf{init\_var\_xx.f90}  \newline

Most variables are accessed and modified only by the model itself, user knows 
nothing about them when the model is running. But before the model runs, some 
parameters must be specified, so that the model knows where to find the 
climatology data, how communicate with the OS, how long it should run and what 
data format user would like to use. These parameters are declared in 
{\bf init\_var\_xx.f90} and classified 
into corresponding namelist. Some of them are newly defined, while others 
refer to other modules defined in {\bf grid\_var\_xx.f90} or 
{\bf control\_xx\_var.f90}.
The user-configurable parameters are set in a plain text file {\bf precfg.txt}
which is provided along with the source files. More detail on it would be given 
in the Appendix.

% Appendix reference required!!

\subsection{Preprocessor}
\label{sec:preprocessor}

Input data to \textbf{prep.f90} is generated by programs \textbf{metgen.f90},
\textbf{inmets.f90} and \textbf{indata.f90}. \newline

\noindent \textbf{metgen.f90} \newline

Grid coordinate is generated through this subroutine on which the governing
equations are applied. Figure ~\ref{fig:Compute_Grid} shows how grid 
parameters are laid out in $X$ and $Y$ directions, 
which are important for grid generation:

\begin{itemize}

\item DX and DY: The lateral grid increments across grid lines.

\item DXV and DYV: The lateral grid increments across grid centers.

\item Y and YV: Latitude cell centers and latitude grid lines.

\item YDEG and YVDEG: Location of latitude cell centers and grid lines in
degree.

\item X0DEG and Y0DEG: Describe the westmost and southmost location of the
modeled region in degree.

\item X1DEG and Y1DEG: Describe the eastmost and northmost location of the
modeled region in degree. 

\end{itemize}

{\bf metgen.f90} generates grids in $X$ and $Y$ directions by itself. 
As for the $Z$ direction, it calls different subroutines in {\bf metgen\_sub.f90}
according to user's setting in the pre-configuration file {\bf precfg.txt}. 
Detail about grid generation in the $Z$ direction will be 
given in Chapter~\ref{sec:z_grid_coordinate}.  \newline \newpage

\begin{figure}[htbp]
\centerline{\scalebox{.35}{\includegraphics{FIGS/Compute_Grid.eps}}}
\caption{In horizontal directions, ghost zones are used for convenient 
specification of inflows and outflows at open boundary points. The conservation 
equations are applied only to interior cells. Colored-grids represent the
ghost zones.}
\label{fig:Compute_Grid}
\end{figure}

%The following important grid parameters are used in this file:

%\begin{itemize}

%\item Resolution parameter: $I0$, $J0$, $K0$. They are the cell number
%in the $X$, $Y$ and $Z$ direction, respectively, as mentioned in 
%Chapter~\ref{sec:computational_grid}.

%\item In the $Z$ direction: $TLZ$, the total depth of the modeled domain. 

%\item In $X$-$Y$ direction: $DY$ and $DX$, $DXV$ and $DYV$the lateral grid 
%increments across grid line and cell-centers, respectively.
%There are $D0DEG$ and $Y0DEG$ as well, 
%which describes the southmost and westmost location of the model region.

%\end{itemize}

\noindent \textbf{inmets.f90}   \newline

This subroutine read and interpolate Etopo5 bathymetric data onto the model grid
generated by {\bf metgen.f90}.
In general, resolution and data position of bathymetry will not match that of 
model grid as shown in Figure~\ref{fig:data_interpolation},
then interpolation is required. \newline

\begin{figure}[htbp]
\centerline{\scalebox{.5}{\includegraphics{FIGS/Itopo5-grid.eps}}}
\caption{Generally the data layout of bathymetric data (dash line) and the one of 
the model region (solid line) are different. Interpolation is done
to distribute bathymetric data properly onto the model grid.}
\label{fig:data_interpolation}
\end{figure}

\noindent \textbf{indata.f90}    \newline

The same as {\bf inmets.f90}, but this one deals with the
Levitus and Hellerman data, which record salinity/temperature and 
wind stress data. Interpolation is done as well. \newline

\noindent \textbf{prep.f90}    \newline

Gridded salinity, temperature, wind stress and bathymetry data are not read and 
applied by the model directly. For computation convenience, some assumption and 
simplification are applied, such as boundary condition.
The code \textbf{prep.f90} creates the following files for input to TIMCOM: a 
two-
dimensional topography array $KB$($I$,$J$) giving the number of layers that are 
water- 
occupied at a given horizontal grid point ($I$,$J$), and other masking and 
geometry arrays; files for initial temperature and salinity (not needed when 
restarting from previous model results); surface wind forcing and restoring data 
files; open lateral boundary files; and matrix data for the EVP surface pressure 
solver. 
One may write alternate program(s) to replace \textbf{prep.f90}, or create the 
depth, intial and boundary data that \textbf{prep.f90} reads and processes for 
the input of TIMCOM. 

\subsection{Main Model Program}

The main model of TIMCOM is made up with the following source files: \newline

\noindent \textbf{timcom.f90}  \newline

This is the main program. It controls the program flow by calling
other subroutines for each phase of TIMCOM such as preprocessing and
main computation.      \newline

\noindent \textbf{initfsxx.f90}  \newline

This subroutine reads input data, grid coordinate,
salinity, temperature, wind stress, boundary conditions and setup
configurations prepared by the preprocessors to feed into the TIMCOM model. \newline                

\noindent \textbf{fsxx.f90}   \newline

The main computation of TIMCOM is done with this file, including the
discretization of governing equations and equations of state, 
conservation laws, open boundary treatments, computation of internal 
longitudinal and latitudinal fluxes,
hydrostatic and non-hydrostatic pressure gradient, incorporation of the
climatological effects such as wind stress and heating, and the EVP solver
for pressure correction.
In other words, it performs the required operations for each corresponding
time step. \newline

\noindent \textbf{climatxx.f90}   \newline

The module defines variables and subroutines corresponding to external driving 
such as heating and wind stress prepared by the preprocessors.  \newline

The iterative computation procedures of the main model are illustrated in
Figure~\ref{fig:main_flowchart}

\begin{figure}[htbp]
\centerline{\scalebox{.5}{\includegraphics{FIGS/Main_Flowchart.eps}}}
\caption{The computation procedures of TIMCOM.}
\label{fig:main_flowchart}
\end{figure}

\subsection{Postprocessors}

TIMCOM containes only the following two postprocessors for
creating visual data according to the result given by the
main model program. \newline  

\noindent \textbf{graphicxx.f90} \newline

Contains subroutines creating graphic data for viewing with 3rd party
programs.  \newline

\noindent \textbf{windows.f90}  \newline   

Defines subroutines for extractomg parts of model output for specific 
sub-domain.      

\section{Z Coordinate System}
\label{sec:z_grid_coordinate}

Parameters related to the description of the $Z$ coordinates,
including $Z$ coordinate position and total number of layers, 
are shown in Figure~\ref{fig:Z_Direction}.   

Although Cartesian coordinate system with uniformly distrubuted
grids is the simplest system for formulation and computation, 
it may not be suitable for everycase. For the case of ocean circulation 
simulation, the bathymetric map is complex, thus finer grid my be required
for higher accuracy near the bottom of the model region, 
which is approximated with stair-step grid. 

In order to satisfy different simulation requirement,
one important feature of TIMCOM is that different
$Z$-coordinate can be chosen: 
linear exponential coordinate, stretched coordinate and 
user-customized coordinate. This is done by setting the option flag {\bf ZOP}
in the pre-configuration file {\bf precfg.txt} to 1, 2 or 3, 
corresponfing to the three aforementioned coordinate systems. \newline

\begin{figure}[htbp]
\centerline{\scalebox{.45}{\includegraphics{FIGS/Z_Direction.eps}}}
\caption{In $z$ direction, there are $K0-1$ layers in total. $K0$ indicate the 
number of interfaces of the layers. $TLZ$ stands for the max depth and  
$Z(k)(k=1,2,\ldots 2 \times K0-1)$ specifies $Z$-level coordinate. 
For computational convenience, the depth value of layer centers are stored
as well, 
so we get $2 \times K0-1$ nodes describing $Z$ coordinate which also equals
to $K0+K1$. No ghost zones are used in $Z$ direction.}
\label{fig:Z_Direction}
\end{figure}

\noindent \textbf{a) Linear Exponential $Z$ Coordinate}\newline

By setting {\bf ZOP} in {\bf precfg.txt} to 1, linear exponential coordinate
will be utilized. The following equations are used to transform $Z$ into $zz$,
as shown in Figure~\ref{fig:Exp_Z_Direction}. 
\newline

\noindent $Z(zz)=A+B\times Exp(C\times zz)+D\times zz,$ \qquad $0 < zz < TLZ$ \newline

Where $zz$ is the linear (unstretched) space coordinate
for expanding coordinate, $0 < d < 1$ \newline

\noindent At $zz=0$, \qquad $Z=A+B \times Exp(0)=0$, \qquad so $A+B=0$ \newline
At $zz=TLZ$, \qquad $Z=A+B \times Exp(C \times TLZ)+D \times TLZ=TLZ$  \newline
\newline
so $A \times (1-Exp(C \times TLZ))=TLZ-D \times TLZ$ \newline 
then $A=TLZ \times (1-D)/(1-Exp(C \times TLZ)), \qquad B=-A$ \newline  

If parameters $C$ and $D$ are provided by the user, which is also done
through editing {\bf precgf.txt}, then $A$ and $B$ be determined.

\begin{figure}[htbp]
\centerline{\scalebox{.3}{\includegraphics{FIGS/Exp_Z_Direction.eps}}}
\caption{Linear exponential $Z$ coordinate. $ODZ$ is the inversed layer
thickness, the distance between two layers; $ODZW$ is the inversed distance
between layer pressure levels.}
\label{fig:Exp_Z_Direction}
\end{figure}

According to Figure ~\ref{fig:Exp_Z_Direction}, the following parameters
could be derived by the program automatically:  \newline

\noindent$DZ=0.01 \times TLZ /(2.0 \times K1)$ \newline
$ODZ(K)=1/(Z(2 \times K+1)-Z(2 \times K-1))$ \qquad ,$K=1,2, \ldots K1$ \newline
$ODZW(K)=1/(Z(2 \times K)-Z(2 \times K-2))$ \qquad, $K=2,3, \ldots K1$
\newline

\noindent\textbf{b) Stretched $Z$ Coordinate}

\begin{figure}[htbp]
\centerline{\scalebox{.3}{\includegraphics{FIGS/Stretch_Z_Direction.eps}}}
\caption{Stretched $Z$ coordinate, where $ZTOP$ is the thickness of the top
layer, whose value is defined in {\bf precfg.txt} by the user.}
\label{fig:Stretch_Z_Direction}
\end{figure}

Stretched coordinate is simpler than the linear exponential one. By specifying
the thickness of the top layer: $ZTOP$, any further grid will be stretched
(or contracted) along the $Z$ direction as the depth getting deeper, by a 
given ratio: $GX$, which is defined in the subroutine {\bf SUB\_Z\_STRETCH of}
{\bf metgen\_sub.f90}. That is, the transformed coordinate $zz$ could be
derived as: \newline
 
\noindent $zz(2)=ZTOP$,\newline 
$zz(K)=zz(K-1)+ZTOP\times GX^{(N-2)}, \qquad K=3,4,\ldots K1$ \newline

\noindent\textbf{c) User Defined $Z$ Coordinate}  \newline

User defined coordinate is fully supported, but this required edit 
the subroutine {\bf SUB\_Z\_USRDEF} in source file {\bf metgen\_sub.f90}
directly, thus recompiling of programe is necessary. The program stores
the $Z$ coordinate position with the array: $Z$, whose dimension
is defined as $K0+K1$. So if users would like to
define their own coordinate, variable declaration and assignment
must be done carfully. 
Since subroutine {\bf SUB\_Z\_USRDEF} is called in the preprocessor: 
{\bf metgen.f90}, passing and managment of parameters should be treated with
caution as well. In general, the subroutines {\bf SUB\_Z\_LINEAR}
and {\bf SUB\_Z\_STRETCH} in {\bf metgen\_sub.f90} responsible for liner
exponential and stretched coordinates could be taken as reference,
   

%\begin{figure}[ht]
%\centerline{\scalebox{.35}{\includegraphics{FIGS/User_Z_Direction.eps}}}
%\caption{In this case, user is able to  define the depth of each layer: Z(K),  
%where $K=1,2, \ldots 2 \times K0-1$.}
%\label{fig:User_Z_Direction}
%\end{figure}


\newpage