New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Chap_LBC.tex in tags/nemo_v3_2/nemo_v3_2/DOC/TexFiles/Chapters – NEMO

source: tags/nemo_v3_2/nemo_v3_2/DOC/TexFiles/Chapters/Chap_LBC.tex @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm

File size: 44.1 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{shlat})}
18\label{LBC_coast}
19%--------------------------------------------nam_lbc-------------------------------------------------------
20\namdisplay{nam_lbc} 
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
27The discrete representation of a domain with complex boundaries (coastlines and
28bottom topography) leads to arrays that include large portions where a computation
29is not required as the model variables remain at zero. Nevertheless, vectorial
30supercomputers are far more efficient when computing over a whole array, and the
31readability of a code is greatly improved when boundary conditions are applied in
32an automatic way rather than by a specific computation before or after each
33computational loop. An efficient way to work over the whole domain while specifying
34the boundary conditions, is to use multiplication by mask arrays in the computation.
35A mask array is a matrix whose elements are $1$ in the ocean domain and $0$ 
36elsewhere. A simple multiplication of a variable by its own mask ensures that it will
37remain zero over land areas. Since most of the boundary conditions consist of a
38zero flux across the solid boundaries, they can be simply applied by multiplying
39variables by the correct mask arrays, $i.e.$ the mask array of the grid point where
40the flux is evaluated. For example, the heat flux in the \textbf{i}-direction is evaluated
41at $u$-points. Evaluating this quantity as,
42
43\begin{equation} \label{Eq_lbc_aaaa}
44\frac{A^{lT} }{e_1 }\frac{\partial T}{\partial i}\equiv \frac{A_u^{lT} 
45}{e_{1u} } \; \delta _{i+1 / 2} \left[ T \right]\;\;mask_u
46\end{equation}
47(where mask$_{u}$ is the mask array at a $u$-point) ensures that the heat flux is
48zero inside land and at the boundaries, since mask$_{u}$ is zero at solid boundaries
49which in this case are defined at $u$-points (normal velocity $u$ remains zero at
50the coast) (Fig.~\ref{Fig_LBC_uv}).
51
52%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
53\begin{figure}[!t] \label{Fig_LBC_uv}  \begin{center}
54\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_LBC_uv.pdf}
55\caption {Lateral boundary (thick line) at T-level. The velocity normal to the
56       boundary is set to zero.}
57\end{center}   \end{figure}
58%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
59
60For momentum the situation is a bit more complex as two boundary conditions
61must be provided along the coast (one each for the normal and tangential velocities).
62The boundary of the ocean in the C-grid is defined by the velocity-faces.
63For example, at a given $T$-level, the lateral boundary (a coastline or an intersection
64with the bottom topography) is made of segments joining $f$-points, and normal
65velocity points are located between two $f-$points (Fig.~\ref{Fig_LBC_uv}).
66The boundary condition on the normal velocity (no flux through solid boundaries)
67can thus be easily implemented using the mask system. The boundary condition
68on the tangential velocity requires a more specific treatment. This boundary
69condition influences the relative vorticity and momentum diffusive trends, and is
70required in order to compute the vorticity at the coast. Four different types of
71lateral boundary condition are available, controlled by the value of the \np{shlat} 
72namelist parameter. (The value of the mask$_{f}$ array along the coastline is set
73equal to this parameter.) These are:
74
75%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
76\begin{figure}[!p] \label{Fig_LBC_shlat}  \begin{center}
77\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_LBC_shlat.pdf}
78\caption {lateral boundary condition (a) free-slip ($shlat=0$) ; (b) no-slip ($shlat=2$)
79; (c) "partial" free-slip ($0<shlat<2$) and (d) "strong" no-slip ($2<shlat$).
80Implied "ghost" velocity inside land area is display in grey. }
81\end{center}   \end{figure}
82%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
83
84\begin{description}
85
86\item[free-slip boundary condition (\np{shlat}=0): ]  the tangential velocity at the
87coastline is equal to the offshore velocity, $i.e.$ the normal derivative of the
88tangential velocity is zero at the coast, so the vorticity: mask$_{f}$ array is set
89to zero inside the land and just at the coast (Fig.~\ref{Fig_LBC_shlat}-a).
90
91\item[no-slip boundary condition (\np{shlat}=2): ] the tangential velocity vanishes
92at the coastline. Assuming that the tangential velocity decreases linearly from
93the closest ocean velocity grid point to the coastline, the normal derivative is
94evaluated as if the velocities at the closest land velocity gridpoint and the closest
95ocean velocity gridpoint were of the same magnitude but in the opposite direction
96(Fig.~\ref{Fig_LBC_shlat}-b). Therefore, the vorticity along the coastlines is given by:
97
98\begin{equation*}
99\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) \ ,
100\end{equation*}
101where $u$ and $v$ are masked fields. Setting the mask$_{f}$ array to $2$ along
102the coastline provides a vorticity field computed with the no-slip boundary condition,
103simply by multiplying it by the mask$_{f}$ :
104\begin{equation} \label{Eq_lbc_bbbb}
105\zeta \equiv \frac{1}{e_{1f} {\kern 1pt}e_{2f} }\left( {\delta _{i+1/2} 
106\left[ {e_{2v} \,v} \right]-\delta _{j+1/2} \left[ {e_{1u} \,u} \right]} 
107\right)\;\mbox{mask}_f
108\end{equation}
109
110\item["partial" free-slip boundary condition (0$<$\np{shlat}$<$2): ] the tangential
111velocity at the coastline is smaller than the offshore velocity, $i.e.$ there is a lateral
112friction but not strong enough to make the tangential velocity at the coast vanish
113(Fig.~\ref{Fig_LBC_shlat}-c). This can be selected by providing a value of mask$_{f}$ 
114strictly inbetween $0$ and $2$.
115
116\item["strong" no-slip boundary condition (2$<$\np{shlat}): ] the viscous boundary
117layer is assumed to be smaller than half the grid size (Fig.~\ref{Fig_LBC_shlat}-d).
118The friction is thus larger than in the no-slip case.
119
120\end{description}
121
122Note that when the bottom topography is entirely represented by the $s$-coor-dinates
123(pure $s$-coordinate), the lateral boundary condition on tangential velocity is of much
124less importance as it is only applied next to the coast where the minimum water depth
125can be quite shallow.
126
127The alternative numerical implementation of the no-slip boundary conditions for an
128arbitrary coast line of \citet{Shchepetkin1996} is also available through the
129\key{noslip\_accurate} CPP key. It is based on a fourth order evaluation of the shear at the
130coast which, in turn, allows a true second order scheme in the interior of the domain
131($i.e.$ the numerical boundary scheme simulates the truncation error of the numerical
132scheme used in the interior of the domain). \citet{Shchepetkin1996} found that such a
133technique considerably improves the quality of the numerical solution. In \NEMO, such
134spectacular improvements have not been found in the half-degree global ocean
135(ORCA05), but significant reductions of numerically induced coastal upwellings were
136found in an eddy resolving simulation of the Alboran Sea \citep{OlivierPh2001}.
137Nevertheless, since a no-slip boundary condition is not recommended in an eddy
138permitting or resolving simulation \citep{Penduff2007}, the use of this option is also
139not recommended.
140
141In practice, the no-slip accurate option changes the way the curl is evaluated at the
142coast (see \mdl{divcur} module), and requires the nature of each coastline grid point
143(convex or concave corners, straight north-south or east-west coast) to be specified. 
144This is performed in routine \rou{dom\_msk\_nsa} in the \mdl{domask} module.
145
146% ================================================================
147% Boundary Condition around the Model Domain
148% ================================================================
149\section{Model Domain Boundary Condition (\jp{jperio})}
150\label{LBC_jperio}
151
152At the model domain boundaries several choices are offered: closed, cyclic east-west,
153south symmetric across the equator, a north-fold, and combination closed-north fold
154or cyclic-north-fold. The north-fold boundary condition is associated with the 3-pole ORCA mesh.
155
156% -------------------------------------------------------------------------------------------------------------
157%        Closed, cyclic, south symmetric (\jp{jperio} = 0, 1 or 2)
158% -------------------------------------------------------------------------------------------------------------
159\subsection{Closed, cyclic, south symmetric (\jp{jperio} = 0, 1 or 2)}
160\label{LBC_jperio012}
161
162The choice of closed, cyclic or symmetric model domain boundary condition is made
163by setting \jp{jperio} to 0, 1 or 2 in file \mdl{par\_oce}. Each time such a boundary
164condition is needed, it is set by a call to routine \mdl{lbclnk}. The computation of
165momentum and tracer trends proceeds from $i=2$ to $i=jpi-1$ and from $j=2$ to
166$j=jpj-1$, $i.e.$ in the model interior. To choose a lateral model boundary condition
167is to specify the first and last rows and columns of the model variables.
168
169\begin{description}
170
171\item[For closed boundary (\textit{jperio=0})], solid walls are imposed at all model
172boundaries: first and last rows and columns are set to zero.
173
174\item[For cyclic east-west boundary (\textit{jperio=1})], first and last rows are set
175to zero (closed) whilst the first column is set to the value of the last-but-one column
176and the last column to the value of the second one (Fig.~\ref{Fig_LBC_jperio}-a).
177Whatever flows out of the eastern (western) end of the basin enters the western
178(eastern) end. Note that there is no option for north-south cyclic or for doubly
179cyclic cases.
180
181\item[For symmetric boundary condition across the equator (\textit{jperio=2})],
182last rows, and first and last columns are set to zero (closed). The row of symmetry
183is chosen to be the $u$- and $T-$points equator line ($j=2$, i.e. at the southern
184end of the domain). For arrays defined at $u-$ or $T-$points, the first row is set
185to the value of the third row while for most of $v$- and $f$-point arrays ($v$, $\zeta$,
186$j\psi$, but \gmcomment{not sure why this is "but"} scalar arrays such as eddy coefficients)
187the first row is set to minus the value of the second row (Fig.~\ref{Fig_LBC_jperio}-b).
188Note that this boundary condition is not yet available for the case of a massively
189parallel computer (\textbf{key{\_}mpp} defined).
190
191\end{description}
192
193%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
194\begin{figure}[!t] \label{Fig_LBC_jperio}  \begin{center}
195\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_LBC_jperio.pdf}
196\caption {setting of (a) east-west cyclic  (b) symmetric across the equator boundary conditions.}
197\end{center}   \end{figure}
198%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
199
200% -------------------------------------------------------------------------------------------------------------
201%        North fold (\textit{jperio = 3 }to $6)$
202% -------------------------------------------------------------------------------------------------------------
203\subsection{North-fold (\textit{jperio = 3 }to $6)$ }
204\label{LBC_north_fold}
205
206The north fold boundary condition has been introduced in order to handle the north
207boundary of a three-polar ORCA grid. Such a grid has two poles in the northern hemisphere.
208\colorbox{yellow}{to be completed...}
209
210%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
211\begin{figure}[!t] \label{Fig_North_Fold_T}  \begin{center}
212\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_North_Fold_T.pdf}
213\caption {North fold boundary with a $T$-point pivot and cyclic east-west boundary condition
214($jperio=4$), as used in ORCA 2, 1/4, and 1/12. Pink shaded area corresponds to the inner
215domain mask (see text). }
216\end{center}   \end{figure}
217%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
218
219% ====================================================================
220% Exchange with neighbouring processors
221% ====================================================================
222\section  [Exchange with neighbouring processors (\textit{lbclnk}, \textit{lib\_mpp})]
223      {Exchange with neighbouring processors (\mdl{lbclnk}, \mdl{lib\_mpp})}
224\label{LBC_mpp}
225
226For massively parallel processing (mpp), a domain decomposition method is used.
227The basic idea of the method is to split the large computation domain of a numerical
228experiment into several smaller domains and solve the set of equations by addressing
229independent local problems. Each processor has its own local memory and computes
230the model equation over a subdomain of the whole model domain. The subdomain
231boundary conditions are specified through communications between processors
232which are organized by explicit statements (message passing method).
233
234A big advantage is that the method does not need many modifications of the initial
235FORTRAN code. From the modeller's point of view, each sub domain running on
236a processor is identical to the "mono-domain" code. In addition, the programmer
237manages the communications between subdomains, and the code is faster when
238the number of processors is increased. The porting of OPA code on an iPSC860
239was achieved during Guyon's PhD [Guyon et al. 1994, 1995] in collaboration with
240CETIIS and ONERA. The implementation in the operational context and the studies
241of performance on a T3D and T3E Cray computers have been made in collaboration
242with IDRIS and CNRS. The present implementation is largely inspired by Guyon's
243work  [Guyon 1995].
244
245The parallelization strategy is defined by the physical characteristics of the
246ocean model. Second order finite difference schemes lead to local discrete
247operators that depend at the very most on one neighbouring point. The only
248non-local computations concern the vertical physics (implicit diffusion, 1.5
249turbulent closure scheme, ...) (delocalization over the whole water column),
250and the solving of the elliptic equation associated with the surface pressure
251gradient computation (delocalization over the whole horizontal domain).
252Therefore, a pencil strategy is used for the data sub-structuration
253\gmcomment{no idea what this means!}
254: the 3D initial domain is laid out on local processor
255memories following a 2D horizontal topological splitting. Each sub-domain
256computes its own surface and bottom boundary conditions and has a side
257wall overlapping interface which defines the lateral boundary conditions for
258computations in the inner sub-domain. The overlapping area consists of the
259two rows at each edge of the sub-domain. After a computation, a communication
260phase starts: each processor sends to its neighbouring processors the update
261values of the points corresponding to the interior overlapping area to its
262neighbouring sub-domain (i.e. the innermost of the two overlapping rows).
263The communication is done through message passing. Usually the parallel virtual
264language, PVM, is used as it is a standard language available on  nearly  all
265MPP computers. More specific languages (i.e. computer dependant languages)
266can be easily used to speed up the communication, such as SHEM on a T3E
267computer. The data exchanges between processors are required at the very
268place where lateral domain boundary conditions are set in the mono-domain
269computation (\S III.10-c): the lbc\_lnk routine which manages such conditions
270is substituted by mpplnk.F or mpplnk2.F routine when running on an MPP
271computer (\key{mpp\_mpi} defined). It has to be pointed out that when using
272the MPP version of the model, the east-west cyclic boundary condition is done
273implicitly, whilst the south-symmetric boundary condition option is not available.
274
275%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
276\begin{figure}[!t] \label{Fig_mpp}  \begin{center}
277\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_mpp.pdf}
278\caption {Positioning of a sub-domain when massively parallel processing is used. }
279\end{center}   \end{figure}
280%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
281
282In the standard version of the OPA model, the splitting is regular and arithmetic.
283 the i-axis is divided by \jp{jpni} and the j-axis by \jp{jpnj} for a number of processors
284 \jp{jpnij} most often equal to $jpni \times jpnj$ (model parameters set in
285 \mdl{par\_oce}). Each processor is independent and without message passing
286 or synchronous process
287 \gmcomment{how does a synchronous process relate to this?},
288 programs run alone and access just its own local memory. For this reason, the
289 main model dimensions are now the local dimensions of the subdomain (pencil)
290 that are named \jp{jpi}, \jp{jpj}, \jp{jpk}. These dimensions include the internal
291 domain and the overlapping rows. The number of rows to exchange (known as
292 the halo) is usually set to one (\jp{jpreci}=1, in \mdl{par\_oce}). The whole domain
293 dimensions are named \jp{jpiglo}, \jp{jpjglo} and \jp{jpk}. The relationship between
294 the whole domain and a sub-domain is:
295\begin{eqnarray} 
296      jpi & = & ( jpiglo-2*jpreci + (jpni-1) ) / jpni + 2*jpreci  \nonumber \\
297      jpj & = & ( jpjglo-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj  \label{Eq_lbc_jpi}
298\end{eqnarray}
299where \jp{jpni}, \jp{jpnj} are the number of processors following the i- and j-axis.
300
301\colorbox{yellow}{Figure IV.3: example of a domain splitting with 9 processors and
302no east-west cyclic boundary conditions.}
303
304One also defines variables nldi and nlei which correspond to the internal
305domain bounds, and the variables nimpp and njmpp which are the position
306of the (1,1) grid-point in the global domain. An element of $T_{l}$, a local array
307(subdomain) corresponds to an element of $T_{g}$, a global array
308(whole domain) by the relationship:
309\begin{equation} \label{Eq_lbc_nimpp}
310T_{g} (i+nimpp-1,j+njmpp-1,k) = T_{l} (i,j,k),
311\end{equation}
312with  $1 \leq i \leq jpi$, $1  \leq j \leq jpj $ , and  $1  \leq k \leq jpk$.
313
314Processors are numbered from 0 to $jpnij-1$, the number is saved in the variable
315nproc. In the standard version, a processor has no more than four neighbouring
316processors named nono (for north), noea (east), noso (south) and nowe (west)
317and two variables, nbondi and nbondj, indicate the relative position of the processor
318\colorbox{yellow}{(see Fig.IV.3)}:
319\begin{itemize}
320\item       nbondi = -1    an east neighbour, no west processor,
321\item       nbondi =  0 an east neighbour, a west neighbour,
322\item       nbondi =  1    no east processor, a west neighbour,
323\item       nbondi =  2    no splitting following the i-axis.
324\end{itemize}
325During the simulation, processors exchange data with their neighbours.
326If there is effectively a neighbour, the processor receives variables from this
327processor on its overlapping row, and sends the data issued from internal
328domain corresponding to the overlapping row of the other processor.
329       
330\colorbox{yellow}{Figure IV.4: pencil splitting with the additional outer halos }
331
332
333The \NEMO model computes equation terms with the help of mask arrays (0 on land
334points and 1 on sea points). It is easily readable and very efficient in the context of
335a computer with vectorial architecture. However, in the case of a scalar processor,
336computations over the land regions become more expensive in terms of CPU time.
337It is worse when we use a complex configuration with a realistic bathymetry like the
338global ocean where more than 50 \% of points are land points. For this reason, a
339pre-processing tool can be used to choose the mpp domain decomposition with a
340maximum number of only land points processors, which can then be eliminated.
341(For example, the mpp\_optimiz tools, available from the DRAKKAR web site.)
342This optimisation is dependent on the specific bathymetry employed. The user
343then chooses optimal parameters \jp{jpni}, \jp{jpnj} and \jp{jpnij} with
344$jpnij < jpni \times jpnj$, leading to the elimination of $jpni \times jpnj - jpnij$ 
345land processors. When those parameters are specified in module \mdl{par\_oce},
346the algorithm in the \rou{inimpp2} routine sets each processor's parameters (nbound,
347nono, noea,...) so that the land-only processors are not taken into account.
348
349\colorbox{yellow}{Note that the inimpp2 routine is general so that the original inimpp
350routine should be suppressed from the code.}
351
352When land processors are eliminated, the value corresponding to these locations in
353the model output files is zero. Note that this is a problem for a mesh output file written
354by such a model configuration, because model users often divide by the scale factors
355($e1t$, $e2t$, etc) and do not expect the grid size to be zero, even on land. It may be
356best not to eliminate land processors when running the model especially to write the
357mesh files as outputs (when \np{nmsh} namelist parameter differs from 0).
358\gmcomment{Steven : dont understand this, no land processor means no output file
359covering this part of globe; its only when files are stitched together into one that you
360can leave a hole}
361
362%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
363\begin{figure}[!ht] \label{Fig_mppini2}  \begin{center}
364\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_mppini2.pdf}
365\caption {Example of Atlantic domain defined for the CLIPPER projet. Initial grid is
366composed of 773 x 1236 horizontal points. (a) the domain is split onto 9 \time 20
367subdomains (jpni=9, jpnj=20). 52 subdomains are land areas. (b) 52 subdomains
368are eliminated (white rectangles) and the resulting number of processors really
369used during the computation is jpnij=128.}
370\end{center}   \end{figure}
371%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
372
373
374% ================================================================
375% Open Boundary Conditions
376% ================================================================
377\section{Open Boundary Conditions (\key{obc})}
378\label{LBC_obc}
379%-----------------------------------------nam_obc  -------------------------------------------
380%-    nobc_dta    =    0     !  = 0 the obc data are equal to the initial state
381%-                           !  = 1 the obc data are read in 'obc   .dta' files
382%-    rdpein      =    1.    !  ???
383%-    rdpwin      =    1.    !  ???
384%-    rdpnin      =   30.    !  ???
385%-    rdpsin      =    1.    !  ???
386%-    rdpeob      = 1500.    !  time relaxation (days) for the east  open boundary
387%-    rdpwob      =   15.    !    "        "           for the west  open boundary
388%-    rdpnob      =  150.    !    "        "           for the north open boundary
389%-    rdpsob      =   15.    !    "        "           for the south open boundary
390%-    zbsic1      =  140.e+6 !  barotropic stream function on first  isolated coastline
391%-    zbsic2      =    1.e+6 !    "                   "    on second isolated coastline
392%-    zbsic3      =    0.    !    "                   "    on thrid  isolated coastline
393%-    ln_obc_clim = .true.   !  climatological obc data files (default T)
394%-    ln_vol_cst  = .true.   !  total volume conserved
395\namdisplay{namobc} 
396
397It is often necessary to implement a model configuration limited to an oceanic
398region or a basin, which communicates with the rest of the global ocean through
399''open boundaries''. As stated by \citet{Roed1986}, an open boundary is a
400computational border where the aim of the calculations is to allow the perturbations
401generated inside the computational domain to leave it without deterioration of the
402inner model solution. However, an open boundary also has to let information from
403the outer ocean enter the model and should support inflow and outflow conditions.
404
405The open boundary package OBC is the first open boundary option developed in
406NEMO (originally in OPA8.2). It allows the user to
407\begin{itemize}
408\item tell the model that a boundary is ''open'' and not closed by a wall, for example
409by modifying the calculation of the divergence of velocity there;
410\item impose values of tracers and velocities at that boundary (values which may
411be taken from a climatology): this is the``fixed OBC'' option.
412\item calculate boundary values by a sophisticated algorithm combining radiation
413and relaxation (``radiative OBC'' option)
414\end{itemize}
415
416The package resides in the OBC directory. It is described here in four parts: the
417boundary geometry (parameters to be set in \mdl{obc\_par}), the forcing data at
418the boundaries (module \mdl{obcdta}),  the radiation algorithm involving the
419namelist and module \mdl{obcrad}, and a brief presentation of boundary update
420and restart files.
421
422%----------------------------------------------
423\subsection{Boundary geometry}
424\label{OBC_geom}
425%
426First one has to realize that open boundaries may not necessarily be located
427at the extremities of the computational domain. They may exist in the middle
428of the domain, for example at Gibraltar Straits if one wants to avoid including
429the Mediterranean in an Atlantic domain. This flexibility has been found necessary
430for the CLIPPER project \citep{Treguier2001}. Because of the complexity of the
431geometry of ocean basins, it may even be necessary to have more than one
432''west'' open boundary, more than one ''north'', etc. This is not possible with
433the OBC option: only one open boundary of each kind, west, east, south and
434north is allowed; these names refer to the grid geometry (not to the direction
435of the geographical ''west'', ''east'', etc).
436
437The open boundary geometry is set by a series of parameters in the module
438\mdl{obc\_par}. For an eastern open boundary, parameters are \jp{lp\_obc\_east} 
439(true if an east open boundary exists), \jp{jpieob} the $i$-index along which
440the eastern open boundary (eob) is located, \jp{jpjed} the $j$-index at which
441it starts, and \jp{jpjef} the $j$-index where it ends (note $d$ is for ''d\'{e}but''
442and $f$ for ''fin'' in French). Similar parameters exist for the west, south and
443north cases (Table~\ref{Tab_obc_param}).
444
445
446%--------------------------------------------------TABLE--------------------------------------------------
447
448\begin{table}[htbp]  \label{Tab_obc_param}
449\begin{center}
450\begin{tabular}{|l|c|c|c|}
451\hline
452Boundary and  & Constant index  & Starting index (d\'{e}but) & Ending index (fin) \\
453Logical flag  &                 &                            &                     \\
454\hline
455West          & \jp{jpiwob} $>= 2$         &  \jp{jpjwd}$>= 2$          &  \jp{jpjwf}<= \jp{jpjglo}-1 \\
456lp\_obc\_west & $i$-index of a $u$ point   & $j$ of a $T$ point   &$j$ of a $T$ point \\
457\hline
458East            & \jp{jpieob}$<=$\jp{jpiglo}-2&\jp{jpjed} $>= 2$         & \jp{jpjef}$<=$ \jp{jpjglo}-1 \\
459 lp\_obc\_east  & $i$-index of a $u$ point    & $j$ of a $T$ point & $j$ of a $T$ point \\
460\hline
461South           & \jp{jpjsob} $>= 2$         & \jp{jpisd} $>= 2$          & \jp{jpisf}$<=$\jp{jpiglo}-1 \\
462lp\_obc\_south  & $j$-index of a $v$ point   & $i$ of a $T$ point   & $i$ of a $T$ point \\
463\hline
464North           & \jp{jpjnob} $<=$ \jp{jpjglo}-2& \jp{jpind} $>= 2$        & \jp{jpinf}$<=$\jp{jpiglo}-1 \\
465lp\_obc\_north  & $j$-index of a $v$ point      & $i$  of a $T$ point & $i$ of a $T$ point \\
466\hline
467\end{tabular}
468\end{center}
469\caption{Names of different indices relating to the open boundaries. In the case
470of a completely open ocean domain with four ocean boundaries, the parameters
471take exactly the values indicated.}
472\end{table}
473
474The open boundaries must be along coordinate lines. On the C-grid, the boundary
475itself is along a line of normal velocity points: $v$ points for a zonal open boundary
476(the south or north one), and $u$ points for a meridional open boundary (the west
477or east one). Another constraint is that there still must be a row of masked points
478all around the domain, as if the domain were a closed basin (unless periodic conditions
479are used together with open boundary conditions). Therefore, an open boundary
480cannot be located at the first/last index, namely, 1, \jp{jpiglo} or \jp{jpjglo}. Also,
481the open boundary algorithm involves calculating the normal velocity points situated
482just on the boundary, as well as the tangential velocity and temperature and salinity
483just outside the boundary. This means that for a west/south boundary, normal
484velocities and temperature are calculated at the same index \jp{jpiwob} and
485\jp{jpjsob}, respectively. For an east/north boundary, the normal velocity is
486calculated at index \jp{jpieob} and \jp{jpjnob}, but the ``outside'' temperature is
487at index \jp{jpieob}+1 and \jp{jpjnob}+1. This means that \jp{jpieob}, \jp{jpjnob} 
488cannot be bigger than \jp{jpiglo}-2, \jp{jpjglo}-2.
489
490
491The starting and ending indices are to be thought of as $T$ point indices: in many
492cases they indicate the first land $T$-point, at the extremity of an open boundary
493(the coast line follows the $f$ grid points, see Fig.~\ref{Fig_obc_north} for an example
494of a northern open boundary). All indices are relative to the global domain. In the
495free surface case it is possible to have ``ocean corners'', that is, an open boundary
496starting and ending in the ocean.
497
498%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
499\begin{figure}[!t] \label{Fig_obc_north}  \begin{center}
500\includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_obc_north.pdf}
501\caption {Localization of the North open boundary points.}
502\end{center} 
503\end{figure}
504%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
505
506Although not compulsory, it is highly recommended that the bathymetry in the
507vicinity of an open boundary follows the following rule: in the direction perpendicular
508to the open line, the water depth should be constant for 4 grid points. This is in
509order to ensure that the radiation condition, which involves model variables next
510to the boundary, is calculated in a consistent way. On Fig.\ref{Fig_obc_north} we
511indicate by an $=$ symbol, the points which should have the same depth. It means
512that at the 4 points near the boundary, the bathymetry is cylindrical \gmcomment{not sure
513why cylindrical}. The line behind the open $T$-line must be 0 in the bathymetry file
514(as shown on Fig.\ref{Fig_obc_north} for example).
515
516%----------------------------------------------
517\subsection{Boundary data}
518\label{OBC_data}
519
520It is necessary to provide information at the boundaries. The simplest case is
521when this information does not change in time and is equal to the initial conditions
522(namelist variable \np{nobc\_dta}=0). This is the case for the standard configuration
523EEL5 with open boundaries. When (\np{nobc\_dta}=1), open boundary information
524is read from netcdf files. For convenience the input files are supposed to be similar
525to the ''history'' NEMO output files, for dimension names and variable names.
526Open boundary arrays must be dimensioned according to the parameters of table~
527\ref{Tab_obc_param}: for example, at the western boundary, arrays have a
528dimension of \jp{jpwf}-\jp{jpwd}+1 in the horizontal and \jp{jpk} in the vertical.
529
530When ocean observations are used to generate the boundary data (a hydrographic
531section for example, as in \citet{Treguier2001}) it happens often that only the velocity
532normal to the boundary is known, which is the reason why the initial OBC code
533assumes that only $T$, $S$, and the normal velocity ($u$ or $v$) needs to be
534specified. As more and more global model solutions and ocean analysis products
535become available, it will be possible to provide information about all the variables
536(including the tangential velocity) so that the specification of four variables at each
537boundaries will become standard. For the sea surface height, one must distinguish
538between the filtered free surface case and the time-splitting or explicit treatment of
539the free surface.
540 In the first case, it is assumed that the user does not wish to represent high
541 frequency motions such as tides. The boundary condition is thus one of zero
542 normal gradient of sea surface height at the open boundaries, following \citet{Marchesiello2001}.
543No information other than the total velocity needs to be provided at the open
544boundaries in that case. In the other two cases (time splitting or explicit free surface),
545the user must provide barotropic information (sea surface height and barotropic
546velocities) and the use of the Flather algorithm for barotropic variables is
547recommanded. However, this algorithm has not yet been fully tested and bugs
548remain in NEMO v2.3. Users should read the code carefully before using it. Finally,
549in the case of the rigid lid approximation the barotropic streamfunction must be
550provided, as documented in \citet{Treguier2001}). This option is no longer
551recommended but remains in NEMO V2.3.
552
553One frequently encountered case is when an open boundary domain is constructed
554from a global or larger scale NEMO configuration. Assuming the domain corresponds
555to indices $ib:ie$, $jb:je$ of the global domain, the bathymetry and forcing of the
556small domain can be created by using the following netcdf utility on the global files:
557ncks -F $-d\;x,ib,ie$ $-d\;y,jb,je$ (part of the nco series of utilities, see http://nco.sourceforge.net).
558The open boundary files can be constructed using ncks
559commands, following table~\ref{Tab_obc_ind}.
560
561%--------------------------------------------------TABLE--------------------------------------------------
562
563\begin{table}[htbp]  \label{Tab_obc_ind}
564\begin{center}
565\begin{tabular}{|l|c|c|c|c|c|}
566\hline
567OBC  & Variable   & file name      & Index  & Start  & end  \\
568West &  T,S       &   obcwest\_TS.nc &  $ib$+1     &   $jb$+1 &  $je-1$  \\
569     &    U       &   obcwest\_U.nc  &  $ib$+1     &   $jb$+1 &  $je-1$  \\ 
570     &    V       &   obcwest\_V.nc  &  $ib$+1     &   $jb$+1 &  $je-1$  \\       
571\hline
572East &  T,S       &   obceast\_TS.nc &  $ie$-1     &   $jb$+1 &  $je-1$  \\
573     &    U       &   obceast\_U.nc  &  $ie$-2     &   $jb$+1 &  $je-1$  \\ 
574     &    V       &   obceast\_V.nc  &  $ie$-1     &   $jb$+1 &  $je-1$  \\       
575\hline         
576South &  T,S      &   obcsouth\_TS.nc &  $jb$+1     &  $ib$+1 &  $ie-1$  \\
577      &    U      &   obcsouth\_U.nc  &  $jb$+1     &  $ib$+1 &  $ie-1$  \\ 
578      &    V      &   obcsouth\_V.nc  &  $jb$+1     &  $ib$+1 &  $ie-1$  \\   
579\hline
580North &  T,S      &   obcnorth\_TS.nc &  $je$-1     &  $ib$+1 &  $ie-1$  \\
581      &    U      &   obcnorth\_U.nc  &  $je$-1     &  $ib$+1 &  $ie-1$  \\ 
582      &    V      &   obcnorth\_V.nc  &  $je$-2     &  $ib$+1 &  $ie-1$  \\ 
583\hline
584\end{tabular}
585\end{center}
586\caption{Requirements for creating open boundary files from a global configuration,
587appropriate for the subdomain of indices $ib:ie$, $jb:je$. ``Index'' designates the
588$i$ or $j$ index along which the $u$ of $v$ boundary point is situated in the global
589configuration, starting and ending with the $j$ or $i$ indices indicated.
590For example, to generate file obcnorth\_V.nc, use the command ncks
591$-F$ $-d\;y,je-2$  $-d\;x,ib+1,ie-1$ } 
592\end{table}
593
594It is assumed that the open boundary files contain the variables for the period of
595the model integration. If the boundary files contain one time frame, the boundary
596data is held fixed in time. If the files contain 12 values, it is assumed that the input
597is a climatology for a repeated annual cycle (corresponding to the case \np{ln\_obc\_clim} 
598= .True.). The case of an arbitrary number of time frames is not yet implemented
599correctly; the user is required to write his own code in the module \mdl{obc\_dta} 
600to deal with this situation.
601
602\subsection{Radiation algorithm}
603\label{OBC_rad}
604
605The art of open boundary management consists in applying a constraint strong
606enough that the inner domain "feels" the rest of the ocean, but weak enough
607that perturbations are allowed to leave the domain with minimum false reflections
608of energy. The constraints are specified separately at each boundary as time
609scales for ''inflow'' and ''outflow'' as defined below. The time scales are set (in days)
610by namelist parameters such as \np{rdpein}, \np{rdpeob} for the eastern open
611boundary for example. When both time scales are zero for a given boundary
612($e.g.$ for the western boundary, \jp{lp\_obc\_west}=.True., \np{rdpwob}=0 and
613\np{rdpwin}=0) this means that the boundary in question is a ''fixed '' boundary
614where the solution is set exactly by the boundary data. This is not recommended,
615except in combination with increased viscosity in a ''sponge'' layer next to the
616boundary in order to avoid spurious reflections. 
617
618
619The radiation\/relaxation \gmcomment{the / doesnt seem to appear in the output} 
620algorithm is applied when either relaxation time (for ''inflow'' or ''outflow'') is
621non-zero. It has been developed and tested in the SPEM model and its
622successor ROMS \citep{Barnier1996, Marchesiello2001}, which is an
623$s$-coordinate model on an Arakawa C-grid. Although the algorithm has
624been numerically successful in the CLIPPER Atlantic models, the physics
625do not work as expected \citep{Treguier2001}. Users are invited to consider
626open boundary conditions (OBC hereafter) with some scepticism
627\citep{Durran2001, Blayo2005}.
628
629The first part of the algorithm calculates a phase velocity to determine
630whether perturbations tend to propagate toward, or away from, the
631boundary. Let us consider a model variable $\phi$.
632The phase velocities ($C_{\phi x}$,$C_{\phi y}$) for the variable $\phi$,
633in the directions normal and tangential to the boundary are
634\begin{equation} \label{Eq_obc_cphi}
635C_{\phi x} = \frac{ -\phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{x} 
636\;\;\;\;\; \;\;\; 
637C_{\phi y} = \frac{ -\phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{y}.
638\end{equation}
639Following \citet{Treguier2001} and \citet{Marchesiello2001} we retain only
640the normal component of the velocity, $C_{\phi x}$, setting $C_{\phi y} =0$ 
641(but unlike the original Orlanski radiation algorithm we retain $\phi_{y}$ in
642the expression for $C_{\phi x}$). 
643
644The discrete form of (\ref{Eq_obc_cphi}), described by \citet{Barnier1998},
645takes into account the two rows of grid points situated inside the domain
646next to the boundary, and the three previous time steps ($n$, $n-1$,
647and $n-2$). The same equation can then be discretized at the boundary at
648time steps $n-1$, $n$ and $n+1$ \gmcomment{since the original was three time-level} 
649in order to extrapolate for the new boundary value $\phi^{n+1}$.
650
651In the open boundary algorithm as implemented in NEMO v2.3, the new boundary
652values are updated differently depending on the sign of $C_{\phi x}$. Let us take
653an eastern boundary as an example. The solution for variable $\phi$ at the
654boundary is given by a generalized wave equation with phase velocity $C_{\phi}$,
655with the addition of a relaxation term, as:
656\begin{eqnarray}
657\phi_{t} &  =  & -C_{\phi x} \phi_{x} + \frac{1}{\tau_{o}} (\phi_{c}-\phi)
658                        \;\;\; \;\;\; \;\;\; (C_{\phi x} > 0), \label{Eq_obc_rado} \\
659\phi_{t} &  =  & \frac{1}{\tau_{i}} (\phi_{c}-\phi)
660\;\;\; \;\;\; \;\;\;\;\;\; (C_{\phi x} < 0), \label{Eq_obc_radi}
661\end{eqnarray}
662where $\phi_{c}$ is the estimate of $\phi$ at the boundary, provided as boundary
663data. Note that in (\ref{Eq_obc_rado}), $C_{\phi x}$ is bounded by the ratio
664$\delta x/\delta t$ for stability reasons. When $C_{\phi x}$ is eastward (outward
665propagation), the radiation condition (\ref{Eq_obc_rado}) is used.
666When  $C_{\phi x}$ is westward (inward propagation), (\ref{Eq_obc_radi}) is
667used with a strong relaxation to climatology (usually $\tau_{i}=\np{rdpein}=$1~day).
668Equation (\ref{Eq_obc_radi}) is solved with a Euler time-stepping scheme. As a
669consequence, setting $\tau_{i}$ smaller than, or equal to the time step is equivalent
670to a fixed boundary condition. A time scale of one day is usually a good compromise
671which guarantees that the inflow conditions remain close to climatology while ensuring
672numerical stability.
673
674In  the case of a western boundary located in the Eastern Atlantic, \citet{Penduff2000} 
675have been able to implement the radiation algorithm without any boundary data,
676using persistence from the previous time step instead. This solution has not worked
677in other cases \citep{Treguier2001}, so that the use of boundary data is recommended.
678Even in the outflow condition (\ref{Eq_obc_rado}), we have found it desirable to
679maintain a weak relaxation to climatology. The time step is usually chosen so as to
680be larger than typical turbulent scales (of order 1000~days \gmcomment{or maybe seconds?}).
681
682The radiation condition is applied to the model variables: temperature, salinity,
683tangential and normal velocities. For normal and tangential velocities, $u$ and $v$,
684radiation is applied with phase velocities calculated from $u$ and $v$ respectively. 
685For the radiation of tracers, we use the phase velocity calculated from the tangential
686velocity in order to avoid calculating too many independent radiation velocities and
687because tangential velocities and tracers have the same position along the boundary
688on a C-grid. 
689
690\subsection{Domain decomposition (\key{mpp\_mpi})}
691\label{OBC_mpp}
692When \key{mpp\_mpi} is active in the code, the computational domain is divided
693into rectangles that are attributed each to a different processor. The open boundary
694code is ``mpp-compatible'' up to a certain point. The radiation algorithm will not
695work if there is an mpp subdomain boundary parallel to the open boundary at the
696index of the boundary, or the grid point after (outside), or three grid points before
697(inside). On the other hand, there is no problem if an mpp subdomain boundary
698cuts the open boundary perpendicularly. These geometrical limitations must be
699checked for by the user (there is no safeguard in the code). 
700The general principle for the open boundary mpp code is that loops over the open
701boundaries {not sure what this means} are performed on local indices (nie0,
702nie1, nje0, nje1 for an eastern boundary for instance) that are initialized in module
703\mdl{obc\_ini}. Those indices have relevant values on the processors that contain
704a segment of an open boundary. For processors that do not include an open
705boundary segment, the indices are such that the calculations within the loops are
706not performed.
707\gmcomment{I dont understand most of the last few sentences}
708 
709Arrays of climatological data that are read from files are seen by all processors
710and have the same dimensions for all (for instance, for the eastern boundary,
711uedta(jpjglo,jpk,2)). On the other hand, the arrays for the calculation of radiation
712are local to each processor (uebnd(jpj,jpk,3,3) for instance).  This allowed the
713CLIPPER model for example, to save on memory where the eastern boundary
714crossed 8 processors so that \jp{jpj} was much smaller than (\jp{jpjef}-\jp{jpjed}+1).
715
716\subsection{Volume conservation}
717\label{OBC_vol}
718
719It is necessary to control the volume inside a domain when using open boundaries.
720With fixed boundaries, it is enough to ensure that the total inflow/outflow has
721reasonable values (either zero or a value compatible with an observed volume
722balance). When using radiative boundary conditions it is necessary to have a
723volume constraint because each open boundary works independently from the
724others. The methodology used to control this volume is identical to the one
725coded in the ROMS model \citep{Marchesiello2001}.
726
727
728%---------------------------------------- EXTRAS
729\colorbox{yellow}{Explain obc\_vol{\ldots}}
730
731\colorbox{yellow}{OBC algorithm for update, OBC restart, list of routines where obc key appears{\ldots}}
732
733\colorbox{yellow}{OBC rigid lid? {\ldots}}
734
735
736
737
738% ====================================================================
739% Flow Relaxation Scheme
740% ====================================================================
741\section{Flow Relaxation Scheme (???)}
742\label{LBC_bdy}
743
744%gm% to be updated by Met Office
Note: See TracBrowser for help on using the repository browser.