source: trunk/DOC/BETA/Chapters/Chap_LBC.tex @ 781

Last change on this file since 781 was 781, checked in by smasson, 13 years ago

doc update, see ticket:1

  • Property svn:executable set to *
File size: 42.7 KB
1% ================================================================
2% Chapter Ñ Lateral Boundary Condition (LBC)
3% ================================================================
4\chapter{Lateral Boundary Condition (LBC) }
8%gm% add here introduction to this chapter
10% ================================================================
11% Boundary Condition at the Coast
12% ================================================================
13\section{Boundary Condition at the Coast (\np{shlat})}
19%The lateral ocean boundary condition continuous to coastlines are Neumann conditions for heat and salt (no flux across boundaries) and Dirichlet conditions for momentum (from free-slip to "strong" no-slip). They are handled automatically by the mask system (see \S\ref{DOM_msk}).
21%OPA works with land and topography grid points in the computational domain due to the presence of continents or islands,and to the use of full or partial step representation of bottom topography. The computation is performed over the whole domain, i.e. we did not try to restrict the computation to ocean point only areas. This choice has two motivations. First, working on ocean only grid point overload the code and harms the code readability. Second,, and more importantly, it drastically reduce the vector portion of the computation, leading to a dramatic increase of CPU time requirement on vector computers. The process of defining which areas are to be masked is described in \S\ref{DOM_msk}, while this section describes how the masking affects the computation of the various terms of the equations with respect to the boundary condition at solid walls.
23The discrete representation of a domain with complex boundaries (coastlines and bottom topography) leads to arrays that include large portions where a computation is not required as the model variables remain at zero. Nevertheless, vectorial supercomputers are far more efficient when computing over a whole array, and the readability of a code is greatly improved when boundary conditions are applied in an automatic way rather than by a specific computation before or after each do loop. An efficient way to work over the whole domain while specifying the boundary conditions is to use the multiplication by mask arrays in the computation. A mask array is a matrix
24which elements are $1 $in the ocean domain and $0$ elsewhere. A simple multiplication of a variable by its own mask ensures that it will remain zero over land areas. Since most of the boundary conditions consist of a zero flux across the solid boundaries, they can be simply settled by
25multiplying variables by the right mask arrays, i.e. the mask array of the grid point where the flux is evaluated. For example, the heat flux in the \textbf{i}-direction is evaluated at $u$-points. Evaluating this quantity as,
27\begin{equation} \label{Eq_lbc_aaaa}
28\frac{A^{lT} }{e_1 }\frac{\partial T}{\partial i}\equiv \frac{A_u^{lT} 
29}{e_{1u} }\delta _{i+1 / 2} \left[ T \right]\;\;mask_u
33\begin{figure}[!t] \label{Fig_LBC_uv}  \begin{center}
35\caption {Lateral boundary (thick line) at T-level. The velocity normal to the boundary is set to zero.}
36\end{center}   \end{figure}
39where mask$_{u}$ is the mask array at $u$-point, ensures that the heat flux is
40zero inside land and at the boundaries as mask$_{u}$ is zero at solid
41boundaries defined at $u$-points in this case (normal velocity $u$ remains zero at
42the coast) (Fig.~\ref{Fig_LBC_uv}).
44On momentum the situation is a bit more complex as two boundary conditions must be provided along the coast (one on the normal velocity and the other on the tangential velocity). The boundary of the ocean in C-grid is defined by the velocity-faces. For example, at a given $T$-level, the lateral boundary (coastline or intersection with the bottom topography) is made of segments
45joining $f$-points, and normal velocity points are located between two $f-$points (Fig.~\ref{Fig_LBC_uv}). The boundary condition on the normal velocity (no flux through solid boundaries) can thus be easily settled by the mask system. The boundary condition on the tangential velocity requires a more specific treatment. It influences the relative vorticity and momentum diffusive trends, and is only required to compute the vorticity at the coast. Four different type of the lateral boundary condition are available, controlled by the value of \np{shlat} namelist parameter (which is equal to the value of the mask$_{f}$ array along the coastline):
48\begin{figure}[!p] \label{Fig_LBC_shlat}  \begin{center}
50\caption {lateral boundary condition (a) free-slip ($shlat=0$) ; (b) no-slip ($shlat=2$) ; (c) "partial" free-slip ($0<shlat<2$) and (d) "strong" no-slip ($2<shlat$). Implied "ghost" velocity inside land area is display in grey. }
51\end{center}   \end{figure}
56\item[free-slip boundary condition (\np{shlat}=0): ]  the tangential velocity at the coastline is equal to the offshore velocity, $i.e.$ the normal derivative of the tangential velocity is zero at the coast, so the vorticity: mask$_{f}$ array is set to zero inside the land and just at the coast (Fig.~\ref{Fig_LBC_shlat}-a).
58\item[no-slip boundary condition (\np{shlat}=2): ] the tangential velocity vanishes at the coastline. Assuming that the tangential velocity decreases linearly from the closest ocean velocity grid point to the coastline, the normal derivative is evaluated as if the closest land velocity gridpoint were of the same magnitude as the closest ocean velocity gridpoint but in the opposite direction (Fig.~\ref{Fig_LBC_shlat}-b). Therefore, the vorticity along the coastlines is given by:
60\zeta \equiv 2 \left(\delta_{i+1/2} \left[e_{2v} v \right] - \delta_{j+1/2} \left[e_{1u} u \right] \right) / \left(e_{1f} e_{2f} \right) \ ,
62where $u$ and $v$ are masked fields. Setting the mask$_{f}$ array to $2$ along the coastline allows to provide a vorticity field computed with the no-slip boundary condition simply by multiplying it by the mask$_{f}$ :
63\begin{equation} \label{Eq_lbc_bbbb}
64\zeta \equiv \frac{1}{e_{1f} {\kern 1pt}e_{2f} }\left( {\delta _{i+1/2} 
65\left[ {e_{2v} \,v} \right]-\delta _{j+1/2} \left[ {e_{1u} \,u} \right]} 
69\item["partial" free-slip boundary condition (0$<$\np{shlat}$<$2): ] the tangential velocity at the coastline is smaller than the offshore velocity, $i.e.$ there is a lateral friction but not strong enough to vanish the tangential velocity at the coast (Fig.~\ref{Fig_LBC_shlat}-c). This can be settled by providing a value of mask$_{f}$ strictly inbetween $0$ and $2$.
71\item["strong" no-slip boundary condition (2$<$\np{shlat}): ] the viscous boundary layer is assumed to be smaller than half the grid size (Fig.~\ref{Fig_LBC_shlat}-d). The friction is thus larger than in the no-slip case.
75Note that when the bottom topography is entirely represented by the $s$-coordinates (pure $s$-coordinate), the lateral boundary condition on momentum tangential velocity is of much  little importance as it is only applied next to the coast where the minimum water depth can be quite shallow.
77The alternative numerical implementation of the no-slip boundary conditions for an arbitrary coast line of \citet{Shchepetkin1996} is also available through the activation of \key{noslip\_accurate}. It is based on a fourth order evaluation of the shear at the coast which, in turn, allows a true second order scheme in the interior of the domain ($i.e.$ the numerical boundary scheme simulates the truncation error of the numerical scheme used in the interior of the domain). \citet{Shchepetkin1996} found that such a technique considerably improves the quality of the numerical solution. In \NEMO, the improvement have not been found so spectacular in the half-degree global ocean (ORCA05), but significant reduction of numerically induced coast upwellings were found in eddy resolving simulation of the Alboran Sea \citep{OlivierPh2001}. Nevertheless, as no-slip boundary condition is not recommended in eddy permitting or resolving simulation \citep{Penduff2007}, the use of this option is not recommended.
79In practice, the no-slip accurate option changes the way the curl is evaluated at the coast (see \mdl{divcur} module), and requires to qualify the nature of coastline grid point (convex or concave corners, straight north-south or east-west coast) which is performed in \mdl{domask} module, \rou{dom\_msk\_nsa} routine.
81% ================================================================
82% Boundary Condition around the Model Domain
83% ================================================================
84\section{Model Domain Boundary Condition (\jp{jperio})}
87At the model domain boundaries several choices are offered: closed, cyclic east-west, south symmetric across the equator, a north-fold, and combination closed-north fold or cyclic-north-fold. The north-fold boundary condition is associated with the 3-pole ORCA mesh.
89% -------------------------------------------------------------------------------------------------------------
90%        Closed, cyclic, south symmetric (\jp{jperio} = 0, 1 or 2)
91% -------------------------------------------------------------------------------------------------------------
92\subsection{Closed, cyclic, south symmetric (\jp{jperio} = 0, 1 or 2)}
95The choice of closed, cyclic or symmetric model domain boundary condition is made by setting \jp{jperio} to 0, 1 or 2 in \mdl{par\_oce} file. Each time such a boundary condition is needed, it is set by a call to \mdl{lbclnk} routine. The computation of momentum and tracer trends proceed from $i=2$ to $i=jpi-1$ and from $j=2$ to $j=jpj-1$, $i.e.$in the inner model domain. To choose a lateral model boundary condition is to specify the first and last rows and columns of the model variables.
97- For closed boundary (\textit{jperio=0}), solid walls are imposed at all model boundaries:
98first and last rows and columns are set to zero.
100- For cyclic east-west boundary (\textit{jperio=1}), first and last rows are set to zero
101(closed) while first column is set to the value of the before last column
102and last column to the value of the second one (Fig.~\ref{Fig_LBC_jperio}-a). Whatever flows
103out of the eastern (western) end of the basin enters the western (eastern)
104end. Note that there is neither option for north-south cyclic nor doubly
105cyclic cases.
107- For symmetric boundary condition across the equator (\textit{jperio=2}), last rows, and
108first and last columns are set to zero (closed). The row of symmetry is
109chosen to be the $u$- and $T-$points equator line ($j=2$, i.e. at the southern end
110of the domain). For arrays defined at $u-$ or $T-$points, the first row is set to
111the value of the third row while for most of $v$- and $f$-points arrays (v, $\zeta$,
112$j\psi$, but scalar arrays such as eddy coefficients) the first row is set to
113minus the value of the second row (Fig.~\ref{Fig_LBC_jperio}-b). Note that this boundary
114condition is not yet available on massively parallel computer
115(\textbf{key{\_}mpp} defined).
118\begin{figure}[!t] \label{Fig_LBC_jperio}  \begin{center}
120\caption {setting of (a) east-west cyclic  (b) symmetric across the equator boundary conditions.}
121\end{center}   \end{figure}
124% -------------------------------------------------------------------------------------------------------------
125%        North fold (\textit{jperio = 3 }to $6)$
126% -------------------------------------------------------------------------------------------------------------
127\subsection{North-fold (\textit{jperio = 3 }to $6)$ }
130The north fold boundary condition has been introduced in order to handle the north boundary of a three-polar ORCA grid. Such a grid has two poles in the northern hemisphere. \colorbox{yellow}{to be completed...}
133\begin{figure}[!t] \label{Fig_North_Fold_T}  \begin{center}
135\caption {North fold boundary with a $T$-point pivot and cyclic east-west boundary condition ($jperio=4$), as used in ORCA 2, 1/4, and 1/12. Pink shaded area corresponds to the inner domain mask (see text). }
136\end{center}   \end{figure}
139% ====================================================================
140% Sub-Domain:Exchange with neighbouring processors
141% ====================================================================
142\section{Exchange with neighbouring processors (\mdl{lbclnk}, \mdl{lib\_mpp})}
145For massively parallel processing (mpp), a domain decomposition method is used. The basis of the method consists in splitting the large computation domain of a numerical experiment into several smaller domains and solving the set of equations by addressing independent local problems. Each processor has its own local memory and computes the model equation over a subdomain of the whole model domain. The subdomain boundary conditions are specified through communications between processors which are explicitly organized by specific statements (message passing method).
147A big advantage is that the method does not need many modifications of the initial FORTRAN code. For the modeller's point of view, each sub domain running on a processor is identical to the "mono-domain" code. In addition, the programmer manages the communications between subdomains, and the code presents more scalability when the number of processors is increased. The porting of OPA code on a iPSC860 was achieved during Guyon's PhD [Guyon et al. 1994, 1995] in collaboration with CETIIS and ONERA. The implementation in the operational context and the studies of performances on a T3D and T3E Cray computers have been made in collaboration with IDRIS and CNRS. The present implementation is largely inspired from Guyon's work  [Guyon 1995].
149   The parallelization strategy is defined by the physical characteristics of the ocean model. Second order finite difference schemes leads to local discrete operators that depend at the very most on one neighbouring point. The only non-local computations concerne the vertical physics (implicit diffusion, 1.5 turbulent closure scheme, ...) (delocalization over the whole water column), and the solving of the elliptic equation associated with the surface pressure gradient computation (delocalization over the whole horizontal domain). Therefore, a pencil strategy is used for the data sub-structuration: the 3D initial domain is laid out on local processor memories following a 2D horizontal topological splitting. Each sub-domain computes its own surface and bottom boundary conditions and has a side wall overlapping interface which stocks lateral boundary conditions for computations in the inner sub-domain. The overlapping area is reduced to one row. After a computation, a communication phase starts: each processor sends to its neighbouring processors the update values of the point corresponding to the overlapping area of its neighbouring sub-domain. The communication is done through message passing. Usually the parallel virtual language, PVM, is used as it is a standard language available on  nearly  all MPP cumputers. More specific languages (i.e. computer dependant languages) can be easily used to speed up the communication, such as SHEM on T3E computer. The data exchanges between processors are required at the very place where lateral domain boundary conditions are set in the mono-domain computation (\S III.10-c): the lbc\_lnk routine which manages such conditions is substituted by mpplnk.F or mpplnk2.F routine when running on MPP computer (\key{mpp\_mpi} defined). It has to be noticed that when using MPP version of the model, the east-west cyclic boundary condition is implicitly done, while the south-symmetric boundary condition option is not available.
152\begin{figure}[!t] \label{Fig_mpp}  \begin{center}
154\caption {Positioning of a sub-domain when massively parallel processing is used. }
155\end{center}   \end{figure}
158   In the standard version of the OPA model, the splitting is regular and arithmetic. the i-axis is divided by \jp{jpni} and the j-axis by \jp{jpnj} for a number of processors \jp{jpnij} most often equal to $jpni \times jpnj$ (model parameters set in \mdl{par\_oce}). Each processor is independent and without message passing or synchronous process, programs run alone and access just at its own local memory. For this reason, the main model dimensions are now the local dimensions of the subdomain (pencil) that are noted \jp{jpi}, \jp{jpj}, \jp{jpk}. These dimensions include the internal domain and the overlapping rows. The number of overlapping rows is usually set to one (\jp{jpreci}=1, in \mdl{par\_oce}). The whole domain dimensions are named \jp{jpiglo}, \jp{jpjglo} and \jp{jpk}. The relationship between the whole domain and a sub-domain is:
160      jpi & = & ( jpiglo-2*jpreci + (jpni-1) ) / jpni + 2*jpreci  \nonumber \\
161      jpj & = & ( jpjglo-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj  \label{Eq_lbc_jpi}
163where \jp{jpni}, \jp{jpnj} are the number of processors following the i- and j-axis.
165\colorbox{yellow}{Figure IV.3: example of a domain splitting with 9 processors and no east-west cyclic boundary conditions.}
167   One defines also variables nldi and nlei which correspond to the internal domain bounds, and the variables nimpp and njmpp which are the position of the (1,1) grid-point in the global domain. An element of $T_{l}$, a local array (subdomain) corresponds to an element of $T_{g}$, a global array (whole domain) by the relationship:
168\begin{equation} \label{Eq_lbc_nimpp}
169T_{g} (i+nimpp-1,j+njmpp-1,k) = T_{l} (i,j,k),
171with  $1 \leq i \leq jpi$, $1  \leq j \leq jpj $ , and  $1  \leq k \leq jpk$.
173   Processors are numbered from 0 to $jpnij-1$, the number is saved in the variable nproc. In the standard version, a processor has no more than four neighbouring processors named nono (for north), noea (east), noso (south) and nowe (west) and two variables, nbondi and nbondj, indicate the situation of the processor \colorbox{yellow}{(see Fig.IV.3)}:
175\item       nbondi = -1    an east neighbour, no west processor,
176\item       nbondi =  0 an east neighbour, a west neighbour,
177\item       nbondi =  1    no east processor, a west neighbour,
178\item       nbondi =  2    no splitting following the i-axis.
180   During the simulation, processors exchange data with their neighbours. If there is effectively a neighbour, the processor receives variables from this processor on its overlapping row, and sends the data issued from internal domain corresponding to the overlapping row of the other processor.
183\colorbox{yellow}{Figure IV.4: pencil splitting with the additional outer halos }
187   The OPA model computes equation terms with the help of mask arrays (0 onto land points and 1 onto sea points). It is easily readable and very efficient in the context of the vectorial architecture. But in the case of scalar processor, computations over the land regions becomes more expensive in term of CPU time. It is all the more so when we use a complex configuration with a realistic bathymetry like the global ocean where more than 50 \% of points are land points. For this reason, a pre-processing tool allows to search in the mpp domain decomposition strategy if a splitting can be found with a maximum number of only land points processors which could be eliminated: the mpp\_optimiz tools (available from the DRAKKAR web site). This optimisation is made with the knowledge of the specific bathymetry. The user chooses optimal parameters \jp{jpni}, \jp{jpnj} and \jp{jpnij} with $jpnij < jpni \times jpnj$, leading to the elimination of $jpni \times jpnj - jpnij$ land processors. When those parameters are specified in module \mdl{par\_oce}, the algorithm in the \rou{inimpp2} routine set each processor name and neighbour parameters (nbound, nono, noea,...) so that the land processors are not taken into account.
189\colorbox{yellow}{Note that the inimpp2 routine is general so that the original inimpp routine should be suppressed from the code.}
191When land processors are eliminated, the value corresponding to these locations in the model output files is zero. Note that this is a problem for a mesh output file written by such a model configuration, because model users often divide by the scale factors ($e1t$, $e2t$, etc) and do not expect the grid size to be zero, even on land. It may be best not to eliminate land processors when running the model especially to write the mesh files as outputs (when \np{nmsh} namelist parameter differs from 0).
194\begin{figure}[!ht] \label{Fig_mppini2}  \begin{center}
196\caption {Example of Atlantic domain defined for the CLIPPER projet. Initial grid is composed by 773 x 1236 horizontal points. (a) the domain is splitting onto 9 \time 20 subdomains (jpni=9, jpnj=20). 52 subdomains are land areas. (b) 52 subdomains are eliminated (white rectangles) and the resulting number of processors really used during the computation is jpnij=128.}
197\end{center}   \end{figure}
201% ================================================================
202% Open Boundary Conditions
203% ================================================================
204\section{Open Boundary Conditions (\key{obc})}
206%-----------------------------------------nam_obc  -------------------------------------------
207%-    nobc_dta    =    0     !  = 0 the obc data are equal to the initial state
208%-                           !  = 1 the obc data are read in 'obc   .dta' files
209%-    rdpein      =    1.    !  ???
210%-    rdpwin      =    1.    !  ???
211%-    rdpnin      =   30.    !  ???
212%-    rdpsin      =    1.    !  ???
213%-    rdpeob      = 1500.    !  time relaxation (days) for the east  open boundary
214%-    rdpwob      =   15.    !    "        "           for the west  open boundary
215%-    rdpnob      =  150.    !    "        "           for the north open boundary
216%-    rdpsob      =   15.    !    "        "           for the south open boundary
217%-    zbsic1      =  140.e+6 !  barotropic stream function on first  isolated coastline
218%-    zbsic2      =    1.e+6 !    "                   "    on second isolated coastline
219%-    zbsic3      =    0.    !    "                   "    on thrid  isolated coastline
220%-    ln_obc_clim = .true.   !  climatological obc data files (default T)
221%-    ln_vol_cst  = .true.   !  total volume conserved
224It is often necessary to implement a model configuration limited to an oceanic region or a basin, which communicates with the rest of the global ocean through  ``open boundaries''. As stated by \citet{Roed1986}, an open boundary is a computational border where the aim of the calculations is to allow the perturbations generated inside the computational domain to leave it without deterioration of the inner model solution. However, an open boundary also has to let information from the outer oceans enter the model and should support inflow and outflow conditions.
226The open boundary package OBC is the first open boundary option developed in NEMO (originally in OPA8.2). It allows the user to
228\item tell the model that a boundary is ``open'' and not closed by a wall, for example by modifying the calculation of the divergence of velocity there, and other necessary modifications;
229\item impose values of tracers and velocities at that boundary (values which may be taken from a climatology): this is the``fixed OBC'' option.
230\item calculate boundary values by a sophisticated algorithm combining radiation and relaxation (``radiative OBC'' option)
233The package resides in the OBC directory. It is described here in four parts: the boundary geometry (parameters to be set in \mdl{obc\_par}), the forcing data at the boundaries (module \mdl{obcdta}),  the radiation algorithm involving the namelist and module \mdl{obcrad}, and a brief presentation of boundary update and restart files.
236\subsection{Boundary geometry}
239First one has to realize that open boundaries may not necessarily be located at the extremities of the computational domain. They may exist in the middle of the domain, for example at Gibraltar Straits if one wants to avoid including the Mediterranean in an Atlantic domain. This flexibility has been found necessary for the CLIPPER project \citep{Treguier2001}. Because of the complexity of the geometry of ocean basins, it may even be necessary to have more than one ``west'' open boundary, more than one ``north'', etc. This is not possible with the OBC option: only one open boundary of each kind, west, east, south and north is allowed; those names refer to the grid geometry (not to the direction of the geographical ``west'', ``east'', etc).
241The open boundary geometry is set by a series of parameters in the module \mdl{obc\_par}. For an eastern open boundary, parameters are \jp{lp\_obc\_east} (true if an east open boundary exists), \jp{jpieob} the $i$-index along which the eastern open boundary (eob) is located, \jp{jpjed} the $j$-index at which it starts, and \jp{jpjef} the $j$-index where it ends (note $d$ is for ``d\'{e}but'' and $f$ for ``fin'' in French). Similar parameters exist for the west, south and north cases (Table~\ref{Tab_obc_param}).
246\begin{table}[htbp]  \label{Tab_obc_param}
250Boundary and  & Constant index  & Starting index (d\'{e}but) & Ending index (fin) \\
251Logical flag  &                 &                            &                     \\
253West          & \jp{jpiwob} $>= 2$         &  \jp{jpjwd}$>= 2$          &  \jp{jpjwf}<= \jp{jpjglo}-1 \\
254lp\_obc\_west & $i$-index of a $u$ point   & $j$ of a $T$ point   &$j$ of a $T$ point \\
256East            & \jp{jpieob}$<=$\jp{jpiglo}-2&\jp{jpjed} $>= 2$         & \jp{jpjef}$<=$ \jp{jpjglo}-1 \\
257 lp\_obc\_east  & $i$-index of a $u$ point    & $j$ of a $T$ point & $j$ of a $T$ point \\
259South           & \jp{jpjsob} $>= 2$         & \jp{jpisd} $>= 2$          & \jp{jpisf}$<=$\jp{jpiglo}-1 \\
260lp\_obc\_south  & $j$-index of a $v$ point   & $i$ of a $T$ point   & $i$ of a $T$ point \\
262North           & \jp{jpjnob} $<=$ \jp{jpjglo}-2& \jp{jpind} $>= 2$        & \jp{jpinf}$<=$\jp{jpiglo}-1 \\
263lp\_obc\_north  & $j$-index of a $v$ point      & $i$  of a $T$ point & $i$ of a $T$ point \\
267\caption{Names of different indices concerning the open boundaries. In the case of a completely open ocean domain with four ocean boundaries, the parameters take exactly the values indicated.}
270The open boundaries must be along coordinate lines. On the C-grid, the boundary itself is along a line of normal velocity points: $v$ points for a zonal open boundary (the south or north one), and $u$ points for a meridional open boundary (the west or east one). Another constraint is that there still must be a row of masked points all around the domain, as if the domain were a closed basin (unless periodic conditions are used together with open boundary conditions). Therefore, an open boundary cannot be located at the first/last index, namely, 1 or \jp{jpiglo}, \jp{jpjglo}. Also, the open boundary algorithm involves calculating the normal velocity points situated just on the boundary, as well as the tangential velocity and temperature, salinity just outside the boundary. This means that for a west/south boundary, normal velocities and temperature are calculated at the same index \jp{jpiwob} and \jp{jpjsob}, respectively. For an east/north boundary, the normal velocity is calculated at index \jp{jpieob} and \jp{jpjnob}, but the ``outside'' temperature is at index \jp{jpieob}+1 and \jp{jpjnob}+1. This means that \jp{jpieob}, \jp{jpjnob} cannot be bigger than \jp{jpiglo}-2, \jp{jpjglo}-2.
273The starting and ending indices are to be thought as $T$ point indices: in many cases they indicate the first land $T$-point, at the extremity of an open boundary (the coast line follows the $f$ grid points, see Fig.~\ref{Fig_obc_north} for the example for a northern open boundary. All indices are relative to the global domain. In the free surface case it is possible to have ``ocean corners'', that is, an open boundary starting and ending in the ocean.
276\begin{figure}[!t] \label{Fig_obc_north}  \begin{center}
278\caption {Localization of the North open boundary points.}
283Although not absolutely compulsory, it is highly recommended that the bathymetry in the vicinity of an open boundary follows the following rule: in the direction perpendicular to the open line, the water depth should be constant for 4 grid points. This is in order to ensure that the radiation condition, which involves mdoel variables next to the boundary, is calculated in a consistent way. On Fig.\ref{Fig_obc_north} we indicate by a $=$ symbol, the points which should have the same depth. It means that in the 4 points near the boundary, the bathymetry is cylindrical. The line behind the open T-line must be 0 in the bathymetry file (as shown on Fig.\ref{Fig_obc_north} for example).
286\subsection{Boundary data}
289It is necessary to provide information at the boundaries. The simple case happens when this information does not change in time and is equal to the initial conditions (namelist variable \np{nobc\_dta}=0). This is the case of the standard configuration EEL5 with open boundaries. When (\np{nobc\_dta}=1), open boundary information is read from netcdf files. For convenience the input files are supposed to be similar to the ''history'' NEMO output files, for dimension names and variable names.
290Open boundary arrays must be dimensioned according to the parameters of table~\ref{Tab_obc_param}: for example, at the western boundary arrays have a dimension of \jp{jpwf}-\jp{jpwd}+1 on the horizontal and \jp{jpk} on the vertical.
292When ocean observations are used to generate the boundary data (a hydrographic section for example, as in \citet{Treguier2001}) it happens often that only the velocity normal to the boundary is known, which is the reason why the initial OBC code assumes that only $T$, $S$, and the normal velocity ($u$ or $v$) needs to be specified. As more and more global model solutions and ocean analysis products are available, it is possible to provide information about all the variables (including the tangential velocity) so that the specification of four variables at each boundaries will soon become standard. Regarding the sea surface height, one must distinguish the filtered free surface case from the time-splitting or explicit treatment of the free surface.
293 In the first case, it is assumed that the user does not wish to represent high frequency motions such as tides. The boundary condition is thus one of zero normal gradient of sea surface height at the open boundaries, following \citet{Marchesiello2001}.
294No information other than the total velocity needs to be provided at the open boundaries in that case. In the other two cases (time stplitting or explicit free surface), the user must provides barotropic information (sea surface height and barotropic velocities) and the use of the Flather algorithm for barotropic variables is recommanded. However, this algorithm has not yet been fully tested and bugs remain in NEMO v2.3. Users should read the code carefully before using it. Finally, in the case of the rigid lid approximation the barotropic streamfunction must be provided, as documented in \citet{Treguier2001}). This option is no longer used but remains in NEMO V2.3.
296One frequently encountered case is when an open boundary domain is constructed from a global or larger scale NEMO configuration. Assuming the domain corresponds to indices $ib:ie$, $jb:je$ of the global domain, the bathymetry and forcing of the small domain can be created by using the following command on the global files: ncks -F $-d\;x,ib,ie$ $-d\;y,jb,je$.
297The open boundary files can be constructed using ncks commands, following table~\ref{Tab_obc_ind}.
301\begin{table}[htbp]  \label{Tab_obc_ind}
305OBC  & Variable   & file name      & Index  & Start  & end  \\
306West &  T,S       &   obcwest\ &  $ib$+1     &   $jb$+1 &  $je-1$  \\
307     &    U       &   obcwest\  &  $ib$+1     &   $jb$+1 &  $je-1$  \\ 
308     &    V       &   obcwest\  &  $ib$+1     &   $jb$+1 &  $je-1$  \\       
310East &  T,S       &   obceast\ &  $ie$-1     &   $jb$+1 &  $je-1$  \\
311     &    U       &   obceast\  &  $ie$-2     &   $jb$+1 &  $je-1$  \\ 
312     &    V       &   obceast\  &  $ie$-1     &   $jb$+1 &  $je-1$  \\       
314South &  T,S      &   obcsouth\ &  $jb$+1     &  $ib$+1 &  $ie-1$  \\
315      &    U      &   obcsouth\  &  $jb$+1     &  $ib$+1 &  $ie-1$  \\ 
316      &    V      &   obcsouth\  &  $jb$+1     &  $ib$+1 &  $ie-1$  \\   
318North &  T,S      &   obcnorth\ &  $je$-1     &  $ib$+1 &  $ie-1$  \\
319      &    U      &   obcnorth\  &  $je$-1     &  $ib$+1 &  $ie-1$  \\ 
320      &    V      &   obcnorth\  &  $je$-2     &  $ib$+1 &  $ie-1$  \\ 
324\caption{Indications for creating open boundary files from a global configuration, appropriate for the subdomain of indices $ib:ie$, $jb:je$. ``Index'' designates the $i$ or $j$ index along which the $u$ of $v$ boundary point is situated in the global configuration, starting and ending with the $j$ or $i$ indices indicated. For example, to generate file obcnorth\, use the command ncks $-F$ $-d\;y,je-2$  $-d\;x,ib+1,ie-1$ } 
327It is assumed that the open boundary files contain the variables for the period of the model integration. If the boundary files contain one time frame, the boundary data is held fixed in time. If the files contain 12 values, it is assumed that the input is a climatology for a repeated annual cycle (corresponding to the case \np{ln\_obc\_clim} = .T.). The case of an arbitrary number of time frames is not yet implemented correctly; the user is supposed to write his own code in the module \mdl{obc\_dta} to deal with this situation.
329\subsection{Radiation algorithm}
332The art of open boundary management consists in applying a constraint strong enough so that the inner domain "feels" the rest of the ocean, but not too strong
333in order to allow perturbations to leave the domain with minimum false reflexions of energy. The cosntraint to the specified boundary data is set by time scales
334specified separately for each boundary and for ``inflow'' and `outflow'' conditions as defined below. These time scales are set (in days) by namelist parameters such as \np{rdpein}, \np{rdpeob} for the eastern open boundary, for example. When both time scales are zero for a given boundary, for example
335\jp{lp\_obc\_west}=.T., \np{rdpwob}=0, \np{rdpwin}=0, this means that the boundary in question (western boundary here) is a ``fixed '' boundary where the solution is set exactly by the boundary data. This is not recommanded, excepted in compination with increased viscosity in a ``sponge'' layer next to the boundary in order to avois spurious reflexions. 
338The radiation\/relaxation algorithm is applied when either relxation time (for ``inflow'' or `outflow'') is nonzero. It has been developed and tested in the SPEM model and its successor ROMS \citep{Barnier1996, Marchesiello2001}, a $s$-coordinate model on an Arakawa C-grid. Although the algorithm has been numerically successful in the CLIPPER Atlantic models, the physics do not work as expected \citep{Treguier2001}. Users are invited to consider open boundary conditions (OBC hereafter) with some skepticism \citep{Durran2001, Blayo2005}.
340The first part of the algorithm consists in calculating a phase
341velocity to determine whether perturbations tend to propagate toward,
342or away from, the boundary.
343Let us consider a model variable $\phi$.
344The phase velocity ($C_{\phi x}$,$C_{\phi y}$) for the variable
345$\phi$, in the directions normal and tangential to
346the boundary is
347\begin{equation} \label{Eq_obc_cphi}
348C_{\phi x} = \frac{ -\phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{x} 
349\;\;\;\;\; \;\;\; 
350C_{\phi y} = \frac{ -\phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{y}.
352Following \citet{Treguier2001} and \citet{Marchesiello2001} retain only the normal
353projection of the total velocity, $C_{\phi x}$, setting
354$C_{\phi y} =0$ (but unlike the original Orlanski radiation algorithm we retain $\phi_{y}$ in the expression
355for $C_{\phi x}$). 
357The discrete form of (\ref{Eq_obc_cphi}), described by \citet{Barnier1998},
358takes into account the two rows of grid points situated inside the domain
359next to the boundary, and the three previous time steps ($n$, $n-1$,
360and $n-2$).
361The same equation can then be discretized at the boundary at
362time steps $n$ and $n+1$ in order to extrapolate the new boundary
363value $\phi^{n+1}$.
365In the present open boundary algorithm, the new boundary values are
366updated differently according to the sign of $C_{\phi x}$.
367Let us take an East boundary as an example. The solution for variable $\phi$ at the boundary is given from a generalized wave equation
368with phase velocity $C_{\phi}$, with the addition of  a relaxation term:
370\phi_{t} &  =  & -C_{\phi x} \phi_{x} + \frac{1}{\tau_{o}} (\phi_{c}-\phi)
371                        \;\;\; \;\;\; \;\;\; (C_{\phi x} > 0), \label{Eq_obc_rado} \\
372\phi_{t} &  =  & \frac{1}{\tau_{i}} (\phi_{c}-\phi)
373\;\;\; \;\;\; \;\;\;\;\;\; (C_{\phi x} < 0), \label{Eq_obc_radi}
375where $\phi_{c}$ is the estimate of $\phi$ at the boundary, provided as boundary data.
376Note that in (\ref{Eq_obc_rado}), $C_{\phi x}$
377is bounded by the ratio $\delta x/\delta t$ for stability reasons.
378When $C_{\phi x}$ is eastward (outward propagation),
379the radiation condition (\ref{Eq_obc_rado}) is used.
380When  $C_{\phi x}$ is westward (inward propagation), (\ref{Eq_obc_radi}) is used with
381a strong relaxation to climatology (usually $\tau_{i}=\np{rdpein}=$1~day).
382The time derivative in (\ref{Eq_obc_radi}) is calculated with a Euler time-stepping
383scheme. It results from this choice that setting
384$\tau_{i}$ smaller than, or equal to the time step is in fact equivalent to a fixed boundary condition;
385a time scale of one day is usually a good compromise which guarantees that the inflow conditions remain close to
386climatology while ensuring numerical stability.
388In  the case of a west boundary located in the Eastern Atlantic, \citet{Penduff2000} have been able to implement the radiation algorithm
389without any boundary data, using persistence from the previous time step instead. This solution did not work in other cases \citep{Treguier2001} so that the use of boundary data is recommended. 
390Even in the outflow condition (\ref{Eq_obc_rado}), we have found desireable to maintain a weak relaxation to climatology.
391The time step is usually chosen so as to be larger than typical turbulent scales (of order 1000~days).
393The radiation condition is applied to the different model variables: temperature, salinity, tangential and normal velocities.
394For normal and tangential velocities $u$ and $v$ radiation is applied with phase velocities calculated from $u$ and $v$ respectively.
395For the radiation of tracers, we use the phase velocity calculated from the tangential velocity, to avoid calculating too many independent
396radiation velocities and because tangential velocities and tracers have the same position along the boundary on a C grid. 
398\subsection{Domain decomposition (\key{mpp\_mpi})}
400When \key{mpp\_mpi} is active in the code, the computational domain is divided into rectangles that are attributed each to a different processor. The open boundary code is ``mpp-compatible'' up to a certain point. The radiation algorithm will not work if there is a mpp subdomain boundary parallel to the open boundary at the index of the boundary, or the grid point after (outside), or three grid points before (inside). On the other hand, there is no problem if a mpp subdomain boundary cuts the open boundary perpendicularly. Those geometry constraints must be checked by the user (there is no safeguard in the code). 
401The general principle for the open boundary mpp code is that loops over the open boundaries are performed on local indices (nie0, nie1, nje0, nje1 for the eastern boundary for instance) that are initialized in module \mdl{obc\_ini}. Those indices have relevant values on the processors that contain a segment of an open boundary. For processors that do not include an open boundary segment, the indices are such that the calculations within the loops are not performed.
403Arrays of climatological data that are read in files are seen by all processors and have the same dimensions for all (for instance, for the eastern boundary, uedta(jpjglo,jpk,2)). On the other hand, the arrays for the calculation of radiation are local to each processor (uebnd(jpj,jpk,3,3) for instance).  This allowed to spare memory in the CLIPPER model, where the eastern boundary crossed 8 processors so that \jp{jpj} was must smaller than (\jp{jpjef}-\jp{jpjed}+1).
405\subsection{Volume conservation}
408It is necessary to control the volume inside a domain with open boundaries. With fixed boundaries, it is enough to ensure that the total inflow/outflow has reasonable values (either zero or a value compatible with an observed volume balance). When using radiative boundary conditions it is necessary to have a volume constraint because each open boundary works independently from the others. The methodology used to control this volume is identical to the one coded in the ROMS model \citep{Marchesiello2001}.
411%---------------------------------------- EXTRAS
412\colorbox{yellow}{Explain obc\_vol{\ldots}}
414\colorbox{yellow}{OBC algorithm for update, OBC restart, list of routines where obc key appears{\ldots}}
416\colorbox{yellow}{OBC rigid lid? {\ldots}}
421% ====================================================================
422% Flow Relaxation Scheme
423% ====================================================================
424\section{Flow Relaxation Scheme (???)}
427%gm% to be updated by Met Office
Note: See TracBrowser for help on using the repository browser.