source: branches/UKMO/v3_6_extra_CMIP6_diagnostics/DOC/TexFiles/Chapters/Chap_LBC.tex @ 7492

Last change on this file since 7492 was 7492, checked in by timgraham, 5 years ago

Final merge of head of v3.6 STABLE branch into this branch

  • Property svn:keywords set to Id
File size: 35.7 KB
Line 
1\documentclass[NEMO_book]{subfiles}
2\begin{document}
3% ================================================================
4% Chapter — Lateral Boundary Condition (LBC)
5% ================================================================
6\chapter{Lateral Boundary Condition (LBC) }
7\label{LBC}
8\minitoc
9
10\newpage
11$\ $\newline    % force a new ligne
12
13
14%gm% add here introduction to this chapter
15
16% ================================================================
17% Boundary Condition at the Coast
18% ================================================================
19\section{Boundary Condition at the Coast (\np{rn\_shlat})}
20\label{LBC_coast}
21%--------------------------------------------nam_lbc-------------------------------------------------------
22\namdisplay{namlbc} 
23%--------------------------------------------------------------------------------------------------------------
24
25%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}).
26
27%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}.
28
29Options are defined through the \ngn{namlbc} namelist variables.
30The discrete representation of a domain with complex boundaries (coastlines and
31bottom topography) leads to arrays that include large portions where a computation
32is not required as the model variables remain at zero. Nevertheless, vectorial
33supercomputers are far more efficient when computing over a whole array, and the
34readability of a code is greatly improved when boundary conditions are applied in
35an automatic way rather than by a specific computation before or after each
36computational loop. An efficient way to work over the whole domain while specifying
37the boundary conditions, is to use multiplication by mask arrays in the computation.
38A mask array is a matrix whose elements are $1$ in the ocean domain and $0$ 
39elsewhere. A simple multiplication of a variable by its own mask ensures that it will
40remain zero over land areas. Since most of the boundary conditions consist of a
41zero flux across the solid boundaries, they can be simply applied by multiplying
42variables by the correct mask arrays, $i.e.$ the mask array of the grid point where
43the flux is evaluated. For example, the heat flux in the \textbf{i}-direction is evaluated
44at $u$-points. Evaluating this quantity as,
45
46\begin{equation} \label{Eq_lbc_aaaa}
47\frac{A^{lT} }{e_1 }\frac{\partial T}{\partial i}\equiv \frac{A_u^{lT} 
48}{e_{1u} } \; \delta _{i+1 / 2} \left[ T \right]\;\;mask_u
49\end{equation}
50(where mask$_{u}$ is the mask array at a $u$-point) ensures that the heat flux is
51zero inside land and at the boundaries, since mask$_{u}$ is zero at solid boundaries
52which in this case are defined at $u$-points (normal velocity $u$ remains zero at
53the coast) (Fig.~\ref{Fig_LBC_uv}).
54
55%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
56\begin{figure}[!t]     \begin{center}
57\includegraphics[width=0.90\textwidth]{Fig_LBC_uv}
58\caption{  \label{Fig_LBC_uv}
59Lateral boundary (thick line) at T-level. The velocity normal to the boundary is set to zero.}
60\end{center}   \end{figure}
61%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
62
63For momentum the situation is a bit more complex as two boundary conditions
64must be provided along the coast (one each for the normal and tangential velocities).
65The boundary of the ocean in the C-grid is defined by the velocity-faces.
66For example, at a given $T$-level, the lateral boundary (a coastline or an intersection
67with the bottom topography) is made of segments joining $f$-points, and normal
68velocity points are located between two $f-$points (Fig.~\ref{Fig_LBC_uv}).
69The boundary condition on the normal velocity (no flux through solid boundaries)
70can thus be easily implemented using the mask system. The boundary condition
71on the tangential velocity requires a more specific treatment. This boundary
72condition influences the relative vorticity and momentum diffusive trends, and is
73required in order to compute the vorticity at the coast. Four different types of
74lateral boundary condition are available, controlled by the value of the \np{rn\_shlat} 
75namelist parameter. (The value of the mask$_{f}$ array along the coastline is set
76equal to this parameter.) These are:
77
78%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
79\begin{figure}[!p] \begin{center}
80\includegraphics[width=0.90\textwidth]{Fig_LBC_shlat}
81\caption{     \label{Fig_LBC_shlat} 
82lateral boundary condition (a) free-slip ($rn\_shlat=0$) ; (b) no-slip ($rn\_shlat=2$)
83; (c) "partial" free-slip ($0<rn\_shlat<2$) and (d) "strong" no-slip ($2<rn\_shlat$).
84Implied "ghost" velocity inside land area is display in grey. }
85\end{center}    \end{figure}
86%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
87
88\begin{description}
89
90\item[free-slip boundary condition (\np{rn\_shlat}=0): ]  the tangential velocity at the
91coastline is equal to the offshore velocity, $i.e.$ the normal derivative of the
92tangential velocity is zero at the coast, so the vorticity: mask$_{f}$ array is set
93to zero inside the land and just at the coast (Fig.~\ref{Fig_LBC_shlat}-a).
94
95\item[no-slip boundary condition (\np{rn\_shlat}=2): ] the tangential velocity vanishes
96at the coastline. Assuming that the tangential velocity decreases linearly from
97the closest ocean velocity grid point to the coastline, the normal derivative is
98evaluated as if the velocities at the closest land velocity gridpoint and the closest
99ocean velocity gridpoint were of the same magnitude but in the opposite direction
100(Fig.~\ref{Fig_LBC_shlat}-b). Therefore, the vorticity along the coastlines is given by:
101
102\begin{equation*}
103\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) \ ,
104\end{equation*}
105where $u$ and $v$ are masked fields. Setting the mask$_{f}$ array to $2$ along
106the coastline provides a vorticity field computed with the no-slip boundary condition,
107simply by multiplying it by the mask$_{f}$ :
108\begin{equation} \label{Eq_lbc_bbbb}
109\zeta \equiv \frac{1}{e_{1f} {\kern 1pt}e_{2f} }\left( {\delta _{i+1/2} 
110\left[ {e_{2v} \,v} \right]-\delta _{j+1/2} \left[ {e_{1u} \,u} \right]} 
111\right)\;\mbox{mask}_f
112\end{equation}
113
114\item["partial" free-slip boundary condition (0$<$\np{rn\_shlat}$<$2): ] the tangential
115velocity at the coastline is smaller than the offshore velocity, $i.e.$ there is a lateral
116friction but not strong enough to make the tangential velocity at the coast vanish
117(Fig.~\ref{Fig_LBC_shlat}-c). This can be selected by providing a value of mask$_{f}$ 
118strictly inbetween $0$ and $2$.
119
120\item["strong" no-slip boundary condition (2$<$\np{rn\_shlat}): ] the viscous boundary
121layer is assumed to be smaller than half the grid size (Fig.~\ref{Fig_LBC_shlat}-d).
122The friction is thus larger than in the no-slip case.
123
124\end{description}
125
126Note that when the bottom topography is entirely represented by the $s$-coor-dinates
127(pure $s$-coordinate), the lateral boundary condition on tangential velocity is of much
128less importance as it is only applied next to the coast where the minimum water depth
129can be quite shallow.
130
131The alternative numerical implementation of the no-slip boundary conditions for an
132arbitrary coast line of \citet{Shchepetkin1996} is also available through the
133\key{noslip\_accurate} CPP key. It is based on a fourth order evaluation of the shear at the
134coast which, in turn, allows a true second order scheme in the interior of the domain
135($i.e.$ the numerical boundary scheme simulates the truncation error of the numerical
136scheme used in the interior of the domain). \citet{Shchepetkin1996} found that such a
137technique considerably improves the quality of the numerical solution. In \NEMO, such
138spectacular improvements have not been found in the half-degree global ocean
139(ORCA05), but significant reductions of numerically induced coastal upwellings were
140found in an eddy resolving simulation of the Alboran Sea \citep{Olivier_PhD01}.
141Nevertheless, since a no-slip boundary condition is not recommended in an eddy
142permitting or resolving simulation \citep{Penduff_al_OS07}, the use of this option is also
143not recommended.
144
145In practice, the no-slip accurate option changes the way the curl is evaluated at the
146coast (see \mdl{divcur} module), and requires the nature of each coastline grid point
147(convex or concave corners, straight north-south or east-west coast) to be specified. 
148This is performed in routine \rou{dom\_msk\_nsa} in the \mdl{domask} module.
149
150% ================================================================
151% Boundary Condition around the Model Domain
152% ================================================================
153\section{Model Domain Boundary Condition (\np{jperio})}
154\label{LBC_jperio}
155
156At the model domain boundaries several choices are offered: closed, cyclic east-west,
157south symmetric across the equator, a north-fold, and combination closed-north fold
158or cyclic-north-fold. The north-fold boundary condition is associated with the 3-pole ORCA mesh.
159
160% -------------------------------------------------------------------------------------------------------------
161%        Closed, cyclic, south symmetric (\np{jperio} = 0, 1 or 2)
162% -------------------------------------------------------------------------------------------------------------
163\subsection{Closed, cyclic, south symmetric (\np{jperio} = 0, 1 or 2)}
164\label{LBC_jperio012}
165
166The choice of closed, cyclic or symmetric model domain boundary condition is made
167by setting \np{jperio} to 0, 1 or 2 in namelist \ngn{namcfg}. Each time such a boundary
168condition is needed, it is set by a call to routine \mdl{lbclnk}. The computation of
169momentum and tracer trends proceeds from $i=2$ to $i=jpi-1$ and from $j=2$ to
170$j=jpj-1$, $i.e.$ in the model interior. To choose a lateral model boundary condition
171is to specify the first and last rows and columns of the model variables.
172
173\begin{description}
174
175\item[For closed boundary (\textit{jperio=0})], solid walls are imposed at all model
176boundaries: first and last rows and columns are set to zero.
177
178\item[For cyclic east-west boundary (\textit{jperio=1})], first and last rows are set
179to zero (closed) whilst the first column is set to the value of the last-but-one column
180and the last column to the value of the second one (Fig.~\ref{Fig_LBC_jperio}-a).
181Whatever flows out of the eastern (western) end of the basin enters the western
182(eastern) end. Note that there is no option for north-south cyclic or for doubly
183cyclic cases.
184
185\item[For symmetric boundary condition across the equator (\textit{jperio=2})],
186last rows, and first and last columns are set to zero (closed). The row of symmetry
187is chosen to be the $u$- and $T-$points equator line ($j=2$, i.e. at the southern
188end of the domain). For arrays defined at $u-$ or $T-$points, the first row is set
189to the value of the third row while for most of $v$- and $f$-point arrays ($v$, $\zeta$,
190$j\psi$, but \gmcomment{not sure why this is "but"} scalar arrays such as eddy coefficients)
191the first row is set to minus the value of the second row (Fig.~\ref{Fig_LBC_jperio}-b).
192Note that this boundary condition is not yet available for the case of a massively
193parallel computer (\textbf{key{\_}mpp} defined).
194
195\end{description}
196
197%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
198\begin{figure}[!t]     \begin{center}
199\includegraphics[width=1.0\textwidth]{Fig_LBC_jperio}
200\caption{    \label{Fig_LBC_jperio}
201setting of (a) east-west cyclic  (b) symmetric across the equator boundary conditions.}
202\end{center}   \end{figure}
203%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
204
205% -------------------------------------------------------------------------------------------------------------
206%        North fold (\textit{jperio = 3 }to $6)$
207% -------------------------------------------------------------------------------------------------------------
208\subsection{North-fold (\textit{jperio = 3 }to $6$) }
209\label{LBC_north_fold}
210
211The north fold boundary condition has been introduced in order to handle the north
212boundary of a three-polar ORCA grid. Such a grid has two poles in the northern hemisphere
213(Fig.\ref{Fig_MISC_ORCA_msh}, and thus requires a specific treatment illustrated in Fig.\ref{Fig_North_Fold_T}.
214Further information can be found in \mdl{lbcnfd} module which applies the north fold boundary condition.
215
216%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
217\begin{figure}[!t]    \begin{center}
218\includegraphics[width=0.90\textwidth]{Fig_North_Fold_T}
219\caption{    \label{Fig_North_Fold_T} 
220North fold boundary with a $T$-point pivot and cyclic east-west boundary condition
221($jperio=4$), as used in ORCA 2, 1/4, and 1/12. Pink shaded area corresponds
222to the inner domain mask (see text). }
223\end{center}   \end{figure}
224%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
225
226% ====================================================================
227% Exchange with neighbouring processors
228% ====================================================================
229\section  [Exchange with neighbouring processors (\textit{lbclnk}, \textit{lib\_mpp})]
230      {Exchange with neighbouring processors (\mdl{lbclnk}, \mdl{lib\_mpp})}
231\label{LBC_mpp}
232
233For massively parallel processing (mpp), a domain decomposition method is used.
234The basic idea of the method is to split the large computation domain of a numerical
235experiment into several smaller domains and solve the set of equations by addressing
236independent local problems. Each processor has its own local memory and computes
237the model equation over a subdomain of the whole model domain. The subdomain
238boundary conditions are specified through communications between processors
239which are organized by explicit statements (message passing method).
240
241A big advantage is that the method does not need many modifications of the initial
242FORTRAN code. From the modeller's point of view, each sub domain running on
243a processor is identical to the "mono-domain" code. In addition, the programmer
244manages the communications between subdomains, and the code is faster when
245the number of processors is increased. The porting of OPA code on an iPSC860
246was achieved during Guyon's PhD [Guyon et al. 1994, 1995] in collaboration with
247CETIIS and ONERA. The implementation in the operational context and the studies
248of performance on a T3D and T3E Cray computers have been made in collaboration
249with IDRIS and CNRS. The present implementation is largely inspired by Guyon's
250work  [Guyon 1995].
251
252The parallelization strategy is defined by the physical characteristics of the
253ocean model. Second order finite difference schemes lead to local discrete
254operators that depend at the very most on one neighbouring point. The only
255non-local computations concern the vertical physics (implicit diffusion,
256turbulent closure scheme, ...) (delocalization over the whole water column),
257and the solving of the elliptic equation associated with the surface pressure
258gradient computation (delocalization over the whole horizontal domain).
259Therefore, a pencil strategy is used for the data sub-structuration
260: the 3D initial domain is laid out on local processor
261memories following a 2D horizontal topological splitting. Each sub-domain
262computes its own surface and bottom boundary conditions and has a side
263wall overlapping interface which defines the lateral boundary conditions for
264computations in the inner sub-domain. The overlapping area consists of the
265two rows at each edge of the sub-domain. After a computation, a communication
266phase starts: each processor sends to its neighbouring processors the update
267values of the points corresponding to the interior overlapping area to its
268neighbouring sub-domain ($i.e.$ the innermost of the two overlapping rows).
269The communication is done through the Message Passing Interface (MPI).
270The data exchanges between processors are required at the very
271place where lateral domain boundary conditions are set in the mono-domain
272computation : the \rou{lbc\_lnk} routine (found in \mdl{lbclnk} module)
273which manages such conditions is interfaced with routines found in \mdl{lib\_mpp} module
274when running on an MPP computer ($i.e.$ when \key{mpp\_mpi} defined).
275It has to be pointed out that when using the MPP version of the model,
276the east-west cyclic boundary condition is done implicitly,
277whilst the south-symmetric boundary condition option is not available.
278
279%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
280\begin{figure}[!t]    \begin{center}
281\includegraphics[width=0.90\textwidth]{Fig_mpp}
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 \NEMO, the splitting is regular and arithmetic.
288The 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$ (parameters set in
290 \ngn{nammpp} namelist). Each processor is independent and without message passing
291 or synchronous process, programs run alone and access just its own local memory.
292 For this reason, the main model dimensions are now the local dimensions of the subdomain (pencil)
293 that are named \jp{jpi}, \jp{jpj}, \jp{jpk}. These dimensions include the internal
294 domain and the overlapping rows. The number of rows to exchange (known as
295 the halo) is usually set to one (\jp{jpreci}=1, in \mdl{par\_oce}). The whole domain
296 dimensions are named \np{jpiglo}, \np{jpjglo} and \jp{jpk}. The relationship between
297 the whole domain and a sub-domain is:
298\begin{eqnarray} 
299      jpi & = & ( jpiglo-2*jpreci + (jpni-1) ) / jpni + 2*jpreci  \nonumber \\
300      jpj & = & ( jpjglo-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj  \label{Eq_lbc_jpi}
301\end{eqnarray}
302where \jp{jpni}, \jp{jpnj} are the number of processors following the i- and j-axis.
303
304One also defines variables nldi and nlei which correspond to the internal domain bounds,
305and the variables nimpp and njmpp which are the position of the (1,1) grid-point in the global domain.
306An element of $T_{l}$, a local array (subdomain) corresponds to an element of $T_{g}$,
307a global array (whole domain) by the relationship:
308\begin{equation} \label{Eq_lbc_nimpp}
309T_{g} (i+nimpp-1,j+njmpp-1,k) = T_{l} (i,j,k),
310\end{equation}
311with  $1 \leq i \leq jpi$, $1  \leq j \leq jpj $ , and  $1  \leq k \leq jpk$.
312
313Processors are numbered from 0 to $jpnij-1$, the number is saved in the variable
314nproc. In the standard version, a processor has no more than four neighbouring
315processors named nono (for north), noea (east), noso (south) and nowe (west)
316and two variables, nbondi and nbondj, indicate the relative position of the processor :
317\begin{itemize}
318\item       nbondi = -1    an east neighbour, no west processor,
319\item       nbondi =  0 an east neighbour, a west neighbour,
320\item       nbondi =  1    no east processor, a west neighbour,
321\item       nbondi =  2    no splitting following the i-axis.
322\end{itemize}
323During the simulation, processors exchange data with their neighbours.
324If there is effectively a neighbour, the processor receives variables from this
325processor on its overlapping row, and sends the data issued from internal
326domain corresponding to the overlapping row of the other processor.
327
328
329The \NEMO model computes equation terms with the help of mask arrays (0 on land
330points and 1 on sea points). It is easily readable and very efficient in the context of
331a computer with vectorial architecture. However, in the case of a scalar processor,
332computations over the land regions become more expensive in terms of CPU time.
333It is worse when we use a complex configuration with a realistic bathymetry like the
334global ocean where more than 50 \% of points are land points. For this reason, a
335pre-processing tool can be used to choose the mpp domain decomposition with a
336maximum number of only land points processors, which can then be eliminated (Fig. \ref{Fig_mppini2})
337(For example, the mpp\_optimiz tools, available from the DRAKKAR web site).
338This optimisation is dependent on the specific bathymetry employed. The user
339then chooses optimal parameters \jp{jpni}, \jp{jpnj} and \jp{jpnij} with
340$jpnij < jpni \times jpnj$, leading to the elimination of $jpni \times jpnj - jpnij$ 
341land processors. When those parameters are specified in \ngn{nammpp} namelist,
342the algorithm in the \rou{inimpp2} routine sets each processor's parameters (nbound,
343nono, noea,...) so that the land-only processors are not taken into account.
344
345\gmcomment{Note that the inimpp2 routine is general so that the original inimpp
346routine should be suppressed from the code.}
347
348When land processors are eliminated, the value corresponding to these locations in
349the model output files is undefined. Note that this is a problem for the meshmask file
350which requires to be defined over the whole domain. Therefore, user should not eliminate
351land processors when creating a meshmask file ($i.e.$ when setting a non-zero value to \np{nn\_msh}).
352
353%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
354\begin{figure}[!ht]     \begin{center}
355\includegraphics[width=0.90\textwidth]{Fig_mppini2}
356\caption {    \label{Fig_mppini2}
357Example of Atlantic domain defined for the CLIPPER projet. Initial grid is
358composed of 773 x 1236 horizontal points.
359(a) the domain is split onto 9 \time 20 subdomains (jpni=9, jpnj=20).
36052 subdomains are land areas.
361(b) 52 subdomains are eliminated (white rectangles) and the resulting number
362of processors really used during the computation is jpnij=128.}
363\end{center}   \end{figure}
364%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
365
366
367% ====================================================================
368% Unstructured open boundaries BDY
369% ====================================================================
370\section{Unstructured Open Boundary Conditions (\key{bdy}) (BDY)}
371\label{LBC_bdy}
372
373%-----------------------------------------nambdy--------------------------------------------
374\namdisplay{nambdy} 
375%-----------------------------------------------------------------------------------------------
376%-----------------------------------------nambdy_index--------------------------------------------
377\namdisplay{nambdy_index} 
378%-----------------------------------------------------------------------------------------------
379%-----------------------------------------nambdy_dta--------------------------------------------
380\namdisplay{nambdy_dta} 
381%-----------------------------------------------------------------------------------------------
382%-----------------------------------------nambdy_dta--------------------------------------------
383\namdisplay{nambdy_dta2} 
384%-----------------------------------------------------------------------------------------------
385
386Options are defined through the \ngn{nambdy} \ngn{nambdy\_index} 
387\ngn{nambdy\_dta} \ngn{nambdy\_dta2} namelist variables.
388The BDY module is an alternative implementation of open boundary
389conditions for regional configurations. It implements the Flow
390Relaxation Scheme algorithm for temperature, salinity, velocities and
391ice fields, and the Flather radiation condition for the depth-mean
392transports. The specification of the location of the open boundary is
393completely flexible and allows for example the open boundary to follow
394an isobath or other irregular contour.
395
396The BDY module was modelled on the OBC module and shares many features
397and a similar coding structure \citep{Chanut2005}.
398
399The BDY module is completely rewritten at NEMO 3.4 and there is a new
400set of namelists. Boundary data files used with earlier versions of
401NEMO may need to be re-ordered to work with this version. See the
402section on the Input Boundary Data Files for details.
403
404%----------------------------------------------
405\subsection{The namelists}
406\label{BDY_namelist}
407
408It is possible to define more than one boundary ``set'' and apply
409different boundary conditions to each set. The number of boundary
410sets is defined by \np{nb\_bdy}.  Each boundary set may be defined
411as a set of straight line segments in a namelist
412(\np{ln\_coords\_file}=.false.) or read in from a file
413(\np{ln\_coords\_file}=.true.). If the set is defined in a namelist,
414then the namelists nambdy\_index must be included separately, one for
415each set. If the set is defined by a file, then a
416``coordinates.bdy.nc'' file must be provided. The coordinates.bdy file
417is analagous to the usual NEMO ``coordinates.nc'' file. In the example
418above, there are two boundary sets, the first of which is defined via
419a file and the second is defined in a namelist. For more details of
420the definition of the boundary geometry see section
421\ref{BDY_geometry}.
422
423For each boundary set a boundary
424condition has to be chosen for the barotropic solution (``u2d'':
425sea-surface height and barotropic velocities), for the baroclinic
426velocities (``u3d''), and for the active tracers\footnote{The BDY
427  module does not deal with passive tracers at this version}
428(``tra''). For each set of variables there is a choice of algorithm
429and a choice for the data, eg. for the active tracers the algorithm is
430set by \np{nn\_tra} and the choice of data is set by
431\np{nn\_tra\_dta}.
432
433The choice of algorithm is currently as follows:
434
435\mbox{}
436
437\begin{itemize}
438\item[0.] No boundary condition applied. So the solution will ``see''
439  the land points around the edge of the edge of the domain.
440\item[1.] Flow Relaxation Scheme (FRS) available for all variables.
441\item[2.] Flather radiation scheme for the barotropic variables. The
442  Flather scheme is not compatible with the filtered free surface
443  ({\it dynspg\_ts}).
444\end{itemize}
445
446\mbox{}
447
448The main choice for the boundary data is
449to use initial conditions as boundary data (\np{nn\_tra\_dta}=0) or to
450use external data from a file (\np{nn\_tra\_dta}=1). For the
451barotropic solution there is also the option to use tidal
452harmonic forcing either by itself or in addition to other external
453data.
454
455If external boundary data is required then the nambdy\_dta namelist
456must be defined. One nambdy\_dta namelist is required for each boundary
457set in the order in which the boundary sets are defined in nambdy. In
458the example given, two boundary sets have been defined and so there
459are two nambdy\_dta namelists. The boundary data is read in using the
460fldread module, so the nambdy\_dta namelist is in the format required
461for fldread. For each variable required, the filename, the frequency
462of the files and the frequency of the data in the files is given. Also
463whether or not time-interpolation is required and whether the data is
464climatological (time-cyclic) data. Note that on-the-fly spatial
465interpolation of boundary data is not available at this version.
466
467In the example namelists given, two boundary sets are defined. The
468first set is defined via a file and applies FRS conditions to
469temperature and salinity and Flather conditions to the barotropic
470variables. External data is provided in daily files (from a
471large-scale model). Tidal harmonic forcing is also used. The second
472set is defined in a namelist. FRS conditions are applied on
473temperature and salinity and climatological data is read from external
474files.
475
476%----------------------------------------------
477\subsection{The Flow Relaxation Scheme}
478\label{BDY_FRS_scheme}
479
480The Flow Relaxation Scheme (FRS) \citep{Davies_QJRMS76,Engerdahl_Tel95},
481applies a simple relaxation of the model fields to
482externally-specified values over a zone next to the edge of the model
483domain. Given a model prognostic variable $\Phi$ 
484\begin{equation}  \label{Eq_bdy_frs1}
485\Phi(d) = \alpha(d)\Phi_{e}(d) + (1-\alpha(d))\Phi_{m}(d)\;\;\;\;\; d=1,N
486\end{equation}
487where $\Phi_{m}$ is the model solution and $\Phi_{e}$ is the specified
488external field, $d$ gives the discrete distance from the model
489boundary  and $\alpha$ is a parameter that varies from $1$ at $d=1$ to
490a small value at $d=N$. It can be shown that this scheme is equivalent
491to adding a relaxation term to the prognostic equation for $\Phi$ of
492the form:
493\begin{equation}  \label{Eq_bdy_frs2}
494-\frac{1}{\tau}\left(\Phi - \Phi_{e}\right)
495\end{equation}
496where the relaxation time scale $\tau$ is given by a function of
497$\alpha$ and the model time step $\Delta t$:
498\begin{equation}  \label{Eq_bdy_frs3}
499\tau = \frac{1-\alpha}{\alpha}  \,\rdt
500\end{equation}
501Thus the model solution is completely prescribed by the external
502conditions at the edge of the model domain and is relaxed towards the
503external conditions over the rest of the FRS zone. The application of
504a relaxation zone helps to prevent spurious reflection of outgoing
505signals from the model boundary.
506
507The function $\alpha$ is specified as a $tanh$ function:
508\begin{equation}  \label{Eq_bdy_frs4}
509\alpha(d) = 1 - \tanh\left(\frac{d-1}{2}\right),       \quad d=1,N
510\end{equation}
511The width of the FRS zone is specified in the namelist as
512\np{nn\_rimwidth}. This is typically set to a value between 8 and 10.
513
514%----------------------------------------------
515\subsection{The Flather radiation scheme}
516\label{BDY_flather_scheme}
517
518The \citet{Flather_JPO94} scheme is a radiation condition on the normal, depth-mean
519transport across the open boundary. It takes the form
520\begin{equation}  \label{Eq_bdy_fla1}
521U = U_{e} + \frac{c}{h}\left(\eta - \eta_{e}\right),
522\end{equation}
523where $U$ is the depth-mean velocity normal to the boundary and $\eta$
524is the sea surface height, both from the model. The subscript $e$
525indicates the same fields from external sources. The speed of external
526gravity waves is given by $c = \sqrt{gh}$, and $h$ is the depth of the
527water column. The depth-mean normal velocity along the edge of the
528model domain is set equal to the
529external depth-mean normal velocity, plus a correction term that
530allows gravity waves generated internally to exit the model boundary.
531Note that the sea-surface height gradient in \eqref{Eq_bdy_fla1}
532is a spatial gradient across the model boundary, so that $\eta_{e}$ is
533defined on the $T$ points with $nbr=1$ and $\eta$ is defined on the
534$T$ points with $nbr=2$. $U$ and $U_{e}$ are defined on the $U$ or
535$V$ points with $nbr=1$, $i.e.$ between the two $T$ grid points.
536
537%----------------------------------------------
538\subsection{Boundary geometry}
539\label{BDY_geometry}
540
541Each open boundary set is defined as a list of points. The information
542is stored in the arrays $nbi$, $nbj$, and $nbr$ in the $idx\_bdy$
543structure.  The $nbi$ and $nbj$ arrays
544define the local $(i,j)$ indices of each point in the boundary zone
545and the $nbr$ array defines the discrete distance from the boundary
546with $nbr=1$ meaning that the point is next to the edge of the
547model domain and $nbr>1$ showing that the point is increasingly
548further away from the edge of the model domain. A set of $nbi$, $nbj$,
549and $nbr$ arrays is defined for each of the $T$, $U$ and $V$
550grids. Figure \ref{Fig_LBC_bdy_geom} shows an example of an irregular
551boundary.
552
553The boundary geometry for each set may be defined in a namelist
554nambdy\_index or by reading in a ``coordinates.bdy.nc'' file. The
555nambdy\_index namelist defines a series of straight-line segments for
556north, east, south and west boundaries. For the northern boundary,
557\np{nbdysegn} gives the number of segments, \np{jpjnob} gives the $j$
558index for each segment and \np{jpindt} and \np{jpinft} give the start
559and end $i$ indices for each segment with similar for the other
560boundaries. These segments define a list of $T$ grid points along the
561outermost row of the boundary ($nbr\,=\, 1$). The code deduces the $U$ and
562$V$ points and also the points for $nbr\,>\, 1$ if
563$nn\_rimwidth\,>\,1$.
564
565The boundary geometry may also be defined from a
566``coordinates.bdy.nc'' file. Figure \ref{Fig_LBC_nc_header}
567gives an example of the header information from such a file. The file
568should contain the index arrays for each of the $T$, $U$ and $V$
569grids. The arrays must be in order of increasing $nbr$. Note that the
570$nbi$, $nbj$ values in the file are global values and are converted to
571local values in the code. Typically this file will be used to generate
572external boundary data via interpolation and so will also contain the
573latitudes and longitudes of each point as shown. However, this is not
574necessary to run the model.
575
576For some choices of irregular boundary the model domain may contain
577areas of ocean which are not part of the computational domain. For
578example if an open boundary is defined along an isobath, say at the
579shelf break, then the areas of ocean outside of this boundary will
580need to be masked out. This can be done by reading a mask file defined
581as \np{cn\_mask\_file} in the nam\_bdy namelist. Only one mask file is
582used even if multiple boundary sets are defined.
583
584%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
585\begin{figure}[!t]      \begin{center}
586\includegraphics[width=1.0\textwidth]{Fig_LBC_bdy_geom}
587\caption {      \label{Fig_LBC_bdy_geom}
588Example of geometry of unstructured open boundary}
589\end{center}   \end{figure}
590%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
591
592%----------------------------------------------
593\subsection{Input boundary data files}
594\label{BDY_data}
595
596The data files contain the data arrays
597in the order in which the points are defined in the $nbi$ and $nbj$
598arrays. The data arrays are dimensioned on: a time dimension;
599$xb$ which is the index of the boundary data point in the horizontal;
600and $yb$ which is a degenerate dimension of 1 to enable the file to be
601read by the standard NEMO I/O routines. The 3D fields also have a
602depth dimension.
603
604At Version 3.4 there are new restrictions on the order in which the
605boundary points are defined (and therefore restrictions on the order
606of the data in the file). In particular:
607
608\mbox{}
609
610\begin{enumerate}
611\item The data points must be in order of increasing $nbr$, ie. all
612  the $nbr=1$ points, then all the $nbr=2$ points etc.
613\item All the data for a particular boundary set must be in the same
614  order. (Prior to 3.4 it was possible to define barotropic data in a
615  different order to the data for tracers and baroclinic velocities).
616\end{enumerate}
617
618\mbox{}
619
620These restrictions mean that data files used with previous versions of
621the model may not work with version 3.4. A fortran utility
622{\it bdy\_reorder} exists in the TOOLS directory which will re-order the
623data in old BDY data files.
624
625%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
626\begin{figure}[!t]     \begin{center}
627\includegraphics[width=1.0\textwidth]{Fig_LBC_nc_header}
628\caption {     \label{Fig_LBC_nc_header} 
629Example of the header for a coordinates.bdy.nc file}
630\end{center}   \end{figure}
631%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
632
633%----------------------------------------------
634\subsection{Volume correction}
635\label{BDY_vol_corr}
636
637There is an option to force the total volume in the regional model to be constant,
638similar to the option in the OBC module. This is controlled  by the \np{nn\_volctl} 
639parameter in the namelist. A value of \np{nn\_volctl}~=~0 indicates that this option is not used.
640If  \np{nn\_volctl}~=~1 then a correction is applied to the normal velocities
641around the boundary at each timestep to ensure that the integrated volume flow
642through the boundary is zero. If \np{nn\_volctl}~=~2 then the calculation of
643the volume change on the timestep includes the change due to the freshwater
644flux across the surface and the correction velocity corrects for this as well.
645
646If more than one boundary set is used then volume correction is
647applied to all boundaries at once.
648
649\newpage
650%----------------------------------------------
651\subsection{Tidal harmonic forcing}
652\label{BDY_tides}
653
654%-----------------------------------------nambdy_tide--------------------------------------------
655\namdisplay{nambdy_tide} 
656%-----------------------------------------------------------------------------------------------
657
658Options are defined through the  \ngn{nambdy\_tide} namelist variables.
659 To be written....
660
661
662
663
664\end{document}
Note: See TracBrowser for help on using the repository browser.