source: branches/2013/dev_r3853_CNRS9_ConfSetting/DOC/TexFiles/Chapters/Chap_LBC.tex @ 3989

Last change on this file since 3989 was 3989, checked in by clevy, 8 years ago

Configuration setting/Step3 and doc, see ticket:#1074

File size: 58.0 KB
Line 
1% ================================================================
2% Chapter � Lateral Boundary Condition (LBC)
3% ================================================================
4\chapter{Lateral Boundary Condition (LBC) }
5\label{LBC}
6\minitoc
7
8\newpage
9$\ $\newline    % force a new ligne
10
11
12%gm% add here introduction to this chapter
13
14% ================================================================
15% Boundary Condition at the Coast
16% ================================================================
17\section{Boundary Condition at the Coast (\np{rn\_shlat})}
18\label{LBC_coast}
19%--------------------------------------------nam_lbc-------------------------------------------------------
20\namdisplay{namlbc} 
21%--------------------------------------------------------------------------------------------------------------
22
23%The lateral ocean boundary conditions contiguous to coastlines are Neumann conditions for heat and salt (no flux across boundaries) and Dirichlet conditions for momentum (ranging from free-slip to "strong" no-slip). They are handled automatically by the mask system (see \S\ref{DOM_msk}).
24
25%OPA allows land and topography grid points in the computational domain due to the presence of continents or islands, and includes the use of a full or partial step representation of bottom topography. The computation is performed over the whole domain, i.e. we do not try to restrict the computation to ocean-only points. This choice has two motivations. Firstly, working on ocean only grid points overloads the code and harms the code readability. Secondly, and more importantly, it drastically reduces the vector portion of the computation, leading to a dramatic increase of CPU time requirement on vector computers.  The current section describes how the masking affects the computation of the various terms of the equations with respect to the boundary condition at solid walls. The process of defining which areas are to be masked is described in \S\ref{DOM_msk}.
26
27Options are defined through the \ngn{namlbc} namelist variables.
28The discrete representation of a domain with complex boundaries (coastlines and
29bottom topography) leads to arrays that include large portions where a computation
30is not required as the model variables remain at zero. Nevertheless, vectorial
31supercomputers are far more efficient when computing over a whole array, and the
32readability of a code is greatly improved when boundary conditions are applied in
33an automatic way rather than by a specific computation before or after each
34computational loop. An efficient way to work over the whole domain while specifying
35the boundary conditions, is to use multiplication by mask arrays in the computation.
36A mask array is a matrix whose elements are $1$ in the ocean domain and $0$ 
37elsewhere. A simple multiplication of a variable by its own mask ensures that it will
38remain zero over land areas. Since most of the boundary conditions consist of a
39zero flux across the solid boundaries, they can be simply applied by multiplying
40variables by the correct mask arrays, $i.e.$ the mask array of the grid point where
41the flux is evaluated. For example, the heat flux in the \textbf{i}-direction is evaluated
42at $u$-points. Evaluating this quantity as,
43
44\begin{equation} \label{Eq_lbc_aaaa}
45\frac{A^{lT} }{e_1 }\frac{\partial T}{\partial i}\equiv \frac{A_u^{lT} 
46}{e_{1u} } \; \delta _{i+1 / 2} \left[ T \right]\;\;mask_u
47\end{equation}
48(where mask$_{u}$ is the mask array at a $u$-point) ensures that the heat flux is
49zero inside land and at the boundaries, since mask$_{u}$ is zero at solid boundaries
50which in this case are defined at $u$-points (normal velocity $u$ remains zero at
51the coast) (Fig.~\ref{Fig_LBC_uv}).
52
53%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
54\begin{figure}[!t]     \begin{center}
55\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_LBC_uv.pdf}
56\caption{  \label{Fig_LBC_uv}
57Lateral boundary (thick line) at T-level. The velocity normal to the boundary is set to zero.}
58\end{center}   \end{figure}
59%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
60
61For momentum the situation is a bit more complex as two boundary conditions
62must be provided along the coast (one each for the normal and tangential velocities).
63The boundary of the ocean in the C-grid is defined by the velocity-faces.
64For example, at a given $T$-level, the lateral boundary (a coastline or an intersection
65with the bottom topography) is made of segments joining $f$-points, and normal
66velocity points are located between two $f-$points (Fig.~\ref{Fig_LBC_uv}).
67The boundary condition on the normal velocity (no flux through solid boundaries)
68can thus be easily implemented using the mask system. The boundary condition
69on the tangential velocity requires a more specific treatment. This boundary
70condition influences the relative vorticity and momentum diffusive trends, and is
71required in order to compute the vorticity at the coast. Four different types of
72lateral boundary condition are available, controlled by the value of the \np{rn\_shlat} 
73namelist parameter. (The value of the mask$_{f}$ array along the coastline is set
74equal to this parameter.) These are:
75
76%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
77\begin{figure}[!p] \begin{center}
78\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_LBC_shlat.pdf}
79\caption{     \label{Fig_LBC_shlat} 
80lateral boundary condition (a) free-slip ($rn\_shlat=0$) ; (b) no-slip ($rn\_shlat=2$)
81; (c) "partial" free-slip ($0<rn\_shlat<2$) and (d) "strong" no-slip ($2<rn\_shlat$).
82Implied "ghost" velocity inside land area is display in grey. }
83\end{center}    \end{figure}
84%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
85
86\begin{description}
87
88\item[free-slip boundary condition (\np{rn\_shlat}=0): ]  the tangential velocity at the
89coastline is equal to the offshore velocity, $i.e.$ the normal derivative of the
90tangential velocity is zero at the coast, so the vorticity: mask$_{f}$ array is set
91to zero inside the land and just at the coast (Fig.~\ref{Fig_LBC_shlat}-a).
92
93\item[no-slip boundary condition (\np{rn\_shlat}=2): ] the tangential velocity vanishes
94at the coastline. Assuming that the tangential velocity decreases linearly from
95the closest ocean velocity grid point to the coastline, the normal derivative is
96evaluated as if the velocities at the closest land velocity gridpoint and the closest
97ocean velocity gridpoint were of the same magnitude but in the opposite direction
98(Fig.~\ref{Fig_LBC_shlat}-b). Therefore, the vorticity along the coastlines is given by:
99
100\begin{equation*}
101\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) \ ,
102\end{equation*}
103where $u$ and $v$ are masked fields. Setting the mask$_{f}$ array to $2$ along
104the coastline provides a vorticity field computed with the no-slip boundary condition,
105simply by multiplying it by the mask$_{f}$ :
106\begin{equation} \label{Eq_lbc_bbbb}
107\zeta \equiv \frac{1}{e_{1f} {\kern 1pt}e_{2f} }\left( {\delta _{i+1/2} 
108\left[ {e_{2v} \,v} \right]-\delta _{j+1/2} \left[ {e_{1u} \,u} \right]} 
109\right)\;\mbox{mask}_f
110\end{equation}
111
112\item["partial" free-slip boundary condition (0$<$\np{rn\_shlat}$<$2): ] the tangential
113velocity at the coastline is smaller than the offshore velocity, $i.e.$ there is a lateral
114friction but not strong enough to make the tangential velocity at the coast vanish
115(Fig.~\ref{Fig_LBC_shlat}-c). This can be selected by providing a value of mask$_{f}$ 
116strictly inbetween $0$ and $2$.
117
118\item["strong" no-slip boundary condition (2$<$\np{rn\_shlat}): ] the viscous boundary
119layer is assumed to be smaller than half the grid size (Fig.~\ref{Fig_LBC_shlat}-d).
120The friction is thus larger than in the no-slip case.
121
122\end{description}
123
124Note that when the bottom topography is entirely represented by the $s$-coor-dinates
125(pure $s$-coordinate), the lateral boundary condition on tangential velocity is of much
126less importance as it is only applied next to the coast where the minimum water depth
127can be quite shallow.
128
129The alternative numerical implementation of the no-slip boundary conditions for an
130arbitrary coast line of \citet{Shchepetkin1996} is also available through the
131\key{noslip\_accurate} CPP key. It is based on a fourth order evaluation of the shear at the
132coast which, in turn, allows a true second order scheme in the interior of the domain
133($i.e.$ the numerical boundary scheme simulates the truncation error of the numerical
134scheme used in the interior of the domain). \citet{Shchepetkin1996} found that such a
135technique considerably improves the quality of the numerical solution. In \NEMO, such
136spectacular improvements have not been found in the half-degree global ocean
137(ORCA05), but significant reductions of numerically induced coastal upwellings were
138found in an eddy resolving simulation of the Alboran Sea \citep{Olivier_PhD01}.
139Nevertheless, since a no-slip boundary condition is not recommended in an eddy
140permitting or resolving simulation \citep{Penduff_al_OS07}, the use of this option is also
141not recommended.
142
143In practice, the no-slip accurate option changes the way the curl is evaluated at the
144coast (see \mdl{divcur} module), and requires the nature of each coastline grid point
145(convex or concave corners, straight north-south or east-west coast) to be specified. 
146This is performed in routine \rou{dom\_msk\_nsa} in the \mdl{domask} module.
147
148% ================================================================
149% Boundary Condition around the Model Domain
150% ================================================================
151\section{Model Domain Boundary Condition (\np{jperio})}
152\label{LBC_jperio}
153
154At the model domain boundaries several choices are offered: closed, cyclic east-west,
155south symmetric across the equator, a north-fold, and combination closed-north fold
156or cyclic-north-fold. The north-fold boundary condition is associated with the 3-pole ORCA mesh.
157
158% -------------------------------------------------------------------------------------------------------------
159%        Closed, cyclic, south symmetric (\np{jperio} = 0, 1 or 2)
160% -------------------------------------------------------------------------------------------------------------
161\subsection{Closed, cyclic, south symmetric (\np{jperio} = 0, 1 or 2)}
162\label{LBC_jperio012}
163
164The choice of closed, cyclic or symmetric model domain boundary condition is made
165by setting \np{jperio} to 0, 1 or 2 in namelist \ngn{namcfg}. Each time such a boundary
166condition is needed, it is set by a call to routine \mdl{lbclnk}. The computation of
167momentum and tracer trends proceeds from $i=2$ to $i=jpi-1$ and from $j=2$ to
168$j=jpj-1$, $i.e.$ in the model interior. To choose a lateral model boundary condition
169is to specify the first and last rows and columns of the model variables.
170
171\begin{description}
172
173\item[For closed boundary (\textit{jperio=0})], solid walls are imposed at all model
174boundaries: first and last rows and columns are set to zero.
175
176\item[For cyclic east-west boundary (\textit{jperio=1})], first and last rows are set
177to zero (closed) whilst the first column is set to the value of the last-but-one column
178and the last column to the value of the second one (Fig.~\ref{Fig_LBC_jperio}-a).
179Whatever flows out of the eastern (western) end of the basin enters the western
180(eastern) end. Note that there is no option for north-south cyclic or for doubly
181cyclic cases.
182
183\item[For symmetric boundary condition across the equator (\textit{jperio=2})],
184last rows, and first and last columns are set to zero (closed). The row of symmetry
185is chosen to be the $u$- and $T-$points equator line ($j=2$, i.e. at the southern
186end of the domain). For arrays defined at $u-$ or $T-$points, the first row is set
187to the value of the third row while for most of $v$- and $f$-point arrays ($v$, $\zeta$,
188$j\psi$, but \gmcomment{not sure why this is "but"} scalar arrays such as eddy coefficients)
189the first row is set to minus the value of the second row (Fig.~\ref{Fig_LBC_jperio}-b).
190Note that this boundary condition is not yet available for the case of a massively
191parallel computer (\textbf{key{\_}mpp} defined).
192
193\end{description}
194
195%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
196\begin{figure}[!t]     \begin{center}
197\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_LBC_jperio.pdf}
198\caption{    \label{Fig_LBC_jperio}
199setting of (a) east-west cyclic  (b) symmetric across the equator boundary conditions.}
200\end{center}   \end{figure}
201%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
202
203% -------------------------------------------------------------------------------------------------------------
204%        North fold (\textit{jperio = 3 }to $6)$
205% -------------------------------------------------------------------------------------------------------------
206\subsection{North-fold (\textit{jperio = 3 }to $6)$ }
207\label{LBC_north_fold}
208
209The north fold boundary condition has been introduced in order to handle the north
210boundary of a three-polar ORCA grid. Such a grid has two poles in the northern hemisphere.
211\colorbox{yellow}{to be completed...}
212
213%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
214\begin{figure}[!t]    \begin{center}
215\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_North_Fold_T.pdf}
216\caption{    \label{Fig_North_Fold_T} 
217North fold boundary with a $T$-point pivot and cyclic east-west boundary condition
218($jperio=4$), as used in ORCA 2, 1/4, and 1/12. Pink shaded area corresponds
219to the inner domain mask (see text). }
220\end{center}   \end{figure}
221%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
222
223% ====================================================================
224% Exchange with neighbouring processors
225% ====================================================================
226\section  [Exchange with neighbouring processors (\textit{lbclnk}, \textit{lib\_mpp})]
227      {Exchange with neighbouring processors (\mdl{lbclnk}, \mdl{lib\_mpp})}
228\label{LBC_mpp}
229
230For massively parallel processing (mpp), a domain decomposition method is used.
231The basic idea of the method is to split the large computation domain of a numerical
232experiment into several smaller domains and solve the set of equations by addressing
233independent local problems. Each processor has its own local memory and computes
234the model equation over a subdomain of the whole model domain. The subdomain
235boundary conditions are specified through communications between processors
236which are organized by explicit statements (message passing method).
237
238A big advantage is that the method does not need many modifications of the initial
239FORTRAN code. From the modeller's point of view, each sub domain running on
240a processor is identical to the "mono-domain" code. In addition, the programmer
241manages the communications between subdomains, and the code is faster when
242the number of processors is increased. The porting of OPA code on an iPSC860
243was achieved during Guyon's PhD [Guyon et al. 1994, 1995] in collaboration with
244CETIIS and ONERA. The implementation in the operational context and the studies
245of performance on a T3D and T3E Cray computers have been made in collaboration
246with IDRIS and CNRS. The present implementation is largely inspired by Guyon's
247work  [Guyon 1995].
248
249The parallelization strategy is defined by the physical characteristics of the
250ocean model. Second order finite difference schemes lead to local discrete
251operators that depend at the very most on one neighbouring point. The only
252non-local computations concern the vertical physics (implicit diffusion, 1.5
253turbulent closure scheme, ...) (delocalization over the whole water column),
254and the solving of the elliptic equation associated with the surface pressure
255gradient computation (delocalization over the whole horizontal domain).
256Therefore, a pencil strategy is used for the data sub-structuration
257\gmcomment{no idea what this means!}
258: the 3D initial domain is laid out on local processor
259memories following a 2D horizontal topological splitting. Each sub-domain
260computes its own surface and bottom boundary conditions and has a side
261wall overlapping interface which defines the lateral boundary conditions for
262computations in the inner sub-domain. The overlapping area consists of the
263two rows at each edge of the sub-domain. After a computation, a communication
264phase starts: each processor sends to its neighbouring processors the update
265values of the points corresponding to the interior overlapping area to its
266neighbouring sub-domain (i.e. the innermost of the two overlapping rows).
267The communication is done through message passing. Usually the parallel virtual
268language, PVM, is used as it is a standard language available on  nearly  all
269MPP computers. More specific languages (i.e. computer dependant languages)
270can be easily used to speed up the communication, such as SHEM on a T3E
271computer. The data exchanges between processors are required at the very
272place where lateral domain boundary conditions are set in the mono-domain
273computation (\S III.10-c): the lbc\_lnk routine which manages such conditions
274is substituted by mpplnk.F or mpplnk2.F routine when running on an MPP
275computer (\key{mpp\_mpi} defined). It has to be pointed out that when using
276the MPP version of the model, the east-west cyclic boundary condition is done
277implicitly, whilst the south-symmetric boundary condition option is not available.
278
279%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
280\begin{figure}[!t]    \begin{center}
281\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_mpp.pdf}
282\caption{   \label{Fig_mpp} 
283Positioning of a sub-domain when massively parallel processing is used. }
284\end{center}   \end{figure}
285%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
286
287In the standard version of the OPA model, the splitting is regular and arithmetic.
288 the i-axis is divided by \jp{jpni} and the j-axis by \jp{jpnj} for a number of processors
289 \jp{jpnij} most often equal to $jpni \times jpnj$ (model parameters set in
290 \mdl{par\_oce}). Each processor is independent and without message passing
291 or synchronous process
292 \gmcomment{how does a synchronous process relate to this?},
293 programs run alone and access just its own local memory. For this reason, the
294 main model dimensions are now the local dimensions of the subdomain (pencil)
295 that are named \jp{jpi}, \jp{jpj}, \jp{jpk}. These dimensions include the internal
296 domain and the overlapping rows. The number of rows to exchange (known as
297 the halo) is usually set to one (\jp{jpreci}=1, in \mdl{par\_oce}). The whole domain
298 dimensions are named \np{jpiglo}, \np{jpjglo} and \jp{jpk}. The relationship between
299 the whole domain and a sub-domain is:
300\begin{eqnarray} 
301      jpi & = & ( jpiglo-2*jpreci + (jpni-1) ) / jpni + 2*jpreci  \nonumber \\
302      jpj & = & ( jpjglo-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj  \label{Eq_lbc_jpi}
303\end{eqnarray}
304where \jp{jpni}, \jp{jpnj} are the number of processors following the i- and j-axis.
305
306\colorbox{yellow}{Figure IV.3: example of a domain splitting with 9 processors and
307no east-west cyclic boundary conditions.}
308
309One also defines variables nldi and nlei which correspond to the internal
310domain bounds, and the variables nimpp and njmpp which are the position
311of the (1,1) grid-point in the global domain. An element of $T_{l}$, a local array
312(subdomain) corresponds to an element of $T_{g}$, a global array
313(whole domain) by the relationship:
314\begin{equation} \label{Eq_lbc_nimpp}
315T_{g} (i+nimpp-1,j+njmpp-1,k) = T_{l} (i,j,k),
316\end{equation}
317with  $1 \leq i \leq jpi$, $1  \leq j \leq jpj $ , and  $1  \leq k \leq jpk$.
318
319Processors are numbered from 0 to $jpnij-1$, the number is saved in the variable
320nproc. In the standard version, a processor has no more than four neighbouring
321processors named nono (for north), noea (east), noso (south) and nowe (west)
322and two variables, nbondi and nbondj, indicate the relative position of the processor
323\colorbox{yellow}{(see Fig.IV.3)}:
324\begin{itemize}
325\item       nbondi = -1    an east neighbour, no west processor,
326\item       nbondi =  0 an east neighbour, a west neighbour,
327\item       nbondi =  1    no east processor, a west neighbour,
328\item       nbondi =  2    no splitting following the i-axis.
329\end{itemize}
330During the simulation, processors exchange data with their neighbours.
331If there is effectively a neighbour, the processor receives variables from this
332processor on its overlapping row, and sends the data issued from internal
333domain corresponding to the overlapping row of the other processor.
334       
335\colorbox{yellow}{Figure IV.4: pencil splitting with the additional outer halos }
336
337
338The \NEMO model computes equation terms with the help of mask arrays (0 on land
339points and 1 on sea points). It is easily readable and very efficient in the context of
340a computer with vectorial architecture. However, in the case of a scalar processor,
341computations over the land regions become more expensive in terms of CPU time.
342It is worse when we use a complex configuration with a realistic bathymetry like the
343global ocean where more than 50 \% of points are land points. For this reason, a
344pre-processing tool can be used to choose the mpp domain decomposition with a
345maximum number of only land points processors, which can then be eliminated.
346(For example, the mpp\_optimiz tools, available from the DRAKKAR web site.)
347This optimisation is dependent on the specific bathymetry employed. The user
348then chooses optimal parameters \jp{jpni}, \jp{jpnj} and \jp{jpnij} with
349$jpnij < jpni \times jpnj$, leading to the elimination of $jpni \times jpnj - jpnij$ 
350land processors. When those parameters are specified in module \mdl{par\_oce},
351the algorithm in the \rou{inimpp2} routine sets each processor's parameters (nbound,
352nono, noea,...) so that the land-only processors are not taken into account.
353
354\colorbox{yellow}{Note that the inimpp2 routine is general so that the original inimpp
355routine should be suppressed from the code.}
356
357When land processors are eliminated, the value corresponding to these locations in
358the model output files is zero. Note that this is a problem for a mesh output file written
359by such a model configuration, because model users often divide by the scale factors
360($e1t$, $e2t$, etc) and do not expect the grid size to be zero, even on land. It may be
361best not to eliminate land processors when running the model especially to write the
362mesh files as outputs (when \np{nn\_msh} namelist parameter differs from 0).
363%%
364\gmcomment{Steven : dont understand this, no land processor means no output file
365covering this part of globe; its only when files are stitched together into one that you
366can leave a hole}
367%%
368
369%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
370\begin{figure}[!ht]     \begin{center}
371\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_mppini2.pdf}
372\caption {    \label{Fig_mppini2}
373Example of Atlantic domain defined for the CLIPPER projet. Initial grid is
374composed of 773 x 1236 horizontal points.
375(a) the domain is split onto 9 \time 20 subdomains (jpni=9, jpnj=20).
37652 subdomains are land areas.
377(b) 52 subdomains are eliminated (white rectangles) and the resulting number
378of processors really used during the computation is jpnij=128.}
379\end{center}   \end{figure}
380%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
381
382
383% ================================================================
384% Open Boundary Conditions
385% ================================================================
386\section{Open Boundary Conditions (\key{obc}) (OBC)}
387\label{LBC_obc}
388%-----------------------------------------nam_obc  -------------------------------------------
389%-    nobc_dta    =    0     !  = 0 the obc data are equal to the initial state
390%-                           !  = 1 the obc data are read in 'obc   .dta' files
391%-    rn_dpein      =    1.    !  ???
392%-    rn_dpwin      =    1.    !  ???
393%-    rn_dpnin      =   30.    !  ???
394%-    rn_dpsin      =    1.    !  ???
395%-    rn_dpeob      = 1500.    !  time relaxation (days) for the east  open boundary
396%-    rn_dpwob      =   15.    !    "        "           for the west  open boundary
397%-    rn_dpnob      =  150.    !    "        "           for the north open boundary
398%-    rn_dpsob      =   15.    !    "        "           for the south open boundary
399%-    ln_obc_clim = .true.   !  climatological obc data files (default T)
400%-    ln_vol_cst  = .true.   !  total volume conserved
401\namdisplay{namobc} 
402
403It is often necessary to implement a model configuration limited to an oceanic
404region or a basin, which communicates with the rest of the global ocean through
405''open boundaries''. As stated by \citet{Roed1986}, an open boundary is a
406computational border where the aim of the calculations is to allow the perturbations
407generated inside the computational domain to leave it without deterioration of the
408inner model solution. However, an open boundary also has to let information from
409the outer ocean enter the model and should support inflow and outflow conditions.
410
411The open boundary package OBC is the first open boundary option developed in
412NEMO (originally in OPA8.2). It allows the user to
413\begin{itemize}
414\item tell the model that a boundary is ''open'' and not closed by a wall, for example
415by modifying the calculation of the divergence of velocity there;
416\item impose values of tracers and velocities at that boundary (values which may
417be taken from a climatology): this is the``fixed OBC'' option.
418\item calculate boundary values by a sophisticated algorithm combining radiation
419and relaxation (``radiative OBC'' option)
420\end{itemize}
421
422Options are defined through the \ngn{namobc} namelist variables.
423The package resides in the OBC directory. It is described here in four parts: the
424boundary geometry (parameters to be set in \mdl{obc\_par}), the forcing data at
425the boundaries (module \mdl{obcdta}),  the radiation algorithm involving the
426namelist and module \mdl{obcrad}, and a brief presentation of boundary update
427and restart files.
428
429%----------------------------------------------
430\subsection{Boundary geometry}
431\label{OBC_geom}
432%
433First one has to realize that open boundaries may not necessarily be located
434at the extremities of the computational domain. They may exist in the middle
435of the domain, for example at Gibraltar Straits if one wants to avoid including
436the Mediterranean in an Atlantic domain. This flexibility has been found necessary
437for the CLIPPER project \citep{Treguier_al_JGR01}. Because of the complexity of the
438geometry of ocean basins, it may even be necessary to have more than one
439''west'' open boundary, more than one ''north'', etc. This is not possible with
440the OBC option: only one open boundary of each kind, west, east, south and
441north is allowed; these names refer to the grid geometry (not to the direction
442of the geographical ''west'', ''east'', etc).
443
444The open boundary geometry is set by a series of parameters in the module
445\mdl{obc\_par}. For an eastern open boundary, parameters are \jp{lp\_obc\_east} 
446(true if an east open boundary exists), \jp{jpieob} the $i$-index along which
447the eastern open boundary (eob) is located, \jp{jpjed} the $j$-index at which
448it starts, and \jp{jpjef} the $j$-index where it ends (note $d$ is for ''d\'{e}but''
449and $f$ for ''fin'' in French). Similar parameters exist for the west, south and
450north cases (Table~\ref{Tab_obc_param}).
451
452
453%--------------------------------------------------TABLE--------------------------------------------------
454\begin{table}[htbp]     \begin{center}    \begin{tabular}{|l|c|c|c|}
455\hline
456Boundary and  & Constant index  & Starting index (d\'{e}but) & Ending index (fin) \\
457Logical flag  &                 &                            &                     \\
458\hline
459West          & \jp{jpiwob} $>= 2$         &  \jp{jpjwd}$>= 2$          &  \jp{jpjwf}<= \np{jpjglo}-1 \\
460lp\_obc\_west & $i$-index of a $u$ point   & $j$ of a $T$ point   &$j$ of a $T$ point \\
461\hline
462East            & \jp{jpieob}$<=$\np{jpiglo}-2&\jp{jpjed} $>= 2$         & \jp{jpjef}$<=$ \np{jpjglo}-1 \\
463 lp\_obc\_east  & $i$-index of a $u$ point    & $j$ of a $T$ point & $j$ of a $T$ point \\
464\hline
465South           & \jp{jpjsob} $>= 2$         & \jp{jpisd} $>= 2$          & \jp{jpisf}$<=$\np{jpiglo}-1 \\
466lp\_obc\_south  & $j$-index of a $v$ point   & $i$ of a $T$ point   & $i$ of a $T$ point \\
467\hline
468North           & \jp{jpjnob} $<=$ \np{jpjglo}-2& \jp{jpind} $>= 2$        & \jp{jpinf}$<=$\np{jpiglo}-1 \\
469lp\_obc\_north  & $j$-index of a $v$ point      & $i$  of a $T$ point & $i$ of a $T$ point \\
470\hline
471\end{tabular}    \end{center}
472\caption{     \label{Tab_obc_param}
473Names of different indices relating to the open boundaries. In the case
474of a completely open ocean domain with four ocean boundaries, the parameters
475take exactly the values indicated.}
476\end{table}
477%------------------------------------------------------------------------------------------------------------
478
479The open boundaries must be along coordinate lines. On the C-grid, the boundary
480itself is along a line of normal velocity points: $v$ points for a zonal open boundary
481(the south or north one), and $u$ points for a meridional open boundary (the west
482or east one). Another constraint is that there still must be a row of masked points
483all around the domain, as if the domain were a closed basin (unless periodic conditions
484are used together with open boundary conditions). Therefore, an open boundary
485cannot be located at the first/last index, namely, 1, \jp{jpiglo} or \jp{jpjglo}. Also,
486the open boundary algorithm involves calculating the normal velocity points situated
487just on the boundary, as well as the tangential velocity and temperature and salinity
488just outside the boundary. This means that for a west/south boundary, normal
489velocities and temperature are calculated at the same index \jp{jpiwob} and
490\jp{jpjsob}, respectively. For an east/north boundary, the normal velocity is
491calculated at index \jp{jpieob} and \jp{jpjnob}, but the ``outside'' temperature is
492at index \jp{jpieob}+1 and \jp{jpjnob}+1. This means that \jp{jpieob}, \jp{jpjnob} 
493cannot be bigger than \jp{jpiglo}-2, \jp{jpjglo}-2.
494
495
496The starting and ending indices are to be thought of as $T$ point indices: in many
497cases they indicate the first land $T$-point, at the extremity of an open boundary
498(the coast line follows the $f$ grid points, see Fig.~\ref{Fig_obc_north} for an example
499of a northern open boundary). All indices are relative to the global domain. In the
500free surface case it is possible to have ``ocean corners'', that is, an open boundary
501starting and ending in the ocean.
502
503%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
504\begin{figure}[!t]     \begin{center}
505\includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_obc_north.pdf}
506\caption{    \label{Fig_obc_north}
507Localization of the North open boundary points.}
508\end{center}     \end{figure}
509%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
510
511Although not compulsory, it is highly recommended that the bathymetry in the
512vicinity of an open boundary follows the following rule: in the direction perpendicular
513to the open line, the water depth should be constant for 4 grid points. This is in
514order to ensure that the radiation condition, which involves model variables next
515to the boundary, is calculated in a consistent way. On Fig.\ref{Fig_obc_north} we
516indicate by an $=$ symbol, the points which should have the same depth. It means
517that at the 4 points near the boundary, the bathymetry is cylindrical \gmcomment{not sure
518why cylindrical}. The line behind the open $T$-line must be 0 in the bathymetry file
519(as shown on Fig.\ref{Fig_obc_north} for example).
520
521%----------------------------------------------
522\subsection{Boundary data}
523\label{OBC_data}
524
525It is necessary to provide information at the boundaries. The simplest case is
526when this information does not change in time and is equal to the initial conditions
527(namelist variable \np{nn\_obcdta}=0). This is the case for the standard configuration
528EEL5 with open boundaries. When (\np{nn\_obcdta}=1), open boundary information
529is read from netcdf files. For convenience the input files are supposed to be similar
530to the ''history'' NEMO output files, for dimension names and variable names.
531Open boundary arrays must be dimensioned according to the parameters of table~
532\ref{Tab_obc_param}: for example, at the western boundary, arrays have a
533dimension of \jp{jpwf}-\jp{jpwd}+1 in the horizontal and \jp{jpk} in the vertical.
534
535When ocean observations are used to generate the boundary data (a hydrographic
536section for example, as in \citet{Treguier_al_JGR01}) it happens often that only the velocity
537normal to the boundary is known, which is the reason why the initial OBC code
538assumes that only $T$, $S$, and the normal velocity ($u$ or $v$) needs to be
539specified. As more and more global model solutions and ocean analysis products
540become available, it will be possible to provide information about all the variables
541(including the tangential velocity) so that the specification of four variables at each
542boundaries will become standard. For the sea surface height, one must distinguish
543between the filtered free surface case and the time-splitting or explicit treatment of
544the free surface.
545 In the first case, it is assumed that the user does not wish to represent high
546 frequency motions such as tides. The boundary condition is thus one of zero
547 normal gradient of sea surface height at the open boundaries, following \citet{Marchesiello2001}.
548No information other than the total velocity needs to be provided at the open
549boundaries in that case. In the other two cases (time splitting or explicit free surface),
550the user must provide barotropic information (sea surface height and barotropic
551velocities) and the use of the Flather algorithm for barotropic variables is
552recommanded. However, this algorithm has not yet been fully tested and bugs
553remain in NEMO v2.3. Users should read the code carefully before using it. Finally,
554in the case of the rigid lid approximation the barotropic streamfunction must be
555provided, as documented in \citet{Treguier_al_JGR01}). This option is no longer
556recommended but remains in NEMO V2.3.
557
558One frequently encountered case is when an open boundary domain is constructed
559from a global or larger scale NEMO configuration. Assuming the domain corresponds
560to indices $ib:ie$, $jb:je$ of the global domain, the bathymetry and forcing of the
561small domain can be created by using the following netcdf utility on the global files:
562ncks -F $-d\;x,ib,ie$ $-d\;y,jb,je$ (part of the nco series of utilities,
563see their \href{http://nco.sourceforge.net}{website}).
564The open boundary files can be constructed using ncks
565commands, following table~\ref{Tab_obc_ind}.
566
567%--------------------------------------------------TABLE--------------------------------------------------
568\begin{table}[htbp]     \begin{center}      \begin{tabular}{|l|c|c|c|c|c|}
569\hline
570OBC  & Variable   & file name      & Index  & Start  & end  \\
571West &  T,S       &   obcwest\_TS.nc &  $ib$+1     &   $jb$+1 &  $je-1$  \\
572     &    U       &   obcwest\_U.nc  &  $ib$+1     &   $jb$+1 &  $je-1$  \\ 
573     &    V       &   obcwest\_V.nc  &  $ib$+1     &   $jb$+1 &  $je-1$  \\       
574\hline
575East &  T,S       &   obceast\_TS.nc &  $ie$-1     &   $jb$+1 &  $je-1$  \\
576     &    U       &   obceast\_U.nc  &  $ie$-2     &   $jb$+1 &  $je-1$  \\ 
577     &    V       &   obceast\_V.nc  &  $ie$-1     &   $jb$+1 &  $je-1$  \\       
578\hline         
579South &  T,S      &   obcsouth\_TS.nc &  $jb$+1     &  $ib$+1 &  $ie-1$  \\
580      &    U      &   obcsouth\_U.nc  &  $jb$+1     &  $ib$+1 &  $ie-1$  \\ 
581      &    V      &   obcsouth\_V.nc  &  $jb$+1     &  $ib$+1 &  $ie-1$  \\   
582\hline
583North &  T,S      &   obcnorth\_TS.nc &  $je$-1     &  $ib$+1 &  $ie-1$  \\
584      &    U      &   obcnorth\_U.nc  &  $je$-1     &  $ib$+1 &  $ie-1$  \\ 
585      &    V      &   obcnorth\_V.nc  &  $je$-2     &  $ib$+1 &  $ie-1$  \\ 
586\hline
587\end{tabular}     \end{center}
588\caption{    \label{Tab_obc_ind}
589Requirements for creating open boundary files from a global configuration,
590appropriate for the subdomain of indices $ib:ie$, $jb:je$. ``Index'' designates the
591$i$ or $j$ index along which the $u$ of $v$ boundary point is situated in the global
592configuration, starting and ending with the $j$ or $i$ indices indicated.
593For example, to generate file obcnorth\_V.nc, use the command ncks
594$-F$ $-d\;y,je-2$  $-d\;x,ib+1,ie-1$ } 
595\end{table}
596%-----------------------------------------------------------------------------------------------------------
597
598It is assumed that the open boundary files contain the variables for the period of
599the model integration. If the boundary files contain one time frame, the boundary
600data is held fixed in time. If the files contain 12 values, it is assumed that the input
601is a climatology for a repeated annual cycle (corresponding to the case \np{ln\_obc\_clim} 
602=true). The case of an arbitrary number of time frames is not yet implemented
603correctly; the user is required to write his own code in the module \mdl{obc\_dta} 
604to deal with this situation.
605
606\subsection{Radiation algorithm}
607\label{OBC_rad}
608
609The art of open boundary management consists in applying a constraint strong
610enough that the inner domain "feels" the rest of the ocean, but weak enough
611that perturbations are allowed to leave the domain with minimum false reflections
612of energy. The constraints are specified separately at each boundary as time
613scales for ''inflow'' and ''outflow'' as defined below. The time scales are set (in days)
614by namelist parameters such as \np{rn\_dpein}, \np{rn\_dpeob} for the eastern open
615boundary for example. When both time scales are zero for a given boundary
616($e.g.$ for the western boundary, \jp{lp\_obc\_west}=true, \np{rn\_dpwob}=0 and
617\np{rn\_dpwin}=0) this means that the boundary in question is a ''fixed '' boundary
618where the solution is set exactly by the boundary data. This is not recommended,
619except in combination with increased viscosity in a ''sponge'' layer next to the
620boundary in order to avoid spurious reflections. 
621
622
623The radiation\/relaxation \gmcomment{the / doesnt seem to appear in the output} 
624algorithm is applied when either relaxation time (for ''inflow'' or ''outflow'') is
625non-zero. It has been developed and tested in the SPEM model and its
626successor ROMS \citep{Barnier1996, Marchesiello2001}, which is an
627$s$-coordinate model on an Arakawa C-grid. Although the algorithm has
628been numerically successful in the CLIPPER Atlantic models, the physics
629do not work as expected \citep{Treguier_al_JGR01}. Users are invited to consider
630open boundary conditions (OBC hereafter) with some scepticism
631\citep{Durran2001, Blayo2005}.
632
633The first part of the algorithm calculates a phase velocity to determine
634whether perturbations tend to propagate toward, or away from, the
635boundary. Let us consider a model variable $\phi$.
636The phase velocities ($C_{\phi x}$,$C_{\phi y}$) for the variable $\phi$,
637in the directions normal and tangential to the boundary are
638\begin{equation} \label{Eq_obc_cphi}
639C_{\phi x} = \frac{ -\phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{x} 
640\;\;\;\;\; \;\;\; 
641C_{\phi y} = \frac{ -\phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{y}.
642\end{equation}
643Following \citet{Treguier_al_JGR01} and \citet{Marchesiello2001} we retain only
644the normal component of the velocity, $C_{\phi x}$, setting $C_{\phi y} =0$ 
645(but unlike the original Orlanski radiation algorithm we retain $\phi_{y}$ in
646the expression for $C_{\phi x}$). 
647
648The discrete form of (\ref{Eq_obc_cphi}), described by \citet{Barnier1998},
649takes into account the two rows of grid points situated inside the domain
650next to the boundary, and the three previous time steps ($n$, $n-1$,
651and $n-2$). The same equation can then be discretized at the boundary at
652time steps $n-1$, $n$ and $n+1$ \gmcomment{since the original was three time-level} 
653in order to extrapolate for the new boundary value $\phi^{n+1}$.
654
655In the open boundary algorithm as implemented in NEMO v2.3, the new boundary
656values are updated differently depending on the sign of $C_{\phi x}$. Let us take
657an eastern boundary as an example. The solution for variable $\phi$ at the
658boundary is given by a generalized wave equation with phase velocity $C_{\phi}$,
659with the addition of a relaxation term, as:
660\begin{eqnarray}
661\phi_{t} &  =  & -C_{\phi x} \phi_{x} + \frac{1}{\tau_{o}} (\phi_{c}-\phi)
662                        \;\;\; \;\;\; \;\;\; (C_{\phi x} > 0), \label{Eq_obc_rado} \\
663\phi_{t} &  =  & \frac{1}{\tau_{i}} (\phi_{c}-\phi)
664\;\;\; \;\;\; \;\;\;\;\;\; (C_{\phi x} < 0), \label{Eq_obc_radi}
665\end{eqnarray}
666where $\phi_{c}$ is the estimate of $\phi$ at the boundary, provided as boundary
667data. Note that in (\ref{Eq_obc_rado}), $C_{\phi x}$ is bounded by the ratio
668$\delta x/\delta t$ for stability reasons. When $C_{\phi x}$ is eastward (outward
669propagation), the radiation condition (\ref{Eq_obc_rado}) is used.
670When  $C_{\phi x}$ is westward (inward propagation), (\ref{Eq_obc_radi}) is
671used with a strong relaxation to climatology (usually $\tau_{i}=\np{rn\_dpein}=$1~day).
672Equation (\ref{Eq_obc_radi}) is solved with a Euler time-stepping scheme. As a
673consequence, setting $\tau_{i}$ smaller than, or equal to the time step is equivalent
674to a fixed boundary condition. A time scale of one day is usually a good compromise
675which guarantees that the inflow conditions remain close to climatology while ensuring
676numerical stability.
677
678In  the case of a western boundary located in the Eastern Atlantic, \citet{Penduff_al_JGR00} 
679have been able to implement the radiation algorithm without any boundary data,
680using persistence from the previous time step instead. This solution has not worked
681in other cases \citep{Treguier_al_JGR01}, so that the use of boundary data is recommended.
682Even in the outflow condition (\ref{Eq_obc_rado}), we have found it desirable to
683maintain a weak relaxation to climatology. The time step is usually chosen so as to
684be larger than typical turbulent scales (of order 1000~days \gmcomment{or maybe seconds?}).
685
686The radiation condition is applied to the model variables: temperature, salinity,
687tangential and normal velocities. For normal and tangential velocities, $u$ and $v$,
688radiation is applied with phase velocities calculated from $u$ and $v$ respectively. 
689For the radiation of tracers, we use the phase velocity calculated from the tangential
690velocity in order to avoid calculating too many independent radiation velocities and
691because tangential velocities and tracers have the same position along the boundary
692on a C-grid. 
693
694\subsection{Domain decomposition (\key{mpp\_mpi})}
695\label{OBC_mpp}
696When \key{mpp\_mpi} is active in the code, the computational domain is divided
697into rectangles that are attributed each to a different processor. The open boundary
698code is ``mpp-compatible'' up to a certain point. The radiation algorithm will not
699work if there is an mpp subdomain boundary parallel to the open boundary at the
700index of the boundary, or the grid point after (outside), or three grid points before
701(inside). On the other hand, there is no problem if an mpp subdomain boundary
702cuts the open boundary perpendicularly. These geometrical limitations must be
703checked for by the user (there is no safeguard in the code). 
704The general principle for the open boundary mpp code is that loops over the open
705boundaries {not sure what this means} are performed on local indices (nie0,
706nie1, nje0, nje1 for an eastern boundary for instance) that are initialized in module
707\mdl{obc\_ini}. Those indices have relevant values on the processors that contain
708a segment of an open boundary. For processors that do not include an open
709boundary segment, the indices are such that the calculations within the loops are
710not performed.
711\gmcomment{I dont understand most of the last few sentences}
712 
713Arrays of climatological data that are read from files are seen by all processors
714and have the same dimensions for all (for instance, for the eastern boundary,
715uedta(jpjglo,jpk,2)). On the other hand, the arrays for the calculation of radiation
716are local to each processor (uebnd(jpj,jpk,3,3) for instance).  This allowed the
717CLIPPER model for example, to save on memory where the eastern boundary
718crossed 8 processors so that \jp{jpj} was much smaller than (\jp{jpjef}-\jp{jpjed}+1).
719
720\subsection{Volume conservation}
721\label{OBC_vol}
722
723It is necessary to control the volume inside a domain when using open boundaries.
724With fixed boundaries, it is enough to ensure that the total inflow/outflow has
725reasonable values (either zero or a value compatible with an observed volume
726balance). When using radiative boundary conditions it is necessary to have a
727volume constraint because each open boundary works independently from the
728others. The methodology used to control this volume is identical to the one
729coded in the ROMS model \citep{Marchesiello2001}.
730
731
732%---------------------------------------- EXTRAS
733\colorbox{yellow}{Explain obc\_vol{\ldots}}
734
735\colorbox{yellow}{OBC algorithm for update, OBC restart, list of routines where obc key appears{\ldots}}
736
737\colorbox{yellow}{OBC rigid lid? {\ldots}}
738
739% ====================================================================
740% Unstructured open boundaries BDY
741% ====================================================================
742\section{Unstructured Open Boundary Conditions (\key{bdy}) (BDY)}
743\label{LBC_bdy}
744
745%-----------------------------------------nambdy--------------------------------------------
746\namdisplay{nambdy} 
747%-----------------------------------------------------------------------------------------------
748%-----------------------------------------nambdy_index--------------------------------------------
749\namdisplay{nambdy_index} 
750%-----------------------------------------------------------------------------------------------
751%-----------------------------------------nambdy_dta--------------------------------------------
752\namdisplay{nambdy_dta} 
753%-----------------------------------------------------------------------------------------------
754%-----------------------------------------nambdy_dta--------------------------------------------
755\namdisplay{nambdy_dta2} 
756%-----------------------------------------------------------------------------------------------
757
758Options are defined through the \ngn{nambdy} \ngn{nambdy\_index} 
759\ngn{nambdy\_dta} \ngn{nambdy\_dta2} namelist variables.
760The BDY module is an alternative implementation of open boundary
761conditions for regional configurations. It implements the Flow
762Relaxation Scheme algorithm for temperature, salinity, velocities and
763ice fields, and the Flather radiation condition for the depth-mean
764transports. The specification of the location of the open boundary is
765completely flexible and allows for example the open boundary to follow
766an isobath or other irregular contour.
767
768The BDY module was modelled on the OBC module and shares many features
769and a similar coding structure \citep{Chanut2005}.
770
771The BDY module is completely rewritten at NEMO 3.4 and there is a new
772set of namelists. Boundary data files used with earlier versions of
773NEMO may need to be re-ordered to work with this version. See the
774section on the Input Boundary Data Files for details.
775
776%----------------------------------------------
777\subsection{The namelists}
778\label{BDY_namelist}
779
780It is possible to define more than one boundary ``set'' and apply
781different boundary conditions to each set. The number of boundary
782sets is defined by \np{nb\_bdy}.  Each boundary set may be defined
783as a set of straight line segments in a namelist
784(\np{ln\_coords\_file}=.false.) or read in from a file
785(\np{ln\_coords\_file}=.true.). If the set is defined in a namelist,
786then the namelists nambdy\_index must be included separately, one for
787each set. If the set is defined by a file, then a
788``coordinates.bdy.nc'' file must be provided. The coordinates.bdy file
789is analagous to the usual NEMO ``coordinates.nc'' file. In the example
790above, there are two boundary sets, the first of which is defined via
791a file and the second is defined in a namelist. For more details of
792the definition of the boundary geometry see section
793\ref{BDY_geometry}.
794
795For each boundary set a boundary
796condition has to be chosen for the barotropic solution (``u2d'':
797sea-surface height and barotropic velocities), for the baroclinic
798velocities (``u3d''), and for the active tracers\footnote{The BDY
799  module does not deal with passive tracers at this version}
800(``tra''). For each set of variables there is a choice of algorithm
801and a choice for the data, eg. for the active tracers the algorithm is
802set by \np{nn\_tra} and the choice of data is set by
803\np{nn\_tra\_dta}.
804
805The choice of algorithm is currently as follows:
806
807\mbox{}
808
809\begin{itemize}
810\item[0.] No boundary condition applied. So the solution will ``see''
811  the land points around the edge of the edge of the domain.
812\item[1.] Flow Relaxation Scheme (FRS) available for all variables.
813\item[2.] Flather radiation scheme for the barotropic variables. The
814  Flather scheme is not compatible with the filtered free surface
815  ({\it dynspg\_ts}).
816\end{itemize}
817
818\mbox{}
819
820The main choice for the boundary data is
821to use initial conditions as boundary data (\np{nn\_tra\_dta}=0) or to
822use external data from a file (\np{nn\_tra\_dta}=1). For the
823barotropic solution there is also the option to use tidal
824harmonic forcing either by itself or in addition to other external
825data.
826
827If external boundary data is required then the nambdy\_dta namelist
828must be defined. One nambdy\_dta namelist is required for each boundary
829set in the order in which the boundary sets are defined in nambdy. In
830the example given, two boundary sets have been defined and so there
831are two nambdy\_dta namelists. The boundary data is read in using the
832fldread module, so the nambdy\_dta namelist is in the format required
833for fldread. For each variable required, the filename, the frequency
834of the files and the frequency of the data in the files is given. Also
835whether or not time-interpolation is required and whether the data is
836climatological (time-cyclic) data. Note that on-the-fly spatial
837interpolation of boundary data is not available at this version.
838
839In the example namelists given, two boundary sets are defined. The
840first set is defined via a file and applies FRS conditions to
841temperature and salinity and Flather conditions to the barotropic
842variables. External data is provided in daily files (from a
843large-scale model). Tidal harmonic forcing is also used. The second
844set is defined in a namelist. FRS conditions are applied on
845temperature and salinity and climatological data is read from external
846files.
847
848%----------------------------------------------
849\subsection{The Flow Relaxation Scheme}
850\label{BDY_FRS_scheme}
851
852The Flow Relaxation Scheme (FRS) \citep{Davies_QJRMS76,Engerdahl_Tel95},
853applies a simple relaxation of the model fields to
854externally-specified values over a zone next to the edge of the model
855domain. Given a model prognostic variable $\Phi$ 
856\begin{equation}  \label{Eq_bdy_frs1}
857\Phi(d) = \alpha(d)\Phi_{e}(d) + (1-\alpha(d))\Phi_{m}(d)\;\;\;\;\; d=1,N
858\end{equation}
859where $\Phi_{m}$ is the model solution and $\Phi_{e}$ is the specified
860external field, $d$ gives the discrete distance from the model
861boundary  and $\alpha$ is a parameter that varies from $1$ at $d=1$ to
862a small value at $d=N$. It can be shown that this scheme is equivalent
863to adding a relaxation term to the prognostic equation for $\Phi$ of
864the form:
865\begin{equation}  \label{Eq_bdy_frs2}
866-\frac{1}{\tau}\left(\Phi - \Phi_{e}\right)
867\end{equation}
868where the relaxation time scale $\tau$ is given by a function of
869$\alpha$ and the model time step $\Delta t$:
870\begin{equation}  \label{Eq_bdy_frs3}
871\tau = \frac{1-\alpha}{\alpha}  \,\rdt
872\end{equation}
873Thus the model solution is completely prescribed by the external
874conditions at the edge of the model domain and is relaxed towards the
875external conditions over the rest of the FRS zone. The application of
876a relaxation zone helps to prevent spurious reflection of outgoing
877signals from the model boundary.
878
879The function $\alpha$ is specified as a $tanh$ function:
880\begin{equation}  \label{Eq_bdy_frs4}
881\alpha(d) = 1 - \tanh\left(\frac{d-1}{2}\right),       \quad d=1,N
882\end{equation}
883The width of the FRS zone is specified in the namelist as
884\np{nn\_rimwidth}. This is typically set to a value between 8 and 10.
885
886%----------------------------------------------
887\subsection{The Flather radiation scheme}
888\label{BDY_flather_scheme}
889
890The \citet{Flather_JPO94} scheme is a radiation condition on the normal, depth-mean
891transport across the open boundary. It takes the form
892\begin{equation}  \label{Eq_bdy_fla1}
893U = U_{e} + \frac{c}{h}\left(\eta - \eta_{e}\right),
894\end{equation}
895where $U$ is the depth-mean velocity normal to the boundary and $\eta$
896is the sea surface height, both from the model. The subscript $e$
897indicates the same fields from external sources. The speed of external
898gravity waves is given by $c = \sqrt{gh}$, and $h$ is the depth of the
899water column. The depth-mean normal velocity along the edge of the
900model domain is set equal to the
901external depth-mean normal velocity, plus a correction term that
902allows gravity waves generated internally to exit the model boundary.
903Note that the sea-surface height gradient in \eqref{Eq_bdy_fla1}
904is a spatial gradient across the model boundary, so that $\eta_{e}$ is
905defined on the $T$ points with $nbr=1$ and $\eta$ is defined on the
906$T$ points with $nbr=2$. $U$ and $U_{e}$ are defined on the $U$ or
907$V$ points with $nbr=1$, $i.e.$ between the two $T$ grid points.
908
909%----------------------------------------------
910\subsection{Boundary geometry}
911\label{BDY_geometry}
912
913Each open boundary set is defined as a list of points. The information
914is stored in the arrays $nbi$, $nbj$, and $nbr$ in the $idx\_bdy$
915structure.  The $nbi$ and $nbj$ arrays
916define the local $(i,j)$ indices of each point in the boundary zone
917and the $nbr$ array defines the discrete distance from the boundary
918with $nbr=1$ meaning that the point is next to the edge of the
919model domain and $nbr>1$ showing that the point is increasingly
920further away from the edge of the model domain. A set of $nbi$, $nbj$,
921and $nbr$ arrays is defined for each of the $T$, $U$ and $V$
922grids. Figure \ref{Fig_LBC_bdy_geom} shows an example of an irregular
923boundary.
924
925The boundary geometry for each set may be defined in a namelist
926nambdy\_index or by reading in a ``coordinates.bdy.nc'' file. The
927nambdy\_index namelist defines a series of straight-line segments for
928north, east, south and west boundaries. For the northern boundary,
929\np{nbdysegn} gives the number of segments, \np{jpjnob} gives the $j$
930index for each segment and \np{jpindt} and \np{jpinft} give the start
931and end $i$ indices for each segment with similar for the other
932boundaries. These segments define a list of $T$ grid points along the
933outermost row of the boundary ($nbr\,=\, 1$). The code deduces the $U$ and
934$V$ points and also the points for $nbr\,>\, 1$ if
935$nn\_rimwidth\,>\,1$.
936
937The boundary geometry may also be defined from a
938``coordinates.bdy.nc'' file. Figure \ref{Fig_LBC_nc_header}
939gives an example of the header information from such a file. The file
940should contain the index arrays for each of the $T$, $U$ and $V$
941grids. The arrays must be in order of increasing $nbr$. Note that the
942$nbi$, $nbj$ values in the file are global values and are converted to
943local values in the code. Typically this file will be used to generate
944external boundary data via interpolation and so will also contain the
945latitudes and longitudes of each point as shown. However, this is not
946necessary to run the model.
947
948For some choices of irregular boundary the model domain may contain
949areas of ocean which are not part of the computational domain. For
950example if an open boundary is defined along an isobath, say at the
951shelf break, then the areas of ocean outside of this boundary will
952need to be masked out. This can be done by reading a mask file defined
953as \np{cn\_mask\_file} in the nam\_bdy namelist. Only one mask file is
954used even if multiple boundary sets are defined.
955
956%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
957\begin{figure}[!t]      \begin{center}
958\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_LBC_bdy_geom.pdf}
959\caption {      \label{Fig_LBC_bdy_geom}
960Example of geometry of unstructured open boundary}
961\end{center}   \end{figure}
962%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
963
964%----------------------------------------------
965\subsection{Input boundary data files}
966\label{BDY_data}
967
968The data files contain the data arrays
969in the order in which the points are defined in the $nbi$ and $nbj$
970arrays. The data arrays are dimensioned on: a time dimension;
971$xb$ which is the index of the boundary data point in the horizontal;
972and $yb$ which is a degenerate dimension of 1 to enable the file to be
973read by the standard NEMO I/O routines. The 3D fields also have a
974depth dimension.
975
976At Version 3.4 there are new restrictions on the order in which the
977boundary points are defined (and therefore restrictions on the order
978of the data in the file). In particular:
979
980\mbox{}
981
982\begin{enumerate}
983\item The data points must be in order of increasing $nbr$, ie. all
984  the $nbr=1$ points, then all the $nbr=2$ points etc.
985\item All the data for a particular boundary set must be in the same
986  order. (Prior to 3.4 it was possible to define barotropic data in a
987  different order to the data for tracers and baroclinic velocities).
988\end{enumerate}
989
990\mbox{}
991
992These restrictions mean that data files used with previous versions of
993the model may not work with version 3.4. A fortran utility
994{\it bdy\_reorder} exists in the TOOLS directory which will re-order the
995data in old BDY data files.
996
997%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
998\begin{figure}[!t]     \begin{center}
999\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_LBC_nc_header.pdf}
1000\caption {     \label{Fig_LBC_nc_header} 
1001Example of the header for a coordinates.bdy.nc file}
1002\end{center}   \end{figure}
1003%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1004
1005%----------------------------------------------
1006\subsection{Volume correction}
1007\label{BDY_vol_corr}
1008
1009There is an option to force the total volume in the regional model to be constant,
1010similar to the option in the OBC module. This is controlled  by the \np{nn\_volctl} 
1011parameter in the namelist. A value of \np{nn\_volctl}~=~0 indicates that this option is not used.
1012If  \np{nn\_volctl}~=~1 then a correction is applied to the normal velocities
1013around the boundary at each timestep to ensure that the integrated volume flow
1014through the boundary is zero. If \np{nn\_volctl}~=~2 then the calculation of
1015the volume change on the timestep includes the change due to the freshwater
1016flux across the surface and the correction velocity corrects for this as well.
1017
1018If more than one boundary set is used then volume correction is
1019applied to all boundaries at once.
1020
1021\newpage
1022%----------------------------------------------
1023\subsection{Tidal harmonic forcing}
1024\label{BDY_tides}
1025
1026%-----------------------------------------nambdy_tide--------------------------------------------
1027\namdisplay{nambdy_tide} 
1028%-----------------------------------------------------------------------------------------------
1029
1030Options are defined through the  \ngn{nambdy\_tide} namelist variables.
1031 To be written....
1032
1033
1034
1035
Note: See TracBrowser for help on using the repository browser.