source: NEMO/branches/2020/ticket_2444/doc/latex/NEMO/subfiles/chap_misc.tex @ 12769

Last change on this file since 12769 was 12769, checked in by mathiot, 11 months ago

#2444: update of the documentation (prior to changes suggested by Dave)

File size: 27.6 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   X.X} & {\em Pierre Mathiot} & {update of the closed sea section}
19    {\em   4.0} & {\em ...} & {\em ...} \\
20    {\em   3.6} & {\em ...} & {\em ...} \\
21    {\em   3.4} & {\em ...} & {\em ...} \\
22    {\em <=3.4} & {\em ...} & {\em ...}
23  \end{tabularx}
28%% =================================================================================================
29\section{Representation of unresolved straits}
32In climate modeling, it often occurs that a crucial connections between water masses is broken as
33the grid mesh is too coarse to resolve narrow straits.
34For example, coarse grid spacing typically closes off the Mediterranean from the Atlantic at
35the Strait of Gibraltar.
36In this case, it is important for climate models to include the effects of salty water entering the Atlantic from
37the Mediterranean.
38Likewise, it is important for the Mediterranean to replenish its supply of water from the Atlantic to
39balance the net evaporation occurring over the Mediterranean region.
40This problem occurs even in eddy permitting simulations.
41For example, in ORCA 1/4\deg\ several straits of the Indonesian archipelago (Ombai, Lombok...)
42are much narrow than even a single ocean grid-point.
44We describe briefly here the two methods that can be used in \NEMO\ to handle such
45improperly resolved straits. The methods consist of opening the strait while ensuring
46that the mass exchanges through the strait are not too large by either artificially
47reducing the cross-sectional area of the strait grid-cells or, locally increasing the
48lateral friction.
50%% =================================================================================================
51\subsection{Hand made geometry changes}
54The first method involves reducing the scale factor in the cross-strait direction to a
55value in better agreement with the true mean width of the strait
56(\autoref{fig:MISC_strait_hand}).  This technique is sometime called "partially open face"
57or "partially closed cells".  The key issue here is only to reduce the faces of $T$-cell
58(\ie\ change the value of the horizontal scale factors at $u$- or $v$-point) but not the
59volume of the $T$-cell.  Indeed, reducing the volume of strait $T$-cell can easily produce
60a numerical instability at that grid point which would require a reduction of the model
61time step.  Thus to instigate a local change in the width of a Strait requires two steps:
65\item Add \texttt{e1e2u} and \texttt{e1e2v} arrays to the \np{cn_domcfg}{cn\_domcfg} file. These 2D
66arrays should contain the products of the unaltered values of: $\texttt{e1u}*\texttt{e2u}$
67and $\texttt{e1u}*\texttt{e2v}$ respectively. That is the original surface areas of $u$-
68and $v$- cells respectively.  These areas are usually defined by the corresponding product
69within the \NEMO\ code but the presence of \texttt{e1e2u} and \texttt{e1e2v} in the
70\np{cn_domcfg}{cn\_domcfg} file will suppress this calculation and use the supplied fields instead.
71If the model domain is provided by user-supplied code in \mdl{usrdef\_hgr}, then this
72routine should also return \texttt{e1e2u} and \texttt{e1e2v} and set the integer return
73argument \texttt{ie1e2u\_v} to a non-zero value. Values other than 0 for this argument
74will suppress the calculation of the areas.
76\item Change values of \texttt{e2u} or \texttt{e1v} (either in the \np{cn_domcfg}{cn\_domcfg} file or
77via code in  \mdl{usrdef\_hgr}), whereever a Strait reduction is required. The choice of
78whether to alter \texttt{e2u} or \texttt{e1v} depends. respectively,  on whether the
79Strait in question is North-South orientated (\eg\ Gibraltar) or East-West orientated (\eg
84The second method is to increase the viscous boundary layer thickness by a local increase
85of the fmask value at the coast. This method can also be effective in wider passages.  The
86concept is illustarted in the second part of  \autoref{fig:MISC_strait_hand} and changes
87to specific locations can be coded in \mdl{usrdef\_fmask}. The \forcode{usr_def_fmask}
88routine is always called after \texttt{fmask} has been defined according to the choice of
89lateral boundary condition as discussed in \autoref{sec:LBC_coast}. The default version of
90\mdl{usrdef\_fmask} contains settings specific to ORCA2 and ORCA1 configurations. These are
91meant as examples only; it is up to the user to verify settings and provide alternatives
92for their own configurations. The default \forcode{usr_def_fmask} makes no changes to
93\texttt{fmask} for any other configuration.
96  \centering
97  \includegraphics[width=0.66\textwidth]{MISC_Gibraltar}
98  \includegraphics[width=0.66\textwidth]{MISC_Gibraltar2}
99  \caption[Two methods to defined the Gibraltar strait]{
100    Example of the Gibraltar strait defined in a 1\deg\ $\times$ 1\deg\ mesh.
101    \textit{Top}: using partially open cells.
102    The meridional scale factor at $v$-point is reduced on both sides of the strait to
103    account for the real width of the strait (about 20 km).
104    Note that the scale factors of the strait $T$-point remains unchanged.
105    \textit{Bottom}: using viscous boundary layers.
106    The four fmask parameters along the strait coastlines are set to a value larger than 4,
107    \ie\ "strong" no-slip case (see \autoref{fig:LBC_shlat}) creating a large viscous boundary layer
108    that allows a reduced transport through the strait.}
109  \label{fig:MISC_strait_hand}
112%% =================================================================================================
113\section[Closed seas (\textit{closea.F90})]{Closed seas (\protect\mdl{closea})}
117  \nlst{namclo}
118  \caption{\forcode{&namclo}}
119  \label{lst:namclo}
122Some configurations include inland seas and lakes as ocean
123points. This is particularly the case for configurations that are
124coupled to an atmosphere model where one might want to include inland
125seas and lakes as ocean model points in order to provide a better
126bottom boundary condition for the atmosphere. However there is no
127route for freshwater to run off from the lakes to the ocean and this
128can lead to large drifts in the sea surface height over the lakes. The
129closea module provides options to either fill in closed seas and lakes
130at run time, or to set the net surface freshwater flux for each lake
131to zero and put the residual flux into the ocean.
133The inland seas and lakes are defined using mask fields in the
134domain configuration file. Special treatment of the closed sea (redistribution of net freshwater or mask those), are defined in \autoref{lst:namclo} and
135can be trigger by \np{ln_closea}{ln\_closea}\forcode{=.true.} in namelist namcfg.
137The options available are the following:
139\item[\np{ln_maskcs}{ln\_maskcs}\forcode{ = .true.}] All the closed seas are masked using \textit{mask\_opensea} variable.
140\item[\np{ln_maskcs}{ln\_maskcs}\forcode{ = .false.}] The net surface flux over each inland sea or group of
141inland seas is set to zero each timestep and the residual flux is
142distributed over a target area.
145When \np{ln_maskcs}{ln\_maskcs}\forcode{ = .false.},
1463 options are available for the redistribution (set up of these options is done in the tool DOMAINcfg):
147\begin{description}[font=$\bullet$ ]
148\item[ glo]: The residual flux is redistributed globally.
149\item[ emp]: The residual flux is redistributed as emp in a river outflow.
150\item[ rnf]: The residual flux is redistributed as rnf in a river outflow if negative. If there is a net evaporation, the residual flux is redistributed globally.
153For each case, 2 masks are needed (\autoref{fig:MISC_closea_mask_example}):
155\item $\bullet$ one describing the 'sources' (ie the closed seas concerned by each options) called \textit{mask\_csglo}, \textit{mask\_csrnf}, \textit{mask\_csemp}.
156\item $\bullet$ one describing each group of inland seas (the Great Lakes for example) and the target area (river outflow or world ocean) for each group of inland seas (St Laurence for the Great Lakes for example) called
157\textit{mask\_csgrpglo}, \textit{mask\_csgrprnf}, \textit{mask\_csgrpemp}.
161  \centering
162  \includegraphics[width=0.66\textwidth]{MISC_closea_mask_example}
163  \caption[Mask fields for the \protect\mdl{closea} module]{
164    Example of mask fields for the \protect\mdl{closea} module.
165    \textit{Left}: a \textit{mask\_csrnf} field;
166    \textit{Right}: a \textit{mask\_csgrprnf} field.
167    In this example, if \protect\np{ln_closea}{ln\_closea} is set to \forcode{.true.},
168    the mean freshwater flux over each of the American Great Lakes will be set to zero,
169    and the total residual for all the lakes, if negative, will be put into
170    the St Laurence Seaway in the area shown.}
171  \label{fig:MISC_closea_mask_example}
174Closed sea not defined (because too small, issue in the bathymetry definition ...) are defined in \textit{mask\_csundef}.
175These points can be masked using the namelist option \np{ln_mask_csundef}{ln\_mask\_csundef}\forcode{= .true.} or used to correct the bathymetry input file.\\
177The masks needed for the closed sea can be created using the DOMAINcfg tool in the utils/tools/DOMAINcfg directory.
178See \autoref{sec:clocfg} for details on the usage of definition of the closed sea masks.
180%% =================================================================================================
181\section{Sub-domain functionality}
184%% =================================================================================================
185\subsection{Simple subsetting of input files via NetCDF attributes}
187The extended grids for use with the under-shelf ice cavities will result in redundant rows
188around Antarctica if the ice cavities are not active.  A simple mechanism for subsetting
189input files associated with the extended domains has been implemented to avoid the need to
190maintain different sets of input fields for use with or without active ice cavities.  This
191subsetting operates for the j-direction only and works by optionally looking for and using
192a global file attribute (named: \np{open_ocean_jstart}{open\_ocean\_jstart}) to determine the starting j-row
193for input.  The use of this option is best explained with an example:
196\noindent Consider an ORCA1
197configuration using the extended grid domain configuration file: \ifile{eORCA1\}
198This file define a horizontal domain of 362x332.  The first row with
199open ocean wet points in the non-isf bathymetry for this set is row 42 (\fortran\ indexing)
200then the formally correct setting for \np{open_ocean_jstart}{open\_ocean\_jstart} is 41.  Using this value as
201the first row to be read will result in a 362x292 domain which is the same size as the
202original ORCA1 domain.  Thus the extended domain configuration file can be used with all
203the original input files for ORCA1 if the ice cavities are not active (\np{ln\_isfcav =
204.false.}).  Full instructions for achieving this are:
207\item Add the new attribute to any input files requiring a j-row offset, i.e:
209ncatted  -a open_ocean_jstart,global,a,d,41
212\item Add the logical switch \np{ln_use_jattr}{ln\_use\_jattr} to \nam{cfg}{cfg} in the configuration
213namelist (if it is not already there) and set \forcode{.true.}
216\noindent Note that with this option, the j-size of the global domain is (extended
217j-size minus \np{open_ocean_jstart}{open\_ocean\_jstart} + 1 ) and this must match the \texttt{jpjglo} value
218for the configuration. This means an alternative version of \ifile{eORCA1\} must
219be created for when \np{ln_use_jattr}{ln\_use\_jattr} is active. The \texttt{ncap2} tool provides a
220convenient way of achieving this:
223ncap2 -s 'jpjglo=292'
226The domain configuration file is unique in this respect since it also contains the value of \jp{jpjglo}
227that is read and used by the model.
228Any other global, 2D and 3D, netcdf, input field can be prepared for use in a reduced domain by adding the
229\texttt{open\_ocean\_jstart} attribute to the file's global attributes.
230In particular this is true for any field that is read by \NEMO\ using the following optional argument to
231the appropriate call to \np{iom_get}{iom\_get}.
237Currently, only the domain configuration variables make use of this optional argument so
238this facility is of little practical use except for tests where no other external input
239files are needed or you wish to use an extended domain configuration with inputs from
240earlier, non-extended configurations. Alternatively, it should be possible to exclude
241empty rows for extended domain, forced ocean runs using interpolation on the fly, by
242adding the optional argument to \texttt{iom\_get} calls for the weights and initial
243conditions. Experimenting with this remains an exercise for the user.
245%% =================================================================================================
246\section[Accuracy and reproducibility (\textit{lib\_fortran.F90})]{Accuracy and reproducibility (\protect\mdl{lib\_fortran})}
249%% =================================================================================================
250\subsection[Issues with intrinsinc SIGN function (\texttt{\textbf{key\_nosignedzero}})]{Issues with intrinsinc SIGN function (\protect\key{nosignedzero})}
253The SIGN(A, B) is the \fortran\ intrinsic function delivers the magnitude of A with the sign of B.
254For example, SIGN(-3.0,2.0) has the value 3.0.
255The problematic case is when the second argument is zero, because, on platforms that support IEEE arithmetic,
256zero is actually a signed number.
257There is a positive zero and a negative zero.
259In \fninety, the processor was required always to deliver a positive result for SIGN(A, B) if B was zero.
260Nevertheless, in \fninety, the processor is allowed to do the correct thing and deliver ABS(A) when
261B is a positive zero and -ABS(A) when B is a negative zero.
262This change in the specification becomes apparent only when B is of type real, and is zero,
263and the processor is capable of distinguishing between positive and negative zero,
264and B is negative real zero.
265Then SIGN delivers a negative result where, under \fninety\ rules, it used to return a positive result.
266This change may be especially sensitive for the ice model,
267so we overwrite the intrinsinc function with our own function simply performing :   \\
268\verb?   IF( B >= 0.e0 ) THEN   ;   SIGN(A,B) = ABS(A)  ?    \\
269\verb?   ELSE                   ;   SIGN(A,B) =-ABS(A)     ?  \\
270\verb?   ENDIF    ? \\
271This feature can be found in \mdl{lib\_fortran} module and is effective when \key{nosignedzero} is defined.
272We use a CPP key as the overwritting of a intrinsic function can present performance issues with
273some computers/compilers.
275%% =================================================================================================
276\subsection{MPP reproducibility}
279The numerical reproducibility of simulations on distributed memory parallel computers is a critical issue.
280In particular, within \NEMO\ global summation of distributed arrays is most susceptible to rounding errors,
281and their propagation and accumulation cause uncertainty in final simulation reproducibility on
282different numbers of processors.
283To avoid so, based on \citet{he.ding_JS01} review of different technics,
284we use a so called self-compensated summation method.
285The idea is to estimate the roundoff error, store it in a buffer, and then add it back in the next addition.
287Suppose we need to calculate $b = a_1 + a_2 + a_3$.
288The following algorithm will allow to split the sum in two
289($sum_1 = a_{1} + a_{2}$ and $b = sum_2 = sum_1 + a_3$) with exactly the same rounding errors as
290the sum performed all at once.
292   sum_1 \ \  &= a_1 + a_2 \\
293   error_1     &= a_2 + ( a_1 - sum_1 ) \\
294   sum_2 \ \  &= sum_1 + a_3 + error_1 \\
295   error_2     &= a_3 + error_1 + ( sum_1 - sum_2 ) \\
296   b \qquad \ &= sum_2 \\
298An example of this feature can be found in \mdl{lib\_fortran} module.
299It is systematicallt used in glob\_sum function (summation over the entire basin excluding duplicated rows and
300columns due to cyclic or north fold boundary condition as well as overlap MPP areas).
301The self-compensated summation method should be used in all summation in i- and/or j-direction.
302See \mdl{closea} module for an example.
303Note also that this implementation may be sensitive to the optimization level.
305%% =================================================================================================
306\subsection{MPP scalability}
309The default method of communicating values across the north-fold in distributed memory applications (\key{mpp\_mpi})
310uses a \textsc{MPI\_ALLGATHER} function to exchange values from each processing region in
311the northern row with every other processing region in the northern row.
312This enables a global width array containing the top 4 rows to be collated on every northern row processor and then
313folded with a simple algorithm.
314Although conceptually simple, this "All to All" communication will hamper performance scalability for
315large numbers of northern row processors.
316From version 3.4 onwards an alternative method is available which only performs direct "Peer to Peer" communications
317between each processor and its immediate "neighbours" across the fold line.
318This is achieved by using the default \textsc{MPI\_ALLGATHER} method during initialisation to
319help identify the "active" neighbours.
320Stored lists of these neighbours are then used in all subsequent north-fold exchanges to
321restrict exchanges to those between associated regions.
322The collated global width array for each region is thus only partially filled but is guaranteed to
323be set at all the locations actually required by each individual for the fold operation.
324This alternative method should give identical results to the default \textsc{ALLGATHER} method and
325is recommended for large values of \np{jpni}{jpni}.
326The new method is activated by setting \np{ln_nnogather}{ln\_nnogather} to be true (\nam{mpp}{mpp}).
327The reproducibility of results using the two methods should be confirmed for each new,
328non-reference configuration.
330%% =================================================================================================
331\section{Model optimisation, control print and benchmark}
335  \nlst{namctl}
336  \caption{\forcode{&namctl}}
337  \label{lst:namctl}
340Options are defined through the  \nam{ctl}{ctl} namelist variables.
342%% =================================================================================================
343\subsection{Vector optimisation}
345\key{vectopt\_loop} enables the internal loops to collapse.
346This is very a very efficient way to increase the length of vector calculations and thus
347to speed up the model on vector computers.
349% Add here also one word on NPROMA technique that has been found useless, since compiler have made significant progress during the last decade.
351% Add also one word on NEC specific optimisation (Novercheck option for example)
353%% =================================================================================================
354\subsection{Status and debugging information output}
357NEMO can produce a range of text information output either: in the main output
358file (ocean.output) written by the normal reporting processor (narea == 1) or various
359specialist output files (e.g. layout.dat, run.stat, tracer.stat etc.). Some, for example
360run.stat and tracer.stat, contain globally collected values for which a single file is
361sufficient. Others, however, contain information that could, potentially, be different
362for each processing region. For computational efficiency, the default volume of text
363information produced is reduced to just a few files from the narea=1 processor.
365When more information is required for monitoring or debugging purposes, the various
366forms of output can be selected via the \np{sn\_cfctl} structure. As well as simple
367on-off switches this structure also allows selection of a range of processors for
368individual reporting (where appropriate) and a time-increment option to restrict
369globally collected values to specified time-step increments.
371Most options within the structure are influenced by the top-level switches shown here
372with their default settings:
375   sn_cfctl%l_allon  = .FALSE.    ! IF T activate all options. If F deactivate all unless l_config is T
376     sn_cfctl%l_config = .TRUE.     ! IF .true. then control which reports are written with the following
379The first switch is a convenience option which can be used to switch on and off all
380sub-options. However, if it is false then switching off all sub-options is only done
381if \texttt{sn_cfctl%l\_config} is also false. Specifically, the logic is:
384  IF ( sn_cfctl%l_allon ) THEN
385    set all suboptions .TRUE.
386    and set procmin, procmax and procincr so that all regions are selected ([0,10000000,1], respectively)
387  ELSEIF ( sn_cfctl%l_config ) THEN
388    honour individual settings of the suboptions from the namelist
389  ELSE
390    set all suboptions .FALSE.
391  ENDIF
394Details of the suboptions follow but first an explanation of the stand-alone option:
395\texttt{sn_cfctl%l_glochk}.  This option modifies the action of the early warning checks
396carried out in \textt{stpctl.F90}. These checks detect probable numerical instabilites
397by searching for excessive sea surface heights or velocities and salinity values
398outside a sensible physical range. If breaches are detected then the default behaviour
399is to locate and report the local indices of the grid-point in breach. These indices
400are included in the error message that precedes the model shutdown. When true,
401\texttt{sn_cfctl%l_glochk} modifies this action by performing a global location of
402the various minimum and maximum values and the global indices are reported. This has
403some value in locating the most severe error in cases where the first detected error
404may not be the worst culprit.
406\subsubsection{Control print suboptions}
408The options that can be individually selected fall into three categories:
410\begin{enumerate} \item{Time step progress information} This category includes
411\texttt{run.stat} and \texttt{tracer.stat} files which record certain physical and
412passive tracer metrics (respectively). Typical contents of \texttt{run.stat} include
413global maximums of ssh, velocity; and global minimums and maximums of temperature
414and salinity.  A netCDF version of \texttt{run.stat} (\texttt{}) is also
415produced with the same time-series data and this can easily be expanded to include
416extra monitoring information.  \texttt{tracer.stat} contains the volume-weighted
417average tracer value for each passive tracer. Collecting these metrics involves
418global communications and will impact on model efficiency so both these options are
419disabled by default by setting the respective options, \texttt{sn\_cfctl%runstat} and
420\texttt{sn\_cfctl%trcstat} to false. A compromise can be made by activating either or
421both of these options and setting the \texttt{sn\_cfctl%timincr} entry to an integer
422value greater than one. This increment determines the time-step frequency at which
423the global metrics are collected and reported.  This increment also applies to the
424time.step file which is otherwise updated every timestep.
425\item{One-time configuration information/progress logs}
427Some run-time configuration information and limited progress information is always
428produced by the first ocean process. This includes the \texttt{ocean.output} file
429which reports on all the namelist options read by the model and remains open to catch
430any warning or error messages generated during execution. A \texttt{layout.dat}
431file is also produced which details the MPI-decomposition used by the model. The
432suboptions: \texttt{sn\_cfctl%oceout} and \texttt{sn\_cfctl%layout} can be used
433to activate the creation of these files by all ocean processes.  For example,
434when \texttt{sn\_cfctl%oceout} is true all processors produce their own version of
435\texttt{ocean.output}.  All files, beyond the the normal reporting processor (narea == 1), are
436named with a \_XXXX extension to their name, where XXXX is a 4-digit area number (with
437leading zeros, if required). This is useful as a debugging aid since all processes can
438report their local conditions. Note though that these files are buffered on most UNIX
439systems so bug-hunting efforts using this facility should also utilise the \fortran:
442   CALL FLUSH(numout)
445statement after any additional write statements to ensure that file contents reflect
446the last model state. Associated with the \texttt{sn\_cfctl%oceout} option is the
447additional \texttt{sn\_cfctl%oasout} suboption. This does not activate its own output
448file but rather activates the writing of addition information regarding the OASIS
449configuration when coupling via oasis and the sbccpl routine. This information is
450written to any active \texttt{ocean.output} files.
451\item{Control sums of trends for debugging}
453NEMO includes an option for debugging reproducibility differences between
454a MPP and mono-processor runs.  This is somewhat dated and clearly only
455useful for this purpose when dealing with configurations that can be run
456on a single processor. The full details can be found in this report: \href{
458control print option in NEMO} The switches to activate production of the control sums
459of trends for either the physics or passive tracers are the \texttt{sn\_cfctl%prtctl}
460and \texttt{sn\_cfctl%prttrc} suboptions, respectively. Although, perhaps, of limited use for its
461original intention, the ability to produce these control sums of trends in specific
462areas provides another tool for diagnosing model behaviour.  If only the output from a
463select few regions is required then additional options are available to activate options
464for only a simple subset of processing regions. These are: \texttt{sn\_cfctl%procmin},
465\texttt{sn\_cfctl%procmax} and \texttt{sn\_cfctl%procincr} which can be used to specify
466the minimum and maximum active areas and the increment. The default values are set
467such that all regions will be active. Note this subsetting can also be used to limit
468which additional \texttt{ocean.output} and \texttt{layout.dat} files are produced if
469those suboptions are active.
474   sn_cfctl%l_glochk = .FALSE.    ! Range sanity checks are local (F) or global (T). Set T for debugging only
475   sn_cfctl%l_allon  = .FALSE.    ! IF T activate all options. If F deactivate all unless l_config is T
476     sn_cfctl%l_config = .TRUE.     ! IF .true. then control which reports are written with the following
477       sn_cfctl%l_runstat = .FALSE. ! switches and which areas produce reports with the proc integer settings.
478       sn_cfctl%l_trcstat = .FALSE. ! The default settings for the proc integers should ensure
479       sn_cfctl%l_oceout  = .FALSE. ! that  all areas report.
480       sn_cfctl%l_layout  = .FALSE. !
481       sn_cfctl%l_prtctl  = .FALSE. !
482       sn_cfctl%l_prttrc  = .FALSE. !
483       sn_cfctl%l_oasout  = .FALSE. !
484       sn_cfctl%procmin   = 0       ! Minimum area number for reporting [default:0]
485       sn_cfctl%procmax   = 1000000 ! Maximum area number for reporting [default:1000000]
486       sn_cfctl%procincr  = 1       ! Increment for optional subsetting of areas [default:1]
487       sn_cfctl%ptimincr  = 1       ! Timestep increment for writing time step progress info
Note: See TracBrowser for help on using the repository browser.