Changeset 11422


Ignore:
Timestamp:
2019-08-08T15:40:47+02:00 (14 months ago)
Author:
jchanut
Message:

#1791, merge with trunk

Location:
NEMO/branches/2019/fix_vvl_ticket1791
Files:
33 deleted
139 edited
32 copied

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/fix_vvl_ticket1791/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg

    r10075 r11422  
    2222   cn_exp      = "Pacific" !  experience name 
    2323   nn_it000    =       1   !  first time step 
    24    nn_itend    =      75   !  last  time step (std 5475) 
     24   nn_itend    =      80   !  last  time step (std 5475) 
    2525   nn_istate   =       0   !  output the initial state (1) or not (0) 
    2626/ 
     
    3030   ln_linssh   = .false.   !  =T  linear free surface  ==>>  model level are fixed in time 
    3131   ! 
    32    rn_rdt      = 5760.     !  time step for the dynamics and tracer 
     32   rn_rdt      = 5400.     !  time step for the dynamics and tracer 
    3333/ 
    3434!----------------------------------------------------------------------- 
  • NEMO/branches/2019/fix_vvl_ticket1791/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg

    r10072 r11422  
    2222   cn_exp      = "Nordic"  !  experience name  
    2323   nn_it000    =    1      !  first time step 
    24    nn_itend    =  300      !  last  time step 
     24   nn_itend    =  320      !  last  time step 
    2525   nn_istate   =      0    !  output the initial state (1) or not (0) 
    2626   ln_clobber  = .true.    !  clobber (overwrite) an existing file 
     
    3131   ln_linssh   = .false.   !  =T  linear free surface  ==>>  model level are fixed in time 
    3232   ! 
    33    rn_rdt      = 1440.     !  time step for the dynamics (and tracer if nn_acc=0) 
     33   rn_rdt      = 1350.     !  time step for the dynamics (and tracer if nn_acc=0) 
    3434/ 
    3535!----------------------------------------------------------------------- 
  • NEMO/branches/2019/fix_vvl_ticket1791/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg

    r10072 r11422  
    2222   cn_exp      = "Nordic"  !  experience name  
    2323   nn_it000    =    1      !  first time step 
    24    nn_itend    =  900      !  last  time step 
     24   nn_itend    =  960      !  last  time step 
    2525   nn_istate   =      0    !  output the initial state (1) or not (0) 
    2626   ln_clobber  = .true.    !  clobber (overwrite) an existing file 
     
    3131   ln_linssh   = .false.   !  =T  linear free surface  ==>>  model level are fixed in time 
    3232   ! 
    33    rn_rdt      = 480.     !  time step for the dynamics (and tracer if nn_acc=0) 
     33   rn_rdt      = 450.     !  time step for the dynamics (and tracer if nn_acc=0) 
    3434/ 
    3535!----------------------------------------------------------------------- 
  • NEMO/branches/2019/fix_vvl_ticket1791/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg

    r10075 r11422  
    2222   cn_exp      =  "ORCA2"  !  experience name 
    2323   nn_it000    =       1   !  first time step 
    24    nn_itend    =      75   !  last  time step (std 5475) 
     24   nn_itend    =      80   !  last  time step (std 5475) 
    2525   nn_istate   =       0   !  output the initial state (1) or not (0) 
    2626/ 
     
    3030   ln_linssh   = .false.   !  =T  linear free surface  ==>>  model level are fixed in time 
    3131   ! 
    32    rn_rdt      = 5760.     !  time step for the dynamics and tracer 
     32   rn_rdt      = 5400.     !  time step for the dynamics and tracer 
    3333/ 
    3434!----------------------------------------------------------------------- 
  • NEMO/branches/2019/fix_vvl_ticket1791/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg

    r10190 r11422  
    2222   cn_exp      =  "ORCA2"  !  experience name 
    2323   nn_it000    =       1   !  first time step 
    24    nn_itend    =    5475   !  last  time step (std 5475) 
     24   nn_itend    =    5840   !  last  time step (std 5475) 
    2525   nn_istate   =       0   !  output the initial state (1) or not (0) 
    2626/ 
     
    2828&namdom        !   time and space domain 
    2929!----------------------------------------------------------------------- 
    30    rn_rdt      = 5760.     !  time step for the dynamics and tracer 
     30   rn_rdt      = 5400.     !  time step for the dynamics and tracer 
    3131/ 
    3232!----------------------------------------------------------------------- 
     
    7575&namsbc        !   Surface Boundary Condition manager                   (default: NO selection) 
    7676!----------------------------------------------------------------------- 
    77    nn_fsbc     = 3         !  frequency of SBC module call 
     77   nn_fsbc     = 4         !  frequency of SBC module call 
    7878                           !     (also = the frequency of sea-ice & iceberg model call) 
    7979                     ! Type of air-sea fluxes  
  • NEMO/branches/2019/fix_vvl_ticket1791/cfgs/SPITZ12/EXPREF/namelist_ice_cfg

    r10911 r11422  
    102102&namdia         !   Diagnostics 
    103103!------------------------------------------------------------------------------ 
     104   ln_icediachk     = .true.         !  check online the heat, mass & salt budgets (T) or not (F) 
    104105/ 
  • NEMO/branches/2019/fix_vvl_ticket1791/doc

    • Property svn:ignore deleted
  • NEMO/branches/2019/fix_vvl_ticket1791/doc/latex

    • Property svn:ignore
      •  

        old new  
        1 *.aux 
        2 *.bbl 
        3 *.blg 
        4 *.dvi 
        5 *.fdb* 
        6 *.fls 
        7 *.idx 
        8 *.ilg 
        9 *.ind 
        10 *.log 
        11 *.maf 
        12 *.mtc* 
        13 *.out 
        14 *.pdf 
        15 *.toc 
        16 _minted-* 
         1figures 
  • NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/HTML_latex2html.sh

    r10146 r11422  
    11#!/bin/bash 
    22 
    3 ./inc/clean.sh 
    4 ./inc/build.sh 
     3#./inc/clean.sh 
     4#./inc/build.sh 
    55 
    6 cd tex_main 
    7 sed -i -e 's#\\documentclass#%\\documentclass#' -e '/{document}/ s/^/%/' \ 
    8    ../tex_sub/*.tex 
    9 sed -i    's#\\subfile{#\\include{#g' \ 
    10    NEMO_manual.tex 
    11 latex2html -local_icons -no_footnode -split 4 -link 2 -mkdir -dir ../html_LaTeX2HTML   \ 
    12             $*                                                                         \ 
    13    NEMO_manual 
    14 sed -i -e 's#%\\documentclass#\\documentclass#' -e '/{document}/ s/^%//' \ 
    15    ../tex_sub/*.tex 
    16 sed -i    's#\\include{#\\subfile{#g' \ 
    17    NEMO_manual.tex 
     6sed -i -e 's#utf8#latin1#'                                 \ 
     7       -e 's#\[outputdir=../build\]{minted}#\[\]{minted}#' \ 
     8       -e '/graphicspath/ s#{../#{../../#g'                \ 
     9          global/packages.tex 
     10 
     11cd ./NEMO/main 
     12sed -i -e 's#\\documentclass#%\\documentclass#' -e '/{document}/ s#^#%#'   ../subfiles/*.tex 
     13sed -i    's#\\subfile{#\\input{#'                                         chapters.tex appendices.tex 
     14 
     15#latex2html                                   -noimages -local_icons -no_footnode -split 4 -link 2 -dir ../html_LaTeX2HTML $* NEMO_manual 
     16latex2html -debug -noreuse -init_file ../../l2hconf.pm -local_icons                               -dir ../build/html               NEMO_manual 
     17 
     18sed -i -e 's#%\\documentclass#\\documentclass#' -e '/{document}/ s#^%##'   ../subfiles/*.tex 
     19sed -i    's#\\input{#\\subfile{#'                                         chapters.tex appendices.tex 
    1820cd - 
    1921 
     22sed -i -e 's#latin1#utf8#'                                 \ 
     23       -e 's#\[\]{minted}#\[outputdir=../build\]{minted}#' \ 
     24       -e '/graphicspath/ s#{../../#{../#g'                \ 
     25     global/packages.tex 
     26 
    2027exit 0 
  • NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO

    • Property svn:ignore deleted
  • NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/main

    • Property svn:ignore
      •  

        old new  
        1 *.aux 
        2 *.bbl 
        3 *.blg 
        4 *.dvi 
        5 *.fdb* 
        6 *.fls 
        7 *.idx 
        81*.ilg 
        92*.ind 
        10 *.log 
        11 *.maf 
        12 *.mtc* 
        13 *.out 
        14 *.pdf 
        15 *.toc 
        16 _minted-* 
    • Property svn:externals set to
      ^/utils/dev/bibtool.rsc bibtool.rsc
  • NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/annex_A.tex

    r10442 r11422  
    1010 
    1111\minitoc 
     12 
     13\vfill 
     14\begin{figure}[b] 
     15\subsubsection*{Changes record} 
     16\begin{tabular}{l||l|m{0.65\linewidth}} 
     17    Release   & Author        & Modifications \\ 
     18    {\em 4.0} & {\em Mike Bell} & {\em review}  \\ 
     19    {\em 3.x} & {\em Gurvan Madec} & {\em original}  \\ 
     20\end{tabular} 
     21\end{figure} 
     22 
    1223 
    1324\newpage 
     
    2940\begin{equation} 
    3041  \label{apdx:A_s_slope} 
    31   \sigma_1 =\frac{1}{e_1 }\;\left. {\frac{\partial z}{\partial i}} \right|_s 
     42  \sigma_1 =\frac{1}{e_1 } \; \left. {\frac{\partial z}{\partial i}} \right|_s 
    3243  \quad \text{and} \quad 
    33   \sigma_2 =\frac{1}{e_2 }\;\left. {\frac{\partial z}{\partial j}} \right|_s 
    34 \end{equation} 
    35  
    36 The chain rule to establish the model equations in the curvilinear $s-$coordinate system is: 
     44  \sigma_2 =\frac{1}{e_2 } \; \left. {\frac{\partial z}{\partial j}} \right|_s . 
     45\end{equation} 
     46 
     47The model fields (e.g. pressure $p$) can be viewed as functions of $(i,j,z,t)$ (e.g. $p(i,j,z,t)$) or as 
     48functions of $(i,j,s,t)$ (e.g. $p(i,j,s,t)$). The symbol $\bullet$ will be used to represent any one of  
     49these fields.  Any ``infinitesimal'' change in $\bullet$ can be written in two forms:  
     50\begin{equation} 
     51  \label{apdx:A_s_infin_changes} 
     52  \begin{aligned} 
     53    & \delta \bullet =  \delta i \left. \frac{ \partial \bullet }{\partial i} \right|_{j,s,t}  
     54                + \delta j \left. \frac{ \partial \bullet }{\partial i} \right|_{i,s,t}  
     55                + \delta s \left. \frac{ \partial \bullet }{\partial s} \right|_{i,j,t}  
     56                + \delta t \left. \frac{ \partial \bullet }{\partial t} \right|_{i,j,s} , \\ 
     57    & \delta \bullet =  \delta i \left. \frac{ \partial \bullet }{\partial i} \right|_{j,z,t}  
     58                + \delta j \left. \frac{ \partial \bullet }{\partial i} \right|_{i,z,t}  
     59                + \delta z \left. \frac{ \partial \bullet }{\partial z} \right|_{i,j,t}  
     60                + \delta t \left. \frac{ \partial \bullet }{\partial t} \right|_{i,j,z} . 
     61  \end{aligned} 
     62\end{equation} 
     63Using the first form and considering a change $\delta i$ with $j, z$ and $t$ held constant, shows that 
     64\begin{equation} 
     65  \label{apdx:A_s_chain_rule} 
     66      \left. {\frac{\partial \bullet }{\partial i}} \right|_{j,z,t}  = 
     67      \left. {\frac{\partial \bullet }{\partial i}} \right|_{j,s,t} 
     68    + \left. {\frac{\partial s       }{\partial i}} \right|_{j,z,t} \;  
     69      \left. {\frac{\partial \bullet }{\partial s}} \right|_{i,j,t} .     
     70\end{equation} 
     71The term $\left. \partial s / \partial i \right|_{j,z,t}$ can be related to the slope of constant $s$ surfaces,  
     72(\autoref{apdx:A_s_slope}), by applying the second of (\autoref{apdx:A_s_infin_changes}) with $\bullet$ set to  
     73$s$ and $j, t$ held constant 
     74\begin{equation} 
     75\label{apdx:a_delta_s} 
     76\delta s|_{j,t} =  
     77         \delta i \left. \frac{ \partial s }{\partial i} \right|_{j,z,t}  
     78       + \delta z \left. \frac{ \partial s }{\partial z} \right|_{i,j,t} . 
     79\end{equation} 
     80Choosing to look at a direction in the $(i,z)$ plane in which $\delta s = 0$ and using 
     81(\autoref{apdx:A_s_slope}) we obtain  
     82\begin{equation} 
     83\left. \frac{ \partial s }{\partial i} \right|_{j,z,t} =   
     84         -  \left. \frac{ \partial z }{\partial i} \right|_{j,s,t} \; 
     85            \left. \frac{ \partial s }{\partial z} \right|_{i,j,t} 
     86    = - \frac{e_1 }{e_3 }\sigma_1  . 
     87\label{apdx:a_ds_di_z} 
     88\end{equation} 
     89Another identity, similar in form to (\autoref{apdx:a_ds_di_z}), can be derived 
     90by choosing $\bullet$ to be $s$ and using the second form of (\autoref{apdx:A_s_infin_changes}) to consider 
     91changes in which $i , j$ and $s$ are constant. This shows that 
     92\begin{equation} 
     93\label{apdx:A_w_in_s} 
     94w_s = \left. \frac{ \partial z }{\partial t} \right|_{i,j,s} =   
     95- \left. \frac{ \partial z }{\partial s} \right|_{i,j,t} 
     96  \left. \frac{ \partial s }{\partial t} \right|_{i,j,z}  
     97  = - e_3 \left. \frac{ \partial s }{\partial t} \right|_{i,j,z} .  
     98\end{equation} 
     99 
     100In what follows, for brevity, indication of the constancy of the $i, j$ and $t$ indices is  
     101usually omitted. Using the arguments outlined above one can show that the chain rules needed to establish  
     102the model equations in the curvilinear $s-$coordinate system are: 
    37103\begin{equation} 
    38104  \label{apdx:A_s_chain_rule} 
    39105  \begin{aligned} 
    40106    &\left. {\frac{\partial \bullet }{\partial t}} \right|_z  = 
    41     \left. {\frac{\partial \bullet }{\partial t}} \right|_s 
    42     -\frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial t} \\ 
     107    \left. {\frac{\partial \bullet }{\partial t}} \right|_s  
     108    + \frac{\partial \bullet }{\partial s}\; \frac{\partial s}{\partial t} , \\ 
    43109    &\left. {\frac{\partial \bullet }{\partial i}} \right|_z  = 
    44110    \left. {\frac{\partial \bullet }{\partial i}} \right|_s 
    45     -\frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial i}= 
    46     \left. {\frac{\partial \bullet }{\partial i}} \right|_s 
    47     -\frac{e_1 }{e_3 }\sigma_1 \frac{\partial \bullet }{\partial s} \\ 
     111    +\frac{\partial \bullet }{\partial s}\; \frac{\partial s}{\partial i}= 
     112    \left. {\frac{\partial \bullet }{\partial i}} \right|_s  
     113    -\frac{e_1 }{e_3 }\sigma_1 \frac{\partial \bullet }{\partial s} , \\ 
    48114    &\left. {\frac{\partial \bullet }{\partial j}} \right|_z  = 
    49     \left. {\frac{\partial \bullet }{\partial j}} \right|_s 
    50     - \frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial j}= 
    51     \left. {\frac{\partial \bullet }{\partial j}} \right|_s 
    52     - \frac{e_2 }{e_3 }\sigma_2 \frac{\partial \bullet }{\partial s} \\ 
    53     &\;\frac{\partial \bullet }{\partial z}  \;\; = \frac{1}{e_3 }\frac{\partial \bullet }{\partial s} 
     115    \left. {\frac{\partial \bullet }{\partial j}} \right|_s  
     116    + \frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial j}= 
     117    \left. {\frac{\partial \bullet }{\partial j}} \right|_s  
     118    - \frac{e_2 }{e_3 }\sigma_2 \frac{\partial \bullet }{\partial s} , \\ 
     119    &\;\frac{\partial \bullet }{\partial z}  \;\; = \frac{1}{e_3 }\frac{\partial \bullet }{\partial s} . 
    54120  \end{aligned} 
    55121\end{equation} 
    56122 
    57 In particular applying the time derivative chain rule to $z$ provides the expression for $w_s$, 
    58 the vertical velocity of the $s-$surfaces referenced to a fix z-coordinate: 
    59 \begin{equation} 
    60   \label{apdx:A_w_in_s} 
    61   w_s   =  \left.   \frac{\partial z }{\partial t}   \right|_s 
    62   = \frac{\partial z}{\partial s} \; \frac{\partial s}{\partial t} 
    63   = e_3 \, \frac{\partial s}{\partial t} 
    64 \end{equation} 
    65123 
    66124% ================================================================ 
     
    79137    { 
    80138    \begin{array}{*{20}l} 
    81       \nabla \cdot {\rm {\bf U}} 
     139      \nabla \cdot {\mathrm {\mathbf U}} 
    82140      &= \frac{1}{e_1 \,e_2 }  \left[ \left. {\frac{\partial (e_2 \,u)}{\partial i}} \right|_z 
    83141        +\left. {\frac{\partial(e_1 \,v)}{\partial j}} \right|_z  \right] 
     
    115173      $, it becomes:} 
    116174    % 
    117       \nabla \cdot {\rm {\bf U}} 
     175      \nabla \cdot {\mathrm {\mathbf U}} 
    118176      & = \frac{1}{e_1 \,e_2 \,e_3 }  \left[ 
    119177        \left.  \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s 
     
    131189\end{subequations} 
    132190 
    133 Here, $w$ is the vertical velocity relative to the $z-$coordinate system. 
    134 Introducing the dia-surface velocity component, 
    135 $\omega $, defined as the volume flux across the moving $s$-surfaces per unit horizontal area: 
     191Here, $w$ is the vertical velocity relative to the $z-$coordinate system.  
     192Using the first form of (\autoref{apdx:A_s_infin_changes})  
     193and the definitions (\autoref{apdx:A_s_slope}) and (\autoref{apdx:A_w_in_s}) for $\sigma_1$, $\sigma_2$ and  $w_s$, 
     194one can show that the vertical velocity, $w_p$ of a point 
     195moving with the horizontal velocity of the fluid along an $s$ surface is given by  
     196\begin{equation} 
     197\label{apdx:A_w_p} 
     198\begin{split} 
     199w_p  = & \left. \frac{ \partial z }{\partial t} \right|_s 
     200     + \frac{u}{e_1} \left. \frac{ \partial z }{\partial i} \right|_s  
     201     + \frac{v}{e_2} \left. \frac{ \partial z }{\partial j} \right|_s \\ 
     202     = & w_s + u \sigma_1 + v \sigma_2 . 
     203\end{split}      
     204\end{equation} 
     205 The vertical velocity across this surface is denoted by 
    136206\begin{equation} 
    137207  \label{apdx:A_w_s} 
    138   \omega  = w - w_s - \sigma_1 \,u - \sigma_2 \,v    \\ 
    139 \end{equation} 
    140 with $w_s$ given by \autoref{apdx:A_w_in_s}, 
    141 we obtain the expression for the divergence of the velocity in the curvilinear $s-$coordinate system: 
    142 \begin{subequations} 
    143   \begin{align*} 
    144     { 
    145     \begin{array}{*{20}l} 
    146       \nabla \cdot {\rm {\bf U}} 
    147       &= \frac{1}{e_1 \,e_2 \,e_3 }    \left[ 
     208  \omega  = w - w_p = w - ( w_s + \sigma_1 \,u + \sigma_2 \,v )  .  
     209\end{equation} 
     210Hence  
     211\begin{equation} 
     212\frac{1}{e_3 } \frac{\partial}{\partial s}   \left[  w -  u\;\sigma_1  - v\;\sigma_2  \right] =  
     213\frac{1}{e_3 } \frac{\partial}{\partial s} \left[  \omega + w_s \right] =  
     214   \frac{1}{e_3 } \left[ \frac{\partial \omega}{\partial s}  
     215 + \left. \frac{ \partial }{\partial t} \right|_s \frac{\partial z}{\partial s} \right] =  
     216   \frac{1}{e_3 } \frac{\partial \omega}{\partial s} + \frac{1}{e_3 } \left. \frac{ \partial e_3}{\partial t} . \right|_s  
     217\end{equation} 
     218 
     219Using (\autoref{apdx:A_w_s}) in our expression for $\nabla \cdot {\mathrm {\mathbf U}}$ we obtain  
     220our final expression for the divergence of the velocity in the curvilinear $s-$coordinate system: 
     221\begin{equation} 
     222      \nabla \cdot {\mathrm {\mathbf U}} = 
     223         \frac{1}{e_1 \,e_2 \,e_3 }    \left[ 
    148224        \left.  \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s 
    149225        +\left.  \frac{\partial (e_1 \,e_3 \,v)}{\partial j} \right|_s        \right] 
    150226        + \frac{1}{e_3 } \frac{\partial \omega }{\partial s} 
    151         + \frac{1}{e_3 } \frac{\partial w_s       }{\partial s} \\ \\ 
    152       &= \frac{1}{e_1 \,e_2 \,e_3 }    \left[ 
    153         \left.  \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s 
    154         +\left.  \frac{\partial (e_1 \,e_3 \,v)}{\partial j} \right|_s        \right] 
    155         + \frac{1}{e_3 } \frac{\partial \omega }{\partial s} 
    156         + \frac{1}{e_3 } \frac{\partial}{\partial s}  \left(  e_3 \; \frac{\partial s}{\partial t}   \right) \\ \\ 
    157       &= \frac{1}{e_1 \,e_2 \,e_3 }    \left[ 
    158         \left.  \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s 
    159         +\left.  \frac{\partial (e_1 \,e_3 \,v)}{\partial j} \right|_s        \right] 
    160         + \frac{1}{e_3 } \frac{\partial \omega }{\partial s} 
    161         + \frac{\partial}{\partial s} \frac{\partial s}{\partial t} 
    162         + \frac{1}{e_3 } \frac{\partial s}{\partial t} \frac{\partial e_3}{\partial s} \\ \\ 
    163       &= \frac{1}{e_1 \,e_2 \,e_3 }    \left[ 
    164         \left.  \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s 
    165         +\left.  \frac{\partial (e_1 \,e_3 \,v)}{\partial j} \right|_s        \right] 
    166         + \frac{1}{e_3 } \frac{\partial \omega }{\partial s} 
    167         + \frac{1}{e_3 } \frac{\partial e_3}{\partial t} 
    168     \end{array} 
    169         } 
    170   \end{align*} 
    171 \end{subequations} 
     227        + \frac{1}{e_3 } \left. \frac{\partial e_3}{\partial t} \right|_s . 
     228\end{equation} 
    172229 
    173230As a result, the continuity equation \autoref{eq:PE_continuity} in the $s-$coordinates is: 
     
    178235    {\left. {\frac{\partial (e_2 \,e_3 \,u)}{\partial i}} \right|_s 
    179236      +  \left. {\frac{\partial (e_1 \,e_3 \,v)}{\partial j}} \right|_s } \right] 
    180   +\frac{1}{e_3 }\frac{\partial \omega }{\partial s} = 0 
    181 \end{equation} 
    182 A additional term has appeared that take into account 
     237  +\frac{1}{e_3 }\frac{\partial \omega }{\partial s} = 0 . 
     238\end{equation} 
     239An additional term has appeared that takes into account 
    183240the contribution of the time variation of the vertical coordinate to the volume budget. 
    184241 
     
    210267        + w \;\frac{\partial u}{\partial z} \\ \\ 
    211268      &= \left. {\frac{\partial u }{\partial t}} \right|_z 
    212         - \left. \zeta \right|_z v 
    213         +  \frac{1}{e_1 \,e_2 }\left[ { \left.{ \frac{\partial (e_2 \,v)}{\partial i} }\right|_z 
     269        -  \frac{1}{e_1 \,e_2 }\left[ { \left.{ \frac{\partial (e_2 \,v)}{\partial i} }\right|_z 
    214270        -\left.{ \frac{\partial (e_1 \,u)}{\partial j} }\right|_z } \right] \; v 
    215271        +  \frac{1}{2e_1} \left.{ \frac{\partial (u^2+v^2)}{\partial i} } \right|_z 
     
    230286        } \\ \\ 
    231287      &= \left. {\frac{\partial u }{\partial t}} \right|_z 
    232         + \left. \zeta \right|_s \;v 
     288        - \left. \zeta \right|_s \;v 
    233289        + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s \\ 
    234290      &\qquad \qquad \qquad \quad 
    235291        + \frac{w}{e_3 } \;\frac{\partial u}{\partial s} 
    236         - \left[   {\frac{\sigma_1 }{e_3 }\frac{\partial v}{\partial s} 
     292        + \left[   {\frac{\sigma_1 }{e_3 }\frac{\partial v}{\partial s} 
    237293        - \frac{\sigma_2 }{e_3 }\frac{\partial u}{\partial s}}   \right]\;v 
    238294        - \frac{\sigma_1 }{2e_3 }\frac{\partial (u^2+v^2)}{\partial s} \\ \\ 
    239295      &= \left. {\frac{\partial u }{\partial t}} \right|_z 
    240         + \left. \zeta \right|_s \;v 
     296        - \left. \zeta \right|_s \;v 
    241297        + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s \\ 
    242298      &\qquad \qquad \qquad \quad 
     
    245301        - \sigma_1 u\frac{\partial u}{\partial s} - \sigma_1 v\frac{\partial v}{\partial s}} \right] \\ \\ 
    246302      &= \left. {\frac{\partial u }{\partial t}} \right|_z 
    247         + \left. \zeta \right|_s \;v 
     303        - \left. \zeta \right|_s \;v 
    248304        + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s 
    249305        + \frac{1}{e_3} \left[  w - \sigma_2 v - \sigma_1 u  \right] 
    250         \; \frac{\partial u}{\partial s}   \\ 
     306        \; \frac{\partial u}{\partial s} .  \\ 
    251307        % 
    252       \intertext{Introducing $\omega$, the dia-a-surface velocity given by (\autoref{apdx:A_w_s}) } 
     308      \intertext{Introducing $\omega$, the dia-s-surface velocity given by (\autoref{apdx:A_w_s}) } 
    253309      % 
    254310      &= \left. {\frac{\partial u }{\partial t}} \right|_z 
    255         + \left. \zeta \right|_s \;v 
     311        - \left. \zeta \right|_s \;v 
    256312        + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s 
    257         + \frac{1}{e_3 } \left( \omega - w_s \right) \frac{\partial u}{\partial s}   \\ 
     313        + \frac{1}{e_3 } \left( \omega + w_s \right) \frac{\partial u}{\partial s}   \\ 
    258314    \end{array} 
    259315    } 
     
    266322  { 
    267323    \begin{array}{*{20}l} 
    268       w_s  \;\frac{\partial u}{\partial s} 
    269       = \frac{\partial s}{\partial t} \;  \frac{\partial u }{\partial s} 
    270       = \left. {\frac{\partial u }{\partial t}} \right|_s  - \left. {\frac{\partial u }{\partial t}} \right|_z \quad , 
     324      \frac{w_s}{e_3}  \;\frac{\partial u}{\partial s} 
     325      = - \left. \frac{\partial s}{\partial t} \right|_z \;  \frac{\partial u }{\partial s} 
     326      = \left. {\frac{\partial u }{\partial t}} \right|_s  - \left. {\frac{\partial u }{\partial t}} \right|_z \ . 
    271327    \end{array} 
    272328  } 
    273329\] 
    274 leads to the $s-$coordinate formulation of the total $z-$coordinate time derivative, 
     330This leads to the $s-$coordinate formulation of the total $z-$coordinate time derivative, 
    275331\ie the total $s-$coordinate time derivative : 
    276332\begin{align} 
     
    278334  \left. \frac{D u}{D t} \right|_s 
    279335  = \left. {\frac{\partial u }{\partial t}} \right|_s 
    280   + \left. \zeta \right|_s \;v 
     336  - \left. \zeta \right|_s \;v 
    281337  + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s 
    282   + \frac{1}{e_3 } \omega \;\frac{\partial u}{\partial s} 
     338  + \frac{1}{e_3 } \omega \;\frac{\partial u}{\partial s} .  
    283339\end{align} 
    284340Therefore, the vector invariant form of the total time derivative has exactly the same mathematical form in 
     
    306362                                         + \frac{1}{e_3}        \frac{\partial \omega}{\partial s}                       \right] \\ \\ 
    307363                                      &&- \frac{v}{e_1 e_2 }\left(    v \;\frac{\partial e_2 }{\partial i} 
    308                                          -u  \;\frac{\partial e_1 }{\partial j}  \right) \\ 
     364                                         -u  \;\frac{\partial e_1 }{\partial j}  \right) . \\ 
    309365  \end{array} 
    310366  } 
     
    328384       -  e_2 u \;\frac{\partial e_3 }{\partial i} 
    329385       -  e_1 v \;\frac{\partial e_3 }{\partial j}   \right) 
    330        -\frac{1}{e_3}        \frac{\partial \omega}{\partial s}                       \right] \\ \\ 
     386       + \frac{1}{e_3}        \frac{\partial \omega}{\partial s}                       \right] \\ \\ 
    331387    && - \frac{v}{e_1 e_2 }\left(   v  \;\frac{\partial e_2 }{\partial i} 
    332388       -u  \;\frac{\partial e_1 }{\partial j}   \right) \\ \\ 
     
    337393    && - \,u \left[  \frac{1}{e_1 e_2 e_3} \left(   \frac{\partial(e_2 e_3 \, u)}{\partial i} 
    338394       + \frac{\partial(e_1 e_3 \, v)}{\partial j}  \right) 
    339        -\frac{1}{e_3}        \frac{\partial \omega}{\partial s}                       \right] 
     395       + \frac{1}{e_3}        \frac{\partial \omega}{\partial s}                       \right] 
    340396       - \frac{v}{e_1 e_2 }\left(   v   \;\frac{\partial e_2 }{\partial i} 
    341        -u   \;\frac{\partial e_1 }{\partial j}  \right)                  \\ 
     397       -u   \;\frac{\partial e_1 }{\partial j}  \right)     .             \\ 
    342398     % 
    343399    \intertext {Introducing a more compact form for the divergence of the momentum fluxes, 
     
    346402  % 
    347403    &= \left. {\frac{\partial u }{\partial t}} \right|_s 
    348     &+ \left.  \nabla \cdot \left(   {{\rm {\bf U}}\,u}   \right)    \right|_s 
     404    &+ \left.  \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,u}   \right)    \right|_s 
    349405      + \,u \frac{1}{e_3 } \frac{\partial e_3}{\partial t} 
    350406      - \frac{v}{e_1 e_2 }\left(    v  \;\frac{\partial e_2 }{\partial i} 
     
    359415  \label{apdx:A_sco_Dt_flux} 
    360416  \left. \frac{D u}{D t} \right|_s   = \frac{1}{e_3}  \left. \frac{\partial ( e_3\,u)}{\partial t} \right|_s 
    361   + \left.  \nabla \cdot \left(   {{\rm {\bf U}}\,u}   \right)    \right|_s 
     417  + \left.  \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,u}   \right)    \right|_s 
    362418  - \frac{v}{e_1 e_2 }\left(    v  \;\frac{\partial e_2 }{\partial i} 
    363     -u  \;\frac{\partial e_1 }{\partial j}            \right) 
     419    -u  \;\frac{\partial e_1 }{\partial j}            \right). 
    364420\end{flalign} 
    365421which is the total time derivative expressed in the curvilinear $s-$coordinate system. 
     
    377433    & =-\frac{1}{\rho_o e_1 }\left[ {\left. {\frac{\partial p}{\partial i}} \right|_s -\frac{e_1 }{e_3 }\sigma_1 \frac{\partial p}{\partial s}} \right] \\ 
    378434    & =-\frac{1}{\rho_o \,e_1 }\left. {\frac{\partial p}{\partial i}} \right|_s +\frac{\sigma_1 }{\rho_o \,e_3 }\left( {-g\;\rho \;e_3 } \right) \\ 
    379     &=-\frac{1}{\rho_o \,e_1 }\left. {\frac{\partial p}{\partial i}} \right|_s -\frac{g\;\rho }{\rho_o }\sigma_1 
     435    &=-\frac{1}{\rho_o \,e_1 }\left. {\frac{\partial p}{\partial i}} \right|_s -\frac{g\;\rho }{\rho_o }\sigma_1 . 
    380436  \end{split} 
    381437\] 
    382438Applying similar manipulation to the second component and 
    383 replacing $\sigma_1$ and $\sigma_2$ by their expression \autoref{apdx:A_s_slope}, it comes: 
     439replacing $\sigma_1$ and $\sigma_2$ by their expression \autoref{apdx:A_s_slope}, it becomes: 
    384440\begin{equation} 
    385441  \label{apdx:A_grad_p_1} 
     
    391447    -\frac{1}{\rho_o \, e_2 }\left. {\frac{\partial p}{\partial j}} \right|_z 
    392448    &=-\frac{1}{\rho_o \,e_2 } \left(    \left.               {\frac{\partial p}{\partial j}} \right|_s 
    393       + g\;\rho \;\left. {\frac{\partial z}{\partial j}} \right|_s   \right) \\ 
     449      + g\;\rho \;\left. {\frac{\partial z}{\partial j}} \right|_s   \right) . \\ 
    394450  \end{split} 
    395451\end{equation} 
     
    399455 
    400456As in $z$-coordinate, 
    401 the horizontal pressure gradient can be split in two parts following \citet{Marsaleix_al_OM08}. 
     457the horizontal pressure gradient can be split in two parts following \citet{marsaleix.auclair.ea_OM08}. 
    402458Let defined a density anomaly, $d$, by $d=(\rho - \rho_o)/ \rho_o$, 
    403459and a hydrostatic pressure anomaly, $p_h'$, by $p_h'= g \; \int_z^\eta d \; e_3 \; dk$. 
     
    405461\[ 
    406462  \begin{split} 
    407     p &= g\; \int_z^\eta \rho \; e_3 \; dk = g\; \int_z^\eta \left(  \rho_o \, d + 1 \right) \; e_3 \; dk   \\ 
    408     &= g \, \rho_o \; \int_z^\eta d \; e_3 \; dk + g \, \int_z^\eta e_3 \; dk 
     463    p &= g\; \int_z^\eta \rho \; e_3 \; dk = g\; \int_z^\eta \rho_o \left( d + 1 \right) \; e_3 \; dk   \\ 
     464    &= g \, \rho_o \; \int_z^\eta d \; e_3 \; dk + \rho_o g \, \int_z^\eta e_3 \; dk . 
    409465  \end{split} 
    410466\] 
     
    412468\begin{equation} 
    413469  \label{apdx:A_pressure} 
    414   p = \rho_o \; p_h' + g \, ( z + \eta ) 
     470  p = \rho_o \; p_h' + \rho_o \, g \, ( \eta - z ) 
    415471\end{equation} 
    416472and the hydrostatic pressure balance expressed in terms of $p_h'$ and $d$ is: 
    417473\[ 
    418   \frac{\partial p_h'}{\partial k} = - d \, g \, e_3 
     474  \frac{\partial p_h'}{\partial k} = - d \, g \, e_3 . 
    419475\] 
    420476 
    421477Substituing \autoref{apdx:A_pressure} in \autoref{apdx:A_grad_p_1} and 
    422 using the definition of the density anomaly it comes the expression in two parts: 
     478using the definition of the density anomaly it becomes an expression in two parts: 
    423479\begin{equation} 
    424480  \label{apdx:A_grad_p_2} 
     
    426482    -\frac{1}{\rho_o \, e_1 } \left. {\frac{\partial p}{\partial i}} \right|_z 
    427483    &=-\frac{1}{e_1 } \left(     \left.              {\frac{\partial p_h'}{\partial i}} \right|_s 
    428       + g\; d  \;\left. {\frac{\partial z}{\partial i}} \right|_s    \right)  - \frac{g}{e_1 } \frac{\partial \eta}{\partial i} \\ 
     484      + g\; d  \;\left. {\frac{\partial z}{\partial i}} \right|_s    \right)  - \frac{g}{e_1 } \frac{\partial \eta}{\partial i} \\ 
    429485             % 
    430486    -\frac{1}{\rho_o \, e_2 }\left. {\frac{\partial p}{\partial j}} \right|_z 
    431487    &=-\frac{1}{e_2 } \left(    \left.               {\frac{\partial p_h'}{\partial j}} \right|_s 
    432       + g\; d \;\left. {\frac{\partial z}{\partial j}} \right|_s   \right)  - \frac{g}{e_2 } \frac{\partial \eta}{\partial j}\\ 
     488      + g\; d \;\left. {\frac{\partial z}{\partial j}} \right|_s   \right)  - \frac{g}{e_2 } \frac{\partial \eta}{\partial j} . \\ 
    433489  \end{split} 
    434490\end{equation} 
     
    463519    -   \frac{1}{e_1 } \left(    \frac{\partial p_h'}{\partial i} + g\; d  \; \frac{\partial z}{\partial i}    \right) 
    464520    -   \frac{g}{e_1 } \frac{\partial \eta}{\partial i} 
    465     +   D_u^{\vect{U}}  +   F_u^{\vect{U}} 
     521    +   D_u^{\vect{U}}  +   F_u^{\vect{U}} , 
    466522  \end{multline} 
    467523  \begin{multline} 
     
    473529    -   \frac{1}{e_2 } \left(    \frac{\partial p_h'}{\partial j} + g\; d  \; \frac{\partial z}{\partial j}    \right) 
    474530    -   \frac{g}{e_2 } \frac{\partial \eta}{\partial j} 
    475     +  D_v^{\vect{U}}  +   F_v^{\vect{U}} 
     531    +  D_v^{\vect{U}}  +   F_v^{\vect{U}} . 
    476532  \end{multline} 
    477533\end{subequations} 
     
    483539    \label{apdx:A_PE_dyn_flux_u} 
    484540    \frac{1}{e_3} \frac{\partial \left(  e_3\,u  \right) }{\partial t} = 
    485     \nabla \cdot \left(   {{\rm {\bf U}}\,u}   \right) 
     541    - \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,u}   \right) 
    486542    +   \left\{ {f + \frac{1}{e_1 e_2 }\left(    v  \;\frac{\partial e_2 }{\partial i} 
    487543          -u  \;\frac{\partial e_1 }{\partial j}            \right)} \right\} \,v     \\ 
    488544    -   \frac{1}{e_1 } \left(    \frac{\partial p_h'}{\partial i} + g\; d  \; \frac{\partial z}{\partial i}    \right) 
    489545    -   \frac{g}{e_1 } \frac{\partial \eta}{\partial i} 
    490     +   D_u^{\vect{U}}  +   F_u^{\vect{U}} 
     546    +   D_u^{\vect{U}}  +   F_u^{\vect{U}} , 
    491547  \end{multline} 
    492548  \begin{multline} 
    493549    \label{apdx:A_dyn_flux_v} 
    494550    \frac{1}{e_3}\frac{\partial \left(  e_3\,v  \right) }{\partial t}= 
    495     -  \nabla \cdot \left(   {{\rm {\bf U}}\,v}   \right) 
    496     +   \left\{ {f + \frac{1}{e_1 e_2 }\left(    v  \;\frac{\partial e_2 }{\partial i} 
     551    -  \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,v}   \right) 
     552    -   \left\{ {f + \frac{1}{e_1 e_2 }\left(    v  \;\frac{\partial e_2 }{\partial i} 
    497553          -u  \;\frac{\partial e_1 }{\partial j}            \right)} \right\} \,u     \\ 
    498554    -   \frac{1}{e_2 } \left(    \frac{\partial p_h'}{\partial j} + g\; d  \; \frac{\partial z}{\partial j}    \right) 
    499555    -   \frac{g}{e_2 } \frac{\partial \eta}{\partial j} 
    500     +  D_v^{\vect{U}}  +   F_v^{\vect{U}} 
     556    +  D_v^{\vect{U}}  +   F_v^{\vect{U}} .  
    501557  \end{multline} 
    502558\end{subequations} 
     
    505561\begin{equation} 
    506562  \label{apdx:A_dyn_zph} 
    507   \frac{\partial p_h'}{\partial k} = - d \, g \, e_3 
     563  \frac{\partial p_h'}{\partial k} = - d \, g \, e_3 . 
    508564\end{equation} 
    509565 
     
    531587  \left[           \frac{\partial }{\partial i} \left( {e_2 \,e_3 \;Tu} \right) 
    532588    +   \frac{\partial }{\partial j} \left( {e_1 \,e_3 \;Tv} \right)               \right]       \\ 
    533   +  \frac{1}{e_3}  \frac{\partial }{\partial k} \left(                   Tw  \right) 
     589  -  \frac{1}{e_3}  \frac{\partial }{\partial k} \left(                   Tw  \right) 
    534590  +  D^{T} +F^{T} 
    535591\end{multline} 
    536592 
    537 The expression for the advection term is a straight consequence of (A.4), 
     593The expression for the advection term is a straight consequence of (\autoref{apdx:A_sco_Continuity}), 
    538594the expression of the 3D divergence in the $s-$coordinates established above.  
    539595 
  • NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/annex_B.tex

    r10442 r11422  
    5353  { 
    5454  \begin{array}{*{20}l} 
    55     D^T=& \frac{1}{e_1\,e_2\,e_3 }\;\left[ {\ \ \ \ e_2\,e_3\,A^{lT} \;\left. 
    56           {\frac{\partial }{\partial i}\left( {\frac{1}{e_1}\;\left. {\frac{\partial T}{\partial i}} \right|_s -\frac{\sigma_1 }{e_3 }\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right.  \\ 
    57         &\qquad \quad \ \ \ +e_1\,e_3\,A^{lT} \;\left. {\frac{\partial }{\partial j}\left( {\frac{1}{e_2 }\;\left. {\frac{\partial T}{\partial j}} \right|_s -\frac{\sigma_2 }{e_3 }\;\frac{\partial T}{\partial s}} \right)} \right|_s \\ 
    58         &\qquad \quad \ \ \ + e_1\,e_2\,A^{lT} \;\frac{\partial }{\partial s}\left( {-\frac{\sigma_1 }{e_1 }\;\left. {\frac{\partial T}{\partial i}} \right|_s -\frac{\sigma_2 }{e_2 }\;\left. {\frac{\partial T}{\partial j}} \right|_s } \right. 
    59           \left. {\left. {+\left( {\varepsilon +\sigma_1^2+\sigma_2 ^2} \right)\;\frac{1}{e_3 }\;\frac{\partial T}{\partial s}} \right)\;\;} \right] 
     55    D^T= \frac{1}{e_1\,e_2\,e_3 } & \left\{ \quad \quad \frac{\partial }{\partial i}  \left. \left[  e_2\,e_3 \, A^{lT} 
     56                               \left( \  \frac{1}{e_1}\; \left. \frac{\partial T}{\partial i} \right|_s  
     57                                       -\frac{\sigma_1 }{e_3 } \; \frac{\partial T}{\partial s} \right) \right]  \right|_s  \right. \\ 
     58        &  \quad \  +   \            \left.   \frac{\partial }{\partial j}  \left. \left[  e_1\,e_3 \, A^{lT} 
     59                               \left( \ \frac{1}{e_2 }\; \left. \frac{\partial T}{\partial j} \right|_s  
     60                                       -\frac{\sigma_2 }{e_3 } \; \frac{\partial T}{\partial s} \right) \right]  \right|_s  \right. \\ 
     61        &  \quad \  +   \           \left.  e_1\,e_2\, \frac{\partial }{\partial s}  \left[ A^{lT} \; \left(  
     62                     -\frac{\sigma_1 }{e_1 } \; \left. \frac{\partial T}{\partial i} \right|_s  
     63                     -\frac{\sigma_2 }{e_2 } \; \left. \frac{\partial T}{\partial j} \right|_s  
     64                          +\left( \varepsilon +\sigma_1^2+\sigma_2 ^2 \right) \; \frac{1}{e_3 } \; \frac{\partial T}{\partial s} \right) \; \right] \;  \right\} . 
    6065  \end{array} 
    6166          } 
    6267\end{align*} 
    6368 
    64 Equation \autoref{apdx:B2} is obtained from \autoref{apdx:B1} without any additional assumption. 
     69\autoref{apdx:B2} is obtained from \autoref{apdx:B1} without any additional assumption. 
    6570Indeed, for the special case $k=z$ and thus $e_3 =1$, 
    6671we introduce an arbitrary vertical coordinate $s = s (i,j,z)$ as in \autoref{apdx:A} and 
     
    9095  { 
    9196  \begin{array}{*{20}l} 
    92     \intertext{Noting that $\frac{1}{e_1} \left. \frac{\partial e_3 }{\partial i} \right|_s = \frac{\partial \sigma_1 }{\partial s}$, it becomes:} 
     97    \intertext{Noting that $\frac{1}{e_1} \left. \frac{\partial e_3 }{\partial i} \right|_s = \frac{\partial \sigma_1 }{\partial s}$, this becomes:} 
    9398    % 
    94     & =\frac{1}{e_1\,e_2\,e_3 }\left[ {\left. {\;\;\;\frac{\partial }{\partial i}\left( {\frac{e_2\,e_3 }{e_1}\,A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)} \right|_s \left. -\, {e_3 \frac{\partial }{\partial i}\left( {\frac{e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\ 
     99    D^T & =\frac{1}{e_1\,e_2\,e_3 }\left[ {\left. {\;\;\;\frac{\partial }{\partial i}\left( {\frac{e_2\,e_3 }{e_1}\,A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)} \right|_s \left. -\, {e_3 \frac{\partial }{\partial i}\left( {\frac{e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\ 
    95100    & \qquad \qquad \quad -e_2 A^{lT}\;\frac{\partial \sigma_1 }{\partial s}\left. {\frac{\partial T}{\partial i}} \right|_s -e_1 \,\sigma_1 \frac{\partial }{\partial s}\left( {\frac{e_2 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right) \\ 
    96     & \qquad \qquad \quad\shoveright{ \left. { +e_1 \,\sigma_1 \frac{\partial }{\partial s}\left( {\frac{e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial z}} \right)\;\;\;} \right] }\\ 
     101    & \qquad \qquad \quad\shoveright{ \left. { +e_1 \,\sigma_1 \frac{\partial }{\partial s}\left( {\frac{e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right] }\\ 
    97102    \\ 
    98103    &=\frac{1}{e_1 \,e_2 \,e_3 } \left[ {\left. {\;\;\;\frac{\partial }{\partial i} \left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)} \right|_s \left. {-\frac{\partial }{\partial i}\left( {e_2 \,\sigma_1 A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\ 
    99104    & \qquad \qquad \quad \left. {+\frac{e_2 \,\sigma_1 }{e_3}A^{lT}\;\frac{\partial T}{\partial s} \;\frac{\partial e_3 }{\partial i}}  \right|_s -e_2 A^{lT}\;\frac{\partial \sigma_1 }{\partial s}\left. {\frac{\partial T}{\partial i}} \right|_s \\ 
    100105    & \qquad \qquad \quad-e_2 \,\sigma_1 \frac{\partial}{\partial s}\left( {A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 \,\sigma_1 ^2}{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right) \\ 
    101     & \qquad \qquad \quad\shoveright{ \left. {-\frac{\partial \left( {e_1 \,e_2 \,\sigma_1 } \right)}{\partial s} \left( {\frac{\sigma_1 }{e_3}A^{lT}\;\frac{\partial T}{\partial s}} \right) + \frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right]} 
     106    & \qquad \qquad \quad\shoveright{ \left. {-\frac{\partial \left( {e_1 \,e_2 \,\sigma_1 } \right)}{\partial s} \left( {\frac{\sigma_1 }{e_3}A^{lT}\;\frac{\partial T}{\partial s}} \right) + \frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right]} . 
    102107  \end{array} 
    103108      } \\ 
     
    105110  \begin{array}{*{20}l} 
    106111    % 
    107     \intertext{using the same remark as just above, it becomes:} 
     112    \intertext{Using the same remark as just above, $D^T$ becomes:} 
    108113    % 
    109     &= \frac{1}{e_1 \,e_2 \,e_3 } \left[ {\left. {\;\;\;\frac{\partial }{\partial i} \left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s -e_2 \,\sigma_1 A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right.\;\;\; \\ 
     114   D^T &= \frac{1}{e_1 \,e_2 \,e_3 } \left[ {\left. {\;\;\;\frac{\partial }{\partial i} \left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s -e_2 \,\sigma_1 A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right.\;\;\; \\ 
    110115    & \qquad \qquad \quad+\frac{e_1 \,e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}\;\frac{\partial \sigma_1 }{\partial s} - \frac {\sigma_1 }{e_3} A^{lT} \;\frac{\partial \left( {e_1 \,e_2 \,\sigma_1 } \right)}{\partial s}\;\frac{\partial T}{\partial s} \\ 
    111116    & \qquad \qquad \quad-e_2 \left( {A^{lT}\;\frac{\partial \sigma_1 }{\partial s}\left. {\frac{\partial T}{\partial i}} \right|_s +\frac{\partial }{\partial s}\left( {\sigma_1 A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)-\frac{\partial \sigma_1 }{\partial s}\;A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right) \\ 
    112     & \qquad \qquad \quad\shoveright{\left. {+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 \,\sigma_1 ^2}{e_3 }A^{lT}\;\frac{\partial T}{\partial s}+\frac{e_1 \,e_2}{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right] } 
     117    & \qquad \qquad \quad\shoveright{\left. {+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 \,\sigma_1 ^2}{e_3 }A^{lT}\;\frac{\partial T}{\partial s}+\frac{e_1 \,e_2}{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right] . } 
    113118  \end{array} 
    114119      } \\ 
     
    117122    % 
    118123    \intertext{Since the horizontal scale factors do not depend on the vertical coordinate, 
    119     the last term of the first line and the first term of the last line cancel, while 
    120     the second line reduces to a single vertical derivative, so it becomes:} 
     124    the two terms on the second line cancel, while 
     125    the third line reduces to a single vertical derivative, so it becomes:} 
    121126  % 
    122     & =\frac{1}{e_1 \,e_2 \,e_3 }\left[ {\left. {\;\;\;\frac{\partial }{\partial i}\left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s -e_2 \,\sigma_1 \,A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\ 
     127    D^T & =\frac{1}{e_1 \,e_2 \,e_3 }\left[ {\left. {\;\;\;\frac{\partial }{\partial i}\left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s -e_2 \,\sigma_1 \,A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\ 
    123128    & \qquad \qquad \quad \shoveright{ \left. {+\frac{\partial }{\partial s}\left( {-e_2 \,\sigma_1 \,A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s +A^{lT}\frac{e_1 \,e_2 }{e_3 }\;\left( {\varepsilon +\sigma_1 ^2} \right)\frac{\partial T}{\partial s}} \right)\;\;\;} \right]} \\ 
    124129    % 
    125     \intertext{in other words, the horizontal/vertical Laplacian operator in the ($i$,$s$) plane takes the following form:} 
     130    \intertext{In other words, the horizontal/vertical Laplacian operator in the ($i$,$s$) plane takes the following form:} 
    126131  \end{array} 
    127132  } \\ 
     
    162167the ($i$,$j$,$k$) curvilinear coordinate system in which 
    163168the equations of the ocean circulation model are formulated, 
    164 takes the following form \citep{Redi_JPO82}: 
     169takes the following form \citep{redi_JPO82}: 
    165170 
    166171\begin{equation} 
     
    169174  \left[ {{ 
    170175        \begin{array}{*{20}c} 
    171           {1+a_1 ^2} \hfill & {-a_1 a_2 } \hfill & {-a_1 } \hfill \\ 
    172           {-a_1 a_2 } \hfill & {1+a_2 ^2} \hfill & {-a_2 } \hfill \\ 
    173           {-a_1 } \hfill & {-a_2 } \hfill & {\varepsilon +a_1 ^2+a_2 ^2} \hfill \\ 
     176          {1+a_2 ^2 +\varepsilon a_1 ^2} \hfill & {-a_1 a_2 (1-\varepsilon)} \hfill & {-a_1 (1-\varepsilon) } \hfill \\ 
     177          {-a_1 a_2 (1-\varepsilon) } \hfill & {1+a_1 ^2 +\varepsilon a_2 ^2} \hfill & {-a_2 (1-\varepsilon)} \hfill \\ 
     178          {-a_1 (1-\varepsilon)} \hfill & {-a_2 (1-\varepsilon)} \hfill & {\varepsilon +a_1 ^2+a_2 ^2} \hfill \\ 
    174179        \end{array} 
    175180      }} \right] 
    176181\end{equation} 
    177 where ($a_1$, $a_2$) are the isopycnal slopes in ($\textbf{i}$, $\textbf{j}$) directions, relative to geopotentials: 
     182where ($a_1$, $a_2$) are $(-1) \times$ the isopycnal slopes in 
     183($\textbf{i}$, $\textbf{j}$) directions, relative to geopotentials (or 
     184equivalently the slopes of the geopotential surfaces in the isopycnal 
     185coordinate framework): 
    178186\[ 
    179187  a_1 =\frac{e_3 }{e_1 }\left( {\frac{\partial \rho }{\partial i}} \right)\left( {\frac{\partial \rho }{\partial k}} \right)^{-1} 
     
    182190  \right)\left( {\frac{\partial \rho }{\partial k}} \right)^{-1} 
    183191\] 
    184  
    185 In practice, isopycnal slopes are generally less than $10^{-2}$ in the ocean, 
    186 so $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{Cox1987}: 
     192and, as before, $\epsilon = A^{vT} / A^{lT}$. 
     193 
     194In practice, $\epsilon$ is small and isopycnal slopes are generally less than $10^{-2}$ in the ocean, 
     195so $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{cox_OM87}. Keeping leading order terms\footnote{Apart from the (1,0)  
     196and (0,1) elements which are set to zero. See \citet{griffies_bk04}, section 14.1.4.1 for a discussion of this point.}: 
    187197\begin{subequations} 
    188198  \label{apdx:B4} 
     
    226236The isopycnal diffusion operator \autoref{apdx:B4}, 
    227237\autoref{apdx:B_ldfiso} conserves tracer quantity and dissipates its square. 
    228 The demonstration of the first property is trivial as \autoref{apdx:B4} is the divergence of fluxes. 
    229 Let us demonstrate the second one: 
     238As \autoref{apdx:B4} is the divergence of a flux, the demonstration of the first property is trivial, providing that the flux normal to the boundary is zero  
     239(as it is when $A_h$ is zero at the boundary). Let us demonstrate the second one: 
    230240\[ 
    231241  \iiint\limits_D T\;\nabla .\left( {\textbf{A}}_{\textbf{I}} \nabla T \right)\,dv 
     
    236246  { 
    237247  \begin{array}{*{20}l} 
    238     \nabla T\;.\left( {{\rm {\bf A}}_{\rm {\bf I}} \nabla T} 
     248    \nabla T\;.\left( {{\mathrm {\mathbf A}}_{\mathrm {\mathbf I}} \nabla T} 
    239249    \right)&=A^{lT}\left[ {\left( {\frac{\partial T}{\partial i}} \right)^2-2a_1 
    240250             \frac{\partial T}{\partial i}\frac{\partial T}{\partial k}+\left( 
     
    246256             j}-a_2 \frac{\partial T}{\partial k}} \right)^2} 
    247257             +\varepsilon \left(\frac{\partial T}{\partial k}\right) ^2\right]      \\ 
    248            & \geq 0 
     258           & \geq 0 .  
    249259  \end{array} 
    250260             } 
     
    275285\end{equation} 
    276286 
    277 where ($r_1$, $r_2$) are the isopycnal slopes in ($\textbf{i}$, $\textbf{j}$) directions, 
    278 relative to $s$-coordinate surfaces: 
     287where ($r_1$, $r_2$) are $(-1)\times$ the isopycnal slopes in ($\textbf{i}$, $\textbf{j}$) directions, 
     288relative to $s$-coordinate surfaces (or equivalently the slopes of the 
     289$s$-coordinate surfaces in the isopycnal coordinate framework): 
    279290\[ 
    280291  r_1 =\frac{e_3 }{e_1 }\left( {\frac{\partial \rho }{\partial i}} \right)\left( {\frac{\partial \rho }{\partial s}} \right)^{-1} 
     
    284295\] 
    285296 
    286 To prove \autoref{apdx:B5} by direct re-expression of \autoref{apdx:B_ldfiso} is straightforward, but laborious. 
     297To prove \autoref{apdx:B_ldfiso_s} by direct re-expression of \autoref{apdx:B_ldfiso} is straightforward, but laborious. 
    287298An easier way is first to note (by reversing the derivation of \autoref{sec:B_2} from \autoref{sec:B_1} ) that 
    288299the weak-slope operator may be \emph{exactly} reexpressed in non-orthogonal $i,j,\rho$-coordinates as 
     
    306317the (rotated,orthogonal) isoneutral axes to the non-orthogonal $i,j,\rho$-coordinates. 
    307318The further transformation into $i,j,s$-coordinates is exact, whatever the steepness of the $s$-surfaces, 
    308 in the same way as the transformation of horizontal/vertical Laplacian diffusion in $z$-coordinates, 
     319in the same way as the transformation of horizontal/vertical Laplacian diffusion in $z$-coordinates in 
    309320\autoref{sec:B_1} onto $s$-coordinates is exact, however steep the $s$-surfaces. 
    310321 
     
    316327\label{sec:B_3} 
    317328 
    318 The second order momentum diffusion operator (Laplacian) in the $z$-coordinate is found by 
     329The second order momentum diffusion operator (Laplacian) in $z$-coordinates is found by 
    319330applying \autoref{eq:PE_lap_vector}, the expression for the Laplacian of a vector, 
    320331to the horizontal velocity vector: 
     
    361372\end{align*} 
    362373Using \autoref{eq:PE_div}, the definition of the horizontal divergence, 
    363 the third componant of the second vector is obviously zero and thus : 
     374the third component of the second vector is obviously zero and thus : 
    364375\[ 
    365   \Delta {\textbf{U}}_h = \nabla _h \left( \chi \right) - \nabla _h \times \left( \zeta \right) + \frac {1}{e_3 } \frac {\partial }{\partial k} \left( {\frac {1}{e_3 } \frac{\partial {\textbf{ U}}_h }{\partial k}} \right) 
     376  \Delta {\textbf{U}}_h = \nabla _h \left( \chi \right) - \nabla _h \times \left( \zeta \textbf{k} \right) + \frac {1}{e_3 } \frac {\partial }{\partial k} \left( {\frac {1}{e_3 } \frac{\partial {\textbf{ U}}_h }{\partial k}} \right) .  
    366377\] 
    367378 
     
    379390  - \nabla _h \times \left( {A^{lm}\;\zeta \;{\textbf{k}}} \right) 
    380391  + \frac{1}{e_3 }\frac{\partial }{\partial k}\left( {\frac{A^{vm}\;}{e_3 } 
    381       \frac{\partial {\rm {\bf U}}_h }{\partial k}} \right) \\ 
     392      \frac{\partial {\mathrm {\mathbf U}}_h }{\partial k}} \right) , \\ 
    382393\end{equation} 
    383394that is, in expanded form: 
     
    386397  & = \frac{1}{e_1} \frac{\partial \left( {A^{lm}\chi   } \right)}{\partial i} 
    387398    -\frac{1}{e_2} \frac{\partial \left( {A^{lm}\zeta } \right)}{\partial j} 
    388     +\frac{1}{e_3} \frac{\partial u}{\partial k}      \\ 
     399    +\frac{1}{e_3} \frac{\partial }{\partial k} \left( \frac{A^{vm}}{e_3} \frac{\partial u}{\partial k} \right)   ,   \\ 
    389400  D^{\textbf{U}}_v 
    390401  & = \frac{1}{e_2 }\frac{\partial \left( {A^{lm}\chi   } \right)}{\partial j} 
    391402    +\frac{1}{e_1 }\frac{\partial \left( {A^{lm}\zeta } \right)}{\partial i} 
    392     +\frac{1}{e_3} \frac{\partial v}{\partial k} 
     403    +\frac{1}{e_3} \frac{\partial }{\partial k} \left( \frac{A^{vm}}{e_3} \frac{\partial v}{\partial k} \right) . 
    393404\end{align*} 
    394405 
  • NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/annex_E.tex

    r10442 r11422  
    4949 
    5050This results in a dissipatively dominant (\ie hyper-diffusive) truncation error 
    51 \citep{Shchepetkin_McWilliams_OM05}. 
    52 The overall performance of the advection scheme is similar to that reported in \cite{Farrow1995}. 
     51\citep{shchepetkin.mcwilliams_OM05}. 
     52The overall performance of the advection scheme is similar to that reported in \cite{farrow.stevens_JPO95}. 
    5353It is a relatively good compromise between accuracy and smoothness. 
    5454It is not a \emph{positive} scheme meaning false extrema are permitted but 
     
    6565the second term which is the diffusive part of the scheme, is evaluated using the \textit{before} velocity 
    6666(forward in time). 
    67 This is discussed by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme. 
     67This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme. 
    6868UBS and QUICK schemes only differ by one coefficient. 
    69 Substituting 1/6 with 1/8 in (\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. 
     69Substituting 1/6 with 1/8 in (\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}. 
    7070This option is not available through a namelist parameter, since the 1/6 coefficient is hard coded. 
    7171Nevertheless it is quite easy to make the substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme. 
     
    8080$\tau_w^{ubs}$ will be evaluated using either \textit{(a)} a centered $2^{nd}$ order scheme, 
    8181or \textit{(b)} a TVD scheme, or \textit{(c)} an interpolation based on conservative parabolic splines following 
    82 \citet{Shchepetkin_McWilliams_OM05} implementation of UBS in ROMS, or \textit{(d)} an UBS. 
     82\citet{shchepetkin.mcwilliams_OM05} implementation of UBS in ROMS, or \textit{(d)} an UBS. 
    8383The $3^{rd}$ case has dispersion properties similar to an eight-order accurate conventional scheme. 
    8484 
     
    255255\subsection{Griffies iso-neutral diffusion operator} 
    256256 
    257 Let try to define a scheme that get its inspiration from the \citet{Griffies_al_JPO98} scheme, 
     257Let try to define a scheme that get its inspiration from the \citet{griffies.gnanadesikan.ea_JPO98} scheme, 
    258258but is formulated within the \NEMO framework 
    259259(\ie using scale factors rather than grid-size and having a position of $T$-points that 
     
    272272Nevertheless, this technique works fine for $T$ and $S$ as they are active tracers 
    273273(\ie they enter the computation of density), but it does not work for a passive tracer. 
    274 \citep{Griffies_al_JPO98} introduce a different way to discretise the off-diagonal terms that 
     274\citep{griffies.gnanadesikan.ea_JPO98} introduce a different way to discretise the off-diagonal terms that 
    275275nicely solve the problem. 
    276276The idea is to get rid of combinations of an averaged in one direction combined with 
     
    308308\begin{figure}[!ht] 
    309309  \begin{center} 
    310     \includegraphics[width=0.70\textwidth]{Fig_ISO_triad} 
     310    \includegraphics[width=\textwidth]{Fig_ISO_triad} 
    311311    \caption{ 
    312312      \protect\label{fig:ISO_triad} 
     
    508508\] 
    509509 
    510 \citep{Griffies_JPO98} introduces another way to implement the eddy induced advection, the so-called skew form. 
     510\citep{griffies_JPO98} introduces another way to implement the eddy induced advection, the so-called skew form. 
    511511It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity. 
    512512For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be transformed as follows: 
     
    574574The horizontal component reduces to the one use for an horizontal laplacian operator and 
    575575the vertical one keeps the same complexity, but not more. 
    576 This property has been used to reduce the computational time \citep{Griffies_JPO98}, 
     576This property has been used to reduce the computational time \citep{griffies_JPO98}, 
    577577but it is not of practical use as usually $A \neq A_e$. 
    578578Nevertheless this property can be used to choose a discret form of \autoref{eq:eiv_skew_continuous} which 
  • NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/annex_iso.tex

    r10442 r11422  
    44\newcommand{\rML}[1][i]{\ensuremath{_{\mathrm{ML}\,#1}}} 
    55\newcommand{\rMLt}[1][i]{\tilde{r}_{\mathrm{ML}\,#1}} 
    6 \newcommand{\triad}[6][]{\ensuremath{{}_{#2}^{#3}{\mathbb{#4}_{#1}}_{#5}^{\,#6}}} 
     6%% Move to ../../global/new_cmds.tex to avoid error with \listoffigures 
     7%\newcommand{\triad}[6][]{\ensuremath{{}_{#2}^{#3}{\mathbb{#4}_{#1}}_{#5}^{\,#6}} 
    78\newcommand{\triadd}[5]{\ensuremath{{}_{#1}^{#2}{\mathbb{#3}}_{#4}^{\,#5}}} 
    89\newcommand{\triadt}[5]{\ensuremath{{}_{#1}^{#2}{\tilde{\mathbb{#3}}}_{#4}^{\,#5}}} 
     
    5253  the vertical skew flux is further reduced to ensure no vertical buoyancy flux, 
    5354  giving an almost pure horizontal diffusive tracer flux within the mixed layer. 
    54   This is similar to the tapering suggested by \citet{Gerdes1991}. See \autoref{subsec:Gerdes-taper} 
     55  This is similar to the tapering suggested by \citet{gerdes.koberle.ea_CD91}. See \autoref{subsec:Gerdes-taper} 
    5556\item[\np{ln\_botmix\_triad}] 
    5657  See \autoref{sec:iso_bdry}.  
     
    7172\label{sec:iso} 
    7273 
    73 We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98}, 
     74We have implemented into \NEMO a scheme inspired by \citet{griffies.gnanadesikan.ea_JPO98}, 
    7475but formulated within the \NEMO framework, using scale factors rather than grid-sizes. 
    7576 
     
    194195\subsection{Expression of the skew-flux in terms of triad slopes} 
    195196 
    196 \citep{Griffies_al_JPO98} introduce a different discretization of the off-diagonal terms that 
     197\citep{griffies.gnanadesikan.ea_JPO98} introduce a different discretization of the off-diagonal terms that 
    197198nicely solves the problem. 
    198199% Instead of multiplying the mean slope calculated at the $u$-point by 
     
    201202\begin{figure}[tb] 
    202203  \begin{center} 
    203     \includegraphics[width=1.05\textwidth]{Fig_GRIFF_triad_fluxes} 
     204    \includegraphics[width=\textwidth]{Fig_GRIFF_triad_fluxes} 
    204205    \caption{ 
    205206      \protect\label{fig:ISO_triad} 
     
    265266\begin{figure}[tb] 
    266267  \begin{center} 
    267     \includegraphics[width=0.80\textwidth]{Fig_GRIFF_qcells} 
     268    \includegraphics[width=\textwidth]{Fig_GRIFF_qcells} 
    268269    \caption{ 
    269270      \protect\label{fig:qcells} 
     
    473474 
    474475To complete the discretization we now need only specify the triad volumes $_i^k\mathbb{V}_{i_p}^{k_p}$. 
    475 \citet{Griffies_al_JPO98} identifies these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter cells, 
     476\citet{griffies.gnanadesikan.ea_JPO98} identifies these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter cells, 
    476477defined in terms of the distances between $T$, $u$,$f$ and $w$-points. 
    477478This is the natural discretization of \autoref{eq:cts-var}. 
     
    658659\begin{figure}[h] 
    659660  \begin{center} 
    660     \includegraphics[width=0.60\textwidth]{Fig_GRIFF_bdry_triads} 
     661    \includegraphics[width=\textwidth]{Fig_GRIFF_bdry_triads} 
    661662    \caption{ 
    662663      \protect\label{fig:bdry_triads} 
     
    685686As discussed in \autoref{subsec:LDF_slp_iso}, 
    686687iso-neutral slopes relative to geopotentials must be bounded everywhere, 
    687 both for consistency with the small-slope approximation and for numerical stability \citep{Cox1987, Griffies_Bk04}. 
     688both for consistency with the small-slope approximation and for numerical stability \citep{cox_OM87, griffies_bk04}. 
    688689The bound chosen in \NEMO is applied to each component of the slope separately and 
    689690has a value of $1/100$ in the ocean interior. 
     
    732733\[ 
    733734  % \label{eq:iso_tensor_ML} 
    734   D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 
     735  D^{lT}=\nabla {\mathrm {\mathbf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 
    735736  \mbox{with}\quad \;\;\Re =\left( {{ 
    736737        \begin{array}{*{20}c} 
     
    829830    (\eg the green triad $i_p=1/2,k_p=-1/2$) are tapered to the appropriate basal triad.} 
    830831  % } 
    831   \includegraphics[width=0.60\textwidth]{Fig_GRIFF_MLB_triads} 
     832  \includegraphics[width=\textwidth]{Fig_GRIFF_MLB_triads} 
    832833\end{figure} 
    833834% >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    847848\[ 
    848849  % \label{eq:iso_tensor_ML2} 
    849   D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 
     850  D^{lT}=\nabla {\mathrm {\mathbf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 
    850851  \mbox{with}\quad \;\;\Re =\left( {{ 
    851852        \begin{array}{*{20}c} 
     
    859860\footnote{ 
    860861  To ensure good behaviour where horizontal density gradients are weak, 
    861   we in fact follow \citet{Gerdes1991} and 
     862  we in fact follow \citet{gerdes.koberle.ea_CD91} and 
    862863  set $\rML^*=\mathrm{sgn}(\tilde{r}_i)\min(|\rMLt^2/\tilde{r}_i|,|\tilde{r}_i|)-\sigma_i$. 
    863864} 
     
    865866This approach is similar to multiplying the iso-neutral diffusion coefficient by 
    866867$\tilde{r}_{\mathrm{max}}^{-2}\tilde{r}_i^{-2}$ for steep slopes, 
    867 as suggested by \citet{Gerdes1991} (see also \citet{Griffies_Bk04}). 
     868as suggested by \citet{gerdes.koberle.ea_CD91} (see also \citet{griffies_bk04}). 
    868869Again it is applied separately to each triad $_i^k\mathbb{R}_{i_p}^{k_p}$ 
    869870 
     
    925926 
    926927However, when \np{ln\_traldf\_triad} is set true, 
    927 \NEMO instead implements eddy induced advection according to the so-called skew form \citep{Griffies_JPO98}. 
     928\NEMO instead implements eddy induced advection according to the so-called skew form \citep{griffies_JPO98}. 
    928929It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity. 
    929930For example in the (\textbf{i},\textbf{k}) plane, 
     
    11391140it is equivalent to a horizontal eiv (eddy-induced velocity) that is uniform within the mixed layer 
    11401141\autoref{eq:eiv_v}. 
    1141 This ensures that the eiv velocities do not restratify the mixed layer \citep{Treguier1997,Danabasoglu_al_2008}. 
     1142This ensures that the eiv velocities do not restratify the mixed layer \citep{treguier.held.ea_JPO97,danabasoglu.ferrari.ea_JC08}. 
    11421143Equivantly, in terms of the skew-flux formulation we use here, 
    11431144the linear slope tapering within the mixed-layer gives a linearly varying vertical flux, 
     
    11531154$uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$ (integer $i$, integer +1/2 $j$, integer +1/2 $k$) 
    11541155points (see Table \autoref{tab:cell}) respectively. 
    1155 We follow \citep{Griffies_Bk04} and calculate the streamfunction at a given $uw$-point from 
     1156We follow \citep{griffies_bk04} and calculate the streamfunction at a given $uw$-point from 
    11561157the surrounding four triads according to: 
    11571158\[ 
  • NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_ASM.tex

    r10442 r11422  
    88\label{chap:ASM} 
    99 
    10 Authors: D. Lea,  M. Martin, K. Mogensen, A. Weaver, ...   % do we keep 
     10\minitoc 
    1111 
    12 \minitoc 
     12\vfill 
     13\begin{figure}[b] 
     14\subsubsection*{Changes record} 
     15\begin{tabular}{l||l|m{0.65\linewidth}} 
     16    Release   & Author        & Modifications \\ 
     17    {\em 4.0} & {\em D. J. Lea} & {\em \NEMO 4.0 updates}  \\ 
     18    {\em 3.4} & {\em D. J. Lea, M. Martin, K. Mogensen, A. Weaver} & {\em Initial version}  \\ 
     19\end{tabular} 
     20\end{figure} 
    1321 
    1422\newpage 
    1523 
    1624The ASM code adds the functionality to apply increments to the model variables: temperature, salinity, 
    17 sea surface height, velocity and sea ice concentration.  
     25sea surface height, velocity and sea ice concentration. 
    1826These are read into the model from a NetCDF file which may be produced by separate data assimilation code. 
    1927The code can also output model background fields which are used as an input to data assimilation code. 
     
    3745it may be preferable to introduce the increment gradually into the ocean model in order to 
    3846minimize spurious adjustment processes. 
    39 This technique is referred to as Incremental Analysis Updates (IAU) \citep{Bloom_al_MWR96}. 
     47This technique is referred to as Incremental Analysis Updates (IAU) \citep{bloom.takacs.ea_MWR96}. 
    4048IAU is a common technique used with 3D assimilation methods such as 3D-Var or OI. 
    4149IAU is used when \np{ln\_asmiau} is set to true. 
    4250 
    43 With IAU, the model state trajectory ${\bf x}$ in the assimilation window ($t_{0} \leq t_{i} \leq t_{N}$) 
     51With IAU, the model state trajectory ${\mathbf x}$ in the assimilation window ($t_{0} \leq t_{i} \leq t_{N}$) 
    4452is corrected by adding the analysis increments for temperature, salinity, horizontal velocity and SSH as 
    4553additional tendency terms to the prognostic equations: 
    4654\begin{align*} 
    4755  % \label{eq:wa_traj_iau} 
    48   {\bf x}^{a}(t_{i}) = M(t_{i}, t_{0})[{\bf x}^{b}(t_{0})] \; + \; F_{i} \delta \tilde{\bf x}^{a} 
     56  {\mathbf x}^{a}(t_{i}) = M(t_{i}, t_{0})[{\mathbf x}^{b}(t_{0})] \; + \; F_{i} \delta \tilde{\mathbf x}^{a} 
    4957\end{align*} 
    50 where $F_{i}$ is a weighting function for applying the increments $\delta\tilde{\bf x}^{a}$ defined such that 
     58where $F_{i}$ is a weighting function for applying the increments $\delta\tilde{\mathbf x}^{a}$ defined such that 
    5159$\sum_{i=1}^{N} F_{i}=1$. 
    52 ${\bf x}^b$ denotes the model initial state and ${\bf x}^a$ is the model state after the increments are applied. 
     60${\mathbf x}^b$ denotes the model initial state and ${\mathbf x}^a$ is the model state after the increments are applied. 
    5361To control the adjustment time of the model to the increment, 
    5462the increment can be applied over an arbitrary sub-window, $t_{m} \leq t_{i} \leq t_{n}$, 
     
    5664Typically the increments are spread evenly over the full window. 
    5765In addition, two different weighting functions have been implemented. 
    58 The first function employs constant weights,  
     66The first function (namelist option \np{niaufn} = 0) employs constant weights, 
    5967\begin{align} 
    6068  \label{eq:F1_i} 
     
    6270  =\left\{ 
    6371  \begin{array}{ll} 
    64     0     &    {\rm if} \; \; \; t_{i} < t_{m}                \\ 
    65     1/M &    {\rm if} \; \; \; t_{m} < t_{i} \leq t_{n} \\ 
    66     0     &    {\rm if} \; \; \; t_{i} > t_{n} 
     72    0     &    {\mathrm if} \; \; \; t_{i} < t_{m}                \\ 
     73    1/M &    {\mathrm if} \; \; \; t_{m} < t_{i} \leq t_{n} \\ 
     74    0     &    {\mathrm if} \; \; \; t_{i} > t_{n} 
    6775  \end{array} 
    68             \right.  
     76            \right. 
    6977\end{align} 
    7078where $M = m-n$. 
    71 The second function employs peaked hat-like weights in order to give maximum weight in the centre of the sub-window, 
     79The second function (namelist option \np{niaufn} = 1) employs peaked hat-like weights in order to give maximum weight in the centre of the sub-window, 
    7280with the weighting reduced linearly to a small value at the window end-points: 
    7381\begin{align} 
     
    7684  =\left\{ 
    7785  \begin{array}{ll} 
    78     0                           &    {\rm if} \; \; \; t_{i}       <     t_{m}                        \\ 
    79     \alpha \, i               &    {\rm if} \; \; \; t_{m}    \leq t_{i}    \leq   t_{M/2}   \\ 
    80     \alpha \, (M - i +1) &    {\rm if} \; \; \; t_{M/2}  <    t_{i}    \leq   t_{n}       \\ 
    81     0                            &   {\rm if} \; \; \; t_{i}        >    t_{n} 
     86    0                           &    {\mathrm if} \; \; \; t_{i}       <     t_{m}                        \\ 
     87    \alpha \, i               &    {\mathrm if} \; \; \; t_{m}    \leq t_{i}    \leq   t_{M/2}   \\ 
     88    \alpha \, (M - i +1) &    {\mathrm if} \; \; \; t_{M/2}  <    t_{i}    \leq   t_{n}       \\ 
     89    0                            &   {\mathrm if} \; \; \; t_{i}        >    t_{n} 
    8290  \end{array} 
    8391                                   \right. 
    8492\end{align} 
    85 where $\alpha^{-1} = \sum_{i=1}^{M/2} 2i$ and $M$ is assumed to be even.  
     93where $\alpha^{-1} = \sum_{i=1}^{M/2} 2i$ and $M$ is assumed to be even. 
    8694The weights described by \autoref{eq:F2_i} provide a smoother transition of the analysis trajectory from 
    8795one assimilation cycle to the next than that described by \autoref{eq:F1_i}. 
     
    92100\label{sec:ASM_div_dmp} 
    93101 
    94 The velocity increments may be initialized by the iterative application of a divergence damping operator. 
    95 In iteration step $n$ new estimates of velocity increments $u^{n}_I$ and $v^{n}_I$ are updated by: 
     102It is quite challenging for data assimilation systems to provide non-divergent velocity increments. 
     103Applying 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}. 
     104 
     105In iteration step $n$ (starting at $n=1$) new estimates of velocity increments $u^{n}_I$ and $v^{n}_I$ are updated by: 
     106 
    96107\begin{equation} 
    97108  \label{eq:asm_dmp} 
     
    105116  \right., 
    106117\end{equation} 
    107 where 
     118 
     119where the divergence is defined as 
     120 
    108121\[ 
    109122  % \label{eq:asm_div} 
     
    112125      +\delta_j \left[ {e_{1v}\,e_{3v}\,v^{n-1}_I} \right]} \right). 
    113126\] 
    114 By the application of \autoref{eq:asm_dmp} and \autoref{eq:asm_dmp} the divergence is filtered in each iteration, 
     127 
     128By the application of \autoref{eq:asm_dmp} the divergence is filtered in each iteration, 
    115129and the vorticity is left unchanged. 
    116130In the presence of coastal boundaries with zero velocity increments perpendicular to the coast 
     
    118132This type of the initialisation reduces the vertical velocity magnitude and 
    119133alleviates the problem of the excessive unphysical vertical mixing in the first steps of the model integration 
    120 \citep{Talagrand_JAS72, Dobricic_al_OS07}. 
     134\citep{talagrand_JAS72, dobricic.pinardi.ea_OS07}. 
    121135Diffusion coefficients are defined as $A_D = \alpha e_{1t} e_{2t}$, where $\alpha = 0.2$. 
    122136The divergence damping is activated by assigning to \np{nn\_divdmp} in the \textit{nam\_asminc} namelist 
    123137a value greater than zero. 
    124 By choosing this value to be of the order of 100 the increments in 
    125 the vertical velocity will be significantly reduced. 
     138This 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. 
    126139 
    127140 
     
    131144\label{sec:ASM_details} 
    132145 
    133 Here we show an example \ngn{namasm} namelist and the header of an example assimilation increments file on 
     146Here we show an example \ngn{nam\_asminc} namelist and the header of an example assimilation increments file on 
    134147the ORCA2 grid. 
    135148 
    136 %------------------------------------------namasm----------------------------------------------------- 
     149%------------------------------------------nam_asminc----------------------------------------------------- 
    137150% 
    138151\nlst{nam_asminc} 
  • NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_CONFIG.tex

    r10442 r11422  
    1818\label{sec:CFG_intro} 
    1919 
    20 The purpose of this part of the manual is to introduce the \NEMO reference configurations. 
     20The purpose of this part of the manual is to introduce the NEMO reference configurations. 
    2121These configurations are offered as means to explore various numerical and physical options, 
    2222thus allowing the user to verify that the code is performing in a manner consistent with that we are running. 
     
    2424The reference configurations also provide a sense for some of the options available in the code, 
    2525though by no means are all options exercised in the reference configurations. 
     26Configuration is defined manually through the \textit{namcfg} namelist variables. 
    2627 
    2728%------------------------------------------namcfg---------------------------------------------------- 
     
    3334% 1D model configuration 
    3435% ================================================================ 
    35 \section{C1D: 1D Water column model (\protect\key{c1d}) } 
     36\section[C1D: 1D Water column model (\texttt{\textbf{key\_c1d}})] 
     37{C1D: 1D Water column model (\protect\key{c1d})} 
    3638\label{sec:CFG_c1d} 
    3739 
    38 BE careful: to be re-written according to suppression of jpizoom and jpjzoom !!!! 
    39  
    40 The 1D model option simulates a stand alone water column within the 3D \NEMO system. 
     40The 1D model option simulates a stand alone water column within the 3D NEMO system. 
    4141It can be applied to the ocean alone or to the ocean-ice system and can include passive tracers or a biogeochemical model. 
    4242It is set up by defining the position of the 1D water column in the grid 
    43 (see \textit{CONFIG/SHARED/namelist\_ref} ).  
     43(see \textit{cfgs/SHARED/namelist\_ref}).  
    4444The 1D model is a very useful tool 
    4545\textit{(a)} to learn about the physics and numerical treatment of vertical mixing processes; 
     
    5050\textit{(d)} to produce extra diagnostics, without the large memory requirement of the full 3D model. 
    5151 
    52 The methodology is based on the use of the zoom functionality over the smallest possible domain: 
    53 a 3x3 domain centered on the grid point of interest, with some extra routines. 
    54 There is no need to define a new mesh, bathymetry, initial state or forcing, 
    55 since the 1D model will use those of the configuration it is a zoom of. 
    56 The chosen grid point is set in \textit{\ngn{namcfg}} namelist by 
    57 setting the \np{jpizoom} and \np{jpjzoom} parameters to the indices of the location of the chosen grid point. 
     52The methodology is based on the configuration of the smallest possible domain: 
     53a 3x3 domain with 75 vertical levels. 
    5854 
    5955The 1D model has some specifies. First, all the horizontal derivatives are assumed to be zero, 
    6056and second, the two components of the velocity are moved on a $T$-point.  
    61 Therefore, defining \key{c1d} changes five main things in the code behaviour:  
     57Therefore, defining \key{c1d} changes some things in the code behaviour:  
    6258\begin{description} 
    6359\item[(1)] 
    64   the lateral boundary condition routine (\rou{lbc\_lnk}) set the value of the central column of 
    65   the 3x3 domain is imposed over the whole domain; 
    66 \item[(3)] 
    67   a call to \rou{lbc\_lnk} is systematically done when reading input data (\ie in \mdl{iom}); 
    68 \item[(3)] 
    6960  a simplified \rou{stp} routine is used (\rou{stp\_c1d}, see \mdl{step\_c1d} module) in which 
    7061  both lateral tendancy terms and lateral physics are not called; 
    71 \item[(4)] 
     62\item[(2)] 
    7263  the vertical velocity is zero 
    7364  (so far, no attempt at introducing a Ekman pumping velocity has been made); 
    74 \item[(5)] 
     65\item[(3)] 
    7566  a simplified treatment of the Coriolis term is performed as $U$- and $V$-points are the same 
    7667  (see \mdl{dyncor\_c1d}). 
    7768\end{description} 
    78 All the relevant \textit{\_c1d} modules can be found in the NEMOGCM/NEMO/OPA\_SRC/C1D directory of 
    79 the \NEMO distribution. 
     69All the relevant \textit{\_c1d} modules can be found in the src/OCE/C1D directory of 
     70the NEMO distribution. 
    8071 
    8172% to be added:  a test case on the yearlong Ocean Weather Station (OWS) Papa dataset of Martin (1985) 
     
    8879 
    8980The ORCA family is a series of global ocean configurations that are run together with 
    90 the LIM sea-ice model (ORCA-LIM) and possibly with PISCES biogeochemical model (ORCA-LIM-PISCES), 
    91 using various resolutions. 
    92 An appropriate namelist is available in \path{CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg} for ORCA2. 
     81the SI3 model (ORCA-ICE) and possibly with PISCES biogeochemical model (ORCA-ICE-PISCES). 
     82An appropriate namelist is available in \path{cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg} for ORCA2. 
    9383The domain of ORCA2 configuration is defined in \ifile{ORCA\_R2\_zps\_domcfg} file, 
    94 this file is available in tar file in the wiki of NEMO: \\ 
    95 https://forge.ipsl.jussieu.fr/nemo/wiki/Users/ReferenceConfigurations/ORCA2\_LIM3\_PISCES \\ 
     84this file is available in tar file on the NEMO community zenodo platform: \\ 
     85https://doi.org/10.5281/zenodo.2640723 
     86 
    9687In this namelist\_cfg the name of domain input file is set in \ngn{namcfg} block of namelist.  
    9788 
     
    9990\begin{figure}[!t] 
    10091  \begin{center} 
    101     \includegraphics[width=0.98\textwidth]{Fig_ORCA_NH_mesh} 
     92    \includegraphics[width=\textwidth]{Fig_ORCA_NH_mesh} 
    10293    \caption{ 
    10394      \protect\label{fig:MISC_ORCA_msh} 
     
    10697      The two "north pole" are the foci of a series of embedded ellipses (blue curves) which 
    10798      are determined analytically and form the i-lines of the ORCA mesh (pseudo latitudes). 
    108       Then, following \citet{Madec_Imbard_CD96}, the normal to the series of ellipses (red curves) is computed which 
     99      Then, following \citet{madec.imbard_CD96}, the normal to the series of ellipses (red curves) is computed which 
    109100      provides the j-lines of the mesh (pseudo longitudes). 
    110101    } 
     
    119110\label{subsec:CFG_orca_grid} 
    120111 
    121 The ORCA grid is a tripolar is based on the semi-analytical method of \citet{Madec_Imbard_CD96}. 
     112The ORCA grid is a tripolar grid based on the semi-analytical method of \citet{madec.imbard_CD96}. 
    122113It allows to construct a global orthogonal curvilinear ocean mesh which has no singularity point inside 
    123114the computational domain since two north mesh poles are introduced and placed on lands. 
     
    131122\begin{figure}[!tbp] 
    132123  \begin{center} 
    133     \includegraphics[width=1.0\textwidth]{Fig_ORCA_NH_msh05_e1_e2} 
    134     \includegraphics[width=0.80\textwidth]{Fig_ORCA_aniso} 
     124    \includegraphics[width=\textwidth]{Fig_ORCA_NH_msh05_e1_e2} 
     125    \includegraphics[width=\textwidth]{Fig_ORCA_aniso} 
    135126    \caption { 
    136127      \protect\label{fig:MISC_ORCA_e1e2} 
     
    158149 
    159150% ------------------------------------------------------------------------------------------------------------- 
    160 %       ORCA-LIM(-PISCES) configurations 
     151%       ORCA-ICE(-PISCES) configurations 
    161152% ------------------------------------------------------------------------------------------------------------- 
    162153\subsection{ORCA pre-defined resolution} 
     
    199190The ORCA\_R2 configuration has the following specificity: starting from a 2\deg~ORCA mesh, 
    200191local mesh refinements were applied to the Mediterranean, Red, Black and Caspian Seas, 
    201 so that the resolution is 1\deg \time 1\deg there. 
     192so that the resolution is 1\deg~ there. 
    202193A local transformation were also applied with in the Tropics in order to refine the meridional resolution up to 
    203 0.5\deg at the Equator. 
     1940.5\deg~ at the Equator. 
    204195 
    205196The ORCA\_R1 configuration has only a local tropical transformation to refine the meridional resolution up to 
     
    211202For ORCA\_R1 and R025, setting the configuration key to 75 allows to use 75 vertical levels, otherwise 46 are used. 
    212203In the other ORCA configurations, 31 levels are used 
    213 (see \autoref{tab:orca_zgr} %\sfcomment{HERE I need to put new table for ORCA2 values} and \autoref{fig:zgr}). 
    214  
    215 Only the ORCA\_R2 is provided with all its input files in the \NEMO distribution. 
    216 It is very similar to that used as part of the climate model developed at IPSL for the 4th IPCC assessment of 
    217 climate change (Marti et al., 2009). 
    218 It is also the basis for the \NEMO contribution to the Coordinate Ocean-ice Reference Experiments (COREs) 
    219 documented in \citet{Griffies_al_OM09}.  
     204(see \autoref{tab:orca_zgr}). %\sfcomment{HERE I need to put new table for ORCA2 values} and \autoref{fig:zgr}). 
     205 
     206Only the ORCA\_R2 is provided with all its input files in the NEMO distribution. 
     207%It is very similar to that used as part of the climate model developed at IPSL for the 4th IPCC assessment of 
     208%climate change (Marti et al., 2009). 
     209%It is also the basis for the \NEMO contribution to the Coordinate Ocean-ice Reference Experiments (COREs) 
     210%documented in \citet{griffies.biastoch.ea_OM09}.  
    220211 
    221212This version of ORCA\_R2 has 31 levels in the vertical, with the highest resolution (10m) in the upper 150m 
    222213(see \autoref{tab:orca_zgr} and \autoref{fig:zgr}).  
    223214The bottom topography and the coastlines are derived from the global atlas of Smith and Sandwell (1997).  
    224 The default forcing uses the boundary forcing from \citet{Large_Yeager_Rep04} (see \autoref{subsec:SBC_blk_core}),  
     215The default forcing uses the boundary forcing from \citet{large.yeager_rpt04} (see \autoref{subsec:SBC_blk_core}),  
    225216which was developed for the purpose of running global coupled ocean-ice simulations without 
    226217an interactive atmosphere. 
    227 This \citet{Large_Yeager_Rep04} dataset is available through 
     218This \citet{large.yeager_rpt04} dataset is available through 
    228219the \href{http://nomads.gfdl.noaa.gov/nomads/forms/mom4/CORE.html}{GFDL web site}. 
    229 The "normal year" of \citet{Large_Yeager_Rep04} has been chosen of the \NEMO distribution since release v3.3.  
    230  
    231 ORCA\_R2 pre-defined configuration can also be run with an AGRIF zoom over the Agulhas current area 
    232 (\key{agrif} defined) and, by setting the appropriate variables, see \path{CONFIG/SHARED/namelist_ref}. 
     220The "normal year" of \citet{large.yeager_rpt04} has been chosen of the NEMO distribution since release v3.3.  
     221 
     222ORCA\_R2 pre-defined configuration can also be run with multiply online nested zooms (\ie with AGRIF, \key{agrif} defined). This is available as the AGRIF\_DEMO configuration that can be found in the \path{cfgs/AGRIF_DEMO/} directory. 
     223 
    233224A regional Arctic or peri-Antarctic configuration is extracted from an ORCA\_R2 or R05 configurations using 
    234225sponge layers at open boundaries.  
     
    237228%       GYRE family: double gyre basin 
    238229% ------------------------------------------------------------------------------------------------------------- 
    239 \section{GYRE family: double gyre basin } 
     230\section{GYRE family: double gyre basin} 
    240231\label{sec:CFG_gyre} 
    241232 
    242 The GYRE configuration \citep{Levy_al_OM10} has been built to 
     233The GYRE configuration \citep{levy.klein.ea_OM10} has been built to 
    243234simulate the seasonal cycle of a double-gyre box model. 
    244 It consists in an idealized domain similar to that used in the studies of \citet{Drijfhout_JPO94} and 
    245 \citet{Hazeleger_Drijfhout_JPO98, Hazeleger_Drijfhout_JPO99, Hazeleger_Drijfhout_JGR00, Hazeleger_Drijfhout_JPO00}, 
     235It consists in an idealized domain similar to that used in the studies of \citet{drijfhout_JPO94} and 
     236\citet{hazeleger.drijfhout_JPO98, hazeleger.drijfhout_JPO99, hazeleger.drijfhout_JGR00, hazeleger.drijfhout_JPO00}, 
    246237over which an analytical seasonal forcing is applied. 
    247238This allows to investigate the spontaneous generation of a large number of interacting, transient mesoscale eddies  
    248239and their contribution to the large scale circulation.  
    249240 
     241The GYRE configuration run together with the PISCES biogeochemical model (GYRE-PISCES). 
    250242The domain geometry is a closed rectangular basin on the $\beta$-plane centred at $\sim$ 30\deg{N} and 
    251243rotated by 45\deg, 3180~km long, 2120~km wide and 4~km deep (\autoref{fig:MISC_strait_hand}). 
     
    253245The configuration is meant to represent an idealized North Atlantic or North Pacific basin. 
    254246The circulation is forced by analytical profiles of wind and buoyancy fluxes. 
    255 The applied forcings vary seasonally in a sinusoidal manner between winter and summer extrema \citep{Levy_al_OM10}.  
     247The applied forcings vary seasonally in a sinusoidal manner between winter and summer extrema \citep{levy.klein.ea_OM10}.  
    256248The wind stress is zonal and its curl changes sign at 22\deg{N} and 36\deg{N}. 
    257249It forces a subpolar gyre in the north, a subtropical gyre in the wider part of the domain and 
     
    266258The GYRE configuration is set like an analytical configuration. 
    267259Through \np{ln\_read\_cfg}\forcode{ = .false.} in \textit{namcfg} namelist defined in 
    268 the reference configuration \path{CONFIG/GYRE/EXP00/namelist_cfg} 
     260the reference configuration \path{cfgs/GYRE_PISCES/EXPREF/namelist_cfg} 
    269261analytical definition of grid in GYRE is done in usrdef\_hrg, usrdef\_zgr routines. 
    270262Its horizontal resolution (and thus the size of the domain) is determined by 
    271263setting \np{nn\_GYRE} in \ngn{namusr\_def}: \\ 
     264 
    272265\np{jpiglo} $= 30 \times$ \np{nn\_GYRE} + 2   \\ 
     266 
    273267\np{jpjglo} $= 20 \times$ \np{nn\_GYRE} + 2   \\ 
     268 
    274269Obviously, the namelist parameters have to be adjusted to the chosen resolution, 
    275 see the Configurations pages on the NEMO web site (Using NEMO\/Configurations). 
     270see the Configurations pages on the NEMO web site (NEMO Configurations). 
    276271In the vertical, GYRE uses the default 30 ocean levels (\jp{jpk}\forcode{ = 31}) (\autoref{fig:zgr}). 
    277272 
     
    281276even though the physical integrity of the solution can be compromised. 
    282277Benchmark is activate via \np{ln\_bench}\forcode{ = .true.} in \ngn{namusr\_def} in 
    283 namelist \path{CONFIG/GYRE/EXP00/namelist_cfg}. 
     278namelist \path{cfgs/GYRE_PISCES/EXPREF/namelist_cfg}. 
    284279 
    285280%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    286281\begin{figure}[!t] 
    287282  \begin{center} 
    288     \includegraphics[width=1.0\textwidth]{Fig_GYRE} 
     283    \includegraphics[width=\textwidth]{Fig_GYRE} 
    289284    \caption{ 
    290285      \protect\label{fig:GYRE} 
    291286      Snapshot of relative vorticity at the surface of the model domain in GYRE R9, R27 and R54. 
    292       From \citet{Levy_al_OM10}. 
     287      From \citet{levy.klein.ea_OM10}. 
    293288    } 
    294289  \end{center} 
     
    304299The AMM, Atlantic Margins Model, is a regional model covering the Northwest European Shelf domain on 
    305300a regular lat-lon grid at approximately 12km horizontal resolution. 
    306 The appropriate \textit{\&namcfg} namelist  is available in \textit{CONFIG/AMM12/EXP00/namelist\_cfg}. 
     301The appropriate \textit{\&namcfg} namelist  is available in \textit{cfgs/AMM12/EXPREF/namelist\_cfg}. 
    307302It is used to build the correct dimensions of the AMM domain. 
    308303 
    309304This configuration tests several features of NEMO functionality specific to the shelf seas. 
    310 In particular, the AMM uses $S$-coordinates in the vertical rather than $z$-coordinates and 
    311 is forced with tidal lateral boundary conditions using a flather boundary condition from the BDY module. 
    312 The AMM configuration uses the GLS (\key{zdfgls}) turbulence scheme, 
    313 the VVL non-linear free surface(\key{vvl}) and time-splitting (\key{dynspg\_ts}). 
     305In particular, the AMM uses $s$-coordinates in the vertical rather than $z$-coordinates and 
     306is forced with tidal lateral boundary conditions using a Flather boundary condition from the BDY module. 
     307Also specific to the AMM configuration is the use of the GLS turbulence scheme (\np{ln\_zdfgls} \forcode{= .true.}). 
    314308 
    315309In addition to the tidal boundary condition the model may also take open boundary conditions from 
  • NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_DIA.tex

    r10509 r11422  
    1010\minitoc 
    1111 
     12\vfill 
     13\begin{figure}[b] 
     14\subsubsection*{Changes record} 
     15\begin{tabular}{l||l|m{0.65\linewidth}} 
     16    Release   & Author        & Modifications \\ 
     17    {\em 4.0} & {\em Mirek Andrejczuk, Massimiliano Drudi} & {\em }  \\ 
     18    {\em }      & {\em Dorotea Iovino, Nicolas Martin} & {\em }  \\ 
     19    {\em 3.6} & {\em Gurvan Madec, Sebastien Masson } & {\em }  \\ 
     20    {\em 3.4} & {\em Gurvan Madec, Rachid Benshila, Andrew Coward } & {\em }  \\  
     21    {\em }      & {\em Christian Ethe, Sebastien Masson } & {\em }  \\  
     22\end{tabular} 
     23\end{figure}  
     24 
    1225\newpage 
    1326 
     
    1528%       Old Model Output  
    1629% ================================================================ 
    17 \section{Old model output (default)} 
     30\section{Model output} 
    1831\label{sec:DIA_io_old} 
    1932 
     
    2538the same run performed in one step. 
    2639It should be noted that this requires that the restart file contains two consecutive time steps for 
    27 all the prognostic variables, and that it is saved in the same binary format as the one used by the computer that 
    28 is to read it (in particular, 32 bits binary IEEE format must not be used for this file). 
     40all the prognostic variables. 
    2941 
    3042The output listing and file(s) are predefined but should be checked and eventually adapted to the user's needs. 
     
    3850the writing work is distributed over the processors in massively parallel computing. 
    3951A complete description of the use of this I/O server is presented in the next section. 
    40  
    41 By default, \key{iomput} is not defined, 
    42 NEMO produces NetCDF with the old IOIPSL library which has been kept for compatibility and its easy installation. 
    43 However, the IOIPSL library is quite inefficient on parallel machines and, since version 3.2, 
    44 many diagnostic options have been added presuming the use of \key{iomput}. 
    45 The usefulness of the default IOIPSL-based option is expected to reduce with each new release. 
    46 If \key{iomput} is not defined, output files and content are defined in the \mdl{diawri} module and 
    47 contain mean (or instantaneous if \key{diainstant} is defined) values over a regular period of 
    48 nn\_write time-steps (namelist parameter).  
    4952 
    5053%\gmcomment{                    % start of gmcomment 
     
    9194in a very easy way. 
    9295All details of iomput functionalities are listed in the following subsections. 
    93 Examples of the XML files that control the outputs can be found in: \path{NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef.xml}, 
    94 \path{NEMOGCM/CONFIG/SHARED/field_def.xml} and \path{NEMOGCM/CONFIG/SHARED/domain_def.xml}. \\ 
     96Examples of the XML files that control the outputs can be found in:  
     97\path{cfgs/ORCA2_ICE_PISCES/EXPREF/iodef.xml}, 
     98\path{cfgs/SHARED/field_def_nemo-oce.xml}, 
     99\path{cfgs/SHARED/field_def_nemo-pisces.xml}, 
     100\path{cfgs/SHARED/field_def_nemo-ice.xml} and \path{cfgs/SHARED/domain_def_nemo.xml}. \\ 
    95101 
    96102The second functionality targets output performance when running in parallel (\key{mpp\_mpi}). 
     
    101107 
    102108In version 3.6, the iom\_put interface depends on 
    103 an external code called \href{https://forge.ipsl.jussieu.fr/ioserver/browser/XIOS/branchs/xios-1.0}{XIOS-1.0}  
     109an external code called \href{https://forge.ipsl.jussieu.fr/ioserver/browser/XIOS/branchs/xios-2.5}{XIOS-2.5}  
    104110(use of revision 618 or higher is required). 
    105111This new IO server can take advantage of the parallel I/O functionality of NetCDF4 to 
     
    168174\xmlline|<variable id="using_server" type="bool"></variable>| 
    169175 
    170 The {\tt using\_server} setting determines whether or not the server will be used in \textit{attached mode} 
    171 (as a library) [{\tt> false <}] or in \textit{detached mode} 
    172 (as an external executable on N additional, dedicated cpus) [{\tt > true <}]. 
     176The {\ttfamily using\_server} setting determines whether or not the server will be used in \textit{attached mode} 
     177(as a library) [{\ttfamily> false <}] or in \textit{detached mode} 
     178(as an external executable on N additional, dedicated cpus) [{\ttfamily > true <}]. 
    173179The \textit{attached mode} is simpler to use but much less efficient for massively parallel applications. 
    174180The type of each file can be either ''multiple\_file'' or ''one\_file''. 
     
    207213\subsubsection{Control of XIOS: the context in iodef.xml} 
    208214 
    209 As well as the {\tt using\_server} flag, other controls on the use of XIOS are set in the XIOS context in iodef.xml. 
     215As well as the {\ttfamily using\_server} flag, other controls on the use of XIOS are set in the XIOS context in iodef.xml. 
    210216See the XML basics section below for more details on XML syntax and rules. 
    211217 
     
    257263See the installation guide on the \href{http://forge.ipsl.jussieu.fr/ioserver/wiki}{XIOS} wiki for help and guidance. 
    258264NEMO will need to link to the compiled XIOS library. 
    259 The \href{https://forge.ipsl.jussieu.fr/nemo/wiki/Users/ModelInterfacing/InputsOutputs#Inputs-OutputsusingXIOS} 
    260 {XIOS with NEMO} guide provides an example illustration of how this can be achieved. 
     265The \href{https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/guide/html/install.html#extract-and-install-xios} 
     266{Extract and install XIOS} guide provides an example illustration of how this can be achieved. 
    261267 
    262268\subsubsection{Add your own outputs} 
     
    269275\begin{enumerate} 
    270276\item[1.] 
    271   in NEMO code, add a \forcode{CALL iom\_put( 'identifier', array )} where you want to output a 2D or 3D array. 
     277  in NEMO code, add a \forcode{CALL iom_put( 'identifier', array )} where you want to output a 2D or 3D array. 
    272278\item[2.] 
    273279  If necessary, add \forcode{USE iom ! I/O manager library} to the list of used modules in 
     
    442448\xmlline|<context src="./nemo_def.xml" />| 
    443449  
    444 \noindent In NEMO, by default, the field and domain definition is done in 2 separate files: 
    445 \path{NEMOGCM/CONFIG/SHARED/field_def.xml} and \path{NEMOGCM/CONFIG/SHARED/domain_def.xml} that 
     450\noindent In NEMO, by default, the field definition is done in 3 separate files ( 
     451\path{cfgs/SHARED/field_def_nemo-oce.xml}, 
     452\path{cfgs/SHARED/field_def_nemo-pisces.xml} and 
     453\path{cfgs/SHARED/field_def_nemo-ice.xml} ) and the  domain definition is done in another file ( \path{cfgs/SHARED/domain_def_nemo.xml} ) 
     454that 
    446455are included in the main iodef.xml file through the following commands: 
    447456\begin{xmllines} 
    448 <field_definition src="./field_def.xml" /> 
    449 <domain_definition src="./domain_def.xml" /> 
     457<context id="nemo" src="./context_nemo.xml"/>  
    450458\end{xmllines} 
    451459 
     
    508516 
    509517Secondly, the group can be used to replace a list of elements. 
    510 Several examples of groups of fields are proposed at the end of the file \path{CONFIG/SHARED/field_def.xml}. 
     518Several examples of groups of fields are proposed at the end of the XML field files ( 
     519\path{cfgs/SHARED/field_def_nemo-oce.xml}, 
     520\path{cfgs/SHARED/field_def_nemo-pisces.xml} and 
     521\path{cfgs/SHARED/field_def_nemo-ice.xml} ) . 
    511522For example, a short list of the usual variables related to the U grid: 
    512523 
     
    514525<field_group id="groupU" > 
    515526   <field field_ref="uoce"  /> 
    516    <field field_ref="suoce" /> 
     527   <field field_ref="ssu" /> 
    517528   <field field_ref="utau"  /> 
    518529</field_group> 
     
    538549the tag family domain. 
    539550It must therefore be done in the domain part of the XML file.  
    540 For example, in \path{CONFIG/SHARED/domain_def.xml}, we provide the following example of a definition of  
     551For example, in \path{cfgs/SHARED/domain_def.xml}, we provide the following example of a definition of  
    541552a 5 by 5 box with the bottom left corner at point (10,10). 
    542553 
    543554\begin{xmllines} 
    544 <domain_group id="grid_T"> 
    545    <domain id="myzoom" zoom_ibegin="10" zoom_jbegin="10" zoom_ni="5" zoom_nj="5" /> 
     555<domain id="myzoomT" domain_ref="grid_T"> 
     556   <zoom_domain ibegin="10" jbegin="10" ni="5" nj="5" /> 
    546557\end{xmllines} 
    547558 
     
    551562\begin{xmllines} 
    552563<file id="myfile_vzoom" output_freq="1d" > 
    553    <field field_ref="toce" domain_ref="myzoom"/> 
     564   <field field_ref="toce" domain_ref="myzoomT"/> 
    554565</file> 
    555566\end{xmllines} 
     
    576587\subsubsection{Define vertical zooms} 
    577588 
    578 Vertical zooms are defined through the attributs zoom\_begin and zoom\_end of the tag family axis. 
     589Vertical zooms are defined through the attributs zoom\_begin and zoom\_n of the tag family axis. 
    579590It must therefore be done in the axis part of the XML file. 
    580 For example, in \path{NEMOGCM/CONFIG/ORCA2_LIM/iodef_demo.xml}, we provide the following example: 
    581  
    582 \begin{xmllines} 
    583 <axis_group id="deptht" long_name="Vertical T levels" unit="m" positive="down" > 
    584    <axis id="deptht" /> 
    585    <axis id="deptht_myzoom" zoom_begin="1" zoom_end="10" /> 
     591For example, in \path{cfgs/ORCA2_ICE_PISCES/EXPREF/iodef_demo.xml}, we provide the following example: 
     592 
     593\begin{xmllines} 
     594<axis_definition> 
     595   <axis id="deptht" long_name="Vertical T levels" unit="m" positive="down" /> 
     596   <axis id="deptht_zoom" azix_ref="deptht" > 
     597      <zoom_axis zoom_begin="1" zoom_n="10" /> 
     598   </axis> 
    586599\end{xmllines} 
    587600 
     
    765778\end{xmllines} 
    766779 
    767 Note that, then the code is crashing, writting real4 variables forces a numerical convection from  
     780Note that, then the code is crashing, writting real4 variables forces a numerical conversion from  
    768781real8 to real4 which will create an internal error in NetCDF and will avoid the creation of the output files. 
    769782Forcing double precision outputs with prec="8" (for example in the field\_definition) will avoid this problem. 
     
    938951    \hline 
    939952  \end{tabularx} 
    940   \caption{Field tags ("\tt{field\_*}")} 
     953  \caption{Field tags ("\ttfamily{field\_*}")} 
    941954\end{table} 
    942955 
     
    974987    \hline 
    975988  \end{tabularx} 
    976   \caption{File tags ("\tt{file\_*}")} 
     989  \caption{File tags ("\ttfamily{file\_*}")} 
    977990\end{table} 
    978991 
     
    10071020    \hline 
    10081021  \end{tabularx} 
    1009   \caption{Axis tags ("\tt{axis\_*}")} 
     1022  \caption{Axis tags ("\ttfamily{axis\_*}")} 
    10101023\end{table} 
    10111024 
     
    10401053    \hline 
    10411054  \end{tabularx} 
    1042   \caption{Domain tags ("\tt{domain\_*)}"} 
     1055  \caption{Domain tags ("\ttfamily{domain\_*)}"} 
    10431056\end{table} 
    10441057 
     
    10731086    \hline 
    10741087  \end{tabularx} 
    1075   \caption{Grid tags ("\tt{grid\_*}")} 
     1088  \caption{Grid tags ("\ttfamily{grid\_*}")} 
    10761089\end{table} 
    10771090 
     
    11141127    \hline 
    11151128  \end{tabularx} 
    1116   \caption{Reference attributes ("\tt{*\_ref}")} 
     1129  \caption{Reference attributes ("\ttfamily{*\_ref}")} 
    11171130\end{table} 
    11181131 
     
    11501163    \hline 
    11511164  \end{tabularx} 
    1152   \caption{Domain attributes ("\tt{zoom\_*}")} 
     1165  \caption{Domain attributes ("\ttfamily{zoom\_*}")} 
    11531166\end{table} 
    11541167 
     
    13181331\subsection{CF metadata standard compliance} 
    13191332 
    1320 Output from the XIOS-1.0 IO server is compliant with  
     1333Output from the XIOS IO server is compliant with  
    13211334\href{http://cfconventions.org/Data/cf-conventions/cf-conventions-1.5/build/cf-conventions.html}{version 1.5} of 
    13221335the CF metadata standard.  
     
    13321345%       NetCDF4 support 
    13331346% ================================================================ 
    1334 \section{NetCDF4 support (\protect\key{netcdf4})} 
     1347\section[NetCDF4 support (\texttt{\textbf{key\_netcdf4}})] 
     1348{NetCDF4 support (\protect\key{netcdf4})} 
    13351349\label{sec:DIA_nc4} 
    13361350 
     
    13401354Chunking and compression can lead to significant reductions in file sizes for a small runtime overhead. 
    13411355For a fuller discussion on chunking and other performance issues the reader is referred to 
    1342 the NetCDF4 documentation found \href{http://www.unidata.ucar.edu/software/netcdf/docs/netcdf.html#Chunking}{here}. 
     1356the NetCDF4 documentation found \href{https://www.unidata.ucar.edu/software/netcdf/docs/netcdf_perf_chunking.html}{here}. 
    13431357 
    13441358The new features are only available when the code has been linked with a NetCDF4 library 
     
    13891403\end{forlines} 
    13901404 
    1391 \noindent for a standard ORCA2\_LIM configuration gives chunksizes of {\small\tt 46x38x1} respectively in 
    1392 the mono-processor case (\ie global domain of {\small\tt 182x149x31}). 
     1405\noindent for a standard ORCA2\_LIM configuration gives chunksizes of {\small\ttfamily 46x38x1} respectively in 
     1406the mono-processor case (\ie global domain of {\small\ttfamily 182x149x31}). 
    13931407An illustration of the potential space savings that NetCDF4 chunking and compression provides is given in  
    13941408table \autoref{tab:NC4} which compares the results of two short runs of the ORCA2\_LIM reference configuration with 
     
    14501464%       Tracer/Dynamics Trends 
    14511465% ------------------------------------------------------------------------------------------------------------- 
    1452 \section{Tracer/Dynamics trends  (\protect\ngn{namtrd})} 
     1466\section[Tracer/Dynamics trends (\texttt{namtrd})] 
     1467{Tracer/Dynamics trends (\protect\ngn{namtrd})} 
    14531468\label{sec:DIA_trd} 
    14541469 
     
    14621477(\ie at the end of each $dyn\cdots.F90$ and/or $tra\cdots.F90$ routines). 
    14631478This capability is controlled by options offered in \ngn{namtrd} namelist. 
    1464 Note that the output are done with xIOS, and therefore the \key{IOM} is required. 
     1479Note that the output are done with XIOS, and therefore the \key{iomput} is required. 
    14651480 
    14661481What is done depends on the \ngn{namtrd} logical set to \forcode{.true.}: 
     
    14881503 
    14891504Note that the mixed layer tendency diagnostic can also be used on biogeochemical models via  
    1490 the \key{trdtrc} and \key{trdmld\_trc} CPP keys. 
     1505the \key{trdtrc} and \key{trdmxl\_trc} CPP keys. 
    14911506 
    14921507\textbf{Note that} in the current version (v3.6), many changes has been introduced but not fully tested. 
    14931508In particular, options associated with \np{ln\_dyn\_mxl}, \np{ln\_vor\_trd}, and \np{ln\_tra\_mxl} are not working, 
    1494 and none of the options have been tested with variable volume (\ie \key{vvl} defined). 
     1509and none of the options have been tested with variable volume (\ie \np{ln\_linssh}\forcode{ = .true.}). 
    14951510 
    14961511% ------------------------------------------------------------------------------------------------------------- 
    14971512%       On-line Floats trajectories 
    14981513% ------------------------------------------------------------------------------------------------------------- 
    1499 \section{FLO: On-Line Floats trajectories (\protect\key{floats})} 
     1514\section[FLO: On-Line Floats trajectories (\texttt{\textbf{key\_floats}})] 
     1515{FLO: On-Line Floats trajectories (\protect\key{floats})} 
    15001516\label{sec:FLO} 
    15011517%--------------------------------------------namflo------------------------------------------------------- 
     
    15061522The on-line computation of floats advected either by the three dimensional velocity field or constraint to 
    15071523remain at a given depth ($w = 0$ in the computation) have been introduced in the system during the CLIPPER project. 
    1508 Options are defined by \ngn{namflo} namelis variables. 
    1509 The algorithm used is based either on the work of \cite{Blanke_Raynaud_JPO97} (default option), 
     1524Options are defined by \ngn{namflo} namelist variables. 
     1525The algorithm used is based either on the work of \cite{blanke.raynaud_JPO97} (default option), 
    15101526or on a $4^th$ Runge-Hutta algorithm (\np{ln\_flork4}\forcode{ = .true.}). 
    1511 Note that the \cite{Blanke_Raynaud_JPO97} algorithm have the advantage of providing trajectories which 
     1527Note that the \cite{blanke.raynaud_JPO97} algorithm have the advantage of providing trajectories which 
    15121528are consistent with the numeric of the code, so that the trajectories never intercept the bathymetry. 
    15131529 
     
    15191535In case of Ariane convention, input filename is \np{init\_float\_ariane}. 
    15201536Its format is: \\ 
    1521 {\scriptsize \texttt{I J K nisobfl itrash itrash}} 
     1537{\scriptsize \texttt{I J K nisobfl itrash}} 
    15221538 
    15231539\noindent with: 
     
    15771593In that case, output filename is trajec\_float. 
    15781594 
    1579 Another possiblity of writing format is Netcdf (\np{ln\_flo\_ascii}\forcode{ = .false.}). 
    1580 There are 2 possibilities: 
    1581  
    1582 - if (\key{iomput}) is used, outputs are selected in  iodef.xml. 
     1595Another possiblity of writing format is Netcdf (\np{ln\_flo\_ascii}\forcode{ = .false.}) with 
     1596\key{iomput} and outputs selected in iodef.xml. 
    15831597Here it is an example of specification to put in files description section: 
    15841598 
     
    15971611\end{xmllines} 
    15981612 
    1599  -  if (\key{iomput}) is not used, a file called \ifile{trajec\_float} will be created by IOIPSL library. 
    1600  
    1601  See also \href{http://stockage.univ-brest.fr/~grima/Ariane/}{here} the web site describing the off-line use of 
    1602  this marvellous diagnostic tool. 
    16031613 
    16041614% ------------------------------------------------------------------------------------------------------------- 
    16051615%       Harmonic analysis of tidal constituents 
    16061616% ------------------------------------------------------------------------------------------------------------- 
    1607 \section{Harmonic analysis of tidal constituents (\protect\key{diaharm}) } 
     1617\section[Harmonic analysis of tidal constituents (\texttt{\textbf{key\_diaharm}})] 
     1618{Harmonic analysis of tidal constituents (\protect\key{diaharm})} 
    16081619\label{sec:DIA_diag_harm} 
    16091620 
    1610 %------------------------------------------namdia_harm---------------------------------------------------- 
     1621%------------------------------------------nam_diaharm---------------------------------------------------- 
    16111622% 
    16121623\nlst{nam_diaharm} 
     
    16161627This on-line Harmonic analysis is actived with \key{diaharm}. 
    16171628 
    1618 Some parameters are available in namelist \ngn{namdia\_harm}: 
     1629Some parameters are available in namelist \ngn{nam\_diaharm}: 
    16191630 
    16201631 - \np{nit000\_han} is the first time step used for harmonic analysis 
     
    16521663%       Sections transports 
    16531664% ------------------------------------------------------------------------------------------------------------- 
    1654 \section{Transports across sections (\protect\key{diadct}) } 
     1665\section[Transports across sections (\texttt{\textbf{key\_diadct}})] 
     1666{Transports across sections (\protect\key{diadct})} 
    16551667\label{sec:DIA_diag_dct} 
    16561668 
     
    16641676 
    16651677Each section is defined by the coordinates of its 2 extremities. 
    1666 The pathways between them are contructed using tools which can be found in \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT} 
    1667 and are written in a binary file \texttt{section\_ijglobal.diadct\_ORCA2\_LIM} which is later read in by 
     1678The pathways between them are contructed using tools which can be found in \texttt{tools/SECTIONS\_DIADCT} 
     1679and are written in a binary file \texttt{section\_ijglobal.diadct} which is later read in by 
    16681680NEMO to compute on-line transports. 
    16691681 
     
    16841696\subsubsection{Creating a binary file containing the pathway of each section} 
    16851697 
    1686 In \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT/run}, 
     1698In \texttt{tools/SECTIONS\_DIADCT/run}, 
    16871699the file \textit{ {list\_sections.ascii\_global}} contains a list of all the sections that are to be computed 
    16881700(this list of sections is based on MERSEA project metrics). 
     
    17331745   
    17341746 The script \texttt{job.ksh} computes the pathway for each section and creates a binary file 
    1735  \texttt{section\_ijglobal.diadct\_ORCA2\_LIM} which is read by NEMO. \\ 
     1747 \texttt{section\_ijglobal.diadct} which is read by NEMO. \\ 
    17361748 
    17371749 It is possible to use this tools for new configuations: \texttt{job.ksh} has to be updated with 
     
    18091821The steric effect is therefore not explicitely represented. 
    18101822This approximation does not represent a serious error with respect to the flow field calculated by the model 
    1811 \citep{Greatbatch_JGR94}, but extra attention is required when investigating sea level, 
     1823\citep{greatbatch_JGR94}, but extra attention is required when investigating sea level, 
    18121824as steric changes are an important contribution to local changes in sea level on seasonal and climatic time scales. 
    18131825This is especially true for investigation into sea level rise due to global warming. 
    18141826 
    18151827Fortunately, the steric contribution to the sea level consists of a spatially uniform component that 
    1816 can be diagnosed by considering the mass budget of the world ocean \citep{Greatbatch_JGR94}. 
     1828can be diagnosed by considering the mass budget of the world ocean \citep{greatbatch_JGR94}. 
    18171829In order to better understand how global mean sea level evolves and thus how the steric sea level can be diagnosed, 
    18181830we compare, in the following, the non-Boussinesq and Boussinesq cases. 
     
    18881900the ocean surface, not by changes in mean mass of the ocean: the steric effect is missing in a Boussinesq fluid. 
    18891901 
    1890 Nevertheless, following \citep{Greatbatch_JGR94}, the steric effect on the volume can be diagnosed by 
     1902Nevertheless, following \citep{greatbatch_JGR94}, the steric effect on the volume can be diagnosed by 
    18911903considering the mass budget of the ocean.  
    18921904The apparent changes in $\mathcal{M}$, mass of the ocean, which are not induced by surface mass flux 
    18931905must be compensated by a spatially uniform change in the mean sea level due to expansion/contraction of the ocean 
    1894 \citep{Greatbatch_JGR94}. 
     1906\citep{greatbatch_JGR94}. 
    18951907In others words, the Boussinesq mass, $\mathcal{M}_o$, can be related to $\mathcal{M}$, 
    18961908the total mass of the ocean seen by the Boussinesq model, via the steric contribution to the sea level, 
     
    19241936This value is a sensible choice for the reference density used in a Boussinesq ocean climate model since, 
    19251937with the exception of only a small percentage of the ocean, density in the World Ocean varies by no more than 
    1926 2$\%$ from this value (\cite{Gill1982}, page 47). 
     19382$\%$ from this value (\cite{gill_bk82}, page 47). 
    19271939 
    19281940Second, we have assumed here that the total ocean surface, $\mathcal{A}$, 
     
    19311943   
    19321944Third, the discretisation of \autoref{eq:steric_Bq} depends on the type of free surface which is considered. 
    1933 In the non linear free surface case, \ie \key{vvl} defined, it is given by 
     1945In the non linear free surface case, \ie \np{ln\_linssh}\forcode{ = .true.}, it is given by 
    19341946 
    19351947\[ 
     
    19541966so that there are no associated ocean currents. 
    19551967Hence, the dynamically relevant sea level is the effective sea level, 
    1956 \ie the sea level as if sea ice (and snow) were converted to liquid seawater \citep{Campin_al_OM08}. 
     1968\ie the sea level as if sea ice (and snow) were converted to liquid seawater \citep{campin.marshall.ea_OM08}. 
    19571969However, in the current version of \NEMO the sea-ice is levitating above the ocean without mass exchanges between 
    19581970ice and ocean. 
     
    19701982where $S_o$ and $p_o$ are the initial salinity and pressure, respectively. 
    19711983 
    1972 Both steric and thermosteric sea level are computed in \mdl{diaar5} which needs the \key{diaar5} defined to 
    1973 be called. 
     1984Both steric and thermosteric sea level are computed in \mdl{diaar5}. 
    19741985 
    19751986% ------------------------------------------------------------------------------------------------------------- 
    19761987%       Other Diagnostics 
    19771988% ------------------------------------------------------------------------------------------------------------- 
    1978 \section{Other diagnostics (\protect\key{diahth}, \protect\key{diaar5})} 
     1989\section[Other diagnostics] 
     1990{Other diagnostics} 
    19791991\label{sec:DIA_diag_others} 
    19801992 
     
    19821994The available ready-to-add diagnostics modules can be found in directory DIA. 
    19831995 
    1984 \subsection{Depth of various quantities (\protect\mdl{diahth})} 
     1996\subsection[Depth of various quantities (\textit{diahth.F90})] 
     1997{Depth of various quantities (\protect\mdl{diahth})} 
    19851998 
    19861999Among the available diagnostics the following ones are obtained when defining the \key{diahth} CPP key: 
    19872000 
    1988 - the mixed layer depth (based on a density criterion \citep{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth}) 
     2001- the mixed layer depth (based on a density criterion \citep{de-boyer-montegut.madec.ea_JGR04}) (\mdl{diahth}) 
    19892002 
    19902003- the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diahth}) 
     
    19942007- the depth of the thermocline (maximum of the vertical temperature gradient) (\mdl{diahth}) 
    19952008 
    1996 % ----------------------------------------------------------- 
    1997 %     Poleward heat and salt transports 
    1998 % ----------------------------------------------------------- 
    1999  
    2000 \subsection{Poleward heat and salt transports (\protect\mdl{diaptr})} 
    2001  
    2002 %------------------------------------------namptr----------------------------------------- 
    2003  
    2004 \nlst{namptr}  
    2005 %----------------------------------------------------------------------------------------- 
    2006  
    2007 The poleward heat and salt transports, their advective and diffusive component, 
    2008 and the meriodional stream function can be computed on-line in \mdl{diaptr} \np{ln\_diaptr} to true 
    2009 (see the \textit{\ngn{namptr} } namelist below). 
    2010 When \np{ln\_subbas}\forcode{ = .true.}, transports and stream function are computed for the Atlantic, Indian, 
    2011 Pacific and Indo-Pacific Oceans (defined north of 30\deg{S}) as well as for the World Ocean. 
    2012 The sub-basin decomposition requires an input file (\ifile{subbasins}) which contains three 2D mask arrays, 
    2013 the Indo-Pacific mask been deduced from the sum of the Indian and Pacific mask (\autoref{fig:mask_subasins}). 
    20142009 
    20152010%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    20162011\begin{figure}[!t] 
    20172012  \begin{center} 
    2018     \includegraphics[width=1.0\textwidth]{Fig_mask_subasins} 
     2013    \includegraphics[width=\textwidth]{Fig_mask_subasins} 
    20192014    \caption{ 
    20202015      \protect\label{fig:mask_subasins} 
     
    20322027%       CMIP specific diagnostics  
    20332028% ----------------------------------------------------------- 
    2034 \subsection{CMIP specific diagnostics (\protect\mdl{diaar5})} 
    2035  
    2036 A series of diagnostics has been added in the \mdl{diaar5}. 
    2037 They corresponds to outputs that are required for AR5 simulations (CMIP5) 
     2029\subsection[CMIP specific diagnostics (\textit{diaar5.F90}, \textit{diaptr.F90})] 
     2030{CMIP specific diagnostics (\protect\mdl{diaar5})} 
     2031 
     2032A series of diagnostics has been added in the \mdl{diaar5} and \mdl{diaptr}. 
     2033In \mdl{diaar5} they correspond to outputs that are required for AR5 simulations (CMIP5) 
    20382034(see also \autoref{sec:DIA_steric} for one of them). 
    2039 Activating those outputs requires to define the \key{diaar5} CPP key. 
     2035The 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). 
     2036 
     2037In \mdl{diaptr} when \np{ln\_diaptr}\forcode{ = .true.}  
     2038(see the \textit{\ngn{namptr} } namelist below) can be computed on-line the poleward heat and salt transports, their advective and diffusive component, and the meriodional stream function . 
     2039When \np{ln\_subbas}\forcode{ = .true.}, transports and stream function are computed for the Atlantic, Indian, 
     2040Pacific and Indo-Pacific Oceans (defined north of 30\deg{S}) as well as for the World Ocean. 
     2041The sub-basin decomposition requires an input file (\ifile{subbasins}) which contains three 2D mask arrays, 
     2042the Indo-Pacific mask been deduced from the sum of the Indian and Pacific mask (\autoref{fig:mask_subasins}). 
     2043 
     2044%------------------------------------------namptr----------------------------------------- 
     2045 
     2046\nlst{namptr}  
     2047%----------------------------------------------------------------------------------------- 
    20402048 
    20412049% ----------------------------------------------------------- 
  • NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_DIU.tex

    r10442 r11422  
    3333(\ie from the temperature of the top few model levels) or from some other source.   
    3434It must be noted that both the cool skin and warm layer models produce estimates of the change in temperature 
    35 ($\Delta T_{\rm{cs}}$ and $\Delta T_{\rm{wl}}$) and 
     35($\Delta T_{\mathrm{cs}}$ and $\Delta T_{\mathrm{wl}}$) and 
    3636both must be added to a foundation SST to obtain the true skin temperature. 
    3737 
     
    6060%=============================================================== 
    6161 
    62 The warm layer is calculated using the model of \citet{Takaya_al_JGR10} (TAKAYA10 model hereafter). 
     62The warm layer is calculated using the model of \citet{takaya.bidlot.ea_JGR10} (TAKAYA10 model hereafter). 
    6363This is a simple flux based model that is defined by the equations 
    6464\begin{align} 
    65 \frac{\partial{\Delta T_{\rm{wl}}}}{\partial{t}}&=&\frac{Q(\nu+1)}{D_T\rho_w c_p 
     65\frac{\partial{\Delta T_{\mathrm{wl}}}}{\partial{t}}&=&\frac{Q(\nu+1)}{D_T\rho_w c_p 
    6666\nu}-\frac{(\nu+1)ku^*_{w}f(L_a)\Delta T}{D_T\Phi\!\left(\frac{D_T}{L}\right)} \mbox{,} 
    6767\label{eq:ecmwf1} \\ 
    6868L&=&\frac{\rho_w c_p u^{*^3}_{w}}{\kappa g \alpha_w Q }\mbox{,}\label{eq:ecmwf2} 
    6969\end{align} 
    70 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. 
     70where $\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. 
    7171In equation (\autoref{eq:ecmwf1}) $\alpha_w=2\times10^{-4}$ is the thermal expansion coefficient of water, 
    7272$\kappa=0.4$ is von K\'{a}rm\'{a}n's constant, $c_p$ is the heat capacity at constant pressure of sea water, 
    7373$\rho_w$ is the water density, and $L$ is the Monin-Obukhov length. 
    7474The tunable variable $\nu$ is a shape parameter that defines the expected subskin temperature profile via 
    75 $T(z) = T(0) - \left( \frac{z}{D_T} \right)^\nu \Delta T_{\rm{wl}}$, 
     75$T(z) = T(0) - \left( \frac{z}{D_T} \right)^\nu \Delta T_{\mathrm{wl}}$, 
    7676where $T$ is the absolute temperature and $z\le D_T$ is the depth below the top of the warm layer. 
    7777The influence of wind on TAKAYA10 comes through the magnitude of the friction velocity of the water $u^*_{w}$, 
     
    8282the diurnal layer, \ie 
    8383\[ 
    84   Q = Q_{\rm{sol}} + Q_{\rm{lw}} + Q_{\rm{h}}\mbox{,} 
     84  Q = Q_{\mathrm{sol}} + Q_{\mathrm{lw}} + Q_{\mathrm{h}}\mbox{,} 
    8585  % \label{eq:e_flux_eqn} 
    8686\] 
    87 where $Q_{\rm{h}}$ is the sensible and latent heat flux, $Q_{\rm{lw}}$ is the long wave flux, 
    88 and $Q_{\rm{sol}}$ is the solar flux absorbed within the diurnal warm layer. 
    89 For $Q_{\rm{sol}}$ the 9 term representation of \citet{Gentemann_al_JGR09} is used. 
     87where $Q_{\mathrm{h}}$ is the sensible and latent heat flux, $Q_{\mathrm{lw}}$ is the long wave flux, 
     88and $Q_{\mathrm{sol}}$ is the solar flux absorbed within the diurnal warm layer. 
     89For $Q_{\mathrm{sol}}$ the 9 term representation of \citet{gentemann.minnett.ea_JGR09} is used. 
    9090In equation \autoref{eq:ecmwf1} the function $f(L_a)=\max(1,L_a^{\frac{2}{3}})$, 
    9191where $L_a=0.3$\footnote{ 
     
    118118%=============================================================== 
    119119 
    120 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}$. 
    121 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 
     120The 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}$. 
     121As 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 
    122122\[ 
    123123  % \label{eq:sunders_eqn} 
    124   \Delta T_{\rm{cs}}=\frac{Q_{\rm{ns}}\delta}{k_t} \mbox{,} 
     124  \Delta T_{\mathrm{cs}}=\frac{Q_{\mathrm{ns}}\delta}{k_t} \mbox{,} 
    125125\] 
    126 where $Q_{\rm{ns}}$ is the, usually negative, non-solar heat flux into the ocean and 
     126where $Q_{\mathrm{ns}}$ is the, usually negative, non-solar heat flux into the ocean and 
    127127$k_t$ is the thermal conductivity of sea water. 
    128128$\delta$ is the thickness of the skin layer and is given by 
     
    132132\end{equation} 
    133133where $\mu$ is the kinematic viscosity of sea water and $\lambda$ is a constant of proportionality which 
    134 \citet{Saunders_JAS82} suggested varied between 5 and 10. 
     134\citet{saunders_JAS67} suggested varied between 5 and 10. 
    135135 
    136 The value of $\lambda$ used in equation (\autoref{eq:sunders_thick_eqn}) is that of \citet{Artale_al_JGR02}, 
    137 which is shown in \citet{Tu_Tsuang_GRL05} to outperform a number of other parametrisations at 
     136The value of $\lambda$ used in equation (\autoref{eq:sunders_thick_eqn}) is that of \citet{artale.iudicone.ea_JGR02}, 
     137which is shown in \citet{tu.tsuang_GRL05} to outperform a number of other parametrisations at 
    138138both low and high wind speeds. 
    139139Specifically, 
  • NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_DOM.tex

    r10502 r11422  
    1818%     - domclo:  closed sea and lakes.... management of closea sea area : specific to global configuration, both forced and coupled 
    1919 
     20\vfill 
     21\begin{figure}[b] 
     22\subsubsection*{Changes record} 
     23\begin{tabular}{m{0.08\linewidth}||m{0.32\linewidth}|m{0.6\linewidth}} 
     24    Release   & Author(s)     & Modifications \\ 
     25\hline 
     26    {\em 4.0} & {\em Simon M{\"u}ller \& Andrew Coward} & {\em Compatibility changes for v4.0. Major simplication has moved many of the options to external domain configuration tools. For now this information has been retained in \autoref{apdx:DOMAINcfg} }  \\ 
     27    {\em 3.x} & {\em Sebastien Masson, Gurvan Madec \& Rashid Benshila } & {\em }  \\ 
     28\end{tabular} 
     29\end{figure} 
     30 
    2031\newpage 
    2132 
    2233Having defined the continuous equations in \autoref{chap:PE} and chosen a time discretization \autoref{chap:STP}, 
    23 we need to choose a discretization on a grid, and numerical algorithms. 
     34we need to choose a grid for spatial discretization and related numerical algorithms. 
    2435In the present chapter, we provide a general description of the staggered grid used in \NEMO, 
    25 and other information relevant to the main directory routines as well as the DOM (DOMain) directory. 
     36and other relevant information about the DOM (DOMain) source-code modules . 
    2637 
    2738% ================================================================ 
     
    4051\begin{figure}[!tb] 
    4152  \begin{center} 
    42     \includegraphics[]{Fig_cell} 
     53    \includegraphics[width=\textwidth]{Fig_cell} 
    4354    \caption{ 
    4455      \protect\label{fig:cell} 
     
    5566The numerical techniques used to solve the Primitive Equations in this model are based on the traditional, 
    5667centred second-order finite difference approximation. 
    57 Special attention has been given to the homogeneity of the solution in the three space directions. 
     68Special attention has been given to the homogeneity of the solution in the three spatial directions. 
    5869The arrangement of variables is the same in all directions. 
    5970It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in 
    6071the centre of each face of the cells (\autoref{fig:cell}). 
    6172This is the generalisation to three dimensions of the well-known ``C'' grid in Arakawa's classification 
    62 \citep{Mesinger_Arakawa_Bk76}. 
     73\citep{mesinger.arakawa_bk76}. 
    6374The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and 
    6475the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points. 
     
    7182Each scale factor is defined as the local analytical value provided by \autoref{eq:scale_factors}. 
    7283As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and 
    73 $\pd[]{z}$ are evaluated in a uniform mesh with a grid size of unity. 
     84$\pd[]{z}$ are evaluated is a uniform mesh with a grid size of unity. 
    7485Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation 
    7586while the scale factors are chosen equal to their local analytical value. 
     
    7990the continuous properties (see \autoref{apdx:C}). 
    8091A similar, related remark can be made about the domain size: 
    81 when needed, an area, volume, or the total ocean depth must be evaluated as the sum of the relevant scale factors 
     92when needed, an area, volume, or the total ocean depth must be evaluated as the product or sum of the relevant scale factors 
    8293(see \autoref{eq:DOM_bar} in the next section). 
    8394 
     
    8798    \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|} 
    8899      \hline 
    89       T  & $i      $ & $j      $ & $k      $ \\ 
     100      t  & $i      $ & $j      $ & $k      $ \\ 
    90101      \hline 
    91102      u  & $i + 1/2$ & $j      $ & $k      $ \\ 
     
    107118      \protect\label{tab:cell} 
    108119      Location of grid-points as a function of integer or integer and a half value of the column, line or level. 
    109       This indexing is only used for the writing of the semi -discrete equation. 
    110       In the code, the indexing uses integer values only and has a reverse direction in the vertical 
     120      This indexing is only used for the writing of the semi -discrete equations. 
     121      In the code, the indexing uses integer values only and is positive downwards in the vertical with $k=1$ at the surface. 
    111122      (see \autoref{subsec:DOM_Num_Index}) 
    112123    } 
    113124  \end{center} 
    114125\end{table} 
     126%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     127 
     128Note that the definition of the scale factors 
     129(\ie as the analytical first derivative of the transformation that 
     130results in $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) 
     131is specific to the \NEMO model \citep{marti.madec.ea_JGR92}. 
     132As an example, a scale factor in the $i$ direction is defined locally at a $t$-point, 
     133whereas many other models on a C grid choose to define such a scale factor as 
     134the distance between the $u$-points on each side of the $t$-point. 
     135Relying on an analytical transformation has two advantages: 
     136firstly, there is no ambiguity in the scale factors appearing in the discrete equations, 
     137since they are first introduced in the continuous equations; 
     138secondly, analytical transformations encourage good practice by the definition of smoothly varying grids 
     139(rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{treguier.dukowicz.ea_JGR96}. 
     140An example of the effect of such a choice is shown in \autoref{fig:zgr_e3}. 
     141%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     142\begin{figure}[!t] 
     143  \begin{center} 
     144    \includegraphics[width=\textwidth]{Fig_zgr_e3} 
     145    \caption{ 
     146      \protect\label{fig:zgr_e3} 
     147      Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, 
     148      and (b) analytically derived grid-point position and scale factors. 
     149      For both grids here, the same $w$-point depth has been chosen but 
     150      in (a) the $t$-points are set half way between $w$-points while 
     151      in (b) they are defined from an analytical function: 
     152      $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$. 
     153      Note the resulting difference between the value of the grid-size $\Delta_k$ and 
     154      those of the scale factor $e_k$. 
     155    } 
     156  \end{center} 
     157\end{figure} 
    115158%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    116159 
     
    132175Following \autoref{eq:PE_grad} and \autoref{eq:PE_lap}, the gradient of a variable $q$ defined at 
    133176a $t$-point has its three components defined at $u$-, $v$- and $w$-points while 
    134 its Laplacian is defined at $t$-point. 
     177its Laplacian is defined at the $t$-point. 
    135178These operators have the following discrete forms in the curvilinear $s$-coordinates system: 
    136179\[ 
     
    171214\end{equation} 
    172215 
    173 The vertical average over the whole water column denoted by an overbar becomes for a quantity $q$ which 
    174 is a masked field (i.e. equal to zero inside solid area): 
     216The vertical average over the whole water column is denoted by an overbar and is for 
     217a masked field $q$ (\ie a quantity that is equal to zero inside solid areas): 
    175218\begin{equation} 
    176219  \label{eq:DOM_bar} 
     
    178221\end{equation} 
    179222where $H_q$  is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points, 
    180 $k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $k^o$ refers to a summation over 
     223$k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $\sum \limits_k$ refers to a summation over 
    181224all grid points of the same type in the direction indicated by the subscript (here $k$). 
    182225 
     
    193236vector points $(u,v,w)$. 
    194237 
    195 Let $a$ and $b$ be two fields defined on the mesh, with value zero inside continental area. 
    196 Using integration by parts it can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) 
     238Let $a$ and $b$ be two fields defined on the mesh, with a value of zero inside continental areas. 
     239It can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) 
    197240are skew-symmetric linear operators, and further that the averaging operators $\overline{\cdots}^{\, i}$, 
    198241$\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie 
     
    218261\begin{figure}[!tb] 
    219262  \begin{center} 
    220     \includegraphics[]{Fig_index_hor} 
     263    \includegraphics[width=\textwidth]{Fig_index_hor} 
    221264    \caption{ 
    222265      \protect\label{fig:index_hor} 
     
    228271%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    229272 
    230 The array representation used in the \fortran code requires an integer indexing while 
    231 the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of 
    232 integer values for $t$-points and both integer and integer and a half values for all the other points. 
    233 Therefore a specific integer indexing must be defined for points other than $t$-points 
     273The array representation used in the \fortran code requires an integer indexing. 
     274However, the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of 
     275integer values for $t$-points only while all the other points involve integer and a half values. 
     276Therefore, a specific integer indexing has been defined for points other than $t$-points 
    234277(\ie velocity and vorticity grid-points). 
    235 Furthermore, the direction of the vertical indexing has been changed so that the surface level is at $k = 1$. 
     278Furthermore, the direction of the vertical indexing has been reversed and the surface level set at $k = 1$. 
    236279 
    237280% ----------------------------------- 
     
    253296\label{subsec:DOM_Num_Index_vertical} 
    254297 
    255 In the vertical, the chosen indexing requires special attention since the $k$-axis is re-orientated downward in 
    256 the \fortran code compared to the indexing used in the semi -discrete equations and 
     298In the vertical, the chosen indexing requires special attention since the direction of the $k$-axis in 
     299the \fortran code is the reverse of that used in the semi -discrete equations and 
    257300given in \autoref{subsec:DOM_cell}. 
    258 The sea surface corresponds to the $w$-level $k = 1$ which is the same index as $t$-level just below 
     301The sea surface corresponds to the $w$-level $k = 1$, which is the same index as the $t$-level just below 
    259302(\autoref{fig:index_vert}). 
    260 The last $w$-level ($k = jpk$) either corresponds to the ocean floor or is inside the bathymetry while 
    261 the last $t$-level is always inside the bathymetry (\autoref{fig:index_vert}). 
    262 Note that for an increasing $k$ index, a $w$-point and the $t$-point just below have the same $k$ index, 
    263 in opposition to what is done in the horizontal plane where 
    264 it is the $t$-point and the nearest velocity points in the direction of the horizontal axis that 
    265 have the same $i$ or $j$ index 
     303The last $w$-level ($k = jpk$) either corresponds to or is below the ocean floor while 
     304the last $t$-level is always outside the ocean domain (\autoref{fig:index_vert}). 
     305Note that a $w$-point and the directly underlaying $t$-point have a common $k$ index (\ie $t$-points and their 
     306nearest $w$-point neighbour in negative index direction), in contrast to the indexing on the horizontal plane where 
     307the $t$-point has the same index as the nearest velocity points in the positive direction of the respective horizontal axis index 
    266308(compare the dashed area in \autoref{fig:index_hor} and \autoref{fig:index_vert}). 
    267309Since the scale factors are chosen to be strictly positive, 
    268 a \textit{minus sign} appears in the \fortran code \textit{before all the vertical derivatives} of 
    269 the discrete equations given in this documentation. 
     310a \textit{minus sign} is included in the \fortran implementations of \textit{all the vertical derivatives} of 
     311the discrete equations given in this manual in order to accommodate the opposing vertical index directions in implementation and documentation. 
    270312 
    271313%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    272314\begin{figure}[!pt] 
    273315  \begin{center} 
    274     \includegraphics[]{Fig_index_vert} 
     316    \includegraphics[width=\textwidth]{Fig_index_vert} 
    275317    \caption{ 
    276318      \protect\label{fig:index_vert} 
    277319      Vertical integer indexing used in the \fortran code. 
    278       Note that the $k$-axis is orientated downward. 
    279       The dashed area indicates the cell in which variables contained in arrays have the same $k$-index. 
     320      Note that the $k$-axis is oriented downward. 
     321      The dashed area indicates the cell in which variables contained in arrays have a common $k$-index. 
    280322    } 
    281323  \end{center} 
     
    283325%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    284326 
     327% ------------------------------------------------------------------------------------------------------------- 
     328%        Domain configuration 
     329% ------------------------------------------------------------------------------------------------------------- 
     330\section{Spatial domain configuration} 
     331\label{subsec:DOM_config} 
     332 
     333\nlst{namcfg} 
     334 
     335Two typical methods are available to specify the spatial domain 
     336configuration; they can be selected using parameter \np{ln\_read\_cfg} 
     337parameter in namelist \ngn{namcfg}.  
     338 
     339If \np{ln\_read\_cfg} is set to \forcode{.true.}, the domain-specific parameters 
     340and fields are read from a netCDF input file, whose name (without its .nc 
     341suffix) can be specified as the value of the \np{cn\_domcfg} parameter in 
     342namelist \ngn{namcfg}. 
     343 
     344If \np{ln\_read\_cfg} is set to \forcode{.false.}, the domain-specific 
     345parameters and fields can be provided (\eg analytically computed) by subroutines 
     346\mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. These subroutines can be supplied in 
     347the \path{MY_SRC} directory of the configuration, and default versions that 
     348configure the spatial domain for the GYRE reference configuration are present in 
     349the \path{src/OCE/USR} directory. 
     350 
     351In version 4.0 there are no longer any options for reading complex bathmetries and  
     352performing a vertical discretization at run-time. Whilst it is occasionally convenient 
     353to have a common bathymetry file and, for example, to run similar models with and 
     354without partial bottom boxes and/or sigma-coordinates, supporting such choices leads to 
     355overly complex code. Worse still is the difficulty of ensuring the model configurations  
     356intended to be identical are indeed so when the model domain itself can be altered by runtime 
     357selections. The code previously used to perform vertical discretization has be incorporated  
     358into an external tool (\path{tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMAINcfg}. 
     359 
     360The next subsections summarise the parameter and fields related to the 
     361configuration of the whole model domain. These represent the minimum information 
     362that must be provided either via the \np{cn\_domcfg} file or set by code 
     363inserted into user-supplied versions of the \mdl{usrdef\_*} subroutines. The 
     364requirements are presented in three sections: the domain size 
     365(\autoref{subsec:DOM_size}), the horizontal mesh 
     366(\autoref{subsec:DOM_hgr}), and the vertical grid 
     367(\autoref{subsec:DOM_zgr}). 
     368 
    285369% ----------------------------------- 
    286370%        Domain Size 
    287371% ----------------------------------- 
    288 \subsubsection{Domain size} 
     372\subsection{Domain size} 
    289373\label{subsec:DOM_size} 
    290374 
    291 The total size of the computational domain is set by the parameters \np{jpiglo}, 
    292 \np{jpjglo} and \np{jpkglo} in the $i$, $j$ and $k$ directions respectively. 
    293 Parameters $jpi$ and $jpj$ refer to the size of each processor subdomain when 
    294 the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined, 
    295 see \autoref{sec:LBC_mpp}). 
    296  
    297 % ================================================================ 
    298 % Domain: List of fields needed 
    299 % ================================================================ 
    300 \section{Needed fields} 
    301 \label{sec:DOM_fields} 
    302 The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that 
    303 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 
    304 The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}. 
    305 The associated scale factors are defined using the analytical first derivative of the transformation 
    306 \autoref{eq:scale_factors}. 
    307 Necessary fields for configuration definition are: 
    308  
    309 \begin{itemize} 
    310 \item 
    311   Geographic position: 
    312   longitude with \texttt{glamt}, \texttt{glamu}, \texttt{glamv}, \texttt{glamf} and 
    313   latitude  with \texttt{gphit}, \texttt{gphiu}, \texttt{gphiv}, \texttt{gphif} 
    314   (all respectively at T, U, V and F point) 
    315 \item 
    316   Coriolis parameter (if domain not on the sphere): \texttt{ff\_f} and \texttt{ff\_t} 
    317   (at T and F point) 
    318 \item 
    319   Scale factors: 
    320   \texttt{e1t}, \texttt{e1u}, \texttt{e1v} and \texttt{e1f} (on i direction), 
    321   \texttt{e2t}, \texttt{e2u}, \texttt{e2v} and \texttt{e2f} (on j direction) and 
    322   \texttt{ie1e2u\_v}, \texttt{e1e2u}, \texttt{e1e2v}. \\ 
    323   \texttt{e1e2u}, \texttt{e1e2v} are u and v surfaces (if gridsize reduction in some straits),  
    324   \texttt{ie1e2u\_v} is to flag set u and v surfaces are neither read nor computed. 
    325 \end{itemize} 
    326   
    327 These fields can be read in an domain input file which name is setted in \np{cn\_domcfg} parameter specified in 
    328 \ngn{namcfg}. 
    329  
    330 \nlst{namcfg} 
    331  
    332 Or they can be defined in an analytical way in \path{MY_SRC} directory of the configuration. 
    333 For Reference Configurations of NEMO input domain files are supplied by NEMO System Team. 
    334 For analytical definition of input fields two routines are supplied: \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. 
    335 They are an example of GYRE configuration parameters, and they are available in \path{src/OCE/USR} directory, 
    336 they provide the horizontal and vertical mesh. 
    337 % ------------------------------------------------------------------------------------------------------------- 
    338 %        Needed fields  
    339 % ------------------------------------------------------------------------------------------------------------- 
    340 %\subsection{List of needed fields to build DOMAIN} 
    341 %\label{subsec:DOM_fields_list} 
    342  
     375The total size of the computational domain is set by the parameters 
     376\np{jpiglo}, \np{jpjglo} and \np{jpkglo} for the $i$, $j$ and $k$ 
     377directions, respectively. Note, that the variables \forcode{jpi} and \forcode{jpj} 
     378refer to the size of each processor subdomain when the code is run in 
     379parallel using domain decomposition (\key{mpp\_mpi} defined, see 
     380\autoref{sec:LBC_mpp}). 
     381 
     382The name of the configuration is set through parameter \np{cn\_cfg}, 
     383and the nominal resolution through parameter \np{nn\_cfg} (unless in 
     384the input file both of variables \forcode{ORCA} and \forcode{ORCA_index} 
     385are present, in which case \np{cn\_cfg} and \np{nn\_cfg} are set from these 
     386values accordingly). 
     387 
     388The global lateral boundary condition type is selected from 8 options 
     389using parameter \np{jperio}. See \autoref{sec:LBC_jperio} for 
     390details on the available options and the corresponding values for 
     391\np{jperio}. 
    343392 
    344393% ================================================================ 
    345394% Domain: Horizontal Grid (mesh)  
    346395% ================================================================ 
    347 \section{Horizontal grid mesh (\protect\mdl{domhgr})} 
    348 \label{sec:DOM_hgr} 
    349  
    350 % ------------------------------------------------------------------------------------------------------------- 
    351 %        Coordinates and scale factors  
    352 % ------------------------------------------------------------------------------------------------------------- 
    353 \subsection{Coordinates and scale factors} 
    354 \label{subsec:DOM_hgr_coord_e} 
    355  
    356 The ocean mesh (\ie the position of all the scalar and vector points) is defined by 
    357 the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 
    358 The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}. 
    359 The associated scale factors are defined using the analytical first derivative of the transformation 
    360 \autoref{eq:scale_factors}. 
    361 These definitions are done in two modules, \mdl{domhgr} and \mdl{domzgr}, 
    362 which provide the horizontal and vertical meshes, respectively. 
    363 This section deals with the horizontal mesh parameters. 
    364  
    365 In a horizontal plane, the location of all the model grid points is defined from 
    366 the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$. 
    367 The horizontal scale factors are calculated using \autoref{eq:scale_factors}. 
    368 For example, when the longitude and latitude are function of a single value 
    369 ($i$ and $j$, respectively) (geographical configuration of the mesh), 
    370 the horizontal mesh definition reduces to define the wanted $\lambda(i)$, $\varphi(j)$, 
    371 and their derivatives $\lambda'(i) \ \varphi'(j)$ in the \mdl{domhgr} module. 
    372 The model computes the grid-point positions and scale factors in the horizontal plane as follows: 
    373 \begin{align*} 
    374    \lambda_t &\equiv \text{glamt} =      \lambda (i      ) 
    375   &\varphi_t &\equiv \text{gphit} =      \varphi (j      ) \\ 
    376    \lambda_u &\equiv \text{glamu} =      \lambda (i + 1/2) 
    377   &\varphi_u &\equiv \text{gphiu} =      \varphi (j      ) \\ 
    378    \lambda_v &\equiv \text{glamv} =      \lambda (i      ) 
    379   &\varphi_v &\equiv \text{gphiv} =      \varphi (j + 1/2) \\ 
    380    \lambda_f &\equiv \text{glamf} =      \lambda (i + 1/2) 
    381   &\varphi_f &\equiv \text{gphif} =      \varphi (j + 1/2) \\ 
    382    e_{1t}    &\equiv \text{e1t}   = r_a |\lambda'(i      ) \; \cos\varphi(j      ) | 
    383   &e_{2t}    &\equiv \text{e2t}   = r_a |\varphi'(j      )                         | \\ 
    384    e_{1u}    &\equiv \text{e1t}   = r_a |\lambda'(i + 1/2) \; \cos\varphi(j      ) | 
    385   &e_{2u}    &\equiv \text{e2t}   = r_a |\varphi'(j      )                         | \\ 
    386    e_{1v}    &\equiv \text{e1t}   = r_a |\lambda'(i      ) \; \cos\varphi(j + 1/2) | 
    387   &e_{2v}    &\equiv \text{e2t}   = r_a |\varphi'(j + 1/2)                         | \\ 
    388    e_{1f}    &\equiv \text{e1t}   = r_a |\lambda'(i + 1/2) \; \cos\varphi(j + 1/2) | 
    389   &e_{2f}    &\equiv \text{e2t}   = r_a |\varphi'(j + 1/2)                         | 
    390 \end{align*} 
    391 where the last letter of each computational name indicates the grid point considered and 
    392 $r_a$ is the earth radius (defined in \mdl{phycst} along with all universal constants). 
    393 Note that the horizontal position of and scale factors at $w$-points are exactly equal to those of $t$-points, 
    394 thus no specific arrays are defined at $w$-points. 
    395  
    396 Note that the definition of the scale factors 
    397 (\ie as the analytical first derivative of the transformation that 
    398 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) 
    399 is specific to the \NEMO model \citep{Marti_al_JGR92}. 
    400 As an example, $e_{1t}$ is defined locally at a $t$-point, 
    401 whereas many other models on a C grid choose to define such a scale factor as 
    402 the distance between the $U$-points on each side of the $t$-point. 
    403 Relying on an analytical transformation has two advantages: 
    404 firstly, there is no ambiguity in the scale factors appearing in the discrete equations, 
    405 since they are first introduced in the continuous equations; 
    406 secondly, analytical transformations encourage good practice by the definition of smoothly varying grids 
    407 (rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{Treguier1996}. 
    408 An example of the effect of such a choice is shown in \autoref{fig:zgr_e3}. 
    409 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    410 \begin{figure}[!t] 
    411   \begin{center} 
    412     \includegraphics[]{Fig_zgr_e3} 
    413     \caption{ 
    414       \protect\label{fig:zgr_e3} 
    415       Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, 
    416       and (b) analytically derived grid-point position and scale factors. 
    417       For both grids here, the same $w$-point depth has been chosen but 
    418       in (a) the $t$-points are set half way between $w$-points while 
    419       in (b) they are defined from an analytical function: 
    420       $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$. 
    421       Note the resulting difference between the value of the grid-size $\Delta_k$ and 
    422       those of the scale factor $e_k$. 
    423     } 
    424   \end{center} 
    425 \end{figure} 
    426 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    427  
    428 % ------------------------------------------------------------------------------------------------------------- 
    429 %        Choice of horizontal grid 
    430 % ------------------------------------------------------------------------------------------------------------- 
    431 \subsection{Choice of horizontal grid} 
    432 \label{subsec:DOM_hgr_msh_choice} 
    433  
    434 % ------------------------------------------------------------------------------------------------------------- 
    435 %        Grid files 
    436 % ------------------------------------------------------------------------------------------------------------- 
    437 \subsection{Output grid files} 
    438 \label{subsec:DOM_hgr_files} 
    439  
    440 All the arrays relating to a particular ocean model configuration (grid-point position, scale factors, masks) 
    441 can be saved in files if \np{nn\_msh} $\not = 0$ (namelist variable in \ngn{namdom}). 
    442 This can be particularly useful for plots and off-line diagnostics. 
    443 In some cases, the user may choose to make a local modification of a scale factor in the code. 
    444 This is the case in global configurations when restricting the width of a specific strait 
    445 (usually a one-grid-point strait that happens to be too wide due to insufficient model resolution). 
    446 An example is Gibraltar Strait in the ORCA2 configuration. 
    447 When such modifications are done, 
    448 the output grid written when \np{nn\_msh} $\not = 0$ is no more equal to the input grid. 
     396\subsection{Horizontal grid mesh (\protect\mdl{domhgr})} 
     397\label{subsec:DOM_hgr} 
     398 
     399% ================================================================ 
     400% Domain: List of hgr-related fields needed 
     401% ================================================================ 
     402\subsubsection{Required fields} 
     403\label{sec:DOM_hgr_fields} 
     404The explicit specification of a range of mesh-related fields are required for the definition of a configuration. These include: 
     405 
     406\begin{Verbatim}[fontsize=\tiny] 
     407int    jpiglo, jpjglo, jpkglo            /* global domain sizes                                          */ 
     408int    jperio                            /* lateral global domain b.c.                                   */ 
     409double glamt, glamu, glamv, glamf        /* geographic longitude (t,u,v and f points respectively)       */ 
     410double gphit, gphiu, gphiv, gphif        /* geographic latitude                                          */ 
     411double e1t, e1u, e1v, e1f                /* horizontal scale factors                                     */ 
     412double e2t, e2u, e2v, e2f                /* horizontal scale factors                                     */ 
     413\end{Verbatim} 
     414 
     415The 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 Table \autoref{tab: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$. 
     416 
     417\subsubsection{Optional fields} 
     418\begin{Verbatim}[fontsize=\tiny] 
     419                                         /* Optional:                                                    */ 
     420int    ORCA, ORCA_index                  /* configuration name, configuration resolution                 */ 
     421double e1e2u, e1e2v                      /* U and V surfaces (if grid size reduction in some straits)    */ 
     422double ff_f, ff_t                        /* Coriolis parameter (if not on the sphere)                    */ 
     423\end{Verbatim} 
     424 
     425NEMO can support the local reduction of key strait widths by altering individual values of 
     426e2u or e1v at the appropriate locations. This is particularly useful for locations such as 
     427Gibraltar or Indonesian Throughflow pinch-points (see \autoref{sec:MISC_strait} for 
     428illustrated examples). The key is to reduce the faces of $T$-cell (\ie change the value of 
     429the horizontal scale factors at $u$- or $v$-point) but not the volume of the cells. Doing 
     430otherwise can lead to numerical instability issues.  In normal operation the surface areas 
     431are computed from $\texttt{e1u} * \texttt{e2u}$ and $\texttt{e1v} * \texttt{e2v}$ but in 
     432cases where a gridsize reduction is required, the unaltered surface areas at $u$ and $v$ 
     433grid points (\texttt{e1e2u} and \texttt{e1e2v}, respectively) must be read or pre-computed 
     434in \mdl{usrdef\_hgr}. If these arrays are present in the \np{cn\_domcfg} file they are 
     435read and the internal computation is suppressed. Versions of \mdl{usrdef\_hgr} which set 
     436their own values of \texttt{e1e2u} and \texttt{e1e2v} should set the surface-area 
     437computation flag: \texttt{ie1e2u\_v} to a non-zero value to suppress their re-computation. 
     438 
     439\smallskip 
     440Similar logic applies to the other optional fields: \texttt{ff\_f} and \texttt{ff\_t} 
     441which can be used to provide the Coriolis parameter at F- and T-points respectively if the 
     442mesh is not on a sphere. If present these fields will be read and used and the normal 
     443calculation ($2*\Omega*\sin(\varphi)$) suppressed. Versions of \mdl{usrdef\_hgr} which set 
     444their own values of \texttt{ff\_f} and \texttt{ff\_t} should set the Coriolis computation 
     445flag: \texttt{iff} to a non-zero value to suppress their re-computation. 
     446 
     447Note that longitudes, latitudes, and scale factors at $w$ points are exactly 
     448equal to those of $t$ points, thus no specific arrays are defined at $w$ points. 
     449 
    449450 
    450451% ================================================================ 
    451452% Domain: Vertical Grid (domzgr) 
    452453% ================================================================ 
    453 \section{Vertical grid (\protect\mdl{domzgr})} 
    454 \label{sec:DOM_zgr} 
    455 %-----------------------------------------nam_zgr & namdom------------------------------------------- 
    456 % 
    457 %\nlst{namzgr}  
    458  
    459 \nlst{namdom}  
     454\subsection[Vertical grid (\textit{domzgr.F90})] 
     455{Vertical grid (\protect\mdl{domzgr})} 
     456\label{subsec:DOM_zgr} 
     457%-----------------------------------------namdom------------------------------------------- 
     458\nlst{namdom} 
    460459%------------------------------------------------------------------------------------------------------------- 
    461460 
    462 Variables are defined through the \ngn{namzgr} and \ngn{namdom} namelists. 
    463461In the vertical, the model mesh is determined by four things:  
    464 (1) the bathymetry given in meters;  
    465 (2) the number of levels of the model (\jp{jpk});  
    466 (3) the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and 
    467 (4) the masking system, \ie the number of wet model levels at each  
    468 $(i,j)$ column of points. 
     462\begin{enumerate} 
     463  \item the bathymetry given in meters;  
     464  \item the number of levels of the model (\jp{jpk});  
     465  \item the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and 
     466  \item the masking system, \ie the number of wet model levels at each  
     467$(i,j)$ location of the horizontal grid. 
     468\end{enumerate} 
    469469 
    470470%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    471471\begin{figure}[!tb] 
    472472  \begin{center} 
    473     \includegraphics[]{Fig_z_zps_s_sps} 
     473    \includegraphics[width=\textwidth]{Fig_z_zps_s_sps} 
    474474    \caption{ 
    475475      \protect\label{fig:z_zps_s_sps} 
     
    480480      (d) hybrid $s-z$ coordinate, 
    481481      (e) hybrid $s-z$ coordinate with partial step, and 
    482       (f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh}~\forcode{= .false.}). 
     482      (f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh}\forcode{ = .false.}). 
    483483      Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e). 
    484484    } 
     
    487487%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    488488 
    489 The choice of a vertical coordinate, even if it is made through \ngn{namzgr} namelist parameters,  
    490 must be done once of all at the beginning of an experiment. 
    491 It is not intended as an option which can be enabled or disabled in the middle of an experiment. 
    492 Three main choices are offered (\autoref{fig:z_zps_s_sps}): 
    493 $z$-coordinate with full step bathymetry (\np{ln\_zco}~\forcode{= .true.}), 
    494 $z$-coordinate with partial step bathymetry (\np{ln\_zps}~\forcode{= .true.}), 
    495 or generalized, $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}). 
    496 Hybridation of the three main coordinates are available: 
    497 $s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps} and \autoref{fig:z_zps_s_sps}). 
    498 By default a non-linear free surface is used: the coordinate follow the time-variation of the free surface so that 
    499 the transformation is time dependent: $z(i,j,k,t)$ (\autoref{fig:z_zps_s_sps}). 
    500 When a linear free surface is assumed (\np{ln\_linssh}~\forcode{= .true.}), 
    501 the vertical coordinate are fixed in time, but the seawater can move up and down across the $z_0$ surface 
    502 (in other words, the top of the ocean in not a rigid-lid). 
    503 The last choice in terms of vertical coordinate concerns the presence (or not) in 
    504 the model domain of ocean cavities beneath ice shelves. 
    505 Setting \np{ln\_isfcav} to true allows to manage ocean cavities, otherwise they are filled in. 
    506 This option is currently only available in $z$- or $zps$-coordinate, 
    507 and partial step are also applied at the ocean/ice shelf interface. 
    508  
    509 Contrary to the horizontal grid, the vertical grid is computed in the code and no provision is made for 
    510 reading it from a file. 
    511 The only input file is the bathymetry (in meters) (\ifile{bathy\_meter}) 
    512 \footnote{ 
    513   N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the \ifile{bathy\_meter} file, 
    514   so that the computation of the number of wet ocean point in each water column is by-passed}. 
    515 If \np{ln\_isfcav}~\forcode{= .true.}, an extra file input file (\ifile{isf\_draft\_meter}) describing 
    516 the ice shelf draft (in meters) is needed. 
    517  
    518 After reading the bathymetry, the algorithm for vertical grid definition differs between the different options: 
    519 \begin{description} 
    520 \item[\textit{zco}] 
    521   set a reference coordinate transformation $z_0(k)$, and set $z(i,j,k,t) = z_0(k)$. 
    522 \item[\textit{zps}] 
    523   set a reference coordinate transformation $z_0(k)$, and calculate the thickness of the deepest level at 
    524   each $(i,j)$ point using the bathymetry, to obtain the final three-dimensional depth and scale factor arrays. 
    525 \item[\textit{sco}] 
    526   smooth the bathymetry to fulfill the hydrostatic consistency criteria and 
    527   set the three-dimensional transformation. 
    528 \item[\textit{s-z} and \textit{s-zps}] 
    529   smooth the bathymetry to fulfill the hydrostatic consistency criteria and 
    530   set the three-dimensional transformation $z(i,j,k)$, 
    531   and possibly introduce masking of extra land points to better fit the original bathymetry file. 
    532 \end{description} 
    533 %%% 
    534 \gmcomment{   add the description of the smoothing:  envelop topography...} 
    535 %%% 
    536  
    537 Unless a linear free surface is used (\np{ln\_linssh}~\forcode{= .false.}), 
    538 the arrays describing the grid point depths and vertical scale factors are three set of 
    539 three dimensional arrays $(i,j,k)$ defined at \textit{before}, \textit{now} and \textit{after} time step. 
    540 The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively. 
    541 They are updated at each model time step using a fixed reference coordinate system which 
    542 computer names have a $\_0$ suffix. 
    543 When the linear free surface option is used (\np{ln\_linssh}~\forcode{= .true.}), \textit{before}, 
    544 \textit{now} and \textit{after} arrays are simply set one for all to their reference counterpart. 
    545  
    546 % ------------------------------------------------------------------------------------------------------------- 
    547 %        Meter Bathymetry 
    548 % ------------------------------------------------------------------------------------------------------------- 
    549 \subsection{Meter bathymetry} 
    550 \label{subsec:DOM_bathy} 
    551  
    552 Three options are possible for defining the bathymetry, according to the namelist variable \np{nn\_bathy} 
    553 (found in \ngn{namdom} namelist):  
    554 \begin{description} 
    555 \item[\np{nn\_bathy}~\forcode{= 0}]: 
    556   a flat-bottom domain is defined. 
    557   The total depth $z_w (jpk)$ is given by the coordinate transformation. 
    558   The domain can either be a closed basin or a periodic channel depending on the parameter \np{jperio}. 
    559 \item[\np{nn\_bathy}~\forcode{= -1}]: 
    560   a domain with a bump of topography one third of the domain width at the central latitude. 
    561   This is meant for the "EEL-R5" configuration, a periodic or open boundary channel with a seamount. 
    562 \item[\np{nn\_bathy}~\forcode{= 1}]: 
    563   read a bathymetry and ice shelf draft (if needed). 
    564   The \ifile{bathy\_meter} file (Netcdf format) provides the ocean depth (positive, in meters) at 
    565   each grid point of the model grid. 
    566   The bathymetry is usually built by interpolating a standard bathymetry product (\eg ETOPO2) onto 
    567   the horizontal ocean mesh. 
    568   Defining the bathymetry also defines the coastline: where the bathymetry is zero, 
    569   no model levels are defined (all levels are masked). 
    570  
    571   The \ifile{isfdraft\_meter} file (Netcdf format) provides the ice shelf draft (positive, in meters) at 
    572   each grid point of the model grid. 
    573   This file is only needed if \np{ln\_isfcav}~\forcode{= .true.}. 
    574   Defining the ice shelf draft will also define the ice shelf edge and the grounding line position. 
    575 \end{description} 
    576  
    577 When a global ocean is coupled to an atmospheric model it is better to represent all large water bodies 
    578 (\eg great lakes, Caspian sea...) even if the model resolution does not allow their communication with 
    579 the rest of the ocean. 
    580 This is unnecessary when the ocean is forced by fixed atmospheric conditions, 
    581 so these seas can be removed from the ocean domain. 
    582 The user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}), 
    583 but the code has to be adapted to the user's configuration. 
    584  
    585 % ------------------------------------------------------------------------------------------------------------- 
    586 %        z-coordinate  and reference coordinate transformation 
    587 % ------------------------------------------------------------------------------------------------------------- 
    588 \subsection[$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and ref. coordinate] 
    589             {$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and reference coordinate} 
    590 \label{subsec:DOM_zco} 
    591  
    592 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    593 \begin{figure}[!tb] 
    594   \begin{center} 
    595     \includegraphics[]{Fig_zgr} 
    596     \caption{ 
    597       \protect\label{fig:zgr} 
    598       Default vertical mesh for ORCA2: 30 ocean levels (L30). 
    599       Vertical level functions for (a) T-point depth and (b) the associated scale factor as computed from 
    600       \autoref{eq:DOM_zgr_ana_1} using \autoref{eq:DOM_zgr_coef} in $z$-coordinate. 
    601     } 
    602   \end{center} 
    603 \end{figure} 
    604 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    605  
    606 The reference coordinate transformation $z_0(k)$ defines the arrays $gdept_0$ and $gdepw_0$ for $t$- and $w$-points, 
    607 respectively. 
    608 As indicated on \autoref{fig:index_vert} \jp{jpk} is the number of $w$-levels. 
    609 $gdepw_0(1)$ is the ocean surface. 
    610 There are at most \jp{jpk}-1 $t$-points inside the ocean, 
    611 the additional $t$-point at $jk = jpk$ is below the sea floor and is not used. 
    612 The vertical location of $w$- and $t$-levels is defined from the analytic expression of the depth $z_0(k)$ whose 
    613 analytical derivative with respect to $k$ provides the vertical scale factors. 
    614 The user must provide the analytical expression of both $z_0$ and its first derivative with respect to $k$. 
    615 This is done in routine \mdl{domzgr} through statement functions, 
    616 using parameters provided in the \ngn{namcfg} namelist. 
    617  
    618 It is possible to define a simple regular vertical grid by giving zero stretching (\np{ppacr}~\forcode{= 0}). 
    619 In that case, the parameters \jp{jpk} (number of $w$-levels) and 
    620 \np{pphmax} (total ocean depth in meters) fully define the grid. 
    621  
    622 For climate-related studies it is often desirable to concentrate the vertical resolution near the ocean surface. 
    623 The following function is proposed as a standard for a $z$-coordinate (with either full or partial steps):  
    624 \begin{gather} 
    625   \label{eq:DOM_zgr_ana_1} 
    626     z_0  (k) = h_{sur} - h_0 \; k - \; h_1 \; \log  \big[ \cosh ((k - h_{th}) / h_{cr}) \big] \\ 
    627     e_3^0(k) = \lt|    - h_0      -    h_1 \; \tanh \big[        (k - h_{th}) / h_{cr}  \big] \rt| 
    628 \end{gather} 
    629 where $k = 1$ to \jp{jpk} for $w$-levels and $k = 1$ to $k = 1$ for $T-$levels. 
    630 Such an expression allows us to define a nearly uniform vertical location of levels at the ocean top and bottom with 
    631 a smooth hyperbolic tangent transition in between (\autoref{fig:zgr}). 
    632  
    633 If the ice shelf cavities are opened (\np{ln\_isfcav}~\forcode{= .true.}), the definition of $z_0$ is the same. 
    634 However, definition of $e_3^0$ at $t$- and $w$-points is respectively changed to: 
    635 \begin{equation} 
    636   \label{eq:DOM_zgr_ana_2} 
    637   \begin{split} 
    638     e_3^T(k) &= z_W (k + 1) - z_W (k    ) \\ 
    639     e_3^W(k) &= z_T (k    ) - z_T (k - 1) 
    640   \end{split} 
    641 \end{equation} 
    642 This formulation decrease the self-generated circulation into the ice shelf cavity  
    643 (which can, in extreme case, leads to blow up).\\ 
    644   
    645 The most used vertical grid for ORCA2 has $10~m$ ($500~m$) resolution in the surface (bottom) layers and 
    646 a depth which varies from 0 at the sea surface to a minimum of $-5000~m$. 
    647 This leads to the following conditions: 
    648 \begin{equation} 
    649   \label{eq:DOM_zgr_coef} 
    650   \begin{array}{ll} 
    651     e_3 (1   + 1/2) =  10. & z(1  ) =     0. \\ 
    652     e_3 (jpk - 1/2) = 500. & z(jpk) = -5000. 
    653   \end{array} 
    654 \end{equation} 
    655  
    656 With the choice of the stretching $h_{cr} = 3$ and the number of levels \jp{jpk}~$= 31$, 
    657 the four coefficients $h_{sur}$, $h_0$, $h_1$, and $h_{th}$ in 
    658 \autoref{eq:DOM_zgr_ana_2} have been determined such that 
    659 \autoref{eq:DOM_zgr_coef} is satisfied, through an optimisation procedure using a bisection method. 
    660 For the first standard ORCA2 vertical grid this led to the following values: 
    661 $h_{sur} = 4762.96$, $h_0 = 255.58, h_1 = 245.5813$, and $h_{th} = 21.43336$. 
    662 The resulting depths and scale factors as a function of the model levels are shown in 
    663 \autoref{fig:zgr} and given in \autoref{tab:orca_zgr}. 
    664 Those values correspond to the parameters \np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist. 
    665  
    666 Rather than entering parameters $h_{sur}$, $h_0$, and $h_1$ directly, it is possible to recalculate them. 
    667 In that case the user sets \np{ppsur}~$=$~\np{ppa0}~$=$~\np{ppa1}~$= 999999$., 
    668 in \ngn{namcfg} namelist, and specifies instead the four following parameters: 
     489The choice of a vertical coordinate is made when setting up the configuration; 
     490it is not intended to be an option which can be changed in the middle of an 
     491experiment. The one exception to this statement being the choice of linear or 
     492non-linear free surface. In v4.0 the linear free surface option is implemented 
     493as a special case of the non-linear free surface. This is computationally 
     494wasteful since it uses the structures for time-varying 3D metrics for fields 
     495that (in the linear free surface case) are fixed. However, the linear 
     496free-surface is rarely used and implementing it this way means a single configuration 
     497file can support both options. 
     498 
     499By default a non-linear free surface is used (\np{ln\_linssh} set to \forcode{ = 
     500.false.} in \ngn{namdom}): the coordinate follow the time-variation of the free 
     501surface so that the transformation is time dependent: $z(i,j,k,t)$ 
     502(\eg \autoref{fig:z_zps_s_sps}f).  When a linear free surface is assumed 
     503(\np{ln\_linssh} set to \forcode{ = .true.} in \ngn{namdom}), the vertical 
     504coordinates are fixed in time, but the seawater can move up and down across the 
     505$z_0$ surface (in other words, the top of the ocean in not a rigid lid). 
     506 
     507Note that settings: \np{ln\_zco}, \np{ln\_zps}, \np{ln\_sco} and \np{ln\_isfcav} mentioned 
     508in the following sections appear to be namelist options but they are no longer truly 
     509namelist options for NEMO. Their value is written to and read from the domain configuration file 
     510and they should be treated as fixed parameters for a particular configuration. They are 
     511namelist options for the \forcode{DOMAINcfg} tool that can be used to build the 
     512configuration file and serve both to provide a record of the choices made whilst building the 
     513configuration and to trigger appropriate code blocks within NEMO. 
     514These values should not be altered in the \np{cn\_domcfg} file. 
     515 
     516\medskip 
     517The decision on these choices must be made when the \np{cn\_domcfg} file is constructed. 
     518Three main choices are offered (\autoref{fig:z_zps_s_sps}a-c): 
     519 
    669520\begin{itemize} 
    670 \item 
    671   \np{ppacr}~$= h_{cr}$: stretching factor (nondimensional). 
    672   The larger \np{ppacr}, the smaller the stretching. 
    673   Values from $3$ to $10$ are usual. 
    674 \item 
    675   \np{ppkth}~$= h_{th}$: is approximately the model level at which maximum stretching occurs 
    676   (nondimensional, usually of order 1/2 or 2/3 of \jp{jpk}) 
    677 \item 
    678   \np{ppdzmin}: minimum thickness for the top layer (in meters). 
    679 \item 
    680   \np{pphmax}: total depth of the ocean (meters). 
     521\item $z$-coordinate with full step bathymetry (\np{ln\_zco}\forcode{ = .true.}), 
     522\item $z$-coordinate with partial step ($zps$) bathymetry (\np{ln\_zps}\forcode{ = .true.}), 
     523\item Generalized, $s$-coordinate (\np{ln\_sco}\forcode{ = .true.}). 
    681524\end{itemize} 
    682 As an example, for the $45$ layers used in the DRAKKAR configuration those parameters are: 
    683 \jp{jpk}~$= 46$, \np{ppacr}~$= 9$, \np{ppkth}~$= 23.563$, \np{ppdzmin}~$= 6~m$, \np{pphmax}~$= 5750~m$. 
    684  
    685 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    686 \begin{table} 
    687   \begin{center} 
    688     \begin{tabular}{c||r|r|r|r} 
    689       \hline 
    690       \textbf{LEVEL} & \textbf{gdept\_1d} & \textbf{gdepw\_1d} & \textbf{e3t\_1d } & \textbf{e3w\_1d} \\ 
    691       \hline 
    692       1              & \textbf{     5.00} &               0.00 & \textbf{   10.00} &            10.00 \\ 
    693       \hline 
    694       2              & \textbf{    15.00} &              10.00 & \textbf{   10.00} &            10.00 \\ 
    695       \hline 
    696       3              & \textbf{    25.00} &              20.00 & \textbf{   10.00} &            10.00 \\ 
    697       \hline 
    698       4              & \textbf{    35.01} &              30.00 & \textbf{   10.01} &            10.00 \\ 
    699       \hline 
    700       5              & \textbf{    45.01} &              40.01 & \textbf{   10.01} &            10.01 \\ 
    701       \hline 
    702       6              & \textbf{    55.03} &              50.02 & \textbf{   10.02} &            10.02 \\ 
    703       \hline 
    704       7              & \textbf{    65.06} &              60.04 & \textbf{   10.04} &            10.03 \\ 
    705       \hline 
    706       8              & \textbf{    75.13} &              70.09 & \textbf{   10.09} &            10.06 \\ 
    707       \hline 
    708       9              & \textbf{    85.25} &              80.18 & \textbf{   10.17} &            10.12 \\ 
    709       \hline 
    710       10             & \textbf{    95.49} &              90.35 & \textbf{   10.33} &            10.24 \\ 
    711       \hline 
    712       11             & \textbf{   105.97} &             100.69 & \textbf{   10.65} &            10.47 \\ 
    713       \hline 
    714       12             & \textbf{   116.90} &             111.36 & \textbf{   11.27} &            10.91 \\ 
    715       \hline 
    716       13             & \textbf{   128.70} &             122.65 & \textbf{   12.47} &            11.77 \\ 
    717       \hline 
    718       14             & \textbf{   142.20} &             135.16 & \textbf{   14.78} &            13.43 \\ 
    719       \hline 
    720       15             & \textbf{   158.96} &             150.03 & \textbf{   19.23} &            16.65 \\ 
    721       \hline 
    722       16             & \textbf{   181.96} &             169.42 & \textbf{   27.66} &            22.78 \\ 
    723       \hline 
    724       17             & \textbf{   216.65} &             197.37 & \textbf{   43.26} &            34.30 \\ 
    725       \hline 
    726       18             & \textbf{   272.48} &             241.13 & \textbf{   70.88} &            55.21 \\ 
    727       \hline 
    728       19             & \textbf{   364.30} &             312.74 & \textbf{  116.11} &            90.99 \\ 
    729       \hline 
    730       20             & \textbf{   511.53} &             429.72 & \textbf{  181.55} &           146.43 \\ 
    731       \hline 
    732       21             & \textbf{   732.20} &             611.89 & \textbf{  261.03} &           220.35 \\ 
    733       \hline 
    734       22             & \textbf{  1033.22} &             872.87 & \textbf{  339.39} &           301.42 \\ 
    735       \hline 
    736       23             & \textbf{  1405.70} &            1211.59 & \textbf{  402.26} &           373.31 \\ 
    737       \hline 
    738       24             & \textbf{  1830.89} &            1612.98 & \textbf{  444.87} &           426.00 \\ 
    739       \hline 
    740       25             & \textbf{  2289.77} &            2057.13 & \textbf{  470.55} &           459.47 \\ 
    741       \hline 
    742       26             & \textbf{  2768.24} &            2527.22 & \textbf{  484.95} &           478.83 \\ 
    743       \hline 
    744       27             & \textbf{  3257.48} &            3011.90 & \textbf{  492.70} &           489.44 \\ 
    745       \hline 
    746       28             & \textbf{  3752.44} &            3504.46 & \textbf{  496.78} &           495.07 \\ 
    747       \hline 
    748       29             & \textbf{  4250.40} &            4001.16 & \textbf{  498.90} &           498.02 \\ 
    749       \hline 
    750       30             & \textbf{  4749.91} &            4500.02 & \textbf{  500.00} &           499.54 \\ 
    751       \hline 
    752       31             & \textbf{  5250.23} &            5000.00 & \textbf{  500.56} &           500.33 \\ 
    753       \hline 
    754     \end{tabular} 
    755   \end{center} 
    756   \caption{ 
    757     \protect\label{tab:orca_zgr} 
    758     Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration as computed from 
    759     \autoref{eq:DOM_zgr_ana_2} using the coefficients given in \autoref{eq:DOM_zgr_coef} 
    760   } 
    761 \end{table} 
    762 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    763  
    764 % ------------------------------------------------------------------------------------------------------------- 
    765 %        z-coordinate with partial step 
    766 % ------------------------------------------------------------------------------------------------------------- 
    767 \subsection{$Z$-coordinate with partial step (\protect\np{ln\_zps}~\forcode{= .true.})} 
    768 \label{subsec:DOM_zps} 
    769 %--------------------------------------------namdom------------------------------------------------------- 
    770  
    771 \nlst{namdom}  
    772 %-------------------------------------------------------------------------------------------------------------- 
    773  
    774 In $z$-coordinate partial step, 
    775 the depths of the model levels are defined by the reference analytical function $z_0(k)$ as described in 
    776 the previous section, \textit{except} in the bottom layer. 
    777 The thickness of the bottom layer is allowed to vary as a function of geographical location $(\lambda,\varphi)$ to 
    778 allow a better representation of the bathymetry, especially in the case of small slopes 
    779 (where the bathymetry varies by less than one level thickness from one grid point to the next). 
    780 The reference layer thicknesses $e_{3t}^0$ have been defined in the absence of bathymetry. 
    781 With partial steps, layers from 1 to \jp{jpk}-2 can have a thickness smaller than $e_{3t}(jk)$. 
    782 The model deepest layer (\jp{jpk}-1) is allowed to have either a smaller or larger thickness than $e_{3t}(jpk)$: 
    783 the maximum thickness allowed is $2*e_{3t}(jpk - 1)$. 
    784 This has to be kept in mind when specifying values in \ngn{namdom} namelist, 
    785 as the maximum depth \np{pphmax} in partial steps: 
    786 for example, with \np{pphmax}~$= 5750~m$ for the DRAKKAR 45 layer grid, 
    787 the maximum ocean depth allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk - 1)$ being $250~m$). 
    788 Two variables in the namdom namelist are used to define the partial step vertical grid. 
    789 The mimimum water thickness (in meters) allowed for a cell partially filled with bathymetry at level jk is 
    790 the minimum of \np{rn\_e3zps\_min} (thickness in meters, usually $20~m$) or $e_{3t}(jk)*$\np{rn\_e3zps\_rat} 
    791 (a fraction, usually 10\%, of the default thickness $e_{3t}(jk)$). 
    792  
    793 \gmcomment{ \colorbox{yellow}{Add a figure here of pstep especially at last ocean level }  } 
    794  
    795 % ------------------------------------------------------------------------------------------------------------- 
    796 %        s-coordinate 
    797 % ------------------------------------------------------------------------------------------------------------- 
    798 \subsection{$S$-coordinate (\protect\np{ln\_sco}~\forcode{= .true.})} 
    799 \label{subsec:DOM_sco} 
    800 %------------------------------------------nam_zgr_sco--------------------------------------------------- 
    801 % 
    802 %\nlst{namzgr_sco}  
    803 %-------------------------------------------------------------------------------------------------------------- 
    804 Options are defined in \ngn{namzgr\_sco}. 
    805 In $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}), the depth and thickness of the model levels are defined from 
    806 the product of a depth field and either a stretching function or its derivative, respectively: 
    807  
    808 \begin{align*} 
    809   % \label{eq:DOM_sco_ana} 
    810   z(k)   &= h(i,j) \; z_0 (k) \\ 
    811   e_3(k) &= h(i,j) \; z_0'(k) 
    812 \end{align*} 
    813  
    814 where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point location in the horizontal and 
    815 $z_0(k)$ is a function which varies from $0$ at the sea surface to $1$ at the ocean bottom. 
    816 The depth field $h$ is not necessary the ocean depth, 
    817 since a mixed step-like and bottom-following representation of the topography can be used 
    818 (\autoref{fig:z_zps_s_sps}) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}). 
    819 The namelist parameter \np{rn\_rmax} determines the slope at which 
    820 the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate. 
    821 The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} as 
    822 the minimum and maximum depths at which the terrain-following vertical coordinate is calculated. 
    823  
    824 Options for stretching the coordinate are provided as examples, 
    825 but care must be taken to ensure that the vertical stretch used is appropriate for the application. 
    826  
    827 The original default NEMO s-coordinate stretching is available if neither of the other options are specified as true 
    828 (\np{ln\_s\_SH94}~\forcode{= .false.} and \np{ln\_s\_SF12}~\forcode{= .false.}). 
    829 This uses a depth independent $\tanh$ function for the stretching \citep{Madec_al_JPO96}: 
    830  
    831 \[ 
    832   z = s_{min} + C (s) (H - s_{min}) 
    833   % \label{eq:SH94_1} 
    834 \] 
    835  
    836 where $s_{min}$ is the depth at which the $s$-coordinate stretching starts and 
    837 allows a $z$-coordinate to placed on top of the stretched coordinate, 
    838 and $z$ is the depth (negative down from the asea surface). 
    839 \begin{gather*} 
    840   s = - \frac{k}{n - 1} \quad \text{and} \quad 0 \leq k \leq n - 1 
    841   % \label{eq:DOM_s} 
    842  \\ 
    843   % \label{eq:DOM_sco_function} 
    844   C(s) = \frac{[\tanh(\theta \, (s + b)) - \tanh(\theta \, b)]}{2 \; \sinh(\theta)} 
    845 \end{gather*} 
    846  
    847 A stretching function, 
    848 modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_s\_SH94}~\forcode{= .true.}), 
    849 is also available and is more commonly used for shelf seas modelling: 
    850  
    851 \[ 
    852   C(s) =   (1 - b) \frac{\sinh(\theta s)}{\sinh(\theta)} 
    853          + b       \frac{\tanh \lt[ \theta \lt(s + \frac{1}{2} \rt) \rt] -   \tanh \lt( \frac{\theta}{2} \rt)} 
    854                         {                                                  2 \tanh \lt( \frac{\theta}{2} \rt)} 
    855   % \label{eq:SH94_2} 
    856 \] 
    857  
    858 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    859 \begin{figure}[!ht] 
    860   \begin{center} 
    861     \includegraphics[]{Fig_sco_function} 
    862     \caption{ 
    863       \protect\label{fig:sco_function} 
    864       Examples of the stretching function applied to a seamount; 
    865       from left to right: surface, surface and bottom, and bottom intensified resolutions 
    866     } 
    867   \end{center} 
    868 \end{figure} 
    869 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    870  
    871 where $H_c$ is the critical depth (\np{rn\_hc}) at which the coordinate transitions from pure $\sigma$ to 
    872 the stretched coordinate, and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and 
    873 bottom control parameters such that $0 \leqslant \theta \leqslant 20$, and $0 \leqslant b \leqslant 1$. 
    874 $b$ has been designed to allow surface and/or bottom increase of the vertical resolution 
    875 (\autoref{fig:sco_function}). 
    876  
    877 Another example has been provided at version 3.5 (\np{ln\_s\_SF12}) that allows a fixed surface resolution in 
    878 an analytical terrain-following stretching \citet{Siddorn_Furner_OM12}. 
    879 In this case the a stretching function $\gamma$ is defined such that: 
    880  
    881 \begin{equation} 
    882   z = - \gamma h \quad \text{with} \quad 0 \leq \gamma \leq 1 
    883   % \label{eq:z} 
    884 \end{equation} 
    885  
    886 The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate: 
    887  
    888 \begin{gather*} 
    889   % \label{eq:DOM_gamma_deriv} 
    890   \gamma =   A \lt( \sigma   - \frac{1}{2} (\sigma^2     + f (\sigma)) \rt) 
    891            + B \lt( \sigma^3 - f           (\sigma) \rt) + f (\sigma)       \\ 
    892   \intertext{Where:} 
    893   % \label{eq:DOM_gamma} 
    894   f(\sigma) = (\alpha + 2) \sigma^{\alpha + 1} - (\alpha + 1) \sigma^{\alpha + 2} 
    895   \quad \text{and} \quad \sigma = \frac{k}{n - 1} 
    896 \end{gather*} 
    897  
    898 This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of 
    899 the user prescribed stretching parameter $\alpha$ (\np{rn\_alpha}) that stretches towards 
    900 the surface ($\alpha > 1.0$) or the bottom ($\alpha < 1.0$) and 
    901 user prescribed surface (\np{rn\_zs}) and bottom depths. 
    902 The bottom cell depth in this example is given as a function of water depth: 
    903  
    904 \[ 
    905   % \label{eq:DOM_zb} 
    906   Z_b = h a + b 
    907 \] 
    908  
    909 where the namelist parameters \np{rn\_zb\_a} and \np{rn\_zb\_b} are $a$ and $b$ respectively. 
    910  
    911 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    912 \begin{figure}[!ht] 
    913   \includegraphics[]{Fig_DOM_compare_coordinates_surface} 
    914   \caption{ 
    915     A comparison of the \citet{Song_Haidvogel_JCP94} $S$-coordinate (solid lines), 
    916     a 50 level $Z$-coordinate (contoured surfaces) and 
    917     the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in the surface $100~m$ for 
    918     a idealised bathymetry that goes from $50~m$ to $5500~m$ depth. 
    919     For clarity every third coordinate surface is shown. 
    920   } 
    921   \label{fig:fig_compare_coordinates_surface} 
    922 \end{figure} 
    923  % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    924  
    925 This gives a smooth analytical stretching in computational space that is constrained to 
    926 given specified surface and bottom grid cell thicknesses in real space. 
    927 This is not to be confused with the hybrid schemes that 
    928 superimpose geopotential coordinates on terrain following coordinates thus 
    929 creating a non-analytical vertical coordinate that 
    930 therefore may suffer from large gradients in the vertical resolutions. 
    931 This stretching is less straightforward to implement than the \citet{Song_Haidvogel_JCP94} stretching, 
    932 but has the advantage of resolving diurnal processes in deep water and has generally flatter slopes. 
    933  
    934 As with the \citet{Song_Haidvogel_JCP94} stretching the stretch is only applied at depths greater than 
    935 the critical depth $h_c$. 
    936 In this example two options are available in depths shallower than $h_c$, 
    937 with pure sigma being applied if the \np{ln\_sigcrit} is true and pure z-coordinates if it is false 
    938 (the z-coordinate being equal to the depths of the stretched coordinate at $h_c$). 
    939  
    940 Minimising the horizontal slope of the vertical coordinate is important in terrain-following systems as 
    941 large slopes lead to hydrostatic consistency. 
    942 A hydrostatic consistency parameter diagnostic following \citet{Haney1991} has been implemented, 
    943 and is output as part of the model mesh file at the start of the run. 
    944  
    945 % ------------------------------------------------------------------------------------------------------------- 
    946 %        z*- or s*-coordinate 
    947 % ------------------------------------------------------------------------------------------------------------- 
    948 \subsection{\zstar- or \sstar-coordinate (\protect\np{ln\_linssh}~\forcode{= .false.})} 
    949 \label{subsec:DOM_zgr_star} 
    950  
    951 This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site. 
    952  
    953 %gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances 
     525 
     526Additionally, hybrid combinations of the three main coordinates are available: 
     527$s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps}d and \autoref{fig:z_zps_s_sps}e). 
     528 
     529A further choice related to vertical coordinate concerns the presence (or not) of ocean 
     530cavities beneath ice shelves within the model domain.  A setting of \np{ln\_isfcav} as 
     531\forcode{.true.} indicates that the domain contains  ocean cavities, otherwise the top, 
     532wet layer of the ocean will always be at the ocean surface.  This option is currently only 
     533available for $z$- or $zps$-coordinates. In the latter case, partial steps are also applied 
     534at the ocean/ice shelf interface. 
     535 
     536Within the model, the arrays describing the grid point depths and vertical scale factors 
     537are three set of three dimensional arrays $(i,j,k)$ defined at \textit{before}, 
     538\textit{now} and \textit{after} time step.  The time at which they are defined is 
     539indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively.  They are updated at each 
     540model time step. The initial fixed reference coordinate system is held in variable names 
     541with a $\_0$ suffix.  When the linear free surface option is used 
     542(\np{ln\_linssh}\forcode{ = .true.}), \textit{before}, \textit{now} and \textit{after} 
     543arrays are initially set to their reference counterpart and remain fixed. 
     544 
     545\subsubsection{Required fields} 
     546\label{sec:DOM_zgr_fields} 
     547The explicit specification of a range of fields related to the vertical grid are required for the definition of a configuration. These include: 
     548 
     549\begin{Verbatim}[fontsize=\tiny] 
     550int    ln_zco, ln_zps, ln_sco            /* flags for z-coord, z-coord with partial steps and s-coord    */ 
     551int    ln_isfcav                         /* flag  for ice shelf cavities                                 */ 
     552double e3t_1d, e3w_1d                    /* reference vertical scale factors at T and W points           */ 
     553double e3t_0, e3u_0, e3v_0, e3f_0, e3w_0 /* vertical scale factors 3D coordinate at T,U,V,F and W points */ 
     554double e3uw_0, e3vw_0                    /* vertical scale factors 3D coordinate at UW and VW points     */ 
     555int    bottom_level, top_level           /* last wet T-points, 1st wet T-points (for ice shelf cavities) */ 
     556                                         /* For reference:                                               */ 
     557float  bathy_metry                       /* bathymetry used in setting top and bottom levels             */ 
     558\end{Verbatim} 
     559 
     560This set of vertical metrics is sufficient to describe the initial depth and thickness of 
     561every gridcell in the model regardless of the choice of vertical coordinate. With constant 
     562z-levels, e3 metrics will be uniform across each horizontal level. In the partial step 
     563case each e3 at the \np{bottom\_level} (and, possibly, \np{top\_level} if ice cavities are 
     564present) may vary from its horizontal neighbours. And, in s-coordinates, variations can 
     565occur throughout the water column. With the non-linear free-surface, all the coordinates 
     566behave more like the s-coordinate in that variations occurr throughout the water column 
     567with displacements related to the sea surface height. These variations are typically much 
     568smaller than those arising from bottom fitted coordinates. The values for vertical metrics 
     569supplied in the domain configuration file can be considered as those arising from a flat 
     570sea surface with zero elevation. 
     571 
     572The \np{bottom\_level} and \np{top\_level} 2D arrays define the \np{bottom\_level} and top 
     573wet levels in each grid column. Without ice cavities, \np{top\_level} is essentially a land 
     574mask (0 on land; 1 everywhere else). With ice cavities, \np{top\_level} determines the 
     575first wet point below the overlying ice shelf. 
     576 
     577 
    954578 
    955579% ------------------------------------------------------------------------------------------------------------- 
    956580%        level bathymetry and mask  
    957581% ------------------------------------------------------------------------------------------------------------- 
    958 \subsection{Level bathymetry and mask} 
     582\subsubsection{Level bathymetry and mask} 
    959583\label{subsec:DOM_msk} 
    960584 
    961 Whatever the vertical coordinate used, the model offers the possibility of representing the bottom topography with 
    962 steps that follow the face of the model cells (step like topography) \citep{Madec_al_JPO96}. 
    963 The distribution of the steps in the horizontal is defined in a 2D integer array, mbathy, which 
    964 gives the number of ocean levels (\ie those that are not masked) at each $t$-point. 
    965 mbathy is computed from the meter bathymetry using the definiton of gdept as the number of $t$-points which 
    966 gdept $\leq$ bathy. 
    967  
    968 Modifications of the model bathymetry are performed in the \textit{bat\_ctl} routine (see \mdl{domzgr} module) after 
    969 mbathy is computed. 
    970 Isolated grid points that do not communicate with another ocean point at the same level are eliminated. 
    971  
    972 As for the representation of bathymetry, a 2D integer array, misfdep, is created. 
    973 misfdep defines the level of the first wet $t$-point. 
    974 All the cells between $k = 1$ and $misfdep(i,j) - 1$ are masked. 
    975 By default, $misfdep(:,:) = 1$ and no cells are masked. 
    976  
    977 In case of ice shelf cavities, modifications of the model bathymetry and ice shelf draft into  
    978 the cavities are performed in the \textit{zgr\_isf} routine. 
    979 The compatibility between ice shelf draft and bathymetry is checked. 
    980 All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded (\ie masked). 
    981 If only one cell on the water column is opened at $t$-, $u$- or $v$-points, 
    982 the bathymetry or the ice shelf draft is dug to fit this constrain. 
    983 If the incompatibility is too strong (need to dig more than 1 cell), the cell is masked. 
    984  
    985 From the \textit{mbathy} and \textit{misfdep} array, the mask fields are defined as follows: 
     585 
     586From \np{top\_level} and \np{bottom\_level} fields, the mask fields are defined as follows: 
    986587\begin{alignat*}{2} 
    987588  tmask(i,j,k) &= &  & 
    988589    \begin{cases} 
    989                   0 &\text{if $                  k  <    misfdep(i,j)$} \\ 
    990                   1 &\text{if $misfdep(i,j) \leq k \leq   mbathy(i,j)$} \\ 
    991                   0 &\text{if $                  k  >     mbathy(i,j)$} 
     590                  0 &\text{if $                  k  <    top\_level(i,j)$} \\ 
     591                  1 &\text{if $bottom\_level(i,j) \leq k \leq   top\_level(i,j)$} \\ 
     592                  0 &\text{if $                  k  >     bottom\_level(i,j)$} 
    992593    \end{cases} 
    993594  \\ 
     
    1005606exactly in the same way as for the bottom boundary. 
    1006607 
    1007 The specification of closed lateral boundaries requires that at least 
    1008 the first and last rows and columns of the \textit{mbathy} array are set to zero. 
    1009 In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to 
    1010 the second one and its first column equal to the last but one (and so too the mask arrays) 
    1011 (see \autoref{fig:LBC_jperio}). 
     608%% The specification of closed lateral boundaries requires that at least 
     609%% the first and last rows and columns of the \textit{mbathy} array are set to zero. 
     610%% In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to 
     611%% the second one and its first column equal to the last but one (and so too the mask arrays) 
     612%% (see \autoref{fig:LBC_jperio}). 
     613 
     614 
     615%------------------------------------------------------------------------------------------------- 
     616%        Closed seas  
     617%------------------------------------------------------------------------------------------------- 
     618\subsection{Closed seas} \label{subsec:DOM_closea}  
     619 
     620When a global ocean is coupled to an atmospheric model it is better to represent all large 
     621water bodies (\eg great lakes, Caspian sea...) even if the model resolution does not allow 
     622their communication with the rest of the ocean.  This is unnecessary when the ocean is 
     623forced by fixed atmospheric conditions, so these seas can be removed from the ocean 
     624domain.  The user has the option to set the bathymetry in closed seas to zero (see 
     625\autoref{sec:MISC_closea}) and to optionally decide on the fate of any freshwater 
     626imbalance over the area. The options are explained in \autoref{sec:MISC_closea} but it 
     627should be noted here that a successful use of these options requires appropriate mask 
     628fields to be present in the domain configuration file. Among the possibilities are: 
     629 
     630\begin{Verbatim}[fontsize=\tiny] 
     631int    closea_mask          /* non-zero values in closed sea areas for optional masking                  */ 
     632int    closea_mask_rnf      /* non-zero values in closed sea areas with runoff locations (precip only)   */ 
     633int    closea_mask_emp      /* non-zero values in closed sea areas with runoff locations (total emp)     */ 
     634\end{Verbatim} 
     635 
     636% ------------------------------------------------------------------------------------------------------------- 
     637%        Grid files 
     638% ------------------------------------------------------------------------------------------------------------- 
     639\subsection{Output grid files} 
     640\label{subsec:DOM_meshmask} 
     641 
     642\nlst{namcfg} 
     643 
     644Most of the arrays relating to a particular ocean model configuration dicussed in this 
     645chapter (grid-point position, scale factors) can be saved in a file if namelist parameter 
     646\np{ln\_write\_cfg} (namelist \ngn{namcfg}) is set to \forcode{.true.}; the output 
     647filename is set thorugh parameter \np{cn\_domcfg\_out}. This is only really useful 
     648if the fields are computed in subroutines \mdl{usrdef\_hgr} or \mdl{usrdef\_zgr} and 
     649checking or confirmation is required. 
     650 
     651\nlst{namdom} 
     652 
     653Alternatively, all the arrays relating to a particular ocean model configuration 
     654(grid-point position, scale factors, depths and masks) can be saved in a file called 
     655\texttt{mesh\_mask} if namelist parameter \np{ln\_meshmask} (namelist \ngn{namdom}) is set 
     656to \forcode{.true.}. This file contains additional fields that can be useful for 
     657post-processing applications 
    1012658 
    1013659% ================================================================ 
    1014660% Domain: Initial State (dtatsd & istate) 
    1015661% ================================================================ 
    1016 \section{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} 
     662\section[Initial state (\textit{istate.F90} and \textit{dtatsd.F90})] 
     663{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} 
    1017664\label{sec:DTA_tsd} 
    1018665%-----------------------------------------namtsd------------------------------------------- 
    1019  
    1020666\nlst{namtsd}  
    1021667%------------------------------------------------------------------------------------------ 
    1022668 
    1023 Options are defined in \ngn{namtsd}. 
    1024 By default, the ocean start from rest (the velocity field is set to zero) and the initialization of temperature and 
    1025 salinity fields is controlled through the \np{ln\_tsd\_ini} namelist parameter. 
     669Basic initial state options are defined in \ngn{namtsd}.  By default, the ocean starts 
     670from rest (the velocity field is set to zero) and the initialization of temperature and 
     671salinity fields is controlled through the \np{ln\_tsd\_init} namelist parameter. 
     672 
    1026673\begin{description} 
    1027 \item[\np{ln\_tsd\_init}~\forcode{= .true.}] 
    1028   use a T and S input files that can be given on the model grid itself or on their native input data grid. 
    1029   In the latter case, 
    1030   the data will be interpolated on-the-fly both in the horizontal and the vertical to the model grid 
    1031   (see \autoref{subsec:SBC_iof}). 
    1032   The information relative to the input files are given in the \np{sn\_tem} and \np{sn\_sal} structures. 
    1033   The computation is done in the \mdl{dtatsd} module. 
    1034 \item[\np{ln\_tsd\_init}~\forcode{= .false.}] 
    1035   use constant salinity value of $35.5~psu$ and an analytical profile of temperature 
    1036   (typical of the tropical ocean), see \rou{istate\_t\_s} subroutine called from \mdl{istate} module. 
     674\item[\np{ln\_tsd\_init}\forcode{= .true.}] 
     675  Use T and S input files that can be given on the model grid itself or on their native 
     676  input data grids.  In the latter case, the data will be interpolated on-the-fly both in 
     677  the horizontal and the vertical to the model grid (see \autoref{subsec:SBC_iof}).  The 
     678  information relating to the input files are specified in the \np{sn\_tem} and 
     679  \np{sn\_sal} structures.  The computation is done in the \mdl{dtatsd} module. 
     680\item[\np{ln\_tsd\_init}\forcode{= .false.}] 
     681  Initial values for T and S are set via a user supplied \rou{usr\_def\_istate} routine 
     682  contained in \mdl{userdef\_istate}. The default version sets horizontally uniform T and 
     683  profiles as used in the  GYRE configuration (see \autoref{sec:CFG_gyre}). 
    1037684\end{description} 
    1038685 
  • NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_DYN.tex

    r10499 r11422  
    6565%           Horizontal divergence and relative vorticity 
    6666%-------------------------------------------------------------------------------------------------------------- 
    67 \subsection{Horizontal divergence and relative vorticity (\protect\mdl{divcur})} 
     67\subsection[Horizontal divergence and relative vorticity (\textit{divcur.F90})] 
     68{Horizontal divergence and relative vorticity (\protect\mdl{divcur})} 
    6869\label{subsec:DYN_divcur} 
    6970 
     
    101102%           Sea Surface Height evolution 
    102103%-------------------------------------------------------------------------------------------------------------- 
    103 \subsection{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} 
     104\subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})] 
     105{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} 
    104106\label{subsec:DYN_sshwzv} 
    105107 
     
    127129Replacing $T$ by the number $1$ in the tracer equation and summing over the water column must lead to 
    128130the sea surface height equation otherwise tracer content will not be conserved 
    129 \citep{Griffies_al_MWR01, Leclair_Madec_OM09}. 
     131\citep{griffies.pacanowski.ea_MWR01, leclair.madec_OM09}. 
    130132 
    131133The vertical velocity is computed by an upward integration of the horizontal divergence starting at the bottom, 
     
    181183%        Vorticity term  
    182184% ------------------------------------------------------------------------------------------------------------- 
    183 \subsection{Vorticity term (\protect\mdl{dynvor})} 
     185\subsection[Vorticity term (\textit{dynvor.F90})] 
     186{Vorticity term (\protect\mdl{dynvor})} 
    184187\label{subsec:DYN_vor} 
    185188%------------------------------------------nam_dynvor---------------------------------------------------- 
     
    203206%                 enstrophy conserving scheme 
    204207%------------------------------------------------------------- 
    205 \subsubsection{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})} 
     208\subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens = .true.})] 
     209{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})} 
    206210\label{subsec:DYN_vor_ens} 
    207211 
     
    226230%                 energy conserving scheme 
    227231%------------------------------------------------------------- 
    228 \subsubsection{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})} 
     232\subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene = .true.})] 
     233{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})} 
    229234\label{subsec:DYN_vor_ene} 
    230235 
     
    246251%                 mix energy/enstrophy conserving scheme 
    247252%------------------------------------------------------------- 
    248 \subsubsection{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.}) } 
     253\subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix = .true.})] 
     254{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.})} 
    249255\label{subsec:DYN_vor_mix} 
    250256 
     
    271277%                 energy and enstrophy conserving scheme 
    272278%------------------------------------------------------------- 
    273 \subsubsection{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.}) } 
     279\subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een = .true.})] 
     280{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.})} 
    274281\label{subsec:DYN_vor_een} 
    275282 
     
    287294Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves....} 
    288295 
    289 A very nice solution to the problem of double averaging was proposed by \citet{Arakawa_Hsu_MWR90}. 
     296A very nice solution to the problem of double averaging was proposed by \citet{arakawa.hsu_MWR90}. 
    290297The idea is to get rid of the double averaging by considering triad combinations of vorticity. 
    291298It is noteworthy that this solution is conceptually quite similar to the one proposed by 
    292 \citep{Griffies_al_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}). 
    293  
    294 The \citet{Arakawa_Hsu_MWR90} vorticity advection scheme for a single layer is modified  
    295 for spherical coordinates as described by \citet{Arakawa_Lamb_MWR81} to obtain the EEN scheme.  
     299\citep{griffies.gnanadesikan.ea_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}). 
     300 
     301The \citet{arakawa.hsu_MWR90} vorticity advection scheme for a single layer is modified  
     302for spherical coordinates as described by \citet{arakawa.lamb_MWR81} to obtain the EEN scheme.  
    296303First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point:  
    297304\[ 
     
    309316\begin{figure}[!ht] 
    310317  \begin{center} 
    311     \includegraphics[width=0.70\textwidth]{Fig_DYN_een_triad} 
     318    \includegraphics[width=\textwidth]{Fig_DYN_een_triad} 
    312319    \caption{ 
    313320      \protect\label{fig:DYN_een_triad} 
     
    327334(with a systematic reduction of $e_{3f}$ when a model level intercept the bathymetry) 
    328335that tends to reinforce the topostrophy of the flow 
    329 (\ie the tendency of the flow to follow the isobaths) \citep{Penduff_al_OS07}.  
     336(\ie the tendency of the flow to follow the isobaths) \citep{penduff.le-sommer.ea_OS07}.  
    330337 
    331338Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as 
     
    356363(\ie $\chi$=$0$) (see \autoref{subsec:C_vorEEN}).  
    357364Applied to a realistic ocean configuration, it has been shown that it leads to a significant reduction of 
    358 the noise in the vertical velocity field \citep{Le_Sommer_al_OM09}. 
     365the noise in the vertical velocity field \citep{le-sommer.penduff.ea_OM09}. 
    359366Furthermore, used in combination with a partial steps representation of bottom topography, 
    360367it improves the interaction between current and topography, 
    361 leading to a larger topostrophy of the flow \citep{Barnier_al_OD06, Penduff_al_OS07}.  
     368leading to a larger topostrophy of the flow \citep{barnier.madec.ea_OD06, penduff.le-sommer.ea_OS07}.  
    362369 
    363370%-------------------------------------------------------------------------------------------------------------- 
    364371%           Kinetic Energy Gradient term 
    365372%-------------------------------------------------------------------------------------------------------------- 
    366 \subsection{Kinetic energy gradient term (\protect\mdl{dynkeg})} 
     373\subsection[Kinetic energy gradient term (\textit{dynkeg.F90})] 
     374{Kinetic energy gradient term (\protect\mdl{dynkeg})} 
    367375\label{subsec:DYN_keg} 
    368376 
     
    384392%           Vertical advection term 
    385393%-------------------------------------------------------------------------------------------------------------- 
    386 \subsection{Vertical advection term (\protect\mdl{dynzad}) } 
     394\subsection[Vertical advection term (\textit{dynzad.F90})] 
     395{Vertical advection term (\protect\mdl{dynzad})} 
    387396\label{subsec:DYN_zad} 
    388397 
     
    403412When \np{ln\_dynzad\_zts}\forcode{ = .true.}, 
    404413a split-explicit time stepping with 5 sub-timesteps is used on the vertical advection term. 
    405 This option can be useful when the value of the timestep is limited by vertical advection \citep{Lemarie_OM2015}.  
     414This option can be useful when the value of the timestep is limited by vertical advection \citep{lemarie.debreu.ea_OM15}.  
    406415Note that in this case, 
    407416a similar split-explicit time stepping should be used on vertical advection of tracer to ensure a better stability, 
     
    430439%           Coriolis plus curvature metric terms 
    431440%-------------------------------------------------------------------------------------------------------------- 
    432 \subsection{Coriolis plus curvature metric terms (\protect\mdl{dynvor}) } 
     441\subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})] 
     442{Coriolis plus curvature metric terms (\protect\mdl{dynvor})} 
    433443\label{subsec:DYN_cor_flux} 
    434444 
     
    451461%           Flux form Advection term 
    452462%-------------------------------------------------------------------------------------------------------------- 
    453 \subsection{Flux form advection term (\protect\mdl{dynadv}) } 
     463\subsection[Flux form advection term (\textit{dynadv.F90})] 
     464{Flux form advection term (\protect\mdl{dynadv})} 
    454465\label{subsec:DYN_adv_flux} 
    455466 
     
    475486a $2^{nd}$ order centered finite difference scheme, CEN2, 
    476487or a $3^{rd}$ order upstream biased scheme, UBS. 
    477 The latter is described in \citet{Shchepetkin_McWilliams_OM05}. 
     488The latter is described in \citet{shchepetkin.mcwilliams_OM05}. 
    478489The schemes are selected using the namelist logicals \np{ln\_dynadv\_cen2} and \np{ln\_dynadv\_ubs}.  
    479490In flux form, the schemes differ by the choice of a space and time interpolation to define the value of 
     
    484495%                 2nd order centred scheme 
    485496%------------------------------------------------------------- 
    486 \subsubsection{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})} 
     497\subsubsection[CEN2: $2^{nd}$ order centred scheme (\forcode{ln_dynadv_cen2 = .true.})] 
     498{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})} 
    487499\label{subsec:DYN_adv_cen2} 
    488500 
     
    507519%                 UBS scheme 
    508520%------------------------------------------------------------- 
    509 \subsubsection{UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})} 
     521\subsubsection[UBS: Upstream Biased Scheme (\forcode{ln_dynadv_ubs = .true.})] 
     522{UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})} 
    510523\label{subsec:DYN_adv_ubs} 
    511524 
     
    523536where $u"_{i+1/2} =\delta_{i+1/2} \left[ {\delta_i \left[ u \right]} \right]$. 
    524537This results in a dissipatively dominant (\ie hyper-diffusive) truncation error 
    525 \citep{Shchepetkin_McWilliams_OM05}. 
    526 The overall performance of the advection scheme is similar to that reported in \citet{Farrow1995}. 
     538\citep{shchepetkin.mcwilliams_OM05}. 
     539The overall performance of the advection scheme is similar to that reported in \citet{farrow.stevens_JPO95}. 
    527540It is a relatively good compromise between accuracy and smoothness. 
    528541It is not a \emph{positive} scheme, meaning that false extrema are permitted. 
     
    542555while the second term, which is the diffusion part of the scheme, 
    543556is evaluated using the \textit{before} velocity (forward in time). 
    544 This is discussed by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme. 
     557This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme. 
    545558 
    546559Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) schemes only differ by 
    547560one coefficient. 
    548 Replacing $1/6$ by $1/8$ in (\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. 
     561Replacing $1/6$ by $1/8$ in (\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}. 
    549562This option is not available through a namelist parameter, since the $1/6$ coefficient is hard coded. 
    550563Nevertheless it is quite easy to make the substitution in the \mdl{dynadv\_ubs} module and obtain a QUICK scheme. 
     
    560573%           Hydrostatic pressure gradient term 
    561574% ================================================================ 
    562 \section{Hydrostatic pressure gradient (\protect\mdl{dynhpg})} 
     575\section[Hydrostatic pressure gradient (\textit{dynhpg.F90})] 
     576{Hydrostatic pressure gradient (\protect\mdl{dynhpg})} 
    563577\label{sec:DYN_hpg} 
    564578%------------------------------------------nam_dynhpg--------------------------------------------------- 
     
    582596%           z-coordinate with full step 
    583597%-------------------------------------------------------------------------------------------------------------- 
    584 \subsection{Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{ = .true.})} 
     598\subsection[Full step $Z$-coordinate (\forcode{ln_dynhpg_zco = .true.})] 
     599{Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{ = .true.})} 
    585600\label{subsec:DYN_hpg_zco} 
    586601 
     
    627642%           z-coordinate with partial step 
    628643%-------------------------------------------------------------------------------------------------------------- 
    629 \subsection{Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{ = .true.})} 
     644\subsection[Partial step $Z$-coordinate (\forcode{ln_dynhpg_zps = .true.})] 
     645{Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{ = .true.})} 
    630646\label{subsec:DYN_hpg_zps} 
    631647 
     
    652668 
    653669Pressure gradient formulations in an $s$-coordinate have been the subject of a vast number of papers 
    654 (\eg, \citet{Song1998, Shchepetkin_McWilliams_OM05}).  
     670(\eg, \citet{song_MWR98, shchepetkin.mcwilliams_OM05}).  
    655671A number of different pressure gradient options are coded but the ROMS-like, 
    656672density Jacobian with cubic polynomial method is currently disabled whilst known bugs are under investigation. 
    657673 
    658 $\bullet$ Traditional coding (see for example \citet{Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{ = .true.}) 
     674$\bullet$ Traditional coding (see for example \citet{madec.delecluse.ea_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{ = .true.}) 
    659675\begin{equation} 
    660676  \label{eq:dynhpg_sco} 
     
    679695$\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}\forcode{ = .true.}) 
    680696 
    681 $\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{Shchepetkin_McWilliams_OM05}  
     697$\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{shchepetkin.mcwilliams_OM05}  
    682698(\np{ln\_dynhpg\_djc}\forcode{ = .true.}) (currently disabled; under development) 
    683699 
    684700Note that expression \autoref{eq:dynhpg_sco} is commonly used when the variable volume formulation is activated 
    685701(\key{vvl}) because in that case, even with a flat bottom, 
    686 the coordinate surfaces are not horizontal but follow the free surface \citep{Levier2007}. 
     702the coordinate surfaces are not horizontal but follow the free surface \citep{levier.treguier.ea_rpt07}. 
    687703The pressure jacobian scheme (\np{ln\_dynhpg\_prj}\forcode{ = .true.}) is available as 
    688704an improved option to \np{ln\_dynhpg\_sco}\forcode{ = .true.} when \key{vvl} is active. 
     
    704720corresponds to the water replaced by the ice shelf. 
    705721This top pressure is constant over time. 
    706 A detailed description of this method is described in \citet{Losch2008}.\\ 
     722A detailed description of this method is described in \citet{losch_JGR08}.\\ 
    707723 
    708724The pressure gradient due to ocean load is computed using the expression \autoref{eq:dynhpg_sco} described in 
     
    712728%           Time-scheme 
    713729%-------------------------------------------------------------------------------------------------------------- 
    714 \subsection{Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{ = .true./.false.})} 
     730\subsection[Time-scheme (\forcode{ln_dynhpg_imp = .{true,false}.})] 
     731{Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{ = .\{true,false\}}.)} 
    715732\label{subsec:DYN_hpg_imp} 
    716733 
     
    722739the physical phenomenon that controls the time-step is internal gravity waves (IGWs). 
    723740A semi-implicit scheme for doubling the stability limit associated with IGWs can be used 
    724 \citep{Brown_Campana_MWR78, Maltrud1998}. 
     741\citep{brown.campana_MWR78, maltrud.smith.ea_JGR98}. 
    725742It involves the evaluation of the hydrostatic pressure gradient as 
    726743an average over the three time levels $t-\rdt$, $t$, and $t+\rdt$ 
     
    773790% Surface Pressure Gradient 
    774791% ================================================================ 
    775 \section{Surface pressure gradient (\protect\mdl{dynspg})} 
     792\section[Surface pressure gradient (\textit{dynspg.F90})] 
     793{Surface pressure gradient (\protect\mdl{dynspg})} 
    776794\label{sec:DYN_spg} 
    777795%-----------------------------------------nam_dynspg---------------------------------------------------- 
     
    790808which imposes a very small time step when an explicit time stepping is used. 
    791809Two methods are proposed to allow a longer time step for the three-dimensional equations:  
    792 the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:PE_flt}),  
     810the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:PE_flt?}),  
    793811and the split-explicit free surface described below. 
    794812The extra term introduced in the filtered method is calculated implicitly,  
     
    811829% Explicit free surface formulation 
    812830%-------------------------------------------------------------------------------------------------------------- 
    813 \subsection{Explicit free surface (\protect\key{dynspg\_exp})} 
     831\subsection[Explicit free surface (\texttt{\textbf{key\_dynspg\_exp}})] 
     832{Explicit free surface (\protect\key{dynspg\_exp})} 
    814833\label{subsec:DYN_spg_exp} 
    815834 
     
    837856% Split-explict free surface formulation 
    838857%-------------------------------------------------------------------------------------------------------------- 
    839 \subsection{Split-explicit free surface (\protect\key{dynspg\_ts})} 
     858\subsection[Split-explicit free surface (\texttt{\textbf{key\_dynspg\_ts}})] 
     859{Split-explicit free surface (\protect\key{dynspg\_ts})} 
    840860\label{subsec:DYN_spg_ts} 
    841861%------------------------------------------namsplit----------------------------------------------------------- 
     
    845865 
    846866The split-explicit free surface formulation used in \NEMO (\key{dynspg\_ts} defined), 
    847 also called the time-splitting formulation, follows the one proposed by \citet{Shchepetkin_McWilliams_OM05}. 
     867also called the time-splitting formulation, follows the one proposed by \citet{shchepetkin.mcwilliams_OM05}. 
    848868The general idea is to solve the free surface equation and the associated barotropic velocity equations with 
    849869a smaller time step than $\rdt$, the time step used for the three dimensional prognostic variables 
     
    862882\begin{equation} 
    863883  \label{eq:BT_dyn} 
    864   \frac{\partial {\rm \overline{{\bf U}}_h} }{\partial t}= 
    865   -f\;{\rm {\bf k}}\times {\rm \overline{{\bf U}}_h} 
    866   -g\nabla _h \eta -\frac{c_b^{\textbf U}}{H+\eta} \rm {\overline{{\bf U}}_h} + \rm {\overline{\bf G}} 
     884  \frac{\partial {\mathrm \overline{{\mathbf U}}_h} }{\partial t}= 
     885  -f\;{\mathrm {\mathbf k}}\times {\mathrm \overline{{\mathbf U}}_h} 
     886  -g\nabla _h \eta -\frac{c_b^{\textbf U}}{H+\eta} \mathrm {\overline{{\mathbf U}}_h} + \mathrm {\overline{\mathbf G}} 
    867887\end{equation} 
    868888\[ 
    869889  % \label{eq:BT_ssh} 
    870   \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\rm{\bf \overline{U}}}_h \,} \right]+P-E 
     890  \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\mathrm{\mathbf \overline{U}}}_h \,} \right]+P-E 
    871891\] 
    872892% \end{subequations} 
    873 where $\rm {\overline{\bf G}}$ is a forcing term held constant, containing coupling term between modes, 
     893where $\mathrm {\overline{\mathbf G}}$ is a forcing term held constant, containing coupling term between modes, 
    874894surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency. 
    875895The third term on the right hand side of \autoref{eq:BT_dyn} represents the bottom stress 
    876896(see section \autoref{sec:ZDF_bfr}), explicitly accounted for at each barotropic iteration. 
    877897Temporal discretization of the system above follows a three-time step Generalized Forward Backward algorithm 
    878 detailed in \citet{Shchepetkin_McWilliams_OM05}. 
     898detailed in \citet{shchepetkin.mcwilliams_OM05}. 
    879899AB3-AM4 coefficients used in \NEMO follow the second-order accurate, 
    880 "multi-purpose" stability compromise as defined in \citet{Shchepetkin_McWilliams_Bk08} 
     900"multi-purpose" stability compromise as defined in \citet{shchepetkin.mcwilliams_ibk09} 
    881901(see their figure 12, lower left).  
    882902 
     
    884904\begin{figure}[!t] 
    885905  \begin{center} 
    886     \includegraphics[width=0.7\textwidth]{Fig_DYN_dynspg_ts} 
     906    \includegraphics[width=\textwidth]{Fig_DYN_dynspg_ts} 
    887907    \caption{ 
    888908      \protect\label{fig:DYN_dynspg_ts} 
     
    936956and time splitting not compatible. 
    937957Advective barotropic velocities are obtained by using a secondary set of filtering weights, 
    938 uniquely defined from the filter coefficients used for the time averaging (\citet{Shchepetkin_McWilliams_OM05}). 
     958uniquely defined from the filter coefficients used for the time averaging (\citet{shchepetkin.mcwilliams_OM05}). 
    939959Consistency between the time averaged continuity equation and the time stepping of tracers is here the key to 
    940960obtain exact conservation. 
     
    953973external gravity waves in idealized or weakly non-linear cases. 
    954974Although the damping is lower than for the filtered free surface, 
    955 it is still significant as shown by \citet{Levier2007} in the case of an analytical barotropic Kelvin wave. 
     975it is still significant as shown by \citet{levier.treguier.ea_rpt07} in the case of an analytical barotropic Kelvin wave. 
    956976 
    957977%>>>>>=============== 
     
    10511071the leap-frog splitting mode in equation \autoref{eq:DYN_spg_ts_ssh}. 
    10521072We have tried various forms of such filtering, 
    1053 with the following method discussed in \cite{Griffies_al_MWR01} chosen due to 
     1073with the following method discussed in \cite{griffies.pacanowski.ea_MWR01} chosen due to 
    10541074its stability and reasonably good maintenance of tracer conservation properties (see ??). 
    10551075 
     
    10811101% Filtered free surface formulation 
    10821102%-------------------------------------------------------------------------------------------------------------- 
    1083 \subsection{Filtered free surface (\protect\key{dynspg\_flt})} 
     1103\subsection[Filtered free surface (\texttt{\textbf{key\_dynspg\_flt}})] 
     1104{Filtered free surface (\protect\key{dynspg\_flt})} 
    10841105\label{subsec:DYN_spg_fltp} 
    10851106 
    1086 The filtered formulation follows the \citet{Roullet_Madec_JGR00} implementation.  
     1107The filtered formulation follows the \citet{roullet.madec_JGR00} implementation.  
    10871108The extra term introduced in the equations (see \autoref{subsec:PE_free_surface}) is solved implicitly.  
    10881109The elliptic solvers available in the code are documented in \autoref{chap:MISC}. 
     
    10921113  \[ 
    10931114    % \label{eq:spg_flt} 
    1094     \frac{\partial {\rm {\bf U}}_h }{\partial t}= {\rm {\bf M}} 
     1115    \frac{\partial {\mathrm {\mathbf U}}_h }{\partial t}= {\mathrm {\mathbf M}} 
    10951116    - g \nabla \left( \tilde{\rho} \ \eta \right) 
    10961117    - g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right) 
     
    10981119  where $T_c$, is a parameter with dimensions of time which characterizes the force, 
    10991120  $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, 
    1100   and $\rm {\bf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient, 
     1121  and $\mathrm {\mathbf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient, 
    11011122  non-linear and viscous terms in \autoref{eq:PE_dyn}. 
    11021123}   %end gmcomment 
     
    11091130% Lateral diffusion term 
    11101131% ================================================================ 
    1111 \section{Lateral diffusion term and operators (\protect\mdl{dynldf})} 
     1132\section[Lateral diffusion term and operators (\textit{dynldf.F90})] 
     1133{Lateral diffusion term and operators (\protect\mdl{dynldf})} 
    11121134\label{sec:DYN_ldf} 
    11131135%------------------------------------------nam_dynldf---------------------------------------------------- 
     
    11431165 
    11441166% ================================================================ 
    1145 \subsection[Iso-level laplacian (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})] 
    1146             {Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})} 
     1167\subsection[Iso-level laplacian (\forcode{ln_dynldf_lap = .true.})] 
     1168{Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})} 
    11471169\label{subsec:DYN_ldf_lap} 
    11481170 
     
    11521174  \left\{ 
    11531175    \begin{aligned} 
    1154       D_u^{l{\rm {\bf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm} 
     1176      D_u^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm} 
    11551177          \;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta_j \left[  
    11561178        {A_f^{lm} \;e_{3f} \zeta } \right] \\ \\ 
    1157       D_v^{l{\rm {\bf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm} 
     1179      D_v^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm} 
    11581180          \;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta_i \left[  
    11591181        {A_f^{lm} \;e_{3f} \zeta } \right] 
     
    11691191%           Rotated laplacian operator 
    11701192%-------------------------------------------------------------------------------------------------------------- 
    1171 \subsection[Rotated laplacian (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})] 
    1172             {Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})} 
     1193\subsection[Rotated laplacian (\forcode{ln_dynldf_iso = .true.})] 
     1194{Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})} 
    11731195\label{subsec:DYN_ldf_iso} 
    11741196 
     
    12281250%           Iso-level bilaplacian operator 
    12291251%-------------------------------------------------------------------------------------------------------------- 
    1230 \subsection[Iso-level bilaplacian (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})] 
    1231             {Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})} 
     1252\subsection[Iso-level bilaplacian (\forcode{ln_dynldf_bilap = .true.})] 
     1253{Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})} 
    12321254\label{subsec:DYN_ldf_bilap} 
    12331255 
     
    12431265%           Vertical diffusion term 
    12441266% ================================================================ 
    1245 \section{Vertical diffusion term (\protect\mdl{dynzdf})} 
     1267\section[Vertical diffusion term (\textit{dynzdf.F90})] 
     1268{Vertical diffusion term (\protect\mdl{dynzdf})} 
    12461269\label{sec:DYN_zdf} 
    12471270%----------------------------------------------namzdf------------------------------------------------------ 
     
    13261349There are two main options for wetting and drying code (wd): 
    13271350(a) an iterative limiter (il) and (b) a directional limiter (dl). 
    1328 The directional limiter is based on the scheme developed by \cite{WarnerEtal13} for RO 
     1351The directional limiter is based on the scheme developed by \cite{warner.defne.ea_CG13} for RO 
    13291352MS 
    1330 which was in turn based on ideas developed for POM by \cite{Oey06}. The iterative 
     1353which was in turn based on ideas developed for POM by \cite{oey_OM06}. The iterative 
    13311354limiter is a new scheme.  The iterative limiter is activated by setting $\mathrm{ln\_wd\_il} = \mathrm{.true.}$ 
    13321355and $\mathrm{ln\_wd\_dl} = \mathrm{.false.}$. The directional limiter is activated 
     
    13721395%   Iterative limiters 
    13731396%----------------------------------------------------------------------------------------- 
    1374 \subsection   [Directional limiter (\textit{wet\_dry})] 
    1375          {Directional limiter (\mdl{wet\_dry})} 
     1397\subsection[Directional limiter (\textit{wet\_dry.F90})] 
     1398{Directional limiter (\mdl{wet\_dry})} 
    13761399\label{subsec:DYN_wd_directional_limiter} 
    13771400The principal idea of the directional limiter is that 
     
    14001423 
    14011424 
    1402 \cite{WarnerEtal13} state that in their scheme the velocity masks at the cell faces for the baroclinic 
     1425\cite{warner.defne.ea_CG13} state that in their scheme the velocity masks at the cell faces for the baroclinic 
    14031426timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic sub-steps is respectively less than 
    14041427or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer 
     
    14121435%----------------------------------------------------------------------------------------- 
    14131436 
    1414 \subsection   [Iterative limiter (\textit{wet\_dry})] 
    1415          {Iterative limiter (\mdl{wet\_dry})} 
     1437\subsection[Iterative limiter (\textit{wet\_dry.F90})] 
     1438{Iterative limiter (\mdl{wet\_dry})} 
    14161439\label{subsec:DYN_wd_iterative_limiter} 
    14171440 
    1418 \subsubsection [Iterative flux limiter (\textit{wet\_dry})] 
    1419          {Iterative flux limiter (\mdl{wet\_dry})} 
     1441\subsubsection[Iterative flux limiter (\textit{wet\_dry.F90})] 
     1442{Iterative flux limiter (\mdl{wet\_dry})} 
    14201443\label{subsubsec:DYN_wd_il_spg_limiter} 
    14211444 
     
    14941517\end{equation} 
    14951518 
    1496 Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\it [Q: Why is 
     1519Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\itshape [Q: Why is 
    14971520this necessary/desirable?]}. Substituting from (\ref{dyn_wd_continuity_coef}) gives an 
    14981521expression for the coefficient needed to multiply the outward flux at this cell in order 
     
    15221545%      Surface pressure gradients 
    15231546%---------------------------------------------------------------------------------------- 
    1524 \subsubsection   [Modification of surface pressure gradients (\textit{dynhpg})] 
    1525          {Modification of surface pressure gradients (\mdl{dynhpg})} 
     1547\subsubsection[Modification of surface pressure gradients (\textit{dynhpg.F90})] 
     1548{Modification of surface pressure gradients (\mdl{dynhpg})} 
    15261549\label{subsubsec:DYN_wd_il_spg} 
    15271550 
     
    15411564%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    15421565\begin{figure}[!ht] \begin{center} 
    1543 \includegraphics[width=0.8\textwidth]{Fig_WAD_dynhpg} 
     1566\includegraphics[width=\textwidth]{Fig_WAD_dynhpg} 
    15441567\caption{ \label{Fig_WAD_dynhpg} 
    15451568Illustrations of the three possible combinations of the logical variables controlling the 
     
    15881611conditions. 
    15891612 
    1590 \subsubsection   [Additional considerations (\textit{usrdef\_zgr})] 
    1591          {Additional considerations (\mdl{usrdef\_zgr})} 
     1613\subsubsection[Additional considerations (\textit{usrdef\_zgr.F90})] 
     1614{Additional considerations (\mdl{usrdef\_zgr})} 
    15921615\label{subsubsec:WAD_additional} 
    15931616 
     
    16031626%      The WAD test cases 
    16041627%---------------------------------------------------------------------------------------- 
    1605 \subsection   [The WAD test cases (\textit{usrdef\_zgr})] 
    1606          {The WAD test cases (\mdl{usrdef\_zgr})} 
     1628\subsection[The WAD test cases (\textit{usrdef\_zgr.F90})] 
     1629{The WAD test cases (\mdl{usrdef\_zgr})} 
    16071630\label{WAD_test_cases} 
    16081631 
     
    16141637% Time evolution term  
    16151638% ================================================================ 
    1616 \section{Time evolution term (\protect\mdl{dynnxt})} 
     1639\section[Time evolution term (\textit{dynnxt.F90})] 
     1640{Time evolution term (\protect\mdl{dynnxt})} 
    16171641\label{sec:DYN_nxt} 
    16181642 
  • NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_LBC.tex

    r10614 r11422  
    1717% Boundary Condition at the Coast 
    1818% ================================================================ 
    19 \section{Boundary condition at the coast (\protect\np{rn\_shlat})} 
     19\section[Boundary condition at the coast (\texttt{rn\_shlat})] 
     20{Boundary condition at the coast (\protect\np{rn\_shlat})} 
    2021\label{sec:LBC_coast} 
    2122%--------------------------------------------nam_lbc------------------------------------------------------- 
     
    5657\begin{figure}[!t] 
    5758  \begin{center} 
    58     \includegraphics[width=0.90\textwidth]{Fig_LBC_uv} 
     59    \includegraphics[width=\textwidth]{Fig_LBC_uv} 
    5960    \caption{ 
    6061      \protect\label{fig:LBC_uv} 
     
    8586\begin{figure}[!p] 
    8687  \begin{center} 
    87     \includegraphics[width=0.90\textwidth]{Fig_LBC_shlat} 
     88    \includegraphics[width=\textwidth]{Fig_LBC_shlat} 
    8889    \caption{ 
    8990      \protect\label{fig:LBC_shlat} 
     
    147148% Boundary Condition around the Model Domain 
    148149% ================================================================ 
    149 \section{Model domain boundary condition (\protect\np{jperio})} 
     150\section[Model domain boundary condition (\texttt{jperio})] 
     151{Model domain boundary condition (\protect\np{jperio})} 
    150152\label{sec:LBC_jperio} 
    151153 
     
    158160%        Closed, cyclic (\np{jperio}\forcode{ = 0..2})  
    159161% ------------------------------------------------------------------------------------------------------------- 
    160 \subsection{Closed, cyclic (\protect\np{jperio}\forcode{= [0127]})} 
     162\subsection[Closed, cyclic (\forcode{jperio = [0127]})] 
     163{Closed, cyclic (\protect\np{jperio}\forcode{ = [0127]})} 
    161164\label{subsec:LBC_jperio012} 
    162165 
     
    194197\begin{figure}[!t] 
    195198  \begin{center} 
    196     \includegraphics[width=1.0\textwidth]{Fig_LBC_jperio} 
     199    \includegraphics[width=\textwidth]{Fig_LBC_jperio} 
    197200    \caption{ 
    198201      \protect\label{fig:LBC_jperio} 
     
    206209%        North fold (\textit{jperio = 3 }to $6)$  
    207210% ------------------------------------------------------------------------------------------------------------- 
    208 \subsection{North-fold (\protect\np{jperio}\forcode{ = 3..6})} 
     211\subsection[North-fold (\forcode{jperio = [3-6]})] 
     212{North-fold (\protect\np{jperio}\forcode{ = [3-6]})} 
    209213\label{subsec:LBC_north_fold} 
    210214 
     
    218222\begin{figure}[!t] 
    219223  \begin{center} 
    220     \includegraphics[width=0.90\textwidth]{Fig_North_Fold_T} 
     224    \includegraphics[width=\textwidth]{Fig_North_Fold_T} 
    221225    \caption{ 
    222226      \protect\label{fig:North_Fold_T} 
     
    232236% Exchange with neighbouring processors  
    233237% ==================================================================== 
    234 \section{Exchange with neighbouring processors (\protect\mdl{lbclnk}, \protect\mdl{lib\_mpp})} 
     238\section[Exchange with neighbouring processors (\textit{lbclnk.F90}, \textit{lib\_mpp.F90})] 
     239{Exchange with neighbouring processors (\protect\mdl{lbclnk}, \protect\mdl{lib\_mpp})} 
    235240\label{sec:LBC_mpp} 
    236241 
     
    280285\begin{figure}[!t] 
    281286  \begin{center} 
    282     \includegraphics[width=0.90\textwidth]{Fig_mpp} 
     287    \includegraphics[width=\textwidth]{Fig_mpp} 
    283288    \caption{ 
    284289      \protect\label{fig:mpp} 
     
    360365\begin{figure}[!ht] 
    361366  \begin{center} 
    362     \includegraphics[width=0.90\textwidth]{Fig_mppini2} 
     367    \includegraphics[width=\textwidth]{Fig_mppini2} 
    363368    \caption { 
    364369      \protect\label{fig:mppini2} 
     
    395400 
    396401The BDY module was modelled on the OBC module (see NEMO 3.4) and shares many features and 
    397 a similar coding structure \citep{Chanut2005}. 
     402a similar coding structure \citep{chanut_rpt05}. 
    398403The specification of the location of the open boundary is completely flexible and 
    399404allows for example the open boundary to follow an isobath or other irregular contour.  
     
    475480\label{subsec:BDY_FRS_scheme} 
    476481 
    477 The Flow Relaxation Scheme (FRS) \citep{Davies_QJRMS76,Engerdahl_Tel95}, 
     482The Flow Relaxation Scheme (FRS) \citep{davies_QJRMS76,engedahl_T95}, 
    478483applies a simple relaxation of the model fields to externally-specified values over 
    479484a zone next to the edge of the model domain. 
     
    514519\label{subsec:BDY_flather_scheme} 
    515520 
    516 The \citet{Flather_JPO94} scheme is a radiation condition on the normal, 
     521The \citet{flather_JPO94} scheme is a radiation condition on the normal, 
    517522depth-mean transport across the open boundary. 
    518523It takes the form 
     
    535540\label{subsec:BDY_orlanski_scheme} 
    536541 
    537 The Orlanski scheme is based on the algorithm described by \citep{Marchesiello2001}, hereafter MMS. 
     542The Orlanski scheme is based on the algorithm described by \citep{marchesiello.mcwilliams.ea_OM01}, hereafter MMS. 
    538543 
    539544The adaptive Orlanski condition solves a wave plus relaxation equation at the boundary: 
     
    636641\begin{figure}[!t] 
    637642  \begin{center} 
    638     \includegraphics[width=1.0\textwidth]{Fig_LBC_bdy_geom} 
     643    \includegraphics[width=\textwidth]{Fig_LBC_bdy_geom} 
    639644    \caption { 
    640645      \protect\label{fig:LBC_bdy_geom} 
     
    670675These restrictions mean that data files used with versions of the 
    671676model prior to Version 3.4 may not work with Version 3.4 onwards. 
    672 A \fortran utility {\it bdy\_reorder} exists in the TOOLS directory which 
     677A \fortran utility {\itshape bdy\_reorder} exists in the TOOLS directory which 
    673678will re-order the data in old BDY data files. 
    674679 
     
    676681\begin{figure}[!t] 
    677682  \begin{center} 
    678     \includegraphics[width=1.0\textwidth]{Fig_LBC_nc_header} 
     683    \includegraphics[width=\textwidth]{Fig_LBC_nc_header} 
    679684    \caption { 
    680685      \protect\label{fig:LBC_nc_header} 
  • NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_LDF.tex

    r10442 r11422  
    2323These three aspects of the lateral diffusion are set through namelist parameters 
    2424(see the \ngn{nam\_traldf} and \ngn{nam\_dynldf} below). 
    25 Note that this chapter describes the standard implementation of iso-neutral tracer mixing, 
    26 and Griffies's implementation, which is used if \np{traldf\_grif}\forcode{ = .true.}, 
    27 is described in Appdx\autoref{apdx:triad} 
    28  
    29 %-----------------------------------nam_traldf - nam_dynldf-------------------------------------------- 
     25Note that this chapter describes the standard implementation of iso-neutral tracer mixing.  
     26Griffies's implementation, which is used if \np{ln\_traldf\_triad}\forcode{ = .true.}, 
     27is described in \autoref{apdx:triad} 
     28 
     29%-----------------------------------namtra_ldf - namdyn_ldf-------------------------------------------- 
    3030 
    3131\nlst{namtra_ldf}  
     
    3434%-------------------------------------------------------------------------------------------------------------- 
    3535 
     36% ================================================================ 
     37% Lateral Mixing Operator 
     38% ================================================================ 
     39\section[Lateral mixing operators] 
     40{Lateral mixing operators} 
     41\label{sec:LDF_op} 
     42We remind here the different lateral mixing operators that can be used. Further details can be found in \autoref{subsec:TRA_ldf_op} and  \autoref{sec:DYN_ldf}. 
     43 
     44\subsection[No lateral mixing (\forcode{ln_traldf_OFF}, \forcode{ln_dynldf_OFF})] 
     45{No lateral mixing (\protect\np{ln\_traldf\_OFF}, \protect\np{ln\_dynldf\_OFF})} 
     46 
     47It is possible to run without explicit lateral diffusion on tracers (\protect\np{ln\_traldf\_OFF}\forcode{ = .true.}) and/or  
     48momentum (\protect\np{ln\_dynldf\_OFF}\forcode{ = .true.}). The latter option is even recommended if using the  
     49UBS advection scheme on momentum (\np{ln\_dynadv\_ubs}\forcode{ = .true.}, 
     50see \autoref{subsec:DYN_adv_ubs}) and can be useful for testing purposes. 
     51 
     52\subsection[Laplacian mixing (\forcode{ln_traldf_lap}, \forcode{ln_dynldf_lap})] 
     53{Laplacian mixing (\protect\np{ln\_traldf\_lap}, \protect\np{ln\_dynldf\_lap})} 
     54Setting \protect\np{ln\_traldf\_lap}\forcode{ = .true.} and/or \protect\np{ln\_dynldf\_lap}\forcode{ = .true.} enables  
     55a second order diffusion on tracers and momentum respectively. Note that in \NEMO 4, one can not combine  
     56Laplacian and Bilaplacian operators for the same variable. 
     57 
     58\subsection[Bilaplacian mixing (\forcode{ln_traldf_blp}, \forcode{ln_dynldf_blp})] 
     59{Bilaplacian mixing (\protect\np{ln\_traldf\_blp}, \protect\np{ln\_dynldf\_blp})} 
     60Setting \protect\np{ln\_traldf\_blp}\forcode{ = .true.} and/or \protect\np{ln\_dynldf\_blp}\forcode{ = .true.} enables  
     61a fourth order diffusion on tracers and momentum respectively. It is implemented by calling the above Laplacian operator twice.  
     62We stress again that from \NEMO 4, the simultaneous use Laplacian and Bilaplacian operators is not allowed. 
    3663 
    3764% ================================================================ 
    3865% Direction of lateral Mixing 
    3966% ================================================================ 
    40 \section{Direction of lateral mixing (\protect\mdl{ldfslp})} 
     67\section[Direction of lateral mixing (\textit{ldfslp.F90})] 
     68{Direction of lateral mixing (\protect\mdl{ldfslp})} 
    4169\label{sec:LDF_slp} 
    4270 
     
    4472\gmcomment{ 
    4573  we should emphasize here that the implementation is a rather old one. 
    46   Better work can be achieved by using \citet{Griffies_al_JPO98, Griffies_Bk04} iso-neutral scheme. 
     74  Better work can be achieved by using \citet{griffies.gnanadesikan.ea_JPO98, griffies_bk04} iso-neutral scheme. 
    4775} 
    4876 
     
    87115%gm%  caution I'm not sure the simplification was a good idea!  
    88116 
    89 These slopes are computed once in \rou{ldfslp\_init} when \np{ln\_sco}\forcode{ = .true.}rue, 
     117These slopes are computed once in \rou{ldf\_slp\_init} when \np{ln\_sco}\forcode{ = .true.}, 
    90118and either \np{ln\_traldf\_hor}\forcode{ = .true.} or \np{ln\_dynldf\_hor}\forcode{ = .true.}.  
    91119 
     
    119147%In practice, \autoref{eq:ldfslp_iso} is of little help in evaluating the neutral surface slopes. Indeed, for an unsimplified equation of state, the density has a strong dependancy on pressure (here approximated as the depth), therefore applying \autoref{eq:ldfslp_iso} using the $in situ$ density, $\rho$, computed at T-points leads to a flattening of slopes as the depth increases. This is due to the strong increase of the $in situ$ density with depth.  
    120148 
    121 %By definition, neutral surfaces are tangent to the local $in situ$ density \citep{McDougall1987}, therefore in \autoref{eq:ldfslp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters). 
     149%By definition, neutral surfaces are tangent to the local $in situ$ density \citep{mcdougall_JPO87}, therefore in \autoref{eq:ldfslp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters). 
    122150 
    123151%In the $z$-coordinate, the derivative of the  \autoref{eq:ldfslp_iso} numerator is evaluated at the same depth \nocite{as what?} ($T$-level, which is the same as the $u$- and $v$-levels), so  the $in situ$ density can be used for its evaluation.  
     
    135163  thus the $in situ$ density can be used. 
    136164  This is not the case for the vertical derivatives: $\delta_{k+1/2}[\rho]$ is replaced by $-\rho N^2/g$, 
    137   where $N^2$ is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following \citet{McDougall1987} 
     165  where $N^2$ is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following \citet{mcdougall_JPO87} 
    138166  (see \autoref{subsec:TRA_bn2}).  
    139167 
     
    144172\item[$s$- or hybrid $s$-$z$- coordinate: ] 
    145173  in the current release of \NEMO, iso-neutral mixing is only employed for $s$-coordinates if 
    146   the Griffies scheme is used (\np{traldf\_grif}\forcode{ = .true.}; 
    147   see Appdx \autoref{apdx:triad}). 
     174  the Griffies scheme is used (\np{ln\_traldf\_triad}\forcode{ = .true.}; 
     175  see \autoref{apdx:triad}). 
    148176  In other words, iso-neutral mixing will only be accurately represented with a linear equation of state 
    149   (\np{nn\_eos}\forcode{ = 1..2}). 
     177  (\np{ln\_seos}\forcode{ = .true.}). 
    150178  In the case of a "true" equation of state, the evaluation of $i$ and $j$ derivatives in \autoref{eq:ldfslp_iso} 
    151179  will include a pressure dependent part, leading to the wrong evaluation of the neutral slopes. 
     
    154182  Note: The solution for $s$-coordinate passes trough the use of different (and better) expression for 
    155183  the constraint on iso-neutral fluxes. 
    156   Following \citet{Griffies_Bk04}, instead of specifying directly that there is a zero neutral diffusive flux of 
     184  Following \citet{griffies_bk04}, instead of specifying directly that there is a zero neutral diffusive flux of 
    157185  locally referenced potential density, we stay in the $T$-$S$ plane and consider the balance between 
    158186  the neutral direction diffusive fluxes of potential temperature and salinity: 
     
    197225 
    198226This implementation is a rather old one. 
    199 It is similar to the one proposed by Cox [1987], except for the background horizontal diffusion. 
    200 Indeed, the Cox implementation of isopycnal diffusion in GFDL-type models requires 
     227It is similar to the one proposed by \citet{cox_OM87}, except for the background horizontal diffusion. 
     228Indeed, the \citet{cox_OM87} implementation of isopycnal diffusion in GFDL-type models requires 
    201229a minimum background horizontal diffusion for numerical stability reasons. 
    202230To overcome this problem, several techniques have been proposed in which the numerical schemes of 
    203 the ocean model are modified \citep{Weaver_Eby_JPO97, Griffies_al_JPO98}. 
    204 Griffies's scheme is now available in \NEMO if \np{traldf\_grif\_iso} is set true; see Appdx \autoref{apdx:triad}. 
    205 Here, another strategy is presented \citep{Lazar_PhD97}: 
     231the ocean model are modified \citep{weaver.eby_JPO97, griffies.gnanadesikan.ea_JPO98}. 
     232Griffies's scheme is now available in \NEMO if \np{ln\_traldf\_triad}=\forcode{= .true.}; see \autoref{apdx:triad}. 
     233Here, another strategy is presented \citep{lazar_phd97}: 
    206234a local filtering of the iso-neutral slopes (made on 9 grid-points) prevents the development of 
    207235grid point noise generated by the iso-neutral diffusion operator (\autoref{fig:LDF_ZDF1}). 
     
    212240 
    213241Nevertheless, this iso-neutral operator does not ensure that variance cannot increase, 
    214 contrary to the \citet{Griffies_al_JPO98} operator which has that property.  
     242contrary to the \citet{griffies.gnanadesikan.ea_JPO98} operator which has that property.  
    215243 
    216244%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    217245\begin{figure}[!ht] 
    218246  \begin{center} 
    219     \includegraphics[width=0.70\textwidth]{Fig_LDF_ZDF1} 
     247    \includegraphics[width=\textwidth]{Fig_LDF_ZDF1} 
    220248    \caption { 
    221249      \protect\label{fig:LDF_ZDF1} 
     
    235263 
    236264 
    237 % In addition and also for numerical stability reasons \citep{Cox1987, Griffies_Bk04},  
     265% In addition and also for numerical stability reasons \citep{cox_OM87, griffies_bk04},  
    238266% the slopes are bounded by $1/100$ everywhere. This limit is decreasing linearly  
    239267% to zero fom $70$ meters depth and the surface (the fact that the eddies "feel" the  
    240268% surface motivates this flattening of isopycnals near the surface). 
    241269 
    242 For numerical stability reasons \citep{Cox1987, Griffies_Bk04}, the slopes must also be bounded by 
    243 $1/100$ everywhere. 
     270For numerical stability reasons \citep{cox_OM87, griffies_bk04}, the slopes must also be bounded by 
     271the namelist scalar \np{rn\_slpmax} (usually $1/100$) everywhere. 
    244272This constraint is applied in a piecewise linear fashion, increasing from zero at the surface to 
    245273$1/100$ at $70$ metres and thereafter decreasing to zero at the bottom of the ocean 
    246274(the fact that the eddies "feel" the surface motivates this flattening of isopycnals near the surface). 
     275\colorbox{yellow}{The way slopes are tapered has be checked. Not sure that this is still what is actually done.} 
    247276 
    248277%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    249278\begin{figure}[!ht] 
    250279  \begin{center} 
    251     \includegraphics[width=0.70\textwidth]{Fig_eiv_slp} 
     280    \includegraphics[width=\textwidth]{Fig_eiv_slp} 
    252281    \caption{ 
    253282      \protect\label{fig:eiv_slp} 
     
    297326(see \autoref{sec:LBC_coast}). 
    298327 
    299  
    300 % ================================================================ 
    301 % Lateral Mixing Operator 
    302 % ================================================================ 
    303 \section{Lateral mixing operators (\protect\mdl{traldf}, \protect\mdl{traldf}) } 
    304 \label{sec:LDF_op} 
    305  
    306  
    307328    
    308329% ================================================================ 
    309330% Lateral Mixing Coefficients 
    310331% ================================================================ 
    311 \section{Lateral mixing coefficient (\protect\mdl{ldftra}, \protect\mdl{ldfdyn}) } 
     332\section[Lateral mixing coefficient (\forcode{nn_aht_ijk_t}, \forcode{nn_ahm_ijk_t})] 
     333{Lateral mixing coefficient (\protect\np{nn\_aht\_ijk\_t}, \protect\np{nn\_ahm\_ijk\_t})} 
    312334\label{sec:LDF_coef} 
    313335 
    314 Introducing a space variation in the lateral eddy mixing coefficients changes the model core memory requirement, 
    315 adding up to four extra three-dimensional arrays for the geopotential or isopycnal second order operator applied to  
    316 momentum. 
    317 Six CPP keys control the space variation of eddy coefficients: three for momentum and three for tracer. 
    318 The three choices allow: 
    319 a space variation in the three space directions (\key{traldf\_c3d},  \key{dynldf\_c3d}), 
    320 in the horizontal plane (\key{traldf\_c2d},  \key{dynldf\_c2d}), 
    321 or in the vertical only (\key{traldf\_c1d},  \key{dynldf\_c1d}). 
    322 The default option is a constant value over the whole ocean on both momentum and tracers.  
    323     
    324 The number of additional arrays that have to be defined and the gridpoint position at which 
    325 they are defined depend on both the space variation chosen and the type of operator used. 
    326 The resulting eddy viscosity and diffusivity coefficients can be a function of more than one variable. 
    327 Changes in the computer code when switching from one option to another have been minimized by 
    328 introducing the eddy coefficients as statement functions 
    329 (include file \textit{ldftra\_substitute.h90} and \textit{ldfdyn\_substitute.h90}). 
    330 The functions are replaced by their actual meaning during the preprocessing step (CPP). 
    331 The specification of the space variation of the coefficient is made in \mdl{ldftra} and \mdl{ldfdyn}, 
    332 or more precisely in include files \textit{traldf\_cNd.h90} and \textit{dynldf\_cNd.h90}, with N=1, 2 or 3. 
    333 The user can modify these include files as he/she wishes. 
    334 The way the mixing coefficient are set in the reference version can be briefly described as follows: 
    335  
    336 \subsubsection{Constant mixing coefficients (default option)} 
    337 When none of the \key{dynldf\_...} and \key{traldf\_...} keys are defined, 
    338 a constant value is used over the whole ocean for momentum and tracers, 
    339 which is specified through the \np{rn\_ahm0} and \np{rn\_aht0} namelist parameters. 
    340  
    341 \subsubsection{Vertically varying mixing coefficients (\protect\key{traldf\_c1d} and \key{dynldf\_c1d})}  
    342 The 1D option is only available when using the $z$-coordinate with full step. 
    343 Indeed in all the other types of vertical coordinate, 
    344 the depth is a 3D function of (\textbf{i},\textbf{j},\textbf{k}) and therefore, 
    345 introducing depth-dependent mixing coefficients will require 3D arrays. 
    346 In the 1D option, a hyperbolic variation of the lateral mixing coefficient is introduced in which 
    347 the surface value is \np{rn\_aht0} (\np{rn\_ahm0}), the bottom value is 1/4 of the surface value, 
    348 and the transition takes place around z=300~m with a width of 300~m 
    349 (\ie both the depth and the width of the inflection point are set to 300~m). 
    350 This profile is hard coded in file \textit{traldf\_c1d.h90}, but can be easily modified by users. 
    351  
    352 \subsubsection{Horizontally varying mixing coefficients (\protect\key{traldf\_c2d} and \protect\key{dynldf\_c2d})} 
    353 By default the horizontal variation of the eddy coefficient depends on the local mesh size and 
     336The specification of the space variation of the coefficient is made in modules \mdl{ldftra} and \mdl{ldfdyn}.  
     337The way the mixing coefficients are set in the reference version can be described as follows: 
     338 
     339\subsection[Mixing coefficients read from file (\forcode{nn_aht_ijk_t = -20, -30}, \forcode{nn_ahm_ijk_t = -20,-30})] 
     340{ Mixing coefficients read from file (\protect\np{nn\_aht\_ijk\_t}\forcode{ = -20, -30}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = -20, -30})} 
     341 
     342Mixing coefficients can be read from file if a particular geographical variation is needed. For example, in the ORCA2 global ocean model,  
     343the laplacian viscosity operator uses $A^l$~= 4.10$^4$ m$^2$/s poleward of 20$^{\circ}$ north and south and 
     344decreases linearly to $A^l$~= 2.10$^3$ m$^2$/s at the equator \citep{madec.delecluse.ea_JPO96, delecluse.madec_icol99}.  
     345Similar modified horizontal variations can be found with the Antarctic or Arctic sub-domain options of ORCA2 and ORCA05.  
     346The provided fields can either be 2d (\np{nn\_aht\_ijk\_t}\forcode{ = -20}, \np{nn\_ahm\_ijk\_t}\forcode{ = -20}) or 3d (\np{nn\_aht\_ijk\_t}\forcode{ = -30},  \np{nn\_ahm\_ijk\_t}\forcode{ = -30}). They must be given at U, V points for tracers and T, F points for momentum (see \autoref{tab:LDF_files}). 
     347 
     348%-------------------------------------------------TABLE--------------------------------------------------- 
     349\begin{table}[htb] 
     350  \begin{center} 
     351    \begin{tabular}{|l|l|l|l|} 
     352      \hline 
     353      Namelist parameter                        & Input filename                               & dimensions & variable names                      \\  \hline 
     354      \np{nn\_ahm\_ijk\_t}\forcode{ = -20}       & \forcode{eddy_viscosity_2D.nc }            &     $(i,j)$         & \forcode{ahmt_2d, ahmf_2d}  \\  \hline 
     355      \np{nn\_aht\_ijk\_t}\forcode{ = -20}           & \forcode{eddy_diffusivity_2D.nc }           &     $(i,j)$          & \forcode{ahtu_2d, ahtv_2d}    \\   \hline 
     356      \np{nn\_ahm\_ijk\_t}\forcode{ = -30}       & \forcode{eddy_viscosity_3D.nc }            &     $(i,j,k)$          & \forcode{ahmt_3d, ahmf_3d}  \\  \hline 
     357      \np{nn\_aht\_ijk\_t}\forcode{ = -30}       & \forcode{eddy_diffusivity_3D.nc }           &     $(i,j,k)$         & \forcode{ahtu_3d, ahtv_3d}    \\   \hline 
     358    \end{tabular} 
     359    \caption{ 
     360      \protect\label{tab:LDF_files} 
     361      Description of expected input files if mixing coefficients are read from NetCDF files. 
     362    } 
     363  \end{center} 
     364\end{table} 
     365%-------------------------------------------------------------------------------------------------------------- 
     366 
     367\subsection[Constant mixing coefficients (\forcode{nn_aht_ijk_t = 0}, \forcode{nn_ahm_ijk_t = 0})] 
     368{ Constant mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 0}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 0})} 
     369 
     370If constant, mixing coefficients are set thanks to a velocity and a length scales ($U_{scl}$, $L_{scl}$) such that: 
     371 
     372\begin{equation} 
     373  \label{eq:constantah} 
     374  A_o^l = \left\{ 
     375    \begin{aligned} 
     376      & \frac{1}{2} U_{scl} L_{scl}            & \text{for laplacian operator } \\ 
     377      & \frac{1}{12} U_{scl} L_{scl}^3                    & \text{for bilaplacian operator } 
     378    \end{aligned} 
     379  \right. 
     380\end{equation} 
     381 
     382 $U_{scl}$ and $L_{scl}$ are given by the namelist parameters \np{rn\_Ud}, \np{rn\_Uv}, \np{rn\_Ld} and \np{rn\_Lv}. 
     383 
     384\subsection[Vertically varying mixing coefficients (\forcode{nn_aht_ijk_t = 10}, \forcode{nn_ahm_ijk_t = 10})] 
     385{Vertically varying mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 10}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 10})} 
     386 
     387In the vertically varying case, a hyperbolic variation of the lateral mixing coefficient is introduced in which 
     388the surface value is given by \autoref{eq:constantah}, the bottom value is 1/4 of the surface value, 
     389and the transition takes place around z=500~m with a width of 200~m. 
     390This profile is hard coded in module \mdl{ldfc1d\_c2d}, but can be easily modified by users. 
     391 
     392\subsection[Mesh size dependent mixing coefficients (\forcode{nn_aht_ijk_t = 20}, \forcode{nn_ahm_ijk_t = 20})] 
     393{Mesh size dependent mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 20}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 20})} 
     394 
     395In that case, the horizontal variation of the eddy coefficient depends on the local mesh size and 
    354396the type of operator used: 
    355397\begin{equation} 
     
    357399  A_l = \left\{ 
    358400    \begin{aligned} 
    359       & \frac{\max(e_1,e_2)}{e_{max}} A_o^l           & \text{for laplacian operator } \\ 
    360       & \frac{\max(e_1,e_2)^{3}}{e_{max}^{3}} A_o^l          & \text{for bilaplacian operator } 
     401      & \frac{1}{2} U_{scl}  \max(e_1,e_2)         & \text{for laplacian operator } \\ 
     402      & \frac{1}{12} U_{scl}  \max(e_1,e_2)^{3}             & \text{for bilaplacian operator } 
    361403    \end{aligned} 
    362404  \right. 
    363405\end{equation} 
    364 where $e_{max}$ is the maximum of $e_1$ and $e_2$ taken over the whole masked ocean domain, 
    365 and $A_o^l$ is the \np{rn\_ahm0} (momentum) or \np{rn\_aht0} (tracer) namelist parameter. 
     406where $U_{scl}$ is a user defined velocity scale (\np{rn\_Ud}, \np{rn\_Uv}). 
    366407This variation is intended to reflect the lesser need for subgrid scale eddy mixing where 
    367408the grid size is smaller in the domain. 
    368 It was introduced in the context of the DYNAMO modelling project \citep{Willebrand_al_PO01}. 
    369 Note that such a grid scale dependance of mixing coefficients significantly increase the range of stability of 
    370 model configurations presenting large changes in grid pacing such as global ocean models. 
     409It was introduced in the context of the DYNAMO modelling project \citep{willebrand.barnier.ea_PO01}. 
     410Note that such a grid scale dependance of mixing coefficients significantly increases the range of stability of 
     411model configurations presenting large changes in grid spacing such as global ocean models. 
    371412Indeed, in such a case, a constant mixing coefficient can lead to a blow up of the model due to 
    372413large coefficient compare to the smallest grid size (see \autoref{sec:STP_forward_imp}), 
    373414especially when using a bilaplacian operator. 
    374415 
    375 Other formulations can be introduced by the user for a given configuration. 
    376 For example, in the ORCA2 global ocean model (see Configurations), 
    377 the laplacian viscosity operator uses \np{rn\_ahm0}~= 4.10$^4$ m$^2$/s poleward of 20$^{\circ}$ north and south and 
    378 decreases linearly to \np{rn\_aht0}~= 2.10$^3$ m$^2$/s at the equator \citep{Madec_al_JPO96, Delecluse_Madec_Bk00}. 
    379 This modification can be found in routine \rou{ldf\_dyn\_c2d\_orca} defined in \mdl{ldfdyn\_c2d}. 
    380 Similar modified horizontal variations can be found with the Antarctic or Arctic sub-domain options of 
    381 ORCA2 and ORCA05 (see \&namcfg namelist). 
    382  
    383 \subsubsection{Space varying mixing coefficients (\protect\key{traldf\_c3d} and \key{dynldf\_c3d})} 
    384  
    385 The 3D space variation of the mixing coefficient is simply the combination of the 1D and 2D cases, 
     416\colorbox{yellow}{CASE \np{nn\_aht\_ijk\_t} = 21 to be added} 
     417 
     418\subsection[Mesh size and depth dependent mixing coefficients (\forcode{nn_aht_ijk_t = 30}, \forcode{nn_ahm_ijk_t = 30})] 
     419{Mesh size and depth dependent mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 30}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 30})} 
     420 
     421The 3D space variation of the mixing coefficient is simply the combination of the 1D and 2D cases above, 
    386422\ie a hyperbolic tangent variation with depth associated with a grid size dependence of 
    387423the magnitude of the coefficient.  
    388424 
    389 \subsubsection{Space and time varying mixing coefficients} 
    390  
    391 There is no default specification of space and time varying mixing coefficient.  
    392 The only case available is specific to the ORCA2 and ORCA05 global ocean configurations. 
    393 It provides only a tracer mixing coefficient for eddy induced velocity (ORCA2) or both iso-neutral and 
    394 eddy induced velocity (ORCA05) that depends on the local growth rate of baroclinic instability. 
    395 This specification is actually used when an ORCA key and both \key{traldf\_eiv} and \key{traldf\_c2d} are defined. 
     425\subsection[Velocity dependent mixing coefficients (\forcode{nn_aht_ijk_t = 31}, \forcode{nn_ahm_ijk_t = 31})] 
     426{Flow dependent mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 31}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 31})} 
     427In that case, the eddy coefficient is proportional to the local velocity magnitude so that the Reynolds number $Re =  \lvert U \rvert  e / A_l$ is constant (and here hardcoded to $12$): 
     428\colorbox{yellow}{JC comment: The Reynolds is effectively set to 12 in the code for both operators but shouldn't it be 2 for Laplacian ?} 
     429 
     430\begin{equation} 
     431  \label{eq:flowah} 
     432  A_l = \left\{ 
     433    \begin{aligned} 
     434      & \frac{1}{12} \lvert U \rvert e          & \text{for laplacian operator } \\ 
     435      & \frac{1}{12} \lvert U \rvert e^3             & \text{for bilaplacian operator }  
     436    \end{aligned} 
     437  \right. 
     438\end{equation} 
     439 
     440\subsection[Deformation rate dependent viscosities (\forcode{nn_ahm_ijk_t = 32})] 
     441{Deformation rate dependent viscosities (\protect\np{nn\_ahm\_ijk\_t}\forcode{ = 32})} 
     442 
     443This option refers to the \citep{smagorinsky_MW63} scheme which is here implemented for momentum only. Smagorinsky chose as a  
     444characteristic time scale $T_{smag}$ the deformation rate and for the lengthscale $L_{smag}$ the maximum wavenumber possible on the horizontal grid, e.g.: 
     445 
     446\begin{equation} 
     447  \label{eq:smag1} 
     448  \begin{split} 
     449    T_{smag}^{-1} & = \sqrt{\left( \partial_x u - \partial_y v\right)^2 + \left( \partial_y u + \partial_x v\right)^2  } \\ 
     450    L_{smag} & = \frac{1}{\pi}\frac{e_1 e_2}{e_1 + e_2} 
     451  \end{split} 
     452\end{equation} 
     453 
     454Introducing a user defined constant $C$ (given in the namelist as \np{rn\_csmc}), one can deduce the mixing coefficients as follows: 
     455 
     456\begin{equation} 
     457  \label{eq:smag2} 
     458  A_{smag} = \left\{ 
     459    \begin{aligned} 
     460      & C^2  T_{smag}^{-1}  L_{smag}^2       & \text{for laplacian operator } \\ 
     461      & \frac{C^2}{8} T_{smag}^{-1} L_{smag}^4        & \text{for bilaplacian operator }  
     462    \end{aligned} 
     463  \right. 
     464\end{equation} 
     465 
     466For stability reasons, upper and lower limits are applied on the resulting coefficient (see \autoref{sec:STP_forward_imp}) so that: 
     467\begin{equation} 
     468  \label{eq:smag3} 
     469    \begin{aligned} 
     470      & C_{min} \frac{1}{2}   \lvert U \rvert  e    < A_{smag} < C_{max} \frac{e^2}{   8\rdt}                 & \text{for laplacian operator } \\ 
     471      & C_{min} \frac{1}{12} \lvert U \rvert  e^3 < A_{smag} < C_{max} \frac{e^4}{64 \rdt}                 & \text{for bilaplacian operator }  
     472    \end{aligned} 
     473\end{equation} 
     474 
     475where $C_{min}$ and $C_{max}$ are adimensional namelist parameters given by \np{rn\_minfac} and \np{rn\_maxfac} respectively. 
     476 
     477\subsection{About space and time varying mixing coefficients} 
    396478 
    397479The following points are relevant when the eddy coefficient varies spatially: 
     
    406488(\autoref{sec:dynldf_properties}). 
    407489 
    408 (3) for isopycnal diffusion on momentum or tracers, an additional purely horizontal background diffusion with 
    409 uniform coefficient can be added by setting a non zero value of \np{rn\_ahmb0} or \np{rn\_ahtb0}, 
    410 a background horizontal eddy viscosity or diffusivity coefficient 
    411 (namelist parameters whose default values are $0$). 
    412 However, the technique used to compute the isopycnal slopes is intended to get rid of such a background diffusion, 
    413 since it introduces spurious diapycnal diffusion (see \autoref{sec:LDF_slp}). 
    414  
    415 (4) when an eddy induced advection term is used (\key{traldf\_eiv}), 
    416 $A^{eiv}$, the eddy induced coefficient has to be defined. 
    417 Its space variations are controlled by the same CPP variable as for the eddy diffusivity coefficient 
    418 (\ie \key{traldf\_cNd}).  
    419  
    420 (5) the eddy coefficient associated with a biharm