Ignore:
Timestamp:
2019-10-29T18:14:49+01:00 (13 months ago)
Message:

Update the branch to r11830 of the trunk!

Location:
Files:
36 deleted
82 edited
92 copied

### Legend:

Unmodified
Removed

• Property svn:ignore deleted
• Property svn:externals set to
^/utils/logos logos

• Property svn:ignore deleted

 r11004 helvetic hyperref idxlayout ifplatform imakeidx koma-script latex latexconfig latex-bin latexconfig latex-fonts latexmk oberdiek parskip preview psnfss subfiles times titlesec titling tools upquote wrapfig xcolor xkeyval xstring zapfchan

 r10146 #!/bin/bash ./inc/clean.sh ./inc/build.sh #./inc/clean.sh #./inc/build.sh cd tex_main sed -i -e 's#\\documentclass#%\\documentclass#' -e '/{document}/ s/^/%/' \ ../tex_sub/*.tex sed -i    's#\\subfile{#\\include{#g' \ NEMO_manual.tex latex2html -local_icons -no_footnode -split 4 -link 2 -mkdir -dir ../html_LaTeX2HTML   \ $* \ NEMO_manual sed -i -e 's#%\\documentclass#\\documentclass#' -e '/{document}/ s/^%//' \ ../tex_sub/*.tex sed -i 's#\\include{#\\subfile{#g' \ NEMO_manual.tex sed -i -e 's#utf8#latin1#' \ -e 's#$outputdir=../build${minted}#{minted}#' \ -e '/graphicspath/ s#{../#{../../#g' \ global/packages.tex cd ./NEMO/main sed -i -e 's#\\documentclass#%\\documentclass#' -e '/{document}/ s#^#%#' ../subfiles/*.tex sed -i 's#\\subfile{#\\input{#' chapters.tex appendices.tex #latex2html -noimages -local_icons -no_footnode -split 4 -link 2 -dir ../html_LaTeX2HTML$* NEMO_manual latex2html -debug -noreuse -init_file ../../l2hconf.pm -local_icons                               -dir ../build/html               NEMO_manual sed -i -e 's#%\\documentclass#\\documentclass#' -e '/{document}/ s#^%##'   ../subfiles/*.tex sed -i    's#\\input{#\\subfile{#'                                         chapters.tex appendices.tex cd - sed -i -e 's#latin1#utf8#'                                 \ -e 's#{minted}#$outputdir=../build${minted}#' \ -e '/graphicspath/ s#{../../#{../#g'                \ global/packages.tex exit 0

• Property svn:ignore deleted
• Property svn:externals set to
^/utils/figures/NEMO figures

• Property svn:ignore
•  old *.aux *.bbl *.blg *.dvi *.fdb* *.fls *.idx *.ilg *.ind *.log *.maf *.mtc* *.out *.pdf *.toc _minted-*

• Property svn:ignore
•  old *.aux *.bbl *.blg *.dvi *.fdb* *.fls *.idx *.ind *.ilg *.ind *.log *.maf *.mtc* *.out *.pdf *.toc _minted-*

 r10442 \begin{document} % ================================================================ % Chapter Assimilation increments (ASM) % ================================================================ \chapter{Apply Assimilation Increments (ASM)} \label{chap:ASM} Authors: D. Lea,  M. Martin, K. Mogensen, A. Weaver, ...   % do we keep %    {\em 4.0} & {\em D. J. Lea} & {\em \NEMO\ 4.0 updates}  \\ %    {\em 3.4} & {\em D. J. Lea, M. Martin, K. Mogensen, A. Weaver} & {\em Initial version}  \\ \minitoc \thispagestyle{plain} \newpage \chaptertoc \paragraph{Changes record} ~\\ {\footnotesize \begin{tabularx}{\textwidth}{l||X|X} Release & Author(s) & Modifications \\ \hline {\em   4.0} & {\em ...} & {\em ...} \\ {\em   3.6} & {\em ...} & {\em ...} \\ {\em   3.4} & {\em ...} & {\em ...} \\ {\em <=3.4} & {\em ...} & {\em ...} \end{tabularx} } \clearpage The ASM code adds the functionality to apply increments to the model variables: temperature, salinity, sea surface height, velocity and sea ice concentration. sea surface height, velocity and sea ice concentration. These are read into the model from a NetCDF file which may be produced by separate data assimilation code. The code can also output model background fields which are used as an input to data assimilation code. This is all controlled by the namelist \textit{\ngn{nam\_asminc} }. This is all controlled by the namelist \nam{_asminc}{\_asminc}. There is a brief description of all the namelist options provided. To build the ASM code \key{asminc} must be set. %=============================================================== %% ================================================================================================= \section{Direct initialization} \label{sec:ASM_DI} Direct initialization (DI) refers to the instantaneous correction of the model background state using the analysis increment. DI is used when \np{ln\_asmdin} is set to true. DI is used when \np{ln_asmdin}{ln\_asmdin} is set to true. %% ================================================================================================= \section{Incremental analysis updates} \label{sec:ASM_IAU} it may be preferable to introduce the increment gradually into the ocean model in order to minimize spurious adjustment processes. This technique is referred to as Incremental Analysis Updates (IAU) \citep{Bloom_al_MWR96}. This technique is referred to as Incremental Analysis Updates (IAU) \citep{bloom.takacs.ea_MWR96}. IAU is a common technique used with 3D assimilation methods such as 3D-Var or OI. IAU is used when \np{ln\_asmiau} is set to true. IAU is used when \np{ln_asmiau}{ln\_asmiau} is set to true. With IAU, the model state trajectory ${\bf x}$ in the assimilation window ($t_{0} \leq t_{i} \leq t_{N}$) With IAU, the model state trajectory ${\mathbf x}$ in the assimilation window ($t_{0} \leq t_{i} \leq t_{N}$) is corrected by adding the analysis increments for temperature, salinity, horizontal velocity and SSH as additional tendency terms to the prognostic equations: \begin{align*} % \label{eq:wa_traj_iau} {\bf x}^{a}(t_{i}) = M(t_{i}, t_{0})[{\bf x}^{b}(t_{0})] \; + \; F_{i} \delta \tilde{\bf x}^{a} % \label{eq:ASM_wa_traj_iau} {\mathbf x}^{a}(t_{i}) = M(t_{i}, t_{0})[{\mathbf x}^{b}(t_{0})] \; + \; F_{i} \delta \tilde{\mathbf x}^{a} \end{align*} where $F_{i}$ is a weighting function for applying the increments $\delta\tilde{\bf x}^{a}$ defined such that where $F_{i}$ is a weighting function for applying the increments $\delta\tilde{\mathbf x}^{a}$ defined such that $\sum_{i=1}^{N} F_{i}=1$. ${\bf x}^b$ denotes the model initial state and ${\bf x}^a$ is the model state after the increments are applied. ${\mathbf x}^b$ denotes the model initial state and ${\mathbf x}^a$ is the model state after the increments are applied. To control the adjustment time of the model to the increment, the increment can be applied over an arbitrary sub-window, $t_{m} \leq t_{i} \leq t_{n}$, Typically the increments are spread evenly over the full window. In addition, two different weighting functions have been implemented. The first function employs constant weights, The first function (namelist option \np{niaufn}{niaufn}=0) employs constant weights, \begin{align} \label{eq:F1_i} \label{eq:ASM_F1_i} F^{(1)}_{i} =\left\{ \begin{array}{ll} 0     &    {\rm if} \; \; \; t_{i} < t_{m}                \\ 1/M &    {\rm if} \; \; \; t_{m} < t_{i} \leq t_{n} \\ 0     &    {\rm if} \; \; \; t_{i} > t_{n} 0     &    {\mathrm if} \; \; \; t_{i} < t_{m}                \\ 1/M &    {\mathrm if} \; \; \; t_{m} < t_{i} \leq t_{n} \\ 0     &    {\mathrm if} \; \; \; t_{i} > t_{n} \end{array} \right. \right. \end{align} where $M = m-n$. The second function employs peaked hat-like weights in order to give maximum weight in the centre of the sub-window, The second function (namelist option \np{niaufn}{niaufn}=1) employs peaked hat-like weights in order to give maximum weight in the centre of the sub-window, with the weighting reduced linearly to a small value at the window end-points: \begin{align} \label{eq:F2_i} \label{eq:ASM_F2_i} F^{(2)}_{i} =\left\{ \begin{array}{ll} 0                           &    {\rm if} \; \; \; t_{i}       <     t_{m}                        \\ \alpha \, i               &    {\rm if} \; \; \; t_{m}    \leq t_{i}    \leq   t_{M/2}   \\ \alpha \, (M - i +1) &    {\rm if} \; \; \; t_{M/2}  <    t_{i}    \leq   t_{n}       \\ 0                            &   {\rm if} \; \; \; t_{i}        >    t_{n} 0                           &    {\mathrm if} \; \; \; t_{i}       <     t_{m}                        \\ \alpha \, i               &    {\mathrm if} \; \; \; t_{m}    \leq t_{i}    \leq   t_{M/2}   \\ \alpha \, (M - i +1) &    {\mathrm if} \; \; \; t_{M/2}  <    t_{i}    \leq   t_{n}       \\ 0                            &   {\mathrm if} \; \; \; t_{i}        >    t_{n} \end{array} \right. \end{align} where $\alpha^{-1} = \sum_{i=1}^{M/2} 2i$ and $M$ is assumed to be even. The weights described by \autoref{eq:F2_i} provide a smoother transition of the analysis trajectory from one assimilation cycle to the next than that described by \autoref{eq:F1_i}. where $\alpha^{-1} = \sum_{i=1}^{M/2} 2i$ and $M$ is assumed to be even. The weights described by \autoref{eq:ASM_F2_i} provide a smoother transition of the analysis trajectory from one assimilation cycle to the next than that described by \autoref{eq:ASM_F1_i}. %========================================================================== % Divergence damping description %%% %% ================================================================================================= \section{Divergence damping initialisation} \label{sec:ASM_div_dmp} The velocity increments may be initialized by the iterative application of a divergence damping operator. In iteration step $n$ new estimates of velocity increments $u^{n}_I$ and $v^{n}_I$ are updated by: It is quite challenging for data assimilation systems to provide non-divergent velocity increments. Applying divergent velocity increments will likely cause spurious vertical velocities in the model. This section describes a method to take velocity increments provided to \NEMO\ ($u^0_I$ and $v^0_I$) and adjust them by the iterative application of a divergence damping operator. The method is also described in \citet{dobricic.pinardi.ea_OS07}. In iteration step $n$ (starting at $n=1$) new estimates of velocity increments $u^{n}_I$ and $v^{n}_I$ are updated by: \label{eq:asm_dmp} \label{eq:ASM_dmp} \left\{ \begin{aligned} \right., where where the divergence is defined as $% \label{eq:asm_div} % \label{eq:ASM_div} \chi^{n-1}_I = \frac{1}{e_{1t}\,e_{2t}\,e_{3t} } \left( {\delta_i \left[ {e_{2u}\,e_{3u}\,u^{n-1}_I} \right] +\delta_j \left[ {e_{1v}\,e_{3v}\,v^{n-1}_I} \right]} \right).$ By the application of \autoref{eq:asm_dmp} and \autoref{eq:asm_dmp} the divergence is filtered in each iteration, By the application of \autoref{eq:ASM_dmp} the divergence is filtered in each iteration, and the vorticity is left unchanged. In the presence of coastal boundaries with zero velocity increments perpendicular to the coast This type of the initialisation reduces the vertical velocity magnitude and alleviates the problem of the excessive unphysical vertical mixing in the first steps of the model integration \citep{Talagrand_JAS72, Dobricic_al_OS07}. \citep{talagrand_JAS72, dobricic.pinardi.ea_OS07}. Diffusion coefficients are defined as $A_D = \alpha e_{1t} e_{2t}$, where $\alpha = 0.2$. The divergence damping is activated by assigning to \np{nn\_divdmp} in the \textit{nam\_asminc} namelist The divergence damping is activated by assigning to \np{nn_divdmp}{nn\_divdmp} in the \nam{_asminc}{\_asminc} namelist a value greater than zero. By choosing this value to be of the order of 100 the increments in the vertical velocity will be significantly reduced. This specifies the number of iterations of the divergence damping. Setting a value of the order of 100 will result in a significant reduction in the vertical velocity induced by the increments. %========================================================================== %% ================================================================================================= \section{Implementation details} \label{sec:ASM_details} Here we show an example \ngn{namasm} namelist and the header of an example assimilation increments file on Here we show an example \nam{_asminc}{\_asminc} namelist and the header of an example assimilation increments file on the ORCA2 grid. %------------------------------------------namasm----------------------------------------------------- % \nlst{nam_asminc} %------------------------------------------------------------------------------------------------------------- \begin{listing} \nlst{nam_asminc} \caption{\forcode{&nam_asminc}} \label{lst:nam_asminc} \end{listing} The header of an assimilation increments file produced using the NetCDF tool \end{clines} \biblio \pindex \subinc{\input{../../global/epilogue}} \end{document}
 r10509 \begin{document} % ================================================================ % Chapter I/O & Diagnostics % ================================================================ \chapter{Output and Diagnostics (IOM, DIA, TRD, FLO)} \label{chap:DIA} \minitoc \newpage % ================================================================ %       Old Model Output % ================================================================ \section{Old model output (default)} %    {\em 4.0} & {\em Mirek Andrejczuk, Massimiliano Drudi} & {\em }  \\ %    {\em }      & {\em Dorotea Iovino, Nicolas Martin} & {\em }  \\ %    {\em 3.6} & {\em Gurvan Madec, Sebastien Masson } & {\em }  \\ %    {\em 3.4} & {\em Gurvan Madec, Rachid Benshila, Andrew Coward } & {\em }  \\ %    {\em }      & {\em Christian Ethe, Sebastien Masson } & {\em }  \\ \thispagestyle{plain} \chaptertoc \paragraph{Changes record} ~\\ {\footnotesize \begin{tabularx}{\textwidth}{l||X|X} Release & Author(s) & Modifications \\ \hline {\em   4.0} & {\em ...} & {\em ...} \\ {\em   3.6} & {\em ...} & {\em ...} \\ {\em   3.4} & {\em ...} & {\em ...} \\ {\em <=3.4} & {\em ...} & {\em ...} \end{tabularx} } \clearpage %% ================================================================================================= \section{Model output} \label{sec:DIA_io_old} the same run performed in one step. It should be noted that this requires that the restart file contains two consecutive time steps for all the prognostic variables, and that it is saved in the same binary format as the one used by the computer that is to read it (in particular, 32 bits binary IEEE format must not be used for this file). all the prognostic variables. The output listing and file(s) are predefined but should be checked and eventually adapted to the user's needs. The output listing is stored in the $ocean.output$ file. The information is printed from within the code on the logical unit $numout$. The output listing is stored in the \textit{ocean.output} file. The information is printed from within the code on the logical unit \texttt{numout}. To locate these prints, use the UNIX command "\textit{grep -i numout}" in the source code directory. A complete description of the use of this I/O server is presented in the next section. By default, \key{iomput} is not defined, NEMO produces NetCDF with the old IOIPSL library which has been kept for compatibility and its easy installation. However, the IOIPSL library is quite inefficient on parallel machines and, since version 3.2, many diagnostic options have been added presuming the use of \key{iomput}. The usefulness of the default IOIPSL-based option is expected to reduce with each new release. If \key{iomput} is not defined, output files and content are defined in the \mdl{diawri} module and contain mean (or instantaneous if \key{diainstant} is defined) values over a regular period of nn\_write time-steps (namelist parameter). %\gmcomment{                    % start of gmcomment % ================================================================ % Diagnostics % ================================================================ %\cmtgm{                    % start of gmcomment %% ================================================================================================= \section{Standard model output (IOM)} \label{sec:DIA_iom} Since version 3.2, iomput is the NEMO output interface of choice. Since version 3.2, iomput is the \NEMO\ output interface of choice. It has been designed to be simple to use, flexible and efficient. The two main purposes of iomput are: The two main purposes of iomput are: \begin{enumerate} \item The complete and flexible control of the output files through external XML files adapted by \item The complete and flexible control of the output files through external XML files adapted by the user from standard templates. \item To achieve high performance and scalable output through the optional distribution of \item To achieve high performance and scalable output through the optional distribution of all diagnostic output related tasks to dedicated processes. \end{enumerate} The first functionality allows the user to specify, without code changes or recompilation, The first functionality allows the user to specify, without code changes or recompilation, aspects of the diagnostic output stream, such as: \begin{itemize} \item The choice of output frequencies that can be different for each file (including real months and years). \item The choice of file contents; includes complete flexibility over which data are written in which files \item The choice of output frequencies that can be different for each file (including real months and years). \item The choice of file contents; includes complete flexibility over which data are written in which files (the same data can be written in different files). \item The possibility to split output files at a chosen frequency. \item The possibility to extract a vertical or an horizontal subdomain. \item The choice of the temporal operation to perform, \eg: average, accumulate, instantaneous, min, max and once. \item Control over metadata via a large XML "database" of possible output fields. \item The possibility to split output files at a chosen frequency. \item The possibility to extract a vertical or an horizontal subdomain. \item The choice of the temporal operation to perform, \eg: average, accumulate, instantaneous, min, max and once. \item Control over metadata via a large XML "database" of possible output fields. \end{itemize} in a very easy way. All details of iomput functionalities are listed in the following subsections. Examples of the XML files that control the outputs can be found in: \path{NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef.xml}, \path{NEMOGCM/CONFIG/SHARED/field_def.xml} and \path{NEMOGCM/CONFIG/SHARED/domain_def.xml}. \\ Examples of the XML files that control the outputs can be found in: \path{cfgs/ORCA2_ICE_PISCES/EXPREF/iodef.xml}, \path{cfgs/SHARED/field_def_nemo-oce.xml}, \path{cfgs/SHARED/field_def_nemo-pisces.xml}, \path{cfgs/SHARED/field_def_nemo-ice.xml} and \path{cfgs/SHARED/domain_def_nemo.xml}. \\ The second functionality targets output performance when running in parallel (\key{mpp\_mpi}). Iomput provides the possibility to specify N dedicated I/O processes (in addition to the NEMO processes) Iomput provides the possibility to specify N dedicated I/O processes (in addition to the \NEMO\ processes) to collect and write the outputs. With an appropriate choice of N by the user, the bottleneck associated with the writing of the output files can be greatly reduced. In version 3.6, the iom\_put interface depends on an external code called \href{https://forge.ipsl.jussieu.fr/ioserver/browser/XIOS/branchs/xios-1.0}{XIOS-1.0} (use of revision 618 or higher is required). In version 3.6, the \rou{iom\_put} interface depends on an external code called \href{https://forge.ipsl.jussieu.fr/ioserver/browser/XIOS/branchs/xios-2.5}{XIOS-2.5} %(use of revision 618 or higher is required). This new IO server can take advantage of the parallel I/O functionality of NetCDF4 to create a single output file and therefore to bypass the rebuilding phase. Note that writing in parallel into the same NetCDF files requires that your NetCDF4 library is linked to an HDF5 library that has been correctly compiled (\ie with the configure option $--$enable-parallel). an HDF5 library that has been correctly compiled (\ie\ with the configure option $--$enable-parallel). Note that the files created by iomput through XIOS are incompatible with NetCDF3. All post-processsing and visualization tools must therefore be compatible with NetCDF4 and not only NetCDF3. Even if not using the parallel I/O functionality of NetCDF4, using N dedicated I/O servers, where N is typically much less than the number of NEMO processors, will reduce the number of output files created. This can greatly reduce the post-processing burden usually associated with using large numbers of NEMO processors. where N is typically much less than the number of \NEMO\ processors, will reduce the number of output files created. This can greatly reduce the post-processing burden usually associated with using large numbers of \NEMO\ processors. Note that for smaller configurations, the rebuilding phase can be avoided, even without a parallel-enabled NetCDF4 library, simply by employing only one dedicated I/O server. %% ================================================================================================= \subsection{XIOS: Reading and writing restart file} XIOS may be used to read single file restart produced by NEMO. Currently only the variables written to file \forcode{numror} can be handled by XIOS. To activate restart reading using XIOS, set \np{ln\_xios\_read}\forcode{ = .true. } in \textit{namelist\_cfg}. This setting will be ignored when multiple restart files are present, and default NEMO functionality will be used for reading. There is no need to change iodef.xml file to use XIOS to read restart, all definitions are done within the NEMO code. For high resolution configuration, however, XIOS may be used to read single file restart produced by \NEMO. Currently only the variables written to file \forcode{numror} can be handled by XIOS. To activate restart reading using XIOS, set \np[=.true. ]{ln_xios_read}{ln\_xios\_read} in \textit{namelist\_cfg}. This setting will be ignored when multiple restart files are present, and default \NEMO functionality will be used for reading. There is no need to change iodef.xml file to use XIOS to read restart, all definitions are done within the \NEMO\ code. For high resolution configuration, however, there may be a need to add the following line in iodef.xml (xios context): \end{xmllines} This variable sets timeout for reading. If XIOS is to be used to read restart from file generated with an earlier NEMO version (3.6 for instance), This variable sets timeout for reading. If XIOS is to be used to read restart from file generated with an earlier \NEMO\ version (3.6 for instance), dimension \forcode{z} defined in restart file must be renamed to \forcode{nav_lev}.\\ XIOS can also be used to write NEMO restart. A namelist parameter \np{nn\_wxios} is used to determine the type of restart NEMO will write. If it is set to 0, default NEMO functionality will be used - each processor writes its own restart file; if it is set to 1 XIOS will write restart into a single file; for \np{nn\_wxios = 2} the restart will be written by XIOS into multiple files, one for each XIOS server. Note, however, that \textbf{NEMO will not read restart generated by XIOS when \np{nn\_wxios = 2}}. The restart will have to be rebuild before continuing the run. This option aims to reduce number of restart files generated by NEMO only, and may be useful when there is a need to change number of processors used to run simulation. XIOS can also be used to write \NEMO\ restart. A namelist parameter \np{nn_wxios}{nn\_wxios} is used to determine the type of restart \NEMO\ will write. If it is set to 0, default \NEMO\ functionality will be used - each processor writes its own restart file; if it is set to 1 XIOS will write restart into a single file; for \np[=2]{nn_wxios}{nn\_wxios} the restart will be written by XIOS into multiple files, one for each XIOS server. Note, however, that \textbf{\NEMO\ will not read restart generated by XIOS when \np[=2]{nn_wxios}{nn\_wxios}}. The restart will have to be rebuild before continuing the run. This option aims to reduce number of restart files generated by \NEMO\ only, and may be useful when there is a need to change number of processors used to run simulation. If an additional variable must be written to a restart file, the following steps are needed: \begin{description} \item[step 1:] add variable name to a list of restart variables (in subroutine \rou{iom\_set\_rst\_vars,} \mdl{iom}) and define correct grid for the variable (\forcode{grid_N_3D} - 3D variable, \forcode{grid_N} - 2D variable, \forcode{grid_vector} - \begin{enumerate} \item Add variable name to a list of restart variables (in subroutine \rou{iom\_set\_rst\_vars,} \mdl{iom}) and define correct grid for the variable (\forcode{grid_N_3D} - 3D variable, \forcode{grid_N} - 2D variable, \forcode{grid_vector} - 1D variable, \forcode{grid_scalar} - scalar), \item[step 2:] add variable to the list of fields written by restart.  This can be done either in subroutine \rou{iom\_set\_rstw\_core} (\mdl{iom}) or by calling  \rou{iom\_set\_rstw\_active} (\mdl{iom}) with the name of a variable as an argument. This convention follows approach for writing restart using iom, where variables are \item Add variable to the list of fields written by restart.  This can be done either in subroutine \rou{iom\_set\_rstw\_core} (\mdl{iom}) or by calling  \rou{iom\_set\_rstw\_active} (\mdl{iom}) with the name of a variable as an argument. This convention follows approach for writing restart using iom, where variables are written either by \rou{rst\_write} or by calling \rou{iom\_rstput} from individual routines. \end{description} \end{enumerate} An older versions of XIOS do not support reading functionality. It's recommended to use at least XIOS2@1451. %% ================================================================================================= \subsection{XIOS: XML Inputs-Outputs Server} %% ================================================================================================= \subsubsection{Attached or detached mode?} \xmlline|| The {\tt using\_server} setting determines whether or not the server will be used in \textit{attached mode} (as a library) [{\tt> false <}] or in \textit{detached mode} (as an external executable on N additional, dedicated cpus) [{\tt > true <}]. The \textit{attached mode} is simpler to use but much less efficient for massively parallel applications. The \texttt{using\_server} setting determines whether or not the server will be used in \textit{attached mode} (as a library) [\texttt{> false <}] or in \textit{detached mode} (as an external executable on N additional, dedicated cpus) [\texttt{ > true <}]. The \textit{attached mode} is simpler to use but much less efficient for massively parallel applications. The type of each file can be either ''multiple\_file'' or ''one\_file''. In \textit{attached mode} and if the type of file is ''multiple\_file'', then each NEMO process will also act as an IO server and produce its own set of output files. then each \NEMO\ process will also act as an IO server and produce its own set of output files. Superficially, this emulates the standard behaviour in previous versions. However, the subdomain written out by each process does not correspond to write to its own set of output files. If the ''one\_file'' option is chosen then all XIOS processes will collect their longitudinal strips and write (in parallel) to a single output file. write (in parallel) to a single output file. Note running in detached mode requires launching a Multiple Process Multiple Data (MPMD) parallel job. The following subsection provides a typical example but the syntax will vary in different MPP environments. %% ================================================================================================= \subsubsection{Number of cpu used by XIOS in detached mode} The number of cores used by the XIOS is specified when launching the model. The number of cores dedicated to XIOS should be from \texttildelow1/10 to \texttildelow1/50 of the number of cores dedicated to NEMO. The number of cores dedicated to XIOS should be from \texttildelow1/10 to \texttildelow1/50 of the number of cores dedicated to \NEMO. Some manufacturers suggest using O($\sqrt{N}$) dedicated IO processors for N processors but this is a general recommendation and not specific to NEMO. this is a general recommendation and not specific to \NEMO. It is difficult to provide precise recommendations because the optimal choice will depend on the particular hardware properties of the target system the particular hardware properties of the target system (parallel filesystem performance, available memory, memory bandwidth etc.) and the volume and frequency of data to be created. \cmd|mpirun -np 62 ./nemo.exe : -np 2 ./xios_server.exe| %% ================================================================================================= \subsubsection{Control of XIOS: the context in iodef.xml} As well as the {\tt using\_server} flag, other controls on the use of XIOS are set in the XIOS context in iodef.xml. As well as the \texttt{using\_server} flag, other controls on the use of XIOS are set in the XIOS context in \textit{iodef.xml}. See the XML basics section below for more details on XML syntax and rules. \begin{table} \scriptsize \begin{tabularx}{\textwidth}{|lXl|} \hline \hline buffer\_size                                                            & buffer size used by XIOS to send data from NEMO to XIOS. buffer size used by XIOS to send data from \NEMO\ to XIOS. Larger is more efficient. Note that needed/used buffer sizes are summarized at the end of the job & \hline buffer\_server\_factor\_size                                            & ratio between NEMO and XIOS buffer size. ratio between \NEMO\ and XIOS buffer size. Should be 2.                                                            & 2        \\ \hline oasis\_codes\_id                                                        & when using oasis, define the identifier of NEMO in the namcouple. when using oasis, define the identifier of \NEMO\ in the namcouple. Note that the identifier of XIOS is xios.x                              & oceanx   \\ \end{table} %% ================================================================================================= \subsection{Practical issues} %% ================================================================================================= \subsubsection{Installation} As mentioned, XIOS is supported separately and must be downloaded and compiled before it can be used with NEMO. As mentioned, XIOS is supported separately and must be downloaded and compiled before it can be used with \NEMO. See the installation guide on the \href{http://forge.ipsl.jussieu.fr/ioserver/wiki}{XIOS} wiki for help and guidance. NEMO will need to link to the compiled XIOS library. The \href{https://forge.ipsl.jussieu.fr/nemo/wiki/Users/ModelInterfacing/InputsOutputs#Inputs-OutputsusingXIOS} {XIOS with NEMO} guide provides an example illustration of how this can be achieved. \NEMO\ will need to link to the compiled XIOS library. The \href{https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/guide/html/install.html#extract-and-install-xios} {Extract and install XIOS} guide provides an example illustration of how this can be achieved. %% ================================================================================================= \subsubsection{Add your own outputs} \begin{enumerate} \item[1.] in NEMO code, add a \forcode{CALL iom\_put( 'identifier', array )} where you want to output a 2D or 3D array. \item[2.] If necessary, add \forcode{USE iom ! I/O manager library} to the list of used modules in \item in \NEMO\ code, add a \forcode{CALL iom_put( 'identifier', array )} where you want to output a 2D or 3D array. \item If necessary, add \forcode{USE iom ! I/O manager library} to the list of used modules in the upper part of your module. \item[3.] in the field\_def.xml file, add the definition of your variable using the same identifier you used in the f90 code \item in the field\_def.xml file, add the definition of your variable using the same identifier you used in the f90 code (see subsequent sections for a details of the XML syntax and rules). For example: \begin{xmllines} ... ... \end{xmllines} \end{xmllines} Note your definition must be added to the field\_group whose reference grid is consistent with the size of the array passed to iomput. (iom\_set\_domain\_attr and iom\_set\_axis\_attr in \mdl{iom}) or defined in the domain\_def.xml file. \eg: \begin{xmllines} \end{xmllines} Note, if your array is computed within the surface module each \np{nn\_fsbc} time\_step, Note, if your array is computed within the surface module each \np{nn_fsbc}{nn\_fsbc} time\_step, add the field definition within the field\_group defined with the id "SBC": \xmlcode{} which has been defined with the correct frequency of operations (iom\_set\_field\_attr in \mdl{iom}) \item[4.] add your field in one of the output files defined in iodef.xml \item add your field in one of the output files defined in iodef.xml (again see subsequent sections for syntax and rules) \begin{xmllines} \begin{xmllines} ... ... \end{xmllines} \end{xmllines} \end{enumerate} %% ================================================================================================= \subsection{XML fundamentals} %% ================================================================================================= \subsubsection{ XML basic rules} See \href{http://www.xmlnews.org/docs/xml-basics.html}{here} for more details. \subsubsection{Structure of the XML file used in NEMO} %% ================================================================================================= \subsubsection{Structure of the XML file used in \NEMO} The XML file used in XIOS is structured by 7 families of tags: \begin{table} \scriptsize \begin{tabular*}{\textwidth}{|p{0.15\textwidth}p{0.4\textwidth}p{0.35\textwidth}|} \hline \begin{table} \scriptsize \begin{tabular}{|p{0.15\textwidth}p{0.4\textwidth}p{0.35\textwidth}|} \hline \xmlcode{}   \\ \hline context nemo    &   context containing IO information for NEMO (mother grid when using AGRIF)  & context nemo    &   context containing IO information for \NEMO\ (mother grid when using AGRIF)  & \xmlcode{}   \\ \hline context 1\_nemo &   context containing IO information for NEMO child grid 1 (when using AGRIF) & context 1\_nemo &   context containing IO information for \NEMO\ child grid 1 (when using AGRIF) & \xmlcode{} \\ \hline context n\_nemo &   context containing IO information for NEMO child grid n (when using AGRIF) & context n\_nemo &   context containing IO information for \NEMO\ child grid n (when using AGRIF) & \xmlcode{} \\ \hline \begin{table} \scriptsize \begin{tabular}{|p{0.15\textwidth}p{0.4\textwidth}p{0.35\textwidth}|} \hline \end{table} \noindent Each context tag related to NEMO (mother or child grids) is divided into 5 parts \noindent Each context tag related to \NEMO\ (mother or child grids) is divided into 5 parts (that can be defined in any order): \begin{table} \scriptsize \begin{tabular}{|p{0.15\textwidth}p{0.4\textwidth}p{0.35\textwidth}|} \hline \end{table} %% ================================================================================================= \subsubsection{Nesting XML files} The inclusion of XML files into the main XML file can be done through the attribute src: \xmlline|| \noindent In NEMO, by default, the field and domain definition is done in 2 separate files: \path{NEMOGCM/CONFIG/SHARED/field_def.xml} and \path{NEMOGCM/CONFIG/SHARED/domain_def.xml} that \noindent In \NEMO, by default, the field definition is done in 3 separate files ( \path{cfgs/SHARED/field_def_nemo-oce.xml}, \path{cfgs/SHARED/field_def_nemo-pisces.xml} and \path{cfgs/SHARED/field_def_nemo-ice.xml} ) and the  domain definition is done in another file ( \path{cfgs/SHARED/domain_def_nemo.xml} ) that are included in the main iodef.xml file through the following commands: \begin{xmllines} \end{xmllines} \end{xmllines} %% ================================================================================================= \subsubsection{Use of inheritance} \begin{xmllines}                 \end{xmllines}     \end{xmllines} Inherit (and overwrite, if needed) the attributes of a tag you are refering to. %% ================================================================================================= \subsubsection{Use of groups} Secondly, the group can be used to replace a list of elements. Several examples of groups of fields are proposed at the end of the file \path{CONFIG/SHARED/field_def.xml}. Several examples of groups of fields are proposed at the end of the XML field files ( \path{cfgs/SHARED/field_def_nemo-oce.xml}, \path{cfgs/SHARED/field_def_nemo-pisces.xml} and \path{cfgs/SHARED/field_def_nemo-ice.xml} ) . For example, a short list of the usual variables related to the U grid:   \end{xmllines} \end{xmllines} %% ================================================================================================= \subsection{Detailed functionalities} the new functionalities offered by the XML interface of XIOS. %% ================================================================================================= \subsubsection{Define horizontal subdomains} Horizontal subdomains are defined through the attributs zoom\_ibegin, zoom\_jbegin, zoom\_ni, zoom\_nj of the tag family domain. It must therefore be done in the domain part of the XML file. For example, in \path{CONFIG/SHARED/domain_def.xml}, we provide the following example of a definition of It must therefore be done in the domain part of the XML file. For example, in \path{cfgs/SHARED/domain_def.xml}, we provide the following example of a definition of a 5 by 5 box with the bottom left corner at point (10,10). \begin{xmllines} \end{xmllines} \begin{xmllines} \end{xmllines} Moorings are seen as an extrem case corresponding to a 1 by 1 subdomain. Moorings are seen as an extrem case corresponding to a 1 by 1 subdomain. The Equatorial section, the TAO, RAMA and PIRATA moorings are already registered in the code and can therefore be outputted without taking care of their (i,j) position in the grid. \end{xmllines} Note that if the domain decomposition used in XIOS cuts the subdomain in several parts and if you use the ''multiple\_file'' type for your output files, you will endup with several files you will need to rebuild using unprovided tools (like ncpdq and ncrcat, Note that if the domain decomposition used in XIOS cuts the subdomain in several parts and if you use the ''multiple\_file'' type for your output files, you will endup with several files you will need to rebuild using unprovided tools (like ncpdq and ncrcat, \href{http://nco.sourceforge.net/nco.html#Concatenation}{see nco manual}). We are therefore advising to use the ''one\_file'' type in this case. %% ================================================================================================= \subsubsection{Define vertical zooms} Vertical zooms are defined through the attributs zoom\_begin and zoom\_end of the tag family axis. Vertical zooms are defined through the attributs zoom\_begin and zoom\_n of the tag family axis. It must therefore be done in the axis part of the XML file. For example, in \path{NEMOGCM/CONFIG/ORCA2_LIM/iodef_demo.xml}, we provide the following example: \begin{xmllines} For example, in \path{cfgs/ORCA2_ICE_PISCES/EXPREF/iodef_demo.xml}, we provide the following example: \begin{xmllines} \end{xmllines} \end{xmllines} %% ================================================================================================= \subsubsection{Control of the output file names} \begin{xmllines} ... \begin{table} \scriptsize \begin{tabularx}{\textwidth}{|lX|} \hline \end{table} \noindent For example, \noindent For example, \xmlline|| \noindent will give the following file name radical: \ifile{myfile\_ORCA2\_19891231\_freq1d} \subsubsection{Other controls of the XML attributes from NEMO} The values of some attributes are defined by subroutine calls within NEMO %% ================================================================================================= \subsubsection{Other controls of the XML attributes from \NEMO} The values of some attributes are defined by subroutine calls within \NEMO (calls to iom\_set\_domain\_attr, iom\_set\_axis\_attr and iom\_set\_field\_attr in \mdl{iom}). Any definition given in the XML file will be overwritten. \begin{table} \scriptsize \begin{tabularx}{\textwidth}{|X|c|c|c|} \begin{tabular}{|l|c|c|} \hline tag ids affected by automatic definition of some of their attributes & name attribute                                                       & attribute value                      \\ attribute value                                                      \\ \hline \hline field\_definition                                                    & freq\_op                                                             & \np{rn\_rdt}                         \\ \np{rn_rdt}{rn\_rdt}                                                         \\ \hline SBC                                                                  & freq\_op                                                             & \np{rn\_rdt} $\times$ \np{nn\_fsbc}  \\ \np{rn_rdt}{rn\_rdt} $\times$ \np{nn_fsbc}{nn\_fsbc}                                  \\ \hline ptrc\_T                                                              & freq\_op                                                             & \np{rn\_rdt} $\times$ \np{nn\_dttrc} \\ \np{rn_rdt}{rn\_rdt} $\times$ \np{nn_dttrc}{nn\_dttrc}                                \\ \hline diad\_T                                                              & freq\_op                                                             & \np{rn\_rdt} $\times$ \np{nn\_dttrc} \\ \np{rn_rdt}{rn\_rdt} $\times$ \np{nn_dttrc}{nn\_dttrc}                                \\ \hline EqT, EqU, EqW                                                        & jbegin, ni,                                                          & according to the grid                \\ & according to the grid                                                \\ & name\_suffix                                                         & \\ \\ \hline TAO, RAMA and PIRATA moorings                                        & zoom\_ibegin, zoom\_jbegin,                                          & according to the grid                \\ & according to the grid                                                \\ & name\_suffix                                                         & \\ \hline \end{tabularx} \\ \hline \end{tabular} \end{table} %% ================================================================================================= \subsubsection{Advanced use of XIOS functionalities} %% ================================================================================================= \subsection{XML reference tables} \label{subsec:IOM_xmlref} \label{subsec:DIA_IOM_xmlref} \begin{enumerate} \item Simple computation: directly define the computation when refering to the variable in the file definition. \item Simple computation: directly define the computation when refering to the variable in the file definition. \begin{xmllines} \end{xmllines} \item Simple computation: define a new variable and use it in the file definition. \item Simple computation: define a new variable and use it in the file definition. in field\_definition: sst2 won't be evaluated. \item Change of variable precision: \item Change of variable precision: \begin{xmllines} \end{xmllines} Note that, then the code is crashing, writting real4 variables forces a numerical convection from Note that, then the code is crashing, writting real4 variables forces a numerical conversion from real8 to real4 which will create an internal error in NetCDF and will avoid the creation of the output files. Forcing double precision outputs with prec="8" (for example in the field\_definition) will avoid this problem. \item add user defined attributes: \begin{xmllines} \item add user defined attributes: \begin{xmllines} blabla_global \end{xmllines} \item use of the @'' function: example 1, weighted temporal average \end{xmllines} \item use of the @'' function: example 1, weighted temporal average - define a new variable in field\_definition \begin{xmllines}     @toce_e3t / @e3t \end{xmllines} Note that in this case, freq\_op must be equal to the file output\_freq. \item use of the @'' function: example 2, monthly SSH standard deviation \item use of the @'' function: example 2, monthly SSH standard deviation - define a new variable in field\_definition \begin{xmllines}     sqrt( @ssh2 - @ssh * @ssh ) \end{xmllines} here we use the default, average. So, in the above case, @ssh2 will do the monthly mean of ssh*ssh. Operation="instant" refers to the temporal operation to be performed on the field ''sqrt( @ssh2 - @ssh * @ssh )'': Operation="instant" refers to the temporal operation to be performed on the field ''sqrt( @ssh2 - @ssh * @ssh )'': here the temporal average is alreday done by the @'' function so we just use instant. field\_ref="ssh" means that attributes not explicitely defined, are inherited from ssh field. Note that in this case, freq\_op must be equal to the file output\_freq. \item use of the @'' function: example 3, monthly average of SST diurnal cycle \item use of the @'' function: example 3, monthly average of SST diurnal cycle - define 2 new variables in field\_definition \begin{xmllines}     \end{xmllines} field\_ref="sst" means that attributes not explicitely defined, are inherited from sst field. %% ================================================================================================= \subsubsection{Tag list per family} \begin{table} \scriptsize \begin{tabularx}{\textwidth}{|l|X|X|l|X|} \hline \hline \end{tabularx} \caption{Context tags} \caption{XIOS: context tags} \end{table} \begin{table} \scriptsize \begin{tabularx}{\textwidth}{|l|X|X|X|l|} \hline \hline \end{tabularx} \caption{Field tags ("\tt{field\_*}")} \caption{XIOS: field tags ("\texttt{field\_*}")} \end{table} \begin{table} \scriptsize \begin{tabularx}{\textwidth}{|l|X|X|X|l|} \hline \hline \end{tabularx} \caption{File tags ("\tt{file\_*}")} \caption{XIOS: file tags ("\texttt{file\_*}")} \end{table} \begin{table} \scriptsize \begin{tabularx}{\textwidth}{|l|X|X|X|X|} \hline \hline \end{tabularx} \caption{Axis tags ("\tt{axis\_*}")} \caption{XIOS: axis tags ("\texttt{axis\_*}")} \end{table} \begin{table} \scriptsize \begin{tabularx}{\textwidth}{|l|X|X|X|X|} \hline \hline \end{tabularx} \caption{Domain tags ("\tt{domain\_*)}"} \caption{XIOS: domain tags ("\texttt{domain\_*)}"} \end{table} \begin{table} \scriptsize \begin{tabularx}{\textwidth}{|l|X|X|X|X|} \hline \hline \end{tabularx} \caption{Grid tags ("\tt{grid\_*}")} \caption{XIOS: grid tags ("\texttt{grid\_*}")} \end{table} %% ================================================================================================= \subsubsection{Attributes list per family} \begin{table} \scriptsize \begin{tabularx}{\textwidth}{|l|X|l|l|} \hline \hline \end{tabularx} \caption{Reference attributes ("\tt{*\_ref}")} \caption{XIOS: reference attributes ("\texttt{*\_ref}")} \end{table} \begin{table} \scriptsize \begin{tabularx}{\textwidth}{|l|X|l|l|} \hline \hline \end{tabularx} \caption{Domain attributes ("\tt{zoom\_*}")} \caption{XIOS: domain attributes ("\texttt{zoom\_*}")} \end{table} \begin{table} \scriptsize \begin{tabularx}{\textwidth}{|l|X|l|l|} \hline \hline \end{tabularx} \caption{File attributes} \caption{XIOS: file attributes} \end{table} \begin{table} \scriptsize \begin{tabularx}{\textwidth}{|l|X|l|l|} \hline \hline \end{tabularx} \caption{Field attributes} \caption{XIOS: field attributes} \end{table} \begin{table} \scriptsize \begin{tabularx}{\textwidth}{|l|X|X|X|} \hline \hline \end{tabularx} \caption{Miscellaneous attributes} \caption{XIOS: miscellaneous attributes} \end{table} %% ================================================================================================= \subsection{CF metadata standard compliance} Output from the XIOS-1.0 IO server is compliant with Output from the XIOS IO server is compliant with \href{http://cfconventions.org/Data/cf-conventions/cf-conventions-1.5/build/cf-conventions.html}{version 1.5} of the CF metadata standard. the CF metadata standard. Therefore while a user may wish to add their own metadata to the output files (as demonstrated in example 4 of section \autoref{subsec:IOM_xmlref}) the metadata should, for the most part, comply with the CF-1.5 standard. section \autoref{subsec:DIA_IOM_xmlref}) the metadata should, for the most part, comply with the CF-1.5 standard. Some metadata that may significantly increase the file size (horizontal cell areas and vertices) are controlled by the namelist parameter \np{ln\_cfmeta} in the \ngn{namrun} namelist. the namelist parameter \np{ln_cfmeta}{ln\_cfmeta} in the \nam{run}{run} namelist. This must be set to true if these metadata are to be included in the output files. % ================================================================ %       NetCDF4 support % ================================================================ \section{NetCDF4 support (\protect\key{netcdf4})} %% ================================================================================================= \section[NetCDF4 support (\texttt{\textbf{key\_netcdf4}})]{NetCDF4 support (\protect\key{netcdf4})} \label{sec:DIA_nc4} Chunking and compression can lead to significant reductions in file sizes for a small runtime overhead. For a fuller discussion on chunking and other performance issues the reader is referred to the NetCDF4 documentation found \href{http://www.unidata.ucar.edu/software/netcdf/docs/netcdf.html#Chunking}{here}. the NetCDF4 documentation found \href{https://www.unidata.ucar.edu/software/netcdf/docs/netcdf_perf_chunking.html}{here}. The new features are only available when the code has been linked with a NetCDF4 library Datasets created with chunking and compression are not backwards compatible with NetCDF3 "classic" format but most analysis codes can be relinked simply with the new libraries and will then read both NetCDF3 and NetCDF4 files. NEMO executables linked with NetCDF4 libraries can be made to produce NetCDF3 files by setting the \np{ln\_nc4zip} logical to false in the \textit{namnc4} namelist: %------------------------------------------namnc4---------------------------------------------------- \nlst{namnc4} %------------------------------------------------------------------------------------------------------------- \NEMO\ executables linked with NetCDF4 libraries can be made to produce NetCDF3 files by setting the \np{ln_nc4zip}{ln\_nc4zip} logical to false in the \nam{nc4}{nc4} namelist: \begin{listing} \nlst{namnc4} \caption{\forcode{&namnc4}} \label{lst:namnc4} \end{listing} If \key{netcdf4} has not been defined, these namelist parameters are not read. In this case, \np{ln\_nc4zip} is set false and dummy routines for a few NetCDF4-specific functions are defined. In this case, \np{ln_nc4zip}{ln\_nc4zip} is set false and dummy routines for a few NetCDF4-specific functions are defined. These functions will not be used but need to be included so that compilation is possible with NetCDF3 libraries. \end{forlines} \noindent for a standard ORCA2\_LIM configuration gives chunksizes of {\small\tt 46x38x1} respectively in the mono-processor case (\ie global domain of {\small\tt 182x149x31}). An illustration of the potential space savings that NetCDF4 chunking and compression provides is given in table \autoref{tab:NC4} which compares the results of two short runs of the ORCA2\_LIM reference configuration with \noindent for a standard ORCA2\_LIM configuration gives chunksizes of {\small\texttt 46x38x1} respectively in the mono-processor case (\ie\ global domain of {\small\texttt 182x149x31}). An illustration of the potential space savings that NetCDF4 chunking and compression provides is given in table \autoref{tab:DIA_NC4} which compares the results of two short runs of the ORCA2\_LIM reference configuration with a 4x2 mpi partitioning. Note the variation in the compression ratio achieved which reflects chiefly the dry to wet volume ratio of each processing region. %------------------------------------------TABLE---------------------------------------------------- \begin{table} \scriptsize \centering \begin{tabular}{lrrr} ORCA2\_2d\_grid\_W\_0007.nc & 4416    & 1368     & 70\%      \\ \end{tabular} \caption{ \protect\label{tab:NC4} Filesize comparison between NetCDF3 and NetCDF4 with chunking and compression } \caption{Filesize comparison between NetCDF3 and NetCDF4 with chunking and compression} \label{tab:DIA_NC4} \end{table} %---------------------------------------------------------------------------------------------------- When \key{iomput} is activated with \key{netcdf4} chunking and compression parameters for fields produced via \np{iom\_put} calls are set via an equivalent and identically named namelist to \textit{namnc4} in \np{xmlio\_server.def}. Typically this namelist serves the mean files whilst the \ngn{ namnc4} in the main namelist file continues to \rou{iom\_put} calls are set via an equivalent and identically named namelist to \nam{nc4}{nc4} in \textit{xmlio\_server.def}. Typically this namelist serves the mean files whilst the \nam{nc4}{nc4} in the main namelist file continues to serve the restart files. This duplication is unfortunate but appropriate since, if using io\_servers, the domain sizes of the individual files produced by the io\_server processes may be different to those produced by the invidual processing regions and different chunking choices may be desired. % ------------------------------------------------------------------------------------------------------------- %       Tracer/Dynamics Trends % ------------------------------------------------------------------------------------------------------------- \section{Tracer/Dynamics trends  (\protect\ngn{namtrd})} %% ================================================================================================= \section[Tracer/Dynamics trends (\forcode{&namtrd})]{Tracer/Dynamics trends (\protect\nam{trd}{trd})} \label{sec:DIA_trd} %------------------------------------------namtrd---------------------------------------------------- \nlst{namtrd} %------------------------------------------------------------------------------------------------------------- \begin{listing} \nlst{namtrd} \caption{\forcode{&namtrd}} \label{lst:namtrd} \end{listing} Each trend of the dynamics and/or temperature and salinity time evolution equations can be send to \mdl{trddyn} and/or \mdl{trdtra} modules (see TRD directory) just after their computation (\ie at the end of each $dyn\cdots.F90$ and/or $tra\cdots.F90$ routines). This capability is controlled by options offered in \ngn{namtrd} namelist. Note that the output are done with xIOS, and therefore the \key{IOM} is required. What is done depends on the \ngn{namtrd} logical set to \forcode{.true.}: (\ie\ at the end of each \textit{dyn....F90} and/or \textit{tra....F90} routines). This capability is controlled by options offered in \nam{trd}{trd} namelist. Note that the output are done with XIOS, and therefore the \key{iomput} is required. What is done depends on the \nam{trd}{trd} logical set to \forcode{.true.}: \begin{description} \item[\np{ln\_glo\_trd}]: at each \np{nn\_trd} time-step a check of the basin averaged properties of \item [{\np{ln_glo_trd}{ln\_glo\_trd}}]: at each \np{nn_trd}{nn\_trd} time-step a check of the basin averaged properties of the momentum and tracer equations is performed. This also includes a check of $T^2$, $S^2$, $\tfrac{1}{2} (u^2+v2)$, and potential energy time evolution equations properties; \item[\np{ln\_dyn\_trd}]: each 3D trend of the evolution of the two momentum components is output; \item[\np{ln\_dyn\_mxl}]: each 3D trend of the evolution of the two momentum components averaged over the mixed layer is output; \item[\np{ln\_vor\_trd}]: a vertical summation of the moment tendencies is performed, \item [{\np{ln_dyn_trd}{ln\_dyn\_trd}}]: each 3D trend of the evolution of the two momentum components is output; \item [{\np{ln_dyn_mxl}{ln\_dyn\_mxl}}]: each 3D trend of the evolution of the two momentum components averaged over the mixed layer is output; \item [{\np{ln_vor_trd}{ln\_vor\_trd}}]: a vertical summation of the moment tendencies is performed, then the curl is computed to obtain the barotropic vorticity tendencies which are output; \item[\np{ln\_KE\_trd}] : each 3D trend of the Kinetic Energy equation is output; \item[\np{ln\_tra\_trd}]: each 3D trend of the evolution of temperature and salinity is output; \item[\np{ln\_tra\_mxl}]: each 2D trend of the evolution of temperature and salinity averaged over the mixed layer is output; \item [{\np{ln_KE_trd}{ln\_KE\_trd}}]  : each 3D trend of the Kinetic Energy equation is output; \item [{\np{ln_tra_trd}{ln\_tra\_trd}}]: each 3D trend of the evolution of temperature and salinity is output; \item [{\np{ln_tra_mxl}{ln\_tra\_mxl}}]: each 2D trend of the evolution of temperature and salinity averaged over the mixed layer is output; \end{description} Note that the mixed layer tendency diagnostic can also be used on biogeochemical models via the \key{trdtrc} and \key{trdmld\_trc} CPP keys. Note that the mixed layer tendency diagnostic can also be used on biogeochemical models via the \key{trdtrc} and \key{trdmxl\_trc} CPP keys. \textbf{Note that} in the current version (v3.6), many changes has been introduced but not fully tested. In particular, options associated with \np{ln\_dyn\_mxl}, \np{ln\_vor\_trd}, and \np{ln\_tra\_mxl} are not working, and none of the options have been tested with variable volume (\ie \key{vvl} defined). % ------------------------------------------------------------------------------------------------------------- %       On-line Floats trajectories % ------------------------------------------------------------------------------------------------------------- \section{FLO: On-Line Floats trajectories (\protect\key{floats})} \label{sec:FLO} %--------------------------------------------namflo------------------------------------------------------- \nlst{namflo} %-------------------------------------------------------------------------------------------------------------- In particular, options associated with \np{ln_dyn_mxl}{ln\_dyn\_mxl}, \np{ln_vor_trd}{ln\_vor\_trd}, and \np{ln_tra_mxl}{ln\_tra\_mxl} are not working, and none of the options have been tested with variable volume (\ie\ \np[=.true.]{ln_linssh}{ln\_linssh}). %% ================================================================================================= \section[FLO: On-Line Floats trajectories (\texttt{\textbf{key\_floats}})]{FLO: On-Line Floats trajectories (\protect\key{floats})} \label{sec:DIA_FLO} \begin{listing} \nlst{namflo} \caption{\forcode{&namflo}} \label{lst:namflo} \end{listing} The on-line computation of floats advected either by the three dimensional velocity field or constraint to remain at a given depth ($w = 0$ in the computation) have been introduced in the system during the CLIPPER project. Options are defined by \ngn{namflo} namelis variables. The algorithm used is based either on the work of \cite{Blanke_Raynaud_JPO97} (default option), or on a $4^th$ Runge-Hutta algorithm (\np{ln\_flork4}\forcode{ = .true.}). Note that the \cite{Blanke_Raynaud_JPO97} algorithm have the advantage of providing trajectories which Options are defined by \nam{flo}{flo} namelist variables. The algorithm used is based either on the work of \cite{blanke.raynaud_JPO97} (default option), or on a $4^th$ Runge-Hutta algorithm (\np[=.true.]{ln_flork4}{ln\_flork4}). Note that the \cite{blanke.raynaud_JPO97} algorithm have the advantage of providing trajectories which are consistent with the numeric of the code, so that the trajectories never intercept the bathymetry. %% ================================================================================================= \subsubsection{Input data: initial coordinates} Initial coordinates can be given with Ariane Tools convention (IJK coordinates, (\np{ln\_ariane}\forcode{ = .true.}) ) or with longitude and latitude. In case of Ariane convention, input filename is \np{init\_float\_ariane}. (IJK coordinates, (\np[=.true.]{ln_ariane}{ln\_ariane}) ) or with longitude and latitude. In case of Ariane convention, input filename is \textit{init\_float\_ariane}. Its format is: \\ {\scriptsize \texttt{I J K nisobfl itrash itrash}} { \texttt{I J K nisobfl itrash}} \noindent with: - I,J,K  : indexes of initial position - nisobfl: 0 for an isobar float, 1 for a float following the w velocity - nisobfl: 0 for an isobar float, 1 for a float following the w velocity - itrash : set to zero; it is a dummy variable to respect Ariane Tools convention \noindent Example: \\ \noindent {\scriptsize { \texttt{ 100.00000  90.00000  -1.50000 1.00000   0.00000   \\ In the other case (longitude and latitude), input filename is init\_float. Its format is: \\ {\scriptsize \texttt{Long Lat depth nisobfl ngrpfl itrash}} { \texttt{Long Lat depth nisobfl ngrpfl itrash}} \noindent with: \noindent Example: \\ \noindent {\scriptsize { \texttt{ 20.0 0.0 0.0 0 1 1    \\ } \\ \np{jpnfl} is the total number of floats during the run. When initial positions are read in a restart file (\np{ln\_rstflo}\forcode{ = .true.} ), \np{jpnflnewflo} can be added in the initialization file. \np{jpnfl}{jpnfl} is the total number of floats during the run. When initial positions are read in a restart file (\np[=.true.]{ln_rstflo}{ln\_rstflo} ), \np{jpnflnewflo}{jpnflnewflo} can be added in the initialization file. %% ================================================================================================= \subsubsection{Output data} \np{nn\_writefl} is the frequency of writing in float output file and \np{nn\_stockfl} is the frequency of \np{nn_writefl}{nn\_writefl} is the frequency of writing in float output file and \np{nn_stockfl}{nn\_stockfl} is the frequency of creation of the float restart file. Output data can be written in ascii files (\np{ln\_flo\_ascii}\forcode{ = .true.}). Output data can be written in ascii files (\np[=.true.]{ln_flo_ascii}{ln\_flo\_ascii}). In that case, output filename is trajec\_float. Another possiblity of writing format is Netcdf (\np{ln\_flo\_ascii}\forcode{ = .false.}). There are 2 possibilities: - if (\key{iomput}) is used, outputs are selected in  iodef.xml. Another possiblity of writing format is Netcdf (\np[=.false.]{ln_flo_ascii}{ln\_flo\_ascii}) with \key{iomput} and outputs selected in iodef.xml. Here it is an example of specification to put in files description section: \end{xmllines} -  if (\key{iomput}) is not used, a file called \ifile{trajec\_float} will be created by IOIPSL library. See also \href{http://stockage.univ-brest.fr/~grima/Ariane/}{here} the web site describing the off-line use of this marvellous diagnostic tool. % ------------------------------------------------------------------------------------------------------------- %       Harmonic analysis of tidal constituents % ------------------------------------------------------------------------------------------------------------- \section{Harmonic analysis of tidal constituents (\protect\key{diaharm}) } %% ================================================================================================= \section[Harmonic analysis of tidal constituents (\texttt{\textbf{key\_diaharm}})]{Harmonic analysis of tidal constituents (\protect\key{diaharm})} \label{sec:DIA_diag_harm} %------------------------------------------namdia_harm---------------------------------------------------- % \nlst{nam_diaharm} %---------------------------------------------------------------------------------------------------------- \begin{listing} \nlst{nam_diaharm} \caption{\forcode{&nam_diaharm}} \label{lst:nam_diaharm} \end{listing} A module is available to compute the amplitude and phase of tidal waves. This on-line Harmonic analysis is actived with \key{diaharm}. Some parameters are available in namelist \ngn{namdia\_harm}: - \np{nit000\_han} is the first time step used for harmonic analysis - \np{nitend\_han} is the  last time step used for harmonic analysis - \np{nstep\_han}  is the  time step frequency for harmonic analysis - \np{nb\_ana}     is the number of harmonics to analyse - \np{tname}       is an array with names of tidal constituents to analyse \np{nit000\_han} and \np{nitend\_han} must be between \np{nit000} and \np{nitend} of the simulation. Some parameters are available in namelist \nam{_diaharm}{\_diaharm}: - \np{nit000_han}{nit000\_han} is the first time step used for harmonic analysis - \np{nitend_han}{nitend\_han} is the  last time step used for harmonic analysis - \np{nstep_han}{nstep\_han}  is the  time step frequency for harmonic analysis % - \np{nb_ana}{nb\_ana}     is the number of harmonics to analyse - \np{tname}{tname}       is an array with names of tidal constituents to analyse \np{nit000_han}{nit000\_han} and \np{nitend_han}{nitend\_han} must be between \np{nit000}{nit000} and \np{nitend}{nitend} of the simulation. The restart capability is not implemented. We obtain in output $C_{j}$ and $S_{j}$ for each tidal wave. % ------------------------------------------------------------------------------------------------------------- %       Sections transports % ------------------------------------------------------------------------------------------------------------- \section{Transports across sections (\protect\key{diadct}) } %% ================================================================================================= \section[Transports across sections (\texttt{\textbf{key\_diadct}})]{Transports across sections (\protect\key{diadct})} \label{sec:DIA_diag_dct} %------------------------------------------namdct---------------------------------------------------- \nlst{namdct} %------------------------------------------------------------------------------------------------------------- \begin{listing} \nlst{nam_diadct} \caption{\forcode{&nam_diadct}} \label{lst:nam_diadct} \end{listing} A module is available to compute the transport of volume, heat and salt through sections. Each section is defined by the coordinates of its 2 extremities. The pathways between them are contructed using tools which can be found in \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT} and are written in a binary file \texttt{section\_ijglobal.diadct\_ORCA2\_LIM} which is later read in by NEMO to compute on-line transports. The pathways between them are contructed using tools which can be found in \texttt{tools/SECTIONS\_DIADCT} and are written in a binary file \texttt{section\_ijglobal.diadct} which is later read in by \NEMO\ to compute on-line transports. The on-line transports module creates three output ascii files: - \texttt{salt\_transport}   for   salt transports (unit: $10^{9}Kg s^{-1}$) \\ Namelist variables in \ngn{namdct} control how frequently the flows are summed and the time scales over which Namelist variables in \nam{_diadct}{\_diadct} control how frequently the flows are summed and the time scales over which they are averaged, as well as the level of output for debugging: \np{nn\_dct}   : frequency of instantaneous transports computing \np{nn\_dctwri}: frequency of writing ( mean of instantaneous transports ) \np{nn\_debug} : debugging of the section \np{nn_dct}{nn\_dct}   : frequency of instantaneous transports computing \np{nn_dctwri}{nn\_dctwri}: frequency of writing ( mean of instantaneous transports ) \np{nn_debug}{nn\_debug} : debugging of the section %% ================================================================================================= \subsubsection{Creating a binary file containing the pathway of each section} In \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT/run}, In \texttt{tools/SECTIONS\_DIADCT/run}, the file \textit{ {list\_sections.ascii\_global}} contains a list of all the sections that are to be computed (this list of sections is based on MERSEA project metrics). Each section is defined by: \\ \noindent {\scriptsize \texttt{long1 lat1 long2 lat2 nclass (ok/no)strpond (no)ice section\_name}} \\ \noindent { \texttt{long1 lat1 long2 lat2 nclass (ok/no)strpond (no)ice section\_name}} \\ with: - \texttt{long2 lat2}, coordinates of the second extremity of the section; - \texttt{nclass}    the number of bounds of your classes (\eg bounds for 2 classes); - \texttt{nclass}    the number of bounds of your classes (\eg\ bounds for 2 classes); - \texttt{okstrpond} to compute    heat and       salt transports, \texttt{nostrpond} if no; \noindent If nclass $\neq$ 0, the next lines contain the class type and the nclass bounds: \\ {\scriptsize { \texttt{ long1 lat1 long2 lat2 nclass (ok/no)strpond (no)ice section\_name \\ - \texttt{zsigp} for potential density classes \\ The script \texttt{job.ksh} computes the pathway for each section and creates a binary file \texttt{section\_ijglobal.diadct\_ORCA2\_LIM} which is read by NEMO. \\ \texttt{section\_ijglobal.diadct} which is read by \NEMO. \\ It is possible to use this tools for new configuations: \texttt{job.ksh} has to be updated with and the ATL\_Cuba\_Florida with 4 temperature clases (5 class bounds), are shown: \\ \noindent {\scriptsize { \texttt{ -68.    -54.5   -60.    -64.7  00 okstrpond noice ACC\_Drake\_Passage \\ } %% ================================================================================================= \subsubsection{To read the output files} The output format is: \\ {\scriptsize { \texttt{ date, time-step number, section number,                \\ \begin{table} \scriptsize \begin{tabular}{|l|l|l|l|l|} \hline \end{table} % ================================================================ % Steric effect in sea surface height % ================================================================ %% ================================================================================================= \section{Diagnosing the steric effect in sea surface height} \label{sec:DIA_steric} Changes in steric sea level are caused when changes in the density of the water column imply an expansion or The steric effect is therefore not explicitely represented. This approximation does not represent a serious error with respect to the flow field calculated by the model \citep{Greatbatch_JGR94}, but extra attention is required when investigating sea level, \citep{greatbatch_JGR94}, but extra attention is required when investigating sea level, as steric changes are an important contribution to local changes in sea level on seasonal and climatic time scales. This is especially true for investigation into sea level rise due to global warming. Fortunately, the steric contribution to the sea level consists of a spatially uniform component that can be diagnosed by considering the mass budget of the world ocean \citep{Greatbatch_JGR94}. can be diagnosed by considering the mass budget of the world ocean \citep{greatbatch_JGR94}. In order to better understand how global mean sea level evolves and thus how the steric sea level can be diagnosed, we compare, in the following, the non-Boussinesq and Boussinesq cases. Let denote $\mathcal{M}$ the total mass    of liquid seawater ($\mathcal{M} = \int_D \rho dv$), $\mathcal{V}$ the total volume  of        seawater      ($\mathcal{V} = \int_D dv$), $\mathcal{A}$ the total surface of       the ocean      ($\mathcal{A} = \int_S ds$), $\bar{\rho}$ the global mean  seawater (\textit{in situ}) density $\mathcal{M}$ the total mass    of liquid seawater ($\mathcal{M} = \int_D \rho dv$), $\mathcal{V}$ the total volume  of        seawater      ($\mathcal{V} = \int_D dv$), $\mathcal{A}$ the total surface of       the ocean      ($\mathcal{A} = \int_S ds$), $\bar{\rho}$ the global mean  seawater (\textit{in situ}) density ($\bar{\rho} = 1/\mathcal{V} \int_D \rho \,dv$), and $\bar{\eta}$ the global mean sea level $\bar{\eta}$ the global mean sea level ($\bar{\eta} = 1/\mathcal{A} \int_S \eta \,ds$). \mathcal{V} &=  \mathcal{A}  \;\bar{\eta} \end{split} \label{eq:MV_nBq} \label{eq:DIA_MV_nBq} \frac{1}{e_3} \partial_t ( e_3\,\rho) + \nabla( \rho \, \textbf{U} ) = \left. \frac{\textit{emp}}{e_3}\right|_\textit{surface} \label{eq:Co_nBq} \label{eq:DIA_Co_nBq} where $\rho$ is the \textit{in situ} density, and \textit{emp} the surface mass exchanges with the other media of the Earth system (atmosphere, sea-ice, land). Its global averaged leads to the total mass change Its global averaged leads to the total mass change \partial_t \mathcal{M} = \mathcal{A} \;\overline{\textit{emp}} \label{eq:Mass_nBq} \label{eq:DIA_Mass_nBq} where $\overline{\textit{emp}} = \int_S \textit{emp}\,ds$ is the net mass flux through the ocean surface. Bringing \autoref{eq:Mass_nBq} and the time derivative of \autoref{eq:MV_nBq} together leads to Bringing \autoref{eq:DIA_Mass_nBq} and the time derivative of \autoref{eq:DIA_MV_nBq} together leads to the evolution equation of the mean sea level \partial_t \bar{\eta} = \frac{\overline{\textit{emp}}}{ \bar{\rho}} - \frac{\mathcal{V}}{\mathcal{A}} \;\frac{\partial_t \bar{\rho} }{\bar{\rho}} \label{eq:ssh_nBq} \label{eq:DIA_ssh_nBq} The first term in equation \autoref{eq:ssh_nBq} alters sea level by adding or subtracting mass from the ocean. The second term arises from temporal changes in the global mean density; \ie from steric effects. The first term in equation \autoref{eq:DIA_ssh_nBq} alters sea level by adding or subtracting mass from the ocean. The second term arises from temporal changes in the global mean density; \ie\ from steric effects. In a Boussinesq fluid, $\rho$ is replaced by $\rho_o$ in all the equation except when $\rho$ appears multiplied by the gravity (\ie in the hydrostatic balance of the primitive Equations). In particular, the mass conservation equation, \autoref{eq:Co_nBq}, degenerates into the incompressibility equation: the gravity (\ie\ in the hydrostatic balance of the primitive Equations). In particular, the mass conservation equation, \autoref{eq:DIA_Co_nBq}, degenerates into the incompressibility equation: $\frac{1}{e_3} \partial_t ( e_3 ) + \nabla( \textbf{U} ) = \left. \frac{\textit{emp}}{\rho_o \,e_3}\right|_ \textit{surface} % \label{eq:Co_Bq} % \label{eq:DIA_Co_Bq}$ $\partial_t \mathcal{V} = \mathcal{A} \;\frac{\overline{\textit{emp}}}{\rho_o} % \label{eq:V_Bq} % \label{eq:DIA_V_Bq}$ the ocean surface, not by changes in mean mass of the ocean: the steric effect is missing in a Boussinesq fluid. Nevertheless, following \citep{Greatbatch_JGR94}, the steric effect on the volume can be diagnosed by considering the mass budget of the ocean. Nevertheless, following \citep{greatbatch_JGR94}, the steric effect on the volume can be diagnosed by considering the mass budget of the ocean. The apparent changes in $\mathcal{M}$, mass of the ocean, which are not induced by surface mass flux must be compensated by a spatially uniform change in the mean sea level due to expansion/contraction of the ocean \citep{Greatbatch_JGR94}. \citep{greatbatch_JGR94}. In others words, the Boussinesq mass, $\mathcal{M}_o$, can be related to $\mathcal{M}$, the total mass of the ocean seen by the Boussinesq model, via the steric contribution to the sea level, \mathcal{M}_o = \mathcal{M} + \rho_o \,\eta_s \,\mathcal{A} \label{eq:M_Bq} \label{eq:DIA_M_Bq} is converted into a mean change in sea level. Introducing the total density anomaly, $\mathcal{D}= \int_D d_a \,dv$, where $d_a = (\rho -\rho_o ) / \rho_o$ is the density anomaly used in \NEMO (cf. \autoref{subsec:TRA_eos}) in \autoref{eq:M_Bq} leads to a very simple form for the steric height: where $d_a = (\rho -\rho_o ) / \rho_o$ is the density anomaly used in \NEMO\ (cf. \autoref{subsec:TRA_eos}) in \autoref{eq:DIA_M_Bq} leads to a very simple form for the steric height: \eta_s = - \frac{1}{\mathcal{A}} \mathcal{D} \label{eq:steric_Bq} \label{eq:DIA_steric_Bq} The above formulation of the steric height of a Boussinesq ocean requires four remarks. First, one can be tempted to define $\rho_o$ as the initial value of $\mathcal{M}/\mathcal{V}$, \ie set $\mathcal{D}_{t=0}=0$, so that the initial steric height is zero. \ie\ set $\mathcal{D}_{t=0}=0$, so that the initial steric height is zero. We do not recommend that. Indeed, in this case $\rho_o$ depends on the initial state of the ocean. This value is a sensible choice for the reference density used in a Boussinesq ocean climate model since, with the exception of only a small percentage of the ocean, density in the World Ocean varies by no more than 2$\%$ from this value (\cite{Gill1982}, page 47). 2$\%$ from this value (\cite{gill_bk82}, page 47). Second, we have assumed here that the total ocean surface, $\mathcal{A}$, does not change when the sea level is changing as it is the case in all global ocean GCMs (wetting and drying of grid point is not allowed). Third, the discretisation of \autoref{eq:steric_Bq} depends on the type of free surface which is considered. In the non linear free surface case, \ie \key{vvl} defined, it is given by Third, the discretisation of \autoref{eq:DIA_steric_Bq} depends on the type of free surface which is considered. In the non linear free surface case, \ie\ \np[=.true.]{ln_linssh}{ln\_linssh}, it is given by $\eta_s = - \frac{ \sum_{i,\,j,\,k} d_a\; e_{1t} e_{2t} e_{3t} }{ \sum_{i,\,j,\,k} e_{1t} e_{2t} e_{3t} } % \label{eq:discrete_steric_Bq_nfs} % \label{eq:DIA_discrete_steric_Bq_nfs}$ \eta_s = - \frac{ \sum_{i,\,j,\,k} d_a\; e_{1t}e_{2t}e_{3t} + \sum_{i,\,j} d_a\; e_{1t}e_{2t} \eta } { \sum_{i,\,j,\,k}       e_{1t}e_{2t}e_{3t} + \sum_{i,\,j}       e_{1t}e_{2t} \eta } % \label{eq:discrete_steric_Bq_fs} % \label{eq:DIA_discrete_steric_Bq_fs} \] so that there are no associated ocean currents. Hence, the dynamically relevant sea level is the effective sea level, \ie the sea level as if sea ice (and snow) were converted to liquid seawater \citep{Campin_al_OM08}. However, in the current version of \NEMO the sea-ice is levitating above the ocean without mass exchanges between \ie\ the sea level as if sea ice (and snow) were converted to liquid seawater \citep{campin.marshall.ea_OM08}. However, in the current version of \NEMO\ the sea-ice is levitating above the ocean without mass exchanges between ice and ocean. Therefore the model effective sea level is always given by $\eta + \eta_s$, whether or not there is sea ice present. $\eta_s = - \frac{1}{\mathcal{A}} \int_D d_a(T,S_o,p_o) \,dv % \label{eq:thermosteric_Bq} % \label{eq:DIA_thermosteric_Bq}$ where $S_o$ and $p_o$ are the initial salinity and pressure, respectively. Both steric and thermosteric sea level are computed in \mdl{diaar5} which needs the \key{diaar5} defined to be called. % ------------------------------------------------------------------------------------------------------------- %       Other Diagnostics % ------------------------------------------------------------------------------------------------------------- \section{Other diagnostics (\protect\key{diahth}, \protect\key{diaar5})} Both steric and thermosteric sea level are computed in \mdl{diaar5}. %% ================================================================================================= \section{Other diagnostics} \label{sec:DIA_diag_others} The available ready-to-add diagnostics modules can be found in directory DIA. \subsection{Depth of various quantities (\protect\mdl{diahth})} %% ================================================================================================= \subsection[Depth of various quantities (\textit{diahth.F90})]{Depth of various quantities (\protect\mdl{diahth})} Among the available diagnostics the following ones are obtained when defining the \key{diahth} CPP key: - the mixed layer depth (based on a density criterion \citep{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth}) - the mixed layer depth (based on a density criterion \citep{de-boyer-montegut.madec.ea_JGR04}) (\mdl{diahth}) - the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diahth}) - the depth of the thermocline (maximum of the vertical temperature gradient) (\mdl{diahth}) % ----------------------------------------------------------- %     Poleward heat and salt transports % ----------------------------------------------------------- \subsection{Poleward heat and salt transports (\protect\mdl{diaptr})} %------------------------------------------namptr----------------------------------------- \nlst{namptr} %----------------------------------------------------------------------------------------- The poleward heat and salt transports, their advective and diffusive component, and the meriodional stream function can be computed on-line in \mdl{diaptr} \np{ln\_diaptr} to true (see the \textit{\ngn{namptr} } namelist below). When \np{ln\_subbas}\forcode{ = .true.}, transports and stream function are computed for the Atlantic, Indian, \begin{figure}[!t] \centering \includegraphics[width=0.66\textwidth]{DIA_mask_subasins} \caption[Decomposition of the World Ocean to compute transports as well as the meridional stream-function]{ Decomposition of the World Ocean (here ORCA2) into sub-basin used in to compute the heat and salt transports as well as the meridional stream-function: Atlantic basin (red), Pacific basin (green), Indian basin (blue), Indo-Pacific basin (blue+green). Note that semi-enclosed seas (Red, Med and Baltic seas) as well as Hudson Bay are removed from the sub-basins. Note also that the Arctic Ocean has been split into Atlantic and Pacific basins along the North fold line. } \label{fig:DIA_mask_subasins} \end{figure} %% ================================================================================================= \subsection[CMIP specific diagnostics (\textit{diaar5.F90}, \textit{diaptr.F90})]{CMIP specific diagnostics (\protect\mdl{diaar5})} A series of diagnostics has been added in the \mdl{diaar5} and \mdl{diaptr}. In \mdl{diaar5} they correspond to outputs that are required for AR5 simulations (CMIP5) (see also \autoref{sec:DIA_steric} for one of them). The module \mdl{diaar5} is active when one of the following outputs is required : global total volume (voltot), global mean ssh (sshtot), global total mass (masstot), global mean temperature (temptot), global mean ssh steric (sshsteric), global mean ssh thermosteric (sshthster), global mean salinity (saltot), sea water pressure at sea floor (botpres), dynamic sea surface height (sshdyn). In \mdl{diaptr} when \np[=.true.]{ln_diaptr}{ln\_diaptr} (see the \nam{ptr}{ptr} namelist below) can be computed on-line the poleward heat and salt transports, their advective and diffusive component, and the meriodional stream function . When \np[=.true.]{ln_subbas}{ln\_subbas}, transports and stream function are computed for the Atlantic, Indian, Pacific and Indo-Pacific Oceans (defined north of 30\deg{S}) as well as for the World Ocean. The sub-basin decomposition requires an input file (\ifile{subbasins}) which contains three 2D mask arrays, the Indo-Pacific mask been deduced from the sum of the Indian and Pacific mask (\autoref{fig:mask_subasins}). %>>>>>>>>>>>>>>>>>>>>>>>>>>>> \begin{figure}[!t] \begin{center} \includegraphics[width=1.0\textwidth]{Fig_mask_subasins} \caption{ \protect\label{fig:mask_subasins} Decomposition of the World Ocean (here ORCA2) into sub-basin used in to compute the heat and salt transports as well as the meridional stream-function: Atlantic basin (red), Pacific basin (green), Indian basin (bleue), Indo-Pacific basin (bleue+green). Note that semi-enclosed seas (Red, Med and Baltic seas) as well as Hudson Bay are removed from the sub-basins. Note also that the Arctic Ocean has been split into Atlantic and Pacific basins along the North fold line. } \end{center} \end{figure} %>>>>>>>>>>>>>>>>>>>>>>>>>>>> % ----------------------------------------------------------- %       CMIP specific diagnostics % ----------------------------------------------------------- \subsection{CMIP specific diagnostics (\protect\mdl{diaar5})} A series of diagnostics has been added in the \mdl{diaar5}. They corresponds to outputs that are required for AR5 simulations (CMIP5) (see also \autoref{sec:DIA_steric} for one of them). Activating those outputs requires to define the \key{diaar5} CPP key. % ----------------------------------------------------------- %       25 hour mean and hourly Surface, Mid and Bed % ----------------------------------------------------------- the Indo-Pacific mask been deduced from the sum of the Indian and Pacific mask (\autoref{fig:DIA_mask_subasins}). \begin{listing} \nlst{namptr} \caption{\forcode{&namptr}} \label{lst:namptr} \end{listing} %% ================================================================================================= \subsection{25 hour mean output for tidal models} %------------------------------------------nam_dia25h------------------------------------- \nlst{nam_dia25h} %----------------------------------------------------------------------------------------- \begin{listing} \nlst{nam_dia25h} \caption{\forcode{&nam_dia25h}} \label{lst:nam_dia25h} \end{listing} A module is available to compute a crudely detided M2 signal by obtaining a 25 hour mean. This diagnostic is actived with the logical $ln\_dia25h$. % ----------------------------------------------------------- %     Top Middle and Bed hourly output % ----------------------------------------------------------- %% ================================================================================================= \subsection{Top middle and bed hourly output} %------------------------------------------nam_diatmb----------------------------------------------------- \nlst{nam_diatmb} %---------------------------------------------------------------------------------------------------------- \begin{listing} \nlst{nam_diatmb} \caption{\forcode{&nam_diatmb}} \label{lst:nam_diatmb} \end{listing} A module is available to output the surface (top), mid water and bed diagnostics of a set of standard variables. This diagnostic is actived with the logical $ln\_diatmb$. % ----------------------------------------------------------- %     Courant numbers % ----------------------------------------------------------- %% ================================================================================================= \subsection{Courant numbers} $C_u = |u|\frac{\rdt}{e_{1u}}, \quad C_v = |v|\frac{\rdt}{e_{2v}}, \quad C_w = |w|\frac{\rdt}{e_{3w}} % \label{eq:CFL} % \label{eq:DIA_CFL}$ Values greater than 1 indicate that information is propagated across more than one grid cell in a single time step. The variables can be activated by setting the \np{nn\_diacfl} namelist parameter to 1 in the \ngn{namctl} namelist. The variables can be activated by setting the \np{nn_diacfl}{nn\_diacfl} namelist parameter to 1 in the \nam{ctl}{ctl} namelist. The diagnostics will be written out to an ascii file named cfl\_diagnostics.ascii. In this file the maximum value of $C_u$, $C_v$, and $C_w$ are printed at each timestep along with the coordinates of At the end of the model run the maximum value of $C_u$, $C_v$, and $C_w$ for the whole model run is printed along with the coordinates of each. The maximum values from the run are also copied to the ocean.output file. % ================================================================ \biblio \pindex The maximum values from the run are also copied to the ocean.output file. \subinc{\input{../../global/epilogue}} \end{document}
 r10442 \begin{document} % ================================================================ % Diurnal SST models (DIU) % Edited by James While % ================================================================ \chapter{Diurnal SST Models (DIU)} \label{chap:DIU} \minitoc \thispagestyle{plain} \chaptertoc \newpage $\$\newline % force a new line \paragraph{Changes record} ~\\ {\footnotesize \begin{tabularx}{\textwidth}{l||X|X} Release & Author(s) & Modifications \\ \hline {\em   4.0} & {\em ...} & {\em ...} \\ {\em   3.6} & {\em ...} & {\em ...} \\ {\em   3.4} & {\em ...} & {\em ...} \\ {\em <=3.4} & {\em ...} & {\em ...} \end{tabularx} } \clearpage Code to produce an estimate of the diurnal warming and cooling of the sea surface skin temperature (skin SST) is found in the DIU directory. temperature (skin SST) is found in the DIU directory. The skin temperature can be split into three parts: \begin{itemize} \item A foundation SST which is free from diurnal warming. \item A warm layer, typically ~3\,m thick, \item A foundation SST which is free from diurnal warming. \item A warm layer, typically ~3\,m thick, where heating from solar radiation can cause a warm stably stratified layer during the daytime \item A cool skin, a thin layer, approximately ~1\, mm thick, \item A cool skin, a thin layer, approximately ~1\, mm thick, where long wave cooling is dominant and cools the immediate ocean surface. \end{itemize} Models are provided for both the warm layer, \mdl{diurnal\_bulk}, and the cool skin, \mdl{cool\_skin}. Foundation SST is not considered as it can be obtained either from the main NEMO model (\ie from the temperature of the top few model levels) or from some other source. Foundation SST is not considered as it can be obtained either from the main \NEMO\ model (\ie\ from the temperature of the top few model levels) or from some other source. It must be noted that both the cool skin and warm layer models produce estimates of the change in temperature ($\Delta T_{\rm{cs}}$ and $\Delta T_{\rm{wl}}$) and ($\Delta T_{\mathrm{cs}}$ and $\Delta T_{\mathrm{wl}}$) and both must be added to a foundation SST to obtain the true skin temperature. Both the cool skin and warm layer models are controlled through the namelist \ngn{namdiu}: Both the cool skin and warm layer models are controlled through the namelist \nam{diu}{diu}: \nlst{namdiu} \begin{listing} \nlst{namdiu} \caption{\forcode{&namdiu}} \label{lst:namdiu} \end{listing} This namelist contains only two variables: \begin{description} \item[\np{ln\_diurnal}] A logical switch for turning on/off both the cool skin and warm layer. \item[\np{ln\_diurnal\_only}] A logical switch which if \forcode{.true.} will run the diurnal model without the other dynamical parts of NEMO. \np{ln\_diurnal\_only} must be \forcode{.false.} if \np{ln\_diurnal} is \forcode{.false.}. \item [{\np{ln_diurnal}{ln\_diurnal}}] A logical switch for turning on/off both the cool skin and warm layer. \item [{\np{ln_diurnal_only}{ln\_diurnal\_only}}] A logical switch which if \forcode{.true.} will run the diurnal model without the other dynamical parts of \NEMO. \np{ln_diurnal_only}{ln\_diurnal\_only} must be \forcode{.false.} if \np{ln_diurnal}{ln\_diurnal} is \forcode{.false.}. \end{description} Initialisation is through the restart file. Specifically the code will expect the presence of the 2-D variable Dsst'' to initialise the warm layer. The cool skin model, which is determined purely by the instantaneous fluxes, has no initialisation variable. The cool skin model, which is determined purely by the instantaneous fluxes, has no initialisation variable. %=============================================================== %% ================================================================================================= \section{Warm layer model} \label{sec:warm_layer_sec} %=============================================================== \label{sec:DIU_warm_layer_sec} The warm layer is calculated using the model of \citet{Takaya_al_JGR10} (TAKAYA10 model hereafter). The warm layer is calculated using the model of \citet{takaya.bidlot.ea_JGR10} (TAKAYA10 model hereafter). This is a simple flux based model that is defined by the equations \begin{align} \frac{\partial{\Delta T_{\rm{wl}}}}{\partial{t}}&=&\frac{Q(\nu+1)}{D_T\rho_w c_p \frac{\partial{\Delta T_{\mathrm{wl}}}}{\partial{t}}&=&\frac{Q(\nu+1)}{D_T\rho_w c_p \nu}-\frac{(\nu+1)ku^*_{w}f(L_a)\Delta T}{D_T\Phi\!\left(\frac{D_T}{L}\right)} \mbox{,} \label{eq:ecmwf1} \\ L&=&\frac{\rho_w c_p u^{*^3}_{w}}{\kappa g \alpha_w Q }\mbox{,}\label{eq:ecmwf2} \label{eq:DIU_ecmwf1} \\ L&=&\frac{\rho_w c_p u^{*^3}_{w}}{\kappa g \alpha_w Q }\mbox{,}\label{eq:DIU_ecmwf2} \end{align} where $\Delta T_{\rm{wl}}$ is the temperature difference between the top of the warm layer and the depth $D_T=3$\,m at which there is assumed to be no diurnal signal. In equation (\autoref{eq:ecmwf1}) $\alpha_w=2\times10^{-4}$ is the thermal expansion coefficient of water, where $\Delta T_{\mathrm{wl}}$ is the temperature difference between the top of the warm layer and the depth $D_T=3$\,m at which there is assumed to be no diurnal signal. In equation (\autoref{eq:DIU_ecmwf1}) $\alpha_w=2\times10^{-4}$ is the thermal expansion coefficient of water, $\kappa=0.4$ is von K\'{a}rm\'{a}n's constant, $c_p$ is the heat capacity at constant pressure of sea water, $\rho_w$ is the water density, and $L$ is the Monin-Obukhov length. The tunable variable $\nu$ is a shape parameter that defines the expected subskin temperature profile via $T(z) = T(0) - \left( \frac{z}{D_T} \right)^\nu \Delta T_{\rm{wl}}$, $T(z) = T(0) - \left( \frac{z}{D_T} \right)^\nu \Delta T_{\mathrm{wl}}$, where $T$ is the absolute temperature and $z\le D_T$ is the depth below the top of the warm layer. The influence of wind on TAKAYA10 comes through the magnitude of the friction velocity of the water $u^*_{w}$, the relationship $u^*_{w} = u_{10}\sqrt{\frac{C_d\rho_a}{\rho_w}}$, where $C_d$ is the drag coefficient, and $\rho_a$ is the density of air. The symbol $Q$ in equation (\autoref{eq:ecmwf1}) is the instantaneous total thermal energy flux into The symbol $Q$ in equation (\autoref{eq:DIU_ecmwf1}) is the instantaneous total thermal energy flux into the diurnal layer, \ie $Q = Q_{\rm{sol}} + Q_{\rm{lw}} + Q_{\rm{h}}\mbox{,} % \label{eq:e_flux_eqn} Q = Q_{\mathrm{sol}} + Q_{\mathrm{lw}} + Q_{\mathrm{h}}\mbox{,} % \label{eq:DIU_e_flux_eqn}$ where $Q_{\rm{h}}$ is the sensible and latent heat flux, $Q_{\rm{lw}}$ is the long wave flux, and $Q_{\rm{sol}}$ is the solar flux absorbed within the diurnal warm layer. For $Q_{\rm{sol}}$ the 9 term representation of \citet{Gentemann_al_JGR09} is used. In equation \autoref{eq:ecmwf1} the function $f(L_a)=\max(1,L_a^{\frac{2}{3}})$, where $Q_{\mathrm{h}}$ is the sensible and latent heat flux, $Q_{\mathrm{lw}}$ is the long wave flux, and $Q_{\mathrm{sol}}$ is the solar flux absorbed within the diurnal warm layer. For $Q_{\mathrm{sol}}$ the 9 term representation of \citet{gentemann.minnett.ea_JGR09} is used. In equation \autoref{eq:DIU_ecmwf1} the function $f(L_a)=\max(1,L_a^{\frac{2}{3}})$, where $L_a=0.3$\footnote{ This is a global average value, more accurately $L_a$ could be computed as $L_a=(u^*_{w}/u_s)^{\frac{1}{2}}$, 4\zeta^2}{1+3\zeta+0.25\zeta^2} &(\zeta \ge 0) \\ (1 - 16\zeta)^{-\frac{1}{2}} & (\zeta < 0) \mbox{,} \end{array} \right. \label{eq:stab_func_eqn} \end{array} \right. \label{eq:DIU_stab_func_eqn} where $\zeta=\frac{D_T}{L}$.  It is clear that the first derivative of (\autoref{eq:stab_func_eqn}), and thus of (\autoref{eq:ecmwf1}), is discontinuous at $\zeta=0$ (\ie $Q\rightarrow0$ in equation (\autoref{eq:ecmwf2})). where $\zeta=\frac{D_T}{L}$.  It is clear that the first derivative of (\autoref{eq:DIU_stab_func_eqn}), and thus of (\autoref{eq:DIU_ecmwf1}), is discontinuous at $\zeta=0$ (\ie\ $Q\rightarrow0$ in equation (\autoref{eq:DIU_ecmwf2})). The two terms on the right hand side of (\autoref{eq:ecmwf1}) represent different processes. The two terms on the right hand side of (\autoref{eq:DIU_ecmwf1}) represent different processes. The first term is simply the diabatic heating or cooling of the diurnal warm layer due to thermal energy fluxes into and out of the layer. In practice the second term acts as a relaxation on the temperature. %=============================================================== %% ================================================================================================= \section{Cool skin model} \label{sec:DIU_cool_skin_sec} \section{Cool skin model} \label{sec:cool_skin_sec} %=============================================================== The cool skin is modelled using the framework of \citet{Saunders_JAS82} who used a formulation of the near surface temperature difference based upon the heat flux and the friction velocity $u^*_{w}$. As the cool skin is so thin (~1\,mm) we ignore the solar flux component to the heat flux and the Saunders equation for the cool skin temperature difference $\Delta T_{\rm{cs}}$ becomes The cool skin is modelled using the framework of \citet{saunders_JAS67} who used a formulation of the near surface temperature difference based upon the heat flux and the friction velocity $u^*_{w}$. As the cool skin is so thin (~1\,mm) we ignore the solar flux component to the heat flux and the Saunders equation for the cool skin temperature difference $\Delta T_{\mathrm{cs}}$ becomes $% \label{eq:sunders_eqn} \Delta T_{\rm{cs}}=\frac{Q_{\rm{ns}}\delta}{k_t} \mbox{,} % \label{eq:DIU_sunders_eqn} \Delta T_{\mathrm{cs}}=\frac{Q_{\mathrm{ns}}\delta}{k_t} \mbox{,}$ where $Q_{\rm{ns}}$ is the, usually negative, non-solar heat flux into the ocean and where $Q_{\mathrm{ns}}$ is the, usually negative, non-solar heat flux into the ocean and $k_t$ is the thermal conductivity of sea water. $\delta$ is the thickness of the skin layer and is given by \label{eq:sunders_thick_eqn} \label{eq:DIU_sunders_thick_eqn} \delta=\frac{\lambda \mu}{u^*_{w}} \mbox{,} where $\mu$ is the kinematic viscosity of sea water and $\lambda$ is a constant of proportionality which \citet{Saunders_JAS82} suggested varied between 5 and 10. \citet{saunders_JAS67} suggested varied between 5 and 10. The value of $\lambda$ used in equation (\autoref{eq:sunders_thick_eqn}) is that of \citet{Artale_al_JGR02}, which is shown in \citet{Tu_Tsuang_GRL05} to outperform a number of other parametrisations at The value of $\lambda$ used in equation (\autoref{eq:DIU_sunders_thick_eqn}) is that of \citet{artale.iudicone.ea_JGR02}, which is shown in \citet{tu.tsuang_GRL05} to outperform a number of other parametrisations at both low and high wind speeds. Specifically, $% \label{eq:artale_lambda_eqn} % \label{eq:DIU_artale_lambda_eqn} \lambda = \frac{ 8.64\times10^4 u^*_{w} k_t }{ \rho c_p h \mu \gamma }\mbox{,}$ $\gamma$ is a dimensionless function of wind speed $u$: $% \label{eq:artale_gamma_eqn} % \label{eq:DIU_artale_gamma_eqn} \gamma = \begin{cases}$ \biblio \pindex \subinc{\input{../../global/epilogue}} \end{document}
 r10502 \begin{document} % ================================================================ % Chapter 2 ——— Space and Time Domain (DOM) % ================================================================ \chapter{Space Domain (DOM)} \label{chap:DOM} \minitoc % Missing things: %  - istate: description of the initial state   ==> this has to be put elsewhere.. %                  perhaps in MISC ?  By the way the initialisation of T S and dynamics %                  should be put outside of DOM routine (better with TRC staff and off-line %                  tracers) %  -geo2ocean:  how to switch from geographic to mesh coordinate %     - domclo:  closed sea and lakes.... management of closea sea area : specific to global configuration, both forced and coupled \newpage Having defined the continuous equations in \autoref{chap:PE} and chosen a time discretization \autoref{chap:STP}, we need to choose a discretization on a grid, and numerical algorithms. % Missing things % -    istate: description of the initial state   ==> this has to be put elsewhere.. %              perhaps in MISC ?  By the way the initialisation of T S and dynamics %              should be put outside of DOM routine (better with TRC staff and off-line %              tracers) % - geo2ocean: how to switch from geographic to mesh coordinate % -    domclo: closed sea and lakes.... %              management of closea sea area: specific to global cfg, both forced and coupled \thispagestyle{plain} \chaptertoc \paragraph{Changes record} ~\\ {\footnotesize \begin{tabularx}{0.8\textwidth}{l||X|X} Release                                                                                 & Author(s)                                                                               & Modifications                                                                           \\ \hline {\em 4.0                                                                              } & {\em Simon M\"{u}ller \& Andrew Coward \newline \newline Simona Flavoni and Tim Graham                                                       } & {\em Compatibility changes: many options moved to external domain configuration tools (see \autoref{apdx:DOMCFG}). \newline Updates                                                                             } \\ {\em 3.6                                                                              } & {\em Rachid Benshila, Christian \'{E}th\'{e}, Pierre Mathiot and Gurvan Madec         } & {\em Updates                                                                          } \\ {\em $\leq$ 3.4                                                                       } & {\em Gurvan Madec and S\'{e}bastien Masson                                            } & {\em First version                                                                    } \end{tabularx} } \clearpage Having defined the continuous equations in \autoref{chap:MB} and chosen a time discretisation \autoref{chap:TD}, we need to choose a grid for spatial discretisation and related numerical algorithms. In the present chapter, we provide a general description of the staggered grid used in \NEMO, and other information relevant to the main directory routines as well as the DOM (DOMain) directory. % ================================================================ % Fundamentals of the Discretisation % ================================================================ and other relevant information about the DOM (DOMain) source code modules. %% ================================================================================================= \section{Fundamentals of the discretisation} \label{sec:DOM_basics} % ------------------------------------------------------------------------------------------------------------- %        Arrangement of Variables % ------------------------------------------------------------------------------------------------------------- %% ================================================================================================= \subsection{Arrangement of variables} \label{subsec:DOM_cell} %>>>>>>>>>>>>>>>>>>>>>>>>>>>> \begin{figure}[!tb] \begin{center} \includegraphics[]{Fig_cell} \caption{ \protect\label{fig:cell} Arrangement of variables. $t$ indicates scalar points where temperature, salinity, density, pressure and horizontal divergence are defined. $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where both relative and planetary vorticities are defined. } \end{center} \begin{figure} \centering \includegraphics[width=0.33\textwidth]{DOM_cell} \caption[Arrangement of variables in the unit cell of space domain]{ Arrangement of variables in the unit cell of space domain. $t$ indicates scalar points where temperature, salinity, density, pressure and horizontal divergence are defined. $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where both relative and planetary vorticities are defined.} \label{fig:DOM_cell} \end{figure} %>>>>>>>>>>>>>>>>>>>>>>>>>>>> The numerical techniques used to solve the Primitive Equations in this model are based on the traditional, centred second-order finite difference approximation. Special attention has been given to the homogeneity of the solution in the three space directions. The numerical techniques used to solve the Primitive Equations in this model are based on the traditional, centred second-order finite difference approximation. Special attention has been given to the homogeneity of the solution in the three spatial directions. The arrangement of variables is the same in all directions. It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in the centre of each face of the cells (\autoref{fig:cell}). This is the generalisation to three dimensions of the well-known C'' grid in Arakawa's classification \citep{Mesinger_Arakawa_Bk76}. The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points. The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:cell}. In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where the scale factors are defined. Each scale factor is defined as the local analytical value provided by \autoref{eq:scale_factors}. It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in the centre of each face of the cells (\autoref{fig:DOM_cell}). This is the generalisation to three dimensions of the well-known C'' grid in Arakawa's classification \citep{mesinger.arakawa_bk76}. The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points. The ocean mesh (\ie\ the position of all the scalar and vector points) is defined by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:DOM_cell}. In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where the scale factors are defined. Each scale factor is defined as the local analytical value provided by \autoref{eq:MB_scale_factors}. As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and $\pd[]{z}$ are evaluated in a uniform mesh with a grid size of unity. Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation while the scale factors are chosen equal to their local analytical value. $\pd[]{z}$ are evaluated is a uniform mesh with a grid size of unity. Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation while the scale factors are chosen equal to their local analytical value. An important point here is that the partial derivative of the scale factors must be evaluated by centred finite difference approximation, not from their analytical expression. This preserves the symmetry of the discrete set of equations and therefore satisfies many of the continuous properties (see \autoref{apdx:C}). This preserves the symmetry of the discrete set of equations and therefore satisfies many of the continuous properties (see \autoref{apdx:INVARIANTS}). A similar, related remark can be made about the domain size: when needed, an area, volume, or the total ocean depth must be evaluated as the sum of the relevant scale factors (see \autoref{eq:DOM_bar} in the next section). %>>>>>>>>>>>>>>>>>>>>>>>>>>>> \begin{table}[!tb] \begin{center} \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|} \hline T  & $i$ & $j$ & $k$ \\ \hline u  & $i + 1/2$ & $j$ & $k$ \\ \hline v  & $i$ & $j + 1/2$ & $k$ \\ \hline w  & $i$ & $j$ & $k + 1/2$ \\ \hline f  & $i + 1/2$ & $j + 1/2$ & $k$ \\ \hline uw & $i + 1/2$ & $j$ & $k + 1/2$ \\ \hline vw & $i$ & $j + 1/2$ & $k + 1/2$ \\ \hline fw & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\ \hline \end{tabular} \caption{ \protect\label{tab:cell} Location of grid-points as a function of integer or integer and a half value of the column, line or level. This indexing is only used for the writing of the semi -discrete equation. In the code, the indexing uses integer values only and has a reverse direction in the vertical (see \autoref{subsec:DOM_Num_Index}) } \end{center} when needed, an area, volume, or the total ocean depth must be evaluated as the product or sum of the relevant scale factors (see \autoref{eq:DOM_bar} in the next section). \begin{table} \centering \begin{tabular}{|l|l|l|l|} \hline t   & $i$ & $j$ & $k$ \\ \hline u   & $i + 1/2$ & $j$ & $k$ \\ \hline v   & $i$ & $j + 1/2$ & $k$ \\ \hline w   & $i$ & $j$ & $k + 1/2$ \\ \hline f   & $i + 1/2$ & $j + 1/2$ & $k$ \\ \hline uw  & $i + 1/2$ & $j$ & $k + 1/2$ \\ \hline vw  & $i$ & $j + 1/2$ & $k + 1/2$ \\ \hline fw  & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\ \hline \end{tabular} \caption[Location of grid-points]{ Location of grid-points as a function of integer or integer and a half value of the column, line or level. This indexing is only used for the writing of the semi-discrete equations. In the code, the indexing uses integer values only and is positive downwards in the vertical with $k=1$ at the surface. (see \autoref{subsec:DOM_Num_Index})} \label{tab:DOM_cell} \end{table} %>>>>>>>>>>>>>>>>>>>>>>>>>>>> % ------------------------------------------------------------------------------------------------------------- %        Vector Invariant Formulation % ------------------------------------------------------------------------------------------------------------- Note that the definition of the scale factors (\ie\ as the analytical first derivative of the transformation that results in $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) is specific to the \NEMO\ model \citep{marti.madec.ea_JGR92}. As an example, a scale factor in the $i$ direction is defined locally at a $t$-point, whereas many other models on a C grid choose to define such a scale factor as the distance between the $u$-points on each side of the $t$-point. Relying on an analytical transformation has two advantages: firstly, there is no ambiguity in the scale factors appearing in the discrete equations, since they are first introduced in the continuous equations; secondly, analytical transformations encourage good practice by the definition of smoothly varying grids (rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{treguier.dukowicz.ea_JGR96}. An example of the effect of such a choice is shown in \autoref{fig:DOM_zgr_e3}. \begin{figure} \centering \includegraphics[width=0.5\textwidth]{DOM_zgr_e3} \caption[Comparison of grid-point position, vertical grid-size and scale factors]{ Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, and (b) analytically derived grid-point position and scale factors. For both grids here, the same $w$-point depth has been chosen but in (a) the $t$-points are set half way between $w$-points while in (b) they are defined from an analytical function: $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$. Note the resulting difference between the value of the grid-size $\Delta_k$ and those of the scale factor $e_k$.} \label{fig:DOM_zgr_e3} \end{figure} %% ================================================================================================= \subsection{Discrete operators} \label{subsec:DOM_operators} Given the values of a variable $q$ at adjacent points, the differencing and averaging operators at the midpoint between them are: Given the values of a variable $q$ at adjacent points, the differencing and averaging operators at the midpoint between them are: \begin{alignat*}{2} % \label{eq:di_mi} % \label{eq:DOM_di_mi} \delta_i [q]      &= &       &q (i + 1/2) - q (i - 1/2) \\ \overline q^{\, i} &= &\big\{ &q (i + 1/2) + q (i - 1/2) \big\} / 2 Similar operators are defined with respect to $i + 1/2$, $j$, $j + 1/2$, $k$, and $k + 1/2$. Following \autoref{eq:PE_grad} and \autoref{eq:PE_lap}, the gradient of a variable $q$ defined at a $t$-point has its three components defined at $u$-, $v$- and $w$-points while its Laplacian is defined at $t$-point. Following \autoref{eq:MB_grad} and \autoref{eq:MB_lap}, the gradient of a variable $q$ defined at a $t$-point has its three components defined at $u$-, $v$- and $w$-points while its Laplacian is defined at the $t$-point. These operators have the following discrete forms in the curvilinear $s$-coordinates system: $\begin{gather*} % \label{eq:DOM_grad} \nabla q \equiv \frac{1}{e_{1u}} \delta_{i + 1/2} [q] \; \, \vect i + \frac{1}{e_{2v}} \delta_{j + 1/2} [q] \; \, \vect j + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k$ \begin{multline*} + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k \\ % \label{eq:DOM_lap} \Delta q \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} \; \lt[   \delta_i \lt( \frac{e_{2u} \, e_{3u}}{e_{1u}} \; \delta_{i + 1/2} [q] \rt) + \delta_j \lt( \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [q] \rt) \; \rt] \\ + \delta_j \lt( \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [q] \rt) \; \rt] + \frac{1}{e_{3t}} \delta_k \lt[ \frac{1              }{e_{3w}} \; \delta_{k + 1/2} [q] \rt] \end{multline*} Following \autoref{eq:PE_curl} and \autoref{eq:PE_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and \end{gather*} Following \autoref{eq:MB_curl} and \autoref{eq:MB_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and its divergence defined at $t$-points: \begin{multline} \begin{multline*} % \label{eq:DOM_curl} \nabla \times \vect A \equiv   \frac{1}{e_{2v} \, e_{3vw}} \Big[   \delta_{i + 1/2} (e_{2v} \, a_2) - \delta_{j + 1/2} (e_{1u} \, a_1) \Big] \vect k \end{multline} \begin{equation} \end{multline*} $% \label{eq:DOM_div} \nabla \cdot \vect A \equiv \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} \Big[ \delta_i (e_{2u} \, e_{3u} \, a_1) + \delta_j (e_{1v} \, e_{3v} \, a_2) \Big] + \frac{1}{e_{3t}} \delta_k (a_3) \end{equation} The vertical average over the whole water column denoted by an overbar becomes for a quantity q which is a masked field (i.e. equal to zero inside solid area):$ The vertical average over the whole water column is denoted by an overbar and is for a masked field $q$ (\ie\ a quantity that is equal to zero inside solid areas): \label{eq:DOM_bar} where $H_q$  is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points, $k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $k^o$ refers to a summation over all grid points of the same type in the direction indicated by the subscript (here $k$). $k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $\sum \limits_k$ refers to a summation over all grid points of the same type in the direction indicated by the subscript (here $k$). In continuous form, the following properties are satisfied: \end{gather} It is straightforward to demonstrate that these properties are verified locally in discrete form as soon as the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at It is straightforward to demonstrate that these properties are verified locally in discrete form as soon as the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at vector points $(u,v,w)$. Let $a$ and $b$ be two fields defined on the mesh, with value zero inside continental area. Using integration by parts it can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) are skew-symmetric linear operators, and further that the averaging operators $\overline{\cdots}^{\, i}$, $\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie \begin{alignat}{4} Let $a$ and $b$ be two fields defined on the mesh, with a value of zero inside continental areas. It can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) are skew-symmetric linear operators, and further that the averaging operators ($\overline{\cdots}^{\, i}$, $\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie \begin{alignat}{5} \label{eq:DOM_di_adj} &\sum \limits_i a_i \; \delta_i [b]      &\equiv &- &&\sum \limits_i \delta      _{   i + 1/2} [a] &b_{i + 1/2} \\ \end{alignat} In other words, the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and In other words, the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and $(\overline{\cdots}^{\, i})^* = \overline{\cdots}^{\, i + 1/2}$, respectively. These two properties will be used extensively in the \autoref{apdx:C} to These two properties will be used extensively in the \autoref{apdx:INVARIANTS} to demonstrate integral conservative properties of the discrete formulation chosen. % ------------------------------------------------------------------------------------------------------------- %        Numerical Indexing % ------------------------------------------------------------------------------------------------------------- %% ================================================================================================= \subsection{Numerical indexing} \label{subsec:DOM_Num_Index} %>>>>>>>>>>>>>>>>>>>>>>>>>>>> \begin{figure}[!tb] \begin{center} \includegraphics[]{Fig_index_hor} \caption{ \protect\label{fig:index_hor} Horizontal integer indexing used in the \fortran code. The dashed area indicates the cell in which variables contained in arrays have the same $i$- and $j$-indices } \end{center} \begin{figure} \centering \includegraphics[width=0.33\textwidth]{DOM_index_hor} \caption[Horizontal integer indexing]{ Horizontal integer indexing used in the \fortran\ code. The dashed area indicates the cell in which variables contained in arrays have the same $i$- and $j$-indices} \label{fig:DOM_index_hor} \end{figure} %>>>>>>>>>>>>>>>>>>>>>>>>>>>> The array representation used in the \fortran code requires an integer indexing while the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of integer values for $t$-points and both integer and integer and a half values for all the other points. Therefore a specific integer indexing must be defined for points other than $t$-points (\ie velocity and vorticity grid-points). Furthermore, the direction of the vertical indexing has been changed so that the surface level is at $k = 1$. % ----------------------------------- %        Horizontal Indexing % ----------------------------------- The array representation used in the \fortran\ code requires an integer indexing. However, the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of integer values for $t$-points only while all the other points involve integer and a half values. Therefore, a specific integer indexing has been defined for points other than $t$-points (\ie\ velocity and vorticity grid-points). Furthermore, the direction of the vertical indexing has been reversed and the surface level set at $k = 1$. %% ================================================================================================= \subsubsection{Horizontal indexing} \label{subsec:DOM_Num_Index_hor} The indexing in the horizontal plane has been chosen as shown in \autoref{fig:index_hor}. The indexing in the horizontal plane has been chosen as shown in \autoref{fig:DOM_index_hor}. For an increasing $i$ index ($j$ index), the $t$-point and the eastward $u$-point (northward $v$-point) have the same index (see the dashed area in \autoref{fig:index_hor}). A $t$-point and its nearest northeast $f$-point have the same $i$-and $j$-indices. % ----------------------------------- %        Vertical indexing % ----------------------------------- (see the dashed area in \autoref{fig:DOM_index_hor}). A $t$-point and its nearest north-east $f$-point have the same $i$-and $j$-indices. %% ================================================================================================= \subsubsection{Vertical indexing} \label{subsec:DOM_Num_Index_vertical} In the vertical, the chosen indexing requires special attention since the $k$-axis is re-orientated downward in the \fortran code compared to the indexing used in the semi -discrete equations and given in \autoref{subsec:DOM_cell}. The sea surface corresponds to the $w$-level $k = 1$ which is the same index as $t$-level just below (\autoref{fig:index_vert}). The last $w$-level ($k = jpk$) either corresponds to the ocean floor or is inside the bathymetry while the last $t$-level is always inside the bathymetry (\autoref{fig:index_vert}). Note that for an increasing $k$ index, a $w$-point and the $t$-point just below have the same $k$ index, in opposition to what is done in the horizontal plane where it is the $t$-point and the nearest velocity points in the direction of the horizontal axis that have the same $i$ or $j$ index (compare the dashed area in \autoref{fig:index_hor} and \autoref{fig:index_vert}). In the vertical, the chosen indexing requires special attention since the direction of the $k$-axis in the \fortran\ code is the reverse of that used in the semi-discrete equations and given in \autoref{subsec:DOM_cell}. The sea surface corresponds to the $w$-level $k = 1$, which is the same index as the $t$-level just below (\autoref{fig:DOM_index_vert}). The last $w$-level ($k = jpk$) either corresponds to or is below the ocean floor while the last $t$-level is always outside the ocean domain (\autoref{fig:DOM_index_vert}). Note that a $w$-point and the directly underlaying $t$-point have a common $k$ index (\ie\ $t$-points and their nearest $w$-point neighbour in negative index direction), in contrast to the indexing on the horizontal plane where the $t$-point has the same index as the nearest velocity points in the positive direction of the respective horizontal axis index (compare the dashed area in \autoref{fig:DOM_index_hor} and \autoref{fig:DOM_index_vert}). Since the scale factors are chosen to be strictly positive, a \textit{minus sign} appears in the \fortran code \textit{before all the vertical derivatives} of the discrete equations given in this documentation. %>>>>>>>>>>>>>>>>>>>>>>>>>>>> \begin{figure}[!pt] \begin{center} \includegraphics[]{Fig_index_vert} \caption{ \protect\label{fig:index_vert} Vertical integer indexing used in the \fortran code. Note that the $k$-axis is orientated downward. The dashed area indicates the cell in which variables contained in arrays have the same $k$-index. } \end{center} a \textit{minus sign} is included in the \fortran\ implementations of \textit{all the vertical derivatives} of the discrete equations given in this manual in order to accommodate the opposing vertical index directions in implementation and documentation. \begin{figure} \centering \includegraphics[width=0.33\textwidth]{DOM_index_vert} \caption[Vertical integer indexing]{ Vertical integer indexing used in the \fortran\ code. Note that the $k$-axis is oriented downward. The dashed area indicates the cell in which variables contained in arrays have a common $k$-index.} \label{fig:DOM_index_vert} \end{figure} %>>>>>>>>>>>>>>>>>>>>>>>>>>>> % ----------------------------------- %        Domain Size % ----------------------------------- \subsubsection{Domain size} %% ================================================================================================= \section{Spatial domain configuration} \label{subsec:DOM_config} Two typical methods are available to specify the spatial domain configuration; they can be selected using parameter \np{ln_read_cfg}{ln\_read\_cfg} parameter in namelist \nam{cfg}{cfg}. If \np{ln_read_cfg}{ln\_read\_cfg} is set to \forcode{.true.}, the domain-specific parameters and fields are read from a NetCDF input file, whose name (without its .nc suffix) can be specified as the value of the \np{cn_domcfg}{cn\_domcfg} parameter in namelist \nam{cfg}{cfg}. If \np{ln_read_cfg}{ln\_read\_cfg} is set to \forcode{.false.}, the domain-specific parameters and fields can be provided (\eg\ analytically computed) by subroutines \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. These subroutines can be supplied in the \path{MY_SRC} directory of the configuration, and default versions that configure the spatial domain for the GYRE reference configuration are present in the \path{./src/OCE/USR} directory. In version 4.0 there are no longer any options for reading complex bathymetries and performing a vertical discretisation at run-time. Whilst it is occasionally convenient to have a common bathymetry file and, for example, to run similar models with and without partial bottom boxes and/or sigma-coordinates, supporting such choices leads to overly complex code. Worse still is the difficulty of ensuring the model configurations intended to be identical are indeed so when the model domain itself can be altered by runtime selections. The code previously used to perform vertical discretisation has been incorporated into an external tool (\path{./tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMCFG}. The next subsections summarise the parameter and fields related to the configuration of the whole model domain. These represent the minimum information that must be provided either via the \np{cn_domcfg}{cn\_domcfg} file or set by code inserted into user-supplied versions of the \texttt{usrdef\_*} subroutines. The requirements are presented in three sections: the domain size (\autoref{subsec:DOM_size}), the horizontal mesh (\autoref{subsec:DOM_hgr}), and the vertical grid (\autoref{subsec:DOM_zgr}). %% ================================================================================================= \subsection{Domain size} \label{subsec:DOM_size} The total size of the computational domain is set by the parameters \np{jpiglo}, \np{jpjglo} and \np{jpkglo} in the $i$, $j$ and $k$ directions respectively. Parameters $jpi$ and $jpj$ refer to the size of each processor subdomain when the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined, see \autoref{sec:LBC_mpp}). % ================================================================ % Domain: List of fields needed % ================================================================ \section{Needed fields} \label{sec:DOM_fields} The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}. The associated scale factors are defined using the analytical first derivative of the transformation \autoref{eq:scale_factors}. Necessary fields for configuration definition are: The total size of the computational domain is set by the parameters \jp{jpiglo}, \jp{jpjglo} and \jp{jpkglo} for the $i$, $j$ and $k$ directions, respectively. Note, that the variables \texttt{jpi} and \texttt{jpj} refer to the size of each processor subdomain when the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined, see \autoref{sec:LBC_mpp}). The name of the configuration is set through parameter \np{cn_cfg}{cn\_cfg}, and the nominal resolution through parameter \np{nn_cfg}{nn\_cfg} (unless in the input file both of variables \texttt{ORCA} and \texttt{ORCA\_index} are present, in which case \np{cn_cfg}{cn\_cfg} and \np{nn_cfg}{nn\_cfg} are set from these values accordingly). The global lateral boundary condition type is selected from 8 options using parameter \jp{jperio}. See \autoref{sec:LBC_jperio} for details on the available options and the corresponding values for \jp{jperio}. %% ================================================================================================= \subsection[Horizontal grid mesh (\textit{domhgr.F90}]{Horizontal grid mesh (\protect\mdl{domhgr})} \label{subsec:DOM_hgr} %% ================================================================================================= \subsubsection{Required fields} \label{sec:DOM_hgr_fields} The explicit specification of a range of mesh-related fields are required for the definition of a configuration. These include: \begin{clines} int    jpiglo, jpjglo, jpkglo     /* global domain sizes                                    */ int    jperio                     /* lateral global domain b.c.                             */ double glamt, glamu, glamv, glamf /* geographic longitude (t,u,v and f points respectively) */ double gphit, gphiu, gphiv, gphif /* geographic latitude                                    */ double e1t, e1u, e1v, e1f         /* horizontal scale factors                               */ double e2t, e2u, e2v, e2f         /* horizontal scale factors                               */ \end{clines} The values of the geographic longitude and latitude arrays at indices $i,j$ correspond to the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$, evaluated at the values as specified in \autoref{tab:DOM_cell} for the respective grid-point position. The calculation of the values of the horizontal scale factor arrays in general additionally involves partial derivatives of $\lambda$ and $\varphi$ with respect to $i$ and $j$, evaluated for the same arguments as $\lambda$ and $\varphi$. %% ================================================================================================= \subsubsection{Optional fields} \begin{clines} /* Optional:                                                 */ int    ORCA, ORCA_index /* configuration name, configuration resolution              */ double e1e2u, e1e2v     /* U and V surfaces (if grid size reduction in some straits) */ double ff_f, ff_t       /* Coriolis parameter (if not on the sphere)                 */ \end{clines} \NEMO\ can support the local reduction of key strait widths by altering individual values of e2u or e1v at the appropriate locations. This is particularly useful for locations such as Gibraltar or Indonesian Throughflow pinch-points (see \autoref{sec:MISC_strait} for illustrated examples). The key is to reduce the faces of $T$-cell (\ie\ change the value of the horizontal scale factors at $u$- or $v$-point) but not the volume of the cells. Doing otherwise can lead to numerical instability issues. In normal operation the surface areas are computed from $e1u * e2u$ and $e1v * e2v$ but in cases where a gridsize reduction is required, the unaltered surface areas at $u$ and $v$ grid points (\texttt{e1e2u} and \texttt{e1e2v}, respectively) must be read or pre-computed in \mdl{usrdef\_hgr}. If these arrays are present in the \np{cn_domcfg}{cn\_domcfg} file they are read and the internal computation is suppressed. Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{e1e2u} and \texttt{e1e2v} should set the surface-area computation flag: \texttt{ie1e2u\_v} to a non-zero value to suppress their re-computation. \smallskip Similar logic applies to the other optional fields: \texttt{ff\_f} and \texttt{ff\_t} which can be used to provide the Coriolis parameter at F- and T-points respectively if the mesh is not on a sphere. If present these fields will be read and used and the normal calculation ($2 * \Omega * \sin(\varphi)$) suppressed. Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{ff\_f} and \texttt{ff\_t} should set the Coriolis computation flag: \texttt{iff} to a non-zero value to suppress their re-computation. Note that longitudes, latitudes, and scale factors at $w$ points are exactly equal to those of $t$ points, thus no specific arrays are defined at $w$ points. %% ================================================================================================= \subsection[Vertical grid (\textit{domzgr.F90})]{Vertical grid (\protect\mdl{domzgr})} \label{subsec:DOM_zgr} \begin{listing} \nlst{namdom} \caption{\forcode{&namdom}} \label{lst:namdom} \end{listing} In the vertical, the model mesh is determined by four things: \begin{enumerate} \item the bathymetry given in meters; \item the number of levels of the model (\jp{jpk}); \item the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and \item the masking system, \ie\ the number of wet model levels at each $(i,j)$ location of the horizontal grid. \end{enumerate} \begin{figure} \centering \includegraphics[width=0.5\textwidth]{DOM_z_zps_s_sps} \caption[Ocean bottom regarding coordinate systems ($z$, $s$ and hybrid $s-z$)]{ The ocean bottom as seen by the model: \begin{enumerate*}[label=(\textit{\alph*})] \item $z$-coordinate with full step, \item $z$-coordinate with partial step, \item $s$-coordinate: terrain following representation, \item hybrid $s-z$ coordinate, \item hybrid $s-z$ coordinate with partial step, and \item same as (e) but in the non-linear free surface (\protect\np[=.false.]{ln_linssh}{ln\_linssh}). \end{enumerate*} Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e).} \label{fig:DOM_z_zps_s_sps} \end{figure} The choice of a vertical coordinate is made when setting up the configuration; it is not intended to be an option which can be changed in the middle of an experiment. The one exception to this statement being the choice of linear or non-linear free surface. In v4.0 the linear free surface option is implemented as a special case of the non-linear free surface. This is computationally wasteful since it uses the structures for time-varying 3D metrics for fields that (in the linear free surface case) are fixed. However, the linear free-surface is rarely used and implementing it this way means a single configuration file can support both options. By default a non-linear free surface is used (\np{ln_linssh}{ln\_linssh} set to \forcode{=.false.} in \nam{dom}{dom}): the coordinate follow the time-variation of the free surface so that the transformation is time dependent: $z(i,j,k,t)$ (\eg\ \autoref{fig:DOM_z_zps_s_sps}f). When a linear free surface is assumed (\np{ln_linssh}{ln\_linssh} set to \forcode{=.true.} in \nam{dom}{dom}), the vertical coordinates are fixed in time, but the seawater can move up and down across the $z_0$ surface (in other words, the top of the ocean in not a rigid lid). Note that settings: \np{ln_zco}{ln\_zco}, \np{ln_zps}{ln\_zps}, \np{ln_sco}{ln\_sco} and \np{ln_isfcav}{ln\_isfcav} mentioned in the following sections appear to be namelist options but they are no longer truly namelist options for \NEMO. Their value is written to and read from the domain configuration file and they should be treated as fixed parameters for a particular configuration. They are namelist options for the \texttt{DOMAINcfg} tool that can be used to build the configuration file and serve both to provide a record of the choices made whilst building the configuration and to trigger appropriate code blocks within \NEMO. These values should not be altered in the \np{cn_domcfg}{cn\_domcfg} file. \medskip The decision on these choices must be made when the \np{cn_domcfg}{cn\_domcfg} file is constructed. Three main choices are offered (\autoref{fig:DOM_z_zps_s_sps}a-c): \begin{itemize} \item Geographic position: longitude with \texttt{glamt}, \texttt{glamu}, \texttt{glamv}, \texttt{glamf} and latitude  with \texttt{gphit}, \texttt{gphiu}, \texttt{gphiv}, \texttt{gphif} (all respectively at T, U, V and F point) \item Coriolis parameter (if domain not on the sphere): \texttt{ff\_f} and \texttt{ff\_t} (at T and F point) \item Scale factors: \texttt{e1t}, \texttt{e1u}, \texttt{e1v} and \texttt{e1f} (on i direction), \texttt{e2t}, \texttt{e2u}, \texttt{e2v} and \texttt{e2f} (on j direction) and \texttt{ie1e2u\_v}, \texttt{e1e2u}, \texttt{e1e2v}. \\ \texttt{e1e2u}, \texttt{e1e2v} are u and v surfaces (if gridsize reduction in some straits), \texttt{ie1e2u\_v} is to flag set u and v surfaces are neither read nor computed. \item $z$-coordinate with full step bathymetry (\np[=.true.]{ln_zco}{ln\_zco}), \item $z$-coordinate with partial step ($zps$) bathymetry (\np[=.true.]{ln_zps}{ln\_zps}), \item Generalized, $s$-coordinate (\np[=.true.]{ln_sco}{ln\_sco}). \end{itemize} These fields can be read in an domain input file which name is setted in \np{cn\_domcfg} parameter specified in \ngn{namcfg}. \nlst{namcfg} Or they can be defined in an analytical way in \path{MY_SRC} directory of the configuration. For Reference Configurations of NEMO input domain files are supplied by NEMO System Team. For analytical definition of input fields two routines are supplied: \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. They are an example of GYRE configuration parameters, and they are available in \path{src/OCE/USR} directory, they provide the horizontal and vertical mesh. % ------------------------------------------------------------------------------------------------------------- %        Needed fields % ------------------------------------------------------------------------------------------------------------- %\subsection{List of needed fields to build DOMAIN} %\label{subsec:DOM_fields_list} % ================================================================ % Domain: Horizontal Grid (mesh) % ================================================================ \section{Horizontal grid mesh (\protect\mdl{domhgr})} \label{sec:DOM_hgr} % ------------------------------------------------------------------------------------------------------------- %        Coordinates and scale factors % ------------------------------------------------------------------------------------------------------------- \subsection{Coordinates and scale factors} \label{subsec:DOM_hgr_coord_e} The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}. The associated scale factors are defined using the analytical first derivative of the transformation \autoref{eq:scale_factors}. These definitions are done in two modules, \mdl{domhgr} and \mdl{domzgr}, which provide the horizontal and vertical meshes, respectively. This section deals with the horizontal mesh parameters. In a horizontal plane, the location of all the model grid points is defined from the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$. The horizontal scale factors are calculated using \autoref{eq:scale_factors}. For example, when the longitude and latitude are function of a single value ($i$ and $j$, respectively) (geographical configuration of the mesh), the horizontal mesh definition reduces to define the wanted $\lambda(i)$, $\varphi(j)$, and their derivatives $\lambda'(i) \ \varphi'(j)$ in the \mdl{domhgr} module. The model computes the grid-point positions and scale factors in the horizontal plane as follows: Additionally, hybrid combinations of the three main coordinates are available: $s-z$ or $s-zps$ coordinate (\autoref{fig:DOM_z_zps_s_sps}d and \autoref{fig:DOM_z_zps_s_sps}e). A further choice related to vertical coordinate concerns the presence (or not) of ocean cavities beneath ice shelves within the model domain. A setting of \np{ln_isfcav}{ln\_isfcav} as \forcode{.true.} indicates that the domain contains ocean cavities, otherwise the top, wet layer of the ocean will always be at the ocean surface. This option is currently only available for $z$- or $zps$-coordinates. In the latter case, partial steps are also applied at the ocean/ice shelf interface. Within the model, the arrays describing the grid point depths and vertical scale factors are three set of three dimensional arrays $(i,j,k)$ defined at \textit{before}, \textit{now} and \textit{after} time step. The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively. They are updated at each model time step. The initial fixed reference coordinate system is held in variable names with a $\_0$ suffix. When the linear free surface option is used (\np[=.true.]{ln_linssh}{ln\_linssh}), \textit{before}, \textit{now} and \textit{after} arrays are initially set to their reference counterpart and remain fixed. %% ================================================================================================= \subsubsection{Required fields} \label{sec:DOM_zgr_fields} The explicit specification of a range of fields related to the vertical grid are required for the definition of a configuration. These include: \begin{clines} int    ln_zco, ln_zps, ln_sco            /* flags for z-coord, z-coord with partial steps and s-coord    */ int    ln_isfcav                         /* flag  for ice shelf cavities                                 */ double e3t_1d, e3w_1d                    /* reference vertical scale factors at T and W points           */ double e3t_0, e3u_0, e3v_0, e3f_0, e3w_0 /* vertical scale factors 3D coordinate at T,U,V,F and W points */ double e3uw_0, e3vw_0                    /* vertical scale factors 3D coordinate at UW and VW points     */ int    bottom_level, top_level           /* last wet T-points, 1st wet T-points (for ice shelf cavities) */ /* For reference:                                               */ float  bathy_metry                       /* bathymetry used in setting top and bottom levels             */ \end{clines} This set of vertical metrics is sufficient to describe the initial depth and thickness of every gridcell in the model regardless of the choice of vertical coordinate. With constant z-levels, e3 metrics will be uniform across each horizontal level. In the partial step case each e3 at the \jp{bottom\_level} (and, possibly, \jp{top\_level} if ice cavities are present) may vary from its horizontal neighbours. And, in s-coordinates, variations can occur throughout the water column. With the non-linear free-surface, all the coordinates behave more like the s-coordinate in that variations occur throughout the water column with displacements related to the sea surface height. These variations are typically much smaller than those arising from bottom fitted coordinates. The values for vertical metrics supplied in the domain configuration file can be considered as those arising from a flat sea surface with zero elevation. The \jp{bottom\_level} and \jp{top\_level} 2D arrays define the \jp{bottom\_level} and top wet levels in each grid column. Without ice cavities, \jp{top\_level} is essentially a land mask (0 on land; 1 everywhere else). With ice cavities, \jp{top\_level} determines the first wet point below the overlying ice shelf. %% ================================================================================================= \subsubsection{Level bathymetry and mask} \label{subsec:DOM_msk} From \jp{top\_level} and \jp{bottom\_level} fields, the mask fields are defined as follows: \begin{align*} \lambda_t &\equiv \text{glamt} =      \lambda (i      ) &\varphi_t &\equiv \text{gphit} =      \varphi (j      ) \\ \lambda_u &\equiv \text{glamu} =      \lambda (i + 1/2) &\varphi_u &\equiv \text{gphiu} =      \varphi (j      ) \\ \lambda_v &\equiv \text{glamv} =      \lambda (i      ) &\varphi_v &\equiv \text{gphiv} =      \varphi (j + 1/2) \\ \lambda_f &\equiv \text{glamf} =      \lambda (i + 1/2) &\varphi_f &\equiv \text{gphif} =      \varphi (j + 1/2) \\ e_{1t}    &\equiv \text{e1t}   = r_a |\lambda'(i      ) \; \cos\varphi(j      ) | &e_{2t}    &\equiv \text{e2t}   = r_a |\varphi'(j      )                         | \\ e_{1u}    &\equiv \text{e1t}   = r_a |\lambda'(i + 1/2) \; \cos\varphi(j      ) | &e_{2u}    &\equiv \text{e2t}   = r_a |\varphi'(j      )                         | \\ e_{1v}    &\equiv \text{e1t}   = r_a |\lambda'(i      ) \; \cos\varphi(j + 1/2) | &e_{2v}    &\equiv \text{e2t}   = r_a |\varphi'(j + 1/2)                         | \\ e_{1f}    &\equiv \text{e1t}   = r_a |\lambda'(i + 1/2) \; \cos\varphi(j + 1/2) | &e_{2f}    &\equiv \text{e2t}   = r_a |\varphi'(j + 1/2)                         | tmask(i,j,k) &= \begin{cases} 0 &\text{if $k < top\_level(i,j)$} \\ 1 &\text{if $bottom\_level(i,j) \leq k \leq top\_level(i,j)$} \\ 0 &\text{if $k > bottom\_level(i,j)$} \end{cases} \\ umask(i,j,k) &= tmask(i,j,k) * tmask(i + 1,j,    k) \\ vmask(i,j,k) &= tmask(i,j,k) * tmask(i    ,j + 1,k) \\ fmask(i,j,k) &= tmask(i,j,k) * tmask(i + 1,j,    k) * tmask(i,j,k) * tmask(i + 1,j,    k) \\ wmask(i,j,k) &= tmask(i,j,k) * tmask(i    ,j,k - 1) \\ \text{with~} wmask(i,j,1) &= tmask(i,j,1) \end{align*} where the last letter of each computational name indicates the grid point considered and $r_a$ is the earth radius (defined in \mdl{phycst} along with all universal constants). Note that the horizontal position of and scale factors at $w$-points are exactly equal to those of $t$-points, thus no specific arrays are defined at $w$-points. Note that the definition of the scale factors (\ie as the analytical first derivative of the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) is specific to the \NEMO model \citep{Marti_al_JGR92}. As an example, $e_{1t}$ is defined locally at a $t$-point, whereas many other models on a C grid choose to define such a scale factor as the distance between the $U$-points on each side of the $t$-point. Relying on an analytical transformation has two advantages: firstly, there is no ambiguity in the scale factors appearing in the discrete equations, since they are first introduced in the continuous equations; secondly, analytical transformations encourage good practice by the definition of smoothly varying grids (rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{Treguier1996}. An example of the effect of such a choice is shown in \autoref{fig:zgr_e3}. %>>>>>>>>>>>>>>>>>>>>>>>>>>>> \begin{figure}[!t] \begin{center} \includegraphics[]{Fig_zgr_e3} \caption{ \protect\label{fig:zgr_e3} Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, and (b) analytically derived grid-point position and scale factors. For both grids here, the same $w$-point depth has been chosen but in (a) the $t$-points are set half way between $w$-points while in (b) they are defined from an analytical function: $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$. Note the resulting difference between the value of the grid-size $\Delta_k$ and those of the scale factor $e_k$. } \end{center} \end{figure} %>>>>>>>>>>>>>>>>>>>>>>>>>>>> % ------------------------------------------------------------------------------------------------------------- %        Choice of horizontal grid % ------------------------------------------------------------------------------------------------------------- \subsection{Choice of horizontal grid} \label{subsec:DOM_hgr_msh_choice} % ------------------------------------------------------------------------------------------------------------- %        Grid files % ------------------------------------------------------------------------------------------------------------- \subsection{Output grid files} \label{subsec:DOM_hgr_files} All the arrays relating to a particular ocean model configuration (grid-point position, scale factors, masks) can be saved in files if \np{nn\_msh} $\not = 0$ (namelist variable in \ngn{namdom}). This can be particularly useful for plots and off-line diagnostics. In some cases, the user may choose to make a local modification of a scale factor in the code. This is the case in global configurations when restricting the width of a specific strait (usually a one-grid-point strait that happens to be too wide due to insufficient model resolution). An example is Gibraltar Strait in the ORCA2 configuration. When such modifications are done, the output grid written when \np{nn\_msh} $\not = 0$ is no more equal to the input grid. % ================================================================ % Domain: Vertical Grid (domzgr) % ================================================================ \section{Vertical grid (\protect\mdl{domzgr})} \label{sec:DOM_zgr} %-----------------------------------------nam_zgr & namdom------------------------------------------- % %\nlst{namzgr} \nlst{namdom} %------------------------------------------------------------------------------------------------------------- Variables are defined through the \ngn{namzgr} and \ngn{namdom} namelists. In the vertical, the model mesh is determined by four things: (1) the bathymetry given in meters; (2) the number of levels of the model (\jp{jpk}); (3) the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and (4) the masking system, \ie the number of wet model levels at each $(i,j)$ column of points. %>>>>>>>>>>>>>>>>>>>>>>>>>>>> \begin{figure}[!tb] \begin{center} \includegraphics[]{Fig_z_zps_s_sps} \caption{ \protect\label{fig:z_zps_s_sps} The ocean bottom as seen by the model: (a) $z$-coordinate with full step, (b) $z$-coordinate with partial step, (c) $s$-coordinate: terrain following representation, (d) hybrid $s-z$ coordinate, (e) hybrid $s-z$ coordinate with partial step, and (f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh}~\forcode{= .false.}). Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e). } \end{center} \end{figure} %>>>>>>>>>>>>>>>>>>>>>>>>>>>> The choice of a vertical coordinate, even if it is made through \ngn{namzgr} namelist parameters, must be done once of all at the beginning of an experiment. It is not intended as an option which can be enabled or disabled in the middle of an experiment. Three main choices are offered (\autoref{fig:z_zps_s_sps}): $z$-coordinate with full step bathymetry (\np{ln\_zco}~\forcode{= .true.}), $z$-coordinate with partial step bathymetry (\np{ln\_zps}~\forcode{= .true.}), or generalized, $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}). Hybridation of the three main coordinates are available: $s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps} and \autoref{fig:z_zps_s_sps}). By default a non-linear free surface is used: the coordinate follow the time-variation of the free surface so that the transformation is time dependent: $z(i,j,k,t)$ (\autoref{fig:z_zps_s_sps}). When a linear free surface is assumed (\np{ln\_linssh}~\forcode{= .true.}), the vertical coordinate are fixed in time, but the seawater can move up and down across the $z_0$ surface (in other words, the top of the ocean in not a rigid-lid). The last choice in terms of vertical coordinate concerns the presence (or not) in the model domain of ocean cavities beneath ice shelves. Setting \np{ln\_isfcav} to true allows to manage ocean cavities, otherwise they are filled in. This option is currently only available in $z$- or $zps$-coordinate, and partial step are also applied at the ocean/ice shelf interface. Contrary to the horizontal grid, the vertical grid is computed in the code and no provision is made for reading it from a file. The only input file is the bathymetry (in meters) (\ifile{bathy\_meter}) \footnote{ N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the \ifile{bathy\_meter} file, so that the computation of the number of wet ocean point in each water column is by-passed}. If \np{ln\_isfcav}~\forcode{= .true.}, an extra file input file (\ifile{isf\_draft\_meter}) describing the ice shelf draft (in meters) is needed. After reading the bathymetry, the algorithm for vertical grid definition differs between the different options: \begin{description} \item[\textit{zco}] set a reference coordinate transformation $z_0(k)$, and set $z(i,j,k,t) = z_0(k)$. \item[\textit{zps}] set a reference coordinate transformation $z_0(k)$, and calculate the thickness of the deepest level at each $(i,j)$ point using the bathymetry, to obtain the final three-dimensional depth and scale factor arrays. \item[\textit{sco}] smooth the bathymetry to fulfill the hydrostatic consistency criteria and set the three-dimensional transformation. \item[\textit{s-z} and \textit{s-zps}] smooth the bathymetry to fulfill the hydrostatic consistency criteria and set the three-dimensional transformation $z(i,j,k)$, and possibly introduce masking of extra land points to better fit the original bathymetry file. \end{description} %%% \gmcomment{   add the description of the smoothing:  envelop topography...} %%% Unless a linear free surface is used (\np{ln\_linssh}~\forcode{= .false.}), the arrays describing the grid point depths and vertical scale factors are three set of three dimensional arrays $(i,j,k)$ defined at \textit{before}, \textit{now} and \textit{after} time step. The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively. They are updated at each model time step using a fixed reference coordinate system which computer names have a $\_0$ suffix. When the linear free surface option is used (\np{ln\_linssh}~\forcode{= .true.}), \textit{before}, \textit{now} and \textit{after} arrays are simply set one for all to their reference counterpart. % ------------------------------------------------------------------------------------------------------------- %        Meter Bathymetry % ------------------------------------------------------------------------------------------------------------- \subsection{Meter bathymetry} \label{subsec:DOM_bathy} Three options are possible for defining the bathymetry, according to the namelist variable \np{nn\_bathy} (found in \ngn{namdom} namelist): \begin{description} \item[\np{nn\_bathy}~\forcode{= 0}]: a flat-bottom domain is defined. The total depth $z_w (jpk)$ is given by the coordinate transformation. The domain can either be a closed basin or a periodic channel depending on the parameter \np{jperio}. \item[\np{nn\_bathy}~\forcode{= -1}]: a domain with a bump of topography one third of the domain width at the central latitude. This is meant for the "EEL-R5" configuration, a periodic or open boundary channel with a seamount. \item[\np{nn\_bathy}~\forcode{= 1}]: read a bathymetry and ice shelf draft (if needed). The \ifile{bathy\_meter} file (Netcdf format) provides the ocean depth (positive, in meters) at each grid point of the model grid. The bathymetry is usually built by interpolating a standard bathymetry product (\eg ETOPO2) onto the horizontal ocean mesh. Defining the bathymetry also defines the coastline: where the bathymetry is zero, no model levels are defined (all levels are masked). The \ifile{isfdraft\_meter} file (Netcdf format) provides the ice shelf draft (positive, in meters) at each grid point of the model grid. This file is only needed if \np{ln\_isfcav}~\forcode{= .true.}. Defining the ice shelf draft will also define the ice shelf edge and the grounding line position. \end{description} When a global ocean is coupled to an atmospheric model it is better to represent all large water bodies (\eg great lakes, Caspian sea...) even if the model resolution does not allow their communication with the rest of the ocean. Note that, without ice shelves cavities, masks at $t-$ and $w-$points are identical with the numerical indexing used (\autoref{subsec:DOM_Num_Index}). Nevertheless, $wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface) exactly in the same way as for the bottom boundary. %% The specification of closed lateral boundaries requires that at least %% the first and last rows and columns of the \textit{mbathy} array are set to zero. %% In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to %% the second one and its first column equal to the last but one (and so too the mask arrays) %% (see \autoref{fig:LBC_jperio}). %        Closed seas %% ================================================================================================= \subsection{Closed seas} \label{subsec:DOM_closea} When a global ocean is coupled to an atmospheric model it is better to represent all large water bodies (\eg\ Great Lakes, Caspian sea, \dots) even if the model resolution does not allow their communication with the rest of the ocean. This is unnecessary when the ocean is forced by fixed atmospheric conditions, so these seas can be removed from the ocean domain. The user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}), but the code has to be adapted to the user's configuration. % ------------------------------------------------------------------------------------------------------------- %        z-coordinate  and reference coordinate transformation % ------------------------------------------------------------------------------------------------------------- \subsection[$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and ref. coordinate] {$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and reference coordinate} \label{subsec:DOM_zco} %>>>>>>>>>>>>>>>>>>>>>>>>>>>> \begin{figure}[!tb] \begin{center} \includegraphics[]{Fig_zgr} \caption{ \protect\label{fig:zgr} Default vertical mesh for ORCA2: 30 ocean levels (L30). Vertical level functions for (a) T-point depth and (b) the associated scale factor as computed from \autoref{eq:DOM_zgr_ana_1} using \autoref{eq:DOM_zgr_coef} in $z$-coordinate. } \end{center} \end{figure} %>>>>>>>>>>>>>>>>>>>>>>>>>>>> The reference coordinate transformation $z_0(k)$ defines the arrays $gdept_0$ and $gdepw_0$ for $t$- and $w$-points, respectively. As indicated on \autoref{fig:index_vert} \jp{jpk} is the number of $w$-levels. $gdepw_0(1)$ is the ocean surface. There are at most \jp{jpk}-1 $t$-points inside the ocean, the additional $t$-point at $jk = jpk$ is below the sea floor and is not used. The vertical location of $w$- and $t$-levels is defined from the analytic expression of the depth $z_0(k)$ whose analytical derivative with respect to $k$ provides the vertical scale factors. The user must provide the analytical expression of both $z_0$ and its first derivative with respect to $k$. This is done in routine \mdl{domzgr} through statement functions, using parameters provided in the \ngn{namcfg} namelist. It is possible to define a simple regular vertical grid by giving zero stretching (\np{ppacr}~\forcode{= 0}). In that case, the parameters \jp{jpk} (number of $w$-levels) and \np{pphmax} (total ocean depth in meters) fully define the grid. For climate-related studies it is often desirable to concentrate the vertical resolution near the ocean surface. The following function is proposed as a standard for a $z$-coordinate (with either full or partial steps): \begin{gather} \label{eq:DOM_zgr_ana_1} z_0  (k) = h_{sur} - h_0 \; k - \; h_1 \; \log  \big[ \cosh ((k - h_{th}) / h_{cr}) \big] \\ e_3^0(k) = \lt|    - h_0      -    h_1 \; \tanh \big[        (k - h_{th}) / h_{cr}  \big] \rt| \end{gather} where $k = 1$ to \jp{jpk} for $w$-levels and $k = 1$ to $k = 1$ for $T-$levels. Such an expression allows us to define a nearly uniform vertical location of levels at the ocean top and bottom with a smooth hyperbolic tangent transition in between (\autoref{fig:zgr}). If the ice shelf cavities are opened (\np{ln\_isfcav}~\forcode{= .true.}), the definition of $z_0$ is the same. However, definition of $e_3^0$ at $t$- and $w$-points is respectively changed to: \label{eq:DOM_zgr_ana_2} \begin{split} e_3^T(k) &= z_W (k + 1) - z_W (k    ) \\ e_3^W(k) &= z_T (k    ) - z_T (k - 1) \end{split} This formulation decrease the self-generated circulation into the ice shelf cavity (which can, in extreme case, leads to blow up).\\ The most used vertical grid for ORCA2 has $10~m$ ($500~m$) resolution in the surface (bottom) layers and a depth which varies from 0 at the sea surface to a minimum of $-5000~m$. This leads to the following conditions: \label{eq:DOM_zgr_coef} \begin{array}{ll} e_3 (1   + 1/2) =  10. & z(1  ) =     0. \\ e_3 (jpk - 1/2) = 500. & z(jpk) = -5000. \end{array} With the choice of the stretching $h_{cr} = 3$ and the number of levels \jp{jpk}~$= 31$, the four coefficients $h_{sur}$, $h_0$, $h_1$, and $h_{th}$ in \autoref{eq:DOM_zgr_ana_2} have been determined such that \autoref{eq:DOM_zgr_coef} is satisfied, through an optimisation procedure using a bisection method. For the first standard ORCA2 vertical grid this led to the following values: $h_{sur} = 4762.96$, $h_0 = 255.58, h_1 = 245.5813$, and $h_{th} = 21.43336$. The resulting depths and scale factors as a function of the model levels are shown in \autoref{fig:zgr} and given in \autoref{tab:orca_zgr}. Those values correspond to the parameters \np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist. Rather than entering parameters $h_{sur}$, $h_0$, and $h_1$ directly, it is possible to recalculate them. In that case the user sets \np{ppsur}~$=$~\np{ppa0}~$=$~\np{ppa1}~$= 999999$., in \ngn{namcfg} namelist, and specifies instead the four following parameters: \begin{itemize} \item \np{ppacr}~$= h_{cr}$: stretching factor (nondimensional). The larger \np{ppacr}, the smaller the stretching. Values from $3$ to $10$ are usual. \item \np{ppkth}~$= h_{th}$: is approximately the model level at which maximum stretching occurs (nondimensional, usually of order 1/2 or 2/3 of \jp{jpk}) \item \np{ppdzmin}: minimum thickness for the top layer (in meters). \item \np{pphmax}: total depth of the ocean (meters). \end{itemize} As an example, for the $45$ layers used in the DRAKKAR configuration those parameters are: \jp{jpk}~$= 46$, \np{ppacr}~$= 9$, \np{ppkth}~$= 23.563$, \np{ppdzmin}~$= 6~m$, \np{pphmax}~$= 5750~m$. %>>>>>>>>>>>>>>>>>>>>>>>>>>>> \begin{table} \begin{center} \begin{tabular}{c||r|r|r|r} \hline \textbf{LEVEL} & \textbf{gdept\_1d} & \textbf{gdepw\_1d} & \textbf{e3t\_1d } & \textbf{e3w\_1d} \\ \hline 1              & \textbf{     5.00} &               0.00 & \textbf{   10.00} &            10.00 \\ \hline 2              & \textbf{    15.00} &              10.00 & \textbf{   10.00} &            10.00 \\ \hline 3              & \textbf{    25.00} &              20.00 & \textbf{   10.00} &            10.00 \\ \hline 4              & \textbf{    35.01} &              30.00 & \textbf{   10.01} &            10.00 \\ \hline 5              & \textbf{    45.01} &              40.01 & \textbf{   10.01} &            10.01 \\ \hline 6              & \textbf{    55.03} &              50.02 & \textbf{   10.02} &            10.02 \\ \hline 7              & \textbf{    65.06} &              60.04 & \textbf{   10.04} &            10.03 \\ \hline 8              & \textbf{    75.13} &              70.09 & \textbf{   10.09} &            10.06 \\ \hline 9              & \textbf{    85.25} &              80.18 & \textbf{   10.17} &            10.12 \\ \hline 10             & \textbf{    95.49} &              90.35 & \textbf{   10.33} &            10.24 \\ \hline 11             & \textbf{   105.97} &             100.69 & \textbf{   10.65} &            10.47 \\ \hline 12             & \textbf{   116.90} &             111.36 & \textbf{   11.27} &            10.91 \\ \hline 13             & \textbf{   128.70} &             122.65 & \textbf{   12.47} &            11.77 \\ \hline 14             & \textbf{   142.20} &             135.16 & \textbf{   14.78} &            13.43 \\ \hline 15             & \textbf{   158.96} &             150.03 & \textbf{   19.23} &            16.65 \\ \hline 16             & \textbf{   181.96} &             169.42 & \textbf{   27.66} &            22.78 \\ \hline 17             & \textbf{   216.65} &             197.37 & \textbf{   43.26} &            34.30 \\ \hline 18             & \textbf{   272.48} &             241.13 & \textbf{   70.88} &            55.21 \\ \hline 19             & \textbf{   364.30} &             312.74 & \textbf{  116.11} &            90.99 \\ \hline 20             & \textbf{   511.53} &             429.72 & \textbf{  181.55} &           146.43 \\ \hline 21             & \textbf{   732.20} &             611.89 & \textbf{  261.03} &           220.35 \\ \hline 22             & \textbf{  1033.22} &             872.87 & \textbf{  339.39} &           301.42 \\ \hline 23             & \textbf{  1405.70} &            1211.59 & \textbf{  402.26} &           373.31 \\ \hline 24             & \textbf{  1830.89} &            1612.98 & \textbf{  444.87} &           426.00 \\ \hline 25             & \textbf{  2289.77} &            2057.13 & \textbf{  470.55} &           459.47 \\ \hline 26             & \textbf{  2768.24} &            2527.22 & \textbf{  484.95} &           478.83 \\ \hline 27             & \textbf{  3257.48} &            3011.90 & \textbf{  492.70} &           489.44 \\ \hline 28             & \textbf{  3752.44} &            3504.46 & \textbf{  496.78} &           495.07 \\ \hline 29             & \textbf{  4250.40} &            4001.16 & \textbf{  498.90} &           498.02 \\ \hline 30             & \textbf{  4749.91} &            4500.02 & \textbf{  500.00} &           499.54 \\ \hline 31             & \textbf{  5250.23} &            5000.00 & \textbf{  500.56} &           500.33 \\ \hline \end{tabular} \end{center} \caption{ \protect\label{tab:orca_zgr} Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration as computed from \autoref{eq:DOM_zgr_ana_2} using the coefficients given in \autoref{eq:DOM_zgr_coef} } \end{table} %>>>>>>>>>>>>>>>>>>>>>>>>>>>> % ------------------------------------------------------------------------------------------------------------- %        z-coordinate with partial step % ------------------------------------------------------------------------------------------------------------- \subsection{$Z$-coordinate with partial step (\protect\np{ln\_zps}~\forcode{= .true.})} \label{subsec:DOM_zps} %--------------------------------------------namdom------------------------------------------------------- \nlst{namdom} %-------------------------------------------------------------------------------------------------------------- In $z$-coordinate partial step, the depths of the model levels are defined by the reference analytical function $z_0(k)$ as described in the previous section, \textit{except} in the bottom layer. The thickness of the bottom layer is allowed to vary as a function of geographical location $(\lambda,\varphi)$ to allow a better representation of the bathymetry, especially in the case of small slopes (where the bathymetry varies by less than one level thickness from one grid point to the next). The reference layer thicknesses $e_{3t}^0$ have been defined in the absence of bathymetry. With partial steps, layers from 1 to \jp{jpk}-2 can have a thickness smaller than $e_{3t}(jk)$. The model deepest layer (\jp{jpk}-1) is allowed to have either a smaller or larger thickness than $e_{3t}(jpk)$: the maximum thickness allowed is $2*e_{3t}(jpk - 1)$. This has to be kept in mind when specifying values in \ngn{namdom} namelist, as the maximum depth \np{pphmax} in partial steps: for example, with \np{pphmax}~$= 5750~m$ for the DRAKKAR 45 layer grid, the maximum ocean depth allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk - 1)$ being $250~m$). Two variables in the namdom namelist are used to define the partial step vertical grid. The mimimum water thickness (in meters) allowed for a cell partially filled with bathymetry at level jk is the minimum of \np{rn\_e3zps\_min} (thickness in meters, usually $20~m$) or $e_{3t}(jk)*$\np{rn\_e3zps\_rat} (a fraction, usually 10\%, of the default thickness $e_{3t}(jk)$). \gmcomment{ \colorbox{yellow}{Add a figure here of pstep especially at last ocean level }  } % ------------------------------------------------------------------------------------------------------------- %        s-coordinate % ------------------------------------------------------------------------------------------------------------- \subsection{$S$-coordinate (\protect\np{ln\_sco}~\forcode{= .true.})} \label{subsec:DOM_sco} %------------------------------------------nam_zgr_sco--------------------------------------------------- % %\nlst{namzgr_sco} %-------------------------------------------------------------------------------------------------------------- Options are defined in \ngn{namzgr\_sco}. In $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}), the depth and thickness of the model levels are defined from the product of a depth field and either a stretching function or its derivative, respectively: \begin{align*} % \label{eq:DOM_sco_ana} z(k)   &= h(i,j) \; z_0 (k) \\ e_3(k) &= h(i,j) \; z_0'(k) \end{align*} where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point location in the horizontal and $z_0(k)$ is a function which varies from $0$ at the sea surface to $1$ at the ocean bottom. The depth field $h$ is not necessary the ocean depth, since a mixed step-like and bottom-following representation of the topography can be used (\autoref{fig:z_zps_s_sps}) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}). The namelist parameter \np{rn\_rmax} determines the slope at which the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate. The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} as the minimum and maximum depths at which the terrain-following vertical coordinate is calculated. Options for stretching the coordinate are provided as examples, but care must be taken to ensure that the vertical stretch used is appropriate for the application. The original default NEMO s-coordinate stretching is available if neither of the other options are specified as true (\np{ln\_s\_SH94}~\forcode{= .false.} and \np{ln\_s\_SF12}~\forcode{= .false.}). This uses a depth independent $\tanh$ function for the stretching \citep{Madec_al_JPO96}: $z = s_{min} + C (s) (H - s_{min}) % \label{eq:SH94_1}$ where $s_{min}$ is the depth at which the $s$-coordinate stretching starts and allows a $z$-coordinate to placed on top of the stretched coordinate, and $z$ is the depth (negative down from the asea surface). \begin{gather*} s = - \frac{k}{n - 1} \quad \text{and} \quad 0 \leq k \leq n - 1 % \label{eq:DOM_s} \\ % \label{eq:DOM_sco_function} C(s) = \frac{[\tanh(\theta \, (s + b)) - \tanh(\theta \, b)]}{2 \; \sinh(\theta)} \end{gather*} A stretching function, modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_s\_SH94}~\forcode{= .true.}), is also available and is more commonly used for shelf seas modelling: $C(s) = (1 - b) \frac{\sinh(\theta s)}{\sinh(\theta)} + b \frac{\tanh \lt[ \theta \lt(s + \frac{1}{2} \rt) \rt] - \tanh \lt( \frac{\theta}{2} \rt)} { 2 \tanh \lt( \frac{\theta}{2} \rt)} % \label{eq:SH94_2}$ %>>>>>>>>>>>>>>>>>>>>>>>>>>>> \begin{figure}[!ht] \begin{center} \includegraphics[]{Fig_sco_function} \caption{ \protect\label{fig:sco_function} Examples of the stretching function applied to a seamount; from left to right: surface, surface and bottom, and bottom intensified resolutions } \end{center} \end{figure} %>>>>>>>>>>>>>>>>>>>>>>>>>>>> where $H_c$ is the critical depth (\np{rn\_hc}) at which the coordinate transitions from pure $\sigma$ to the stretched coordinate, and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and bottom control parameters such that $0 \leqslant \theta \leqslant 20$, and $0 \leqslant b \leqslant 1$. $b$ has been designed to allow surface and/or bottom increase of the vertical resolution (\autoref{fig:sco_function}). Another example has been provided at version 3.5 (\np{ln\_s\_SF12}) that allows a fixed surface resolution in an analytical terrain-following stretching \citet{Siddorn_Furner_OM12}. In this case the a stretching function $\gamma$ is defined such that: z = - \gamma h \quad \text{with} \quad 0 \leq \gamma \leq 1 % \label{eq:z} The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate: \begin{gather*} % \label{eq:DOM_gamma_deriv} \gamma =   A \lt( \sigma   - \frac{1}{2} (\sigma^2     + f (\sigma)) \rt) + B \lt( \sigma^3 - f           (\sigma) \rt) + f (\sigma)       \\ \intertext{Where:} % \label{eq:DOM_gamma} f(\sigma) = (\alpha + 2) \sigma^{\alpha + 1} - (\alpha + 1) \sigma^{\alpha + 2} \quad \text{and} \quad \sigma = \frac{k}{n - 1} \end{gather*} This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of the user prescribed stretching parameter $\alpha$ (\np{rn\_alpha}) that stretches towards the surface ($\alpha > 1.0$) or the bottom ($\alpha < 1.0$) and user prescribed surface (\np{rn\_zs}) and bottom depths. The bottom cell depth in this example is given as a function of water depth: $% \label{eq:DOM_zb} Z_b = h a + b$ where the namelist parameters \np{rn\_zb\_a} and \np{rn\_zb\_b} are $a$ and $b$ respectively. %>>>>>>>>>>>>>>>>>>>>>>>>>>>> \begin{figure}[!ht] \includegraphics[]{Fig_DOM_compare_coordinates_surface} \caption{ A comparison of the \citet{Song_Haidvogel_JCP94} $S$-coordinate (solid lines), a 50 level $Z$-coordinate (contoured surfaces) and the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in the surface $100~m$ for a idealised bathymetry that goes from $50~m$ to $5500~m$ depth. For clarity every third coordinate surface is shown. } \label{fig:fig_compare_coordinates_surface} \end{figure} % >>>>>>>>>>>>>>>>>>>>>>>>>>>> This gives a smooth analytical stretching in computational space that is constrained to given specified surface and bottom grid cell thicknesses in real space. This is not to be confused with the hybrid schemes that superimpose geopotential coordinates on terrain following coordinates thus creating a non-analytical vertical coordinate that therefore may suffer from large gradients in the vertical resolutions. This stretching is less straightforward to implement than the \citet{Song_Haidvogel_JCP94} stretching, but has the advantage of resolving diurnal processes in deep water and has generally flatter slopes. As with the \citet{Song_Haidvogel_JCP94} stretching the stretch is only applied at depths greater than the critical depth $h_c$. In this example two options are available in depths shallower than $h_c$, with pure sigma being applied if the \np{ln\_sigcrit} is true and pure z-coordinates if it is false (the z-coordinate being equal to the depths of the stretched coordinate at $h_c$). Minimising the horizontal slope of the vertical coordinate is important in terrain-following systems as large slopes lead to hydrostatic consistency. A hydrostatic consistency parameter diagnostic following \citet{Haney1991} has been implemented, and is output as part of the model mesh file at the start of the run. % ------------------------------------------------------------------------------------------------------------- %        z*- or s*-coordinate % ------------------------------------------------------------------------------------------------------------- \subsection{\zstar- or \sstar-coordinate (\protect\np{ln\_linssh}~\forcode{= .false.})} \label{subsec:DOM_zgr_star} This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site. %gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances % ------------------------------------------------------------------------------------------------------------- %        level bathymetry and mask % ------------------------------------------------------------------------------------------------------------- \subsection{Level bathymetry and mask} \label{subsec:DOM_msk} Whatever the vertical coordinate used, the model offers the possibility of representing the bottom topography with steps that follow the face of the model cells (step like topography) \citep{Madec_al_JPO96}. The distribution of the steps in the horizontal is defined in a 2D integer array, mbathy, which gives the number of ocean levels (\ie those that are not masked) at each $t$-point. mbathy is computed from the meter bathymetry using the definiton of gdept as the number of $t$-points which gdept $\leq$ bathy. Modifications of the model bathymetry are performed in the \textit{bat\_ctl} routine (see \mdl{domzgr} module) after mbathy is computed. Isolated grid points that do not communicate with another ocean point at the same level are eliminated. As for the representation of bathymetry, a 2D integer array, misfdep, is created. misfdep defines the level of the first wet $t$-point. All the cells between $k = 1$ and $misfdep(i,j) - 1$ are masked. By default, $misfdep(:,:) = 1$ and no cells are masked. In case of ice shelf cavities, modifications of the model bathymetry and ice shelf draft into the cavities are performed in the \textit{zgr\_isf} routine. The compatibility between ice shelf draft and bathymetry is checked. All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded (\ie masked). If only one cell on the water column is opened at $t$-, $u$- or $v$-points, the bathymetry or the ice shelf draft is dug to fit this constrain. If the incompatibility is too strong (need to dig more than 1 cell), the cell is masked. From the \textit{mbathy} and \textit{misfdep} array, the mask fields are defined as follows: \begin{alignat*}{2} tmask(i,j,k) &= &  & \begin{cases} 0 &\text{if $k < misfdep(i,j)$} \\ 1 &\text{if $misfdep(i,j) \leq k \leq mbathy(i,j)$} \\ 0 &\text{if $k > mbathy(i,j)$} \end{cases} \\ umask(i,j,k) &= &  &tmask(i,j,k) * tmask(i + 1,j,    k) \\ vmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i    ,j + 1,k) \\ fmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i + 1,j,    k) \\ &  &* &tmask(i,j,k) * tmask(i + 1,j,    k) \\ wmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i    ,j,k - 1) \\ \text{with~} wmask(i,j,1) &= & &tmask(i,j,1) \end{alignat*} Note that, without ice shelves cavities, masks at $t-$ and $w-$points are identical with the numerical indexing used (\autoref{subsec:DOM_Num_Index}). Nevertheless, $wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface) exactly in the same way as for the bottom boundary. The specification of closed lateral boundaries requires that at least the first and last rows and columns of the \textit{mbathy} array are set to zero. In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to the second one and its first column equal to the last but one (and so too the mask arrays) (see \autoref{fig:LBC_jperio}). % ================================================================ % Domain: Initial State (dtatsd & istate) % ================================================================ \section{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} \label{sec:DTA_tsd} %-----------------------------------------namtsd------------------------------------------- \nlst{namtsd} %------------------------------------------------------------------------------------------ Options are defined in \ngn{namtsd}. By default, the ocean start from rest (the velocity field is set to zero) and the initialization of temperature and salinity fields is controlled through the \np{ln\_tsd\_ini} namelist parameter. The user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}) and to optionally decide on the fate of any freshwater imbalance over the area. The options are explained in \autoref{sec:MISC_closea} but it should be noted here that a successful use of these options requires appropriate mask fields to be present in the domain configuration file. Among the possibilities are: \begin{clines} int closea_mask     /* non-zero values in closed sea areas for optional masking                */ int closea_mask_rnf /* non-zero values in closed sea areas with runoff locations (precip only) */ int closea_mask_emp /* non-zero values in closed sea areas with runoff locations (total emp)   */ \end{clines} %% ================================================================================================= \subsection{Output grid files} \label{subsec:DOM_meshmask} Most of the arrays relating to a particular ocean model configuration discussed in this chapter (grid-point position, scale factors) can be saved in a file if namelist parameter \np{ln_write_cfg}{ln\_write\_cfg} (namelist \nam{cfg}{cfg}) is set to \forcode{.true.}; the output filename is set through parameter \np{cn_domcfg_out}{cn\_domcfg\_out}. This is only really useful if the fields are computed in subroutines \mdl{usrdef\_hgr} or \mdl{usrdef\_zgr} and checking or confirmation is required. Alternatively, all the arrays relating to a particular ocean model configuration (grid-point position, scale factors, depths and masks) can be saved in a file called \texttt{mesh\_mask} if namelist parameter \np{ln_meshmask}{ln\_meshmask} (namelist \nam{dom}{dom}) is set to \forcode{.true.}. This file contains additional fields that can be useful for post-processing applications. %% ================================================================================================= \section[Initial state (\textit{istate.F90} and \textit{dtatsd.F90})]{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} \label{sec:DOM_DTA_tsd} \begin{listing} \nlst{namtsd} \caption{\forcode{&namtsd}} \label{lst:namtsd} \end{listing} Basic initial state options are defined in \nam{tsd}{tsd}. By default, the ocean starts from rest (the velocity field is set to zero) and the initialization of temperature and salinity fields is controlled through the \np{ln_tsd_init}{ln\_tsd\_init} namelist parameter. \begin{description} \item[\np{ln\_tsd\_init}~\forcode{= .true.}] use a T and S input files that can be given on the model grid itself or on their native input data grid. \item [{\np[=.true.]{ln_tsd_init}{ln\_tsd\_init}}] Use T and S input files that can be given on the model grid itself or on their native input data grids. In the latter case, the data will be interpolated on-the-fly both in the horizontal and the vertical to the model grid (see \autoref{subsec:SBC_iof}). The information relative to the input files are given in the \np{sn\_tem} and \np{sn\_sal} structures. The information relating to the input files are specified in the \np{sn_tem}{sn\_tem} and \np{sn_sal}{sn\_sal} structures. The computation is done in the \mdl{dtatsd} module. \item[\np{ln\_tsd\_init}~\forcode{= .false.}] use constant salinity value of $35.5~psu$ and an analytical profile of temperature (typical of the tropical ocean), see \rou{istate\_t\_s} subroutine called from \mdl{istate} module. \item [{\np[=.false.]{ln_tsd_init}{ln\_tsd\_init}}] Initial values for T and S are set via a user supplied \rou{usr\_def\_istate} routine contained in \mdl{userdef\_istate}. The default version sets horizontally uniform T and profiles as used in the GYRE configuration (see \autoref{sec:CFGS_gyre}). \end{description} \biblio \pindex \subinc{\input{../../global/epilogue}} \end{document}
 r10499 \begin{document} % ================================================================ % Chapter ——— Ocean Dynamics (DYN) % ================================================================ \chapter{Ocean Dynamics (DYN)} \label{chap:DYN} \minitoc \thispagestyle{plain} \chaptertoc \paragraph{Changes record} ~\\ {\footnotesize \begin{tabularx}{\textwidth}{l||X|X} Release & Author(s) & Modifications \\ \hline {\em   4.0} & {\em ...} & {\em ...} \\ {\em   3.6} & {\em ...} & {\em ...} \\ {\em   3.4} & {\em ...} & {\em ...} \\ {\em <=3.4} & {\em ...} & {\em ...} \end{tabularx} } \clearpage Using the representation described in \autoref{chap:DOM}, (surface wind stress calculation using bulk formulae, estimation of mixing coefficients) that are carried out in modules SBC, LDF and ZDF and are described in \autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively. \autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively. In the present chapter we also describe the diagnostic equations used to compute the horizontal divergence, curl of the velocities (\emph{divcur} module) and the vertical velocity (\emph{wzvmod} module). The different options available to the user are managed by namelist variables. For term \textit{ttt} in the momentum equations, the logical namelist variables are \textit{ln\_dynttt\_xxx}, The different options available to the user are managed by namelist variables. For term \textit{ttt} in the momentum equations, the logical namelist variables are \textit{ln\_dynttt\_xxx}, where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme. If a CPP key is used for this term its name is \key{ttt}. %If a CPP key is used for this term its name is \key{ttt}. The corresponding code can be found in the \textit{dynttt\_xxx} module in the DYN directory, and it is usually computed in the \textit{dyn\_ttt\_xxx} subroutine. The user has the option of extracting and outputting each tendency term from the 3D momentum equations (\key{trddyn} defined), as described in \autoref{chap:MISC}. Furthermore, the tendency terms associated with the 2D barotropic vorticity balance (when \key{trdvor} is defined) (\texttt{trddyn?} defined), as described in \autoref{chap:MISC}. Furthermore, the tendency terms associated with the 2D barotropic vorticity balance (when \texttt{trdvor?} is defined) can be derived from the 3D terms. %%% \gmcomment{STEVEN: not quite sure I've got the sense of the last sentence. does MISC correspond to "extracting tendency terms" or "vorticity balance"?} % ================================================================ % Sea Surface Height evolution & Diagnostics variables % ================================================================ \cmtgm{STEVEN: not quite sure I've got the sense of the last sentence. Does MISC correspond to "extracting tendency terms" or "vorticity balance"?} %% ================================================================================================= \section{Sea surface height and diagnostic variables ($\eta$, $\zeta$, $\chi$, $w$)} \label{sec:DYN_divcur_wzv} %-------------------------------------------------------------------------------------------------------------- %           Horizontal divergence and relative vorticity %-------------------------------------------------------------------------------------------------------------- \subsection{Horizontal divergence and relative vorticity (\protect\mdl{divcur})} %% ================================================================================================= \subsection[Horizontal divergence and relative vorticity (\textit{divcur.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{divcur})} \label{subsec:DYN_divcur} The vorticity is defined at an $f$-point (\ie corner point) as follows: \label{eq:divcur_cur} The vorticity is defined at an $f$-point (\ie\ corner point) as follows: \label{eq:DYN_divcur_cur} \zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta_{i+1/2} \left[ {e_{2v}\;v} \right] -\delta_{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right) The horizontal divergence is defined at a $T$-point. It is given by: % \label{eq:divcur_div} % \label{eq:DYN_divcur_div} \chi =\frac{1}{e_{1t}\,e_{2t}\,e_{3t} } \left( {\delta_i \left[ {e_{2u}\,e_{3u}\,u} \right] ensure perfect restartability. The vorticity and divergence at the \textit{now} time step are used for the computation of the nonlinear advection and of the vertical velocity respectively. %-------------------------------------------------------------------------------------------------------------- % Sea Surface Height evolution %-------------------------------------------------------------------------------------------------------------- \subsection{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} the nonlinear advection and of the vertical velocity respectively. %% ================================================================================================= \subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} \label{subsec:DYN_sshwzv} The sea surface height is given by: \label{eq:dynspg_ssh} \label{eq:DYN_spg_ssh} \begin{aligned} \frac{\partial \eta }{\partial t} \end{aligned} where \textit{emp} is the surface freshwater budget (evaporation minus precipitation), where \textit{emp} is the surface freshwater budget (evaporation minus precipitation), expressed in Kg/m^2/s (which is equal to mm/s), and \rho_w=1,035~Kg/m^3 is the reference density of sea water (Boussinesq approximation). If river runoff is expressed as a surface freshwater flux (see \autoref{chap:SBC}) then \textit{emp} can be written as the evaporation minus precipitation, minus the river runoff. \textit{emp} can be written as the evaporation minus precipitation, minus the river runoff. The sea-surface height is evaluated using exactly the same time stepping scheme as the tracer equation \autoref{eq:tra_nxt}: the tracer equation \autoref{eq:TRA_nxt}: a leapfrog scheme in combination with an Asselin time filter, \ie the velocity appearing in \autoref{eq:dynspg_ssh} is centred in time (\textit{now} velocity). \ie\ the velocity appearing in \autoref{eq:DYN_spg_ssh} is centred in time (\textit{now} velocity). This is of paramount importance. Replacing T by the number 1 in the tracer equation and summing over the water column must lead to the sea surface height equation otherwise tracer content will not be conserved \citep{Griffies_al_MWR01, Leclair_Madec_OM09}. \citep{griffies.pacanowski.ea_MWR01, leclair.madec_OM09}. The vertical velocity is computed by an upward integration of the horizontal divergence starting at the bottom, taking into account the change of the thickness of the levels: \label{eq:wzv} \label{eq:DYN_wzv} \left\{ \begin{aligned} In the case of a non-linear free surface (\key{vvl}), the top vertical velocity is -\textit{emp}/\rho_w, In the case of a non-linear free surface (\texttt{vvl?}), the top vertical velocity is -\textit{emp}/\rho_w, as changes in the divergence of the barotropic transport are absorbed into the change of the level thicknesses, re-orientated downward. \gmcomment{not sure of this... to be modified with the change in emp setting} In the case of a linear free surface, the time derivative in \autoref{eq:wzv} disappears. \cmtgm{not sure of this... to be modified with the change in emp setting} In the case of a linear free surface, the time derivative in \autoref{eq:DYN_wzv} disappears. The upper boundary condition applies at a fixed level z=0. The top vertical velocity is thus equal to the divergence of the barotropic transport (\ie the first term in the right-hand-side of \autoref{eq:dynspg_ssh}). (\ie\ the first term in the right-hand-side of \autoref{eq:DYN_spg_ssh}). Note also that whereas the vertical velocity has the same discrete expression in z- and s-coordinates, its physical meaning is not the same: in the second case, w is the velocity normal to the s-surfaces. Note also that the k-axis is re-orientated downwards in the \fortran code compared to the indexing used in the semi-discrete equations such as \autoref{eq:wzv} (see \autoref{subsec:DOM_Num_Index_vertical}). % ================================================================ % Coriolis and Advection terms: vector invariant form % ================================================================ Note also that the k-axis is re-orientated downwards in the \fortran\ code compared to the indexing used in the semi-discrete equations such as \autoref{eq:DYN_wzv} (see \autoref{subsec:DOM_Num_Index_vertical}). %% ================================================================================================= \section{Coriolis and advection: vector invariant form} \label{sec:DYN_adv_cor_vect} %-----------------------------------------nam_dynadv---------------------------------------------------- \nlst{namdyn_adv} %------------------------------------------------------------------------------------------------------------- \begin{listing} \nlst{namdyn_adv} \caption{\forcode{&namdyn_adv}} \label{lst:namdyn_adv} \end{listing} The vector invariant form of the momentum equations is the one most often used in applications of the \NEMO ocean model. applications of the \NEMO\ ocean model. The flux form option (see next section) has been present since version 2. Options are defined through the \ngn{namdyn\_adv} namelist variables Coriolis and Options are defined through the \nam{dyn_adv}{dyn\_adv} namelist variables Coriolis and momentum advection terms are evaluated using a leapfrog scheme, \ie the velocity appearing in these expressions is centred in time (\textit{now} velocity). \ie\ the velocity appearing in these expressions is centred in time (\textit{now} velocity). At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied following \autoref{chap:LBC}. % ------------------------------------------------------------------------------------------------------------- % Vorticity term % ------------------------------------------------------------------------------------------------------------- \subsection{Vorticity term (\protect\mdl{dynvor})} %% ================================================================================================= \subsection[Vorticity term (\textit{dynvor.F90})]{Vorticity term (\protect\mdl{dynvor})} \label{subsec:DYN_vor} %------------------------------------------nam_dynvor---------------------------------------------------- \nlst{namdyn_vor} %------------------------------------------------------------------------------------------------------------- Options are defined through the \ngn{namdyn\_vor} namelist variables. Four discretisations of the vorticity term (\np{ln\_dynvor\_xxx}\forcode{ = .true.}) are available: \begin{listing} \nlst{namdyn_vor} \caption{\forcode{&namdyn_vor}} \label{lst:namdyn_vor} \end{listing} Options are defined through the \nam{dyn_vor}{dyn\_vor} namelist variables. Four discretisations of the vorticity term (\texttt{ln\_dynvor\_xxx}\forcode{=.true.}) are available: conserving potential enstrophy of horizontally non-divergent flow (ENS scheme); conserving horizontal kinetic energy (ENE scheme); horizontal kinetic energy for the planetary vorticity term (MIX scheme); or conserving both the potential enstrophy of horizontally non-divergent flow and horizontal kinetic energy (EEN scheme) (see \autoref{subsec:C_vorEEN}). (EEN scheme) (see \autoref{subsec:INVARIANTS_vorEEN}). In the case of ENS, ENE or MIX schemes the land sea mask may be slightly modified to ensure the consistency of vorticity term with analytical equations (\np{ln\_dynvor\_con}\forcode{ = .true.}). vorticity term with analytical equations (\np[=.true.]{ln_dynvor_con}{ln\_dynvor\_con}). The vorticity terms are all computed in dedicated routines that can be found in the \mdl{dynvor} module. %------------------------------------------------------------- % enstrophy conserving scheme %------------------------------------------------------------- \subsubsection{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})} %% ================================================================================================= \subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens})]{Enstrophy conserving scheme (\protect\np{ln_dynvor_ens}{ln\_dynvor\_ens})} \label{subsec:DYN_vor_ens} In the enstrophy conserving case (ENS scheme), the discrete formulation of the vorticity term provides a global conservation of the enstrophy ( [ (\zeta +f ) / e_{3f} ]^2  in s-coordinates) for a horizontally non-divergent flow (\ie \chi=0), ( [ (\zeta +f ) / e_{3f} ]^2  in s-coordinates) for a horizontally non-divergent flow (\ie\ \chi=0), but does not conserve the total kinetic energy. It is given by: \label{eq:dynvor_ens} \label{eq:DYN_vor_ens} \left\{ \begin{aligned} \end{aligned} \right. %------------------------------------------------------------- % energy conserving scheme %------------------------------------------------------------- \subsubsection{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})} %% ================================================================================================= \subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene})]{Energy conserving scheme (\protect\np{ln_dynvor_ene}{ln\_dynvor\_ene})} \label{subsec:DYN_vor_ene} It is given by: \label{eq:dynvor_ene} \label{eq:DYN_vor_ene} \left\{ \begin{aligned} \end{aligned} \right. %------------------------------------------------------------- % mix energy/enstrophy conserving scheme %------------------------------------------------------------- \subsubsection{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.}) } %% ================================================================================================= \subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix})]{Mixed energy/enstrophy conserving scheme (\protect\np{ln_dynvor_mix}{ln\_dynvor\_mix})} \label{subsec:DYN_vor_mix} For the mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the two previous schemes is used. It consists of the ENS scheme (\autoref{eq:dynvor_ens}) for the relative vorticity term, and of the ENE scheme (\autoref{eq:dynvor_ene}) applied to the planetary vorticity term. \[ % \label{eq:dynvor_mix} It consists of the ENS scheme (\autoref{eq:DYN_vor_ens}) for the relative vorticity term, and of the ENE scheme (\autoref{eq:DYN_vor_ene}) applied to the planetary vorticity term. \[ % \label{eq:DYN_vor_mix} \left\{ { \begin{aligned} %------------------------------------------------------------- %                 energy and enstrophy conserving scheme %------------------------------------------------------------- \subsubsection{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.}) } %% ================================================================================================= \subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een})]{Energy and enstrophy conserving scheme (\protect\np{ln_dynvor_een}{ln\_dynvor\_een})} \label{subsec:DYN_vor_een} the presence of grid point oscillation structures that will be invisible to the operator. These structures are \textit{computational modes} that will be at least partly damped by the momentum diffusion operator (\ie the subgrid-scale advection), but not by the resolved advection term. the momentum diffusion operator (\ie\ the subgrid-scale advection), but not by the resolved advection term. The ENS and ENE schemes therefore do not contribute to dump any grid point noise in the horizontal velocity field. Such noise would result in more noise in the vertical velocity field, an undesirable feature. $u$ and $v$ are located at different grid points, a price worth paying to avoid a double averaging in the pressure gradient term as in the $B$-grid. \gmcomment{ To circumvent this, Adcroft (ADD REF HERE) \cmtgm{ To circumvent this, Adcroft (ADD REF HERE) Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves....} A very nice solution to the problem of double averaging was proposed by \citet{Arakawa_Hsu_MWR90}. A very nice solution to the problem of double averaging was proposed by \citet{arakawa.hsu_MWR90}. The idea is to get rid of the double averaging by considering triad combinations of vorticity. It is noteworthy that this solution is conceptually quite similar to the one proposed by \citep{Griffies_al_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}). The \citet{Arakawa_Hsu_MWR90} vorticity advection scheme for a single layer is modified for spherical coordinates as described by \citet{Arakawa_Lamb_MWR81} to obtain the EEN scheme. First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point: $% \label{eq:pot_vor} \citep{griffies.gnanadesikan.ea_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:INVARIANTS}). The \citet{arakawa.hsu_MWR90} vorticity advection scheme for a single layer is modified for spherical coordinates as described by \citet{arakawa.lamb_MWR81} to obtain the EEN scheme. First consider the discrete expression of the potential vorticity, q, defined at an f-point: \[ % \label{eq:DYN_pot_vor} q = \frac{\zeta +f} {e_{3f} }$ where the relative vorticity is defined by (\autoref{eq:divcur_cur}), the Coriolis parameter is given by $f=2 \,\Omega \;\sin \varphi _f$ and the layer thickness at $f$-points is: \label{eq:een_e3f} where the relative vorticity is defined by (\autoref{eq:DYN_divcur_cur}), the Coriolis parameter is given by $f=2 \,\Omega \;\sin \varphi _f$ and the layer thickness at $f$-points is: \label{eq:DYN_een_e3f} e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2} %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> \begin{figure}[!ht] \begin{center} \includegraphics[width=0.70\textwidth]{Fig_DYN_een_triad} \caption{ \protect\label{fig:DYN_een_triad} Triads used in the energy and enstrophy conserving scheme (een) for $u$-component (upper panel) and $v$-component (lower panel). } \end{center} \centering \includegraphics[width=0.66\textwidth]{DYN_een_triad} \caption[Triads used in the energy and enstrophy conserving scheme (EEN)]{ Triads used in the energy and enstrophy conserving scheme (EEN) for $u$-component (upper panel) and $v$-component (lower panel).} \label{fig:DYN_een_triad} \end{figure} % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> A key point in \autoref{eq:een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made. A key point in \autoref{eq:DYN_een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made. It uses the sum of masked t-point vertical scale factor divided either by the sum of the four t-point masks (\np{nn\_een\_e3f}\forcode{ = 1}), or just by $4$ (\np{nn\_een\_e3f}\forcode{ = .true.}). (\np[=1]{nn_een_e3f}{nn\_een\_e3f}), or just by $4$ (\np[=.true.]{nn_een_e3f}{nn\_een\_e3f}). The latter case preserves the continuity of $e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and extends by continuity the value of $e_{3f}$ into the land areas. (with a systematic reduction of $e_{3f}$ when a model level intercept the bathymetry) that tends to reinforce the topostrophy of the flow (\ie the tendency of the flow to follow the isobaths) \citep{Penduff_al_OS07}. (\ie\ the tendency of the flow to follow the isobaths) \citep{penduff.le-sommer.ea_OS07}. Next, the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as the following triad combinations of the neighbouring potential vorticities defined at f-points (\autoref{fig:DYN_een_triad}): \label{eq:Q_triads} (\autoref{fig:DYN_een_triad}): \label{eq:DYN_Q_triads} _i^j \mathbb{Q}^{i_p}_{j_p} = \frac{1}{12} \ \left(   q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p}  \right) where the indices $i_p$ and $k_p$ take the values: $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$. Finally, the vorticity terms are represented as: \label{eq:dynvor_een} where the indices $i_p$ and $k_p$ take the values: $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$. Finally, the vorticity terms are represented as: \label{eq:DYN_vor_een} \left\{ { \begin{aligned} \end{aligned} } \right. This EEN scheme in fact combines the conservation properties of the ENS and ENE schemes. It conserves both total energy and potential enstrophy in the limit of horizontally nondivergent flow (\ie $\chi$=$0$) (see \autoref{subsec:C_vorEEN}). (\ie\ $\chi$=$0$) (see \autoref{subsec:INVARIANTS_vorEEN}). Applied to a realistic ocean configuration, it has been shown that it leads to a significant reduction of the noise in the vertical velocity field \citep{Le_Sommer_al_OM09}. the noise in the vertical velocity field \citep{le-sommer.penduff.ea_OM09}. Furthermore, used in combination with a partial steps representation of bottom topography, it improves the interaction between current and topography, leading to a larger topostrophy of the flow \citep{Barnier_al_OD06, Penduff_al_OS07}. %-------------------------------------------------------------------------------------------------------------- %           Kinetic Energy Gradient term %-------------------------------------------------------------------------------------------------------------- \subsection{Kinetic energy gradient term (\protect\mdl{dynkeg})} leading to a larger topostrophy of the flow \citep{barnier.madec.ea_OD06, penduff.le-sommer.ea_OS07}. %% ================================================================================================= \subsection[Kinetic energy gradient term (\textit{dynkeg.F90})]{Kinetic energy gradient term (\protect\mdl{dynkeg})} \label{subsec:DYN_keg} As demonstrated in \autoref{apdx:C}, As demonstrated in \autoref{apdx:INVARIANTS}, there is a single discrete formulation of the kinetic energy gradient term that, together with the formulation chosen for the vertical advection (see below), conserves the total kinetic energy: % \label{eq:dynkeg} % \label{eq:DYN_keg} \left\{ \begin{aligned} %-------------------------------------------------------------------------------------------------------------- %           Vertical advection term %-------------------------------------------------------------------------------------------------------------- \subsection{Vertical advection term (\protect\mdl{dynzad}) } %% ================================================================================================= \subsection[Vertical advection term (\textit{dynzad.F90})]{Vertical advection term (\protect\mdl{dynzad})} \label{subsec:DYN_zad} conserves the total kinetic energy. Indeed, the change of KE due to the vertical advection is exactly balanced by the change of KE due to the gradient of KE (see \autoref{apdx:C}). % \label{eq:dynzad} the change of KE due to the gradient of KE (see \autoref{apdx:INVARIANTS}). \[ % \label{eq:DYN_zad} \left\{ \begin{aligned} \right. When \np{ln\_dynzad\_zts}\forcode{ = .true.}, When \np[=.true.]{ln_dynzad_zts}{ln\_dynzad\_zts}, a split-explicit time stepping with 5 sub-timesteps is used on the vertical advection term. This option can be useful when the value of the timestep is limited by vertical advection \citep{Lemarie_OM2015}. This option can be useful when the value of the timestep is limited by vertical advection \citep{lemarie.debreu.ea_OM15}. Note that in this case, a similar split-explicit time stepping should be used on vertical advection of tracer to ensure a better stability, an option which is only available with a TVD scheme (see \np{ln\_traadv\_tvd\_zts} in \autoref{subsec:TRA_adv_tvd}). % ================================================================ % Coriolis and Advection : flux form % ================================================================ an option which is only available with a TVD scheme (see \np{ln_traadv_tvd_zts}{ln\_traadv\_tvd\_zts} in \autoref{subsec:TRA_adv_tvd}). %% ================================================================================================= \section{Coriolis and advection: flux form} \label{sec:DYN_adv_cor_flux} %------------------------------------------nam_dynadv---------------------------------------------------- \nlst{namdyn_adv} %------------------------------------------------------------------------------------------------------------- Options are defined through the \ngn{namdyn\_adv} namelist variables. Options are defined through the \nam{dyn_adv}{dyn\_adv} namelist variables. In the flux form (as in the vector invariant form), the Coriolis and momentum advection terms are evaluated using a leapfrog scheme, \ie the velocity appearing in their expressions is centred in time (\textit{now} velocity). \ie\ the velocity appearing in their expressions is centred in time (\textit{now} velocity). At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied following \autoref{chap:LBC}. %-------------------------------------------------------------------------------------------------------------- %           Coriolis plus curvature metric terms %-------------------------------------------------------------------------------------------------------------- \subsection{Coriolis plus curvature metric terms (\protect\mdl{dynvor}) } %% ================================================================================================= \subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})]{Coriolis plus curvature metric terms (\protect\mdl{dynvor})} \label{subsec:DYN_cor_flux} In flux form, the vorticity term reduces to a Coriolis term in which the Coriolis parameter has been modified to account for the "metric" term. This altered Coriolis parameter is thus discretised at $f$-points. It is given by: It is given by: \begin{multline*} % \label{eq:dyncor_metric} % \label{eq:DYN_cor_metric} f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i}  -  u\frac{\partial e_1 }{\partial j}} \right)  \\ \equiv   f + \frac{1}{e_{1f} e_{2f} } \left( { \ \overline v ^{i+1/2}\delta_{i+1/2} \left[ {e_{2u} } \right] -  \overline u ^{j+1/2}\delta_{j+1/2} \left[ {e_{1u} } \right]  }  \ \right) \end{multline*} Any of the (\autoref{eq:dynvor_ens}), (\autoref{eq:dynvor_ene}) and (\autoref{eq:dynvor_een}) schemes can be used to \end{multline*} Any of the (\autoref{eq:DYN_vor_ens}), (\autoref{eq:DYN_vor_ene}) and (\autoref{eq:DYN_vor_een}) schemes can be used to compute the product of the Coriolis parameter and the vorticity. However, the energy-conserving scheme (\autoref{eq:dynvor_een}) has exclusively been used to date. This term is evaluated using a leapfrog scheme, \ie the velocity is centred in time (\textit{now} velocity). %-------------------------------------------------------------------------------------------------------------- %           Flux form Advection term %-------------------------------------------------------------------------------------------------------------- \subsection{Flux form advection term (\protect\mdl{dynadv}) } However, the energy-conserving scheme (\autoref{eq:DYN_vor_een}) has exclusively been used to date. This term is evaluated using a leapfrog scheme, \ie\ the velocity is centred in time (\textit{now} velocity). %% ================================================================================================= \subsection[Flux form advection term (\textit{dynadv.F90})]{Flux form advection term (\protect\mdl{dynadv})} \label{subsec:DYN_adv_flux} The discrete expression of the advection term is given by: \[ % \label{eq:dynadv} % \label{eq:DYN_adv} \left\{ \begin{aligned} a $2^{nd}$ order centered finite difference scheme, CEN2, or a $3^{rd}$ order upstream biased scheme, UBS. The latter is described in \citet{Shchepetkin_McWilliams_OM05}. The schemes are selected using the namelist logicals \np{ln\_dynadv\_cen2} and \np{ln\_dynadv\_ubs}. The latter is described in \citet{shchepetkin.mcwilliams_OM05}. The schemes are selected using the namelist logicals \np{ln_dynadv_cen2}{ln\_dynadv\_cen2} and \np{ln_dynadv_ubs}{ln\_dynadv\_ubs}. In flux form, the schemes differ by the choice of a space and time interpolation to define the value of $u$ and $v$ at the centre of each face of $u$- and $v$-cells, \ie at the $T$-, $f$-, and $uw$-points for $u$ and at the $f$-, $T$- and $vw$-points for $v$. %------------------------------------------------------------- $u$ and $v$ at the centre of each face of $u$- and $v$-cells, \ie\ at the $T$-, $f$-, and $uw$-points for $u$ and at the $f$-, $T$- and $vw$-points for $v$. %                 2nd order centred scheme %------------------------------------------------------------- \subsubsection{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})} %% ================================================================================================= \subsubsection[CEN2: $2^{nd}$ order centred scheme (\forcode{ln_dynadv_cen2})]{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln_dynadv_cen2}{ln\_dynadv\_cen2})} \label{subsec:DYN_adv_cen2} In the centered $2^{nd}$ order formulation, the velocity is evaluated as the mean of the two neighbouring points: \label{eq:dynadv_cen2} \label{eq:DYN_adv_cen2} \left\{ \begin{aligned} \end{aligned} \right. The scheme is non diffusive (\ie conserves the kinetic energy) but dispersive (\ie it may create false extrema). The scheme is non diffusive (\ie\ conserves the kinetic energy) but dispersive (\ie\ it may create false extrema). It is therefore notoriously noisy and must be used in conjunction with an explicit diffusion operator to produce a sensible solution. so $u$ and $v$ are the \emph{now} velocities. %------------------------------------------------------------- %                 UBS scheme %------------------------------------------------------------- \subsubsection{UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})} %% ================================================================================================= \subsubsection[UBS: Upstream Biased Scheme (\forcode{ln_dynadv_ubs})]{UBS: Upstream Biased Scheme (\protect\np{ln_dynadv_ubs}{ln\_dynadv\_ubs})} \label{subsec:DYN_adv_ubs} For example, the evaluation of $u_T^{ubs}$ is done as follows: \label{eq:dynadv_ubs} \label{eq:DYN_adv_ubs} u_T^{ubs} =\overline u ^i-\;\frac{1}{6} \begin{cases} where $u"_{i+1/2} =\delta_{i+1/2} \left[ {\delta_i \left[ u \right]} \right]$. This results in a dissipatively dominant (\ie hyper-diffusive) truncation error \citep{Shchepetkin_McWilliams_OM05}. The overall performance of the advection scheme is similar to that reported in \citet{Farrow1995}. This results in a dissipatively dominant (\ie\ hyper-diffusive) truncation error \citep{shchepetkin.mcwilliams_OM05}. The overall performance of the advection scheme is similar to that reported in \citet{farrow.stevens_JPO95}. It is a relatively good compromise between accuracy and smoothness. It is not a \emph{positive} scheme, meaning that false extrema are permitted. But the amplitudes of the false extrema are significantly reduced over those in the centred second order method. As the scheme already includes a diffusion component, it can be used without explicit lateral diffusion on momentum (\ie \np{ln\_dynldf\_lap}\forcode{ = }\np{ln\_dynldf\_bilap}\forcode{ = .false.}), As the scheme already includes a diffusion component, it can be used without explicit lateral diffusion on momentum (\ie\ \np[=]{ln_dynldf_lap}{ln\_dynldf\_lap}\np[=.false.]{ln_dynldf_bilap}{ln\_dynldf\_bilap}), and it is recommended to do so. The UBS scheme is not used in all directions. In the vertical, the centred $2^{nd}$ order evaluation of the advection is preferred, \ie $u_{uw}^{ubs}$ and $u_{vw}^{ubs}$ in \autoref{eq:dynadv_cen2} are used. UBS is diffusive and is associated with vertical mixing of momentum. \gmcomment{ gm  pursue the In the vertical, the centred $2^{nd}$ order evaluation of the advection is preferred, \ie\ $u_{uw}^{ubs}$ and $u_{vw}^{ubs}$ in \autoref{eq:DYN_adv_cen2} are used. UBS is diffusive and is associated with vertical mixing of momentum. \cmtgm{ gm  pursue the sentence:Since vertical mixing of momentum is a source term of the TKE equation...  } For stability reasons, the first term in (\autoref{eq:dynadv_ubs}), For stability reasons, the first term in (\autoref{eq:DYN_adv_ubs}), which corresponds to a second order centred scheme, is evaluated using the \textit{now} velocity (centred in time), while the second term, which is the diffusion part of the scheme, is evaluated using the \textit{before} velocity (forward in time). This is discussed by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme. This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme. Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) schemes only differ by one coefficient. Replacing $1/6$ by $1/8$ in (\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. Replacing $1/6$ by $1/8$ in (\autoref{eq:DYN_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}. This option is not available through a namelist parameter, since the $1/6$ coefficient is hard coded. Nevertheless it is quite easy to make the substitution in the \mdl{dynadv\_ubs} module and obtain a QUICK scheme. there is also the possibility of using a $4^{th}$ order evaluation of the advective velocity as in ROMS. This is an error and should be suppressed soon. %%%