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 6625 for branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters – NEMO

Ignore:
Timestamp:
2016-05-26T11:08:07+02:00 (8 years ago)
Author:
kingr
Message:

Rolled back to r6613

Location:
branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Abstracts_Foreword.tex

    r6617 r6625  
    1313be a flexible tool for studying the ocean and its interactions with the others components of  
    1414the earth climate system over a wide range of space and time scales.  
    15 Prognostic variables are the three-dimensional velocity field, a non-linear sea surface height,  
    16 the \textit{Conservative} Temperature and the \textit{Absolute} Salinity.  
    17 In the horizontal direction, the model uses a curvilinear orthogonal grid and in the vertical direction,  
    18 a full or partial step $z$-coordinate, or $s$-coordinate, or a mixture of the two.  
    19 The distribution of variables is a three-dimensional Arakawa C-type grid.  
    20 Various physical choices are available to describe ocean physics, including TKE, and GLS vertical physics.  
    21 Within NEMO, the ocean is interfaced with a sea-ice model (LIM or CICE), passive tracer and  
    22 biogeochemical models (TOP) and, via the OASIS coupler, with several atmospheric general circulation models.  
    23 It also support two-way grid embedding via the AGRIF software. 
     15Prognostic variables are the three-dimensional velocity field, a linear  
     16or non-linear sea surface height, the temperature and the salinity. In the horizontal direction,  
     17the model uses a curvilinear orthogonal grid and in the vertical direction, a full or partial step  
     18$z$-coordinate, or $s$-coordinate, or a mixture of the two. The distribution of variables is a  
     19three-dimensional Arakawa C-type grid. Various physical choices are available to describe  
     20ocean physics, including TKE, GLS and KPP vertical physics. Within NEMO, the ocean is  
     21interfaced with a sea-ice model (LIM v2 and v3), passive tracer and biogeochemical models (TOP)  
     22and, via the OASIS coupler, with several atmospheric general circulation models. It also  
     23support two-way grid embedding via the AGRIF software. 
    2424 
    2525% ================================================================ 
    … …  
    3131interactions avec les autres composantes du syst\`{e}me climatique terrestre.  
    3232Les variables pronostiques sont le champ tridimensionnel de vitesse, une hauteur de la mer  
    33 lin\'{e}aire, la Temp\'{e}rature Conservative et la Salinit\'{e} Absolue.  
     33lin\'{e}aire ou non, la temperature et la salinit\'{e}.  
    3434La distribution des variables se fait sur une grille C d'Arakawa tridimensionnelle utilisant une  
    3535coordonn\'{e}e verticale $z$ \`{a} niveaux entiers ou partiels, ou une coordonn\'{e}e s, ou encore  
    3636une combinaison des deux. Diff\'{e}rents choix sont propos\'{e}s pour d\'{e}crire la physique  
    37 oc\'{e}anique, incluant notamment des physiques verticales TKE et GLS. A travers l'infrastructure  
    38 NEMO, l'oc\'{e}an est interfac\'{e} avec des mod\`{e}les de glace de mer (LIM ou CICE),  
    39 de biog\'{e}ochimie marine et de traceurs passifs, et, via le coupleur OASIS, \`{a} plusieurs  
    40 mod\`{e}les de circulation g\'{e}n\'{e}rale atmosph\'{e}rique.  
    41 Il supporte \'{e}galement l'embo\^{i}tement interactif de maillages via le logiciel AGRIF. 
     37oc\'{e}anique, incluant notamment des physiques verticales TKE, GLS et KPP. A travers l'infrastructure  
     38NEMO, l'oc\'{e}an est interfac\'{e} avec des mod\`{e}les de glace de mer, de biog\'{e}ochimie  
     39et de traceurs passifs, et, via le coupleur OASIS, \`{a} plusieurs mod\`{e}les de circulation  
     40g\'{e}n\'{e}rale atmosph\'{e}rique. Il supporte \'{e}galement l'embo\^{i}tement interactif de  
     41maillages via le logiciel AGRIF. 
    4242}  
    4343 

  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Annex_C.tex

    r6617 r6625  
    410410\end{aligned}   } \right. 
    411411\end{equation}  
    412 where the indices $i_p$ and $j_p$ take the following value:  
     412where the indices $i_p$ and $k_p$ take the following value:  
    413413$i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$, 
    414414and the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$, defined at $T$-point, are given by:  
    … …  
    11031103The discrete formulation of the horizontal diffusion of momentum ensures the  
    11041104conservation of potential vorticity and the horizontal divergence, and the  
    1105 dissipation of the square of these quantities ($i.e.$ enstrophy and the  
     1105dissipation of the square of these quantities (i.e. enstrophy and the  
    11061106variance of the horizontal divergence) as well as the dissipation of the  
    11071107horizontal kinetic energy. In particular, when the eddy coefficients are  
    … …  
    11271127&\int \limits_D \frac{1} {e_3 } \textbf{k} \cdot \nabla \times  
    11281128   \Bigl[    \nabla_h  \left( A^{\,lm}\;\chi  \right) 
    1129            - \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right)    \Bigr]\;dv   \\  
    1130 %\end{flalign*} 
     1129             - \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right)    \Bigr]\;dv  = 0 
     1130\end{flalign*} 
    11311131%%%%%%%%%%  recheck here....  (gm) 
    1132 %\begin{flalign*} 
    1133 =& \int \limits_D  -\frac{1} {e_3 } \textbf{k} \cdot \nabla \times  
    1134    \Bigl[ \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right)  \Bigr]\;dv \\  
    1135 %\end{flalign*} 
    1136 %\begin{flalign*} 
     1132\begin{flalign*} 
     1133= \int \limits_D  -\frac{1} {e_3 } \textbf{k} \cdot \nabla \times  
     1134   \Bigl[ \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right)  \Bigr]\;dv &&& \\  
     1135\end{flalign*} 
     1136\begin{flalign*} 
    11371137\equiv& \sum\limits_{i,j} 
    11381138   \left\{ 
    1139      \delta_{i+1/2} \left[  \frac {e_{2v}} {e_{1v}\,e_{3v}}  \delta_i \left[ A_f^{\,lm} e_{3f} \zeta  \right]  \right] 
    1140    + \delta_{j+1/2} \left[  \frac {e_{1u}} {e_{2u}\,e_{3u}}  \delta_j \left[ A_f^{\,lm} e_{3f} \zeta  \right]  \right] 
    1141    \right\}     \\  
     1139   \delta_{i+1/2}  
     1140   \left[  
     1141   \frac {e_{2v}} {e_{1v}\,e_{3v}}  \delta_i 
     1142      \left[ A_f^{\,lm} e_{3f} \zeta  \right] 
     1143    \right] 
     1144   + \delta_{j+1/2}  
     1145   \left[  
     1146   \frac {e_{1u}} {e_{2u}\,e_{3u}} \delta_j  
     1147      \left[ A_f^{\,lm} e_{3f} \zeta  \right] 
     1148   \right] 
     1149   \right\}  
     1150   && \\  
    11421151% 
    11431152\intertext{Using \eqref{DOM_di_adj}, it follows:} 
    … …  
    11451154\equiv& \sum\limits_{i,j,k}  
    11461155   -\,\left\{ 
    1147       \frac{e_{2v}} {e_{1v}\,e_{3v}}  \delta_i  \left[ A_f^{\,lm} e_{3f} \zeta  \right]\;\delta_i \left[ 1\right] 
    1148     + \frac{e_{1u}} {e_{2u}\,e_{3u}}  \delta_j  \left[ A_f^{\,lm} e_{3f} \zeta  \right]\;\delta_j \left[ 1\right] 
     1156      \frac{e_{2v}} {e_{1v}\,e_{3v}}  \delta_i 
     1157      \left[ A_f^{\,lm} e_{3f} \zeta  \right]\;\delta_i \left[ 1\right] 
     1158   + \frac{e_{1u}} {e_{2u}\,e_{3u}} \delta_j  
     1159      \left[ A_f^{\,lm} e_{3f} \zeta  \right]\;\delta_j \left[ 1\right] 
    11491160   \right\} \quad \equiv 0  
    1150     \\  
     1161   && \\  
    11511162\end{flalign*} 
    11521163 
    … …  
    11561167\subsection{Dissipation of Horizontal Kinetic Energy} 
    11571168\label{Apdx_C.3.2} 
     1169 
    11581170 
    11591171The lateral momentum diffusion term dissipates the horizontal kinetic energy: 
    … …  
    12091221\label{Apdx_C.3.3} 
    12101222 
     1223 
    12111224The lateral momentum diffusion term dissipates the enstrophy when the eddy  
    12121225coefficients are horizontally uniform: 
    … …  
    12151228   \left[   \nabla_h \left( A^{\,lm}\;\chi  \right) 
    12161229          - \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right)   \right]\;dv &&&\\ 
    1217 &\quad = A^{\,lm} \int \limits_D \zeta \textbf{k} \cdot \nabla \times  
     1230&= A^{\,lm} \int \limits_D \zeta \textbf{k} \cdot \nabla \times  
    12181231   \left[    \nabla_h \times \left( \zeta \; \textbf{k} \right)   \right]\;dv &&&\\ 
    1219 &\quad \equiv A^{\,lm} \sum\limits_{i,j,k}  \zeta \;e_{3f}  
     1232&\equiv A^{\,lm} \sum\limits_{i,j,k}  \zeta \;e_{3f}  
    12201233   \left\{     \delta_{i+1/2} \left[  \frac{e_{2v}} {e_{1v}\,e_{3v}} \delta_i \left[ e_{3f} \zeta  \right]   \right] 
    12211234             + \delta_{j+1/2} \left[  \frac{e_{1u}} {e_{2u}\,e_{3u}} \delta_j \left[ e_{3f} \zeta  \right]   \right]      \right\}   &&&\\  
    … …  
    12231236\intertext{Using \eqref{DOM_di_adj}, it follows:} 
    12241237% 
    1225 &\quad \equiv  - A^{\,lm} \sum\limits_{i,j,k}  
     1238&\equiv  - A^{\,lm} \sum\limits_{i,j,k}  
    12261239   \left\{    \left(  \frac{1} {e_{1v}\,e_{3v}}  \delta_i \left[ e_{3f} \zeta  \right]  \right)^2   b_v 
    1227             + \left(  \frac{1} {e_{2u}\,e_{3u}}  \delta_j \left[ e_{3f} \zeta  \right] \right)^2   b_u  \right\}  \quad \leq \;0    &&&\\ 
     1240            + \left(  \frac{1} {e_{2u}\,e_{3u}}  \delta_j \left[ e_{3f} \zeta  \right] \right)^2   b_u  \right\}      &&&\\ 
     1241& \leq \;0       &&&\\  
    12281242\end{flalign*} 
    12291243 
    … …  
    12361250When the horizontal divergence of the horizontal diffusion of momentum  
    12371251(discrete sense) is taken, the term associated with the vertical curl of the  
    1238 vorticity is zero locally, due to \eqref{Eq_DOM_div_curl}.  
    1239 The resulting term conserves the $\chi$ and dissipates $\chi^2$  
    1240 when the eddy coefficients are horizontally uniform. 
     1252vorticity is zero locally, due to (!!! II.1.8  !!!!!). The resulting term conserves the  
     1253$\chi$ and dissipates $\chi^2$ when the eddy coefficients are  
     1254horizontally uniform. 
    12411255\begin{flalign*} 
    12421256& \int\limits_D  \nabla_h \cdot  
    12431257   \Bigl[     \nabla_h \left( A^{\,lm}\;\chi \right) 
    12441258             - \nabla_h \times \left( A^{\,lm}\;\zeta \;\textbf{k} \right)    \Bigr]  dv 
    1245 = \int\limits_D  \nabla_h \cdot \nabla_h \left( A^{\,lm}\;\chi  \right)   dv   \\ 
     1259= \int\limits_D  \nabla_h \cdot \nabla_h \left( A^{\,lm}\;\chi  \right)   dv   &&&\\ 
    12461260% 
    12471261&\equiv \sum\limits_{i,j,k}  
    12481262   \left\{   \delta_i \left[ A_u^{\,lm} \frac{e_{2u}\,e_{3u}} {e_{1u}}  \delta_{i+1/2} \left[ \chi \right]  \right] 
    1249            + \delta_j \left[ A_v^{\,lm} \frac{e_{1v}\,e_{3v}} {e_{2v}}  \delta_{j+1/2} \left[ \chi \right]  \right]    \right\}    \\  
     1263           + \delta_j \left[ A_v^{\,lm} \frac{e_{1v}\,e_{3v}} {e_{2v}}  \delta_{j+1/2} \left[ \chi \right]  \right]    \right\}    &&&\\  
    12501264% 
    12511265\intertext{Using \eqref{DOM_di_adj}, it follows:} 
    … …  
    12531267&\equiv \sum\limits_{i,j,k}  
    12541268   - \left\{   \frac{e_{2u}\,e_{3u}} {e_{1u}}  A_u^{\,lm} \delta_{i+1/2} \left[ \chi \right] \delta_{i+1/2} \left[ 1 \right]  
    1255              + \frac{e_{1v}\,e_{3v}} {e_{2v}}  A_v^{\,lm} \delta_{j+1/2} \left[ \chi \right] \delta_{j+1/2} \left[ 1 \right]    \right\}  
    1256    \quad \equiv 0      \\  
     1269             + \frac{e_{1v}\,e_{3v}}  {e_{2v}}  A_v^{\,lm} \delta_{j+1/2} \left[ \chi \right] \delta_{j+1/2} \left[ 1 \right]    \right\}  
     1270   \qquad \equiv 0     &&& \\  
    12571271\end{flalign*} 
    12581272 
    … …  
    12671281   \left[    \nabla_h              \left( A^{\,lm}\;\chi                    \right) 
    12681282           - \nabla_h   \times  \left( A^{\,lm}\;\zeta \;\textbf{k} \right)    \right]\;  dv 
    1269  = A^{\,lm}   \int\limits_D \chi \;\nabla_h \cdot \nabla_h \left( \chi \right)\;  dv    \\  
     1283 = A^{\,lm}   \int\limits_D \chi \;\nabla_h \cdot \nabla_h \left( \chi \right)\;  dv    &&&\\  
    12701284% 
    12711285&\equiv A^{\,lm}  \sum\limits_{i,j,k}  \frac{1} {e_{1t}\,e_{2t}\,e_{3t}}  \chi  
    … …  
    12731287      \delta_i  \left[   \frac{e_{2u}\,e_{3u}} {e_{1u}}  \delta_{i+1/2} \left[ \chi \right]   \right] 
    12741288   + \delta_j  \left[   \frac{e_{1v}\,e_{3v}} {e_{2v}}   \delta_{j+1/2} \left[ \chi \right]   \right] 
    1275    \right\} \;   e_{1t}\,e_{2t}\,e_{3t}    \\  
     1289   \right\} \;   e_{1t}\,e_{2t}\,e_{3t}    &&&\\  
    12761290% 
    12771291\intertext{Using \eqref{DOM_di_adj}, it turns out to be:} 
    … …  
    12791293&\equiv - A^{\,lm} \sum\limits_{i,j,k} 
    12801294   \left\{    \left(  \frac{1} {e_{1u}}  \delta_{i+1/2}  \left[ \chi \right]  \right)^2  b_u 
    1281             + \left(  \frac{1} {e_{2v}}  \delta_{j+1/2}  \left[ \chi \right]  \right)^2  b_v    \right\}     
    1282 \quad \leq 0             \\ 
     1295                 + \left(  \frac{1} {e_{2v}}  \delta_{j+1/2}  \left[ \chi \right]  \right)^2  b_v    \right\} \;    &&&\\ 
     1296% 
     1297&\leq 0              &&&\\ 
    12831298\end{flalign*} 
    12841299 
    … …  
    12881303\section{Conservation Properties on Vertical Momentum Physics} 
    12891304\label{Apdx_C_4} 
     1305 
    12901306 
    12911307As for the lateral momentum physics, the continuous form of the vertical diffusion  
    … …  
    13031319   \left(   \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k}   \right)\; dv    \quad &\leq 0     \\ 
    13041320\end{align*} 
    1305  
    13061321The first property is obvious. The second results from: 
     1322 
    13071323\begin{flalign*} 
    13081324\int\limits_D  
    … …  
    13431359   e_{1f}\,e_{2f}\,e_{3f} \; \equiv 0   && \\ 
    13441360\end{flalign*} 
    1345  
    13461361If the vertical diffusion coefficient is uniform over the whole domain, the  
    13471362enstrophy is dissipated, $i.e.$ 
    … …  
    13511366      \left( \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k}   \right)   \right)\; dv = 0   &&&\\ 
    13521367\end{flalign*} 
    1353  
    13541368This property is only satisfied in $z$-coordinates: 
     1369 
    13551370\begin{flalign*} 
    13561371\int\limits_D \zeta \, \textbf{k} \cdot \nabla \times  
    … …  
    14621477 
    14631478The numerical schemes used for tracer subgridscale physics are written such  
    1464 that the heat and salt contents are conserved (equations in flux form).  
    1465 Since a flux form is used to compute the temperature and salinity,  
    1466 the quadratic form of these quantities ($i.e.$ their variance) globally tends to diminish.  
    1467 As for the advection term, there is conservation of mass only if the Equation Of Seawater is linear.  
     1479that the heat and salt contents are conserved (equations in flux form, second  
     1480order centered finite differences). Since a flux form is used to compute the  
     1481temperature and salinity, the quadratic form of these quantities (i.e. their variance)  
     1482globally tends to diminish. As for the advection term, there is generally no strict  
     1483conservation of mass, even if in practice the mass is conserved to a very high  
     1484accuracy.  
    14681485 
    14691486% ------------------------------------------------------------------------------------------------------------- 

  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Annex_D.tex

    r6617 r6625  
    120120\hline 
    121121public  \par or  \par module variable&  
    122 \textbf{m n} \par \textit{but not} \par \textbf{nn\_ np\_}&  
     122\textbf{m n} \par \textit{but not} \par \textbf{nn\_}&  
    123123\textbf{a b e f g h o q r} \par \textbf{t} \textit{to} \textbf{x} \par but not \par \textbf{fs rn\_}&  
    124124\textbf{l} \par \textit{but not} \par \textbf{lp ld} \par \textbf{ ll ln\_}&  
    … …  
    156156\hline 
    157157parameter&  
    158 \textbf{jp np\_}&  
     158\textbf{jp}&  
    159159\textbf{pp}&  
    160160\textbf{lp}&  
    … …  
    190190%-------------------------------------------------------------------------------------------------------------- 
    191191 
    192 N.B.   Parameter here, in not only parameter in the \textsc{Fortran} acceptation, it is also used for code variables  
    193 that are read in namelist and should never been modified during a simulation.  
    194 It is the case, for example, for the size of a domain (jpi,jpj,jpk). 
    195  
    196192\newpage 
    197193% ================================================================ 

  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Chap_DIA.tex

    r6617 r6625  
    22% Chapter I/O & Diagnostics 
    33% ================================================================ 
    4 \chapter{Output and Diagnostics (IOM, DIA, TRD, FLO)} 
     4\chapter{Ouput and Diagnostics (IOM, DIA, TRD, FLO)} 
    55\label{DIA} 
    66\minitoc 
    77 
    88\newpage 
    9 $\ $\newline    % force a new line 
     9$\ $\newline    % force a new ligne 
    1010 
    1111% ================================================================ 
    … …  
    4848 
    4949 
    50 Since version 3.2, iomput is the NEMO output interface of choice.  
    51 It has been designed to be simple to use, flexible and efficient.  
    52 The two main purposes of iomput are:  
     50Since version 3.2, iomput is the NEMO output interface of choice. It has been designed to be simple to use, flexible and efficient. The two main purposes of iomput are:  
    5351\begin{enumerate} 
    5452\item The complete and flexible control of the output files through external XML files adapted by the user from standard templates.  
    … …  
    11181116% ------------------------------------------------------------------------------------------------------------- 
    11191117\section[Tracer/Dynamics Trends (TRD)] 
    1120                   {Tracer/Dynamics Trends  (\ngn{namtrd})} 
     1118                  {Tracer/Dynamics Trends  (\key{trdtra}, \key{trddyn},    \\  
     1119                                                             \key{trddvor}, \key{trdmld})} 
    11211120\label{DIA_trd} 
    11221121 
    … …  
    11251124%------------------------------------------------------------------------------------------------------------- 
    11261125 
    1127 Each trend of the dynamics and/or temperature and salinity time evolution equations  
    1128 can be send to \mdl{trddyn} and/or \mdl{trdtra} modules (see TRD directory) just after their computation  
    1129 ($i.e.$ at the end of each $dyn\cdots.F90$ and/or $tra\cdots.F90$ routines).  
    1130 This capability is controlled by options offered in \ngn{namtrd} namelist.  
    1131 Note that the output are done with xIOS, and therefore the \key{IOM} is required. 
    1132  
    1133 What is done depends on the \ngn{namtrd} logical set to \textit{true}: 
     1126When \key{trddyn} and/or \key{trddyn} CPP variables are defined, each  
     1127trend of the dynamics and/or temperature and salinity time evolution equations  
     1128is stored in three-dimensional arrays just after their computation ($i.e.$ at the end  
     1129of each $dyn\cdots.F90$ and/or $tra\cdots.F90$ routines). Options are defined by 
     1130\ngn{namtrd} namelist variables. These trends are then  
     1131used in \mdl{trdmod} (see TRD directory) every \textit{nn\_trd } time-steps. 
     1132 
     1133What is done depends on the CPP keys defined: 
    11341134\begin{description} 
    1135 \item[\np{ln\_glo\_trd}] : at each \np{nn\_trd} time-step a check of the basin averaged properties  
    1136 of the momentum and tracer equations is performed. This also includes a check of $T^2$, $S^2$,  
    1137 $\tfrac{1}{2} (u^2+v2)$, and potential energy time evolution equations properties ;  
    1138 \item[\np{ln\_dyn\_trd}] : each 3D trend of the evolution of the two momentum components is output ;  
    1139 \item[\np{ln\_dyn\_mxl}] : each 3D trend of the evolution of the two momentum components averaged  
    1140                            over the mixed layer is output  ;  
    1141 \item[\np{ln\_vor\_trd}] : a vertical summation of the moment tendencies is performed,  
    1142                            then the curl is computed to obtain the barotropic vorticity tendencies which are output ; 
    1143 \item[\np{ln\_KE\_trd}]  : each 3D trend of the Kinetic Energy equation is output ; 
    1144 \item[\np{ln\_tra\_trd}] : each 3D trend of the evolution of temperature and salinity is output ; 
    1145 \item[\np{ln\_tra\_mxl}] : each 2D trend of the evolution of temperature and salinity averaged  
    1146                            over the mixed layer is output ; 
     1135\item[\key{trddyn}, \key{trdtra}] : a check of the basin averaged properties of the momentum  
     1136and/or tracer equations is performed ;  
     1137\item[\key{trdvor}] : a vertical summation of the moment tendencies is performed,  
     1138then the curl is computed to obtain the barotropic vorticity tendencies which are output ; 
     1139\item[\key{trdmld}] : output of the tracer tendencies averaged vertically   
     1140either over the mixed layer (\np{nn\_ctls}=0),  
     1141or       over a fixed number of model levels (\np{nn\_ctls}$>$1 provides the number of level),  
     1142or       over a spatially varying but temporally fixed number of levels (typically the base  
     1143of the winter mixed layer) read in \ifile{ctlsurf\_idx} (\np{nn\_ctls}=1) ; 
    11471144\end{description} 
     1145 
     1146The units in the output file can be changed using the \np{nn\_ucf} namelist parameter.  
     1147For example, in case of salinity tendency the units are given by PSU/s/\np{nn\_ucf}. 
     1148Setting \np{nn\_ucf}=86400 ($i.e.$ the number of second in a day) provides the tendencies in PSU/d. 
     1149 
     1150When \key{trdmld} is defined, two time averaging procedure are proposed. 
     1151Setting \np{ln\_trdmld\_instant} to \textit{true}, a simple time averaging is performed,  
     1152so that the resulting tendency is the contribution to the change of a quantity between  
     1153the two instantaneous values taken at the extremities of the time averaging period. 
     1154Setting \np{ln\_trdmld\_instant} to \textit{false}, a double time averaging is performed,  
     1155so that the resulting tendency is the contribution to the change of a quantity between  
     1156two \textit{time mean} values. The later option requires the use of an extra file, \ifile{restart\_mld}   
     1157(\np{ln\_trdmld\_restart}=true), to restart a run. 
     1158 
    11481159 
    11491160Note that the mixed layer tendency diagnostic can also be used on biogeochemical models  
    11501161via the \key{trdtrc} and \key{trdmld\_trc} CPP keys. 
    1151  
    1152 \textbf{Note that} in the current version (v3.6), many changes has been introduced but not fully tested.  
    1153 In particular, options associated with \np{ln\_dyn\_mxl}, \np{ln\_vor\_trd}, and \np{ln\_tra\_mxl}  
    1154 are not working, and none of the option have been tested with variable volume ($i.e.$ \key{vvl} defined). 
    1155  
    11561162 
    11571163% ------------------------------------------------------------------------------------------------------------- 
    … …  
    12741280\label{DIA_diag_harm} 
    12751281 
     1282A module is available to compute the amplitude and phase for tidal waves.  
     1283This diagnostic is actived with \key{diaharm}. 
     1284 
    12761285%------------------------------------------namdia_harm---------------------------------------------------- 
    12771286\namdisplay{namdia_harm} 
    12781287%---------------------------------------------------------------------------------------------------------- 
    12791288 
    1280 A module is available to compute the amplitude and phase of tidal waves.  
    1281 This on-line Harmonic analysis is actived with \key{diaharm}. 
    1282 Some parameters are available in namelist \ngn{namdia\_harm} : 
    1283  
    1284 - \np{nit000\_han} is the first time step used for harmonic analysis 
    1285  
    1286 - \np{nitend\_han} is the last time step used for harmonic analysis 
    1287  
    1288 - \np{nstep\_han} is the time step frequency for harmonic analysis 
    1289  
    1290 - \np{nb\_ana} is the number of harmonics to analyse 
    1291  
    1292 - \np{tname} is an array with names of tidal constituents to analyse 
    1293  
    1294 \np{nit000\_han} and \np{nitend\_han} must be between \np{nit000} and \np{nitend} of the simulation. 
     1289Concerning the on-line Harmonic analysis, some parameters are available in namelist 
     1290\ngn{namdia\_harm} : 
     1291 
     1292- \texttt{nit000\_han} is the first time step used for harmonic analysis 
     1293 
     1294- \texttt{nitend\_han} is the last time step used for harmonic analysis 
     1295 
     1296- \texttt{nstep\_han} is the time step frequency for harmonic analysis 
     1297 
     1298- \texttt{nb\_ana} is the number of harmonics to analyse 
     1299 
     1300- \texttt{tname} is an array with names of tidal constituents to analyse 
     1301 
     1302\texttt{nit000\_han} and \texttt{nitend\_han} must be between \texttt{nit000} and \texttt{nitend} of the simulation. 
    12951303The restart capability is not implemented. 
    12961304 
    1297 The Harmonic analysis solve the following equation: 
     1305The Harmonic analysis solve this equation: 
    12981306\begin{equation} 
    12991307h_{i} - A_{0} + \sum^{nb\_ana}_{j=1}[A_{j}cos(\nu_{j}t_{j}-\phi_{j})] = e_{i} 
    … …  
    13161324\label{DIA_diag_dct} 
    13171325 
    1318 A module is available to compute the transport of volume, heat and salt through sections.  
    1319 This diagnostic is actived with \key{diadct}. 
     1326A module is available to compute the transport of volume, heat and salt through sections. This diagnostic 
     1327is actived with \key{diadct}. 
    13201328 
    13211329Each section is defined by the coordinates of its 2 extremities. The pathways between them are contructed 
    … …  
    13391347%------------------------------------------------------------------------------------------------------------- 
    13401348 
    1341 \np{nn\_dct}: frequency of instantaneous transports computing 
    1342  
    1343 \np{nn\_dctwri}: frequency of writing ( mean of instantaneous transports ) 
    1344  
    1345 \np{nn\_debug}: debugging of the section 
     1349\texttt{nn\_dct}: frequency of instantaneous transports computing 
     1350 
     1351\texttt{nn\_dctwri}: frequency of writing ( mean of instantaneous transports ) 
     1352 
     1353\texttt{nn\_debug}: debugging of the section 
    13461354 
    13471355\subsubsection{ To create a binary file containing the pathway of each section } 
    … …  
    14741482the \key{diahth} CPP key:  
    14751483 
    1476 - the mixed layer depth (based on a density criterion \citep{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth}) 
     1484- the mixed layer depth (based on a density criterion, \citet{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth}) 
    14771485 
    14781486- the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diahth}) 

  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Chap_DOM.tex

    r6617 r6625  
    11% ================================================================ 
    2 % Chapter 2 ——— Space and Time Domain (DOM) 
     2% Chapter 2 � Space and Time Domain (DOM) 
    33% ================================================================ 
    44\chapter{Space Domain (DOM) } 
    … …  
    138138and $f$-points, and its divergence defined at $t$-points: 
    139139\begin{eqnarray}  \label{Eq_DOM_curl} 
    140  \nabla \times {\rm{\bf A}}\equiv & 
     140 \nabla \times {\rm {\bf A}}\equiv & 
    141141      \frac{1}{e_{2v}\,e_{3vw} } \ \left( \delta_{j +1/2} \left[e_{3w}\,a_3 \right] -\delta_{k+1/2} \left[e_{2v} \,a_2 \right] \right)  &\ \mathbf{i} \\  
    142142 +& \frac{1}{e_{2u}\,e_{3uw}} \ \left( \delta_{k+1/2} \left[e_{1u}\,a_1  \right] -\delta_{i +1/2} \left[e_{3w}\,a_3 \right] \right)  &\ \mathbf{j} \\ 
    … …  
    183183Let $a$ and $b$ be two fields defined on the mesh, with value zero inside  
    184184continental area. Using integration by parts it can be shown that the differencing  
    185 operators ($\delta_i$, $\delta_j$ and $\delta_k$) are skew-symmetric linear operators,  
    186 and further that the averaging operators $\overline{\,\cdot\,}^{\,i}$,  
     185operators ($\delta_i$, $\delta_j$ and $\delta_k$) are anti-symmetric linear  
     186operators, and further that the averaging operators $\overline{\,\cdot\,}^{\,i}$,  
    187187$\overline{\,\cdot\,}^{\,k}$ and $\overline{\,\cdot\,}^{\,k}$) are symmetric linear  
    188188operators, $i.e.$ 
    … …  
    364364For both grids here,  the same $w$-point depth has been chosen but in (a) the  
    365365$t$-points are set half way between $w$-points while in (b) they are defined from  
    366 an analytical function: $z(k)=5\,(k-1/2)^3 - 45\,(k-1/2)^2 + 140\,(k-1/2) - 150$.  
     366an analytical function: $z(k)=5\,(i-1/2)^3 - 45\,(i-1/2)^2 + 140\,(i-1/2) - 150$.  
    367367Note the resulting difference between the value of the grid-size $\Delta_k$ and  
    368368those of the scale factor $e_k$. } 
    … …  
    425425 
    426426The choice of the grid must be consistent with the boundary conditions specified  
    427 by \np{jperio}, a parameter found in \ngn{namcfg} namelist (see {\S\ref{LBC}). 
     427by the parameter \np{jperio} (see {\S\ref{LBC}). 
    428428 
    429429% ------------------------------------------------------------------------------------------------------------- 
    … …  
    481481%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    482482 
    483 The choice of a vertical coordinate, even if it is made through \ngn{namzgr} namelist parameters,  
     483The choice of a vertical coordinate, even if it is made through a namelist parameter,  
    484484must be done once of all at the beginning of an experiment. It is not intended as an  
    485485option which can be enabled or disabled in the middle of an experiment. Three main  
    … …  
    494494bathymetry or $s$-coordinate (hybrid and partial step coordinates have not  
    495495yet been tested in NEMO v2.3). If using $z$-coordinate with partial step bathymetry 
    496 (\np{ln\_zps}~=~true), ocean cavity beneath ice shelves can be open (\np{ln\_isfcav}~=~true)  
    497 and partial step are also applied at the ocean/ice shelf interface.  
     496(\np{ln\_zps}~=~true), ocean cavity beneath ice shelves can be open (\np{ln\_isfcav}~=~true). 
    498497 
    499498Contrary to the horizontal grid, the vertical grid is computed in the code and no  
    500499provision is made for reading it from a file. The only input file is the bathymetry  
    501 (in meters) (\ifile{bathy\_meter}).  
     500(in meters) (\ifile{bathy\_meter})  
    502501\footnote{N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the  
    503502\ifile{bathy\_meter} file, so that the computation of the number of wet ocean point  
    … …  
    541540 
    542541Three options are possible for defining the bathymetry, according to the  
    543 namelist variable \np{nn\_bathy} (found in \ngn{namdom} namelist):  
     542namelist variable \np{nn\_bathy}:  
    544543\begin{description} 
    545544\item[\np{nn\_bathy} = 0] a flat-bottom domain is defined. The total depth $z_w (jpk)$  
    … …  
    549548domain width at the central latitude. This is meant for the "EEL-R5" configuration,  
    550549a periodic or open boundary channel with a seamount.  
    551 \item[\np{nn\_bathy} = 1] read a bathymetry and ice shelf draft (if needed). 
    552  The \ifile{bathy\_meter} file (Netcdf format) provides the ocean depth (positive, in meters) 
    553  at each grid point of the model grid. The bathymetry is usually built by interpolating a standard bathymetry product  
     550\item[\np{nn\_bathy} = 1] read a bathymetry. The \ifile{bathy\_meter} file (Netcdf format)  
     551provides the ocean depth (positive, in meters) at each grid point of the model grid.  
     552The bathymetry is usually built by interpolating a standard bathymetry product  
    554553($e.g.$ ETOPO2) onto the horizontal ocean mesh. Defining the bathymetry also  
    555554defines the coastline: where the bathymetry is zero, no model levels are defined  
    556555(all levels are masked). 
    557  
    558 The \ifile{isfdraft\_meter} file (Netcdf format) provides the ice shelf draft (positive, in meters) 
    559  at each grid point of the model grid. This file is only needed if \np{ln\_isfcav}~=~true.  
    560 Defining the ice shelf draft will also define the ice shelf edge and the grounding line position. 
    561556\end{description} 
    562557 
    … …  
    615610(Fig.~\ref{Fig_zgr}). 
    616611 
    617 If the ice shelf cavities are opened (\np{ln\_isfcav}=~true~}), the definition of $z_0$ is the same.  
    618 However, definition of $e_3^0$ at $t$- and $w$-points is respectively changed to: 
    619 \begin{equation} \label{DOM_zgr_ana} 
    620 \begin{split} 
    621  e_3^T(k) &= z_W (k+1) - z_W (k)   \\ 
    622  e_3^W(k) &= z_T (k)   - z_T (k-1) \\ 
    623 \end{split} 
    624 \end{equation} 
    625 This formulation decrease the self-generated circulation into the ice shelf cavity  
    626 (which can, in extreme case, leads to blow up).\\ 
    627  
    628   
    629612The most used vertical grid for ORCA2 has $10~m$ ($500~m)$ resolution in the  
    630613surface (bottom) layers and a depth which varies from 0 at the sea surface to a  
    … …  
    738721usually 10\%, of the default thickness $e_{3t}(jk)$). 
    739722 
    740 \gmcomment{ \colorbox{yellow}{Add a figure here of pstep especially at last ocean level } } 
     723 \colorbox{yellow}{Add a figure here of pstep especially at last ocean level } 
    741724 
    742725% ------------------------------------------------------------------------------------------------------------- 
    … …  
    877860gives the number of ocean levels ($i.e.$ those that are not masked) at each  
    878861$t$-point. mbathy is computed from the meter bathymetry using the definiton of  
    879 gdept as the number of $t$-points which gdept $\leq$ bathy. 
     862gdept as the number of $t$-points which gdept $\leq$ bathy.  
    880863 
    881864Modifications of the model bathymetry are performed in the \textit{bat\_ctl}  
    882865routine (see \mdl{domzgr} module) after mbathy is computed. Isolated grid points  
    883 that do not communicate with another ocean point at the same level are eliminated.\\ 
    884  
    885 As for the representation of bathymetry, a 2D integer array, misfdep, is created.  
    886 misfdep defines the level of the first wet $t$-point. All the cells between $k=1$ and $misfdep(i,j)-1$ are masked.  
    887 By default, misfdep(:,:)=1 and no cells are masked. 
    888  
    889 In case of ice shelf cavities (\np{ln\_isfcav}~=~true), modifications of the model bathymetry and ice shelf draft in  
    890 the cavities are performed through the \textit{zgr\_isf} routine. The compatibility between ice shelf draft and bathymetry is checked:  
    891 if only one cell on the water column is opened at $t$-, $u$- or $v$-points, the bathymetry or the ice shelf draft is dug to have a 2-level water column  
    892 (i.e. two unmasked levels). If the incompatibility is too strong (i.e. need to dig more than one cell), the entire water column is masked.\\  
     866that do not communicate with another ocean point at the same level are eliminated. 
    893867 
    894868From the \textit{mbathy} array, the mask fields are defined as follows: 
    895869\begin{align*} 
    896 tmask(i,j,k) &= \begin{cases}   \; 0&   \text{ if $k < misfdep(i,j) $ } \\ 
    897                                 \; 1&   \text{ if $misfdep(i,j) \leq k\leq mbathy(i,j)$  }    \\ 
    898                                 \; 0&   \text{ if $k > mbathy(i,j)$  }    \end{cases}     \\ 
     870tmask(i,j,k) &= \begin{cases}   \; 1&   \text{ if $k\leq mbathy(i,j)$  }    \\ 
     871                                                \; 0&   \text{ if $k\leq mbathy(i,j)$  }    \end{cases}     \\ 
    899872umask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i+1,j,k)   \\ 
    900873vmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i,j+1,k)   \\ 
    901874fmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i+1,j,k)   \\ 
    902                    & \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k) \\ 
    903 wmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i,j,k-1) \text{ with } wmask(i,j,1) = tmask(i,j,1)  
     875                   & \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k) 
    904876\end{align*} 
    905877 
    906 Note, wmask is now defined. It allows, in case of ice shelves,  
    907 to deal with the top boundary (ice shelf/ocean interface) exactly in the same way as for the bottom boundary.  
    908  
    909 The specification of closed lateral boundaries requires that at least the first and last  
     878Note that \textit{wmask} is not defined as it is exactly equal to \textit{tmask} with  
     879the numerical indexing used (\S~\ref{DOM_Num_Index}). Moreover, the  
     880specification of closed lateral boundaries requires that at least the first and last  
    910881rows and columns of the \textit{mbathy} array are set to zero. In the particular  
    911882case of an east-west cyclical boundary condition, \textit{mbathy} has its last  

  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Chap_DYN.tex

    r6617 r6625  
    11% ================================================================ 
    2 % Chapter ——— Ocean Dynamics (DYN) 
     2% Chapter � Ocean Dynamics (DYN) 
    33% ================================================================ 
    44\chapter{Ocean Dynamics (DYN)} 
    55\label{DYN} 
    66\minitoc 
     7 
     8% add a figure for  dynvor ens, ene latices 
    79 
    810%\vspace{2.cm} 
    … …  
    163165%------------------------------------------------------------------------------------------------------------- 
    164166 
    165 The vector invariant form of the momentum equations (\np{ln\_dynhpg\_vec}~=~true) is the one most  
    166 often used in applications of the \NEMO ocean model. The flux form option (\np{ln\_dynhpg\_vec}~=false) 
    167 (see next section) has been present since version $2$.  
    168 Options are defined through the \ngn{namdyn\_adv} namelist variables. 
    169 Coriolis and momentum advection terms are evaluated using a leapfrog scheme,  
    170 $i.e.$ the velocity appearing in these expressions is centred in time (\textit{now} velocity).  
     167The vector invariant form of the momentum equations is the one most  
     168often used in applications of the \NEMO ocean model. The flux form option  
     169(see next section) has been present since version $2$. Options are defined 
     170through the \ngn{namdyn\_adv} namelist variables 
     171Coriolis and momentum advection terms are evaluated using a leapfrog  
     172scheme, $i.e.$ the velocity appearing in these expressions is centred in  
     173time (\textit{now} velocity).  
    171174At the lateral boundaries either free slip, no slip or partial slip boundary  
    172175conditions are applied following Chap.\ref{LBC}. 
    … …  
    300303%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    301304 
    302 A key point in \eqref{Eq_een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made.  
    303 It uses the sum of masked t-point vertical scale factor divided either  
    304 by the sum of the four t-point masks (\np{ln\_dynvor\_een\_old}~=~false),  
    305 or  just by $4$ (\np{ln\_dynvor\_een\_old}~=~true). 
    306 The latter case preserves the continuity of $e_{3f}$ when one or more of the neighbouring $e_{3t}$  
    307 tends to zero and extends by continuity the value of $e_{3f}$ into the land areas.  
    308 This case introduces a sub-grid-scale topography at f-points (with a systematic reduction of $e_{3f}$  
    309 when a model level intercept the bathymetry) that tends to reinforce the topostrophy of the flow  
    310 ($i.e.$ the tendency of the flow to follow the isobaths) \citep{Penduff_al_OS07}.  
     305Note that a key point in \eqref{Eq_een_e3f} is that the averaging in the \textbf{i}- and  
     306\textbf{j}- directions uses the masked vertical scale factor but is always divided by  
     307$4$, not by the sum of the masks at the four $T$-points. This preserves the continuity of  
     308$e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and  
     309extends by continuity the value of $e_{3f}$ into the land areas. This feature is essential for  
     310the $z$-coordinate with partial steps. 
    311311 
    312312Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as  
    … …  
    374374\end{aligned}         \right. 
    375375\end{equation}  
    376 When \np{ln\_dynzad\_zts}~=~\textit{true}, a split-explicit time stepping with 5 sub-timesteps is used  
    377 on the vertical advection term. 
    378 This option can be useful when the value of the timestep is limited by vertical advection \citep{Lemarie_OM2015}.  
    379 Note that in this case, a similar split-explicit time stepping should be used on  
    380 vertical advection of tracer to ensure a better stability,  
    381 an option which is only available with a TVD scheme (see \np{ln\_traadv\_tvd\_zts} in \S\ref{TRA_adv_tvd}). 
    382  
    383376 
    384377% ================================================================ 
    … …  
    498491those in the centred second order method. As the scheme already includes  
    499492a diffusion component, it can be used without explicit  lateral diffusion on momentum  
    500 ($i.e.$ setting both \np{ln\_dynldf\_lap} and \np{ln\_dynldf\_bilap} to \textit{false}),  
    501 and it is recommended to do so. 
     493($i.e.$ \np{ln\_dynldf\_lap}=\np{ln\_dynldf\_bilap}=false), and it is recommended to do so. 
    502494 
    503495The UBS scheme is not used in all directions. In the vertical, the centred $2^{nd}$  
    … …  
    637629($e_{3w}$). 
    638630  
     631$\bullet$ Traditional coding with adaptation for ice shelf cavities (\np{ln\_dynhpg\_isf}=true). 
     632This scheme need the activation of ice shelf cavities (\np{ln\_isfcav}=true). 
     633 
    639634$\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}=true) 
    640635 
    … …  
    651646pressure Jacobian method is used to solve the horizontal pressure gradient. This method can provide 
    652647a more accurate calculation of the horizontal pressure gradient than the standard scheme. 
    653  
    654 \subsection{Ice shelf cavity} 
    655 \label{DYN_hpg_isf} 
    656 Beneath an ice shelf, the total pressure gradient is the sum of the pressure gradient due to the ice shelf load and 
    657  the pressure gradient due to the ocean load. If cavities are present (\np{ln\_isfcav}~=~true) these two terms can be 
    658  calculated by setting \np{ln\_dynhpg\_isf}~=~true. No other scheme is working with ice shelves.\\ 
    659  
    660 $\bullet$ The main hypothesis to compute the ice shelf load is that the ice shelf is in isostatic equilibrium. 
    661  The top pressure is computed integrating a reference density profile (prescribed as density of a water at 34.4  
    662 PSU and -1.9$\degres C$) from the sea surface to the ice shelf base, which corresponds to the load of the water 
    663 column in which the ice shelf is floatting. This top pressure is constant over time. A detailed description of  
    664 this method is described in \citet{Losch2008}.\\ 
    665  
    666 $\bullet$ The ocean load is computed using the expression \eqref{Eq_dynhpg_sco} described in \ref{DYN_hpg_sco}.  
    667 A treatment of the top and bottom partial cells similar to the one described in \ref{DYN_hpg_zps} is done  
    668 to reduce the residual circulation generated by the top partial cell.  
    669648 
    670649%-------------------------------------------------------------------------------------------------------------- 
    … …  
    739718$\ $\newline      %force an empty line 
    740719 
     720%%% 
    741721Options are defined through the \ngn{namdyn\_spg} namelist variables. 
    742 The surface pressure gradient term is related to the representation of the free surface (\S\ref{PE_hor_pg}).  
    743 The main distinction is between the fixed volume case (linear free surface) and the variable volume case  
    744 (nonlinear free surface, \key{vvl} is defined). In the linear free surface case (\S\ref{PE_free_surface})  
    745 the vertical scale factors $e_{3}$ are fixed in time, while they are time-dependent in the nonlinear case  
    746 (\S\ref{PE_free_surface}).  
    747 With both linear and nonlinear free surface, external gravity waves are allowed in the equations,  
    748 which imposes a very small time step when an explicit time stepping is used.  
    749 Two methods are proposed to allow a longer time step for the three-dimensional equations:  
    750 the filtered free surface, which is a modification of the continuous equations (see \eqref{Eq_PE_flt}),  
    751 and the split-explicit free surface described below.  
    752 The extra term introduced in the filtered method is calculated implicitly,  
    753 so that the update of the next velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. 
     722The surface pressure gradient term is related to the representation of the free surface (\S\ref{PE_hor_pg}). The main distinction is between the fixed volume case (linear free surface) and the variable volume case (nonlinear free surface, \key{vvl} is defined). In the linear free surface case (\S\ref{PE_free_surface}) the vertical scale factors $e_{3}$ are fixed in time, while they are time-dependent in the nonlinear case (\S\ref{PE_free_surface}). With both linear and nonlinear free surface, external gravity waves are allowed in the equations, which imposes a very small time step when an explicit time stepping is used. Two methods are proposed to allow a longer time step for the three-dimensional equations: the filtered free surface, which is a modification of the continuous equations (see \eqref{Eq_PE_flt}), and the split-explicit free surface described below. The extra term introduced in the filtered method is calculated implicitly, so that the update of the next velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. 
     723 
     724%%% 
    754725 
    755726 
    … …  
    765736implicitly, so that a solver is used to compute it. As a consequence the update of the $next$  
    766737velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. 
     738 
    767739 
    768740 
    … …  
    807779$\rdt_e = \rdt / nn\_baro$. This parameter can be optionally defined automatically (\np{ln\_bt\_nn\_auto}=true)  
    808780considering that the stability of the barotropic system is essentially controled by external waves propagation.  
    809 Maximum Courant number is in that case time independent, and easily computed online from the input bathymetry. 
    810 Therefore, $\rdt_e$ is adjusted so that the Maximum allowed Courant number is smaller than \np{rn\_bt\_cmax}. 
     781Maximum allowed Courant number is in that case time independent, and easily computed online from the input bathymetry. 
    811782 
    812783%%% 
    … …  
    831802Schematic of the split-explicit time stepping scheme for the external  
    832803and internal modes. Time increases to the right. In this particular exemple,  
    833 a boxcar averaging window over $nn\_baro$ barotropic time steps is used ($nn\_bt\_flt=1$) and $nn\_baro=5$. 
     804a boxcar averaging window over $nn\_baro$ barotropic time steps is used ($nn\_bt\_filt=1$) and $nn\_baro=5$. 
    834805Internal mode time steps (which are also the model time steps) are denoted  
    835806by $t-\rdt$, $t$ and $t+\rdt$. Variables with $k$ superscript refer to instantaneous barotropic variables,  
    … …  
    837808The former are used to obtain time filtered quantities at $t+\rdt$ while the latter are used to obtain time averaged  
    838809transports to advect tracers. 
    839 a) Forward time integration: \np{ln\_bt\_fw}=true, \np{ln\_bt\_av}=true.  
    840 b) Centred time integration: \np{ln\_bt\_fw}=false, \np{ln\_bt\_av}=true.  
    841 c) Forward time integration with no time filtering (POM-like scheme): \np{ln\_bt\_fw}=true, \np{ln\_bt\_av}=false. } 
     810a) Forward time integration: \np{ln\_bt\_fw}=true, \np{ln\_bt\_ave}=true.  
     811b) Centred time integration: \np{ln\_bt\_fw}=false, \np{ln\_bt\_ave}=true.  
     812c) Forward time integration with no time filtering (POM-like scheme): \np{ln\_bt\_fw}=true, \np{ln\_bt\_ave}=false. } 
    842813\end{center}    \end{figure} 
    843814%>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   > 
    … …  
    845816In the default case (\np{ln\_bt\_fw}=true), the external mode is integrated  
    846817between \textit{now} and  \textit{after} baroclinic time-steps (Fig.~\ref{Fig_DYN_dynspg_ts}a). To avoid aliasing of fast barotropic motions into three dimensional equations, time filtering is eventually applied on barotropic  
    847 quantities (\np{ln\_bt\_av}=true). In that case, the integration is extended slightly beyond  \textit{after} time step to provide time filtered quantities.  
     818quantities (\np{ln\_bt\_ave}=true). In that case, the integration is extended slightly beyond  \textit{after} time step to provide time filtered quantities.  
    848819These are used for the subsequent initialization of the barotropic mode in the following baroclinic step.  
    849820Since external mode equations written at baroclinic time steps finally follow a forward time stepping scheme,  
    … …  
    866837%%% 
    867838 
    868 One can eventually choose to feedback instantaneous values by not using any time filter (\np{ln\_bt\_av}=false).  
     839One can eventually choose to feedback instantaneous values by not using any time filter (\np{ln\_bt\_ave}=false).  
    869840In that case, external mode equations are continuous in time, ie they are not re-initialized when starting a new  
    870841sub-stepping sequence. This is the method used so far in the POM model, the stability being maintained by refreshing at (almost)  
    … …  
    11871158 
    11881159Besides the surface and bottom stresses (see the above section) which are  
    1189 introduced as boundary conditions on the vertical mixing, three other forcings  
    1190 may enter the dynamical equations by affecting the surface pressure gradient.  
    1191  
    1192 (1) When \np{ln\_apr\_dyn}~=~true (see \S\ref{SBC_apr}), the atmospheric pressure is taken  
    1193 into account when computing the surface pressure gradient. 
    1194  
    1195 (2) When \np{ln\_tide\_pot}~=~true and \key{tide} is defined (see \S\ref{SBC_tide}),  
    1196 the tidal potential is taken into account when computing the surface pressure gradient. 
    1197  
    1198 (3) When \np{nn\_ice\_embd}~=~2 and LIM or CICE is used ($i.e.$ when the sea-ice is embedded in the ocean),  
    1199 the snow-ice mass is taken into account when computing the surface pressure gradient. 
    1200  
    1201  
    1202 \gmcomment{ missing : the lateral boundary condition !!!   another external forcing 
    1203  } 
     1160introduced as boundary conditions on the vertical mixing, two other forcings  
     1161enter the dynamical equations.  
     1162 
     1163One is the effect of atmospheric pressure on the ocean dynamics. 
     1164Another forcing term is the tidal potential. 
     1165Both of which will be introduced into the reference version soon.  
     1166 
     1167\gmcomment{atmospheric pressure is there!!!!    include its description } 
    12041168 
    12051169% ================================================================ 

  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Chap_LBC.tex

    r6617 r6625  
    11% ================================================================ 
    2 % Chapter — Lateral Boundary Condition (LBC)  
     2% Chapter � Lateral Boundary Condition (LBC)  
    33% ================================================================ 
    44\chapter{Lateral Boundary Condition (LBC) } 
    … …  
    204204%        North fold (\textit{jperio = 3 }to $6)$  
    205205% ------------------------------------------------------------------------------------------------------------- 
    206 \subsection{North-fold (\textit{jperio = 3 }to $6$) } 
     206\subsection{North-fold (\textit{jperio = 3 }to $6)$ } 
    207207\label{LBC_north_fold} 
    208208 
    209209The north fold boundary condition has been introduced in order to handle the north  
    210 boundary of a three-polar ORCA grid. Such a grid has two poles in the northern hemisphere  
    211 (Fig.\ref{Fig_MISC_ORCA_msh}, and thus requires a specific treatment illustrated in Fig.\ref{Fig_North_Fold_T}.  
    212 Further information can be found in \mdl{lbcnfd} module which applies the north fold boundary condition. 
     210boundary of a three-polar ORCA grid. Such a grid has two poles in the northern hemisphere.  
     211\colorbox{yellow}{to be completed...} 
    213212 
    214213%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    251250ocean model. Second order finite difference schemes lead to local discrete  
    252251operators that depend at the very most on one neighbouring point. The only  
    253 non-local computations concern the vertical physics (implicit diffusion,  
     252non-local computations concern the vertical physics (implicit diffusion, 1.5  
    254253turbulent closure scheme, ...) (delocalization over the whole water column),  
    255254and the solving of the elliptic equation associated with the surface pressure  
    256255gradient computation (delocalization over the whole horizontal domain).  
    257256Therefore, a pencil strategy is used for the data sub-structuration  
     257\gmcomment{no idea what this means!} 
    258258: the 3D initial domain is laid out on local processor  
    259259memories following a 2D horizontal topological splitting. Each sub-domain  
    … …  
    264264phase starts: each processor sends to its neighbouring processors the update  
    265265values of the points corresponding to the interior overlapping area to its  
    266 neighbouring sub-domain ($i.e.$ the innermost of the two overlapping rows).  
    267 The communication is done through the Message Passing Interface (MPI).  
    268 The data exchanges between processors are required at the very  
     266neighbouring sub-domain (i.e. the innermost of the two overlapping rows).  
     267The communication is done through message passing. Usually the parallel virtual  
     268language, PVM, is used as it is a standard language available on  nearly  all  
     269MPP computers. More specific languages (i.e. computer dependant languages)  
     270can be easily used to speed up the communication, such as SHEM on a T3E  
     271computer. The data exchanges between processors are required at the very  
    269272place where lateral domain boundary conditions are set in the mono-domain  
    270 computation : the \rou{lbc\_lnk} routine (found in \mdl{lbclnk} module)  
    271 which manages such conditions is interfaced with routines found in \mdl{lib\_mpp} module  
    272 when running on an MPP computer ($i.e.$ when \key{mpp\_mpi} defined).  
    273 It has to be pointed out that when using the MPP version of the model,  
    274 the east-west cyclic boundary condition is done implicitly,  
    275 whilst the south-symmetric boundary condition option is not available. 
     273computation (\S III.10-c): the lbc\_lnk routine which manages such conditions  
     274is substituted by mpplnk.F or mpplnk2.F routine when running on an MPP  
     275computer (\key{mpp\_mpi} defined). It has to be pointed out that when using  
     276the MPP version of the model, the east-west cyclic boundary condition is done  
     277implicitly, whilst the south-symmetric boundary condition option is not available. 
    276278 
    277279%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    283285%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    284286 
    285 In the standard version of \NEMO, the splitting is regular and arithmetic. 
    286 The i-axis is divided by \jp{jpni} and the j-axis by \jp{jpnj} for a number of processors  
    287 \jp{jpnij} most often equal to $jpni \times jpnj$ (parameters set in  
    288  \ngn{nammpp} namelist). Each processor is independent and without message passing  
    289  or synchronous process, programs run alone and access just its own local memory.  
    290  For this reason, the main model dimensions are now the local dimensions of the subdomain (pencil)  
     287In the standard version of the OPA model, the splitting is regular and arithmetic. 
     288 the i-axis is divided by \jp{jpni} and the j-axis by \jp{jpnj} for a number of processors  
     289 \jp{jpnij} most often equal to $jpni \times jpnj$ (model parameters set in  
     290 \mdl{par\_oce}). Each processor is independent and without message passing  
     291 or synchronous process  
     292 \gmcomment{how does a synchronous process relate to this?},  
     293 programs run alone and access just its own local memory. For this reason, the  
     294 main model dimensions are now the local dimensions of the subdomain (pencil)  
    291295 that are named \jp{jpi}, \jp{jpj}, \jp{jpk}. These dimensions include the internal  
    292296 domain and the overlapping rows. The number of rows to exchange (known as  
    … …  
    300304where \jp{jpni}, \jp{jpnj} are the number of processors following the i- and j-axis. 
    301305 
    302 One also defines variables nldi and nlei which correspond to the internal domain bounds,  
    303 and the variables nimpp and njmpp which are the position of the (1,1) grid-point in the global domain.  
    304 An element of $T_{l}$, a local array (subdomain) corresponds to an element of $T_{g}$,  
    305 a global array (whole domain) by the relationship:  
     306\colorbox{yellow}{Figure IV.3: example of a domain splitting with 9 processors and  
     307no east-west cyclic boundary conditions.} 
     308 
     309One also defines variables nldi and nlei which correspond to the internal  
     310domain bounds, and the variables nimpp and njmpp which are the position  
     311of the (1,1) grid-point in the global domain. An element of $T_{l}$, a local array  
     312(subdomain) corresponds to an element of $T_{g}$, a global array  
     313(whole domain) by the relationship:  
    306314\begin{equation} \label{Eq_lbc_nimpp} 
    307315T_{g} (i+nimpp-1,j+njmpp-1,k) = T_{l} (i,j,k), 
    … …  
    312320nproc. In the standard version, a processor has no more than four neighbouring  
    313321processors named nono (for north), noea (east), noso (south) and nowe (west)  
    314 and two variables, nbondi and nbondj, indicate the relative position of the processor : 
     322and two variables, nbondi and nbondj, indicate the relative position of the processor  
     323\colorbox{yellow}{(see Fig.IV.3)}: 
    315324\begin{itemize} 
    316325\item       nbondi = -1    an east neighbour, no west processor, 
    … …  
    323332processor on its overlapping row, and sends the data issued from internal  
    324333domain corresponding to the overlapping row of the other processor. 
     334        
     335\colorbox{yellow}{Figure IV.4: pencil splitting with the additional outer halos } 
    325336 
    326337 
    … …  
    332343global ocean where more than 50 \% of points are land points. For this reason, a  
    333344pre-processing tool can be used to choose the mpp domain decomposition with a  
    334 maximum number of only land points processors, which can then be eliminated (Fig. \ref{Fig_mppini2}) 
    335 (For example, the mpp\_optimiz tools, available from the DRAKKAR web site).  
     345maximum number of only land points processors, which can then be eliminated.  
     346(For example, the mpp\_optimiz tools, available from the DRAKKAR web site.)  
    336347This optimisation is dependent on the specific bathymetry employed. The user  
    337348then chooses optimal parameters \jp{jpni}, \jp{jpnj} and \jp{jpnij} with  
    338349$jpnij < jpni \times jpnj$, leading to the elimination of $jpni \times jpnj - jpnij$  
    339 land processors. When those parameters are specified in \ngn{nammpp} namelist,  
     350land processors. When those parameters are specified in module \mdl{par\_oce},  
    340351the algorithm in the \rou{inimpp2} routine sets each processor's parameters (nbound,  
    341352nono, noea,...) so that the land-only processors are not taken into account.  
    342353 
    343 \gmcomment{Note that the inimpp2 routine is general so that the original inimpp  
     354\colorbox{yellow}{Note that the inimpp2 routine is general so that the original inimpp  
    344355routine should be suppressed from the code.} 
    345356 
    346357When land processors are eliminated, the value corresponding to these locations in  
    347 the model output files is undefined. Note that this is a problem for the meshmask file  
    348 which requires to be defined over the whole domain. Therefore, user should not eliminate  
    349 land processors when creating a meshmask file ($i.e.$ when setting a non-zero value to \np{nn\_msh}). 
     358the model output files is zero. Note that this is a problem for a mesh output file written  
     359by such a model configuration, because model users often divide by the scale factors  
     360($e1t$, $e2t$, etc) and do not expect the grid size to be zero, even on land. It may be  
     361best not to eliminate land processors when running the model especially to write the  
     362mesh files as outputs (when \np{nn\_msh} namelist parameter differs from 0). 
     363%% 
     364\gmcomment{Steven : dont understand this, no land processor means no output file  
     365covering this part of globe; its only when files are stitched together into one that you  
     366can leave a hole} 
     367%% 
    350368 
    351369%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    362380%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    363381 
     382 
     383% ================================================================ 
     384% Open Boundary Conditions  
     385% ================================================================ 
     386\section{Open Boundary Conditions (\key{obc}) (OBC)} 
     387\label{LBC_obc} 
     388%-----------------------------------------nam_obc  ------------------------------------------- 
     389%-    nobc_dta    =    0     !  = 0 the obc data are equal to the initial state 
     390%-                           !  = 1 the obc data are read in 'obc   .dta' files 
     391%-    rn_dpein      =    1.    !  ??? 
     392%-    rn_dpwin      =    1.    !  ??? 
     393%-    rn_dpnin      =   30.    !  ??? 
     394%-    rn_dpsin      =    1.    !  ??? 
     395%-    rn_dpeob      = 1500.    !  time relaxation (days) for the east  open boundary 
     396%-    rn_dpwob      =   15.    !    "        "           for the west  open boundary 
     397%-    rn_dpnob      =  150.    !    "        "           for the north open boundary 
     398%-    rn_dpsob      =   15.    !    "        "           for the south open boundary 
     399%-    ln_obc_clim = .true.   !  climatological obc data files (default T) 
     400%-    ln_vol_cst  = .true.   !  total volume conserved 
     401\namdisplay{namobc}  
     402 
     403It is often necessary to implement a model configuration limited to an oceanic  
     404region or a basin, which communicates with the rest of the global ocean through  
     405''open boundaries''. As stated by \citet{Roed1986}, an open boundary is a  
     406computational border where the aim of the calculations is to allow the perturbations  
     407generated inside the computational domain to leave it without deterioration of the  
     408inner model solution. However, an open boundary also has to let information from  
     409the outer ocean enter the model and should support inflow and outflow conditions.  
     410 
     411The open boundary package OBC is the first open boundary option developed in  
     412NEMO (originally in OPA8.2). It allows the user to  
     413\begin{itemize} 
     414\item tell the model that a boundary is ''open'' and not closed by a wall, for example  
     415by modifying the calculation of the divergence of velocity there; 
     416\item impose values of tracers and velocities at that boundary (values which may  
     417be taken from a climatology): this is the``fixed OBC'' option.  
     418\item calculate boundary values by a sophisticated algorithm combining radiation  
     419and relaxation (``radiative OBC'' option) 
     420\end{itemize} 
     421 
     422Options are defined through the \ngn{namobc} namelist variables. 
     423The package resides in the OBC directory. It is described here in four parts: the  
     424boundary geometry (parameters to be set in \mdl{obc\_par}), the forcing data at  
     425the boundaries (module \mdl{obcdta}),  the radiation algorithm involving the  
     426namelist and module \mdl{obcrad}, and a brief presentation of boundary update  
     427and restart files. 
     428 
     429%---------------------------------------------- 
     430\subsection{Boundary geometry} 
     431\label{OBC_geom} 
     432% 
     433First one has to realize that open boundaries may not necessarily be located  
     434at the extremities of the computational domain. They may exist in the middle  
     435of the domain, for example at Gibraltar Straits if one wants to avoid including  
     436the Mediterranean in an Atlantic domain. This flexibility has been found necessary  
     437for the CLIPPER project \citep{Treguier_al_JGR01}. Because of the complexity of the  
     438geometry of ocean basins, it may even be necessary to have more than one  
     439''west'' open boundary, more than one ''north'', etc. This is not possible with  
     440the OBC option: only one open boundary of each kind, west, east, south and  
     441north is allowed; these names refer to the grid geometry (not to the direction  
     442of the geographical ''west'', ''east'', etc). 
     443 
     444The open boundary geometry is set by a series of parameters in the module  
     445\mdl{obc\_par}. For an eastern open boundary, parameters are \jp{lp\_obc\_east}  
     446(true if an east open boundary exists), \jp{jpieob} the $i$-index along which  
     447the eastern open boundary (eob) is located, \jp{jpjed} the $j$-index at which  
     448it starts, and \jp{jpjef} the $j$-index where it ends (note $d$ is for ''d\'{e}but''  
     449and $f$ for ''fin'' in French). Similar parameters exist for the west, south and  
     450north cases (Table~\ref{Tab_obc_param}). 
     451 
     452 
     453%--------------------------------------------------TABLE-------------------------------------------------- 
     454\begin{table}[htbp]     \begin{center}    \begin{tabular}{|l|c|c|c|} 
     455\hline 
     456Boundary and  & Constant index  & Starting index (d\'{e}but) & Ending index (fin) \\ 
     457Logical flag  &                 &                            &                     \\ 
     458\hline 
     459West          & \jp{jpiwob} $>= 2$         &  \jp{jpjwd}$>= 2$          &  \jp{jpjwf}<= \np{jpjglo}-1 \\ 
     460lp\_obc\_west & $i$-index of a $u$ point   & $j$ of a $T$ point   &$j$ of a $T$ point \\ 
     461\hline 
     462East            & \jp{jpieob}$<=$\np{jpiglo}-2&\jp{jpjed} $>= 2$         & \jp{jpjef}$<=$ \np{jpjglo}-1 \\ 
     463 lp\_obc\_east  & $i$-index of a $u$ point    & $j$ of a $T$ point & $j$ of a $T$ point \\ 
     464\hline 
     465South           & \jp{jpjsob} $>= 2$         & \jp{jpisd} $>= 2$          & \jp{jpisf}$<=$\np{jpiglo}-1 \\ 
     466lp\_obc\_south  & $j$-index of a $v$ point   & $i$ of a $T$ point   & $i$ of a $T$ point \\ 
     467\hline 
     468North           & \jp{jpjnob} $<=$ \np{jpjglo}-2& \jp{jpind} $>= 2$        & \jp{jpinf}$<=$\np{jpiglo}-1 \\ 
     469lp\_obc\_north  & $j$-index of a $v$ point      & $i$  of a $T$ point & $i$ of a $T$ point \\ 
     470\hline 
     471\end{tabular}    \end{center} 
     472\caption{     \label{Tab_obc_param} 
     473Names of different indices relating to the open boundaries. In the case  
     474of a completely open ocean domain with four ocean boundaries, the parameters  
     475take exactly the values indicated.} 
     476\end{table} 
     477%------------------------------------------------------------------------------------------------------------ 
     478 
     479The open boundaries must be along coordinate lines. On the C-grid, the boundary  
     480itself is along a line of normal velocity points: $v$ points for a zonal open boundary  
     481(the south or north one), and $u$ points for a meridional open boundary (the west  
     482or east one). Another constraint is that there still must be a row of masked points  
     483all around the domain, as if the domain were a closed basin (unless periodic conditions  
     484are used together with open boundary conditions). Therefore, an open boundary  
     485cannot be located at the first/last index, namely, 1, \jp{jpiglo} or \jp{jpjglo}. Also,  
     486the open boundary algorithm involves calculating the normal velocity points situated  
     487just on the boundary, as well as the tangential velocity and temperature and salinity  
     488just outside the boundary. This means that for a west/south boundary, normal  
     489velocities and temperature are calculated at the same index \jp{jpiwob} and  
     490\jp{jpjsob}, respectively. For an east/north boundary, the normal velocity is  
     491calculated at index \jp{jpieob} and \jp{jpjnob}, but the ``outside'' temperature is  
     492at index \jp{jpieob}+1 and \jp{jpjnob}+1. This means that \jp{jpieob}, \jp{jpjnob}  
     493cannot be bigger than \jp{jpiglo}-2, \jp{jpjglo}-2.  
     494 
     495 
     496The starting and ending indices are to be thought of as $T$ point indices: in many  
     497cases they indicate the first land $T$-point, at the extremity of an open boundary  
     498(the coast line follows the $f$ grid points, see Fig.~\ref{Fig_obc_north} for an example  
     499of a northern open boundary). All indices are relative to the global domain. In the  
     500free surface case it is possible to have ``ocean corners'', that is, an open boundary  
     501starting and ending in the ocean. 
     502 
     503%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     504\begin{figure}[!t]     \begin{center} 
     505\includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_obc_north.pdf} 
     506\caption{    \label{Fig_obc_north} 
     507Localization of the North open boundary points.} 
     508\end{center}     \end{figure} 
     509%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     510 
     511Although not compulsory, it is highly recommended that the bathymetry in the  
     512vicinity of an open boundary follows the following rule: in the direction perpendicular  
     513to the open line, the water depth should be constant for 4 grid points. This is in  
     514order to ensure that the radiation condition, which involves model variables next  
     515to the boundary, is calculated in a consistent way. On Fig.\ref{Fig_obc_north} we  
     516indicate by an $=$ symbol, the points which should have the same depth. It means  
     517that at the 4 points near the boundary, the bathymetry is cylindrical \gmcomment{not sure  
     518why cylindrical}. The line behind the open $T$-line must be 0 in the bathymetry file  
     519(as shown on Fig.\ref{Fig_obc_north} for example). 
     520 
     521%---------------------------------------------- 
     522\subsection{Boundary data} 
     523\label{OBC_data} 
     524 
     525It is necessary to provide information at the boundaries. The simplest case is  
     526when this information does not change in time and is equal to the initial conditions  
     527(namelist variable \np{nn\_obcdta}=0). This is the case for the standard configuration  
     528EEL5 with open boundaries. When (\np{nn\_obcdta}=1), open boundary information  
     529is read from netcdf files. For convenience the input files are supposed to be similar  
     530to the ''history'' NEMO output files, for dimension names and variable names.  
     531Open boundary arrays must be dimensioned according to the parameters of table~ 
     532\ref{Tab_obc_param}: for example, at the western boundary, arrays have a  
     533dimension of \jp{jpwf}-\jp{jpwd}+1 in the horizontal and \jp{jpk} in the vertical.  
     534 
     535When ocean observations are used to generate the boundary data (a hydrographic  
     536section for example, as in \citet{Treguier_al_JGR01}) it happens often that only the velocity  
     537normal to the boundary is known, which is the reason why the initial OBC code  
     538assumes that only $T$, $S$, and the normal velocity ($u$ or $v$) needs to be  
     539specified. As more and more global model solutions and ocean analysis products  
     540become available, it will be possible to provide information about all the variables  
     541(including the tangential velocity) so that the specification of four variables at each  
     542boundaries will become standard. For the sea surface height, one must distinguish  
     543between the filtered free surface case and the time-splitting or explicit treatment of  
     544the free surface. 
     545 In the first case, it is assumed that the user does not wish to represent high  
     546 frequency motions such as tides. The boundary condition is thus one of zero  
     547 normal gradient of sea surface height at the open boundaries, following \citet{Marchesiello2001}.  
     548No information other than the total velocity needs to be provided at the open  
     549boundaries in that case. In the other two cases (time splitting or explicit free surface),  
     550the user must provide barotropic information (sea surface height and barotropic  
     551velocities) and the use of the Flather algorithm for barotropic variables is  
     552recommanded. However, this algorithm has not yet been fully tested and bugs  
     553remain in NEMO v2.3. Users should read the code carefully before using it. Finally,  
     554in the case of the rigid lid approximation the barotropic streamfunction must be  
     555provided, as documented in \citet{Treguier_al_JGR01}). This option is no longer  
     556recommended but remains in NEMO V2.3. 
     557 
     558One frequently encountered case is when an open boundary domain is constructed  
     559from a global or larger scale NEMO configuration. Assuming the domain corresponds  
     560to indices $ib:ie$, $jb:je$ of the global domain, the bathymetry and forcing of the  
     561small domain can be created by using the following netcdf utility on the global files:  
     562ncks -F $-d\;x,ib,ie$ $-d\;y,jb,je$ (part of the nco series of utilities,  
     563see their \href{http://nco.sourceforge.net}{website}).  
     564The open boundary files can be constructed using ncks  
     565commands, following table~\ref{Tab_obc_ind}.  
     566 
     567%--------------------------------------------------TABLE-------------------------------------------------- 
     568\begin{table}[htbp]     \begin{center}      \begin{tabular}{|l|c|c|c|c|c|} 
     569\hline 
     570OBC  & Variable   & file name      & Index  & Start  & end  \\ 
     571West &  T,S       &   obcwest\_TS.nc &  $ib$+1     &   $jb$+1 &  $je-1$  \\ 
     572     &    U       &   obcwest\_U.nc  &  $ib$+1     &   $jb$+1 &  $je-1$  \\  
     573     &    V       &   obcwest\_V.nc  &  $ib$+1     &   $jb$+1 &  $je-1$  \\        
     574\hline 
     575East &  T,S       &   obceast\_TS.nc &  $ie$-1     &   $jb$+1 &  $je-1$  \\ 
     576     &    U       &   obceast\_U.nc  &  $ie$-2     &   $jb$+1 &  $je-1$  \\  
     577     &    V       &   obceast\_V.nc  &  $ie$-1     &   $jb$+1 &  $je-1$  \\        
     578\hline          
     579South &  T,S      &   obcsouth\_TS.nc &  $jb$+1     &  $ib$+1 &  $ie-1$  \\ 
     580      &    U      &   obcsouth\_U.nc  &  $jb$+1     &  $ib$+1 &  $ie-1$  \\  
     581      &    V      &   obcsouth\_V.nc  &  $jb$+1     &  $ib$+1 &  $ie-1$  \\     
     582\hline 
     583North &  T,S      &   obcnorth\_TS.nc &  $je$-1     &  $ib$+1 &  $ie-1$  \\ 
     584      &    U      &   obcnorth\_U.nc  &  $je$-1     &  $ib$+1 &  $ie-1$  \\  
     585      &    V      &   obcnorth\_V.nc  &  $je$-2     &  $ib$+1 &  $ie-1$  \\   
     586\hline 
     587\end{tabular}     \end{center} 
     588\caption{    \label{Tab_obc_ind} 
     589Requirements for creating open boundary files from a global configuration,  
     590appropriate for the subdomain of indices $ib:ie$, $jb:je$. ``Index'' designates the  
     591$i$ or $j$ index along which the $u$ of $v$ boundary point is situated in the global  
     592configuration, starting and ending with the $j$ or $i$ indices indicated.  
     593For example, to generate file obcnorth\_V.nc, use the command ncks  
     594$-F$ $-d\;y,je-2$  $-d\;x,ib+1,ie-1$ }  
     595\end{table} 
     596%----------------------------------------------------------------------------------------------------------- 
     597 
     598It is assumed that the open boundary files contain the variables for the period of  
     599the model integration. If the boundary files contain one time frame, the boundary  
     600data is held fixed in time. If the files contain 12 values, it is assumed that the input  
     601is a climatology for a repeated annual cycle (corresponding to the case \np{ln\_obc\_clim}  
     602=true). The case of an arbitrary number of time frames is not yet implemented  
     603correctly; the user is required to write his own code in the module \mdl{obc\_dta}  
     604to deal with this situation.  
     605 
     606\subsection{Radiation algorithm} 
     607\label{OBC_rad} 
     608 
     609The art of open boundary management consists in applying a constraint strong  
     610enough that the inner domain "feels" the rest of the ocean, but weak enough 
     611that perturbations are allowed to leave the domain with minimum false reflections  
     612of energy. The constraints are specified separately at each boundary as time  
     613scales for ''inflow'' and ''outflow'' as defined below. The time scales are set (in days)  
     614by namelist parameters such as \np{rn\_dpein}, \np{rn\_dpeob} for the eastern open  
     615boundary for example. When both time scales are zero for a given boundary  
     616($e.g.$ for the western boundary, \jp{lp\_obc\_west}=true, \np{rn\_dpwob}=0 and  
     617\np{rn\_dpwin}=0) this means that the boundary in question is a ''fixed '' boundary  
     618where the solution is set exactly by the boundary data. This is not recommended,  
     619except in combination with increased viscosity in a ''sponge'' layer next to the  
     620boundary in order to avoid spurious reflections.   
     621 
     622 
     623The radiation\/relaxation \gmcomment{the / doesnt seem to appear in the output}  
     624algorithm is applied when either relaxation time (for ''inflow'' or ''outflow'') is  
     625non-zero. It has been developed and tested in the SPEM model and its  
     626successor ROMS \citep{Barnier1996, Marchesiello2001}, which is an  
     627$s$-coordinate model on an Arakawa C-grid. Although the algorithm has  
     628been numerically successful in the CLIPPER Atlantic models, the physics  
     629do not work as expected \citep{Treguier_al_JGR01}. Users are invited to consider  
     630open boundary conditions (OBC hereafter) with some scepticism  
     631\citep{Durran2001, Blayo2005}. 
     632 
     633The first part of the algorithm calculates a phase velocity to determine  
     634whether perturbations tend to propagate toward, or away from, the  
     635boundary. Let us consider a model variable $\phi$.  
     636The phase velocities ($C_{\phi x}$,$C_{\phi y}$) for the variable $\phi$,  
     637in the directions normal and tangential to the boundary are 
     638\begin{equation} \label{Eq_obc_cphi} 
     639C_{\phi x} = \frac{ -\phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{x}  
     640\;\;\;\;\; \;\;\;  
     641C_{\phi y} = \frac{ -\phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{y}.  
     642\end{equation} 
     643Following \citet{Treguier_al_JGR01} and \citet{Marchesiello2001} we retain only  
     644the normal component of the velocity, $C_{\phi x}$, setting $C_{\phi y} =0$  
     645(but unlike the original Orlanski radiation algorithm we retain $\phi_{y}$ in  
     646the expression for $C_{\phi x}$).   
     647 
     648The discrete form of (\ref{Eq_obc_cphi}), described by \citet{Barnier1998}, 
     649takes into account the two rows of grid points situated inside the domain  
     650next to the boundary, and the three previous time steps ($n$, $n-1$, 
     651and $n-2$). The same equation can then be discretized at the boundary at 
     652time steps $n-1$, $n$ and $n+1$ \gmcomment{since the original was three time-level}  
     653in order to extrapolate for the new boundary value $\phi^{n+1}$.  
     654 
     655In the open boundary algorithm as implemented in NEMO v2.3, the new boundary  
     656values are updated differently depending on the sign of $C_{\phi x}$. Let us take  
     657an eastern boundary as an example. The solution for variable $\phi$ at the  
     658boundary is given by a generalized wave equation with phase velocity $C_{\phi}$,  
     659with the addition of a relaxation term, as: 
     660\begin{eqnarray} 
     661\phi_{t} &  =  & -C_{\phi x} \phi_{x} + \frac{1}{\tau_{o}} (\phi_{c}-\phi)  
     662                        \;\;\; \;\;\; \;\;\; (C_{\phi x} > 0), \label{Eq_obc_rado} \\ 
     663\phi_{t} &  =  & \frac{1}{\tau_{i}} (\phi_{c}-\phi)  
     664\;\;\; \;\;\; \;\;\;\;\;\; (C_{\phi x} < 0), \label{Eq_obc_radi} 
     665\end{eqnarray} 
     666where $\phi_{c}$ is the estimate of $\phi$ at the boundary, provided as boundary  
     667data. Note that in (\ref{Eq_obc_rado}), $C_{\phi x}$ is bounded by the ratio  
     668$\delta x/\delta t$ for stability reasons. When $C_{\phi x}$ is eastward (outward  
     669propagation), the radiation condition (\ref{Eq_obc_rado}) is used.  
     670When  $C_{\phi x}$ is westward (inward propagation), (\ref{Eq_obc_radi}) is  
     671used with a strong relaxation to climatology (usually $\tau_{i}=\np{rn\_dpein}=$1~day). 
     672Equation (\ref{Eq_obc_radi}) is solved with a Euler time-stepping scheme. As a  
     673consequence, setting $\tau_{i}$ smaller than, or equal to the time step is equivalent  
     674to a fixed boundary condition. A time scale of one day is usually a good compromise  
     675which guarantees that the inflow conditions remain close to climatology while ensuring  
     676numerical stability.  
     677 
     678In  the case of a western boundary located in the Eastern Atlantic, \citet{Penduff_al_JGR00}  
     679have been able to implement the radiation algorithm without any boundary data,  
     680using persistence from the previous time step instead. This solution has not worked  
     681in other cases \citep{Treguier_al_JGR01}, so that the use of boundary data is recommended.  
     682Even in the outflow condition (\ref{Eq_obc_rado}), we have found it desirable to  
     683maintain a weak relaxation to climatology. The time step is usually chosen so as to  
     684be larger than typical turbulent scales (of order 1000~days \gmcomment{or maybe seconds?}). 
     685 
     686The radiation condition is applied to the model variables: temperature, salinity,  
     687tangential and normal velocities. For normal and tangential velocities, $u$ and $v$,  
     688radiation is applied with phase velocities calculated from $u$ and $v$ respectively.   
     689For the radiation of tracers, we use the phase velocity calculated from the tangential  
     690velocity in order to avoid calculating too many independent radiation velocities and  
     691because tangential velocities and tracers have the same position along the boundary  
     692on a C-grid.   
     693 
     694\subsection{Domain decomposition (\key{mpp\_mpi})} 
     695\label{OBC_mpp} 
     696When \key{mpp\_mpi} is active in the code, the computational domain is divided  
     697into rectangles that are attributed each to a different processor. The open boundary  
     698code is ``mpp-compatible'' up to a certain point. The radiation algorithm will not  
     699work if there is an mpp subdomain boundary parallel to the open boundary at the  
     700index of the boundary, or the grid point after (outside), or three grid points before  
     701(inside). On the other hand, there is no problem if an mpp subdomain boundary  
     702cuts the open boundary perpendicularly. These geometrical limitations must be  
     703checked for by the user (there is no safeguard in the code).   
     704The general principle for the open boundary mpp code is that loops over the open  
     705boundaries {not sure what this means} are performed on local indices (nie0,  
     706nie1, nje0, nje1 for an eastern boundary for instance) that are initialized in module  
     707\mdl{obc\_ini}. Those indices have relevant values on the processors that contain  
     708a segment of an open boundary. For processors that do not include an open  
     709boundary segment, the indices are such that the calculations within the loops are  
     710not performed. 
     711\gmcomment{I dont understand most of the last few sentences} 
     712  
     713Arrays of climatological data that are read from files are seen by all processors  
     714and have the same dimensions for all (for instance, for the eastern boundary,  
     715uedta(jpjglo,jpk,2)). On the other hand, the arrays for the calculation of radiation  
     716are local to each processor (uebnd(jpj,jpk,3,3) for instance).  This allowed the  
     717CLIPPER model for example, to save on memory where the eastern boundary  
     718crossed 8 processors so that \jp{jpj} was much smaller than (\jp{jpjef}-\jp{jpjed}+1).  
     719 
     720\subsection{Volume conservation} 
     721\label{OBC_vol} 
     722 
     723It is necessary to control the volume inside a domain when using open boundaries.  
     724With fixed boundaries, it is enough to ensure that the total inflow/outflow has  
     725reasonable values (either zero or a value compatible with an observed volume  
     726balance). When using radiative boundary conditions it is necessary to have a  
     727volume constraint because each open boundary works independently from the  
     728others. The methodology used to control this volume is identical to the one  
     729coded in the ROMS model \citep{Marchesiello2001}. 
     730 
     731 
     732%---------------------------------------- EXTRAS 
     733\colorbox{yellow}{Explain obc\_vol{\ldots}} 
     734 
     735\colorbox{yellow}{OBC algorithm for update, OBC restart, list of routines where obc key appears{\ldots}} 
     736 
     737\colorbox{yellow}{OBC rigid lid? {\ldots}} 
    364738 
    365739% ==================================================================== 

  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Chap_LDF.tex

    r6617 r6625  
    11 
    22% ================================================================ 
    3 % Chapter ——— Lateral Ocean Physics (LDF) 
     3% Chapter � Lateral Ocean Physics (LDF) 
    44% ================================================================ 
    55\chapter{Lateral Ocean Physics (LDF)} 
    … …  
    6868When none of the \textbf{key\_dynldf\_...} and \textbf{key\_traldf\_...} keys are  
    6969defined, a constant value is used over the whole ocean for momentum and  
    70 tracers, which is specified through the \np{rn\_ahm\_0\_lap} and \np{rn\_aht\_0} namelist  
     70tracers, which is specified through the \np{rn\_ahm0} and \np{rn\_aht0} namelist  
    7171parameters. 
    7272 
    … …  
    7777mixing coefficients will require 3D arrays. In the 1D option, a hyperbolic variation  
    7878of the lateral mixing coefficient is introduced in which the surface value is  
    79 \np{rn\_aht\_0} (\np{rn\_ahm\_0\_lap}), the bottom value is 1/4 of the surface value,  
     79\np{rn\_aht0} (\np{rn\_ahm0}), the bottom value is 1/4 of the surface value,  
    8080and the transition takes place around z=300~m with a width of 300~m  
    8181($i.e.$ both the depth and the width of the inflection point are set to 300~m).  
    … …  
    9393\end{equation} 
    9494where $e_{max}$ is the maximum of $e_1$ and $e_2$ taken over the whole masked  
    95 ocean domain, and $A_o^l$ is the \np{rn\_ahm\_0\_lap} (momentum) or \np{rn\_aht\_0} (tracer)  
     95ocean domain, and $A_o^l$ is the \np{rn\_ahm0} (momentum) or \np{rn\_aht0} (tracer)  
    9696namelist parameter. This variation is intended to reflect the lesser need for subgrid  
    9797scale eddy mixing where the grid size is smaller in the domain. It was introduced in  
    … …  
    105105Other formulations can be introduced by the user for a given configuration.  
    106106For example, in the ORCA2 global ocean model (see Configurations), the laplacian  
    107 viscosity operator uses \np{rn\_ahm\_0\_lap}~= 4.10$^4$ m$^2$/s poleward of 20$^{\circ}$  
    108 north and south and decreases linearly to \np{rn\_aht\_0}~= 2.10$^3$ m$^2$/s  
     107viscosity operator uses \np{rn\_ahm0}~= 4.10$^4$ m$^2$/s poleward of 20$^{\circ}$  
     108north and south and decreases linearly to \np{rn\_aht0}~= 2.10$^3$ m$^2$/s  
    109109at the equator \citep{Madec_al_JPO96, Delecluse_Madec_Bk00}. This modification  
    110110can be found in routine \rou{ldf\_dyn\_c2d\_orca} defined in \mdl{ldfdyn\_c2d}.  
    … …  
    120120\subsubsection{Space and Time Varying Mixing Coefficients} 
    121121 
    122 There are no default specifications of space and time varying mixing coefficient.  One 
    123 available case is specific to the ORCA2 and ORCA05 global ocean configurations. It 
    124 provides only a tracer mixing coefficient for eddy induced velocity (ORCA2) or both 
    125 iso-neutral and eddy induced velocity (ORCA05) that depends on the local growth rate of 
    126 baroclinic instability. This specification is actually used when an ORCA key 
     122There is no default specification of space and time varying mixing coefficient.  
     123The only case available is specific to the ORCA2 and ORCA05 global ocean  
     124configurations. It provides only a tracer  
     125mixing coefficient for eddy induced velocity (ORCA2) or both iso-neutral and  
     126eddy induced velocity (ORCA05) that depends on the local growth rate of  
     127baroclinic instability. This specification is actually used when an ORCA key  
    127128and both \key{traldf\_eiv} and \key{traldf\_c2d} are defined. 
    128  
    129 \subsubsection{Smagorinsky viscosity (\key{dynldf\_c3d} and \key{dynldf\_smag})} 
    130  
    131 The \key{dynldf\_smag} key activates a 3D, time-varying viscosity that depends on the 
    132 resolved motions. Following \citep{Smagorinsky_93} the viscosity coefficient is set 
    133 proportional to a local deformation rate based on the horizontal shear and tension, 
    134 namely: 
    135  
    136 \begin{equation} 
    137 A_{m_{Smag}} = \left(\frac{{\sf CM_{Smag}}}{\pi}\right)^2L^2\vert{D}\vert 
    138 \end{equation} 
    139  
    140 \noindent where the deformation rate $\vert{D}\vert$ is given by  
    141  
    142 \begin{equation} 
    143 \vert{D}\vert=\sqrt{\left({\frac{\partial{u}} {\partial{x}}} 
    144                          -{\frac{\partial{v}} {\partial{y}}}\right)^2 
    145                  +  \left({\frac{\partial{u}} {\partial{y}}} 
    146                          +{\frac{\partial{v}} {\partial{x}}}\right)^2}  
    147 \end{equation} 
    148  
    149 \noindent and $L$ is the local gridscale given by: 
    150  
    151 \begin{equation} 
    152 L^2 = \frac{2{e_1}^2 {e_2}^2}{\left ( {e_1}^2 + {e_2}^2 \right )} 
    153 \end{equation} 
    154  
    155 \citep{Griffies_Hallberg_MWR00} suggest values in the range 2.2 to 4.0 of the coefficient 
    156 $\sf CM_{Smag}$ for oceanic flows. This value is set via the \np{rn\_cmsmag\_1} namelist 
    157 parameter. An additional parameter: \np{rn\_cmsh} is included in NEMO for experimenting 
    158 with the contribution of the shear term. A value of 1.0 (the default) calculates the 
    159 deformation rate as above; a value of 0.0 will discard the shear term entirely. 
    160  
    161 For numerical stability, the calculated viscosity is bounded according to the following: 
    162  
    163 \begin{equation} 
    164 {\rm MIN}\left ({ L^2\over {8\Delta{t}}}, rn\_ahm\_m\_lap\right ) \geq A_{m_{Smag}}  
    165                                                                   \geq rn\_ahm\_0\_lap 
    166 \end{equation} 
    167  
    168 \noindent with both parameters for the upper and lower bounds being provided via the 
    169 indicated namelist parameters. 
    170  
    171 \bigskip When $ln\_dynldf\_bilap = .true.$, a biharmonic version of the Smagorinsky 
    172 viscosity is also available which sets a coefficient for the biharmonic viscosity as: 
    173  
    174 \begin{equation} 
    175 B_{m_{Smag}} = - \left(\frac{{\sf CM_{bSmag}}}{\pi}\right)^2 {L^4\over 8}\vert{D}\vert 
    176 \end{equation} 
    177  
    178 \noindent which is bounded according to: 
    179  
    180 \begin{equation} 
    181 {\rm MAX}\left (-{ L^4\over {64\Delta{t}}}, rn\_ahm\_m\_blp\right ) \leq B_{m_{Smag}}  
    182                                                                     \leq rn\_ahm\_0\_blp 
    183 \end{equation} 
    184  
    185 \noindent Note the reversal of the inequalities here because NEMO requires the biharmonic 
    186 coefficients as negative numbers. $\sf CM_{bSmag}$ is set via the \np{rn\_cmsmag\_2} 
    187 namelist parameter and the bounding values have corresponding entries in the namelist too. 
    188  
    189 \bigskip The current implementation in NEMO also allows for 3D, time-varying diffusivities 
    190 to be set using the Smagorinsky approach. Users should note that this option is not 
    191 recommended for many applications since diffusivities will tend to be largest near 
    192 boundaries (where shears are greatest) leading to spurious upwellings 
    193 (\citep{Griffies_Bk04}, chapter 18.3.4). Nevertheless the option is there for those 
    194 wishing to experiment. This choice requires both \key{traldf\_c3d} and \key{traldf\_smag} 
    195 and uses the \np{rn\_chsmag} (${\sf CH_{Smag}}$), \np{rn\_smsh} and \np{rn\_aht\_m} 
    196 namelist parameters in an analogous way to \np{rn\_cmsmag\_1}, \np{rn\_cmsh} and 
    197 \np{rn\_ahm\_m\_lap} (see above) to set the diffusion coefficient: 
    198  
    199 \begin{equation} 
    200 A_{h_{Smag}} = \left(\frac{{\sf CH_{Smag}}}{\pi}\right)^2L^2\vert{D}\vert 
    201 \end{equation} 
    202  
    203   
    204 For numerical stability, the calculated diffusivity is bounded according to the following: 
    205  
    206 \begin{equation} 
    207 {\rm MIN}\left ({ L^2\over {8\Delta{t}}}, rn\_aht\_m\right ) \geq A_{h_{Smag}}  
    208                                                              \geq rn\_aht\_0 
    209 \end{equation} 
    210  
    211  
    212129 
    213130$\ $\newline    % force a new ligne 
    … …  
    227144(3) for isopycnal diffusion on momentum or tracers, an additional purely  
    228145horizontal background diffusion with uniform coefficient can be added by  
    229 setting a non zero value of \np{rn\_ahmb\_0} or \np{rn\_ahtb\_0}, a background horizontal  
     146setting a non zero value of \np{rn\_ahmb0} or \np{rn\_ahtb0}, a background horizontal  
    230147eddy viscosity or diffusivity coefficient (namelist parameters whose default  
    231148values are $0$). However, the technique used to compute the isopycnal  

  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Chap_MISC.tex

    r6617 r6625  
    3434has been made to set them in a generic way. However, examples of how  
    3535they can be set up is given in the ORCA 2\deg and 0.5\deg configurations. For example,  
    36 for details of implementation in ORCA2, search:  
    37 \texttt{ IF( cp\_cfg == "orca" .AND. jp\_cfg == 2 ) } 
     36for details of implementation in ORCA2, search: 
     37\vspace{-10pt}   
     38\begin{alltt}   
     39\tiny     
     40\begin{verbatim} 
     41IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) 
     42\end{verbatim}   
     43\end{alltt} 
    3844 
    3945% ------------------------------------------------------------------------------------------------------------- 
    … …  
    8389%-------------------------------------------------------------------------------------------------------------- 
    8490 
     91\colorbox{yellow}{Add a short description of CLA staff here or in lateral boundary condition chapter?} 
    8592Options are defined through the  \ngn{namcla} namelist variables. 
    86 This option is an obsolescent feature that will be removed in version 3.7 and followings.  
    8793 
    8894%The problem is resolved here by allowing the mixing of tracers and mass/volume between non-adjacent water columns at nominated regions within the model. Momentum is not mixed. The scheme conserves total tracer content, and total volume (the latter in $z*$- or $s*$-coordinate), and maintains compatibility between the tracer and mass/volume budgets.   

  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Chap_Model_Basics.tex

    r6617 r6625  
    247247sufficient to solve a linearized version of (\ref{Eq_PE_ssh}), which still allows  
    248248to take into account freshwater fluxes applied at the ocean surface \citep{Roullet_Madec_JGR00}. 
    249 Nevertheless, with the linearization, an exact conservation of heat and salt contents is lost. 
    250249 
    251250The filtering of EGWs in models with a free surface is usually a matter of discretisation  
    252 of the temporal derivatives, using a split-explicit method \citep{Killworth_al_JPO91, Zhang_Endoh_JGR92}  
    253 or the implicit scheme \citep{Dukowicz1994} or the addition of a filtering force in the momentum equation  
    254 \citep{Roullet_Madec_JGR00}. With the present release, \NEMO offers the choice between  
    255 an explicit free surface (see \S\ref{DYN_spg_exp}) or a split-explicit scheme strongly  
    256 inspired the one proposed by \citet{Shchepetkin_McWilliams_OM05} (see \S\ref{DYN_spg_ts}). 
    257  
    258 %\newpage 
    259 %$\ $\newline    % force a new line 
     251of the temporal derivatives, using the time splitting method \citep{Killworth_al_JPO91, Zhang_Endoh_JGR92}  
     252or the implicit scheme \citep{Dukowicz1994}. In \NEMO, we use a slightly different approach  
     253developed by \citet{Roullet_Madec_JGR00}: the damping of EGWs is ensured by introducing an  
     254additional force in the momentum equation. \eqref{Eq_PE_dyn} becomes:  
     255\begin{equation} \label{Eq_PE_flt} 
     256\frac{\partial {\rm {\bf U}}_h }{\partial t}= {\rm {\bf M}} 
     257- g \nabla \left( \tilde{\rho} \ \eta \right)  
     258- g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right)  
     259\end{equation} 
     260where $T_c$, is a parameter with dimensions of time which characterizes the force,  
     261$\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, and $\rm {\bf M}$  
     262represents the collected contributions of the Coriolis, hydrostatic pressure gradient,  
     263non-linear and viscous terms in \eqref{Eq_PE_dyn}. 
     264 
     265The new force can be interpreted as a diffusion of vertically integrated volume flux divergence.  
     266The time evolution of $D$ is thus governed by a balance of two terms, $-g$ \textbf{A} $\eta$  
     267and $g \, T_c \,$ \textbf{A} $D$, associated with a propagative regime and a diffusive regime  
     268in the temporal spectrum, respectively. In the diffusive regime, the EGWs no longer propagate,  
     269$i.e.$ they are stationary and damped. The diffusion regime applies to the modes shorter than  
     270$T_c$. For longer ones, the diffusion term vanishes. Hence, the temporally unresolved EGWs  
     271can be damped by choosing $T_c > \rdt$. \citet{Roullet_Madec_JGR00} demonstrate that  
     272(\ref{Eq_PE_flt}) can be integrated with a leap frog scheme except the additional term which  
     273has to be computed implicitly. This is not surprising since the use of a large time step has a  
     274necessarily numerical cost. Two gains arise in comparison with the previous formulations.  
     275Firstly, the damping of EGWs can be quantified through the magnitude of the additional term.  
     276Secondly, the numerical scheme does not need any tuning. Numerical stability is ensured as  
     277soon as $T_c > \rdt$. 
     278 
     279When the variations of free surface elevation are small compared to the thickness of the first  
     280model layer, the free surface equation (\ref{Eq_PE_ssh}) can be linearized. As emphasized  
     281by \citet{Roullet_Madec_JGR00} the linearization of (\ref{Eq_PE_ssh}) has consequences on the  
     282conservation of salt in the model. With the nonlinear free surface equation, the time evolution  
     283of the total salt content is  
     284\begin{equation} \label{Eq_PE_salt_content} 
     285    \frac{\partial }{\partial t}\int\limits_{D\eta } {S\;dv}  
     286                        =\int\limits_S {S\;(-\frac{\partial \eta }{\partial t}-D+P-E)\;ds} 
     287\end{equation} 
     288where $S$ is the salinity, and the total salt is integrated over the whole ocean volume  
     289$D_\eta$ bounded by the time-dependent free surface. The right hand side (which is an  
     290integral over the free surface) vanishes when the nonlinear equation (\ref{Eq_PE_ssh})  
     291is satisfied, so that the salt is perfectly conserved. When the free surface equation is  
     292linearized, \citet{Roullet_Madec_JGR00} show that the total salt content integrated in the fixed  
     293volume $D$ (bounded by the surface $z=0$) is no longer conserved: 
     294\begin{equation} \label{Eq_PE_salt_content_linear} 
     295         \frac{\partial }{\partial t}\int\limits_D {S\;dv}  
     296               = - \int\limits_S {S\;\frac{\partial \eta }{\partial t}ds}  
     297\end{equation} 
     298 
     299The right hand side of (\ref{Eq_PE_salt_content_linear}) is small in equilibrium solutions  
     300\citep{Roullet_Madec_JGR00}. It can be significant when the freshwater forcing is not balanced and  
     301the globally averaged free surface is drifting. An increase in sea surface height \textit{$\eta $}  
     302results in a decrease of the salinity in the fixed volume $D$. Even in that case though,  
     303the total salt integrated in the variable volume $D_{\eta}$ varies much less, since  
     304(\ref{Eq_PE_salt_content_linear}) can be rewritten as  
     305\begin{equation} \label{Eq_PE_salt_content_corrected} 
     306\frac{\partial }{\partial t}\int\limits_{D\eta } {S\;dv}  
     307=\frac{\partial}{\partial t} \left[ \;{\int\limits_D {S\;dv} +\int\limits_S {S\eta \;ds} } \right] 
     308=\int\limits_S {\eta \;\frac{\partial S}{\partial t}ds} 
     309\end{equation} 
     310 
     311Although the total salt content is not exactly conserved with the linearized free surface,  
     312its variations are driven by correlations of the time variation of surface salinity with the  
     313sea surface height, which is a negligible term. This situation contrasts with the case of  
     314the rigid lid approximation in which case freshwater forcing is represented by a virtual  
     315salt flux, leading to a spurious source of salt at the ocean surface  
     316\citep{Huang_JPO93, Roullet_Madec_JGR00}. 
     317 
     318\newpage 
     319$\ $\newline    % force a new ligne 
    260320 
    261321% ================================================================ 
    … …  
    713773\end{equation} 
    714774 
    715 The equations solved by the ocean model \eqref{Eq_PE} in $s-$coordinate can be written as follows (see Appendix~\ref{Apdx_A_momentum}): 
     775The equations solved by the ocean model \eqref{Eq_PE} in $s-$coordinate can be written as follows: 
    716776 
    717777 \vspace{0.5cm} 
    718 $\bullet$ Vector invariant form of the momentum equation : 
     778* momentum equation: 
    719779\begin{multline} \label{Eq_PE_sco_u} 
    720 \frac{\partial  u  }{\partial t}= 
     780\frac{1}{e_3} \frac{\partial \left(  e_3\,u  \right) }{\partial t}= 
    721781   +   \left( {\zeta +f} \right)\,v                                     
    722782   -   \frac{1}{2\,e_1} \frac{\partial}{\partial i} \left(  u^2+v^2   \right)  
    … …  
    727787\end{multline} 
    728788\begin{multline} \label{Eq_PE_sco_v} 
    729 \frac{\partial v }{\partial t}= 
     789\frac{1}{e_3} \frac{\partial \left(  e_3\,v  \right) }{\partial t}= 
    730790   -   \left( {\zeta +f} \right)\,u    
    731791   -   \frac{1}{2\,e_2 }\frac{\partial }{\partial j}\left(  u^2+v^2  \right)         
    … …  
    735795   +  D_v^{\vect{U}}  +   F_v^{\vect{U}} \quad 
    736796\end{multline} 
    737  
    738  \vspace{0.5cm} 
    739 $\bullet$ Vector invariant form of the momentum equation : 
    740 \begin{multline} \label{Eq_PE_sco_u} 
    741 \frac{1}{e_3} \frac{\partial \left(  e_3\,u  \right) }{\partial t}= 
    742    +   \left( { f + \frac{1}{e_1 \; e_2 } 
    743                \left(    v \frac{\partial e_2}{\partial i} 
    744                   -u \frac{\partial e_1}{\partial j}  \right)}    \right) \, v    \\ 
    745    - \frac{1}{e_1 \; e_2 \; e_3 }   \left(  
    746                \frac{\partial \left( {e_2 \, e_3 \, u\,u} \right)}{\partial i} 
    747       +        \frac{\partial \left( {e_1 \, e_3 \, v\,u} \right)}{\partial j}   \right) 
    748    - \frac{1}{e_3 }\frac{\partial \left( { \omega\,u} \right)}{\partial k}    \\ 
    749    - \frac{1}{e_1} \frac{\partial}{\partial i} \left( \frac{p_s + p_h}{\rho _o}    \right)     
    750    +  g\frac{\rho }{\rho _o}\sigma _1  
    751    +   D_u^{\vect{U}}  +   F_u^{\vect{U}} \quad 
    752 \end{multline} 
    753 \begin{multline} \label{Eq_PE_sco_v} 
    754 \frac{1}{e_3} \frac{\partial \left(  e_3\,v  \right) }{\partial t}= 
    755    -   \left( { f + \frac{1}{e_1 \; e_2} 
    756                \left(    v \frac{\partial e_2}{\partial i} 
    757                   -u \frac{\partial e_1}{\partial j}  \right)}    \right) \, u   \\ 
    758    - \frac{1}{e_1 \; e_2 \; e_3 }   \left(  
    759                \frac{\partial \left( {e_2 \; e_3  \,u\,v} \right)}{\partial i} 
    760       +        \frac{\partial \left( {e_1 \; e_3  \,v\,v} \right)}{\partial j}   \right) 
    761                  - \frac{1}{e_3 } \frac{\partial \left( { \omega\,v} \right)}{\partial k}    \\ 
    762    -   \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho _o}  \right)  
    763     +  g\frac{\rho }{\rho _o }\sigma _2    
    764    +  D_v^{\vect{U}}  +   F_v^{\vect{U}} \quad 
    765 \end{multline} 
    766  
    767797where the relative vorticity, \textit{$\zeta $}, the surface pressure gradient, and the hydrostatic  
    768798pressure have the same expressions as in $z$-coordinates although they do not represent  
    769799exactly the same quantities. $\omega$ is provided by the continuity equation  
    770800(see Appendix~\ref{Apdx_A}): 
     801 
    771802\begin{equation} \label{Eq_PE_sco_continuity} 
    772803\frac{\partial e_3}{\partial t} + e_3 \; \chi + \frac{\partial \omega }{\partial s} = 0    
    … …  
    778809 
    779810 \vspace{0.5cm} 
    780 $\bullet$ tracer equations: 
     811* tracer equations: 
    781812\begin{multline} \label{Eq_PE_sco_t} 
    782813\frac{1}{e_3} \frac{\partial \left(  e_3\,T  \right) }{\partial t}= 
    … …  
    9921023\label{PE_zco_tilde} 
    9931024 
    994 The $\tilde{z}$-coordinate has been developed by \citet{Leclair_Madec_OM11}. 
    995 It is available in \NEMO since the version 3.4. Nevertheless, it is currently not robust enough  
    996 to be used in all possible configurations. Its use is therefore not recommended. 
    997  
     1025The $\tilde{z}$-coordinate has been developed by \citet{Leclair_Madec_OM10s}. 
     1026It is not available in the current version of \NEMO. 
    9981027 
    9991028\newpage  
    … …  
    11281157operator acting along $s-$surfaces (see \S\ref{LDF}). 
    11291158 
    1130 \subsubsection{Lateral Laplacian tracer diffusive operator} 
    1131  
    1132 The lateral Laplacian tracer diffusive operator is defined by (see Appendix~\ref{Apdx_B}): 
     1159\subsubsection{Lateral second order tracer diffusive operator} 
     1160 
     1161The lateral second order tracer diffusive operator is defined by (see Appendix~\ref{Apdx_B}): 
    11331162\begin{equation} \label{Eq_PE_iso_tensor} 
    11341163D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad  
    … …  
    11511180ocean (see Appendix~\ref{Apdx_B}). 
    11521181 
    1153 For \textit{iso-level} diffusion, $r_1$ and $r_2 $ are zero. $\Re $ reduces to the identity  
    1154 in the horizontal direction, no rotation is applied.  
    1155  
    11561182For \textit{geopotential} diffusion, $r_1$ and $r_2 $ are the slopes between the  
    1157 geopotential and computational surfaces: they are equal to $\sigma _1$ and $\sigma _2$,  
    1158 respectively (see \eqref{Eq_PE_sco_slope} ). 
     1183geopotential and computational surfaces: in $z$-coordinates they are zero  
     1184($r_1 = r_2 = 0$) while in $s$-coordinate (including $\textit{z*}$ case) they are  
     1185equal to $\sigma _1$ and $\sigma _2$, respectively (see \eqref{Eq_PE_sco_slope} ). 
    11591186 
    11601187For \textit{isoneutral} diffusion $r_1$ and $r_2$ are the slopes between the isoneutral  
    … …  
    12041231to zero in the vicinity of the boundaries. The latter strategy is used in \NEMO (cf. Chap.~\ref{LDF}). 
    12051232 
    1206 \subsubsection{Lateral bilaplacian tracer diffusive operator} 
    1207  
    1208 The lateral bilaplacian tracer diffusive operator is defined by: 
     1233\subsubsection{Lateral fourth order tracer diffusive operator} 
     1234 
     1235The lateral fourth order tracer diffusive operator is defined by: 
    12091236\begin{equation} \label{Eq_PE_bilapT} 
    12101237D^{lT}=\Delta \left( {A^{lT}\;\Delta T} \right)  
    12111238\qquad \text{where} \  D^{lT}=\Delta \left( {A^{lT}\;\Delta T} \right) 
    12121239 \end{equation} 
     1240 
    12131241It is the second order operator given by \eqref{Eq_PE_iso_tensor} applied twice with  
    12141242the eddy diffusion coefficient correctly placed.  
    12151243 
    1216 \subsubsection{Lateral Laplacian momentum diffusive operator} 
    1217  
    1218 The Laplacian momentum diffusive operator along $z$- or $s$-surfaces is found by  
     1244 
     1245\subsubsection{Lateral second order momentum diffusive operator} 
     1246 
     1247The second order momentum diffusive operator along $z$- or $s$-surfaces is found by  
    12191248applying \eqref{Eq_PE_lap_vector} to the horizontal velocity vector (see Appendix~\ref{Apdx_B}): 
    12201249\begin{equation} \label{Eq_PE_lapU} 
    … …  
    12501279of the Equator in a geographical coordinate system \citep{Lengaigne_al_JGR03}. 
    12511280 
    1252 \subsubsection{lateral bilaplacian momentum diffusive operator} 
     1281\subsubsection{lateral fourth order momentum diffusive operator} 
    12531282 
    12541283As for tracers, the fourth order momentum diffusive operator along $z$ or $s$-surfaces  

  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Chap_Model_Basics_zstar.tex

    r6617 r6625  
    11% ================================================================ 
    2 % Chapter 1 ——— Model Basics 
     2% Chapter 1 � Model Basics 
    33% ================================================================ 
    44% ================================================================ 

  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Chap_SBC.tex

    r6617 r6625  
    11% ================================================================ 
    2 % Chapter —— Surface Boundary Condition (SBC, ISF, ICB)  
     2% Chapter � Surface Boundary Condition (SBC, ISF, ICB)  
    33% ================================================================ 
    44\chapter{Surface Boundary Condition (SBC, ISF, ICB) } 
    … …  
    1717   \item the two components of the surface ocean stress $\left( {\tau _u \;,\;\tau _v} \right)$ 
    1818   \item the incoming solar and non solar heat fluxes $\left( {Q_{ns} \;,\;Q_{sr} } \right)$ 
    19    \item the surface freshwater budget $\left( {\textit{emp}} \right)$ 
    20    \item the surface salt flux associated with freezing/melting of seawater $\left( {\textit{sfx}} \right)$ 
     19   \item the surface freshwater budget $\left( {\textit{emp},\;\textit{emp}_S } \right)$ 
    2120\end{itemize} 
    2221plus an optional field: 
    … …  
    2827are controlled by namelist \ngn{namsbc} variables: an analytical formulation (\np{ln\_ana}~=~true),  
    2928a flux formulation (\np{ln\_flx}~=~true), a bulk formulae formulation (CORE  
    30 (\np{ln\_blk\_core}~=~true), CLIO (\np{ln\_blk\_clio}~=~true) or MFS 
     29(\np{ln\_core}~=~true), CLIO (\np{ln\_clio}~=~true) or MFS 
    3130\footnote { Note that MFS bulk formulae compute fluxes only for the ocean component} 
    32 (\np{ln\_blk\_mfs}~=~true) bulk formulae) and a coupled or mixed forced/coupled formulation  
    33 (exchanges with a atmospheric model via the OASIS coupler) (\np{ln\_cpl} or \np{ln\_mixcpl}~=~true).  
    34 When used ($i.e.$ \np{ln\_apr\_dyn}~=~true), the atmospheric pressure forces both ocean and ice dynamics. 
    35  
    36 The frequency at which the forcing fields have to be updated is given by the \np{nn\_fsbc} namelist parameter.  
     31(\np{ln\_mfs}~=~true) bulk formulae) and a coupled  
     32formulation (exchanges with a atmospheric model via the OASIS coupler)  
     33(\np{ln\_cpl}~=~true). When used, the atmospheric pressure forces both  
     34ocean and ice dynamics (\np{ln\_apr\_dyn}~=~true). 
     35The frequency at which the six or seven fields have to be updated is the \np{nn\_fsbc}  
     36namelist parameter.  
    3737When the fields are supplied from data files (flux and bulk formulations), the input fields  
    38 need not be supplied on the model grid. Instead a file of coordinates and weights can  
     38need not be supplied on the model grid.  Instead a file of coordinates and weights can  
    3939be supplied which maps the data from the supplied grid to the model points  
    4040(so called "Interpolation on the Fly", see \S\ref{SBC_iof}). 
    … …  
    4242can be masked to avoid spurious results in proximity of the coasts  as large sea-land gradients characterize 
    4343most of the atmospheric variables. 
    44  
    4544In addition, the resulting fields can be further modified using several namelist options.  
    46 These options control  
    47 \begin{itemize} 
    48 \item the rotation of vector components supplied relative to an east-north  
    49 coordinate system onto the local grid directions in the model ;  
    50 \item the addition of a surface restoring term to observed SST and/or SSS (\np{ln\_ssr}~=~true) ;  
    51 \item the modification of fluxes below ice-covered areas (using observed ice-cover or a sea-ice model) (\np{nn\_ice}~=~0,1, 2 or 3) ;  
    52 \item the addition of river runoffs as surface freshwater fluxes or lateral inflow (\np{ln\_rnf}~=~true) ;  
    53 \item the addition of isf melting as lateral inflow (parameterisation) (\np{nn\_isf}~=~2 or 3 and \np{ln\_isfcav}~=~false)  
    54 or as fluxes applied at the land-ice ocean interface (\np{nn\_isf}~=~1 or 4 and \np{ln\_isfcav}~=~true) ;  
    55 \item the addition of a freshwater flux adjustment in order to avoid a mean sea-level drift (\np{nn\_fwb}~=~0,~1~or~2) ;  
    56 \item the transformation of the solar radiation (if provided as daily mean) into a diurnal cycle (\np{ln\_dm2dc}~=~true) ;  
    57 and a neutral drag coefficient can be read from an external wave model (\np{ln\_cdgw}~=~true).  
    58 \end{itemize} 
    59 The latter option is possible only in case core or mfs bulk formulas are selected. 
     45These options control  the rotation of vector components supplied relative to an east-north  
     46coordinate system onto the local grid directions in the model; the addition of a surface  
     47restoring term to observed SST and/or SSS (\np{ln\_ssr}~=~true); the modification of fluxes  
     48below ice-covered areas (using observed ice-cover or a sea-ice model)  
     49(\np{nn\_ice}~=~0,1, 2 or 3); the addition of river runoffs as surface freshwater  
     50fluxes or lateral inflow (\np{ln\_rnf}~=~true); the addition of isf melting as lateral inflow (parameterisation)  
     51(\np{nn\_isf}~=~2 or 3 and \np{ln\_isfcav}~=~false) or as surface flux at the land-ice ocean interface 
     52(\np{nn\_isf}~=~1 or 4 and \np{ln\_isfcav}~=~true);  
     53the addition of a freshwater flux adjustment in order to avoid a mean sea-level drift (\np{nn\_fwb}~=~0,~1~or~2); the  
     54transformation of the solar radiation (if provided as daily mean) into a diurnal  
     55cycle (\np{ln\_dm2dc}~=~true); and a neutral drag coefficient can be read from an external wave  
     56model (\np{ln\_cdgw}~=~true). The latter option is possible only in case core or mfs bulk formulas are selected. 
    6057 
    6158In this chapter, we first discuss where the surface boundary condition appears in the 
    … …  
    7673 
    7774The surface ocean stress is the stress exerted by the wind and the sea-ice  
    78 on the ocean. It is applied in \mdl{dynzdf} module as a surface boundary condition of the  
    79 computation of the momentum vertical mixing trend (see \eqref{Eq_dynzdf_sbc} in \S\ref{DYN_zdf}). 
    80 As such, it has to be provided as a 2D vector interpolated  
    81 onto the horizontal velocity ocean mesh, $i.e.$ resolved onto the model  
    82 (\textbf{i},\textbf{j}) direction at $u$- and $v$-points. 
     75on the ocean. The two components of stress are assumed to be interpolated  
     76onto the ocean mesh, $i.e.$ resolved onto the model (\textbf{i},\textbf{j}) direction  
     77at $u$- and $v$-points They are applied as a surface boundary condition of the  
     78computation of the momentum vertical mixing trend (\mdl{dynzdf} module) : 
     79\begin{equation} \label{Eq_sbc_dynzdf} 
     80\left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1} 
     81    = \frac{1}{\rho _o} \binom{\tau _u}{\tau _v } 
     82\end{equation} 
     83where $(\tau _u ,\;\tau _v )=(utau,vtau)$ are the two components of the wind  
     84stress vector in the $(\textbf{i},\textbf{j})$ coordinate system. 
    8385 
    8486The surface heat flux is decomposed into two parts, a non solar and a solar heat  
    8587flux, $Q_{ns}$ and $Q_{sr}$, respectively. The former is the non penetrative part  
    86 of the heat flux ($i.e.$ the sum of sensible, latent and long wave heat fluxes  
    87 plus the heat content of the mass exchange with the atmosphere and sea-ice).  
    88 It is applied in \mdl{trasbc} module as a surface boundary condition trend of  
    89 the first level temperature time evolution equation (see \eqref{Eq_tra_sbc}  
    90 and \eqref{Eq_tra_sbc_lin} in \S\ref{TRA_sbc}).  
    91 The latter is the penetrative part of the heat flux. It is applied as a 3D  
    92 trends of the temperature equation (\mdl{traqsr} module) when \np{ln\_traqsr}=\textit{true}. 
    93 The way the light penetrates inside the water column is generally a sum of decreasing  
    94 exponentials (see \S\ref{TRA_qsr}).  
    95  
    96 The surface freshwater budget is provided by the \textit{emp} field. 
    97 It represents the mass flux exchanged with the atmosphere (evaporation minus precipitation)  
    98 and possibly with the sea-ice and ice shelves (freezing minus melting of ice).  
    99 It affects both the ocean in two different ways:  
    100 $(i)$   it changes the volume of the ocean and therefore appears in the sea surface height  
    101 equation as a volume flux, and  
    102 $(ii)$  it changes the surface temperature and salinity through the heat and salt contents  
    103 of the mass exchanged with the atmosphere, the sea-ice and the ice shelves.  
    104  
     88of the heat flux ($i.e.$ the sum of sensible, latent and long wave heat fluxes).  
     89It is applied as a surface boundary condition trend of the first level temperature  
     90time evolution equation (\mdl{trasbc} module).  
     91\begin{equation} \label{Eq_sbc_trasbc_q} 
     92\frac{\partial T}{\partial t}\equiv \cdots \;+\;\left. {\frac{Q_{ns} }{\rho  
     93_o \;C_p \;e_{3t} }} \right|_{k=1} \quad 
     94\end{equation} 
     95$Q_{sr}$ is the penetrative part of the heat flux. It is applied as a 3D  
     96trends of the temperature equation (\mdl{traqsr} module) when \np{ln\_traqsr}=True. 
     97 
     98\begin{equation} \label{Eq_sbc_traqsr} 
     99\frac{\partial T}{\partial t}\equiv \cdots \;+\frac{Q_{sr} }{\rho_o C_p \,e_{3t} }\delta _k \left[ {I_w } \right] 
     100\end{equation} 
     101where $I_w$ is a non-dimensional function that describes the way the light  
     102penetrates inside the water column. It is generally a sum of decreasing  
     103exponentials (see \S\ref{TRA_qsr}). 
     104 
     105The surface freshwater budget is provided by fields: \textit{emp} and $\textit{emp}_S$ which  
     106may or may not be identical. Indeed, a surface freshwater flux has two effects:  
     107it changes the volume of the ocean and it changes the surface concentration of  
     108salt (and other tracers). Therefore it appears in the sea surface height as a volume  
     109flux, \textit{emp} (\textit{dynspg\_xxx} modules), and in the salinity time evolution equations  
     110as a concentration/dilution effect,  
     111$\textit{emp}_{S}$ (\mdl{trasbc} module).  
     112\begin{equation} \label{Eq_trasbc_emp} 
     113\begin{aligned} 
     114&\frac{\partial \eta }{\partial t}\equiv \cdots \;+\;\textit{emp}\quad  \\  
     115\\ 
     116 &\frac{\partial S}{\partial t}\equiv \cdots \;+\left. {\frac{\textit{emp}_S \;S}{e_{3t} }} \right|_{k=1} \\  
     117 \end{aligned} 
     118\end{equation}  
     119 
     120In the real ocean, $\textit{emp}=\textit{emp}_S$ and the ocean salt content is conserved,  
     121but it exist several numerical reasons why this equality should be broken.  
     122For example, when the ocean is coupled to a sea-ice model, the water exchanged between  
     123ice and ocean is slightly salty (mean sea-ice salinity is $\sim $\textit{4 psu}). In this case,  
     124$\textit{emp}_{S}$ take into account both concentration/dilution effect associated with  
     125freezing/melting and the salt flux between ice and ocean, while \textit{emp} is  
     126only the volume flux. In addition, in the current version of \NEMO, the sea-ice is  
     127assumed to be above the ocean (the so-called levitating sea-ice). Freezing/melting does  
     128not change the ocean volume (no impact on \textit{emp}) but it modifies the SSS. 
     129%gm  \colorbox{yellow}{(see {\S} on LIM sea-ice model)}. 
     130 
     131Note that SST can also be modified by a freshwater flux. Precipitation (in  
     132particular solid precipitation) may have a temperature significantly different from  
     133the SST. Due to the lack of information about the temperature of  
     134precipitation, we assume it is equal to the SST. Therefore, no  
     135concentration/dilution term appears in the temperature equation. It has to  
     136be emphasised that this absence does not mean that there is no heat flux  
     137associated with precipitation! Precipitation can change the ocean volume and thus the 
     138ocean heat content. It is therefore associated with a heat flux (not yet   
     139diagnosed in the model) \citep{Roullet_Madec_JGR00}). 
    105140 
    106141%\colorbox{yellow}{Miss: } 
    … …  
    117152%Sbcmod manage the ``providing'' (fourniture) to the ocean the 7 fields 
    118153% 
    119 %Fluxes update only each nn{\_}fsbc time step (namsbc) explain relation  
    120 %between nn{\_}fsbc and nf{\_}ice, do we define nf{\_}blk??? ? only one  
    121 %nn{\_}fsbc 
     154%Fluxes update only each nf{\_}sbc time step (namsbc) explain relation  
     155%between nf{\_}sbc and nf{\_}ice, do we define nf{\_}blk??? ? only one  
     156%nf{\_}sbc 
    122157% 
    123158%Explain here all the namlist namsbc variable{\ldots}. 
    124 %  
    125 % explain : use or not of surface currents 
    126159% 
    127160%\colorbox{yellow}{End Miss } 
    128161 
    129 The ocean model provides, at each time step, to the surface module (\mdl{sbcmod})  
    130 the surface currents, temperature and salinity.   
    131 These variables are averaged over \np{nn\_fsbc} time-step (\ref{Tab_ssm}),  
    132 and it is these averaged fields which are used to computes the surface fluxes  
    133 at a frequency of \np{nn\_fsbc} time-step. 
    134  
     162The ocean model provides the surface currents, temperature and salinity  
     163averaged over \np{nf\_sbc} time-step (\ref{Tab_ssm}).The computation of the  
     164mean is done in \mdl{sbcmod} module. 
    135165 
    136166%-------------------------------------------------TABLE--------------------------------------------------- 
    … …  
    145175\caption{  \label{Tab_ssm}    
    146176Ocean variables provided by the ocean to the surface module (SBC).  
    147 The variable are averaged over nn{\_}fsbc time step, $i.e.$ the frequency of  
     177The variable are averaged over nf{\_}sbc time step, $i.e.$ the frequency of  
    148178computation of surface fluxes.} 
    149179\end{center}   \end{table} 
    … …  
    429459%-------------------------------------------------------------------------------------------------------------- 
    430460 
    431 In some circumstances it may be useful to avoid calculating the 3D temperature, salinity and velocity fields  
    432 and simply read them in from a previous run or receive them from OASIS.   
     461In some circumstances it may be useful to avoid calculating the 3D temperature, salinity and velocity fields and simply read them in from  a previous run.   
     462Options are defined through the  \ngn{namsbc\_sas} namelist variables. 
    433463For example: 
    434464 
    435 \begin{itemize} 
    436 \item  Multiple runs of the model are required in code development to see the effect of different algorithms in 
     465\begin{enumerate} 
     466\item  Multiple runs of the model are required in code development to see the affect of different algorithms in 
    437467       the bulk formulae. 
    438468\item  The effect of different parameter sets in the ice model is to be examined. 
    439 \item  Development of sea-ice algorithms or parameterizations. 
    440 \item  spinup of the iceberg floats 
    441 \item  ocean/sea-ice simulation with both media running in parallel (\np{ln\_mixcpl}~=~\textit{true}) 
    442 \end{itemize} 
     469\end{enumerate} 
    443470 
    444471The StandAlone Surface scheme provides this utility. 
    445 Its options are defined through the \ngn{namsbc\_sas} namelist variables. 
    446472A new copy of the model has to be compiled with a configuration based on ORCA2\_SAS\_LIM. 
    447473However no namelist parameters need be changed from the settings of the previous run (except perhaps nn{\_}date0) 
    … …  
    449475Routines replaced are: 
    450476 
    451 \begin{itemize} 
    452 \item \mdl{nemogcm} : This routine initialises the rest of the model and repeatedly calls the stp time stepping routine (step.F90) 
     477\begin{enumerate} 
     478\item  \mdl{nemogcm} 
     479 
     480       This routine initialises the rest of the model and repeatedly calls the stp time stepping routine (step.F90) 
    453481       Since the ocean state is not calculated all associated initialisations have been removed. 
    454 \item  \mdl{step} : The main time stepping routine now only needs to call the sbc routine (and a few utility functions). 
    455 \item  \mdl{sbcmod} : This has been cut down and now only calculates surface forcing and the ice model required.  New surface modules 
     482\item  \mdl{step} 
     483 
     484       The main time stepping routine now only needs to call the sbc routine (and a few utility functions). 
     485\item  \mdl{sbcmod} 
     486 
     487       This has been cut down and now only calculates surface forcing and the ice model required.  New surface modules 
    456488       that can function when only the surface level of the ocean state is defined can also be added (e.g. icebergs). 
    457 \item  \mdl{daymod} : No ocean restarts are read or written (though the ice model restarts are retained), so calls to restart functions 
     489\item  \mdl{daymod} 
     490 
     491       No ocean restarts are read or written (though the ice model restarts are retained), so calls to restart functions 
    458492       have been removed.  This also means that the calendar cannot be controlled by time in a restart file, so the user 
    459493       must make sure that nn{\_}date0 in the model namelist is correct for his or her purposes. 
    460 \item  \mdl{stpctl} : Since there is no free surface solver, references to it have been removed from \rou{stp\_ctl} module. 
    461 \item  \mdl{diawri} : All 3D data have been removed from the output.  The surface temperature, salinity and velocity components (which 
     494\item  \mdl{stpctl} 
     495 
     496       Since there is no free surface solver, references to it have been removed from \rou{stp\_ctl} module. 
     497\item  \mdl{diawri} 
     498 
     499       All 3D data have been removed from the output.  The surface temperature, salinity and velocity components (which 
    462500       have been read in) are written along with relevant forcing and ice data. 
    463 \end{itemize} 
     501\end{enumerate} 
    464502 
    465503One new routine has been added: 
    466504 
    467 \begin{itemize} 
    468 \item  \mdl{sbcsas} : This module initialises the input files needed for reading temperature, salinity and velocity arrays at the surface. 
     505\begin{enumerate} 
     506\item  \mdl{sbcsas} 
     507       This module initialises the input files needed for reading temperature, salinity and velocity arrays at the surface. 
    469508       These filenames are supplied in namelist namsbc{\_}sas.  Unfortunately because of limitations with the \mdl{iom} module, 
    470509       the full 3D fields from the mean files have to be read in and interpolated in time, before using just the top level. 
    471510       Since fldread is used to read in the data, Interpolation on the Fly may be used to change input data resolution. 
    472 \end{itemize} 
    473  
    474  
    475 % Missing the description of the 2 following variables: 
    476 %   ln_3d_uve   = .true.    !  specify whether we are supplying a 3D u,v and e3 field 
    477 %   ln_read_frq = .false.    !  specify whether we must read frq or not 
    478  
    479  
     511\end{enumerate} 
    480512 
    481513% ================================================================ 
    … …  
    558590reanalysis and satellite data. They use an inertial dissipative method to compute  
    559591the turbulent transfer coefficients (momentum, sensible heat and evaporation)  
    560 from the 10 meters wind speed, air temperature and specific humidity. 
     592from the 10 metre wind speed, air temperature and specific humidity. 
    561593This \citet{Large_Yeager_Rep04} dataset is available through the  
    562594\href{http://nomads.gfdl.noaa.gov/nomads/forms/mom4/CORE.html}{GFDL web site}.  
    … …  
    593625or larger than the one of the input atmospheric fields. 
    594626 
    595 The  \np{sn\_wndi}, \np{sn\_wndj}, \np{sn\_qsr}, \np{sn\_qlw}, \np{sn\_tair},\np{sn\_humi},\np{sn\_prec}, \np{sn\_snow}, \np{sn\_tdif} parameters describe the fields and the way they have to be used (spatial and temporal interpolations).  
    596  
    597 \np{cn\_dir} is the directory of location of bulk files 
    598 \np{ln\_taudif} is the flag to specify if we use Hight Frequency (HF) tau information (.true.) or not (.false.) 
    599 \np{rn\_zqt}: is the height of humidity and temperature measurements (m) 
    600 \np{rn\_zu}: is the height of wind measurements (m) 
    601 The multiplicative factors to activate (value is 1) or deactivate (value is 0) :  
    602 \np{rn\_pfac} for precipitations (total and snow) 
    603 \np{rn\_efac} for evaporation  
    604 \np{rn\_vfac} for for ice/ocean velocities in the calculation of wind stress   
    605  
    606627% ------------------------------------------------------------------------------------------------------------- 
    607628%        CLIO Bulk formulea 
    … …  
    699720are sent to the atmospheric component. 
    700721 
    701 A generalised coupled interface has been developed.  
    702 It is currently interfaced with OASIS-3-MCT (\key{oasis3}).  
     722A generalised coupled interface has been developed. It is currently interfaced with OASIS 3 
     723(\key{oasis3}) and does not support OASIS 4 
     724\footnote{The \key{oasis4} exist. It activates portion of the code that are still under development.}.  
    703725It has been successfully used to interface \NEMO to most of the European atmospheric  
    704726GCM (ARPEGE, ECHAM, ECMWF, HadAM, HadGAM, LMDz),  
    … …  
    765787\label{SBC_tide} 
    766788 
    767 %------------------------------------------nam_tide--------------------------------------- 
     789A module is available to use the tidal potential forcing and is activated with with \key{tide}. 
     790 
     791 
     792%------------------------------------------nam_tide---------------------------------------------------- 
    768793\namdisplay{nam_tide} 
    769 %----------------------------------------------------------------------------------------- 
    770  
    771 A module is available to compute the tidal potential and use it in the momentum equation. 
    772 This option is activated when \key{tide} is defined. 
    773  
    774 Some parameters are available in namelist \ngn{nam\_tide}: 
     794%------------------------------------------------------------------------------------------------------------- 
     795 
     796Concerning the tidal potential, some parameters are available in namelist \ngn{nam\_tide}: 
    775797 
    776798- \np{ln\_tide\_pot} activate the tidal potential forcing 
    … …  
    779801 
    780802- \np{clname} is the name of constituent 
     803 
    781804 
    782805The tide is generated by the forces of gravity ot the Earth-Moon and Earth-Sun sytem; 
    … …  
    872895lowest box the river water is being added to (i.e. the total depth that river water is being added to in the model). 
    873896 
    874 %Christian: 
    875 If the depth information is not provide in the NetCDF file, it can be estimate from the runoff input file at the initial time-step, by setting the namelist parameter \np{ln\_rnf\_depth\_ini} to true. 
    876  
    877 This estimation is a simple linear relation between the runoff and a given depth :  
    878 \begin{equation}  
    879 h\_dep  = \frac{rn\_dep\_max} {rn\_rnf\_max}  rnf 
    880 \end{equation} 
    881 where  \np{rn\_dep\_max} is the given maximum depth over which the runoffs is spread,  
    882  \np{rn\_rnf\_max} is the maximum value of the runoff climatologie over the global domain 
    883 and rnf is the maximum value in time of the runoff climatology at each grid cell (computed online). 
    884  
    885 The estimated depth array can be output if needed in a NetCDF file by setting the namelist parameter \np{nn\_rnf\_depth\_file} to 1. 
    886  
    887897The mass/volume addition due to the river runoff is, at each relevant depth level, added to the horizontal divergence  
    888898(\textit{hdivn}) in the subroutine \rou{sbc\_rnf\_div} (called from \mdl{divcur}). 
    … …  
    948958\namdisplay{namsbc_isf} 
    949959%-------------------------------------------------------------------------------------------------------- 
    950 Namelist variable in \ngn{namsbc}, \np{nn\_isf}, controls the ice shelf representation used (Fig. \ref{Fig_SBC_isf}):  
    951  
    952 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    953 \begin{figure}[!h]    \begin{center} 
    954 \includegraphics[width=0.8\textwidth]{./TexFiles/Figures/Fig_SBC_isf.pdf} 
    955 \caption{ \label{Fig_SBC_isf} 
    956 Schematic for all the options available trough \np{nn\_isf}.} 
    957 \end{center}   \end{figure} 
    958 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    959  
     960Namelist variable in \ngn{namsbc}, \np{nn\_isf},  control the kind of ice shelf representation used.  
    960961\begin{description} 
    961 \item[\np{nn\_isf}~=~0] 
    962 The ice shelf routines are not used. The ice shelf melting is not computed or prescribed, the cavity have to be closed.  
    963 If needed, the ice shelf melting should be added to the runoff or the precipitation file. 
    964  
    965962\item[\np{nn\_isf}~=~1] 
    966 The ice shelf cavity is represented. The fwf and heat flux are computed. Two different bulk formula are available: 
    967    \begin{description} 
    968    \item[\np{nn\_isfblk}~=~1] 
    969    The bulk formula used to compute the melt is based the one described in \citet{Hunter2006}. 
    970         This formulation is based on a balance between the upward ocean heat flux and the latent heat flux at the ice shelf base. 
    971  
    972    \item[\np{nn\_isfblk}~=~2]  
    973    The bulk formula used to compute the melt is based the one described in \citet{Jenkins1991}. 
    974         This formulation is based on a 3 equations formulation (a heat flux budget, a salt flux budget and a linearised freezing point temperature equation). 
    975    \end{description} 
    976  
    977 For this 2 bulk formulations, there are 3 different ways to compute the exchange coeficient: 
    978    \begin{description} 
    979         \item[\np{nn\_gammablk~=~0~}] 
    980    The salt and heat exchange coefficients are constant and defined by \np{rn\_gammas0} and \np{rn\_gammat0} 
    981  
    982    \item[\np{nn\_gammablk~=~1~}] 
    983    The salt and heat exchange coefficients are velocity dependent and defined as $\np{rn\_gammas0} \times u_{*}$ and $\np{rn\_gammat0} \times u_{*}$ 
    984         where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn\_hisf\_tbl} meters). 
    985         See \citet{Jenkins2010} for all the details on this formulation. 
    986     
    987    \item[\np{nn\_gammablk~=~2~}] 
    988    The salt and heat exchange coefficients are velocity and stability dependent and defined as  
    989         $\gamma_{T,S} = \frac{u_{*}}{\Gamma_{Turb} + \Gamma^{T,S}_{Mole}}$ 
    990         where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn\_hisf\_tbl} meters),  
    991         $\Gamma_{Turb}$ the contribution of the ocean stability and  
    992         $\Gamma^{T,S}_{Mole}$ the contribution of the molecular diffusion. 
    993         See \citet{Holland1999} for all the details on this formulation. 
    994         \end{description} 
     963The ice shelf cavity is represented. The fwf and heat flux are computed.  
     964Full description, sensitivity and validation in preparation.  
    995965 
    996966\item[\np{nn\_isf}~=~2] 
    … …  
    998968The fwf is distributed along the ice shelf edge between the depth of the average grounding line (GL) 
    999969(\np{sn\_depmax\_isf}) and the base of the ice shelf along the calving front (\np{sn\_depmin\_isf}) as in (\np{nn\_isf}~=~3).  
    1000 Furthermore the fwf and heat flux are computed using the \citet{Beckmann2003} parameterisation of isf melting.  
     970Furthermore the fwf is computed using the \citet{Beckmann2003} parameterisation of isf melting.  
    1001971The effective melting length (\np{sn\_Leff\_isf}) is read from a file. 
    1002972 
    1003973\item[\np{nn\_isf}~=~3] 
    1004974A simple parameterisation of isf is used. The ice shelf cavity is not represented.  
    1005 The fwf (\np{sn\_rnfisf}) is prescribed and distributed along the ice shelf edge between the depth of the average grounding line (GL) 
    1006 (\np{sn\_depmax\_isf}) and the base of the ice shelf along the calving front (\np{sn\_depmin\_isf}).  
    1007 The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$. 
     975The fwf (\np{sn\_rnfisf}) is distributed along the ice shelf edge between the depth of the average grounding line (GL) 
     976(\np{sn\_depmax\_isf}) and the base of the ice shelf along the calving front (\np{sn\_depmin\_isf}). 
     977Full description, sensitivity and validation in preparation. 
    1008978 
    1009979\item[\np{nn\_isf}~=~4] 
    1010 The ice shelf cavity is opened. However, the fwf is not computed but specified from file \np{sn\_fwfisf}).  
    1011 The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$.\\ 
     980The ice shelf cavity is represented. However, the fwf (\np{sn\_fwfisf}) and heat flux (\np{sn\_qisf}) are  
     981not computed but specified from file.  
    1012982\end{description} 
    1013983 
    1014  
    1015 $\bullet$ \np{nn\_isf}~=~1 and \np{nn\_isf}~=~2 compute a melt rate based on the water mass properties, ocean velocities and depth. 
    1016  This flux is thus highly dependent of the model resolution (horizontal and vertical), realism of the water masses onto the shelf ...\\ 
    1017  
    1018 $\bullet$ \np{nn\_isf}~=~3 and \np{nn\_isf}~=~4 read the melt rate from a file. You have total control of the fwf forcing. 
    1019 This can be usefull if the water masses on the shelf are not realistic or the resolution (horizontal/vertical) are too  
    1020 coarse to have realistic melting or for studies where you need to control your heat and fw input.\\  
    1021  
    1022 Two namelist parameters control how the heat and fw fluxes are passed to NEMO: \np{rn\_hisf\_tbl} and \np{ln\_divisf} 
    1023 \begin{description} 
    1024 \item[\np{rn\_hisf\_tbl}] is the top boundary layer thickness as defined in \citet{Losch2008}.  
    1025 This parameter is only used if \np{nn\_isf}~=~1 or \np{nn\_isf}~=~4 
    1026 It allows you to control over which depth you want to spread the heat and fw fluxes.  
    1027  
    1028 If \np{rn\_hisf\_tbl} = 0.0, the fluxes are put in the top level whatever is its tickness.  
    1029  
    1030 If \np{rn\_hisf\_tbl} $>$ 0.0, the fluxes are spread over the first \np{rn\_hisf\_tbl} m (ie over one or several cells). 
    1031  
    1032 \item[\np{ln\_divisf}] is a flag to apply the fw flux as a volume flux or as a salt flux.  
    1033  
    1034 \np{ln\_divisf}~=~true applies the fwf as a volume flux. This volume flux is implemented with in the same way as for the runoff. 
    1035 The fw addition due to the ice shelf melting is, at each relevant depth level, added to the horizontal divergence  
    1036 (\textit{hdivn}) in the subroutine \rou{sbc\_isf\_div}, called from \mdl{divcur}.  
    1037 See the runoff section \ref{SBC_rnf} for all the details about the divergence correction.  
    1038  
    1039 \np{ln\_divisf}~=~false applies the fwf and heat flux directly on the salinity and temperature tendancy. 
    1040  
    1041 \item[\np{ln\_conserve}] is a flag for \np{nn\_isf}~=~1. A conservative boundary layer scheme as described in \citet{Jenkins2001}  
    1042 is used if \np{ln\_conserve}=true. It takes into account the fact that the melt water is at freezing T and needs to be warm up to ocean temperature.  
    1043 It is only relevant for \np{ln\_divisf}~=~false.  
    1044 If \np{ln\_divisf}~=~true, \np{ln\_conserve} has to be set to false to avoid a double counting of the contribution.  
    1045   
    1046 \end{description} 
     984\np{nn\_isf}~=~1 and \np{nn\_isf}~=~2 compute a melt rate based on the water masse properties, ocean velocities and depth. 
     985 This flux is thus highly dependent of the model resolution (horizontal and vertical), realism of the water masse onto the shelf ... 
     986 
     987\np{nn\_isf}~=~3 and \np{nn\_isf}~=~4 read the melt rate and heat flux from a file. You have total control of the fwf scenario. 
     988 
     989 This can be usefull if the water masses on the shelf are not realistic or the resolution (horizontal/vertical) are too  
     990coarse to have realistic melting or for sensitivity studies where you want to control your input.  
     991Full description, sensitivity and validation in preparation.  
     992 
     993There is 2 ways to apply the fwf to NEMO. The first possibility (\np{ln\_divisf}~=~false) applied the fwf 
     994 and heat flux directly on the salinity and temperature tendancy. The second possibility (\np{ln\_divisf}~=~true) 
     995 apply the fwf as for the runoff fwf (see \S\ref{SBC_rnf}). The mass/volume addition due to the ice shelf melting is, 
     996 at each relevant depth level, added to the horizontal divergence (\textit{hdivn}) in the subroutine \rou{sbc\_isf\_div}  
     997(called from \mdl{divcur}).  
    1047998% 
    1048999% ================================================================ 
    10491000%        Handling of icebergs 
    10501001% ================================================================ 
    1051 \section{Handling of icebergs (ICB)} 
     1002\section{ Handling of icebergs (ICB) } 
    10521003\label{ICB_icebergs} 
    10531004%------------------------------------------namberg---------------------------------------------------- 
    … …  
    10551006%------------------------------------------------------------------------------------------------------------- 
    10561007 
    1057 Icebergs are modelled as lagrangian particles in NEMO \citep{Marsh_GMD2015}. 
    1058 Their physical behaviour is controlled by equations as described in \citet{Martin_Adcroft_OM10} ). 
    1059 (Note that the authors kindly provided a copy of their code to act as a basis for implementation in NEMO). 
    1060 Icebergs are initially spawned into one of ten classes which have specific mass and thickness as described  
    1061 in the \ngn{namberg} namelist:  
     1008Icebergs are modelled as lagrangian particles in NEMO. 
     1009Their physical behaviour is controlled by equations as described in  \citet{Martin_Adcroft_OM10} ). 
     1010(Note that the authors kindly provided a copy of their code to act as a basis for implementation in NEMO.) 
     1011Icebergs are initially spawned into one of ten classes which have specific mass and thickness as described in the \ngn{namberg} namelist:  
    10621012\np{rn\_initial\_mass} and \np{rn\_initial\_thickness}. 
    10631013Each class has an associated scaling (\np{rn\_mass\_scaling}), which is an integer representing how many icebergs  
    … …  
    12431193The presence at the sea surface of an ice covered area modifies all the fluxes  
    12441194transmitted to the ocean. There are several way to handle sea-ice in the system  
    1245 depending on the value of the \np{nn\_ice} namelist parameter found in \ngn{namsbc} namelist.   
     1195depending on the value of the \np{nn{\_}ice} namelist parameter.   
    12461196\begin{description} 
    12471197\item[nn{\_}ice = 0]  there will never be sea-ice in the computational domain.  
    … …  
    13181268% ------------------------------------------------------------------------------------------------------------- 
    13191269\subsection   [Neutral drag coefficient from external wave model (\textit{sbcwave})] 
    1320               {Neutral drag coefficient from external wave model (\mdl{sbcwave})} 
     1270                        {Neutral drag coefficient from external wave model (\mdl{sbcwave})} 
    13211271\label{SBC_wave} 
    13221272%------------------------------------------namwave---------------------------------------------------- 
    13231273\namdisplay{namsbc_wave} 
    13241274%------------------------------------------------------------------------------------------------------------- 
    1325  
    1326 In order to read a neutral drag coeff, from an external data source ($i.e.$ a wave model), the  
    1327 logical variable \np{ln\_cdgw} in \ngn{namsbc} namelist must be set to \textit{true}.  
     1275\begin{description} 
     1276 
     1277\item [??] In order to read a neutral drag coeff, from an external data source (i.e. a wave model), the  
     1278logical variable \np{ln\_cdgw} 
     1279 in $namsbc$ namelist must be defined ${.true.}$.  
    13281280The \mdl{sbcwave} module containing the routine \np{sbc\_wave} reads the 
    13291281namelist \ngn{namsbc\_wave} (for external data names, locations, frequency, interpolation and all  
    13301282the miscellanous options allowed by Input Data generic Interface see \S\ref{SBC_input})  
    1331 and a 2D field of neutral drag coefficient.  
    1332 Then using the routine TURB\_CORE\_1Z or TURB\_CORE\_2Z, and starting from the neutral drag coefficent provided,  
    1333 the drag coefficient is computed according to stable/unstable conditions of the air-sea interface following \citet{Large_Yeager_Rep04}. 
    1334  
     1283and a 2D field of neutral drag coefficient. Then using the routine  
     1284TURB\_CORE\_1Z or TURB\_CORE\_2Z, and starting from the neutral drag coefficent provided, the drag coefficient is computed according  
     1285to stable/unstable conditions of the air-sea interface following \citet{Large_Yeager_Rep04}. 
     1286 
     1287\end{description} 
    13351288 
    13361289% Griffies doc: 
    1337 % When running ocean-ice simulations, we are not explicitly representing land processes,  
    1338 % such as rivers, catchment areas, snow accumulation, etc. However, to reduce model drift,  
    1339 % it is important to balance the hydrological cycle in ocean-ice models.  
    1340 % We thus need to prescribe some form of global normalization to the precipitation minus evaporation plus river runoff.  
    1341 % The result of the normalization should be a global integrated zero net water input to the ocean-ice system over  
    1342 % a chosen time scale.  
    1343 %How often the normalization is done is a matter of choice. In mom4p1, we choose to do so at each model time step,  
    1344 % so that there is always a zero net input of water to the ocean-ice system.  
    1345 % Others choose to normalize over an annual cycle, in which case the net imbalance over an annual cycle is used  
    1346 % to alter the subsequent year�s water budget in an attempt to damp the annual water imbalance.  
    1347 % Note that the annual budget approach may be inappropriate with interannually varying precipitation forcing.  
    1348 % When running ocean-ice coupled models, it is incorrect to include the water transport between the ocean  
    1349 % and ice models when aiming to balance the hydrological cycle.  
    1350 % The reason is that it is the sum of the water in the ocean plus ice that should be balanced when running ocean-ice models,  
    1351 % not the water in any one sub-component. As an extreme example to illustrate the issue,  
    1352 % consider an ocean-ice model with zero initial sea ice. As the ocean-ice model spins up,  
    1353 % there should be a net accumulation of water in the growing sea ice, and thus a net loss of water from the ocean.  
    1354 % The total water contained in the ocean plus ice system is constant, but there is an exchange of water between  
    1355 % the subcomponents. This exchange should not be part of the normalization used to balance the hydrological cycle  
    1356 % in ocean-ice models.  
    1357  
    1358  
     1290% When running ocean-ice simulations, we are not explicitly representing land processes, such as rivers, catchment areas, snow accumulation, etc. However, to reduce model drift, it is important to balance the hydrological cycle in ocean-ice models. We thus need to prescribe some form of global normalization to the precipitation minus evaporation plus river runoff. The result of the normalization should be a global integrated zero net water input to the ocean-ice system over a chosen time scale.  
     1291%How often the normalization is done is a matter of choice. In mom4p1, we choose to do so at each model time step, so that there is always a zero net input of water to the ocean-ice system. Others choose to normalize over an annual cycle, in which case the net imbalance over an annual cycle is used to alter the subsequent year�s water budget in an attempt to damp the annual water imbalance. Note that the annual budget approach may be inappropriate with interannually varying precipitation forcing.  
     1292%When running ocean-ice coupled models, it is incorrect to include the water transport between the ocean and ice models when aiming to balance the hydrological cycle. The reason is that it is the sum of the water in the ocean plus ice that should be balanced when running ocean-ice models, not the water in any one sub-component. As an extreme example to illustrate the issue, consider an ocean-ice model with zero initial sea ice. As the ocean-ice model spins up, there should be a net accumulation of water in the growing sea ice, and thus a net loss of water from the ocean. The total water contained in the ocean plus ice system is constant, but there is an exchange of water between the subcomponents. This exchange should not be part of the normalization used to balance the hydrological cycle in ocean-ice models.  
     1293 
     1294 

  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Chap_STO.tex

    r6617 r6625  
    55\label{STO} 
    66 
    7 Authors: P.-A. Bouttier 
     7\minitoc 
    88 
    9 \minitoc 
    109 
    1110\newpage 
    1211$\ $\newline    % force a new line 
    13  
    14 The stochastic parametrization module aims to explicitly simulate uncertainties in the model. More particularly, \cite{Brankart_OM2013} has shown that, because of the nonlinearity of the seawater equation of state, unresolved scales represent a major source of uncertainties in the computation of the large scale horizontal density gradient (from T/S large scale fields), and that the impact of these uncertainties can be simulated by random processes representing unresolved T/S fluctuations. 
    15  
    16 The stochastic formulation of the equation of state can be written as: 
    17 \begin{equation} 
    18  \label{eq:eos_sto} 
    19   \rho = \frac{1}{2} \sum_{i=1}^m\{ \rho[T+\Delta T_i,S+\Delta S_i,p_o(z)] + \rho[T-\Delta T_i,S-\Delta S_i,p_o(z)] \} 
    20 \end{equation} 
    21 where $p_o(z)$ is the reference pressure depending on the depth and $\Delta T_i$ and $\Delta S_i$ are a set of T/S perturbations defined as the scalar product of the respective local T/S gradients with random walks $\mathbf{\xi}$: 
    22 \begin{equation} 
    23  \label{eq:sto_pert} 
    24  \Delta T_i = \mathbf{\xi}_i \cdot \nabla T \qquad \hbox{and} \qquad \Delta S_i = \mathbf{\xi}_i \cdot \nabla S 
    25 \end{equation} 
    26 $\mathbf{\xi}_i$ are produced by a first-order autoregressive processes (AR-1) with a parametrized decorrelation time scale, and horizontal and vertical standard deviations $\sigma_s$. $\mathbf{\xi}$ are uncorrelated over the horizontal and fully correlated along the vertical. 
    27  
    28  
    29 \section{Stochastic processes} 
    30 \label{STO_the_details} 
    31  
    32 The starting point of our implementation of stochastic parameterizations 
    33 in NEMO is to observe that many existing parameterizations are based 
    34 on autoregressive processes, which are used as a basic source of randomness 
    35 to transform a deterministic model into a probabilistic model. 
    36 A generic approach is thus to add one single new module in NEMO, 
    37 generating processes with appropriate statistics 
    38 to simulate each kind of uncertainty in the model 
    39 (see \cite{Brankart_al_GMD2015} for more details). 
    40  
    41 In practice, at every model grid point, independent Gaussian autoregressive 
    42 processes~$\xi^{(i)},\,i=1,\ldots,m$ are first generated 
    43 using the same basic equation: 
    44  
    45 \begin{equation} 
    46 \label{eq:autoreg} 
    47 \xi^{(i)}_{k+1} = a^{(i)} \xi^{(i)}_k + b^{(i)} w^{(i)} + c^{(i)} 
    48 \end{equation} 
    49  
    50 \noindent 
    51 where $k$ is the index of the model timestep; and 
    52 $a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are parameters defining 
    53 the mean ($\mu^{(i)}$) standard deviation ($\sigma^{(i)}$) 
    54 and correlation timescale ($\tau^{(i)}$) of each process: 
    55  
    56 \begin{itemize} 
    57 \item for order~1 processes, $w^{(i)}$ is a Gaussian white noise, 
    58 with zero mean and standard deviation equal to~1, and the parameters 
    59 $a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are given by: 
    60  
    61 \begin{equation} 
    62 \label{eq:ord1} 
    63 \left\{ 
    64 \begin{array}{l} 
    65 a^{(i)} = \varphi \\ 
    66 b^{(i)} = \sigma^{(i)} \sqrt{ 1 - \varphi^2 }  
    67  \qquad\qquad\mbox{with}\qquad\qquad 
    68 \varphi = \exp \left( - 1 / \tau^{(i)} \right) \\ 
    69 c^{(i)} = \mu^{(i)} \left( 1 - \varphi \right) \\ 
    70 \end{array} 
    71 \right. 
    72 \end{equation} 
    73  
    74 \item for order~$n>1$ processes, $w^{(i)}$ is an order~$n-1$ autoregressive process, 
    75 with zero mean, standard deviation equal to~$\sigma^{(i)}$; correlation timescale 
    76 equal to~$\tau^{(i)}$; and the parameters 
    77 $a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are given by: 
    78  
    79 \begin{equation} 
    80 \label{eq:ord2} 
    81 \left\{ 
    82 \begin{array}{l} 
    83 a^{(i)} = \varphi \\ 
    84 b^{(i)} = \frac{n-1}{2(4n-3)} \sqrt{ 1 - \varphi^2 }  
    85  \qquad\qquad\mbox{with}\qquad\qquad 
    86 \varphi = \exp \left( - 1 / \tau^{(i)} \right) \\ 
    87 c^{(i)} = \mu^{(i)} \left( 1 - \varphi \right) \\ 
    88 \end{array} 
    89 \right. 
    90 \end{equation} 
    91  
    92 \end{itemize} 
    93  
    94 \noindent 
    95 In this way, higher order processes can be easily generated recursively using the same piece of code implementing Eq.~(\ref{eq:autoreg}), and using succesively processes from order $0$ to~$n-1$ as~$w^{(i)}$. 
    96 The parameters in Eq.~(\ref{eq:ord2}) are computed so that this recursive application 
    97 of Eq.~(\ref{eq:autoreg}) leads to processes with the required standard deviation 
    98 and correlation timescale, with the additional condition that 
    99 the $n-1$ first derivatives of the autocorrelation function 
    100 are equal to zero at~$t=0$, so that the resulting processes 
    101 become smoother and smoother as $n$ is increased. 
    102  
    103 Overall, this method provides quite a simple and generic way of generating a wide class of stochastic processes. However, this also means that new model parameters are needed to specify each of these stochastic processes. As in any parameterization of lacking physics, a very important issues then to tune these new parameters using either first principles, model simulations, or real-world observations. 
    104  
    105 \section{Implementation details} 
    106 \label{STO_thech_details} 
    107 The computer code implementing stochastic parametrisations is made of one single FORTRAN module, 
    108 with 3 public routines to be called by the model (in our case, NEMO): 
    109  
    110 The first routine ({sto\_par}) is a direct implementation of Eq.~(\ref{eq:autoreg}), 
    111 applied at each model grid point (in 2D or 3D), 
    112 and called at each model time step ($k$) to update 
    113 every autoregressive process ($i=1,\ldots,m$). 
    114 This routine also includes a filtering operator, applied to $w^{(i)}$, 
    115 to introduce a spatial correlation between the stochastic processes. 
    116  
    117 The second routine ({sto\_par\_init}) 
    118 is an initialization routine mainly dedicated 
    119 to the computation of parameters $a^{(i)}, b^{(i)}, c^{(i)}$ 
    120 for each autoregressive process, as a function of the statistical properties 
    121 required by the model user (mean, standard deviation, time correlation, 
    122 order of the process,\ldots). Parameters for the processes can be specified through the following namelist parameters: 
    123 \begin{alltt} 
    124 \tiny 
    125 \begin{verbatim} 
    126    nn_sto_eos = 1                ! number of independent random walks  
    127    rn_eos_stdxy = 1.4            ! random walk horz. standard deviation (in grid points) 
    128    rn_eos_stdz  = 0.7            ! random walk vert. standard deviation (in grid points) 
    129    rn_eos_tcor  = 1440.0         ! random walk time correlation (in timesteps) 
    130    nn_eos_ord  = 1               ! order of autoregressive processes 
    131    nn_eos_flt  = 0               ! passes of Laplacian filter 
    132    rn_eos_lim  = 2.0             ! limitation factor (default = 3.0) 
    133 \end{verbatim} 
    134 \end{alltt} 
    135 This routine also includes the initialization (seeding) 
    136 of the random number generator. 
    137  
    138 The third routine ({sto\_rst\_write}) writes a ``restart file'' 
    139 with the current value of all autoregressive processes 
    140 to allow restarting a simulation from where it has been interrupted. 
    141 This file also contains the current state of the random number generator. 
    142 In case of a restart, this file is then read by the initialization routine 
    143 ({sto\_par\_init}), so that the simulation can continue exactly 
    144 as if it was not interrupted. 
    145 Restart capabilities of the module are driven by the following namelist parameters: 
    146 \begin{alltt} 
    147 \tiny 
    148 \begin{verbatim} 
    149    ln_rststo = .false.           ! start from mean parameter (F) or from restart file (T) 
    150    ln_rstseed = .true.           ! read seed of RNG from restart file 
    151    cn_storst_in  = "restart_sto" !  suffix of stochastic parameter restart file (input) 
    152    cn_storst_out = "restart_sto" !  suffix of stochastic parameter restart file (output) 
    153 \end{verbatim} 
    154 \end{alltt} 
    155  
    156 In the particular case of the stochastic equation of state, there is also an additional module ({sto\_pts}) implementing Eq~\ref{eq:sto_pert} and specific piece of code in the equation of state implementing Eq~\ref{eq:eos_sto}. 
    157  
    158  

  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Chap_TRA.tex

    r6617 r6625  
    11% ================================================================ 
    2 % Chapter 1 ——— Ocean Tracers (TRA) 
     2% Chapter 1 � Ocean Tracers (TRA) 
    33% ================================================================ 
    44\chapter{Ocean Tracers (TRA)} 
    … …  
    3636(BBL) parametrisation, and an internal damping (DMP) term. The terms QSR,  
    3737BBC, BBL and DMP are optional. The external forcings and parameterisations  
    38 require complex inputs and complex calculations ($e.g.$ bulk formulae, estimation  
     38require complex inputs and complex calculations (e.g. bulk formulae, estimation  
    3939of mixing coefficients) that are carried out in the SBC, LDF and ZDF modules and  
    4040described in chapters \S\ref{SBC}, \S\ref{LDF} and  \S\ref{ZDF}, respectively.  
    41 Note that \mdl{tranpc}, the non-penetrative convection module, although  
    42 located in the NEMO/OPA/TRA directory as it directly modifies the tracer fields,  
    43 is described with the model vertical physics (ZDF) together with other available  
    44 parameterization of convection. 
     41Note that \mdl{tranpc}, the non-penetrative convection module,  although  
     42(temporarily) located in the NEMO/OPA/TRA directory, is described with the  
     43model vertical physics (ZDF). 
     44%%% 
     45\gmcomment{change the position of eosbn2 in the reference code} 
     46%%% 
    4547 
    4648In the present chapter we also describe the diagnostic equations used to compute  
    47 the sea-water properties (density, Brunt-V\"{a}is\"{a}l\"{a} frequency, specific heat and  
     49the sea-water properties (density, Brunt-Vais\"{a}l\"{a} frequency, specific heat and  
    4850freezing point with associated modules \mdl{eosbn2} and \mdl{phycst}). 
    4951 
    … …  
    5456found in the \textit{trattt} or \textit{trattt\_xxx} module, in the NEMO/OPA/TRA directory. 
    5557 
    56 The user has the option of extracting each tendency term on the RHS of the tracer  
    57 equation for output (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}~=~true), as described in Chap.~\ref{DIA}. 
     58The user has the option of extracting each tendency term on the rhs of the tracer  
     59equation for output (\key{trdtra} is defined), as described in Chap.~\ref{MISC}. 
    5860 
    5961$\ $\newline    % force a new ligne 
    … …  
    123125\end{description} 
    124126In all cases, this boundary condition retains local conservation of tracer.  
    125 Global conservation is obtained in non-linear free surface case,  
    126 but \textit{not} in the linear free surface case. Nevertheless, in the latter case,  
    127 it is achieved to a good approximation since the non-conservative  
     127Global conservation is obtained in both rigid-lid and non-linear free surface  
     128cases, but not in the linear free surface case. Nevertheless, in the latter 
     129case, it is achieved to a good approximation since the non-conservative  
    128130term is the product of the time derivative of the tracer and the free surface  
    129131height, two quantities that are not correlated (see \S\ref{PE_free_surface},  
    … …  
    131133 
    132134The velocity field that appears in (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_zco})  
    133 is the centred (\textit{now}) \textit{effective} ocean velocity, $i.e.$ the \textit{eulerian} velocity 
    134 (see Chap.~\ref{DYN}) plus the eddy induced velocity (\textit{eiv})  
    135 and/or the mixed layer eddy induced velocity (\textit{eiv}) 
    136 when those parameterisations are used (see Chap.~\ref{LDF}). 
     135is the centred (\textit{now}) \textit{eulerian} ocean velocity (see Chap.~\ref{DYN}).  
     136When eddy induced velocity (\textit{eiv}) parameterisation is used it is the \textit{now}  
     137\textit{effective} velocity ($i.e.$ the sum of the eulerian and eiv velocities) which is used. 
    137138 
    138139The choice of an advection scheme is made in the \textit{\ngn{nam\_traadv}} namelist, by  
    … …  
    145146 
    146147Note that  
    147 (1) cen2 and TVD schemes require an explicit diffusion  
     148(1) cen2, cen4 and TVD schemes require an explicit diffusion  
    148149operator while the other schemes are diffusive enough so that they do not  
    149150require additional diffusion ;  
    150 (2) cen2, MUSCL2, and UBS are not \textit{positive} schemes 
     151(2) cen2, cen4, MUSCL2, and UBS are not \textit{positive} schemes 
    151152\footnote{negative values can appear in an initially strictly positive tracer field  
    152153which is advected} 
    … …  
    188189temperature is close to the freezing point). 
    189190This combined scheme has been included for specific grid points in the ORCA2  
    190 configuration only. This is an obsolescent feature as the recommended  
     191and ORCA4 configurations only. This is an obsolescent feature as the recommended  
    191192advection scheme for the ORCA configuration is TVD (see  \S\ref{TRA_adv_tvd}). 
    192193 
    … …  
    195196have this order of accuracy. \gmcomment{Note also that ... blah, blah} 
    196197 
     198% ------------------------------------------------------------------------------------------------------------- 
     199%        4nd order centred scheme   
     200% ------------------------------------------------------------------------------------------------------------- 
     201\subsection   [$4^{nd}$ order centred scheme (cen4) (\np{ln\_traadv\_cen4})] 
     202           {$4^{nd}$ order centred scheme (cen4) (\np{ln\_traadv\_cen4}=true)} 
     203\label{TRA_adv_cen4} 
     204 
     205In the $4^{th}$ order formulation (to be implemented), tracer values are  
     206evaluated at velocity points as a $4^{th}$ order interpolation, and thus depend on  
     207the four neighbouring $T$-points. For example, in the $i$-direction: 
     208\begin{equation} \label{Eq_tra_adv_cen4} 
     209\tau _u^{cen4}  
     210=\overline{   T - \frac{1}{6}\,\delta _i \left[ \delta_{i+1/2}[T] \,\right]   }^{\,i+1/2} 
     211\end{equation} 
     212 
     213Strictly speaking, the cen4 scheme is not a $4^{th}$ order advection scheme  
     214but a $4^{th}$ order evaluation of advective fluxes, since the divergence of  
     215advective fluxes \eqref{Eq_tra_adv} is kept at $2^{nd}$ order. The phrase ``$4^{th}$  
     216order scheme'' used in oceanographic literature is usually associated  
     217with the scheme presented here. Introducing a \textit{true} $4^{th}$ order advection  
     218scheme is feasible but, for consistency reasons, it requires changes in the  
     219discretisation of the tracer advection together with changes in both the  
     220continuity equation and the momentum advection terms.   
     221 
     222A direct consequence of the pseudo-fourth order nature of the scheme is that  
     223it is not non-diffusive, i.e. the global variance of a tracer is not preserved using  
     224\textit{cen4}. Furthermore, it must be used in conjunction with an explicit  
     225diffusion operator to produce a sensible solution. The time-stepping is also  
     226performed using a leapfrog scheme in conjunction with an Asselin time-filter,  
     227so $T$ in (\ref{Eq_tra_adv_cen4}) is the \textit{now} tracer. 
     228 
     229At a $T$-grid cell adjacent to a boundary (coastline, bottom and surface), an  
     230additional hypothesis must be made to evaluate $\tau _u^{cen4}$. This  
     231hypothesis usually reduces the order of the scheme. Here we choose to set  
     232the gradient of $T$ across the boundary to zero. Alternative conditions can be  
     233specified, such as a reduction to a second order scheme for these near boundary  
     234grid points. 
    197235 
    198236% ------------------------------------------------------------------------------------------------------------- 
    … …  
    232270used for the diffusive part.  
    233271 
    234 An additional option has been added controlled by \np{ln\_traadv\_tvd\_zts}.  
    235 By setting this logical to true, a TVD scheme is used on both horizontal and vertical direction,  
    236 but on the latter, a split-explicit time stepping is used, with 5 sub-timesteps.  
    237 This option can be useful when the value of the timestep is limited by vertical advection \citep{Lemarie_OM2015}.  
    238 Note that in this case, a similar split-explicit time stepping should be used on  
    239 vertical advection of momentum to ensure a better stability (see \np{ln\_dynzad\_zts} in \S\ref{DYN_zad}). 
    240  
    241  
    242272% ------------------------------------------------------------------------------------------------------------- 
    243273%        MUSCL scheme   
    … …  
    266296 
    267297For an ocean grid point adjacent to land and where the ocean velocity is  
    268 directed toward land, two choices are available: an upstream flux (\np{ln\_traadv\_muscl}=true)  
    269 or a second order flux (\np{ln\_traadv\_muscl2}=true).  
    270 Note that the latter choice does not ensure the \textit{positive} character of the scheme.  
    271 Only the former can be used on both active and passive tracers.  
    272 The two MUSCL schemes are implemented in the \mdl{traadv\_tvd} and \mdl{traadv\_tvd2} modules. 
    273  
    274 Note that when using np{ln\_traadv\_msc\_ups}~=~true in addition to \np{ln\_traadv\_muscl}=true,  
    275 the MUSCL fluxes are replaced by upstream fluxes in vicinity of river mouths. 
     298directed toward land, two choices are available: an upstream flux  
     299(\np{ln\_traadv\_muscl}=true) or a second order flux  
     300(\np{ln\_traadv\_muscl2}=true). Note that the latter choice does not ensure  
     301the \textit{positive} character of the scheme. Only the former can be used  
     302on both active and passive tracers. The two MUSCL schemes are implemented  
     303in the \mdl{traadv\_tvd} and \mdl{traadv\_tvd2} modules. 
    276304 
    277305% ------------------------------------------------------------------------------------------------------------- 
    … …  
    388416direction (as for the UBS case) should be implemented to restore this property. 
    389417 
     418 
     419% ------------------------------------------------------------------------------------------------------------- 
     420%        PPM scheme   
     421% ------------------------------------------------------------------------------------------------------------- 
     422\subsection   [Piecewise Parabolic Method (PPM) (\np{ln\_traadv\_ppm})] 
     423         {Piecewise Parabolic Method (PPM) (\np{ln\_traadv\_ppm}=true)} 
     424\label{TRA_adv_ppm} 
     425 
     426The Piecewise Parabolic Method (PPM) proposed by Colella and Woodward (1984)  
     427\sgacomment{reference?} 
     428is based on a quadradic piecewise construction. Like the QCK scheme, it is associated  
     429with the ULTIMATE QUICKEST limiter \citep{Leonard1991}. It has been implemented  
     430in \NEMO by G. Reffray (MERCATOR-ocean) but is not yet offered in the reference  
     431version 3.3. 
    390432 
    391433% ================================================================ 
    … …  
    422464surfaces is given by:  
    423465\begin{equation} \label{Eq_tra_ldf_lap} 
    424 D_T^{lT} =\frac{1}{b_t} \left( \; 
     466D_T^{lT} =\frac{1}{b_tT} \left( \; 
    425467   \delta _{i}\left[ A_u^{lT} \; \frac{e_{2u}\,e_{3u}}{e_{1u}} \;\delta _{i+1/2} [T] \right]  
    426468+ \delta _{j}\left[ A_v^{lT} \;  \frac{e_{1v}\,e_{3v}}{e_{2v}} \;\delta _{j+1/2} [T] \right]  \;\right) 
    … …  
    619661the thickness of the top model layer.  
    620662 
    621 Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components  
    622 ($i.e.$ atmosphere, sea-ice, land), the change in the heat and salt content of the surface layer  
    623 of the ocean is due both to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$)  
    624 and to the heat and salt content of the mass exchange. They are both included directly in $Q_{ns}$,  
    625 the surface heat flux, and $F_{salt}$, the surface salt flux (see \S\ref{SBC} for further details). 
    626 By doing this, the forcing formulation is the same for any tracer (including temperature and salinity). 
     663Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components ($i.e.$ atmosphere, sea-ice, land), 
     664the change in the heat and salt content of the surface layer of the ocean is due both  
     665to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$) 
     666 and to the heat and salt content of the mass exchange. 
     667\sgacomment{ the following does not apply to the release to which this documentation is  
     668attached and so should not be included .... 
     669In a forthcoming release, these two parts, computed in the surface module (SBC), will be included directly 
     670in $Q_{ns}$, the surface heat flux and $F_{salt}$, the surface salt flux. 
     671The specification of these fluxes is further detailed in the SBC chapter (see \S\ref{SBC}).  
     672This change will provide a forcing formulation which is the same for any tracer (including temperature and salinity). 
     673  
     674In the current version, the situation is a little bit more complicated. } 
    627675 
    628676The surface module (\mdl{sbcmod}, see \S\ref{SBC}) provides the following  
    … …  
    631679$\bullet$ $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface  
    632680(i.e. the difference between the total surface heat flux and the fraction of the short wave flux that  
    633 penetrates into the water column, see \S\ref{TRA_qsr}) plus the heat content associated with  
    634 of the mass exchange with the atmosphere and lands. 
    635  
    636 $\bullet$ $\textit{sfx}$, the salt flux resulting from ice-ocean mass exchange (freezing, melting, ridging...) 
    637  
    638 $\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation)  
    639  and possibly with the sea-ice and ice-shelves. 
    640  
    641 $\bullet$ \textit{rnf}, the mass flux associated with runoff  
    642 (see \S\ref{SBC_rnf} for further detail of how it acts on temperature and salinity tendencies) 
    643  
    644 $\bullet$ \textit{fwfisf}, the mass flux associated with ice shelf melt, (see \S\ref{SBC_isf} for further details  
    645 on how the ice shelf melt is computed and applied).\\ 
    646  
    647 In the non-linear free surface case (\key{vvl} is defined), the surface boundary condition  
    648 on temperature and salinity is applied as follows: 
     681penetrates into the water column, see \S\ref{TRA_qsr}) 
     682 
     683$\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation) 
     684 
     685$\bullet$ $\textit{emp}_S$, an equivalent mass flux taking into account the effect of ice-ocean mass exchange 
     686 
     687$\bullet$ \textit{rnf}, the mass flux associated with runoff (see \S\ref{SBC_rnf} for further detail of how it acts on temperature and salinity tendencies) 
     688 
     689The $\textit{emp}_S$ field is not simply the budget of evaporation-precipitation+freezing-melting because  
     690the sea-ice is not currently embedded in the ocean but levitates above it. There is no mass 
     691exchanged between the sea-ice and the ocean. Instead we only take into account the salt 
     692flux associated with the non-zero salinity of sea-ice, and the concentration/dilution effect 
     693due to the freezing/melting (F/M) process. These two parts of the forcing are then converted into  
     694an equivalent mass flux given by $\textit{emp}_S - \textit{emp}$. As a result of this mess,  
     695the surface boundary condition on temperature and salinity is applied as follows: 
     696 
     697In the nonlinear free surface case (\key{vvl} is defined): 
    649698\begin{equation} \label{Eq_tra_sbc} 
    650699\begin{aligned} 
    651  &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }  &\overline{ Q_{ns}       }^t  & \\  
    652 & F^S =\frac{ 1 }{\rho _o  \,      \left. e_{3t} \right|_{k=1} }  &\overline{ \textit{sfx} }^t   & \\    
     700 &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }    
     701           &\overline{ \left( Q_{ns} - \textit{emp}\;C_p\,\left. T \right|_{k=1} \right) }^t  & \\  
     702% 
     703& F^S =\frac{ 1 }{\rho _o \,\left. e_{3t} \right|_{k=1} }  
     704           &\overline{ \left( (\textit{emp}_S - \textit{emp})\;\left. S \right|_{k=1}  \right) }^t   & \\    
     705 \end{aligned} 
     706\end{equation}  
     707 
     708In the linear free surface case (\key{vvl} not defined): 
     709\begin{equation} \label{Eq_tra_sbc_lin} 
     710\begin{aligned} 
     711 &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }  &\overline{ Q_{ns} }^t  & \\  
     712% 
     713& F^S =\frac{ 1 }{\rho _o \,\left. e_{3t} \right|_{k=1} }  
     714           &\overline{ \left( \textit{emp}_S\;\left. S \right|_{k=1}  \right) }^t   & \\    
    653715 \end{aligned} 
    654716\end{equation}  
    … …  
    657719divergence of odd and even time step (see \S\ref{STP}). 
    658720 
    659 In the linear free surface case (\key{vvl} is \textit{not} defined),  
    660 an additional term has to be added on both temperature and salinity.  
    661 On temperature, this term remove the heat content associated with mass exchange 
    662 that has been added to $Q_{ns}$. On salinity, this term mimics the concentration/dilution effect that 
    663 would have resulted from a change in the volume of the first level. 
    664 The resulting surface boundary condition is applied as follows: 
    665 \begin{equation} \label{Eq_tra_sbc_lin} 
    666 \begin{aligned} 
    667  &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }    
    668            &\overline{ \left( Q_{ns} - \textit{emp}\;C_p\,\left. T \right|_{k=1} \right) }^t  & \\  
    669 % 
    670 & F^S =\frac{ 1 }{\rho _o \,\left. e_{3t} \right|_{k=1} }  
    671            &\overline{ \left( \;\textit{sfx} - \textit{emp} \;\left. S \right|_{k=1}  \right) }^t   & \\    
    672  \end{aligned} 
    673 \end{equation}  
    674 Note that an exact conservation of heat and salt content is only achieved with non-linear free surface.  
    675 In the linear free surface case, there is a small imbalance. The imbalance is larger  
     721The two set of equations, \eqref{Eq_tra_sbc} and \eqref{Eq_tra_sbc_lin}, are obtained  
     722by assuming that the temperature of precipitation and evaporation are equal to 
     723the ocean surface temperature and that their salinity is zero. Therefore, the heat content 
     724of the \textit{emp} budget must be added to the temperature equation in the variable volume case,  
     725while it does not appear in the constant volume case. Similarly, the \textit{emp} budget affects  
     726the ocean surface salinity in the constant volume case (through the concentration dilution effect) 
     727while it does not appears explicitly in the variable volume case since salinity change will be 
     728induced by volume change. In both constant and variable volume cases, surface salinity  
     729will change with ice-ocean salt flux and F/M flux (both contained in $\textit{emp}_S - \textit{emp}$) without mass exchanges. 
     730 
     731Note that the concentration/dilution effect due to F/M is computed using 
     732a constant ice salinity as well as a constant ocean salinity.  
     733This approximation suppresses the correlation between \textit{SSS}  
     734and F/M flux, allowing the ice-ocean salt exchanges to be conservative. 
     735Indeed, if this approximation is not made, even if the F/M budget is zero  
     736on average over the whole ocean domain and over the seasonal cycle,  
     737the associated salt flux is not zero, since sea-surface salinity and F/M flux are  
     738intrinsically correlated (high \textit{SSS} are found where freezing is  
     739strong whilst low \textit{SSS} is usually associated with high melting areas). 
     740 
     741Even using this approximation, an exact conservation of heat and salt content  
     742is only achieved in the variable volume case. In the constant volume case,  
     743there is a small imbalance associated with the product $(\partial_t\eta - \textit{emp}) * \textit{SSS}$. 
     744Nevertheless, the salt content variation is quite small and will not induce 
     745a long term drift as there is no physical reason for $(\partial_t\eta - \textit{emp})$  
     746and \textit{SSS} to be correlated \citep{Roullet_Madec_JGR00}.  
     747Note that, while quite small, the imbalance in the constant volume case is larger  
    676748than the imbalance associated with the Asselin time filter \citep{Leclair_Madec_OM09}.  
    677 This is the reason why the modified filter is not applied in the linear free surface case (see \S\ref{STP}). 
     749This is the reason why the modified filter is not applied in the constant volume case. 
    678750 
    679751% ------------------------------------------------------------------------------------------------------------- 
    … …  
    749821($i.e.$ the inverses of the extinction length scales) are tabulated over 61 nonuniform  
    750822chlorophyll classes ranging from 0.01 to 10 g.Chl/L (see the routine \rou{trc\_oce\_rgb}  
    751 in \mdl{trc\_oce} module). Four types of chlorophyll can be chosen in the RGB formulation: 
    752 \begin{description}  
    753 \item[\np{nn\_chdta}=0]  
    754 a constant 0.05 g.Chl/L value everywhere ;  
    755 \item[\np{nn\_chdta}=1]   
    756 an observed time varying chlorophyll deduced from satellite surface ocean color measurement  
    757 spread uniformly in the vertical direction ;  
    758 \item[\np{nn\_chdta}=2]   
    759 same as previous case except that a vertical profile of chlorophyl is used.  
    760 Following \cite{Morel_Berthon_LO89}, the profile is computed from the local surface chlorophyll value ; 
    761 \item[\np{ln\_qsr\_bio}=true]   
    762 simulated time varying chlorophyll by TOP biogeochemical model.  
    763 In this case, the RGB formulation is used to calculate both the phytoplankton  
    764 light limitation in PISCES or LOBSTER and the oceanic heating rate.  
    765 \end{description}  
     823in \mdl{trc\_oce} module). Three types of chlorophyll can be chosen in the RGB formulation: 
     824(1) a constant 0.05 g.Chl/L value everywhere (\np{nn\_chdta}=0) ; (2) an observed  
     825time varying chlorophyll (\np{nn\_chdta}=1) ; (3) simulated time varying chlorophyll 
     826by TOP biogeochemical model (\np{ln\_qsr\_bio}=true). In the latter case, the RGB  
     827formulation is used to calculate both the phytoplankton light limitation in PISCES  
     828or LOBSTER and the oceanic heating rate.  
     829 
    766830The trend in \eqref{Eq_tra_qsr} associated with the penetration of the solar radiation  
    767831is added to the temperature trend, and the surface heat flux is modified in routine \mdl{traqsr}.  
    … …  
    795859\label{TRA_bbc} 
    796860%--------------------------------------------nambbc-------------------------------------------------------- 
    797 \namdisplay{nambbc} 
     861\namdisplay{namtra_bbc} 
    798862%-------------------------------------------------------------------------------------------------------------- 
    799863%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    10391103\subsection[DMP\_TOOLS]{Generating resto.nc using DMP\_TOOLS} 
    10401104 
    1041 DMP\_TOOLS can be used to generate a netcdf file containing the restoration coefficient $\gamma$.  
    1042 Note that in order to maintain bit comparison with previous NEMO versions DMP\_TOOLS must be compiled  
    1043 and run on the same machine as the NEMO model. A mesh\_mask.nc file for the model configuration is required as an input.  
    1044 This can be generated by carrying out a short model run with the namelist parameter \np{nn\_msh} set to 1.  
    1045 The namelist parameter \np{ln\_tradmp} will also need to be set to .false. for this to work.  
    1046 The \nl{nam\_dmp\_create} namelist in the DMP\_TOOLS directory is used to specify options for the restoration coefficient. 
     1105DMP\_TOOLS can be used to generate a netcdf file containing the restoration coefficient $\gamma$. Note that in order to maintain bit comparison with previous NEMO versions DMP\_TOOLS must be compiled and run on the same machine as the NEMO model. A mesh\_mask.nc file for the model configuration is required as an input. This can be generated by carrying out a short model run with the namelist parameter \np{nn\_msh} set to 1. The namelist parameter \np{ln\_tradmp} will also need to be set to .false. for this to work. The \nl{nam\_dmp\_create} namelist in the DMP\_TOOLS directory is used to specify options for the restoration coefficient. 
    10471106 
    10481107%--------------------------------------------nam_dmp_create------------------------------------------------- 
    … …  
    10521111\np{cp\_cfg}, \np{cp\_cpz}, \np{jp\_cfg} and \np{jperio} specify the model configuration being used and should be the same as specified in \nl{namcfg}. The variable \nl{lzoom} is used to specify that the damping is being used as in case \textit{a} above to provide boundary conditions to a zoom configuration. In the case of the arctic or antarctic zoom configurations this includes some specific treatment. Otherwise damping is applied to the 6 grid points along the ocean boundaries. The open boundaries are specified by the variables \np{lzoom\_n}, \np{lzoom\_e}, \np{lzoom\_s}, \np{lzoom\_w} in the \nl{nam\_zoom\_dmp} name list. 
    10531112 
    1054 The remaining switch namelist variables determine the spatial variation of the restoration coefficient in non-zoom configurations.  
    1055 \np{ln\_full\_field} specifies that newtonian damping should be applied to the whole model domain.  
    1056 \np{ln\_med\_red\_seas} specifies grid specific restoration coefficients in the Mediterranean Sea  
    1057 for the ORCA4, ORCA2 and ORCA05 configurations.  
    1058 If \np{ln\_old\_31\_lev\_code} is set then the depth variation of the coeffients will be specified as  
    1059 a function of the model number. This option is included to allow backwards compatability of the ORCA2 reference  
    1060 configurations with previous model versions.  
    1061 \np{ln\_coast} specifies that the restoration coefficient should be reduced near to coastlines.  
    1062 This option only has an effect if \np{ln\_full\_field} is true.  
    1063 \np{ln\_zero\_top\_layer} specifies that the restoration coefficient should be zero in the surface layer.  
    1064 Finally \np{ln\_custom} specifies that the custom module will be called.  
    1065 This module is contained in the file custom.F90 and can be edited by users. For example damping could be applied in a specific region. 
    1066  
    1067 The restoration coefficient can be set to zero in equatorial regions by specifying a positive value of \np{nn\_hdmp}.  
    1068 Equatorward of this latitude the restoration coefficient will be zero with a smooth transition to  
    1069 the full values of a 10$^{\circ}$ latitud band.  
    1070 This is often used because of the short adjustment time scale in the equatorial region  
    1071 \citep{Reverdin1991, Fujio1991, Marti_PhD92}. The time scale associated with the damping depends on the depth as a  
    1072 hyperbolic tangent, with \np{rn\_surf} as surface value, \np{rn\_bot} as bottom value and a transition depth of \np{rn\_dep}.   
     1113The remaining switch namelist variables determine the spatial variation of the restoration coefficient in non-zoom configurations. \np{ln\_full\_field} specifies that newtonian damping should be applied to the whole model domain. \np{ln\_med\_red\_seas} specifies grid specific restoration coefficients in the Mediterranean Sea for the ORCA4, ORCA2 and ORCA05 configurations. If \np{ln\_old\_31\_lev\_code} is set then the depth variation of the coeffients will be specified as a function of the model number. This option is included to allow backwards compatability of the ORCA2 reference configurations with previous model versions. \np{ln\_coast} specifies that the restoration coefficient should be reduced near to coastlines. This option only has an effect if \np{ln\_full\_field} is true. \np{ln\_zero\_top\_layer} specifies that the restoration coefficient should be zero in the surface layer. Finally \np{ln\_custom} specifies that the custom module will be called. This module is contained in the file custom.F90 and can be edited by users. For example damping could be applied in a specific region. 
     1114 
     1115The restoration coefficient can be set to zero in equatorial regions by specifying a positive value of \np{nn\_hdmp}. Equatorward of this latitude the restoration coefficient will be zero with a smooth transition to the full values of a 10$^{\circ}$ latitud band. This is often used because of the short adjustment time scale in the equatorial region \citep{Reverdin1991, Fujio1991, Marti_PhD92}. The time scale associated with the damping depends on the depth as a hyperbolic tangent, with \np{rn\_surf} as surface value, \np{rn\_bot} as bottom value and a transition depth of \np{rn\_dep}.   
    10731116 
    10741117% ================================================================ 
    … …  
    11241167%        Equation of State 
    11251168% ------------------------------------------------------------------------------------------------------------- 
    1126 \subsection{Equation Of Seawater (\np{nn\_eos} = -1, 0, or 1)} 
     1169\subsection{Equation of State (\np{nn\_eos} = 0, 1 or 2)} 
    11271170\label{TRA_eos} 
    11281171 
    1129 The Equation Of Seawater (EOS) is an empirical nonlinear thermodynamic relationship  
    1130 linking seawater density, $\rho$, to a number of state variables,  
    1131 most typically temperature, salinity and pressure.  
    1132 Because density gradients control the pressure gradient force through the hydrostatic balance,  
    1133 the equation of state provides a fundamental bridge between the distribution of active tracers  
    1134 and the fluid dynamics. Nonlinearities of the EOS are of major importance, in particular  
    1135 influencing the circulation through determination of the static stability below the mixed layer,  
    1136 thus controlling rates of exchange between the atmosphere  and the ocean interior \citep{Roquet_JPO2015}.  
    1137 Therefore an accurate EOS based on either the 1980 equation of state (EOS-80, \cite{UNESCO1983})  
    1138 or TEOS-10 \citep{TEOS10} standards should be used anytime a simulation of the real  
    1139 ocean circulation is attempted \citep{Roquet_JPO2015}.  
    1140 The use of TEOS-10 is highly recommended because  
    1141 \textit{(i)} it is the new official EOS,  
    1142 \textit{(ii)} it is more accurate, being based on an updated database of laboratory measurements, and  
    1143 \textit{(iii)} it uses Conservative Temperature and Absolute Salinity (instead of potential temperature  
    1144 and practical salinity for EOS-980, both variables being more suitable for use as model variables  
    1145 \citep{TEOS10, Graham_McDougall_JPO13}.  
    1146 EOS-80 is an obsolescent feature of the NEMO system, kept only for backward compatibility. 
    1147 For process studies, it is often convenient to use an approximation of the EOS. To that purposed,  
    1148 a simplified EOS (S-EOS) inspired by \citet{Vallis06} is also available. 
    1149  
    1150 In the computer code, a density anomaly, $d_a= \rho / \rho_o - 1$,  
    1151 is computed, with $\rho_o$ a reference density. Called \textit{rau0}  
    1152 in the code, $\rho_o$ is set in \mdl{phycst} to a value of $1,026~Kg/m^3$.  
     1172It is necessary to know the equation of state for the ocean very accurately  
     1173to determine stability properties (especially the Brunt-Vais\"{a}l\"{a} frequency),  
     1174particularly in the deep ocean. The ocean seawater volumic mass, $\rho$,  
     1175abusively called density, is a non linear empirical function of \textit{in situ}  
     1176temperature, salinity and pressure. The reference equation of state is that  
     1177defined by the Joint Panel on Oceanographic Tables and Standards  
     1178\citep{UNESCO1983}. It was the standard equation of state used in early  
     1179releases of OPA. However, even though this computation is fully vectorised,  
     1180it is quite time consuming ($15$ to $20${\%} of the total CPU time) since  
     1181it requires the prior computation of the \textit{in situ} temperature from the  
     1182model \textit{potential} temperature using the \citep{Bryden1973} polynomial  
     1183for adiabatic lapse rate and a $4^th$ order Runge-Kutta integration scheme.  
     1184Since OPA6, we have used the \citet{JackMcD1995} equation of state for  
     1185seawater instead. It allows the computation of the \textit{in situ} ocean density  
     1186directly as a function of \textit{potential} temperature relative to the surface  
     1187(an \NEMO variable), the practical salinity (another \NEMO variable) and the  
     1188pressure (assuming no pressure variation along geopotential surfaces, $i.e.$  
     1189the pressure in decibars is approximated by the depth in meters).  
     1190Both the \citet{UNESCO1983} and \citet{JackMcD1995} equations of state  
     1191have exactly the same except that the values of the various coefficients have  
     1192been adjusted by \citet{JackMcD1995} in order to directly use the \textit{potential}  
     1193temperature instead of the \textit{in situ} one. This reduces the CPU time of the  
     1194\textit{in situ} density computation to about $3${\%} of the total CPU time,  
     1195while maintaining a quite accurate equation of state. 
     1196 
     1197In the computer code, a \textit{true} density anomaly, $d_a= \rho / \rho_o - 1$,  
     1198is computed, with $\rho_o$ a reference volumic mass. Called \textit{rau0}  
     1199in the code, $\rho_o$ is defined in \mdl{phycst}, and a value of $1,035~Kg/m^3$.  
    11531200This is a sensible choice for the reference density used in a Boussinesq ocean  
    11541201climate model, as, with the exception of only a small percentage of the ocean,  
    1155 density in the World Ocean varies by no more than 2$\%$ from that value \citep{Gill1982}. 
    1156  
    1157 Options are defined through the  \ngn{nameos} namelist variables, and in particular \np{nn\_eos}  
    1158 which controls the EOS used (=-1 for TEOS10 ; =0 for EOS-80 ; =1 for S-EOS). 
    1159 \begin{description} 
    1160  
    1161 \item[\np{nn\_eos}$=-1$] the polyTEOS10-bsq equation of seawater \citep{Roquet_OM2015} is used.   
    1162 The accuracy of this approximation is comparable to the TEOS-10 rational function approximation,  
    1163 but it is optimized for a boussinesq fluid and the polynomial expressions have simpler  
    1164 and more computationally efficient expressions for their derived quantities  
    1165 which make them more adapted for use in ocean models.  
    1166 Note that a slightly higher precision polynomial form is now used replacement of the TEOS-10  
    1167 rational function approximation for hydrographic data analysis  \citep{TEOS10}.  
    1168 A key point is that conservative state variables are used:  
    1169 Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: $\degres C$, notation: $\Theta$). 
    1170 The pressure in decibars is approximated by the depth in meters.  
    1171 With TEOS10, the specific heat capacity of sea water, $C_p$, is a constant. It is set to  
    1172 $C_p=3991.86795711963~J\,Kg^{-1}\,\degres K^{-1}$, according to \citet{TEOS10}. 
    1173  
    1174 Choosing polyTEOS10-bsq implies that the state variables used by the model are  
    1175 $\Theta$ and $S_A$. In particular, the initial state deined by the user have to be given as  
    1176 \textit{Conservative} Temperature and \textit{Absolute} Salinity.  
    1177 In addition, setting \np{ln\_useCT} to \textit{true} convert the Conservative SST to potential SST  
    1178 prior to either computing the air-sea and ice-sea fluxes (forced mode)  
    1179 or sending the SST field to the atmosphere (coupled mode). 
    1180  
    1181 \item[\np{nn\_eos}$=0$] the polyEOS80-bsq equation of seawater is used. 
    1182 It takes the same polynomial form as the polyTEOS10, but the coefficients have been optimized  
    1183 to accurately fit EOS80 (Roquet, personal comm.). The state variables used in both the EOS80  
    1184 and the ocean model are:  
    1185 the Practical Salinity ((unit: psu, notation: $S_p$)) and Potential Temperature (unit: $\degres C$, notation: $\theta$). 
    1186 The pressure in decibars is approximated by the depth in meters.   
    1187 With thsi EOS, the specific heat capacity of sea water, $C_p$, is a function of temperature,  
    1188 salinity and pressure \citep{UNESCO1983}. Nevertheless, a severe assumption is made in order to  
    1189 have a heat content ($C_p T_p$) which is conserved by the model: $C_p$ is set to a constant  
    1190 value, the TEOS10 value.  
    1191   
    1192 \item[\np{nn\_eos}$=1$] a simplified EOS (S-EOS) inspired by \citet{Vallis06} is chosen,  
    1193 the coefficients of which has been optimized to fit the behavior of TEOS10 (Roquet, personal comm.)  
    1194 (see also \citet{Roquet_JPO2015}). It provides a simplistic linear representation of both  
    1195 cabbeling and thermobaricity effects which is enough for a proper treatment of the EOS  
    1196 in theoretical studies \citep{Roquet_JPO2015}. 
     1202density in the World Ocean varies by no more than 2$\%$ from $1,035~kg/m^3$  
     1203\citep{Gill1982}. 
     1204 
     1205Options are defined through the  \ngn{nameos} namelist variables. 
     1206The default option (namelist parameter \np{nn\_eos}=0) is the \citet{JackMcD1995}  
     1207equation of state. Its use is highly recommended. However, for process studies,  
     1208it is often convenient to use a linear approximation of the density. 
    11971209With such an equation of state there is no longer a distinction between  
    1198 \textit{conservative} and \textit{potential} temperature, as well as between \textit{absolute}  
    1199 and \textit{practical} salinity. 
    1200 S-EOS takes the following expression: 
    1201 \begin{equation} \label{Eq_tra_S-EOS} 
     1210\textit{in situ} and \textit{potential} density and both cabbeling and thermobaric 
     1211effects are removed. 
     1212Two linear formulations are available: a function of $T$ only (\np{nn\_eos}=1)  
     1213and a function of both $T$ and $S$ (\np{nn\_eos}=2): 
     1214\begin{equation} \label{Eq_tra_eos_linear} 
    12021215\begin{split} 
    1203   d_a(T,S,z)  =  ( & - a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * T_a  \\ 
    1204                                 & + b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * S_a  \\ 
    1205                                 & - \nu \; T_a \; S_a \;  ) \; / \; \rho_o                     \\ 
    1206   with \ \  T_a = T-10  \; ;  & \;  S_a = S-35  \; ;\;  \rho_o = 1026~Kg/m^3 
     1216  d_a(T)       &=  \rho (T)      /  \rho_o   - 1     =  \  0.0285         -  \alpha   \;T     \\  
     1217  d_a(T,S)    &=  \rho (T,S)   /  \rho_o   - 1     =  \  \beta \; S       -  \alpha   \;T     
    12071218\end{split} 
    12081219\end{equation}  
    1209 where the computer name of the coefficients as well as their standard value are given in \ref{Tab_SEOS}. 
    1210 In fact, when choosing S-EOS, various approximation of EOS can be specified simply by changing  
    1211 the associated coefficients.  
    1212 Setting to zero the two thermobaric coefficients ($\mu_1$, $\mu_2$) remove thermobaric effect from S-EOS. 
    1213 setting to zero the three cabbeling coefficients ($\lambda_1$, $\lambda_2$, $\nu$) remove cabbeling effect from S-EOS. 
    1214 Keeping non-zero value to $a_0$ and $b_0$ provide a linear EOS function of T and S. 
    1215  
    1216 \end{description} 
    1217  
    1218  
    1219 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    1220 \begin{table}[!tb] 
    1221 \begin{center} \begin{tabular}{|p{26pt}|p{72pt}|p{56pt}|p{136pt}|} 
    1222 \hline 
    1223 coeff.   & computer name   & S-EOS     &  description                      \\ \hline 
    1224 $a_0$       & \np{rn\_a0}     & 1.6550 $10^{-1}$ &  linear thermal expansion coeff.    \\ \hline 
    1225 $b_0$       & \np{rn\_b0}     & 7.6554 $10^{-1}$ &  linear haline  expansion coeff.    \\ \hline 
    1226 $\lambda_1$ & \np{rn\_lambda1}& 5.9520 $10^{-2}$ &  cabbeling coeff. in $T^2$          \\ \hline 
    1227 $\lambda_2$ & \np{rn\_lambda2}& 5.4914 $10^{-4}$ &  cabbeling coeff. in $S^2$       \\ \hline 
    1228 $\nu$       & \np{rn\_nu}     & 2.4341 $10^{-3}$ &  cabbeling coeff. in $T \, S$       \\ \hline 
    1229 $\mu_1$     & \np{rn\_mu1}    & 1.4970 $10^{-4}$ &  thermobaric coeff. in T         \\ \hline 
    1230 $\mu_2$     & \np{rn\_mu2}    & 1.1090 $10^{-5}$ &  thermobaric coeff. in S            \\ \hline 
    1231 \end{tabular} 
    1232 \caption{ \label{Tab_SEOS} 
    1233 Standard value of S-EOS coefficients. } 
    1234 \end{center} 
    1235 \end{table} 
    1236 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    1237  
    1238  
    1239 % ------------------------------------------------------------------------------------------------------------- 
    1240 %        Brunt-V\"{a}is\"{a}l\"{a} Frequency 
    1241 % ------------------------------------------------------------------------------------------------------------- 
    1242 \subsection{Brunt-V\"{a}is\"{a}l\"{a} Frequency (\np{nn\_eos} = 0, 1 or 2)} 
     1220where $\alpha$ and $\beta$ are the thermal and haline expansion  
     1221coefficients, and $\rho_o$, the reference volumic mass, $rau0$.  
     1222($\alpha$ and $\beta$ can be modified through the \np{rn\_alpha} and  
     1223\np{rn\_beta} namelist variables). Note that when $d_a$ is a function  
     1224of $T$ only (\np{nn\_eos}=1), the salinity is a passive tracer and can be  
     1225used as such. 
     1226 
     1227% ------------------------------------------------------------------------------------------------------------- 
     1228%        Brunt-Vais\"{a}l\"{a} Frequency 
     1229% ------------------------------------------------------------------------------------------------------------- 
     1230\subsection{Brunt-Vais\"{a}l\"{a} Frequency (\np{nn\_eos} = 0, 1 or 2)} 
    12431231\label{TRA_bn2} 
    12441232 
    1245 An accurate computation of the ocean stability (i.e. of $N$, the brunt-V\"{a}is\"{a}l\"{a} 
    1246  frequency) is of paramount importance as determine the ocean stratification and  
    1247  is used in several ocean parameterisations (namely TKE, GLS, Richardson number dependent  
    1248  vertical diffusion, enhanced vertical diffusion, non-penetrative convection, tidal mixing  
    1249  parameterisation, iso-neutral diffusion). In particular, $N^2$ has to be computed at the local pressure  
    1250  (pressure in decibar being approximated by the depth in meters). The expression for $N^2$  
    1251  is given by:  
     1233An accurate computation of the ocean stability (i.e. of $N$, the brunt-Vais\"{a}l\"{a} 
     1234 frequency) is of paramount importance as it is used in several ocean  
     1235 parameterisations (namely TKE, KPP, Richardson number dependent  
     1236 vertical diffusion, enhanced vertical diffusion, non-penetrative convection,  
     1237 iso-neutral diffusion). In particular, one must be aware that $N^2$ has to  
     1238 be computed with an \textit{in situ} reference. The expression for $N^2$  
     1239 depends on the type of equation of state used (\np{nn\_eos} namelist parameter). 
     1240 
     1241For \np{nn\_eos}=0 (\citet{JackMcD1995} equation of state), the \citet{McDougall1987}  
     1242polynomial expression is used (with the pressure in decibar approximated by  
     1243the depth in meters):  
    12521244\begin{equation} \label{Eq_tra_bn2} 
     1245N^2 = \frac{g}{e_{3w}} \; \beta   \  
     1246      \left(  \alpha / \beta \ \delta_{k+1/2}[T]     - \delta_{k+1/2}[S]   \right)  
     1247\end{equation}  
     1248where $\alpha$ and $\beta$ are the thermal and haline expansion coefficients.  
     1249They are a function of  $\overline{T}^{\,k+1/2},\widetilde{S}=\overline{S}^{\,k+1/2} - 35.$,  
     1250and  $z_w$, with $T$ the \textit{potential} temperature and $\widetilde{S}$ a salinity anomaly.  
     1251Note that both $\alpha$ and $\beta$ depend on \textit{potential}  
     1252temperature and salinity which are averaged at $w$-points prior  
     1253to the computation instead of being computed at $T$-points and  
     1254then averaged to $w$-points. 
     1255 
     1256When a linear equation of state is used (\np{nn\_eos}=1 or 2,  
     1257\eqref{Eq_tra_bn2} reduces to: 
     1258\begin{equation} \label{Eq_tra_bn2_linear} 
    12531259N^2 = \frac{g}{e_{3w}} \left(   \beta \;\delta_{k+1/2}[S] - \alpha \;\delta_{k+1/2}[T]   \right) 
    12541260\end{equation}  
    1255 where $(T,S) = (\Theta, S_A)$ for TEOS10, $= (\theta, S_p)$ for TEOS-80, or $=(T,S)$ for S-EOS,  
    1256 and, $\alpha$ and $\beta$ are the thermal and haline expansion coefficients.  
    1257 The coefficients are a polynomial function of temperature, salinity and depth which expression  
    1258 depends on the chosen EOS. They are computed through \textit{eos\_rab}, a \textsc{Fortran}  
    1259 function that can be found in \mdl{eosbn2}. 
     1261where $\alpha$ and $\beta $ are the constant coefficients used to  
     1262defined the linear equation of state \eqref{Eq_tra_eos_linear}. 
     1263 
     1264% ------------------------------------------------------------------------------------------------------------- 
     1265%        Specific Heat 
     1266% ------------------------------------------------------------------------------------------------------------- 
     1267\subsection    [Specific Heat (\textit{phycst})] 
     1268         {Specific Heat (\mdl{phycst})} 
     1269\label{TRA_adv_ldf} 
     1270 
     1271The specific heat of sea water, $C_p$, is a function of temperature, salinity  
     1272and pressure \citep{UNESCO1983}. It is only used in the model to convert  
     1273surface heat fluxes into surface temperature increase and so the pressure  
     1274dependence is neglected. The dependence on $T$ and $S$ is weak.  
     1275For example, with $S=35~psu$, $C_p$ increases from $3989$ to $4002$  
     1276when $T$ varies from -2~\degres C to 31~\degres C. Therefore, $C_p$ has  
     1277been chosen as a constant: $C_p=4.10^3~J\,Kg^{-1}\,\degres K^{-1}$.  
     1278Its value is set in \mdl{phycst} module.  
     1279 
    12601280 
    12611281% ------------------------------------------------------------------------------------------------------------- 
    … …  
    12781298sea water ($i.e.$ referenced to the surface $p=0$), thus the pressure dependent  
    12791299terms in \eqref{Eq_tra_eos_fzp} (last term) have been dropped. The freezing 
    1280 point is computed through \textit{eos\_fzp}, a \textsc{Fortran} function that can be found  
     1300point is computed through \textit{tfreez}, a \textsc{Fortran} function that can be found  
    12811301in \mdl{eosbn2}.   
    12821302 
    … …  
    12881308\label{TRA_zpshde} 
    12891309 
    1290 \gmcomment{STEVEN: to be consistent with earlier discussion of differencing and averaging operators,  
    1291                    I've changed "derivative" to "difference" and "mean" to "average"} 
    1292  
    1293 With partial cells (\np{ln\_zps}=true) at bottom and top (\np{ln\_isfcav}=true), in general, tracers in horizontally  
     1310\gmcomment{STEVEN: to be consistent with earlier discussion of differencing and averaging operators, I've changed "derivative" to "difference" and "mean" to "average"} 
     1311 
     1312With partial bottom cells (\np{ln\_zps}=true), in general, tracers in horizontally  
    12941313adjacent cells live at different depths. Horizontal gradients of tracers are needed  
    12951314for horizontal diffusion (\mdl{traldf} module) and for the hydrostatic pressure  
    12961315gradient (\mdl{dynhpg} module) to be active.  
    12971316\gmcomment{STEVEN from gm : question: not sure of  what -to be active- means} 
    1298  
    12991317Before taking horizontal gradients between the tracers next to the bottom, a linear  
    13001318interpolation in the vertical is used to approximate the deeper tracer as if it actually  
    … …  
    13721390\gmcomment{gm :   this last remark has to be done} 
    13731391%%% 
    1374  
    1375 If under ice shelf seas opened (\np{ln\_isfcav}=true), the partial cell properties  
    1376 at the top are computed in the same way as for the bottom. Some extra variables are,  
    1377 however, computed to reduce the flow generated at the top and bottom if $z*$ coordinates activated. 
    1378 The extra variables calculated and used by \S\ref{DYN_hpg_isf} are: 
    1379  
    1380 $\bullet$ $\overline{T}_k^{\,i+1/2}$ as described in \eqref{Eq_zps_hde} 
    1381  
    1382 $\bullet$ $\delta _{i+1/2} Z_{T_k} = \widetilde {Z}^{\,i}_{T_k}-Z^{\,i}_{T_k}$ to compute  
    1383 the pressure gradient correction term used by \eqref{Eq_dynhpg_sco} in \S\ref{DYN_hpg_isf}, 
    1384  with $\widetilde {Z}_{T_k}$ the depth of the point $\widetilde {T}_{k}$ in case of $z^*$ coordinates  
    1385 (this term = 0 in z-coordinates) 

  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Chap_ZDF.tex

    r6617 r6625  
    3333points, respectively (see \S\ref{TRA_zdf} and \S\ref{DYN_zdf}). These  
    3434coefficients can be assumed to be either constant, or a function of the local  
    35 Richardson number, or computed from a turbulent closure model (TKE, GLS or KPP formulation).  
    36 The computation of these coefficients is initialized in the \mdl{zdfini} module  
    37 and performed in the \mdl{zdfric}, \mdl{zdftke}, \mdl{zdfgls} or \mdl{zdfkpp} modules.  
    38 The trends due to the vertical momentum and tracer diffusion, including the surface forcing,  
    39 are computed and added to the general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively.  
     35Richardson number, or computed from a turbulent closure model (either  
     36TKE or KPP formulation). The computation of these coefficients is initialized  
     37in the \mdl{zdfini} module and performed in the \mdl{zdfric}, \mdl{zdftke} or  
     38\mdl{zdfkpp} modules. The trends due to the vertical momentum and tracer  
     39diffusion, including the surface forcing, are computed and added to the  
     40general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively.  
    4041These trends can be computed using either a forward time stepping scheme  
    4142(namelist parameter \np{ln\_zdfexp}=true) or a backward time stepping  
    … …  
    261262\end{equation} 
    262263 
    263 At the ocean surface, a non zero length scale is set through the  \np{rn\_mxl0} namelist  
     264At the ocean surface, a non zero length scale is set through the  \np{rn\_lmin0} namelist  
    264265parameter. Usually the surface scale is given by $l_o = \kappa \,z_o$  
    265266where $\kappa = 0.4$ is von Karman's constant and $z_o$ the roughness  
    266267parameter of the surface. Assuming $z_o=0.1$~m \citep{Craig_Banner_JPO94}  
    267 leads to a 0.04~m, the default value of \np{rn\_mxl0}. In the ocean interior  
     268leads to a 0.04~m, the default value of \np{rn\_lsurf}. In the ocean interior  
    268269a minimum length scale is set to recover the molecular viscosity when $\bar{e}$  
    269270reach its minimum value ($1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ). 
    … …  
    294295As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$,  
    295296with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}~=~67.83 corresponds  
    296 to $\alpha_{CB} = 100$. Further setting  \np{ln\_mxl0} to true applies \eqref{ZDF_Lsbc}  
    297 as surface boundary condition on length scale, with $\beta$ hard coded to the Stacey's value. 
     297to $\alpha_{CB} = 100$. further setting  \np{ln\_lsurf} to true applies \eqref{ZDF_Lsbc}  
     298as surface boundary condition on length scale, with $\beta$ hard coded to the Stacet's value. 
    298299Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters)  
    299300is applied on surface $\bar{e}$ value. 
    … …  
    354355%--------------------------------------------------------------% 
    355356 
    356 Vertical mixing parameterizations commonly used in ocean general circulation models  
    357 tend to produce mixed-layer depths that are too shallow during summer months and windy conditions. 
    358 This bias is particularly acute over the Southern Ocean.  
    359 To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme  \cite{Rodgers_2014}.  
    360 The parameterization is an empirical one, $i.e.$ not derived from theoretical considerations,  
    361 but rather is meant to account for observed processes that affect the density structure of  
    362 the ocean’s planetary boundary layer that are not explicitly captured by default in the TKE scheme  
    363 ($i.e.$ near-inertial oscillations and ocean swells and waves). 
    364  
    365 When using this parameterization ($i.e.$ when \np{nn\_etau}~=~1), the TKE input to the ocean ($S$)  
    366 imposed by the winds in the form of near-inertial oscillations, swell and waves is parameterized  
    367 by \eqref{ZDF_Esbc} the standard TKE surface boundary condition, plus a depth depend one given by: 
    368 \begin{equation}  \label{ZDF_Ehtau} 
    369 S = (1-f_i) \; f_r \; e_s \; e^{-z / h_\tau}  
    370 \end{equation} 
    371 where  
    372 $z$ is the depth,   
    373 $e_s$ is TKE surface boundary condition,  
    374 $f_r$ is the fraction of the surface TKE that penetrate in the ocean,  
    375 $h_\tau$ is a vertical mixing length scale that controls exponential shape of the penetration,  
    376 and $f_i$ is the ice concentration (no penetration if $f_i=1$, that is if the ocean is entirely  
    377 covered by sea-ice). 
    378 The value of $f_r$, usually a few percents, is specified through \np{rn\_efr} namelist parameter.  
    379 The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np{nn\_etau}~=~0)  
    380 or a latitude dependent value (varying from 0.5~m at the Equator to a maximum value of 30~m  
    381 at high latitudes (\np{nn\_etau}~=~1).  
    382  
    383 Note that two other option existe, \np{nn\_etau}~=~2, or 3. They correspond to applying  
    384 \eqref{ZDF_Ehtau} only at the base of the mixed layer, or to using the high frequency part  
    385 of the stress to evaluate the fraction of TKE that penetrate the ocean.  
    386 Those two options are obsolescent features introduced for test purposes. 
    387 They will be removed in the next release.  
    388  
    389  
     357To be add here a description of "penetration of TKE" and the associated namelist parameters 
     358 \np{nn\_etau}, \np{rn\_efr} and \np{nn\_htau}. 
    390359 
    391360% from Burchard et al OM 2008 :  
    392 % the most critical process not reproduced by statistical turbulence models is the activity of  
    393 % internal waves and their interaction with turbulence. After the Reynolds decomposition,  
    394 % internal waves are in principle included in the RANS equations, but later partially  
    395 % excluded by the hydrostatic assumption and the model resolution.  
    396 % Thus far, the representation of internal wave mixing in ocean models has been relatively crude  
    397 % (e.g. Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002). 
     361% the most critical process not reproduced by statistical turbulence models is the activity of internal waves and their interaction with turbulence. After the Reynolds decomposition, internal waves are in principle included in the RANS equations, but later partially excluded by the hydrostatic assumption and the model resolution. Thus far, the representation of internal wave mixing in ocean models has been relatively crude (e.g. Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002). 
    398362 
    399363 
    … …  
    622586Options are defined through the  \ngn{namzdf\_kpp} namelist variables. 
    623587 
    624 Note that KPP is an obsolescent feature of the \NEMO system.  
    625 It will be removed in the next release (v3.7 and followings). 
     588\colorbox{yellow}{Add a description of KPP here.} 
    626589 
    627590 
    … …  
    673636 
    674637Options are defined through the  \ngn{namzdf} namelist variables. 
    675 The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}~=~\textit{true}.  
     638The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}=true.  
    676639It is applied at each \np{nn\_npc} time step and mixes downwards instantaneously  
    677640the statically unstable portion of the water column, but only until the density  
    … …  
    681644(Fig. \ref{Fig_npc}): starting from the top of the ocean, the first instability is  
    682645found. Assume in the following that the instability is located between levels  
    683 $k$ and $k+1$. The temperature and salinity in the two levels are  
     646$k$ and $k+1$. The potential temperature and salinity in the two levels are  
    684647vertically mixed, conserving the heat and salt contents of the water column.  
    685648The new density is then computed by a linear approximation. If the new  
    … …  
    701664\citep{Madec_al_JPO91, Madec_al_DAO91, Madec_Crepon_Bk91}. 
    702665 
    703 The current implementation has been modified in order to deal with any non linear  
    704 equation of seawater (L. Brodeau, personnal communication).  
    705 Two main differences have been introduced compared to the original algorithm:  
    706 $(i)$ the stability is now checked using the Brunt-V\"{a}is\"{a}l\"{a} frequency  
    707 (not the the difference in potential density) ;  
    708 $(ii)$ when two levels are found unstable, their thermal and haline expansion coefficients  
    709 are vertically mixed in the same way their temperature and salinity has been mixed. 
    710 These two modifications allow the algorithm to perform properly and accurately  
    711 with TEOS10 or EOS-80 without having to recompute the expansion coefficients at each  
    712 mixing iteration. 
     666Note that in the current implementation of this algorithm presents several  
     667limitations. First, potential density referenced to the sea surface is used to  
     668check whether the density profile is stable or not. This is a strong  
     669simplification which leads to large errors for realistic ocean simulations.  
     670Indeed, many water masses of the world ocean, especially Antarctic Bottom 
     671Water, are unstable when represented in surface-referenced potential density.  
     672The scheme will erroneously mix them up. Second, the mixing of potential  
     673density is assumed to be linear. This assures the convergence of the algorithm  
     674even when the equation of state is non-linear. Small static instabilities can thus  
     675persist due to cabbeling: they will be treated at the next time step.  
     676Third, temperature and salinity, and thus density, are mixed, but the  
     677corresponding velocity fields remain unchanged. When using a Richardson  
     678Number dependent eddy viscosity, the mixing of momentum is done through  
     679the vertical diffusion: after a static adjustment, the Richardson Number is zero  
     680and thus the eddy viscosity coefficient is at a maximum. When this convective  
     681adjustment algorithm is used with constant vertical eddy viscosity, spurious  
     682solutions can occur since the vertical momentum diffusion remains small even  
     683after a static adjustment. In that case, we recommend the addition of momentum  
     684mixing in a manner that mimics the mixing in temperature and salinity  
     685\citep{Speich_PhD92, Speich_al_JPO96}. 
    713686 
    714687% ------------------------------------------------------------------------------------------------------------- 
    … …  
    716689% ------------------------------------------------------------------------------------------------------------- 
    717690\subsection   [Enhanced Vertical Diffusion (\np{ln\_zdfevd})] 
    718               {Enhanced Vertical Diffusion (\np{ln\_zdfevd}=true)} 
     691         {Enhanced Vertical Diffusion (\np{ln\_zdfevd}=true)} 
    719692\label{ZDF_evd} 
    720693 
    … …  
    857830% Bottom Friction 
    858831% ================================================================ 
    859 \section  [Bottom and Top Friction (\textit{zdfbfr})]   {Bottom and Top Friction (\mdl{zdfbfr} module)} 
     832\section  [Bottom and top Friction (\textit{zdfbfr})]   {Bottom Friction (\mdl{zdfbfr} module)} 
    860833\label{ZDF_bfr} 
    861834 
    … …  
    865838 
    866839Options to define the top and bottom friction are defined through the  \ngn{nambfr} namelist variables. 
    867 The bottom friction represents the friction generated by the bathymetry.  
    868 The top friction represents the friction generated by the ice shelf/ocean interface.  
    869 As the friction processes at the top and bottom are represented similarly, only the bottom friction is described in detail below.\\ 
    870  
     840The top friction is activated only if the ice shelf cavities are opened (\np{ln\_isfcav}~=~true). 
     841As the friction processes at the top and bottom are the represented similarly, only the bottom friction is described in detail. 
    871842 
    872843Both the surface momentum flux (wind stress) and the bottom momentum  
    … …  
    941912$H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$.  
    942913This is the default value used in \NEMO. It corresponds to a decay time scale  
    943 of 115~days. It can be changed by specifying \np{rn\_bfri1} (namelist parameter). 
     914of 115~days. It can be changed by specifying \np{rn\_bfric1} (namelist parameter). 
    944915 
    945916For the linear friction case the coefficients defined in the general  
    … …  
    951922\end{split} 
    952923\end{equation} 
    953 When \np{nn\_botfr}=1, the value of $r$ used is \np{rn\_bfri1}.  
     924When \np{nn\_botfr}=1, the value of $r$ used is \np{rn\_bfric1}.  
    954925Setting \np{nn\_botfr}=0 is equivalent to setting $r=0$ and leads to a free-slip  
    955926bottom boundary condition. These values are assigned in \mdl{zdfbfr}.  
    … …  
    958929in the \ifile{bfr\_coef} input NetCDF file. The mask values should vary from 0 to 1.  
    959930Locations with a non-zero mask value will have the friction coefficient increased  
    960 by $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri1}. 
     931by $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfric1}. 
    961932 
    962933% ------------------------------------------------------------------------------------------------------------- 
    … …  
    978949$e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{Killworth1992}  
    979950uses $C_D = 1.4\;10^{-3}$ and $e_b =2.5\;\;10^{-3}$m$^2$\;s$^{-2}$.  
    980 The CME choices have been set as default values (\np{rn\_bfri2} and \np{rn\_bfeb2}  
     951The CME choices have been set as default values (\np{rn\_bfric2} and \np{rn\_bfeb2}  
    981952namelist parameters). 
    982953 
    … …  
    993964\end{equation} 
    994965 
    995 The coefficients that control the strength of the non-linear bottom friction are 
    996 initialised as namelist parameters: $C_D$= \np{rn\_bfri2}, and $e_b$ =\np{rn\_bfeb2}. 
    997 Note for applications which treat tides explicitly a low or even zero value of 
    998 \np{rn\_bfeb2} is recommended. From v3.2 onwards a local enhancement of $C_D$ is possible 
    999 via an externally defined 2D mask array (\np{ln\_bfr2d}=true).  This works in the same way 
    1000 as for the linear bottom friction case with non-zero masked locations increased by 
    1001 $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri2}. 
    1002  
    1003 % ------------------------------------------------------------------------------------------------------------- 
    1004 %       Bottom Friction Log-layer 
    1005 % ------------------------------------------------------------------------------------------------------------- 
    1006 \subsection{Log-layer Bottom Friction enhancement (\np{nn\_botfr} = 2, \np{ln\_loglayer} = .true.)} 
    1007 \label{ZDF_bfr_loglayer} 
    1008  
    1009 In the non-linear bottom friction case, the drag coefficient, $C_D$, can be optionally 
    1010 enhanced using a "law of the wall" scaling. If  \np{ln\_loglayer} = .true., $C_D$ is no 
    1011 longer constant but is related to the thickness of the last wet layer in each column by: 
    1012  
    1013 \begin{equation} 
    1014 C_D = \left ( {\kappa \over {\rm log}\left ( 0.5e_{3t}/rn\_bfrz0 \right ) } \right )^2 
    1015 \end{equation} 
    1016  
    1017 \noindent where $\kappa$ is the von-Karman constant and \np{rn\_bfrz0} is a roughness 
    1018 length provided via the namelist. 
    1019  
    1020 For stability, the drag coefficient is bounded such that it is kept greater or equal to 
    1021 the base \np{rn\_bfri2} value and it is not allowed to exceed the value of an additional 
    1022 namelist parameter: \np{rn\_bfri2\_max}, i.e.: 
    1023  
    1024 \begin{equation} 
    1025 rn\_bfri2 \leq C_D \leq rn\_bfri2\_max 
    1026 \end{equation} 
    1027  
    1028 \noindent Note also that a log-layer enhancement can also be applied to the top boundary 
    1029 friction if under ice-shelf cavities are in use (\np{ln\_isfcav}=.true.).  In this case, the 
    1030 relevant namelist parameters are \np{rn\_tfrz0}, \np{rn\_tfri2} 
    1031 and \np{rn\_tfri2\_max}. 
     966The coefficients that control the strength of the non-linear bottom friction are  
     967initialised as namelist parameters: $C_D$= \np{rn\_bfri2}, and $e_b$ =\np{rn\_bfeb2}.  
     968Note for applications which treat tides explicitly a low or even zero value of  
     969\np{rn\_bfeb2} is recommended. From v3.2 onwards a local enhancement of $C_D$  
     970is possible via an externally defined 2D mask array (\np{ln\_bfr2d}=true).  
     971See previous section for details. 
    1032972 
    1033973% ------------------------------------------------------------------------------------------------------------- 
    … …  
    13131253 
    13141254% ================================================================ 
    1315 % Internal wave-driven mixing 
    1316 % ================================================================ 
    1317 \section{Internal wave-driven mixing (\key{zdftmx\_new})} 
    1318 \label{ZDF_tmx_new} 
    1319  
    1320 %--------------------------------------------namzdf_tmx_new------------------------------------------ 
    1321 \namdisplay{namzdf_tmx_new} 
    1322 %-------------------------------------------------------------------------------------------------------------- 
    1323  
    1324 The parameterization of mixing induced by breaking internal waves is a generalization  
    1325 of the approach originally proposed by \citet{St_Laurent_al_GRL02}.  
    1326 A three-dimensional field of internal wave energy dissipation $\epsilon(x,y,z)$ is first constructed,  
    1327 and the resulting diffusivity is obtained as  
    1328 \begin{equation} \label{Eq_Kwave} 
    1329 A^{vT}_{wave} =  R_f \,\frac{ \epsilon }{ \rho \, N^2 } 
    1330 \end{equation} 
    1331 where $R_f$ is the mixing efficiency and $\epsilon$ is a specified three dimensional distribution  
    1332 of the energy available for mixing. If the \np{ln\_mevar} namelist parameter is set to false,  
    1333 the mixing efficiency is taken as constant and equal to 1/6 \citep{Osborn_JPO80}.  
    1334 In the opposite (recommended) case, $R_f$ is instead a function of the turbulence intensity parameter  
    1335 $Re_b = \frac{ \epsilon}{\nu \, N^2}$, with $\nu$ the molecular viscosity of seawater,  
    1336 following the model of \cite{Bouffard_Boegman_DAO2013}  
    1337 and the implementation of \cite{de_lavergne_JPO2016_efficiency}. 
    1338 Note that $A^{vT}_{wave}$ is bounded by $10^{-2}\,m^2/s$, a limit that is often reached when the mixing efficiency is constant. 
    1339  
    1340 In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary  
    1341 as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to true, a recommended choice).  
    1342 This parameterization of differential mixing, due to \cite{Jackson_Rehmann_JPO2014},  
    1343 is implemented as in \cite{de_lavergne_JPO2016_efficiency}. 
    1344  
    1345 The three-dimensional distribution of the energy available for mixing, $\epsilon(i,j,k)$, is constructed  
    1346 from three static maps of column-integrated internal wave energy dissipation, $E_{cri}(i,j)$,  
    1347 $E_{pyc}(i,j)$, and $E_{bot}(i,j)$, combined to three corresponding vertical structures  
    1348 (de Lavergne et al., in prep): 
    1349 \begin{align*} 
    1350 F_{cri}(i,j,k) &\propto e^{-h_{ab} / h_{cri} }\\ 
    1351 F_{pyc}(i,j,k) &\propto N^{n\_p}\\ 
    1352 F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} } 
    1353 \end{align*}  
    1354 In the above formula, $h_{ab}$ denotes the height above bottom,  
    1355 $h_{wkb}$ denotes the WKB-stretched height above bottom, defined by 
    1356 \begin{equation*} 
    1357 h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz'  } \; , 
    1358 \end{equation*} 
    1359 The $n_p$ parameter (given by \np{nn\_zpyc} in \ngn{namzdf\_tmx\_new} namelist)  controls the stratification-dependence of the pycnocline-intensified dissipation.  
    1360 It can take values of 1 (recommended) or 2. 
    1361 Finally, the vertical structures $F_{cri}$ and $F_{bot}$ require the specification of  
    1362 the decay scales $h_{cri}(i,j)$ and $h_{bot}(i,j)$, which are defined by two additional input maps.  
    1363 $h_{cri}$ is related to the large-scale topography of the ocean (etopo2)  
    1364 and $h_{bot}$ is a function of the energy flux $E_{bot}$, the characteristic horizontal scale of  
    1365 the abyssal hill topography \citep{Goff_JGR2010} and the latitude. 
    1366  
    1367 % ================================================================ 
    1368  
    1369  
    1370  

  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Introduction.tex

    r6617 r6625  
    2424release 8.2, described in \citet{Madec1998}. This model has been used for a wide  
    2525range of applications, both regional or global, as a forced ocean model and as a  
    26 model coupled with the sea-ice and/or the atmosphere.   
     26model coupled with the atmosphere. A complete list of references is found on the  
     27\NEMO web site.  
    2728 
    2829This manual is organised in as follows. Chapter~\ref{PE} presents the model basics,  
    2930$i.e.$ the equations and their assumptions, the vertical coordinates used, and the  
    3031subgrid scale physics. This part deals with the continuous equations of the model  
    31 (primitive equations, with temperature, salinity and an equation of seawater).  
     32(primitive equations, with potential temperature, salinity and an equation of state).  
    3233The equations are written in a curvilinear coordinate system, with a choice of vertical  
    3334coordinates ($z$ or $s$, with the rescaled height coordinate formulation \textit{z*}, or   
    … …  
    7879space and time variable coefficient \citet{Treguier1997}. The model has vertical harmonic  
    7980viscosity and diffusion with a space and time variable coefficient, with options to compute  
    80 the coefficients with \citet{Blanke1993}, \citet{Pacanowski_Philander_JPO81},  
     81the coefficients with \citet{Blanke1993}, \citet{Large_al_RG94}, \citet{Pacanowski_Philander_JPO81},  
    8182or \citet{Umlauf_Burchard_JMS03} mixing schemes. 
    8283 \vspace{1cm} 
    8384  
    84 %%gm    To be put somewhere else .... 
    85  
     85  
    8686\noindent CPP keys and namelists are used for inputs to the code.  \newline 
    8787 
    … …  
    112112 \vspace{1cm} 
    113113 
    114 %%gm  end 
    115114 
    116115Model outputs management and specific online diagnostics are described in chapters~\ref{DIA}. 
    … …  
    228227\item a deep re-writting and simplification of the off-line tracer component (OFF\_SRC) ;  
    229228\item the merge of passive and active advection and diffusion modules ; 
    230 \item Use of the Flexible Configuration Manager (FCM) to build configurations, generate the Makefile and produce the executable ; 
     229\item  Use of the Flexible Configuration Manager (FCM) to build configurations, generate the Makefile and produce the executable ; 
    231230\item Linear-tangent and Adjoint component (TAM) added, phased with v3.0 
    232231\end{enumerate} 
    … …  
    250249 
    251250 
    252  \vspace{1cm} 
    253 $\bullet$ The main modifications from NEMO/OPA v3.4 and  v3.6 are :\\ 
    254 \begin{enumerate} 
    255 \item I/O management: NEMO in now interfaced with XIOS, a Input/Output server having a versatile xml user interface, and  
    256 allowing I/O to be performed on dedicated processors thus improving scalability and performance on massively parallel platforms.  
    257 \item ICB module \citep{Marsh_GMD2015}: icebergs as lagrangian floats ;  
    258 \item SAS: Stand Alone Surface module allowing testing of forcing set with bulk formulae, to run sea-ice models without ocean, to run ICB icebergs module alone, and to test AGRIF with sea-ice 
    259 \item ISF : Under ice-selves cavities (parametrisation and/or explicit representation) 
    260 \item Coupled interface for next IPCC requirements (multi category sea-ice, calving and iceberg module) 
    261 \item Ocean and ice allowed to be explicitly coupled through OASIS, using StandAlone Surface module) 
    262 \item On line coarsening of ocean I/O 
    263 \item Major evolution of LIM3 sea-ice model \citep{Rousset_GMD2015} 
    264 \item Open boundaries: completion of BDY/OBC merge : BDY is now the only Open boundary module available 
    265 \item re-visit of the specification of heat/salt(tracers)/mass fluxes ; 
    266 \item levitating or fully embedded sea-ice (for LIM and CICE) ; 
    267 \item a new parameterization of mixing induced by breaking internal waves (de Lavergne et al. in prep.) 
    268 And also: 
    269 \item update of AGRIF package and AGRIF compatibility with LIM2 sea-ice model ; 
    270 \item A new vertical sigma coordinate stretching function \citep{Siddorn_Furner_OM12} ; 
    271 \item Smagorinsky eddy coefficients: the \cite{Griffies_Hallberg_MWR00} Smagorinsky type diffusivity/viscosity for lateral mixing has been introduced ; 
    272 \item Standard Fox Kemper parametrisation 
    273 \item Analytical tropical cyclones taken in account using track and magnitude observations (Vincent et al. JGR 2012a,b) ; 
    274 \item OBS: observation operators improved and now available in Standalone mode ; 
    275 \item Log layer option for bottom friction 
    276 \item Faster split-explicit time stepping ;  
    277 \item Z-tilde ALE coordinates \citep{Leclair_Madec_OM11} ;  
    278 \item implicit bottom friction ; 
    279 \item Runoff improved and SBC with BGC 
    280 \item MPP assessment and optimisation 
    281 \item First steps of wave coupling 
    282  
    283 Features becoming obsolete: LIM2 (replaced by LIM3 monocategory) ; IOIPSL (replaced by XIOS) ;  
    284  
    285 Features that has been removed : LOBSTER (now included in PISCES) ; OBC, replaced by BDY ;    
    286  
    287  
    288  
    289 \end{enumerate} 
    290  
    291  
Note: See TracChangeset for help on using the changeset viewer.