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_misc.tex in NEMO/trunk/doc/latex/NEMO/subfiles – NEMO

source: NEMO/trunk/doc/latex/NEMO/subfiles/chap_misc.tex @ 11339

Last change on this file since 11339 was 11339, checked in by acc, 5 years ago

Update of chap_misc.tex for v4.0 compatibility

File size: 22.1 KB
Line 
1\documentclass[../main/NEMO_manual]{subfiles}
2
3\begin{document}
4% ================================================================
5% Chapter --- Miscellaneous Topics
6% ================================================================
7\chapter{Miscellaneous Topics}
8\label{chap:MISC}
9
10\minitoc
11
12\newpage
13
14% ================================================================
15% Representation of Unresolved Straits
16% ================================================================
17\section{Representation of unresolved straits}
18\label{sec:MISC_strait}
19
20In climate modeling, it often occurs that a crucial connections between water masses is broken as
21the grid mesh is too coarse to resolve narrow straits.
22For example, coarse grid spacing typically closes off the Mediterranean from the Atlantic at
23the Strait of Gibraltar.
24In this case, it is important for climate models to include the effects of salty water entering the Atlantic from
25the Mediterranean.
26Likewise, it is important for the Mediterranean to replenish its supply of water from the Atlantic to
27balance the net evaporation occurring over the Mediterranean region.
28This problem occurs even in eddy permitting simulations.
29For example, in ORCA 1/4\deg several straits of the Indonesian archipelago (Ombai, Lombok...)
30are much narrow than even a single ocean grid-point.
31
32We describe briefly here the two methods that can be used in \NEMO to handle such
33improperly resolved straits. The methods consist of opening the strait while ensuring
34that the mass exchanges through the strait are not too large by either artificially
35reducing the cross-sectional area of the strait grid-cells or, locally increasing the
36lateral friction.
37
38% -------------------------------------------------------------------------------------------------------------
39%       Hand made geometry changes
40% -------------------------------------------------------------------------------------------------------------
41\subsection{Hand made geometry changes}
42\label{subsec:MISC_strait_hand}
43
44The first method involves reducing the scale factor in the cross-strait direction to a
45value in better agreement with the true mean width of the strait
46(\autoref{fig:MISC_strait_hand}).  This technique is sometime called "partially open face"
47or "partially closed cells".  The key issue here is only to reduce the faces of $T$-cell
48(\ie change the value of the horizontal scale factors at $u$- or $v$-point) but not the
49volume of the $T$-cell.  Indeed, reducing the volume of strait $T$-cell can easily produce
50a numerical instability at that grid point which would require a reduction of the model
51time step.  Thus to instigate a local change in the width of a Strait requires two steps:
52
53\begin{itemize}
54
55\item Add \texttt{e1e2u} and \texttt{e1e2v} arrays to the \np{cn\_domcfg} file. These 2D
56arrays should contain the products of the unaltered values of: $\texttt{e1u}*\texttt{e2u}$
57and $\texttt{e1u}*\texttt{e2v}$ respectively. That is the original surface areas of $u$-
58and $v$- cells respectively.  These areas are usually defined by the corresponding product
59within the \NEMO code but the presence of \texttt{e1e2u} and \texttt{e1e2v} in the
60\np{cn\_domcfg} file will suppress this calculation and use the supplied fields instead.
61If the model domain is provided by user-supplied code in \mdl{usrdef\_hgr}, then this
62routine should also return \texttt{e1e2u} and \texttt{e1e2v} and set the integer return
63argument \texttt{ie1e2u\_v} to a non-zero value. Values other than 0 for this argument
64will suppress the calculation of the areas.
65
66\item Change values of \texttt{e2u} or \texttt{e1v} (either in the \np{cn\_domcfg} file or
67via code in  \mdl{usrdef\_hgr}), whereever a Strait reduction is required. The choice of
68whether to alter \texttt{e2u} or \texttt{e1v} depends. respectively,  on whether the
69Strait in question is North-South orientated (\eg Gibraltar) or East-West orientated (\eg
70Lombok).
71
72\end{itemize}
73
74
75The second method is to increase the viscous boundary layer thickness by a local increase
76of the fmask value at the coast. This method can also be effective in wider passages.  The
77concept is illustarted in the second part of  \autoref{fig:MISC_strait_hand} and changes
78to specific locations can be coded in \mdl{usrdef\_fmask}. The \forcode{usr_def_fmask}
79routine is always called after \texttt{fmask} has been defined according to the choice of
80lateral boundary condition as discussed in \autoref{sec:LBC_coast}. The default version of
81\mdl{usrdef\_fmask} contains settings specific to ORCA2 and ORCA1 configurations. These are
82meant as examples only; it is up to the user to verify settings and provide alternatives
83for their own configurations. The default \forcode{usr_def_fmask} makes no changes to
84\texttt{fmask} for any other configuration.
85
86%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
87\begin{figure}[!tbp]
88  \begin{center}
89    \includegraphics[width=\textwidth]{Fig_Gibraltar}
90    \includegraphics[width=\textwidth]{Fig_Gibraltar2}
91    \caption{
92      \protect\label{fig:MISC_strait_hand}
93      Example of the Gibraltar strait defined in a $1^{\circ} \times 1^{\circ}$ mesh.
94      \textit{Top}: using partially open cells.
95      The meridional scale factor at $v$-point is reduced on both sides of the strait to account for
96      the real width of the strait (about 20 km).
97      Note that the scale factors of the strait $T$-point remains unchanged.
98      \textit{Bottom}: using viscous boundary layers.
99      The four fmask parameters along the strait coastlines are set to a value larger than 4,
100      \ie "strong" no-slip case (see \autoref{fig:LBC_shlat}) creating a large viscous boundary layer that
101      allows a reduced transport through the strait.
102    }
103  \end{center}
104\end{figure}
105%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
106
107%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
108\begin{figure}[!tbp]
109  \begin{center}
110    \includegraphics[width=\textwidth]{Fig_closea_mask_example}
111    \caption{
112      \protect\label{fig:closea_mask_example}
113      Example of mask fields for the closea module. \textit{Left}: a
114      closea\_mask field; \textit{Right}: a closea\_mask\_rnf
115      field. In this example, if ln\_closea is set to .true., the mean
116      freshwater flux over each of the American Great Lakes will be
117      set to zero, and the total residual for all the lakes, if
118      negative, will be put into the St Laurence Seaway in the area
119      shown.
120    }
121  \end{center}
122\end{figure}
123%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
124
125% ================================================================
126% Closed seas
127% ================================================================
128\section[Closed seas (\textit{closea.F90})]
129{Closed seas (\protect\mdl{closea})}
130\label{sec:MISC_closea}
131
132Some configurations include inland seas and lakes as ocean
133points. This is particularly the case for configurations that are
134coupled to an atmosphere model where one might want to include inland
135seas and lakes as ocean model points in order to provide a better
136bottom boundary condition for the atmosphere. However there is no
137route for freshwater to run off from the lakes to the ocean and this
138can lead to large drifts in the sea surface height over the lakes. The
139closea module provides options to either fill in closed seas and lakes
140at run time, or to set the net surface freshwater flux for each lake
141to zero and put the residual flux into the ocean.
142
143Prior to NEMO 4 the locations of inland seas and lakes was set via
144hardcoded indices for various ORCA configurations. From NEMO 4 onwards
145the inland seas and lakes are defined using mask fields in the
146domain configuration file. The options are as follows.
147
148\begin{enumerate}
149\item{{\bfseries No ``closea\_mask'' field is included in domain configuration
150  file.} In this case the closea module does nothing.}
151
152\item{{\bfseries A field called closea\_mask is included in the domain
153configuration file and ln\_closea=.false. in namelist namcfg.} In this
154case the inland seas defined by the closea\_mask field are filled in
155(turned to land points) at run time. That is every point in
156closea\_mask that is nonzero is set to be a land point.}
157
158\item{{\bfseries A field called closea\_mask is included in the domain
159configuration file and ln\_closea=.true. in namelist namcfg.} Each
160inland sea or group of inland seas is set to a positive integer value
161in the closea\_mask field (see Figure \ref{fig:closea_mask_example}
162for an example). The net surface flux over each inland sea or group of
163inland seas is set to zero each timestep and the residual flux is
164distributed over the global ocean (ie. all ocean points where
165closea\_mask is zero).}
166
167\item{{\bfseries Fields called closea\_mask and closea\_mask\_rnf are
168included in the domain configuration file and ln\_closea=.true. in
169namelist namcfg.} This option works as for option 3, except that if
170the net surface flux over an inland sea is negative (net
171precipitation) it is put into the ocean at specified runoff points. A
172net positive surface flux (net evaporation) is still spread over the
173global ocean. The mapping from inland seas to runoff points is defined
174by the closea\_mask\_rnf field. Each mapping is defined by a positive
175integer value for the inland sea(s) and the corresponding runoff
176points. An example is given in Figure
177\ref{fig:closea_mask_example}. If no mapping is provided for a
178particular inland sea then the residual is spread over the global
179ocean.}
180
181\item{{\bfseries Fields called closea\_mask and closea\_mask\_emp are
182included in the domain configuration file and ln\_closea=.true. in
183namelist namcfg.} This option works the same as option 4 except that
184the nonzero net surface flux is sent to the ocean at the specified
185runoff points regardless of whether it is positive or negative. The
186mapping from inland seas to runoff points in this case is defined by
187the closea\_mask\_emp field.}
188\end{enumerate}
189
190There is a python routine to create the closea\_mask fields and append
191them to the domain configuration file in the utils/tools/DOMAINcfg directory.
192
193% ================================================================
194% Sub-Domain Functionality
195% ================================================================
196\section{Sub-domain functionality}
197\label{sec:MISC_zoom}
198
199\subsection{Simple subsetting of input files via NetCDF attributes}
200
201The extended grids for use with the under-shelf ice cavities will result in redundant rows
202around Antarctica if the ice cavities are not active.  A simple mechanism for subsetting
203input files associated with the extended domains has been implemented to avoid the need to
204maintain different sets of input fields for use with or without active ice cavities.  This
205subsetting operates for the j-direction only and works by optionally looking for and using
206a global file attribute (named: \np{open\_ocean\_jstart}) to determine the starting j-row
207for input.  The use of this option is best explained with an example:
208\medskip
209
210\noindent Consider an ORCA1
211configuration using the extended grid domain configuration file: \ifile{eORCA1\_domcfg.nc}
212This file define a horizontal domain of 362x332.  The first row with
213open ocean wet points in the non-isf bathymetry for this set is row 42 (\fortran indexing)
214then the formally correct setting for \np{open\_ocean\_jstart} is 41.  Using this value as
215the first row to be read will result in a 362x292 domain which is the same size as the
216original ORCA1 domain.  Thus the extended domain configuration file can be used with all
217the original input files for ORCA1 if the ice cavities are not active (\np{ln\_isfcav =
218.false.}).  Full instructions for achieving this are:
219
220\begin{itemize}
221\item  Add the new attribute to any input files requiring a j-row offset, i.e:
222\begin{cmds} 
223ncatted  -a open_ocean_jstart,global,a,d,41 eORCA1_domcfg.nc
224\end{cmds}
225
226\item Add the logical switch \np{ln\_use\_jattr} to \ngn{namcfg} in the configuration
227namelist (if it is not already there) and set \np{.true.}
228\end{itemize}
229
230\noindent Note that with this option, the j-size of the global domain is (extended
231j-size minus \np{open\_ocean\_jstart} + 1 ) and this must match the \texttt{jpjglo} value
232for the configuration. This means an alternative version of \ifile{eORCA1\_domcfg.nc} must
233be created for when \np{ln\_use\_jattr} is active. The \texttt{ncap2} tool provides a
234convenient way of achieving this:
235
236\begin{cmds}
237ncap2 -s 'jpjglo=292' eORCA1_domcfg.nc nORCA1_domcfg.nc
238\end{cmds}
239
240The domain configuration file is unique in this respect since it also contains the value
241of \texttt{jpjglo} that is read and used by the model. Any other global, 2D and 3D,
242netcdf, input field can be prepared for use in a reduced domain by adding the
243\texttt{open\_ocean\_jstart} attribute to the file's  global attributes. In particular
244this is true for any field that is read by \NEMO using the following optional argument to
245the appropriate call to \np{iom\_get}.
246\begin{forlines}
247lrowattr=ln_use_jattr
248\end{forlines}
249
250Currently, only the domain configuration variables make use of this optional argument so
251this facility is of little practical use except for tests where no other external input
252files are needed or you wish to use an extended domain configuration with inputs from
253earlier, non-extended configurations. Alternatively, it should be possible to exclude
254empty rows for extended domain, forced ocean runs using interpolation on the fly, by
255adding the optional argument to \texttt{iom\_get} calls for the weights and initial
256conditions. Experimenting with this remains an exercise for the user.
257
258% ================================================================
259% Accuracy and Reproducibility
260% ================================================================
261\section[Accuracy and reproducibility (\textit{lib\_fortran.F90})]
262{Accuracy and reproducibility (\protect\mdl{lib\_fortran})}
263\label{sec:MISC_fortran}
264
265\subsection[Issues with intrinsinc SIGN function (\texttt{\textbf{key\_nosignedzero}})]
266{Issues with intrinsinc SIGN function (\protect\key{nosignedzero})}
267\label{subsec:MISC_sign}
268
269The SIGN(A, B) is the \fortran intrinsic function delivers the magnitude of A with the sign of B.
270For example, SIGN(-3.0,2.0) has the value 3.0.
271The problematic case is when the second argument is zero, because, on platforms that support IEEE arithmetic,
272zero is actually a signed number.
273There is a positive zero and a negative zero.
274
275In \fninety, the processor was required always to deliver a positive result for SIGN(A, B) if B was zero.
276Nevertheless, in \fninety, the processor is allowed to do the correct thing and deliver ABS(A) when
277B is a positive zero and -ABS(A) when B is a negative zero.
278This change in the specification becomes apparent only when B is of type real, and is zero,
279and the processor is capable of distinguishing between positive and negative zero,
280and B is negative real zero.
281Then SIGN delivers a negative result where, under \fninety rules, it used to return a positive result.
282This change may be especially sensitive for the ice model,
283so we overwrite the intrinsinc function with our own function simply performing :   \\
284\verb?   IF( B >= 0.e0 ) THEN   ;   SIGN(A,B) = ABS(A)  ?    \\
285\verb?   ELSE                   ;   SIGN(A,B) =-ABS(A)     ?  \\
286\verb?   ENDIF    ? \\
287This feature can be found in \mdl{lib\_fortran} module and is effective when \key{nosignedzero} is defined.
288We use a CPP key as the overwritting of a intrinsic function can present performance issues with
289some computers/compilers.
290
291
292\subsection{MPP reproducibility}
293\label{subsec:MISC_glosum}
294
295The numerical reproducibility of simulations on distributed memory parallel computers is a critical issue.
296In particular, within NEMO global summation of distributed arrays is most susceptible to rounding errors,
297and their propagation and accumulation cause uncertainty in final simulation reproducibility on
298different numbers of processors.
299To avoid so, based on \citet{he.ding_JS01} review of different technics,
300we use a so called self-compensated summation method.
301The idea is to estimate the roundoff error, store it in a buffer, and then add it back in the next addition.
302
303Suppose we need to calculate $b = a_1 + a_2 + a_3$.
304The following algorithm will allow to split the sum in two
305($sum_1 = a_{1} + a_{2}$ and $b = sum_2 = sum_1 + a_3$) with exactly the same rounding errors as
306the sum performed all at once.
307\begin{align*}
308   sum_1 \ \  &= a_1 + a_2 \\
309   error_1     &= a_2 + ( a_1 - sum_1 ) \\
310   sum_2 \ \  &= sum_1 + a_3 + error_1 \\
311   error_2     &= a_3 + error_1 + ( sum_1 - sum_2 ) \\
312   b \qquad \ &= sum_2 \\
313\end{align*}
314An example of this feature can be found in \mdl{lib\_fortran} module.
315It is systematicallt used in glob\_sum function (summation over the entire basin excluding duplicated rows and
316columns due to cyclic or north fold boundary condition as well as overlap MPP areas).
317The self-compensated summation method should be used in all summation in i- and/or j-direction.
318See \mdl{closea} module for an example.
319Note also that this implementation may be sensitive to the optimization level.
320
321\subsection{MPP scalability}
322\label{subsec:MISC_mppsca}
323
324The default method of communicating values across the north-fold in distributed memory applications (\key{mpp\_mpi})
325uses a \textsc{MPI\_ALLGATHER} function to exchange values from each processing region in
326the northern row with every other processing region in the northern row.
327This enables a global width array containing the top 4 rows to be collated on every northern row processor and then
328folded with a simple algorithm.
329Although conceptually simple, this "All to All" communication will hamper performance scalability for
330large numbers of northern row processors.
331From version 3.4 onwards an alternative method is available which only performs direct "Peer to Peer" communications
332between each processor and its immediate "neighbours" across the fold line.
333This is achieved by using the default \textsc{MPI\_ALLGATHER} method during initialisation to
334help identify the "active" neighbours.
335Stored lists of these neighbours are then used in all subsequent north-fold exchanges to
336restrict exchanges to those between associated regions.
337The collated global width array for each region is thus only partially filled but is guaranteed to
338be set at all the locations actually required by each individual for the fold operation.
339This alternative method should give identical results to the default \textsc{ALLGATHER} method and
340is recommended for large values of \np{jpni}.
341The new method is activated by setting \np{ln\_nnogather} to be true (\ngn{nammpp}).
342The reproducibility of results using the two methods should be confirmed for each new,
343non-reference configuration.
344
345% ================================================================
346% Model optimisation, Control Print and Benchmark
347% ================================================================
348\section{Model optimisation, control print and benchmark}
349\label{sec:MISC_opt}
350%--------------------------------------------namctl-------------------------------------------------------
351
352\nlst{namctl} 
353%--------------------------------------------------------------------------------------------------------------
354
355Options are defined through the  \ngn{namctl} namelist variables.
356
357\subsection{Vector optimisation}
358
359\key{vectopt\_loop} enables the internal loops to collapse.
360This is very a very efficient way to increase the length of vector calculations and thus
361to speed up the model on vector computers.
362 
363% Add here also one word on NPROMA technique that has been found useless, since compiler have made significant progress during the last decade.
364 
365% Add also one word on NEC specific optimisation (Novercheck option for example)
366 
367\subsection{Control print}
368
369The \np{ln\_ctl} switch was originally used as a debugging option in two modes:
370
371\begin{enumerate}
372\item{\np{ln\_ctl}: compute and print the trends averaged over the interior domain in all TRA, DYN, LDF and
373ZDF modules.
374This option is very helpful when diagnosing the origin of an undesired change in model results. }
375
376\item{also \np{ln\_ctl} but using the nictl and njctl namelist parameters to check the source of differences between
377mono and multi processor runs.}
378\end{enumerate}
379
380However, in recent versions it has also been used to force all processors to assume the
381reporting role. Thus when \np{ln\_ctl} is true all processors produce their own versions
382of files such as: ocean.output, layout.dat, etc.  All such files, beyond the the normal
383reporting processor (narea == 1), are named with a \_XXXX extension to their name, where
384XXXX is a 4-digit area number (with leading zeros, if required). Other reporting files
385such as run.stat (and its netCDF counterpart: run.stat.nc) and tracer.stat contain global
386information and are only ever produced by the reporting master (narea == 1). For version
3874.0 a start has been made to return \np{ln\_ctl} to its original function by introducing
388a new control structure which allows finer control over which files are produced. This
389feature is still evolving but it does already allow the user to: select individually the
390production of run.stat and tracer.stat files and to toggle the production of other files
391on processors other than the reporting master. These other reporters can be a simple
392subset of processors as defined by a minimum, maximum and incremental processor number.
393
394Note, that production of the run.stat and tracer.stat files require global communications.
395For run.stat, these are global min and max operations to find metrics such as the gloabl
396maximum velocity. For tracer.stat these are global sums of tracer fields. To improve model
397performance these operations are disabled by default and, where necessary, any use of the
398global values have been replaced with local calculations. For example, checks on the CFL
399criterion are now done on the local domain and only reported if a breach is detected.
400
401Experienced users may wish to still monitor this information as a check on model progress.
402If so, the best compromise will be to activate the files with:
403
404\begin{verbatim}
405     sn_cfctl%l_config = .TRUE.
406       sn_cfctl%l_runstat = .TRUE.
407       sn_cfctl%l_trcstat = .TRUE.
408\end{verbatim}
409
410and to use the new time increment setting to ensure the values are collected and reported
411at a suitably long interval. For example:
412
413\begin{verbatim}     
414       sn_cfctl%ptimincr  = 25
415\end{verbatim}
416
417will carry out the global communications and write the information every 25 timesteps. This
418increment also applies to the time.step file which is otherwise updated every timestep.
419
420% ================================================================
421\biblio
422
423\pindex
424
425\end{document}
Note: See TracBrowser for help on using the repository browser.