New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 1831 for branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_DYN.tex – NEMO

Ignore:
Timestamp:
2010-04-12T16:59:59+02:00 (14 years ago)
Author:
gm
Message:

cover, namelist, rigid-lid, e3t, appendices, see ticket: #658

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_DYN.tex

    r1224 r1831  
    2424         + \text{HPG} + \text{SPG} + \text{LDF} + \text{ZDF} 
    2525\end{equation*} 
    26  
    2726NXT stands for next, referring to the time-stepping. The first group of terms on  
    2827the rhs of the momentum equations corresponds to the Coriolis and advection  
    2928terms that are decomposed into a vorticity part (VOR), a kinetic energy part (KEG)  
    3029and, a vertical advection part (ZAD) in the vector invariant formulation or a Coriolis  
    31 and advection part(COR+ADV) in the flux formulation. The terms following these  
     30and advection part (COR+ADV) in the flux formulation. The terms following these  
    3231are the pressure gradient contributions (HPG, Hydrostatic Pressure Gradient,  
    3332and SPG, Surface Pressure Gradient); and contributions from lateral diffusion  
    … …  
    6059 
    6160$\ $\newline    % force a new ligne 
     61 
     62% ================================================================ 
     63% Sea Surface Height evolution & Diagnostics variables 
     64% ================================================================ 
     65\section{Sea surface height and diagnostic variables ($\eta$, $\zeta$, $\chi$, $w$)} 
     66\label{DYN_divcur_wzv} 
     67 
     68%-------------------------------------------------------------------------------------------------------------- 
     69%           Horizontal divergence and relative vorticity 
     70%-------------------------------------------------------------------------------------------------------------- 
     71\subsection   [Horizontal divergence and relative vorticity (\textit{divcur})] 
     72         {Horizontal divergence and relative vorticity (\mdl{divcur})} 
     73\label{DYN_divcur} 
     74 
     75The vorticity is defined at an $f$-point ($i.e.$ corner point) as follows: 
     76\begin{equation} \label{Eq_divcur_cur} 
     77\zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta _{i+1/2} \left[ {e_{2v}\;v} \right] 
     78                          -\delta _{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right) 
     79\end{equation}  
     80 
     81The horizontal divergence is defined at a $T$-point. It is given by: 
     82\begin{equation} \label{Eq_divcur_div} 
     83\chi =\frac{1}{e_{1t}\,e_{2t}\,e_{3t} } 
     84      \left( {\delta _i \left[ {e_{2u}\,e_{3u}\,u} \right] 
     85             +\delta _j \left[ {e_{1v}\,e_{3v}\,v} \right]} \right) 
     86\end{equation}  
     87 
     88Note that in the $z$-coordinate with full step (\key{zco} is defined),  
     89$e_{3u}$=$e_{3v}$=$e_{3f}$ so that they cancel in \eqref{Eq_divcur_div}. 
     90 
     91Note also that whereas the vorticity have the same discrete expression in $z$-  
     92and $s$-coordinate, its physical meaning is not identical. $\zeta$ is a pseudo  
     93vorticity along $s$-surfaces (only pseudo because $(u,v)$ are still defined along  
     94geopotential surfaces, but are no more necessary defined at the same depth). 
     95 
     96The vorticity and divergence at the \textit{before} step are used in the computation  
     97of the horizontal diffusion of momentum. Note that because they have been  
     98calculated prior to the Asselin filtering of the \textit{before} velocities, the  
     99\textit{before} vorticity and divergence arrays must be included in the restart file  
     100to ensure perfect restartability. The vorticity and divergence at the \textit{now}  
     101time step are used for the computation of the nonlinear advection and of the  
     102vertical velocity respectively.  
     103 
     104%-------------------------------------------------------------------------------------------------------------- 
     105%           Sea Surface Height evolution 
     106%-------------------------------------------------------------------------------------------------------------- 
     107\subsection   [Sea surface height evolution and vertical velocity (\textit{sshwzv})] 
     108         {Horizontal divergence and relative vorticity (\mdl{sshwzv})} 
     109\label{DYN_sshwzv} 
     110 
     111The sea surface height is given by : 
     112\begin{equation} \label{Eq_dynspg_ssh} 
     113\begin{aligned} 
     114\frac{\partial \eta }{\partial t} 
     115&\equiv    \frac{1}{e_{1t} e_{2t} }\sum\limits_k { \left(  \delta _i \left[ {e_{2u}\,e_{3u}\;u} \right] 
     116                                                                                  +\delta _j \left[ {e_{1v}\,e_{3v}\;v} \right]  \right) }  
     117           -    \frac{\textit{emp}}{\rho _w }   \\ 
     118&\equiv    \sum\limits_k {\chi \ e_{3t}}  -  \frac{\textit{emp}}{\rho _w } 
     119\end{aligned} 
     120\end{equation} 
     121where \textit{emp} is the surface freshwater budget (evaporation minus precipitation),  
     122expressed in Kg/m$^2$/s (which is equal to mm/s), and $\rho _w$=1,000~Kg/m$^3$  
     123is the volumic mass of pure water. If river runoff is expressed as a surface freshwater  
     124flux (see \S\ref{SBC}) then \textit{emp} can be written as the evaporation minus  
     125precipitation, minus the river runoff. The sea-surface height is evaluated  
     126using exactly the same time stepping as the tracer equation \eqref{Eq_tra_nxt}:  
     127a leapfrog scheme in combination with an Asselin time filter, $i.e.$ the velocity appearing  
     128in \eqref{Eq_dynspg_ssh} is centred in time (\textit{now} velocity).  
     129This is of paramount importance. Substituing $T$ by $1$ in the tracer equation and summing 
     130over the water column must lead to the sea surface height equation otherwise tracer content 
     131could not be conserved \ref{Griffies_al_MWR01, LeclairMadec2009}. 
     132 
     133The vertical velocity is computed by an upward integration of the horizontal  
     134divergence from the bottom, taken into account the change of the thickness of the levels : 
     135 
     136\begin{equation} \label{Eq_wzv} 
     137\left\{   \begin{aligned} 
     138&\left. w \right|_{3/2} \quad= 0    \\ 
     139&\left. w \right|_{k+1/2}     = \left. w \right|_{k+1/2}  + e_{3t}\;  \left. \chi \right|_k   
     140                                         - \frac{ e_{3t}^{t+1} - e_{3t}^{t-1} } {2 \Delta t} 
     141\end{aligned}   \right. 
     142\end{equation} 
     143 
     144In case of a non-linear free surface (\key{vvl}), the top vertical velocity is $-\textit{emp}/\rho_w$,  
     145as changes in the the divergence of the barotropic transport is absorbed in the change  
     146of the levels thickness.re-oriented downward co 
     147In case of linear free surface, the time derivative in \eqref{Eq_wzv} cancel out. 
     148The upper boundary condition applies at a fixed level $z=0$. The top vertical velocity  
     149is thus equal to the divergence of the barotropic transport ($i.e.$ the first term in the 
     150right-hand-side of \eqref{Eq_dynspg_ssh}). 
     151 
     152Note also that whereas the vertical velocity has the same discrete  
     153expression in $z$- and $s$-coordinate, its physical meaning is not the same:  
     154in the second case, $w$ is the velocity normal to the $s$-surfaces.  
     155Note also that the $k$-axis is re-oriented downward in the \textsc{fortran} code compare  
     156to the indexing used in the semi-discrete equations such as \eqref{Eq_wzv}  
     157(see  \S\ref{DOM_Num_Index_vertical}).  
     158 
     159 
    62160% ================================================================ 
    63161% Coriolis and Advection terms: vector invariant form 
    … …  
    66164\label{DYN_adv_cor_vect} 
    67165%-----------------------------------------nam_dynadv---------------------------------------------------- 
    68 \namdisplay{nam_dynadv}  
     166\namdisplay{namdyn_adv}  
    69167%------------------------------------------------------------------------------------------------------------- 
    70168 
    … …  
    85183\label{DYN_vor} 
    86184%------------------------------------------nam_dynvor---------------------------------------------------- 
    87 \namdisplay{nam_dynvor}  
     185\namdisplay{namdyn_vor}  
    88186%------------------------------------------------------------------------------------------------------------- 
    89187 
    90 Different discretisations of the vorticity term (\textit{ln\_dynvor\_xxx}=.true.) are  
    91 available: conserving potential enstrophy of horizontally non-divergent flow ;  
    92 conserving horizontal kinetic energy ; or conserving potential enstrophy for the  
    93 relative vorticity term and horizontal kinetic energy for the planetary vorticity  
    94 term (see  Appendix~\ref{Apdx_C}). The vorticity terms are given below for the  
    95 general case, but note that in the full step $z$-coordinate (\key{zco} is defined),  
    96 $e_{3u} =e_{3v} =e_{3f}$ so that the vertical scale factors disappear. They are  
    97 all computed in dedicated routines that can be found in the \mdl{dynvor} module. 
     188Four discretisations of the vorticity term (\textit{ln\_dynvor\_xxx}=.true.) are available:  
     189conserving potential enstrophy of horizontally non-divergent flow (ENS scheme) ;  
     190conserving horizontal kinetic energy (ENE scheme) ; conserving potential enstrophy for  
     191the relative vorticity term and horizontal kinetic energy for the planetary vorticity  
     192term (MIX scheme) ; or conserving both the potential enstrophy of horizontally non-divergent  
     193flow and horizontal kinetic energy (ENE scheme) (see  Appendix~\ref{Apdx_C_vor_zad}).  
     194The vorticity terms are given below for the general case, but note that in the full step  
     195$z$-coordinate (\key{zco} is defined), $e_{3u}$=$e_{3v}$=$e_{3f}$ so that the vertical scale  
     196factors disappear. They are all computed in dedicated routines that can be found in  
     197the \mdl{dynvor} module. 
    98198 
    99199%------------------------------------------------------------- 
    … …  
    106206vorticity term provides a global conservation of the enstrophy  
    107207($ [ (\zeta +f ) / e_{3f} ]^2 $ in $s$-coordinates) for a horizontally non-divergent  
    108 flow ($i.e.$ $\chi=0$), but does not conserve the total kinetic energy. It is given by: 
     208flow ($i.e.$ $\chi$=$0$), but does not conserve the total kinetic energy. It is given by: 
    109209\begin{equation} \label{Eq_dynvor_ens} 
    110210\left\{  
    111211\begin{aligned} 
    112212{+\frac{1}{e_{1u} } } & {\overline {\left( { \frac{\zeta +f}{e_{3f} }} \right)} }^{\,i}  
    113                                 & {\overline{\overline {\left( {e_{1v} e_{3v} v} \right)}} }^{\,i, j+1/2}    \\ 
     213                                & {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i, j+1/2}    \\ 
    114214{- \frac{1}{e_{2v} } } & {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)} }^{\,j}   
    115                                 & {\overline{\overline {\left( {e_{2u} e_{3u} u} \right)}} }^{\,i+1/2, j}   
     215                                & {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2, j}   
    116216\end{aligned}  
    117217 \right. 
    … …  
    129229\left\{   \begin{aligned} 
    130230{+\frac{1}{e_{1u}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right) 
    131                             \;  \overline {\left( {e_{1v} e_{3v} v} \right)} ^{\,i+1/2}} }^{\,j} }    \\ 
     231                            \;  \overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} }    \\ 
    132232{- \frac{1}{e_{2v}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right) 
    133                             \;  \overline {\left( {e_{2u} e_{3u} u} \right)} ^{\,j+1/2}} }^{\,i} } 
     233                            \;  \overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} } 
    134234\end{aligned}    \right. 
    135235\end{equation}  
    … …  
    146246to the planetary vorticity term. 
    147247\begin{equation} \label{Eq_dynvor_mix} 
    148 \left\{ { 
    149 \begin{aligned} 
     248\left\{ {     \begin{aligned} 
    150249 {+\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^{\,i}  
    151  \; {\overline{\overline {\left( {e_{1v} \; e_{3v} \ v} \right)}} }^{\,i,j+1/2} -\frac{1}{e_{1u} } 
     250 \; {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i,j+1/2} -\frac{1}{e_{1u} } 
    152251 \; {\overline {\left( {\frac{f}{e_{3f} }} \right)  
    153  \;\overline {\left( {e_{1v} \; e_{3v} \ v} \right)} ^{\,i+1/2}} }^{\,j} } \\ 
     252 \;\overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} } \\ 
    154253{-\frac{1}{e_{2v} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^j 
    155  \; {\overline{\overline {\left( {e_{2u} \; e_{3u} \ u} \right)}} }^{\,i+1/2,j} +\frac{1}{e_{2v} } 
     254 \; {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2,j} +\frac{1}{e_{2v} } 
    156255 \; {\overline {\left( {\frac{f}{e_{3f} }} \right) 
    157  \;\overline {\left( {e_{2u}\; e_{3u} \ u} \right)} ^{\,j+1/2}} }^{\,i} } \hfill 
    158 \end{aligned}  
    159 } \right. 
     256 \;\overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} } \hfill 
     257\end{aligned}     } \right. 
    160258\end{equation}  
    161259 
    … …  
    166264\label{DYN_vor_een} 
    167265 
    168 In the energy and enstrophy conserving scheme (EEN scheme), the vorticity term  
    169 is  evaluated using the vorticity advection scheme of \citet{Arakawa1990}.  
    170 This scheme conserves both total energy and potential enstrophy in the limit of  
    171 horizontally nondivergent flow ($i.e. \ \chi=0$). While EEN is more complicated  
    172 than ENS or ENE and does not conserve potential enstrophy and total energy in  
    173 general flow, it tolerates arbitrarily thin layers. This feature is essential for  
    174 $z$-coordinate with partial step.  
    175 %%% 
    176 \gmcomment{gm :   it actually conserve kinetic energy  !   show that in appendix C } 
    177 %%% 
    178  
    179 The \citet{Arakawa1990} vorticity advection scheme for a single layer is modified  
    180 for spherical coordinates as described by \citet{Arakawa1981} to obtain the EEN  
    181 scheme. The potential vorticity, defined at an $f$-point, is:  
     266In both the ENS and ENE schemes, it is apparent that the combination of $i$ and $j$  
     267averages of the velocity allows for the presence of grid point oscillation structures  
     268that will be invisible to the operator. These structures are \textit{computational modes}  
     269that will beat least partly damped by the momentum diffusive operator ($i.e.$ the  
     270subgrid-scale advection), but not by the resolved advection term. These two schemes 
     271therefore do not contribute to dump grid point noise in the horizontal velocity field,  
     272which results in more noise in vertical velocity field, an undesired feature. This is a well-known  
     273characteristics of $C$-grid discretization where $u$ and $v$ are located at different grid point, 
     274a price to pay to avoid a double averaging on the pressure gradient term as in $B$-grid.  
     275To circumvent this, Adcroft (ADD REF HERE)  
     276we have proposed to introduce a second velocity ... blahblah....  
     277 
     278Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves.... 
     279 
     280A very nice solution to that problem was proposed by \citet{Arakawa_Hsu_MWR90}. The idea is 
     281to get rid of the double averaging by considering triad combinations of vorticity.  
     282It is noteworthy that this solution is conceptually quite similar to the one proposed by 
     283\citep{Griffies_al_JPO98} for the discretization of the iso-neutral diffusive operator. 
     284 
     285The \citet{Arakawa_Hsu_MWR90} vorticity advection scheme for a single layer is modified  
     286for spherical coordinates as described by \citet{Arakawa_Lamb_MWR81} to obtain the EEN scheme.  
     287Let first provides the discrete expression of the potential vorticity, $q$, defined at an $f$-point:  
    182288\begin{equation} \label{Eq_pot_vor} 
    183 q_f  = \frac{\zeta +f} {e_{3f} } 
     289q  = \frac{\zeta +f} {e_{3f} } 
    184290\end{equation} 
    185291where the relative vorticity is defined by (\ref{Eq_divcur_cur}), the Coriolis parameter  
    … …  
    202308\textbf{j}- directions uses the masked vertical scale factor but is always divided by  
    203309$4$, not by the sum of the mask at $T$-point. This preserves the continuity of  
    204 $e_{3f}$ when one or more of the neighbouring $e_{3T}$ tends to zero and  
    205 extends by continuity the value of $e_{3f}$ in the land areas.  
    206 %%% 
    207 \gmcomment{this has to be further investigate in case of several step topography} 
    208 %%% 
     310$e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and  
     311extends by continuity the value of $e_{3f}$ in the land areas. This feature is essential for  
     312$z$-coordinate with partial step. 
     313 
     314 
     315Then, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at $T$-point as  
     316the following triad combinations of the neighbouring potential vorticities defined at f-point  
     317(Fig.~\ref{Fig_DYN_een_triad}):  
     318\begin{equation} \label{Q_triads} 
     319_i^j \mathbb{Q}^{i_p}_{j_p} 
     320= \frac{1}{12} \ \left(   q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p}  \right) 
     321\end{equation} 
     322where the indices $i_p$ and $k_p$ take the following value: $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$.  
    209323 
    210324The vorticity terms are represented as:  
    … …  
    212326\left\{ { 
    213327\begin{aligned} 
    214  +q\,e_3 \, v  &\equiv +\frac{1}{e_{1u} }  \left[  
    215 {{\begin{array}{*{20}c} 
    216       {\,\ \ a_{j+1/2}^{i   } \left( {e_{1v} e_{3v} \ v} \right)_{j+1}^{i+1/2} }  
    217    {\,+\,b_{j+1/2}^{i   } \left( {e_{1v} e_{3v} \ v} \right)_{j+1}^{i-1/2}  } \\ 
    218  \\ 
    219      {  +\,c_{j-1/2}^{i   }  \left( {e_{1v} e_{3v} \ v} \right)_{j    }^{i+1/2}         }  
    220    {\,+\,d_{j+1/2}^{i   } \left( {e_{1v} e_{3v} \ v} \right)_{j+1}^{i+1/2} } \\ 
    221 \end{array} }} \right] \\  
    222 \\ 
    223 -q\,e_3 \,u       &\equiv -\frac{1}{e_{2v} }  \left[  
    224 {{\begin{array}{*{20}c} 
    225    {\,\ \ a_{j-1/2}^{i   }  \left( {e_{2u} e_{3u} \ u} \right)_{j+1}^{i+1/2} }  
    226    {\,+\,b_{j-1/2}^{i+1}  \left( {e_{2u} e_{3u} \ u} \right)_{j+1/2}^{i+1} } \\ 
    227  \\ 
    228       {  +\,c_{j+1/2}^{i+1} \left( {e_{2u} e_{3u} \ u} \right)_{j+1/2}^{i+1} }  
    229    {\,+\,d_{j+1/2}^{i   }  \left( {e_{2u} e_{3u} \ u} \right)_{j+1/2}^{i   } } \\ 
    230 \end{array} }} \right] 
     328 +q\,e_3 \, v  &\equiv +\frac{1}{e_{1u} }   \sum_{\substack{i_p,\,k_p}}  
     329                         {^{i+1/2-i_p}_j}  \mathbb{Q}^{i_p}_{j_p}  \left( e_{1v}\,e_{3v} \;v  \right)^{i+1/2-i_p}_{j+j_p}   \\ 
     330 - q\,e_3 \, u     &\equiv -\frac{1}{e_{2v} }    \sum_{\substack{i_p,\,k_p}}  
     331                         {^i_{j+1/2-j_p}}  \mathbb{Q}^{i_p}_{j_p}  \left( e_{2u}\,e_{3u} \;u  \right)^{i+i_p}_{j+1/2-j_p}   \\ 
    231332\end{aligned}  
    232333} \right. 
    233334\end{equation}  
    234 where $a$, $b$, $c$ and $d$ are the following triad combinations of the  
    235 neighbouring potential vorticities (Fig.~\ref{Fig_DYN_een_triad}):  
    236 \begin{equation} \label{Eq_een_triads} 
    237 \left\{  
    238 \begin{aligned} 
    239  a_{\,j+1/2}^i & = \frac{1}  {12} \left( q_{j+1/2}^{i+1} + q_{j+1 /2}^i + q_{j-1/2}^i  \right)    \\  
    240  \\ 
    241  b_{\,j+1/2}^i & = \frac{1}  {12} \left( q_{j+1/2}^{i-1} +q_{j+1/2}^i +q_{j-1/2}^i   \right)     \\  
    242 \\ 
    243  c_{\,j+1/2}^i & = \frac{1}  {12} \left( q_{j-1/2}^{i-1} +q_{j+1/2}^i +q_{j-1/2}^i   \right)     \\  
    244 \\ 
    245  d_{\,j+1/2}^i & = \frac{1}  {12} \left( q_{j-1/2}^{i+1} +q_{j+1/2}^i +q_{j-1/2}^i  \right)     \\  
    246 \end{aligned}  
    247 \right. 
    248 \end{equation} 
     335 
     336This EEN scheme in fact combines the conservation properties of ENS and ENE schemes.  
     337It conserves both total energy and potential enstrophy in the limit of horizontally  
     338nondivergent flow ($i.e.$ $\chi$=$0$) (see  Appendix~\ref{Apdx_C_vor_zad}).  
     339Applied to a realistic ocean configuration, it has been shown that its larger mantis 
     340leads to a significant reduction of the noise in the vertical velocity field \citep{Le_Sommer_al_OM09}.  
     341Furthermore, used in combination with partial step representation of bottom topography, 
     342it improves the interaction between current and topography, leading to a larger 
     343topostrophy of the flow  \citep{Barnier_al_OD06, Penduff_al_OS07}.  
    249344 
    250345%-------------------------------------------------------------------------------------------------------------- 
    … …  
    260355\begin{equation} \label{Eq_dynkeg} 
    261356\left\{ \begin{aligned} 
    262  -\frac{1}{2 \; e_{1u} }  
    263  & \ \delta _{i+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right]   \\ 
    264  -\frac{1}{2 \; e_{2v} }  
    265  & \ \delta _{j+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right]     
     357 -\frac{1}{2 \; e_{1u} }  & \ \delta _{i+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right]   \\ 
     358 -\frac{1}{2 \; e_{2v} }  & \ \delta _{j+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right]     
    266359\end{aligned} \right. 
    267360\end{equation}  
    … …  
    280373\begin{equation} \label{Eq_dynzad} 
    281374\left\{     \begin{aligned} 
    282  -\frac{1}  { e_{1u}\,e_{2u}\,e_{3u} }  &   
    283   \ {\overline {\overline{ e_{1T}\,e_{2T}\,w } ^{\,i+1/2}  \;\delta _{k+1/2} \left[ u \right]  }^{\,k}   } \\ 
    284  -\frac{1}  { e_{1v}\,e_{2v}\,e_{3v} }  & 
    285   \ {\overline {\overline{ e_{1T}\,e_{2T}\,w } ^{\,j+1/2}  \;\delta _{k+1/2} \left[ u \right]  }^{\,k}   } 
    286 \end{aligned} \right. 
     375-\frac{1} {e_{1u}\,e_{2u}\,e_{3u}} &\ \overline{\ \overline{ e_{1t}\,e_{2t}\;w } ^{\,i+1/2}  \;\delta _{k+1/2} \left[ u \right]\  }^{\,k}  \\ 
     376-\frac{1} {e_{1v}\,e_{2v}\,e_{3v}}  &\ \overline{\ \overline{ e_{1t}\,e_{2t}\;w } ^{\,j+1/2}  \;\delta _{k+1/2} \left[ u \right]\  }^{\,k}  
     377\end{aligned}         \right. 
    287378\end{equation}  
    288379 
    … …  
    293384\label{DYN_adv_cor_flux} 
    294385%------------------------------------------nam_dynadv---------------------------------------------------- 
    295 \namdisplay{nam_dynadv}  
     386\namdisplay{namdyn_adv}  
    296387%------------------------------------------------------------------------------------------------------------- 
    297388 
    … …  
    315406\begin{multline} \label{Eq_dyncor_metric} 
    316407f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i}  -  u\frac{\partial e_1 }{\partial j}} \right)  \\ 
    317    \equiv   f + \frac{1}{e_{1f} e_{2f} }  
    318    \left( { \ \overline v ^{i+1/2}\delta _{i+1/2} \left[ {e_{2u} } \right]   
    319             -  \overline u ^{j+1/2}\delta _{j+1/2} \left[ {e_{1u} } \right]  }  \ \right) 
     408   \equiv   f + \frac{1}{e_{1f} e_{2f} } \left( { \ \overline v ^{i+1/2}\delta _{i+1/2} \left[ {e_{2u} } \right]   
     409                                                                 -  \overline u ^{j+1/2}\delta _{j+1/2} \left[ {e_{1u} } \right]  }  \ \right) 
    320410\end{multline}  
    321411 
    … …  
    338428\begin{aligned} 
    339429\frac{1}{e_{1u}\,e_{2u}\,e_{3u}}  
    340 \left(      \delta _{i+1/2} \left[ \overline{e_{2u}\,e_{3u}\ u }^{i       }  \ u_T      \right]     
    341           + \delta _{j       } \left[ \overline{e_{1u}\,e_{3u}\ v }^{i+1/2}  \ u_F      \right] \right.  \ \;   \\ 
    342 \left.   + \delta _{k      } \left[ \overline{e_{1w}\,e_{2w} w}^{i+1/2}  \ u_{uw} \right] \right)   \\ 
     430\left(      \delta _{i+1/2} \left[ \overline{e_{2u}\,e_{3u}\;u }^{i       }  \ u_t      \right]     
     431          + \delta _{j       } \left[ \overline{e_{1u}\,e_{3u}\;v }^{i+1/2}  \ u_f      \right] \right.  \ \;   \\ 
     432\left.   + \delta _{k      } \left[ \overline{e_{1w}\,e_{2w}\;w}^{i+1/2}  \ u_{uw} \right] \right)   \\ 
    343433\\ 
    344434\frac{1}{e_{1v}\,e_{2v}\,e_{3v}}  
    345 \left(   \delta _{i      } \left[  \overline{e_{2u}\,e_{3u } \ u }^{j+1/2} \ v_F       \right]  
    346          + \delta _{j+1/2} \left[ \overline{e_{1u}\,e_{3u } \ v }^{i       } \ v_T       \right] \right.  \ \, \\ 
    347 \left.  + \delta _{k     } \left[  \overline{e_{1w}\,e_{2w} \ w}^{j+1/2} \ v_{vw}  \right] \right) \\ 
     435\left(     \delta _{i       } \left[ \overline{e_{2u}\,e_{3u }\;u }^{j+1/2} \ v_f       \right]  
     436         + \delta _{j+1/2} \left[ \overline{e_{1u}\,e_{3u }\;v }^{i       } \ v_t       \right] \right.  \ \, \, \\ 
     437\left.  + \delta _{k      } \left[ \overline{e_{1w}\,e_{2w}\;w}^{j+1/2} \ v_{vw}  \right] \right) \\ 
    348438\end{aligned} 
    349439\right. 
    … …  
    369459\begin{equation} \label{Eq_dynadv_cen2} 
    370460\left\{     \begin{aligned} 
    371  u_T^{cen2} &=\overline u^{i }      \quad &   
    372   u_F^{cen2} &=\overline u^{j+1/2}     \quad & 
    373  u_{uw}^{cen2} &=\overline u^{k+1/2}      \\ 
    374  v_F^{cen2} &=\overline v ^{i+1/2}     \quad & 
    375  v_F^{cen2} &=\overline v^j      \quad & 
    376  v_{vw}^{cen2} &=\overline v ^{k+1/2}      \\ 
    377 \end{aligned} \right. 
     461 u_T^{cen2} &=\overline u^{i }       \quad &  u_F^{cen2} &=\overline u^{j+1/2}  \quad &  u_{uw}^{cen2} &=\overline u^{k+1/2}   \\ 
     462 v_F^{cen2} &=\overline v ^{i+1/2} \quad & v_F^{cen2} &=\overline v^j      \quad &  v_{vw}^{cen2} &=\overline v ^{k+1/2}  \\ 
     463\end{aligned}      \right. 
    378464\end{equation}  
    379465 
    … …  
    417503(centred in time), while the second term, which is the diffusive part of the scheme,  
    418504is evaluated using the \textit{before} velocity (forward in time). This is discussed  
    419 by \citet{Webb1998} in the context of the Quick advection scheme. 
     505by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme. 
    420506 
    421507Note that the UBS and Quadratic Upstream Interpolation for Convective Kinematics  
    422508(QUICK) schemes only differ by one coefficient. Substituting $1/6$ with $1/8$ in  
    423 (\ref{Eq_dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb1998}.  
     509(\ref{Eq_dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.  
    424510This option is not available through a namelist parameter, since the $1/6$ coefficient  
    425511is hard coded. Nevertheless it is quite easy to make the substitution in  
    … …  
    440526\label{DYN_hpg} 
    441527%------------------------------------------nam_dynhpg--------------------------------------------------- 
    442 \namdisplay{nam_dynhpg}  
    443 \namdisplay{namflg}  
     528\namdisplay{namdyn_hpg}  
    444529%------------------------------------------------------------------------------------------------------------- 
    445 %%% 
    446 \gmcomment{Suppress the namflg namelist and incorporate it in the namhpg namelist} 
    447 %%% 
    448530 
    449531The key distinction between the different algorithms used for the hydrostatic  
    … …  
    454536 
    455537The hydrostatic pressure gradient term is evaluated either using a leapfrog scheme,  
    456 $i.e.$ the density appearing in its expression is centred in time (\emph{now} rho), or  
     538$i.e.$ the density appearing in its expression is centred in time (\emph{now} $\rho$), or  
    457539a semi-implcit scheme. At the lateral boundaries either free slip, no slip or partial slip  
    458540boundary conditions are applied. 
    … …  
    495577Note that the $1/2$ factor in (\ref{Eq_dynhpg_zco_surf}) is adequate because of  
    496578the definition of $e_{3w}$ as the vertical derivative of the scale factor at the surface  
    497 level ($z=0)$.  
     579level ($z=0$). Note also that in case of variable volume level (\key{vvl} defined), the  
     580surface pressure gradient is included in \eqref{Eq_dynhpg_zco_surf} and \eqref{Eq_dynhpg_zco}  
     581through the space and time variations of the vertical scale factor $e_{3w}$. 
    498582 
    499583%-------------------------------------------------------------------------------------------------------------- 
    … …  
    528612gradient options are coded, but they are not yet fully documented or tested.  
    529613 
    530 $\bullet$ Traditional coding (see for example \citet{Madec1996}: (\np{ln\_dynhpg\_sco}=.true.,  
     614$\bullet$ Traditional coding (see for example \citet{Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}=.true.,  
    531615\np{ln\_dynhpg\_hel}=.true.) 
    532616\begin{equation} \label{Eq_dynhpg_sco} 
    533617\left\{ \begin{aligned} 
    534618 - \frac{1}                   {\rho_o \, e_{1u}} \;   \delta _{i+1/2} \left[  p^h  \right]  
    535 + \frac{g\; \overline {\rho}^{i+1/2}}  {\rho_o \, e_{1u}} \;   \delta _{i+1/2} \left[  z_T  \right]    \\ 
     619+ \frac{g\; \overline {\rho}^{i+1/2}}  {\rho_o \, e_{1u}} \;   \delta _{i+1/2} \left[  z_t   \right]    \\ 
    536620 - \frac{1}                   {\rho_o \, e_{2v}} \;   \delta _{j+1/2} \left[  p^h  \right]   
    537 + \frac{g\; \overline {\rho}^{j+1/2}}  {\rho_o \, e_{2v}} \;   \delta _{j+1/2} \left[  z_T  \right]    \\ 
     621+ \frac{g\; \overline {\rho}^{j+1/2}}  {\rho_o \, e_{2v}} \;   \delta _{j+1/2} \left[  z_t   \right]    \\ 
    538622\end{aligned} \right. 
    539623\end{equation}  
    … …  
    551635(\np{ln\_dynhpg\_djc}=.true.) 
    552636 
    553 $\bullet$ Rotated axes scheme (rot) \citep{Thiem2006} (\np{ln\_dynhpg\_rot}=.true.) 
     637$\bullet$ Rotated axes scheme (rot) \citep{Thiem_Berntsen_OM06} (\np{ln\_dynhpg\_rot}=.true.) 
    554638 
    555639Note that expression \eqref{Eq_dynhpg_sco} is used when the variable volume  
    … …  
    617701Note that in the semi-implicit case, it is necessary to save the filtered density, an  
    618702extra three-dimensional field, in the restart file to restart the model with exact  
    619 reproducibility. This option is controlled by the namelist parameter  
    620 \np{nn\_dynhpg\_rst}=.true.. 
     703reproducibility. This option is controlled by  \np{nn\_dynhpg\_rst}, a namelist parameter. 
    621704 
    622705% ================================================================ 
    … …  
    627710\label{DYN_hpg_spg} 
    628711%-----------------------------------------nam_dynspg---------------------------------------------------- 
    629 \namdisplay{nam_dynspg}  
     712\namdisplay{namdyn_spg}  
    630713%------------------------------------------------------------------------------------------------------------ 
    631714 
    632 The form of the surface pressure gradient term is dependent on the representation  
    633 of the free surface (\S\ref{PE_hor_pg}). The main distinction is between the fixed  
    634 volume case (linear free surface or rigid lid) and the variable volume case  
    635 (nonlinear free surface, \key{vvl} is defined). In the linear free surface case  
    636 (\S\ref{PE_free_surface}) and the rigid lid case (\S\ref{PE_rigid_lid}), the vertical  
    637 scale factors $e_{3}$ are fixed in time, whilst in the nonlinear case  
    638 (\S\ref{PE_free_surface}) they are time-dependent. With both linear and nonlinear  
    639 free surface, external gravity waves are allowed in the equations, which imposes  
    640 a very small time step when an explicit time stepping is used. Two methods are  
    641 proposed to allow a longer time step for the three-dimensional equations: the  
    642 filtered free surface method, which involves a modification of the continuous  
    643 equations (see \eqref{Eq_PE_flt}), and the split-explicit free surface method  
    644 described below. The extra term introduced in the filtered method is calculated  
    645 implicitly, so that the update of the $next$ velocities is done in module  
    646 \mdl{dynspg\_flt} and not in \mdl{dynnxt}. 
    647  
    648 %-------------------------------------------------------------------------------------------------------------- 
    649 % Linear free surface formulation 
    650 %-------------------------------------------------------------------------------------------------------------- 
    651 \subsection{Linear free surface formulation (\key{exp} or \textbf{\_ts} or \textbf{\_flt})} 
    652 \label{DYN_spg_linear} 
    653  
    654 In the linear free surface formulation, the sea surface height is assumed to be  
    655 small compared to the thickness of the ocean levels, so that $(a)$ the time  
    656 evolution of the sea surface height becomes a linear equation, and $(b)$ the  
    657 thickness of the ocean levels is assumed to be constant with time.  
    658 As mentioned in (\S\ref{PE_free_surface}) the linearization affects the  
    659 conservation of tracers.  
    660  
    661 %------------------------------------------------------------- 
    662 % Explicit 
    663 %------------------------------------------------------------- 
    664 \subsubsection{Explicit (\key{dynspg\_exp})} 
     715$\ $\newline      %force an empty line 
     716 
     717The form of the surface pressure gradient term depends on how the user wants to handle  
     718the fast external gravity waves that are solution of the analytical equation (\S\ref{PE_hor_pg}).  
     719Three formulation are available, all controlled by a CPP key (ln\_dynspg\_xxx): 
     720a explicit formulation which required a small time step ; 
     721a filtered free surface formulation which allows a larger time step by adding a filtering  
     722term in the momentum equation ;  
     723and a plit-explicit free surface formulation, described below, which also allows a larger time step. 
     724 
     725The extra term introduced in the filtered method is calculated  
     726implicitly, so that a solver is used to compute it and that the update of the $next$  
     727velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. 
     728 
     729 
     730 
     731%-------------------------------------------------------------------------------------------------------------- 
     732% Explicit free surface formulation 
     733%-------------------------------------------------------------------------------------------------------------- 
     734\subsection{Explicit free surface (\key{dynspg\_exp})} 
    665735\label{DYN_spg_exp} 
    666736 
    667 In the explicit free surface formulation, the model time step is chosen to be  
    668 small enough to describe the external gravity waves (typically a few tens of  
    669 seconds). The sea surface height is given by : 
    670 \begin{equation} \label{Eq_dynspg_ssh} 
    671 \frac{\partial \eta }{\partial t}\equiv \frac{\text{EMP}}{\rho _w }+\frac{1}{e_{1T}  
    672 e_{2T} }\sum\limits_k {\left( {\delta _i \left[ {e_{2u} e_{3u} u}  
    673 \right]+\delta _j \left[ {e_{1v} e_{3v} v} \right]} \right)}  
    674 \end{equation} 
    675 where EMP is the surface freshwater budget, expressed in Kg/m$^2$/s  
    676 (which is equal to mm/s), and $\rho _w$=1,000~Kg/m$^3$ is the volumic  
    677 mass of pure water. If river runoff is expressed as a surface freshwater flux  
    678 (see \S\ref{SBC}) then EMP can be written as the evaporation minus  
    679 precipitation, minus the river runoff. The sea-surface height is evaluated  
    680 using a leapfrog scheme in combination with an Asselin time filter, $i.e.$  
    681 the velocity appearing in \eqref{Eq_dynspg_ssh} is centred in time  
    682 (\textit{now} velocity).  
     737In the explicit free surface formulation (\key{dynspg\_exp} defined), the model time step  
     738is chosen to be small enough to describe the external gravity waves (typically a few tens of seconds).  
    683739 
    684740The surface pressure gradient, also evaluated using a leap-frog scheme, is  
    … …  
    686742\begin{equation} \label{Eq_dynspg_exp} 
    687743\left\{ \begin{aligned} 
    688  - \frac{1}{e_{1u}} \;  \delta _{i+1/2} \left[ \,\eta\,  \right]    \\ 
    689  - \frac{1}{e_{2v}} \;  \delta _{j+1/2} \left[ \,\eta\,  \right]   
     744 - \frac{1}{e_{1u}\,\rho_o} \;   \delta _{i+1/2} \left[  \,\rho \,\eta\,  \right]   \\ 
     745 - \frac{1}{e_{2v}\,\rho_o} \;   \delta _{j+1/2} \left[  \,\rho \,\eta\,  \right]   
    690746\end{aligned} \right. 
    691747\end{equation}  
    692748 
    693 Consistent with the linearization, a factor of $\left. \rho \right|_{k=1} / \rho _o$  
    694 is omitted in \eqref{Eq_dynspg_exp}.  
    695  
    696 %------------------------------------------------------------- 
    697 % Split-explicit time-stepping 
    698 %------------------------------------------------------------- 
    699 \subsubsection{Split-explicit time-stepping (\key{dynspg\_ts})} 
     749Note that in the non-linear free surface case ($i.e.$ \key{vvl} defined), the surface pressure  
     750gradient is already included in the momentum tendency  through the level thickness variation  
     751when computing the hydrostatic pressure gradient. Thus, nothing is done in the \mdl{dynspg\_exp} module. 
     752 
     753%-------------------------------------------------------------------------------------------------------------- 
     754% Split-explict free surface formulation 
     755%-------------------------------------------------------------------------------------------------------------- 
     756\subsection{Split-Explicit free surface (\key{dynspg\_ts})} 
    700757\label{DYN_spg_ts} 
    701 %--------------------------------------------namdom---------------------------------------------------- 
    702 \namdisplay{namdom}  
    703 %-------------------------------------------------------------------------------------------------------------- 
     758 
     759In the split-explicit free surface formulation, also called time-splitting formulation  
     760(\key{dynspg\_ts} defined) 
     761 
    704762 
    705763The split-explicit free surface formulation used in \NEMO follows the one  
    706 proposed by \citet{Griffies2004}. The general idea is to solve the free surface  
    707 equation with a small time step \np{rdtbt}, while the three dimensional  
    708 prognostic variables are solved with a longer time step that is a multiple of  
    709 \np{rdtbt} (Fig.\ref {Fig_DYN_dynspg_ts}).  
     764proposed by \citet{Griffies_Bk04}. The general idea is to solve the free surface  
     765equation and the associated barotropic velocity equations with a smaller time  
     766step than $\Delta t$, the time step used for the three dimensional prognostic  
     767variables (Fig.\ref {Fig_DYN_dynspg_ts}).  
     768The size of the small time step, $\Delta_e$ (the external mode or barotropic time step) 
     769 is provided through \np{nn\_baro} namelist parameter as:  
     770$\Delta_e = \Delta / nn\_baro$. 
     771  
    710772 
    711773%>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   > 
    … …  
    739801shown by \citet{Levier2007} in the case of an analytical barotropic Kelvin wave.  
    740802 
    741 %------------------------------------------------------------- 
    742 % Filtered formulation  
    743 %------------------------------------------------------------- 
    744 \subsubsection{Filtered formulation (\key{dynspg\_flt})} 
    745 \label{DYN_spg_flt} 
    746  
    747 The filtered formulation follows the \citet{Roullet2000} implementation. The extra  
     803%-------------------------------------------------------------------------------------------------------------- 
     804% Filtered free surface formulation 
     805%-------------------------------------------------------------------------------------------------------------- 
     806\subsection{Filtered free surface (\key{dynspg\_flt})} 
     807\label{DYN_spg_fltp} 
     808 
     809 
     810The filtered formulation follows the \citet{Roullet_Madec_JGR00} implementation. The extra  
    748811term introduced in the equations (see {\S}I.2.2) is solved implicitly. The elliptic  
    749 solvers available in the code are documented in \S\ref{MISC}. The amplitude of  
    750 the extra term is given by the namelist variable \np{rnu}. The default value is 1,  
    751 as recommended by \citet{Roullet2000} 
    752  
    753 \gmcomment{\np{rnu}=1 to be suppressed from namelist !} 
    754  
    755 %------------------------------------------------------------- 
    756 % Non-linear free surface formulation  
    757 %------------------------------------------------------------- 
    758 \subsection{Non-linear free surface formulation (\key{vvl})} 
    759 \label{DYN_spg_vvl} 
    760  
    761 In the non-linear free surface formulation, the variations of volume are fully  
    762 taken into account. This option is presented in a report \citep{Levier2007}  
    763 available on the \NEMO web site. The three time-stepping methods (explicit,  
    764 split-explicit and filtered) are the same as in \S\ref{DYN_spg_linear} except  
    765 that the ocean depth is now time-dependent. In particular, this means that  
    766 in the filtered case, the matrix to be inverted has to be recomputed at each  
    767 time-step. 
    768  
    769 %-------------------------------------------------------------------------------------------------------------- 
    770 %           Rigid-lid formulation 
    771 %-------------------------------------------------------------------------------------------------------------- 
    772 \subsection{Rigid-lid formulation (\key{dynspg\_rl})} 
    773 \label{DYN_spg_rl} 
    774  
    775 With the rigid lid formulation, an elliptic equation has to be solved for  
    776 the barotropic streamfunction. For consistency this equation is obtained by  
    777 taking the discrete curl of the discrete vertical sum of the discrete  
    778 momentum equation:  
    779 \begin{equation}\label{Eq_dynspg_rl} 
    780 \frac{1}{\rho _o }\nabla _h p_s \equiv \left( {{\begin{array}{*{20}c} 
    781  {\overline M_u +\frac{1}{H\;e_2 } \delta_ j \left[ \partial_t \psi \right]}      \\ 
    782  \\ 
    783  {\overline M_v -\frac{1}{H\;e_1 }  \delta_i  \left[ \partial_t \psi \right]}        \\ 
    784 \end{array} }} \right) 
    785 \end{equation} 
    786  
    787 Here ${\rm {\bf M}}= \left( M_u,M_v \right)$ represents the collected  
    788 contributions of nonlinear, viscous and hydrostatic pressure gradient terms in  
    789 \eqref{Eq_PE_dyn} and the overbar indicates a vertical average over the  
    790 whole water column (i.e. from $z=-H$, the ocean bottom, to $z=0$, the rigid-lid).  
    791 The time derivative of $\psi$ is the solution of an elliptic equation: 
    792 \begin{multline} \label{Eq_bsf} 
    793    \delta_{i+1/2} \left[ \frac{e_{2v}}{H_v\;e_{1v}} \delta_{i} \left[  \partial_t \psi \right] \right] 
    794 + \delta_{j+1/2} \left[ \frac{e_{1u}}{H_u\;e_{2u}} \delta_{j} \left[  \partial_t \psi \right] \right] 
    795 \\ = 
    796   \delta_{i+1/2} \left[ e_{2v} M_v  \right] 
    797 - \delta_{j+1/2} \left[ e_{1u} M_u  \right] 
    798 \end{multline} 
    799  
    800 The elliptic solvers available in the code are documented in \S\ref{MISC}).  
    801 The boundary conditions must be given on all separate landmasses (islands),  
    802 which is done by integrating the vorticity equation around each island. This  
    803 requires identifying each island in the bathymetry file, a cumbersome task.  
    804 This explains why the rigid lid option is not recommended with complex  
    805 domains such as the global ocean. Parameters jpisl (number of islands) and  
    806 jpnisl (maximum number of points per island) of the \hf{par\_oce} file are  
    807 related to this option. 
     812solvers available in the code are documented in \S\ref{MISC}. 
    808813 
    809814 
    … …  
    815820\label{DYN_ldf} 
    816821%------------------------------------------nam_dynldf---------------------------------------------------- 
    817 \namdisplay{nam_dynldf}  
     822\namdisplay{namdyn_ldf}  
    818823%------------------------------------------------------------------------------------------------------------- 
    819824 
    … …  
    821826(rotated or not) or biharmonic operators. The coefficients may be constant  
    822827or spatially variable; the description of the coefficients is found in the chapter  
    823 on lateralphysics (Chap.\ref{LDF}). The lateral diffusion of momentum is  
    824 evaluated using a forward scheme, i.e. the velocity appearing in its expression  
     828on lateral physics (Chap.\ref{LDF}). The lateral diffusion of momentum is  
     829evaluated using a forward scheme, $i.e.$ the velocity appearing in its expression  
    825830is the \textit{before} velocity in time, except for the pure vertical component  
    826831that appears when a tensor of rotation is used. This latter term is solved  
    … …  
    850855As explained in \S\ref{PE_ldf}, this formulation (as the gradient of a divergence  
    851856and curl of the vorticity) preserves symmetry and ensures a complete  
    852 separation between the vorticity and divergence parts. Note that in the full step  
    853 $z$-coordinate (\key{zco} is defined), $e_{3u} =e_{3v} =e_{3f}$ so that they  
    854 cancel in the rotational part of \eqref{Eq_dynldf_lap}. 
     857separation between the vorticity and divergence parts of the momentum diffusion.  
     858Note that in the full step $z$-coordinate (\key{zco} is defined), $e_{3u} =e_{3v} =e_{3f}$  
     859so that they cancel in the rotational part of \eqref{Eq_dynldf_lap}. 
    855860 
    856861%-------------------------------------------------------------------------------------------------------------- 
    … …  
    875880 D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\ 
    876881&  \left\{\quad  {\delta _{i+1/2} \left[ {A_T^{lm}  \left(  
    877     {\frac{e_{2T} \; e_{3T} }{e_{1T} } \,\delta _{i}[u] 
    878    -e_{2T} \; r_{1T} \,\overline{\overline {\delta _{k+1/2}[u]}}^{\,i,\,k}} 
     882    {\frac{e_{2t} \; e_{3t} }{e_{1t} } \,\delta _{i}[u] 
     883   -e_{2t} \; r_{1t} \,\overline{\overline {\delta _{k+1/2}[u]}}^{\,i,\,k}} 
    879884 \right)} \right]}   \right. 
    880885\\  
    … …  
    902907 \right)} \right]}   \right. 
    903908\\  
    904 & \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1T}\,e_{3T} }{e_{2T}  
    905 }\,\delta _{j} [v] - e_{1T}\, r_{2T}  
     909& \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1t}\,e_{3t} }{e_{2t}  
     910}\,\delta _{j} [v] - e_{1t}\, r_{2t}  
    906911\,\overline{\overline {\delta _{k+1/2} [v]}} ^{\,j,\,k}}  
    907912\right)} \right]  
    … …  
    954959can be used for the vertical diffusion term : $(a)$ a forward time differencing  
    955960scheme (\np{ln\_zdfexp}=.true.) using a time splitting technique  
    956 (\np{n\_zdfexp} $>$ 1) or $(b)$ a backward (or implicit) time differencing scheme  
     961(\np{nn\_zdfexp} $>$ 1) or $(b)$ a backward (or implicit) time differencing scheme  
    957962(\np{ln\_zdfexp}=.false.) (see \S\ref{DOM_nxt}). Note that namelist variables  
    958 \np{ln\_zdfexp} and \np{n\_zdfexp} apply to both tracers and dynamics.  
     963\np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics.  
    959964 
    960965The formulation of the vertical subgrid scale physics is the same whatever  
    … …  
    10211026The general framework for dynamics time stepping is a leap-frog scheme,  
    10221027$i.e.$ a three level centred time scheme associated with an Asselin time filter  
    1023 (cf. \S\ref{DOM_nxt}) 
    1024 \begin{equation} \label{Eq_dynnxt} 
    1025 \begin{split} 
    1026 &u^{t+\Delta t} = u^{t-\Delta t} + 2 \, \Delta t  \ \text{RHS}_u^t   \\ 
    1027 \\ 
     1028(cf. \S\ref{DOM_nxt}). The scheme is applied to the velocity except when using  
     1029the flux form of momentum advection (cf. \S\ref{DYN_adv_cor_flux}) in variable  
     1030volume level case (\key{vvl} defined), where it has to be applied to the thickness  
     1031weighted velocity (see \S\ref{Apdx_A_momentum})   
     1032 
     1033$\bullet$ vector invariant form or linear free surface (\np{ln\_dynhpg\_vec}=.true. ; not \key{vvl}): 
     1034\begin{equation} \label{Eq_dynnxt_vec} 
     1035\left\{   \begin{aligned} 
     1036&u^{t+\Delta t} = u_f^{t-\Delta t} + 2\Delta t  \ \text{RHS}_u^t     \\ 
    10281037&u_f^t \;\quad = u^t+\gamma \,\left[ {u_f^{t-\Delta t} -2u^t+u^{t+\Delta t}} \right] 
    1029 \end{split} 
     1038\end{aligned}   \right. 
     1039\end{equation}  
     1040 
     1041$\bullet$ flux form and nonlinear free surface (\np{ln\_dynhpg\_vec}=.false. ; \key{vvl} defined): 
     1042\begin{equation} \label{Eq_dynnxt_flux} 
     1043\left\{   \begin{aligned} 
     1044&\left(e_{3u}\,u\right)^{t+\Delta t} = \left(e_{3u}\,u\right)_f^{t-\Delta t} + 2\Delta t \; e_{3u} \;\text{RHS}_u^t     \\ 
     1045&\left(e_{3u}\,u\right)_f^t \;\quad = \left(e_{3u}\,u\right)^t 
     1046  +\gamma \,\left[ {\left(e_{3u}\,u\right)_f^{t-\Delta t} -2\left(e_{3u}\,u\right)^t+\left(e_{3u}\,u\right)^{t+\Delta t}} \right] 
     1047\end{aligned}   \right. 
    10301048\end{equation}  
    10311049where RHS is the right hand side of the momentum equation, the subscript $f$  
    10321050denotes filtered values and $\gamma$ is the Asselin coefficient. $\gamma$ is  
    1033 initialized as \np{atfp} (namelist parameter). Its default value is \np{atfp} = 0.1. 
    1034  
    1035 Note that whith the filtered free surface, the update of the \textit{next} velocities  
     1051initialized as \np{nn\_atfp} (namelist parameter). Its default value is \np{nn\_atfp} = 0.1. 
     1052 
     1053Note that with the filtered free surface, the update of the \textit{next} velocities  
    10361054is done in the \mdl{dynsp\_flt} module, and only the swap of arrays  
    1037 and Asselin filtering is done in \mdl{dynnxt.} 
    1038  
    1039 % ================================================================ 
    1040 % Diagnostic variables 
    1041 % ================================================================ 
    1042 \section{Diagnostic variables ($\zeta$, $\chi$, $w$)} 
    1043 \label{DYN_divcur_wzv} 
    1044  
    1045 %-------------------------------------------------------------------------------------------------------------- 
    1046 %           Horizontal divergence and relative vorticity 
    1047 %-------------------------------------------------------------------------------------------------------------- 
    1048 \subsection   [Horizontal divergence and relative vorticity (\textit{divcur})] 
    1049          {Horizontal divergence and relative vorticity (\mdl{divcur})} 
    1050 \label{DYN_divcur} 
    1051  
    1052 The vorticity is defined at an $f$-point ($i.e.$ corner point) as follows: 
    1053 \begin{equation} \label{Eq_divcur_cur} 
    1054 \zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta _{i+1/2} \left[ {e_{2v}\;v} \right] 
    1055                           -\delta _{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right) 
    1056 \end{equation}  
    1057  
    1058 The horizontal divergence is defined at a $T$-point. It is given by: 
    1059 \begin{equation} \label{Eq_divcur_div} 
    1060 \chi =\frac{1}{e_{1T}\,e_{2T}\,e_{3T} } 
    1061       \left( {\delta _i \left[ {e_{2u}\,e_{3u}\,u} \right] 
    1062            +\delta _j \left[ {e_{1v}\,e_{3v}\,v} \right]} \right) 
    1063 \end{equation}  
    1064  
    1065 Note that in the $z$-coordinate with full step (\key{zco} is defined),  
    1066 $e_{3u} =e_{3v} =e_{3f}$ so that they cancel in \eqref{Eq_divcur_div}. 
    1067  
    1068 Note also that whereas the vorticity have the same discrete expression in $z$-  
    1069 and $s$-coordinate, its physical meaning is not identical. $\zeta$ is a pseudo  
    1070 vorticity along $s$-surfaces (only pseudo because $(u,v)$ are still defined along  
    1071 geopotential surfaces, but are no more necessary defined at the same depth). 
    1072  
    1073 The vorticity and divergence at the \textit{before} step are used in the computation  
    1074 of the horizontal diffusion of momentum. Note that because they have been  
    1075 calculated prior to the Asselin filtering of the \textit{before} velocities, the  
    1076 \textit{before} vorticity and divergence arrays must be included in the restart file  
    1077 to ensure perfect restartability. The vorticity and divergence at the \textit{now}  
    1078 time step are used for the computation of the nonlinear advection and of the  
    1079 vertical velocity respectively.  
    1080  
    1081 %-------------------------------------------------------------------------------------------------------------- 
    1082 %           Vertical Velocity 
    1083 %-------------------------------------------------------------------------------------------------------------- 
    1084 \subsection   [Vertical velocity (\textit{wzvmod})] 
    1085          {Vertical velocity (\mdl{wzvmod})} 
    1086 \label{DYN_wzv} 
    1087  
    1088 The vertical velocity is computed by an upward integration of the horizontal  
    1089 divergence from the bottom : 
    1090  
    1091 \begin{equation} \label{Eq_wzv} 
    1092 \left\{   \begin{aligned} 
    1093 &\left. w \right|_{3/2} \quad= 0    \\ 
    1094 \\ 
    1095 &\left. w \right|_{k+1/2}     = \left. w \right|_{k+1/2}  + e_{3t}\;  \left. \chi \right|_k   
    1096 \end{aligned}   \right. 
    1097 \end{equation}  
    1098  
    1099 With a free surface, the top vertical velocity is non-zero, due to the  
    1100 freshwater forcing and the variations of the free surface elevation. With a  
    1101 linear free surface or with a rigid lid, the upper boundary condition  
    1102 applies at a fixed level $z=0$. Note that in the rigid-lid case (\key{dynspg\_rl}  
    1103 is defined), the surface boundary condition ($\left. w \right|_\text{surface}=0)$ is  
    1104 automatically achieved at least at computer accuracy, due to the the way the  
    1105 surface pressure gradient is expressed in discrete form (Appendix~\ref{Apdx_C}). 
    1106  
    1107 Note also that whereas the vertical velocity has the same discrete  
    1108 expression in $z$- and $s$-coordinate, its physical meaning is not the same:  
    1109 in the second case, $w$ is the velocity normal to the $s$-surfaces.  
    1110  
    1111 With the variable volume option, the calculation of the vertical velocity is  
    1112 modified (see \citet{Levier2007}, report available on the \NEMO web site).  
    1113  
    1114 % ================================================================ 
     1055and Asselin filtering is done in \mdl{dynnxt}. 
     1056 
     1057 
     1058 
     1059% ================================================================ 
Note: See TracChangeset for help on using the changeset viewer.