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 NEMO/trunk/doc/latex/NEMO/subfiles – NEMO

source: NEMO/trunk/doc/latex/NEMO/subfiles/chap_LBC.tex @ 11577

Last change on this file since 11577 was 11577, checked in by nicolasmartin, 5 years ago

New LaTeX commands \nam and \np to mention namelist content
(Partial commit to serve as a backup before other large edits)
In order to benefit of the syntax highlighting and to have a simpler syntax for
citing namelist block (\nam) and parameter (\np) with an optional variable assignment (\forcode{...}),
at this time the only viable solution I found is to require a double marker for
what it looks like the same item:

  1. Marker with the real name: 'tra_adv' block or 'ln_flx' parameter
  2. Marker with underscore character escaping: 'tra\_adv' block or 'ln\_flx' parameter

Despite many searches and attempts, I did not find a workaround to edit on-the-fly one or
the other marker.
In fact, the problem is on one side that the LaTeX index interprets '_' as a switch for lowering like
in math mode while on the other hand the backslash is considered for Pygments as a typo in Fortran
(red box).

For instance, \nam and \np have as of now the aforementioned 2 mandatory arguments in
the previous order (between braces) + an optional argument for \np when the parameter is defined
(between brackets at the first position):

  • \nam: LaTeX code in the \nam{tra_adv}{tra\_adv} -> PDF ' in the &namtra_adv (namelist X.X) ' with syntax highlighting, the hyperlink and the index entry
  • \np: LaTeX code \np[=.true.]{ln_flx}{ln\_flx} -> PDF ln_flux=.true. with syntax highlighting for the whole string and the entry in the 'parameters' index
File size: 43.5 KB
Line 
1\documentclass[../main/NEMO_manual]{subfiles}
2
3\begin{document}
4% ================================================================
5% Chapter — Lateral Boundary Condition (LBC)
6% ================================================================
7\chapter{Lateral Boundary Condition (LBC)}
8\label{chap:LBC}
9
10\chaptertoc
11
12\newpage
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 (\forcode{rn_shlat})]{Boundary condition at the coast (\protect\np{rn_shlat}{rn\_shlat})}
20\label{sec:LBC_coast}
21%--------------------------------------------namlbc-------------------------------------------------------
22
23\begin{listing}
24  \nlst{namlbc}
25  \caption{\forcode{&namlbc}}
26  \label{lst:namlbc}
27\end{listing}
28%--------------------------------------------------------------------------------------------------------------
29
30%The lateral ocean boundary conditions contiguous to coastlines are Neumann conditions for heat and salt
31%(no flux across boundaries) and Dirichlet conditions for momentum (ranging from free-slip to "strong" no-slip).
32%They are handled automatically by the mask system (see \autoref{subsec:DOM_msk}).
33
34%OPA allows land and topography grid points in the computational domain due to the presence of continents or islands,
35%and includes the use of a full or partial step representation of bottom topography.
36%The computation is performed over the whole domain, \ie\ we do not try to restrict the computation to ocean-only points.
37%This choice has two motivations.
38%Firstly, working on ocean only grid points overloads the code and harms the code readability.
39%Secondly, and more importantly, it drastically reduces the vector portion of the computation,
40%leading to a dramatic increase of CPU time requirement on vector computers.
41%The current section describes how the masking affects the computation of the various terms of the equations
42%with respect to the boundary condition at solid walls.
43%The process of defining which areas are to be masked is described in \autoref{subsec:DOM_msk}.
44
45Options are defined through the \nam{lbc} namelist variables.
46The discrete representation of a domain with complex boundaries (coastlines and bottom topography) leads to
47arrays that include large portions where a computation is not required as the model variables remain at zero.
48Nevertheless, vectorial supercomputers are far more efficient when computing over a whole array,
49and the readability of a code is greatly improved when boundary conditions are applied in
50an automatic way rather than by a specific computation before or after each computational loop.
51An efficient way to work over the whole domain while specifying the boundary conditions,
52is to use multiplication by mask arrays in the computation.
53A mask array is a matrix whose elements are $1$ in the ocean domain and $0$ elsewhere.
54A simple multiplication of a variable by its own mask ensures that it will remain zero over land areas.
55Since most of the boundary conditions consist of a zero flux across the solid boundaries,
56they can be simply applied by multiplying variables by the correct mask arrays,
57\ie\ the mask array of the grid point where the flux is evaluated.
58For example, the heat flux in the \textbf{i}-direction is evaluated at $u$-points.
59Evaluating this quantity as,
60
61\[
62  % \label{eq:LBC_aaaa}
63  \frac{A^{lT} }{e_1 }\frac{\partial T}{\partial i}\equiv \frac{A_u^{lT}
64  }{e_{1u} } \; \delta_{i+1 / 2} \left[ T \right]\;\;mask_u
65\]
66(where mask$_{u}$ is the mask array at a $u$-point) ensures that the heat flux is zero inside land and
67at the boundaries, since mask$_{u}$ is zero at solid boundaries which in this case are defined at $u$-points
68(normal velocity $u$ remains zero at the coast) (\autoref{fig:LBC_uv}).
69
70%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
71\begin{figure}[!t]
72  \centering
73  \includegraphics[width=0.66\textwidth]{Fig_LBC_uv}
74  \caption[Lateral boundary at $T$-level]{
75    Lateral boundary (thick line) at T-level.
76    The velocity normal to the boundary is set to zero.}
77  \label{fig:LBC_uv}
78\end{figure}
79%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
80
81For momentum the situation is a bit more complex as two boundary conditions must be provided along the coast
82(one each for the normal and tangential velocities).
83The boundary of the ocean in the C-grid is defined by the velocity-faces.
84For example, at a given $T$-level,
85the lateral boundary (a coastline or an intersection with the bottom topography) is made of
86segments joining $f$-points, and normal velocity points are located between two $f-$points (\autoref{fig:LBC_uv}).
87The boundary condition on the normal velocity (no flux through solid boundaries)
88can thus be easily implemented using the mask system.
89The boundary condition on the tangential velocity requires a more specific treatment.
90This boundary condition influences the relative vorticity and momentum diffusive trends,
91and is required in order to compute the vorticity at the coast.
92Four different types of lateral boundary condition are available,
93controlled by the value of the \np{rn_shlat}{rn\_shlat} namelist parameter
94(The value of the mask$_{f}$ array along the coastline is set equal to this parameter).
95These are:
96
97%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
98\begin{figure}[!p]
99  \centering
100  \includegraphics[width=0.66\textwidth]{Fig_LBC_shlat}
101  \caption[Lateral boundary conditions]{
102    Lateral boundary conditions
103    (a) free-slip                       (\protect\np{rn_shlat}{rn\_shlat}\forcode{=0});
104    (b) no-slip                         (\protect\np{rn_shlat}{rn\_shlat}\forcode{=2});
105    (c) "partial" free-slip (\forcode{0<}\protect\np{rn_shlat}{rn\_shlat}\forcode{<2}) and
106    (d) "strong" no-slip    (\forcode{2<}\protect\np{rn_shlat}{rn\_shlat}).
107    Implied "ghost" velocity inside land area is display in grey.}
108  \label{fig:LBC_shlat}
109\end{figure}
110%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
111
112\begin{description}
113
114\item[free-slip boundary condition (\np{rn_shlat}{rn\_shlat}\forcode{=0}):] the tangential velocity at
115  the coastline is equal to the offshore velocity,
116  \ie\ the normal derivative of the tangential velocity is zero at the coast,
117  so the vorticity: mask$_{f}$ array is set to zero inside the land and just at the coast
118  (\autoref{fig:LBC_shlat}-a).
119
120\item[no-slip boundary condition (\np{rn_shlat}{rn\_shlat}\forcode{=2}):] the tangential velocity vanishes at the coastline.
121  Assuming that the tangential velocity decreases linearly from
122  the closest ocean velocity grid point to the coastline,
123  the normal derivative is evaluated as if the velocities at the closest land velocity gridpoint and
124  the closest ocean velocity gridpoint were of the same magnitude but in the opposite direction
125  (\autoref{fig:LBC_shlat}-b).
126  Therefore, the vorticity along the coastlines is given by:
127
128  \[
129    \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) \ ,
130  \]
131  where $u$ and $v$ are masked fields.
132  Setting the mask$_{f}$ array to $2$ along the coastline provides a vorticity field computed with
133  the no-slip boundary condition, simply by multiplying it by the mask$_{f}$ :
134  \[
135    % \label{eq:LBC_bbbb}
136    \zeta \equiv \frac{1}{e_{1f} {\kern 1pt}e_{2f} }\left( {\delta_{i+1/2}
137        \left[ {e_{2v} \,v} \right]-\delta_{j+1/2} \left[ {e_{1u} \,u} \right]}
138    \right)\;\mbox{mask}_f
139  \]
140
141\item["partial" free-slip boundary condition (0$<$\np{rn_shlat}{rn\_shlat}$<$2):] the tangential velocity at
142  the coastline is smaller than the offshore velocity, \ie\ there is a lateral friction but
143  not strong enough to make the tangential velocity at the coast vanish (\autoref{fig:LBC_shlat}-c).
144  This can be selected by providing a value of mask$_{f}$ strictly inbetween $0$ and $2$.
145
146\item["strong" no-slip boundary condition (2$<$\np{rn_shlat}{rn\_shlat}):] the viscous boundary layer is assumed to
147  be smaller than half the grid size (\autoref{fig:LBC_shlat}-d).
148  The friction is thus larger than in the no-slip case.
149
150\end{description}
151
152Note that when the bottom topography is entirely represented by the $s$-coordinates (pure $s$-coordinate),
153the lateral boundary condition on tangential velocity is of much less importance as
154it is only applied next to the coast where the minimum water depth can be quite shallow.
155
156
157% ================================================================
158% Boundary Condition around the Model Domain
159% ================================================================
160\section[Model domain boundary condition (\forcode{jperio})]{Model domain boundary condition (\protect\jp{jperio})}
161\label{sec:LBC_jperio}
162
163At the model domain boundaries several choices are offered:
164closed, cyclic east-west, cyclic north-south, a north-fold, and combination closed-north fold or
165bi-cyclic east-west and north-fold.
166The north-fold boundary condition is associated with the 3-pole ORCA mesh.
167
168% -------------------------------------------------------------------------------------------------------------
169%        Closed, cyclic (\jp{jperio}\forcode{ = 0..2})
170% -------------------------------------------------------------------------------------------------------------
171\subsection[Closed, cyclic (\forcode{=0,1,2,7})]{Closed, cyclic (\protect\jp{jperio}\forcode{=0,1,2,7})}
172\label{subsec:LBC_jperio012}
173
174The choice of closed or cyclic model domain boundary condition is made by
175setting \jp{jperio} to 0, 1, 2 or 7 in namelist \nam{cfg}.
176Each time such a boundary condition is needed, it is set by a call to routine \mdl{lbclnk}.
177The computation of momentum and tracer trends proceeds from $i=2$ to $i=jpi-1$ and from $j=2$ to $j=jpj-1$,
178\ie\ in the model interior.
179To choose a lateral model boundary condition is to specify the first and last rows and columns of
180the model variables.
181
182\begin{description}
183
184\item[For closed boundary (\jp{jperio}\forcode{=0})],
185  solid walls are imposed at all model boundaries:
186  first and last rows and columns are set to zero.
187
188\item[For cyclic east-west boundary (\jp{jperio}\forcode{=1})],
189  first and last rows are set to zero (closed) whilst the first column is set to
190  the value of the last-but-one column and the last column to the value of the second one
191  (\autoref{fig:LBC_jperio}-a).
192  Whatever flows out of the eastern (western) end of the basin enters the western (eastern) end.
193
194\item[For cyclic north-south boundary (\jp{jperio}\forcode{=2})],
195  first and last columns are set to zero (closed) whilst the first row is set to
196  the value of the last-but-one row and the last row to the value of the second one
197  (\autoref{fig:LBC_jperio}-a).
198  Whatever flows out of the northern (southern) end of the basin enters the southern (northern) end.
199
200\item[Bi-cyclic east-west and north-south boundary (\jp{jperio}\forcode{=7})] combines cases 1 and 2.
201
202\end{description}
203
204%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
205\begin{figure}[!t]
206  \centering
207  \includegraphics[width=0.66\textwidth]{Fig_LBC_jperio}
208  \caption[Setting of east-west cyclic and symmetric across the Equator boundary conditions]{
209    Setting of (a) east-west cyclic (b) symmetric across the Equator boundary conditions}
210  \label{fig:LBC_jperio}
211\end{figure}
212%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
213
214% -------------------------------------------------------------------------------------------------------------
215%        North fold (\textit{jperio = 3 }to $6)$
216% -------------------------------------------------------------------------------------------------------------
217\subsection[North-fold (\forcode{=3,6})]{North-fold (\protect\jp{jperio}\forcode{=3,6})}
218\label{subsec:LBC_north_fold}
219
220The north fold boundary condition has been introduced in order to handle the north boundary of
221a three-polar ORCA grid.
222Such a grid has two poles in the northern hemisphere (\autoref{fig:CFGS_ORCA_msh},
223and thus requires a specific treatment illustrated in \autoref{fig:LBC_North_Fold_T}.
224Further information can be found in \mdl{lbcnfd} module which applies the north fold boundary condition.
225
226%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
227\begin{figure}[!t]
228  \centering
229  \includegraphics[width=0.66\textwidth]{Fig_North_Fold_T}
230  \caption[North fold boundary in ORCA 2\deg, 1/4\deg and 1/12\deg]{
231    North fold boundary with a $T$-point pivot and cyclic east-west boundary condition ($jperio=4$),
232    as used in ORCA 2\deg, 1/4\deg and 1/12\deg.
233    Pink shaded area corresponds to the inner domain mask (see text).}
234  \label{fig:LBC_North_Fold_T}
235\end{figure}
236%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
237
238% ====================================================================
239% Exchange with neighbouring processors
240% ====================================================================
241\section[Exchange with neighbouring processors (\textit{lbclnk.F90}, \textit{lib\_mpp.F90})]{Exchange with neighbouring processors (\protect\mdl{lbclnk}, \protect\mdl{lib\_mpp})}
242\label{sec:LBC_mpp}
243
244%-----------------------------------------nammpp--------------------------------------------
245
246\begin{listing}
247  \nlst{nammpp}
248  \caption{\forcode{&nammpp}}
249  \label{lst:nammpp}
250\end{listing}
251%-----------------------------------------------------------------------------------------------
252
253For massively parallel processing (mpp), a domain decomposition method is used.
254The basic idea of the method is to split the large computation domain of a numerical experiment into several smaller domains and
255solve the set of equations by addressing independent local problems.
256Each processor has its own local memory and computes the model equation over a subdomain of the whole model domain.
257The subdomain boundary conditions are specified through communications between processors which are organized by
258explicit statements (message passing method).
259The present implementation is largely inspired by Guyon's work [Guyon 1995].
260
261The parallelization strategy is defined by the physical characteristics of the ocean model.
262Second order finite difference schemes lead to local discrete operators that
263depend at the very most on one neighbouring point.
264The only non-local computations concern the vertical physics
265(implicit diffusion, turbulent closure scheme, ...).
266Therefore, a pencil strategy is used for the data sub-structuration:
267the 3D initial domain is laid out on local processor memories following a 2D horizontal topological splitting.
268Each sub-domain computes its own surface and bottom boundary conditions and
269has a side wall overlapping interface which defines the lateral boundary conditions for
270computations in the inner sub-domain.
271The overlapping area consists of the two rows at each edge of the sub-domain.
272After a computation, a communication phase starts:
273each processor sends to its neighbouring processors the update values of the points corresponding to
274the interior overlapping area to its neighbouring sub-domain (\ie\ the innermost of the two overlapping rows).
275Communications are first done according to the east-west direction and next according to the north-south direction.
276There is no specific communications for the corners.
277The communication is done through the Message Passing Interface (MPI) and requires \key{mpp\_mpi}.
278Use also \key{mpi2} if MPI3 is not available on your computer.
279The data exchanges between processors are required at the very place where
280lateral domain boundary conditions are set in the mono-domain computation:
281the \rou{lbc\_lnk} routine (found in \mdl{lbclnk} module) which manages such conditions is interfaced with
282routines found in \mdl{lib\_mpp} module.
283The output file \textit{communication\_report.txt} provides the list of which routines do how
284many communications during 1 time step of the model.\\
285
286%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
287\begin{figure}[!t]
288  \centering
289  \includegraphics[width=0.66\textwidth]{Fig_mpp}
290  \caption{Positioning of a sub-domain when massively parallel processing is used}
291  \label{fig:LBC_mpp}
292\end{figure}
293%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
294
295In \NEMO, the splitting is regular and arithmetic.
296The total number of subdomains corresponds to the number of MPI processes allocated to \NEMO\ when the model is launched
297(\ie\ mpirun -np x ./nemo will automatically give x subdomains).
298The i-axis is divided by \np{jpni} and the j-axis by \np{jpnj}.
299These parameters are defined in \nam{mpp} namelist.
300If \np{jpni} and \np{jpnj} are < 1, they will be automatically redefined in the code to give the best domain decomposition
301(see bellow).
302
303Each processor is independent and without message passing or synchronous process, programs run alone and access just its own local memory.
304For this reason,
305the main model dimensions are now the local dimensions of the subdomain (pencil) that are named \jp{jpi}, \jp{jpj}, \jp{jpk}.
306These dimensions include the internal domain and the overlapping rows.
307The number of rows to exchange (known as the halo) is usually set to one (nn\_hls=1, in \mdl{par\_oce},
308and must be kept to one until further notice).
309The whole domain dimensions are named \jp{jpiglo}, \jp{jpjglo} and \jp{jpk}.
310The relationship between the whole domain and a sub-domain is:
311\begin{gather*}
312  jpi = ( jpiglo-2\times nn\_hls + (jpni-1) ) / jpni + 2\times nn\_hls \\
313  jpj = ( jpjglo-2\times nn\_hls + (jpnj-1) ) / jpnj + 2\times nn\_hls
314\end{gather*}
315
316One also defines variables nldi and nlei which correspond to the internal domain bounds, and the variables nimpp and njmpp which are the position of the (1,1) grid-point in the global domain (\autoref{fig:LBC_mpp}). Note that since the version 4, there is no more extra-halo area as defined in \autoref{fig:LBC_mpp} so \jp{jpi} is now always equal to nlci and \jp{jpj} equal to nlcj.
317
318An element of $T_{l}$, a local array (subdomain) corresponds to an element of $T_{g}$,
319a global array (whole domain) by the relationship:
320\[
321  % \label{eq:LBC_nimpp}
322  T_{g} (i+nimpp-1,j+njmpp-1,k) = T_{l} (i,j,k),
323\]
324with $1 \leq i \leq jpi$, $1  \leq j \leq jpj $ , and  $1  \leq k \leq jpk$.
325
326The 1-d arrays $mig(1:\jp{jpi})$ and $mjg(1:\jp{jpj})$, defined in \rou{dom\_glo} routine (\mdl{domain} module), should be used to get global domain indices from local domain indices. The 1-d arrays, $mi0(1:\jp{jpiglo})$, $mi1(1:\jp{jpiglo})$ and $mj0(1:\jp{jpjglo})$, $mj1(1:\jp{jpjglo})$ have the reverse purpose and should be used to define loop indices expressed in global domain indices (see examples in \mdl{dtastd} module).\\
327
328The \NEMO\ model computes equation terms with the help of mask arrays (0 on land points and 1 on sea points). It is therefore possible that an MPI subdomain contains only land points. To save ressources, we try to supress from the computational domain as much land subdomains as possible. For example if $N_{mpi}$ processes are allocated to NEMO, the domain decomposition will be given by the following equation:
329\[
330  N_{mpi} = jpni \times jpnj - N_{land} + N_{useless}
331\]
332$N_{land}$ is the total number of land subdomains in the domain decomposition defined by \np{jpni} and \np{jpnj}. $N_{useless}$ is the number of land subdomains that are kept in the compuational domain in order to make sure that $N_{mpi}$ MPI processes are indeed allocated to a given subdomain. The values of $N_{mpi}$, \np{jpni}, \np{jpnj}$N_{land}$ and $N_{useless}$ are printed in the output file \texttt{ocean.output}. $N_{useless}$ must, of course, be as small as possible to limit the waste of ressources. A warning is issued in  \texttt{ocean.output} if $N_{useless}$ is not zero. Note that non-zero value of $N_{useless}$ is uselly required when using AGRIF as, up to now, the parent grid and each of the child grids must use all the $N_{mpi}$ processes.
333
334If the domain decomposition is automatically defined (when \np{jpni} and \np{jpnj} are < 1), the decomposition chosen by the model will minimise the sub-domain size (defined as $max_{all domains}(jpi \times jpj)$) and maximize the number of eliminated land subdomains. This means that no other domain decomposition (a set of \np{jpni} and \np{jpnj} values) will use less processes than $(jpni  \times  jpnj - N_{land})$ and get a smaller subdomain size.
335In order to specify $N_{mpi}$ properly (minimize $N_{useless}$), you must run the model once with \np{ln_list}{ln\_list} activated. In this case, the model will start the initialisation phase, print the list of optimum decompositions ($N_{mpi}$, \np{jpni} and \np{jpnj}) in \texttt{ocean.output} and directly abort. The maximum value of $N_{mpi}$ tested in this list is given by $max(N_{MPI\_tasks}, jpni \times jpnj)$. For example, run the model on 40 nodes with ln\_list activated and $jpni = 10000$ and $jpnj = 1$, will print the list of optimum domains decomposition from 1 to about 10000.
336
337Processors are numbered from 0 to $N_{mpi} - 1$. Subdomains containning some ocean points are numbered first from 0 to $jpni * jpnj - N_{land} -1$. The remaining $N_{useless}$ land subdomains are numbered next, which means that, for a given (\np{jpni}, \np{jpnj}), the numbers attributed to he ocean subdomains do not vary with $N_{useless}$.
338
339When land processors are eliminated, the value corresponding to these locations in the model output files is undefined. \np{ln_mskland}{ln\_mskland} must be activated in order avoid Not a Number values in output files. Note that it is better to not eliminate land processors when creating a meshmask file (\ie\ when setting a non-zero value to \np{nn_msh}{nn\_msh}).
340
341%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
342\begin{figure}[!ht]
343  \centering
344  \includegraphics[width=0.66\textwidth]{Fig_mppini2}
345  \caption[Atlantic domain defined for the CLIPPER projet]{
346    Example of Atlantic domain defined for the CLIPPER projet.
347    Initial grid is composed of 773 x 1236 horizontal points.
348    (a) the domain is split onto 9 $times$ 20 subdomains (jpni=9, jpnj=20).
349    52 subdomains are land areas.
350    (b) 52 subdomains are eliminated (white rectangles) and
351    the resulting number of processors really used during the computation is jpnij=128.}
352  \label{fig:LBC_mppini2}
353\end{figure}
354%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
355
356
357% ====================================================================
358% Unstructured open boundaries BDY
359% ====================================================================
360\section{Unstructured open boundary conditions (BDY)}
361\label{sec:LBC_bdy}
362
363%-----------------------------------------nambdy--------------------------------------------
364
365\begin{listing}
366  \nlst{nambdy}
367  \caption{\forcode{&nambdy}}
368  \label{lst:nambdy}
369\end{listing}
370%-----------------------------------------------------------------------------------------------
371%-----------------------------------------nambdy_dta--------------------------------------------
372
373\begin{listing}
374  \nlst{nambdy_dta}
375  \caption{\forcode{&nambdy_dta}}
376  \label{lst:nambdy_dta}
377\end{listing}
378%-----------------------------------------------------------------------------------------------
379
380Options are defined through the \nam{bdy} and \nam{bdy_dta}{bdy\_dta} namelist variables.
381The BDY module is the core implementation of open boundary conditions for regional configurations on
382ocean temperature, salinity, barotropic-baroclinic velocities, ice-snow concentration, thicknesses, temperatures, salinity and melt ponds concentration and thickness.
383
384The BDY module was modelled on the OBC module (see \NEMO\ 3.4) and shares many features and
385a similar coding structure \citep{chanut_rpt05}.
386The specification of the location of the open boundary is completely flexible and
387allows any type of setup, from regular boundaries to irregular contour (it includes the possibility to set an open boundary able to follow an isobath).
388Boundary data files used with versions of \NEMO\ prior to Version 3.4 may need to be re-ordered to work with this version.
389See the section on the Input Boundary Data Files for details.
390
391%----------------------------------------------
392\subsection{Namelists}
393\label{subsec:LBC_bdy_namelist}
394
395The BDY module is activated by setting \np{ln_bdy}{ln\_bdy}\forcode{=.true.} .
396It is possible to define more than one boundary ``set'' and apply different boundary conditions to each set.
397The number of boundary sets is defined by \np{nb_bdy}{nb\_bdy}.
398Each boundary set can be either defined as a series of straight line segments directly in the namelist
399(\np{ln_coords_file}{ln\_coords\_file}\forcode{=.false.}, and a namelist block \nam{bdy_index}{bdy\_index} must be included for each set) or read in from a file (\np{ln\_coords\_file}\forcode{=.true.}, and a ``\ifile{coordinates.bdy}'' file must be provided).
400The coordinates.bdy file is analagous to the usual \NEMO\ ``\ifile{coordinates}'' file.
401In the example above, there are two boundary sets, the first of which is defined via a file and
402the second is defined in the namelist.
403For more details of the definition of the boundary geometry see section \autoref{subsec:LBC_bdy_geometry}.
404
405For each boundary set a boundary condition has to be chosen for the barotropic solution
406(``u2d'':sea-surface height and barotropic velocities), for the baroclinic velocities (``u3d''),
407for the active tracers \footnote{The BDY module does not deal with passive tracers at this version} (``tra''), and for sea-ice (``ice'').
408For each set of variables one has to choose an algorithm and the boundary data (set resp. by \np{cn_tra}{cn\_tra} and \np{nn_tra_dta}{nn\_tra\_dta} for tracers).\\
409
410The choice of algorithm is currently as follows:
411
412\begin{description}
413\item[\forcode{'none'}:] No boundary condition applied.
414  So the solution will ``see'' the land points around the edge of the edge of the domain.
415\item[\forcode{'specified'}:] Specified boundary condition applied (only available for baroclinic velocity and tracer variables).
416\item[\forcode{'neumann'}:] Value at the boundary are duplicated (No gradient). Only available for baroclinic velocity and tracer variables.
417\item[\forcode{'frs'}:] Flow Relaxation Scheme (FRS) available for all variables.
418\item[\forcode{'Orlanski'}:] Orlanski radiation scheme (fully oblique) for barotropic, baroclinic and tracer variables.
419\item[\forcode{'Orlanski_npo'}:] Orlanski radiation scheme for barotropic, baroclinic and tracer variables.
420\item[\forcode{'flather'}:] Flather radiation scheme for the barotropic variables only.
421\end{description}
422
423The boundary data is either set to initial conditions
424(\np{nn_tra_dta}{nn\_tra\_dta}\forcode{=0}) or forced with external data from a file (\np{nn\_tra\_dta}\forcode{=1}).
425In case the 3d velocity data contain the total velocity (ie, baroclinic and barotropic velocity),
426the bdy code can derived baroclinic and barotropic velocities by setting \np{ln_full_vel}{ln\_full\_vel}\forcode{=.true.}
427For the barotropic solution there is also the option to use tidal harmonic forcing either by
428itself (\np{nn_dyn2d_dta}{nn\_dyn2d\_dta}\forcode{=2}) or in addition to other external data (\np{nn\_dyn2d\_dta}\forcode{=3}).\\
429If not set to initial conditions, sea-ice salinity, temperatures and melt ponds data at the boundary can either be read in a file or defined as constant (by \np{rn_ice_sal}{rn\_ice\_sal}, \np{rn_ice_tem}{rn\_ice\_tem}, \np{rn_ice_apnd}{rn\_ice\_apnd}, \np{rn_ice_hpnd}{rn\_ice\_hpnd}). Ice age is constant and defined by \np{rn_ice_age}{rn\_ice\_age}.
430
431If external boundary data is required then the \nam{bdy_dta}{bdy\_dta} namelist must be defined.
432One \nam{bdy_dta}{bdy\_dta} namelist is required for each boundary set, adopting the same order of indexes in which the boundary sets are defined in nambdy.
433In the example given, two boundary sets have been defined. The first one is reading data file in the \nam{bdy_dta}{bdy\_dta} namelist shown above
434and the second one is using data from intial condition (no namelist block needed).
435The boundary data is read in using the fldread module,
436so the \nam{bdy_dta}{bdy\_dta} namelist is in the format required for fldread.
437For each required variable, the filename, the frequency of the files and
438the frequency of the data in the files are given.
439Also whether or not time-interpolation is required and whether the data is climatological (time-cyclic) data.
440For sea-ice salinity, temperatures and melt ponds, reading the files are skipped and constant values are used if filenames are defined as {'NOT USED'}.\\
441
442There is currently an option to vertically interpolate the open boundary data onto the native grid at run-time.
443If \np{nn_bdy_jpk}{nn\_bdy\_jpk}$<-1$, it is assumed that the lateral boundary data are already on the native grid.
444However, if \np{nn_bdy_jpk}{nn\_bdy\_jpk} is set to the number of vertical levels present in the boundary data,
445a bilinear interpolation onto the native grid will be triggered at runtime.
446For this to be successful the additional variables: $gdept$, $gdepu$, $gdepv$, $e3t$, $e3u$ and $e3v$, are required to be present in the lateral boundary files.
447These correspond to the depths and scale factors of the input data,
448the latter used to make any adjustment to the velocity fields due to differences in the total water depths between the two vertical grids.\\
449
450In the example of given namelists, two boundary sets are defined.
451The first set is defined via a file and applies FRS conditions to temperature and salinity and
452Flather conditions to the barotropic variables. No condition specified for the baroclinic velocity and sea-ice.
453External data is provided in daily files (from a large-scale model).
454Tidal harmonic forcing is also used.
455The second set is defined in a namelist.
456FRS conditions are applied on temperature and salinity and climatological data is read from initial condition files.
457
458%----------------------------------------------
459\subsection{Flow relaxation scheme}
460\label{subsec:LBC_bdy_FRS_scheme}
461
462The Flow Relaxation Scheme (FRS) \citep{davies_QJRMS76,engedahl_T95},
463applies a simple relaxation of the model fields to externally-specified values over
464a zone next to the edge of the model domain.
465Given a model prognostic variable $\Phi$
466\[
467  % \label{eq:LBC_bdy_frs1}
468  \Phi(d) = \alpha(d)\Phi_{e}(d) + (1-\alpha(d))\Phi_{m}(d)\;\;\;\;\; d=1,N
469\]
470where $\Phi_{m}$ is the model solution and $\Phi_{e}$ is the specified external field,
471$d$ gives the discrete distance from the model boundary and
472$\alpha$ is a parameter that varies from $1$ at $d=1$ to a small value at $d=N$.
473It can be shown that this scheme is equivalent to adding a relaxation term to
474the prognostic equation for $\Phi$ of the form:
475\[
476  % \label{eq:LBC_bdy_frs2}
477  -\frac{1}{\tau}\left(\Phi - \Phi_{e}\right)
478\]
479where the relaxation time scale $\tau$ is given by a function of $\alpha$ and the model time step $\Delta t$:
480\[
481  % \label{eq:LBC_bdy_frs3}
482  \tau = \frac{1-\alpha}{\alpha\,\rdt
483\]
484Thus the model solution is completely prescribed by the external conditions at the edge of the model domain and
485is relaxed towards the external conditions over the rest of the FRS zone.
486The application of a relaxation zone helps to prevent spurious reflection of
487outgoing signals from the model boundary.
488
489The function $\alpha$ is specified as a $tanh$ function:
490\[
491  % \label{eq:LBC_bdy_frs4}
492  \alpha(d) = 1 - \tanh\left(\frac{d-1}{2}\right),       \quad d=1,N
493\]
494The width of the FRS zone is specified in the namelist as \np{nn_rimwidth}{nn\_rimwidth}.
495This is typically set to a value between 8 and 10.
496
497%----------------------------------------------
498\subsection{Flather radiation scheme}
499\label{subsec:LBC_bdy_flather_scheme}
500
501The \citet{flather_JPO94} scheme is a radiation condition on the normal,
502depth-mean transport across the open boundary.
503It takes the form
504\begin{equation}
505  \label{eq:LBC_bdy_fla1}
506  U = U_{e} + \frac{c}{h}\left(\eta - \eta_{e}\right),
507\end{equation}
508where $U$ is the depth-mean velocity normal to the boundary and $\eta$ is the sea surface height,
509both from the model.
510The subscript $e$ indicates the same fields from external sources.
511The speed of external gravity waves is given by $c = \sqrt{gh}$, and $h$ is the depth of the water column.
512The depth-mean normal velocity along the edge of the model domain is set equal to
513the external depth-mean normal velocity,
514plus a correction term that allows gravity waves generated internally to exit the model boundary.
515Note that the sea-surface height gradient in \autoref{eq:LBC_bdy_fla1} is a spatial gradient across the model boundary,
516so that $\eta_{e}$ is defined on the $T$ points with $nbr=1$ and $\eta$ is defined on the $T$ points with $nbr=2$.
517$U$ and $U_{e}$ are defined on the $U$ or $V$ points with $nbr=1$, \ie\ between the two $T$ grid points.
518
519%----------------------------------------------
520\subsection{Orlanski radiation scheme}
521\label{subsec:LBC_bdy_orlanski_scheme}
522
523The Orlanski scheme is based on the algorithm described by \citep{marchesiello.mcwilliams.ea_OM01}, hereafter MMS.
524
525The adaptive Orlanski condition solves a wave plus relaxation equation at the boundary:
526\begin{equation}
527  \label{eq:LBC_wave_continuous}
528  \frac{\partial\phi}{\partial t} + c_x \frac{\partial\phi}{\partial x} + c_y \frac{\partial\phi}{\partial y} =
529  -\frac{1}{\tau}(\phi - \phi^{ext})
530\end{equation}
531
532where $\phi$ is the model field, $x$ and $y$ refer to the normal and tangential directions to the boundary respectively, and the phase
533velocities are diagnosed from the model fields as:
534
535\begin{equation}
536  \label{eq:LBC_cx}
537  c_x = -\frac{\partial\phi}{\partial t}\frac{\partial\phi / \partial x}{(\partial\phi /\partial x)^2 + (\partial\phi /\partial y)^2}
538\end{equation}
539\begin{equation}
540  \label{eq:LBC_cy}
541  c_y = -\frac{\partial\phi}{\partial t}\frac{\partial\phi / \partial y}{(\partial\phi /\partial x)^2 + (\partial\phi /\partial y)^2}
542\end{equation}
543
544(As noted by MMS, this is a circular diagnosis of the phase speeds which only makes sense on a discrete grid).
545Equation (\autoref{eq:LBC_wave_continuous}) is defined adaptively depending on the sign of the phase velocity normal to the boundary $c_x$.
546For $c_x$ outward, we have
547
548\begin{equation}
549\tau = \tau_{out}
550\end{equation}
551
552For $c_x$ inward, the radiation equation is not applied:
553
554\begin{equation}
555  \label{eq:LBC_tau_in}
556  \tau = \tau_{in}\,\,\,;\,\,\, c_x = c_y = 0
557\end{equation}
558
559Generally the relaxation time scale at inward propagation points (\np{rn_time_dmp}{rn\_time\_dmp}) is set much shorter than the time scale at outward propagation
560points (\np{rn_time_dmp_out}{rn\_time\_dmp\_out}) so that the solution is constrained more strongly by the external data at inward propagation points.
561See \autoref{subsec:LBC_bdy_relaxation} for detailed on the spatial shape of the scaling.\\
562The ``normal propagation of oblique radiation'' or NPO approximation (called \forcode{'orlanski_npo'}) involves assuming
563that $c_y$ is zero in equation (\autoref{eq:LBC_wave_continuous}), but including
564this term in the denominator of equation (\autoref{eq:LBC_cx}). Both versions of the scheme are options in BDY. Equations
565(\autoref{eq:LBC_wave_continuous}) - (\autoref{eq:LBC_tau_in}) correspond to equations (13) - (15) and (2) - (3) in MMS.\\
566
567%----------------------------------------------
568\subsection{Relaxation at the boundary}
569\label{subsec:LBC_bdy_relaxation}
570
571In addition to a specific boundary condition specified as \np{cn_tra}{cn\_tra} and \np{cn_dyn3d}{cn\_dyn3d}, relaxation on baroclinic velocities and tracers variables are available.
572It is control by the namelist parameter \np{ln_tra_dmp}{ln\_tra\_dmp} and \np{ln_dyn3d_dmp}{ln\_dyn3d\_dmp} for each boundary set.
573
574The relaxation time scale value (\np{rn_time_dmp}{rn\_time\_dmp} and \np{rn_time_dmp_out}{rn\_time\_dmp\_out}, $\tau$) are defined at the boundaries itself.
575This time scale ($\alpha$) is weighted by the distance ($d$) from the boundary over \np{nn_rimwidth}{nn\_rimwidth} cells ($N$):
576
577\[
578  \alpha = \frac{1}{\tau}(\frac{N+1-d}{N})^2,       \quad d=1,N
579\]
580
581The same scaling is applied in the Orlanski damping.
582
583%----------------------------------------------
584\subsection{Boundary geometry}
585\label{subsec:LBC_bdy_geometry}
586
587Each open boundary set is defined as a list of points.
588The information is stored in the arrays $nbi$, $nbj$, and $nbr$ in the $idx\_bdy$ structure.
589The $nbi$ and $nbj$ arrays define the local $(i,j)$ indexes of each point in the boundary zone and
590the $nbr$ array defines the discrete distance from the boundary: $nbr=1$ means that
591the boundary point is next to the edge of the model domain, while $nbr>1$ means that
592the boundary point is increasingly further away from the edge of the model domain.
593A set of $nbi$, $nbj$, and $nbr$ arrays is defined for each of the $T$, $U$ and $V$ grids.
594\autoref{fig:LBC_bdy_geom} shows an example of an irregular boundary.
595
596The boundary geometry for each set may be defined in a namelist nambdy\_index or
597by reading in a ``\ifile{coordinates.bdy}'' file.
598The nambdy\_index namelist defines a series of straight-line segments for north, east, south and west boundaries.
599One nambdy\_index namelist block is needed for each boundary condition defined by indexes.
600For the northern boundary, \texttt{nbdysegn} gives the number of segments,
601\jp{jpjnob} gives the $j$ index for each segment and \jp{jpindt} and
602\jp{jpinft} give the start and end $i$ indices for each segment with similar for the other boundaries.
603These segments define a list of $T$ grid points along the outermost row of the boundary ($nbr\,=\, 1$).
604The code deduces the $U$ and $V$ points and also the points for $nbr\,>\, 1$ if \np{nn_rimwidth}{nn\_rimwidth}\forcode{>1}.
605
606The boundary geometry may also be defined from a ``\ifile{coordinates.bdy}'' file.
607\autoref{fig:LBC_nc_header} gives an example of the header information from such a file, based on the description of geometrical setup given above.
608The file should contain the index arrays for each of the $T$, $U$ and $V$ grids.
609The arrays must be in order of increasing $nbr$.
610Note that the $nbi$, $nbj$ values in the file are global values and are converted to local values in the code.
611Typically this file will be used to generate external boundary data via interpolation and so
612will also contain the latitudes and longitudes of each point as shown.
613However, this is not necessary to run the model.
614
615For some choices of irregular boundary the model domain may contain areas of ocean which
616are not part of the computational domain.
617For example, if an open boundary is defined along an isobath, say at the shelf break,
618then the areas of ocean outside of this boundary will need to be masked out.
619This can be done by reading a mask file defined as \np{cn_mask_file}{cn\_mask\_file} in the nam\_bdy namelist.
620Only one mask file is used even if multiple boundary sets are defined.
621
622%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
623\begin{figure}[!t]
624  \centering
625  \includegraphics[width=0.66\textwidth]{Fig_LBC_bdy_geom}
626  \caption[Geometry of unstructured open boundary]{Example of geometry of unstructured open boundary}
627  \label{fig:LBC_bdy_geom}
628\end{figure}
629%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
630
631%----------------------------------------------
632\subsection{Input boundary data files}
633\label{subsec:LBC_bdy_data}
634
635The data files contain the data arrays in the order in which the points are defined in the $nbi$ and $nbj$ arrays.
636The data arrays are dimensioned on:
637a time dimension;
638$xb$ which is the index of the boundary data point in the horizontal;
639and $yb$ which is a degenerate dimension of 1 to enable the file to be read by the standard \NEMO\ I/O routines.
640The 3D fields also have a depth dimension.
641
642From Version 3.4 there are new restrictions on the order in which the boundary points are defined
643(and therefore restrictions on the order of the data in the file).
644In particular:
645
646\begin{enumerate}
647\item The data points must be in order of increasing $nbr$,
648  ie. all the $nbr=1$ points, then all the $nbr=2$ points etc.
649\item All the data for a particular boundary set must be in the same order.
650  (Prior to 3.4 it was possible to define barotropic data in a different order to
651  the data for tracers and baroclinic velocities).
652\end{enumerate}
653
654These restrictions mean that data files used with versions of the
655model prior to Version 3.4 may not work with Version 3.4 onwards.
656A \fortran\ utility {\itshape bdy\_reorder} exists in the TOOLS directory which
657will re-order the data in old BDY data files.
658
659%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
660\begin{figure}[!t]
661  \centering
662  \includegraphics[width=0.66\textwidth]{Fig_LBC_nc_header}
663  \caption[Header for a \protect\ifile{coordinates.bdy} file]{
664    Example of the header for a \protect\ifile{coordinates.bdy} file}
665  \label{fig:LBC_nc_header}
666\end{figure}
667%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
668
669%----------------------------------------------
670\subsection{Volume correction}
671\label{subsec:LBC_bdy_vol_corr}
672
673There is an option to force the total volume in the regional model to be constant.
674This is controlled  by the \np{ln_vol}{ln\_vol} parameter in the namelist.
675A value of \np{ln_vol}{ln\_vol}\forcode{=.false.} indicates that this option is not used.
676Two options to control the volume are available (\np{nn_volctl}{nn\_volctl}).
677If \np{nn_volctl}{nn\_volctl}\forcode{=0} then a correction is applied to the normal barotropic velocities around the boundary at
678each timestep to ensure that the integrated volume flow through the boundary is zero.
679If \np{nn_volctl}{nn\_volctl}\forcode{=1} then the calculation of the volume change on
680the timestep includes the change due to the freshwater flux across the surface and
681the correction velocity corrects for this as well.
682
683If more than one boundary set is used then volume correction is
684applied to all boundaries at once.
685
686%----------------------------------------------
687\subsection{Tidal harmonic forcing}
688\label{subsec:LBC_bdy_tides}
689
690%-----------------------------------------nambdy_tide--------------------------------------------
691
692\begin{listing}
693  \nlst{nambdy_tide}
694  \caption{\forcode{&nambdy_tide}}
695  \label{lst:nambdy_tide}
696\end{listing}
697%-----------------------------------------------------------------------------------------------
698
699Tidal forcing at open boundaries requires the activation of surface
700tides (i.e., in \nam{_tide}{\_tide}, \np{ln_tide}{ln\_tide} needs to be set to
701\forcode{.true.} and the required constituents need to be activated by
702including their names in the \np{clname} array; see
703\autoref{sec:SBC_tide}). Specific options related to the reading in of
704the complex harmonic amplitudes of elevation (SSH) and barotropic
705velocity (u,v) at open boundaries are defined through the
706\nam{bdy_tide}{bdy\_tide} namelist parameters.\\
707
708The tidal harmonic data at open boundaries can be specified in two
709different ways, either on a two-dimensional grid covering the entire
710model domain or along open boundary segments; these two variants can
711be selected by setting \np{ln_bdytide_2ddta }{ln\_bdytide\_2ddta } to \forcode{.true.} or
712\forcode{.false.}, respectively. In either case, the real and
713imaginary parts of SSH and the two barotropic velocity components for
714each activated tidal constituent \textit{tcname} have to be provided
715separately: when two-dimensional data is used, variables
716\textit{tcname\_z1} and \textit{tcname\_z2} for real and imaginary SSH,
717respectively, are expected in input file \np{filtide} with suffix
718\ifile{\_grid\_T}, variables \textit{tcname\_u1} and
719\textit{tcname\_u2} for real and imaginary u, respectively, are
720expected in input file \np{filtide} with suffix \ifile{\_grid\_U}, and
721\textit{tcname\_v1} and \textit{tcname\_v2} for real and imaginary v,
722respectively, are expected in input file \np{filtide} with suffix
723\ifile{\_grid\_V}; when data along open boundary segments is used,
724variables \textit{z1} and \textit{z2} (real and imaginary part of SSH)
725are expected to be available from file \np{filtide} with suffix
726\ifile{tcname\_grid\_T}, variables \textit{u1} and \textit{u2} (real
727and imaginary part of u) are expected to be available from file
728\np{filtide} with suffix \ifile{tcname\_grid\_U}, and variables
729\textit{v1} and \textit{v2} (real and imaginary part of v) are
730expected to be available from file \np{filtide} with suffix
731\ifile{tcname\_grid\_V}. If \np{ln_bdytide_conj}{ln\_bdytide\_conj} is set to
732\forcode{.true.}, the data is expected to be in complex conjugate
733form.
734
735Note that the barotropic velocity components are assumed to be defined
736on the native model grid and should be rotated accordingly when they
737are converted from their definition on a different source grid. To do
738so, the u, v amplitudes and phases can be converted into tidal
739ellipses, the grid rotation added to the ellipse inclination, and then
740converted back (care should be taken regarding conventions of the
741direction of rotation). %, e.g. anticlockwise or clockwise.
742
743\biblio
744
745\pindex
746
747\end{document}
Note: See TracBrowser for help on using the repository browser.