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

Last change on this file since 11693 was 11693, checked in by nicolasmartin, 13 months ago

Macros renaming

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