- Timestamp:
- 2015-12-14T09:23:38+01:00 (9 years ago)
- Location:
- branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters
- Files:
-
- 10 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Abstracts_Foreword.tex
r3294 r6040 13 13 be a flexible tool for studying the ocean and its interactions with the others components of 14 14 the earth climate system over a wide range of space and time scales. 15 Prognostic variables are the three-dimensional velocity field, a linear16 or non-linear sea surface height, the temperature and the salinity. In the horizontal direction,17 the model uses a curvilinear orthogonal grid and in the vertical direction, a full or partial step18 $z$-coordinate, or $s$-coordinate, or a mixture of the two. The distribution of variables is a19 three-dimensional Arakawa C-type grid. Various physical choices are available to describe20 ocean physics, including TKE, GLS and KPP vertical physics. Within NEMO, the ocean is21 interfaced with a sea-ice model (LIM v2 and v3), passive tracer and biogeochemical models (TOP)22 and, via the OASIS coupler, with several atmospheric general circulation models. It also23 support two-way grid embedding via the AGRIF software.15 Prognostic variables are the three-dimensional velocity field, a non-linear sea surface height, 16 the \textit{Conservative} Temperature and the \textit{Absolut} 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. 24 24 25 25 % ================================================================ … … 31 31 interactions avec les autres composantes du syst\`{e}me climatique terrestre. 32 32 Les variables pronostiques sont le champ tridimensionnel de vitesse, une hauteur de la mer 33 lin\'{e}aire ou non, la temperature et la salinit\'{e}.33 lin\'{e}aire, la Ttemperature Conservative et la Salinit\'{e} Absolue. 34 34 La distribution des variables se fait sur une grille C d'Arakawa tridimensionnelle utilisant une 35 35 coordonn\'{e}e verticale $z$ \`{a} niveaux entiers ou partiels, ou une coordonn\'{e}e s, ou encore 36 36 une 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 , GLS et KPP. A travers l'infrastructure38 NEMO, l'oc\'{e}an est interfac\'{e} avec des mod\`{e}les de glace de mer , de biog\'{e}ochimie39 et de traceurs passifs, et, via le coupleur OASIS, \`{a} plusieurs mod\`{e}les de circulation40 g\'{e}n\'{e}rale atmosph\'{e}rique. Il supporte \'{e}galement l'embo\^{i}tement interactif de41 maillages via le logiciel AGRIF.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. 42 42 } 43 43 -
branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Annex_C.tex
r3294 r6040 410 410 \end{aligned} } \right. 411 411 \end{equation} 412 where the indices $i_p$ and $ k_p$ take the following value:412 where the indices $i_p$ and $j_p$ take the following value: 413 413 $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$, 414 414 and the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$, defined at $T$-point, are given by: -
branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_DOM.tex
r5120 r6040 1 1 % ================================================================ 2 % Chapter 2 �Space and Time Domain (DOM)2 % Chapter 2 ——— Space and Time Domain (DOM) 3 3 % ================================================================ 4 4 \chapter{Space Domain (DOM) } … … 364 364 For both grids here, the same $w$-point depth has been chosen but in (a) the 365 365 $t$-points are set half way between $w$-points while in (b) they are defined from 366 an analytical function: $z(k)=5\,( i-1/2)^3 - 45\,(i-1/2)^2 + 140\,(i-1/2) - 150$.366 an analytical function: $z(k)=5\,(k-1/2)^3 - 45\,(k-1/2)^2 + 140\,(k-1/2) - 150$. 367 367 Note the resulting difference between the value of the grid-size $\Delta_k$ and 368 368 those of the scale factor $e_k$. } … … 475 475 (d) hybrid $s-z$ coordinate, 476 476 (e) hybrid $s-z$ coordinate with partial step, and 477 (f) same as (e) but with variable volume associated with the non-linear free surface.478 Note that the variable volume option (\key{vvl})can be used with any of the477 (f) same as (e) but in the non-linear free surface (\np{ln\_linssh}=false). 478 Note that the non-linear free surface can be used with any of the 479 479 5 coordinates (a) to (e).} 480 480 \end{center} \end{figure} 481 481 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 482 482 483 The choice of a vertical coordinate, even if it is made through a namelist parameter,483 The choice of a vertical coordinate, even if it is made through \ngn{namzgr} namelist parameters, 484 484 must be done once of all at the beginning of an experiment. It is not intended as an 485 485 option which can be enabled or disabled in the middle of an experiment. Three main … … 488 488 (\np{ln\_zps}~=~true), or generalized, $s$-coordinate (\np{ln\_sco}~=~true). 489 489 Hybridation of the three main coordinates are available: $s-z$ or $s-zps$ coordinate 490 (Fig.~\ref{Fig_z_zps_s_sps}d and \ref{Fig_z_zps_s_sps}e). When using the variable 491 volume option \key{vvl} ($i.e.$ non-linear free surface), the coordinate follow the 492 time-variation of the free surface so that the transformation is time dependent: 493 $z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). This option can be used with full step 494 bathymetry or $s$-coordinate (hybrid and partial step coordinates have not 495 yet 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). 490 (Fig.~\ref{Fig_z_zps_s_sps}d and \ref{Fig_z_zps_s_sps}e). By default a non-linear free surface is used: 491 the coordinate follow the time-variation of the free surface so that the transformation is time dependent: 492 $z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). When a linear free surface is assumed (\np{ln\_linssh}=true), 493 the vertical coordinate are fixed in time, but the seawater can move up and down across the z=0 surface 494 (in other words, the top of the ocean in not a rigid-lid). 495 The last choice in terms of vertical coordinate concerns the presence (or not) in the model domain 496 of ocean cavities beneath ice shelves. Setting \np{ln\_isfcav} to true allows to manage ocean cavities, 497 otherwise they are filled in. 497 498 498 499 Contrary to the horizontal grid, the vertical grid is computed in the code and no … … 519 520 %%% 520 521 521 The arrays describing the grid point depths and vertical scale factors 522 are three dimensional arrays $(i,j,k)$ even in the case of $z$-coordinate with 523 full step bottom topography. In non-linear free surface (\key{vvl}), their knowledge 524 is required at \textit{before}, \textit{now} and \textit{after} time step, while they 525 do not vary in time in linear free surface case. 526 To improve the code readability while providing this flexibility, the vertical coordinate 527 and scale factors are defined as functions of 528 $(i,j,k)$ with "fs" as prefix (examples: \textit{fse3t\_b, fse3t\_n, fse3t\_a,} 529 for the \textit{before}, \textit{now} and \textit{after} scale factors at $t$-point) 530 that can be either three different arrays when \key{vvl} is defined, or a single fixed arrays. 531 These functions are defined in the file \hf{domzgr\_substitute} of the DOM directory. 532 They are used throughout the code, and replaced by the corresponding arrays at 533 the time of pre-processing (CPP capability). 522 Unless a linear free surface is used (\np{ln\_linssh}=false), the arrays describing 523 the grid point depths and vertical scale factors are three set of three dimensional arrays $(i,j,k)$ 524 defined at \textit{before}, \textit{now} and \textit{after} time step. The time at which they are 525 defined is indicated by a suffix:$\_b$, $\_n$, or $\_a$, respectively. They are updated at each model time step 526 using a fixed reference coordinate system which computer names have a $\_0$ suffix. 527 When the linear free surface option is used (\np{ln\_linssh}=true), \textit{before}, \textit{now} 528 and \textit{after} arrays are simply set one for all to their reference counterpart. 529 534 530 535 531 % ------------------------------------------------------------------------------------------------------------- … … 840 836 % z*- or s*-coordinate 841 837 % ------------------------------------------------------------------------------------------------------------- 842 \subsection{$z^*$- or $s^*$-coordinate ( add \key{vvl}) }843 \label{DOM_zgr_ vvl}838 \subsection{$z^*$- or $s^*$-coordinate (\np{ln\_linssh}=false) } 839 \label{DOM_zgr_star} 844 840 845 841 This option is described in the Report by Levier \textit{et al.} (2007), available on -
branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_DYN.tex
r5120 r6040 1 1 % ================================================================ 2 % Chapter �Ocean Dynamics (DYN)2 % Chapter ——— Ocean Dynamics (DYN) 3 3 % ================================================================ 4 4 \chapter{Ocean Dynamics (DYN)} 5 5 \label{DYN} 6 6 \minitoc 7 8 % add a figure for dynvor ens, ene latices9 7 10 8 %\vspace{2.cm} -
branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_LDF.tex
r4147 r6040 1 1 2 2 % ================================================================ 3 % Chapter �Lateral Ocean Physics (LDF)3 % Chapter ——— Lateral Ocean Physics (LDF) 4 4 % ================================================================ 5 5 \chapter{Lateral Ocean Physics (LDF)} … … 15 15 described in \S\ref{PE_zdf} and their discrete formulation in \S\ref{TRA_ldf} 16 16 and \S\ref{DYN_ldf}). In this section we further discuss each lateral physics option. 17 Choosing one lateral physics scheme means for the user defining, (1) the space 18 and time variations of the eddy coefficients ; (2) the direction along which the 19 lateral diffusive fluxes are evaluated (model level, geopotential or isopycnal 20 surfaces); and (3) the type of operator used (harmonic, or biharmonic operators, 21 and for tracers only, eddy induced advection on tracers). These three aspects 22 of the lateral diffusion are set through namelist parameters and CPP keys 23 (see the \textit{\ngn{nam\_traldf}} and \textit{\ngn{nam\_dynldf}} below). Note 24 that this chapter describes the default implementation of iso-neutral 17 Choosing one lateral physics scheme means for the user defining, 18 (1) the type of operator used (laplacian or bilaplacian operators, or no lateral mixing term) ; 19 (2) the direction along which the lateral diffusive fluxes are evaluated (model level, geopotential or isopycnal surfaces) ; and 20 (3) the space and time variations of the eddy coefficients. 21 These three aspects of the lateral diffusion are set through namelist parameters 22 (see the \textit{\ngn{nam\_traldf}} and \textit{\ngn{nam\_dynldf}} below). 23 Note that this chapter describes the standard implementation of iso-neutral 25 24 tracer mixing, and Griffies's implementation, which is used if 26 25 \np{traldf\_grif}=true, is described in Appdx\ref{sec:triad} … … 33 32 34 33 % ================================================================ 34 % Direction of lateral Mixing 35 % ================================================================ 36 \section [Direction of Lateral Mixing (\textit{ldfslp})] 37 {Direction of Lateral Mixing (\mdl{ldfslp})} 38 \label{LDF_slp} 39 40 %%% 41 \gmcomment{ we should emphasize here that the implementation is a rather old one. 42 Better work can be achieved by using \citet{Griffies_al_JPO98, Griffies_Bk04} iso-neutral scheme. } 43 44 A direction for lateral mixing has to be defined when the desired operator does 45 not act along the model levels. This occurs when $(a)$ horizontal mixing is 46 required on tracer or momentum (\np{ln\_traldf\_hor} or \np{ln\_dynldf\_hor}) 47 in $s$- or mixed $s$-$z$- coordinates, and $(b)$ isoneutral mixing is required 48 whatever the vertical coordinate is. This direction of mixing is defined by its 49 slopes in the \textbf{i}- and \textbf{j}-directions at the face of the cell of the 50 quantity to be diffused. For a tracer, this leads to the following four slopes : 51 $r_{1u}$, $r_{1w}$, $r_{2v}$, $r_{2w}$ (see \eqref{Eq_tra_ldf_iso}), while 52 for momentum the slopes are $r_{1t}$, $r_{1uw}$, $r_{2f}$, $r_{2uw}$ for 53 $u$ and $r_{1f}$, $r_{1vw}$, $r_{2t}$, $r_{2vw}$ for $v$. 54 55 %gm% add here afigure of the slope in i-direction 56 57 \subsection{slopes for tracer geopotential mixing in the $s$-coordinate} 58 59 In $s$-coordinates, geopotential mixing ($i.e.$ horizontal mixing) $r_1$ and 60 $r_2$ are the slopes between the geopotential and computational surfaces. 61 Their discrete formulation is found by locally solving \eqref{Eq_tra_ldf_iso} 62 when the diffusive fluxes in the three directions are set to zero and $T$ is 63 assumed to be horizontally uniform, $i.e.$ a linear function of $z_T$, the 64 depth of a $T$-point. 65 %gm { Steven : My version is obviously wrong since I'm left with an arbitrary constant which is the local vertical temperature gradient} 66 67 \begin{equation} \label{Eq_ldfslp_geo} 68 \begin{aligned} 69 r_{1u} &= \frac{e_{3u}}{ \left( e_{1u}\;\overline{\overline{e_{3w}}}^{\,i+1/2,\,k} \right)} 70 \;\delta_{i+1/2}[z_t] 71 &\approx \frac{1}{e_{1u}}\; \delta_{i+1/2}[z_t] 72 \\ 73 r_{2v} &= \frac{e_{3v}}{\left( e_{2v}\;\overline{\overline{e_{3w}}}^{\,j+1/2,\,k} \right)} 74 \;\delta_{j+1/2} [z_t] 75 &\approx \frac{1}{e_{2v}}\; \delta_{j+1/2}[z_t] 76 \\ 77 r_{1w} &= \frac{1}{e_{1w}}\;\overline{\overline{\delta_{i+1/2}[z_t]}}^{\,i,\,k+1/2} 78 &\approx \frac{1}{e_{1w}}\; \delta_{i+1/2}[z_{uw}] 79 \\ 80 r_{2w} &= \frac{1}{e_{2w}}\;\overline{\overline{\delta_{j+1/2}[z_t]}}^{\,j,\,k+1/2} 81 &\approx \frac{1}{e_{2w}}\; \delta_{j+1/2}[z_{vw}] 82 \\ 83 \end{aligned} 84 \end{equation} 85 86 %gm% caution I'm not sure the simplification was a good idea! 87 88 These slopes are computed once in \rou{ldfslp\_init} when \np{ln\_sco}=True, 89 and either \np{ln\_traldf\_hor}=True or \np{ln\_dynldf\_hor}=True. 90 91 \subsection{Slopes for tracer iso-neutral mixing}\label{LDF_slp_iso} 92 In iso-neutral mixing $r_1$ and $r_2$ are the slopes between the iso-neutral 93 and computational surfaces. Their formulation does not depend on the vertical 94 coordinate used. Their discrete formulation is found using the fact that the 95 diffusive fluxes of locally referenced potential density ($i.e.$ $in situ$ density) 96 vanish. So, substituting $T$ by $\rho$ in \eqref{Eq_tra_ldf_iso} and setting the 97 diffusive fluxes in the three directions to zero leads to the following definition for 98 the neutral slopes: 99 100 \begin{equation} \label{Eq_ldfslp_iso} 101 \begin{split} 102 r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac{\delta_{i+1/2}[\rho]} 103 {\overline{\overline{\delta_{k+1/2}[\rho]}}^{\,i+1/2,\,k}} 104 \\ 105 r_{2v} &= \frac{e_{3v}}{e_{2v}}\; \frac{\delta_{j+1/2}\left[\rho \right]} 106 {\overline{\overline{\delta_{k+1/2}[\rho]}}^{\,j+1/2,\,k}} 107 \\ 108 r_{1w} &= \frac{e_{3w}}{e_{1w}}\; 109 \frac{\overline{\overline{\delta_{i+1/2}[\rho]}}^{\,i,\,k+1/2}} 110 {\delta_{k+1/2}[\rho]} 111 \\ 112 r_{2w} &= \frac{e_{3w}}{e_{2w}}\; 113 \frac{\overline{\overline{\delta_{j+1/2}[\rho]}}^{\,j,\,k+1/2}} 114 {\delta_{k+1/2}[\rho]} 115 \\ 116 \end{split} 117 \end{equation} 118 119 %gm% rewrite this as the explanation is not very clear !!! 120 %In practice, \eqref{Eq_ldfslp_iso} is of little help in evaluating the neutral surface slopes. Indeed, for an unsimplified equation of state, the density has a strong dependancy on pressure (here approximated as the depth), therefore applying \eqref{Eq_ldfslp_iso} using the $in situ$ density, $\rho$, computed at T-points leads to a flattening of slopes as the depth increases. This is due to the strong increase of the $in situ$ density with depth. 121 122 %By definition, neutral surfaces are tangent to the local $in situ$ density \citep{McDougall1987}, therefore in \eqref{Eq_ldfslp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters). 123 124 %In the $z$-coordinate, the derivative of the \eqref{Eq_ldfslp_iso} numerator is evaluated at the same depth \nocite{as what?} ($T$-level, which is the same as the $u$- and $v$-levels), so the $in situ$ density can be used for its evaluation. 125 126 As the mixing is performed along neutral surfaces, the gradient of $\rho$ in 127 \eqref{Eq_ldfslp_iso} has to be evaluated at the same local pressure (which, 128 in decibars, is approximated by the depth in meters in the model). Therefore 129 \eqref{Eq_ldfslp_iso} cannot be used as such, but further transformation is 130 needed depending on the vertical coordinate used: 131 132 \begin{description} 133 134 \item[$z$-coordinate with full step : ] in \eqref{Eq_ldfslp_iso} the densities 135 appearing in the $i$ and $j$ derivatives are taken at the same depth, thus 136 the $in situ$ density can be used. This is not the case for the vertical 137 derivatives: $\delta_{k+1/2}[\rho]$ is replaced by $-\rho N^2/g$, where $N^2$ 138 is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following 139 \citet{McDougall1987} (see \S\ref{TRA_bn2}). 140 141 \item[$z$-coordinate with partial step : ] this case is identical to the full step 142 case except that at partial step level, the \emph{horizontal} density gradient 143 is evaluated as described in \S\ref{TRA_zpshde}. 144 145 \item[$s$- or hybrid $s$-$z$- coordinate : ] in the current release of \NEMO, 146 iso-neutral mixing is only employed for $s$-coordinates if the 147 Griffies scheme is used (\np{traldf\_grif}=true; see Appdx \ref{sec:triad}). 148 In other words, iso-neutral mixing will only be accurately represented with a 149 linear equation of state (\np{nn\_eos}=1 or 2). In the case of a "true" equation 150 of state, the evaluation of $i$ and $j$ derivatives in \eqref{Eq_ldfslp_iso} 151 will include a pressure dependent part, leading to the wrong evaluation of 152 the neutral slopes. 153 154 %gm% 155 Note: The solution for $s$-coordinate passes trough the use of different 156 (and better) expression for the constraint on iso-neutral fluxes. Following 157 \citet{Griffies_Bk04}, instead of specifying directly that there is a zero neutral 158 diffusive flux of locally referenced potential density, we stay in the $T$-$S$ 159 plane and consider the balance between the neutral direction diffusive fluxes 160 of potential temperature and salinity: 161 \begin{equation} 162 \alpha \ \textbf{F}(T) = \beta \ \textbf{F}(S) 163 \end{equation} 164 %gm{ where vector F is ....} 165 166 This constraint leads to the following definition for the slopes: 167 168 \begin{equation} \label{Eq_ldfslp_iso2} 169 \begin{split} 170 r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac 171 {\alpha_u \;\delta_{i+1/2}[T] - \beta_u \;\delta_{i+1/2}[S]} 172 {\alpha_u \;\overline{\overline{\delta_{k+1/2}[T]}}^{\,i+1/2,\,k} 173 -\beta_u \;\overline{\overline{\delta_{k+1/2}[S]}}^{\,i+1/2,\,k} } 174 \\ 175 r_{2v} &= \frac{e_{3v}}{e_{2v}}\; \frac 176 {\alpha_v \;\delta_{j+1/2}[T] - \beta_v \;\delta_{j+1/2}[S]} 177 {\alpha_v \;\overline{\overline{\delta_{k+1/2}[T]}}^{\,j+1/2,\,k} 178 -\beta_v \;\overline{\overline{\delta_{k+1/2}[S]}}^{\,j+1/2,\,k} } 179 \\ 180 r_{1w} &= \frac{e_{3w}}{e_{1w}}\; \frac 181 {\alpha_w \;\overline{\overline{\delta_{i+1/2}[T]}}^{\,i,\,k+1/2} 182 -\beta_w \;\overline{\overline{\delta_{i+1/2}[S]}}^{\,i,\,k+1/2} } 183 {\alpha_w \;\delta_{k+1/2}[T] - \beta_w \;\delta_{k+1/2}[S]} 184 \\ 185 r_{2w} &= \frac{e_{3w}}{e_{2w}}\; \frac 186 {\alpha_w \;\overline{\overline{\delta_{j+1/2}[T]}}^{\,j,\,k+1/2} 187 -\beta_w \;\overline{\overline{\delta_{j+1/2}[S]}}^{\,j,\,k+1/2} } 188 {\alpha_w \;\delta_{k+1/2}[T] - \beta_w \;\delta_{k+1/2}[S]} 189 \\ 190 \end{split} 191 \end{equation} 192 where $\alpha$ and $\beta$, the thermal expansion and saline contraction 193 coefficients introduced in \S\ref{TRA_bn2}, have to be evaluated at the three 194 velocity points. In order to save computation time, they should be approximated 195 by the mean of their values at $T$-points (for example in the case of $\alpha$: 196 $\alpha_u=\overline{\alpha_T}^{i+1/2}$, $\alpha_v=\overline{\alpha_T}^{j+1/2}$ 197 and $\alpha_w=\overline{\alpha_T}^{k+1/2}$). 198 199 Note that such a formulation could be also used in the $z$-coordinate and 200 $z$-coordinate with partial steps cases. 201 202 \end{description} 203 204 This implementation is a rather old one. It is similar to the one 205 proposed by Cox [1987], except for the background horizontal 206 diffusion. Indeed, the Cox implementation of isopycnal diffusion in 207 GFDL-type models requires a minimum background horizontal diffusion 208 for numerical stability reasons. To overcome this problem, several 209 techniques have been proposed in which the numerical schemes of the 210 ocean model are modified \citep{Weaver_Eby_JPO97, 211 Griffies_al_JPO98}. Griffies's scheme is now available in \NEMO if 212 \np{traldf\_grif\_iso} is set true; see Appdx \ref{sec:triad}. Here, 213 another strategy is presented \citep{Lazar_PhD97}: a local 214 filtering of the iso-neutral slopes (made on 9 grid-points) prevents 215 the development of grid point noise generated by the iso-neutral 216 diffusion operator (Fig.~\ref{Fig_LDF_ZDF1}). This allows an 217 iso-neutral diffusion scheme without additional background horizontal 218 mixing. This technique can be viewed as a diffusion operator that acts 219 along large-scale (2~$\Delta$x) \gmcomment{2deltax doesnt seem very 220 large scale} iso-neutral surfaces. The diapycnal diffusion required 221 for numerical stability is thus minimized and its net effect on the 222 flow is quite small when compared to the effect of an horizontal 223 background mixing. 224 225 Nevertheless, this iso-neutral operator does not ensure that variance cannot increase, 226 contrary to the \citet{Griffies_al_JPO98} operator which has that property. 227 228 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 229 \begin{figure}[!ht] \begin{center} 230 \includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_LDF_ZDF1.pdf} 231 \caption { \label{Fig_LDF_ZDF1} 232 averaging procedure for isopycnal slope computation.} 233 \end{center} \end{figure} 234 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 235 236 %There are three additional questions about the slope calculation. 237 %First the expression for the rotation tensor has been obtain assuming the "small slope" approximation, so a bound has to be imposed on slopes. 238 %Second, numerical stability issues also require a bound on slopes. 239 %Third, the question of boundary condition specified on slopes... 240 241 %from griffies: chapter 13.1.... 242 243 244 245 % In addition and also for numerical stability reasons \citep{Cox1987, Griffies_Bk04}, 246 % the slopes are bounded by $1/100$ everywhere. This limit is decreasing linearly 247 % to zero fom $70$ meters depth and the surface (the fact that the eddies "feel" the 248 % surface motivates this flattening of isopycnals near the surface). 249 250 For numerical stability reasons \citep{Cox1987, Griffies_Bk04}, the slopes must also 251 be bounded by $1/100$ everywhere. This constraint is applied in a piecewise linear 252 fashion, increasing from zero at the surface to $1/100$ at $70$ metres and thereafter 253 decreasing to zero at the bottom of the ocean. (the fact that the eddies "feel" the 254 surface motivates this flattening of isopycnals near the surface). 255 256 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 257 \begin{figure}[!ht] \begin{center} 258 \includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_eiv_slp.pdf} 259 \caption { \label{Fig_eiv_slp} 260 Vertical profile of the slope used for lateral mixing in the mixed layer : 261 \textit{(a)} in the real ocean the slope is the iso-neutral slope in the ocean interior, 262 which has to be adjusted at the surface boundary (i.e. it must tend to zero at the 263 surface since there is no mixing across the air-sea interface: wall boundary 264 condition). Nevertheless, the profile between the surface zero value and the interior 265 iso-neutral one is unknown, and especially the value at the base of the mixed layer ; 266 \textit{(b)} profile of slope using a linear tapering of the slope near the surface and 267 imposing a maximum slope of 1/100 ; \textit{(c)} profile of slope actually used in 268 \NEMO: a linear decrease of the slope from zero at the surface to its ocean interior 269 value computed just below the mixed layer. Note the huge change in the slope at the 270 base of the mixed layer between \textit{(b)} and \textit{(c)}.} 271 \end{center} \end{figure} 272 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 273 274 \colorbox{yellow}{add here a discussion about the flattening of the slopes, vs tapering the coefficient.} 275 276 \subsection{slopes for momentum iso-neutral mixing} 277 278 The iso-neutral diffusion operator on momentum is the same as the one used on 279 tracers but applied to each component of the velocity separately (see 280 \eqref{Eq_dyn_ldf_iso} in section~\ref{DYN_ldf_iso}). The slopes between the 281 surface along which the diffusion operator acts and the surface of computation 282 ($z$- or $s$-surfaces) are defined at $T$-, $f$-, and \textit{uw}- points for the 283 $u$-component, and $T$-, $f$- and \textit{vw}- points for the $v$-component. 284 They are computed from the slopes used for tracer diffusion, $i.e.$ 285 \eqref{Eq_ldfslp_geo} and \eqref{Eq_ldfslp_iso} : 286 287 \begin{equation} \label{Eq_ldfslp_dyn} 288 \begin{aligned} 289 &r_{1t}\ \ = \overline{r_{1u}}^{\,i} &&& r_{1f}\ \ &= \overline{r_{1u}}^{\,i+1/2} \\ 290 &r_{2f} \ \ = \overline{r_{2v}}^{\,j+1/2} &&& r_{2t}\ &= \overline{r_{2v}}^{\,j} \\ 291 &r_{1uw} = \overline{r_{1w}}^{\,i+1/2} &&\ \ \text{and} \ \ & r_{1vw}&= \overline{r_{1w}}^{\,j+1/2} \\ 292 &r_{2uw}= \overline{r_{2w}}^{\,j+1/2} &&& r_{2vw}&= \overline{r_{2w}}^{\,j+1/2}\\ 293 \end{aligned} 294 \end{equation} 295 296 The major issue remaining is in the specification of the boundary conditions. 297 The same boundary conditions are chosen as those used for lateral 298 diffusion along model level surfaces, i.e. using the shear computed along 299 the model levels and with no additional friction at the ocean bottom (see 300 {\S\ref{LBC_coast}). 301 302 303 % ================================================================ 304 % Lateral Mixing Operator 305 % ================================================================ 306 \section [Lateral Mixing Operators (\textit{ldftra}, \textit{ldfdyn})] 307 {Lateral Mixing Operators (\mdl{traldf}, \mdl{traldf}) } 308 \label{LDF_op} 309 310 311 312 % ================================================================ 35 313 % Lateral Mixing Coefficients 36 314 % ================================================================ … … 38 316 {Lateral Mixing Coefficient (\mdl{ldftra}, \mdl{ldfdyn}) } 39 317 \label{LDF_coef} 40 41 318 42 319 Introducing a space variation in the lateral eddy mixing coefficients changes … … 165 442 166 443 % ================================================================ 167 % Direction of lateral Mixing168 % ================================================================169 \section [Direction of Lateral Mixing (\textit{ldfslp})]170 {Direction of Lateral Mixing (\mdl{ldfslp})}171 \label{LDF_slp}172 173 %%%174 \gmcomment{ we should emphasize here that the implementation is a rather old one.175 Better work can be achieved by using \citet{Griffies_al_JPO98, Griffies_Bk04} iso-neutral scheme. }176 177 A direction for lateral mixing has to be defined when the desired operator does178 not act along the model levels. This occurs when $(a)$ horizontal mixing is179 required on tracer or momentum (\np{ln\_traldf\_hor} or \np{ln\_dynldf\_hor})180 in $s$- or mixed $s$-$z$- coordinates, and $(b)$ isoneutral mixing is required181 whatever the vertical coordinate is. This direction of mixing is defined by its182 slopes in the \textbf{i}- and \textbf{j}-directions at the face of the cell of the183 quantity to be diffused. For a tracer, this leads to the following four slopes :184 $r_{1u}$, $r_{1w}$, $r_{2v}$, $r_{2w}$ (see \eqref{Eq_tra_ldf_iso}), while185 for momentum the slopes are $r_{1t}$, $r_{1uw}$, $r_{2f}$, $r_{2uw}$ for186 $u$ and $r_{1f}$, $r_{1vw}$, $r_{2t}$, $r_{2vw}$ for $v$.187 188 %gm% add here afigure of the slope in i-direction189 190 \subsection{slopes for tracer geopotential mixing in the $s$-coordinate}191 192 In $s$-coordinates, geopotential mixing ($i.e.$ horizontal mixing) $r_1$ and193 $r_2$ are the slopes between the geopotential and computational surfaces.194 Their discrete formulation is found by locally solving \eqref{Eq_tra_ldf_iso}195 when the diffusive fluxes in the three directions are set to zero and $T$ is196 assumed to be horizontally uniform, $i.e.$ a linear function of $z_T$, the197 depth of a $T$-point.198 %gm { Steven : My version is obviously wrong since I'm left with an arbitrary constant which is the local vertical temperature gradient}199 200 \begin{equation} \label{Eq_ldfslp_geo}201 \begin{aligned}202 r_{1u} &= \frac{e_{3u}}{ \left( e_{1u}\;\overline{\overline{e_{3w}}}^{\,i+1/2,\,k} \right)}203 \;\delta_{i+1/2}[z_t]204 &\approx \frac{1}{e_{1u}}\; \delta_{i+1/2}[z_t]205 \\206 r_{2v} &= \frac{e_{3v}}{\left( e_{2v}\;\overline{\overline{e_{3w}}}^{\,j+1/2,\,k} \right)}207 \;\delta_{j+1/2} [z_t]208 &\approx \frac{1}{e_{2v}}\; \delta_{j+1/2}[z_t]209 \\210 r_{1w} &= \frac{1}{e_{1w}}\;\overline{\overline{\delta_{i+1/2}[z_t]}}^{\,i,\,k+1/2}211 &\approx \frac{1}{e_{1w}}\; \delta_{i+1/2}[z_{uw}]212 \\213 r_{2w} &= \frac{1}{e_{2w}}\;\overline{\overline{\delta_{j+1/2}[z_t]}}^{\,j,\,k+1/2}214 &\approx \frac{1}{e_{2w}}\; \delta_{j+1/2}[z_{vw}]215 \\216 \end{aligned}217 \end{equation}218 219 %gm% caution I'm not sure the simplification was a good idea!220 221 These slopes are computed once in \rou{ldfslp\_init} when \np{ln\_sco}=True,222 and either \np{ln\_traldf\_hor}=True or \np{ln\_dynldf\_hor}=True.223 224 \subsection{Slopes for tracer iso-neutral mixing}\label{LDF_slp_iso}225 In iso-neutral mixing $r_1$ and $r_2$ are the slopes between the iso-neutral226 and computational surfaces. Their formulation does not depend on the vertical227 coordinate used. Their discrete formulation is found using the fact that the228 diffusive fluxes of locally referenced potential density ($i.e.$ $in situ$ density)229 vanish. So, substituting $T$ by $\rho$ in \eqref{Eq_tra_ldf_iso} and setting the230 diffusive fluxes in the three directions to zero leads to the following definition for231 the neutral slopes:232 233 \begin{equation} \label{Eq_ldfslp_iso}234 \begin{split}235 r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac{\delta_{i+1/2}[\rho]}236 {\overline{\overline{\delta_{k+1/2}[\rho]}}^{\,i+1/2,\,k}}237 \\238 r_{2v} &= \frac{e_{3v}}{e_{2v}}\; \frac{\delta_{j+1/2}\left[\rho \right]}239 {\overline{\overline{\delta_{k+1/2}[\rho]}}^{\,j+1/2,\,k}}240 \\241 r_{1w} &= \frac{e_{3w}}{e_{1w}}\;242 \frac{\overline{\overline{\delta_{i+1/2}[\rho]}}^{\,i,\,k+1/2}}243 {\delta_{k+1/2}[\rho]}244 \\245 r_{2w} &= \frac{e_{3w}}{e_{2w}}\;246 \frac{\overline{\overline{\delta_{j+1/2}[\rho]}}^{\,j,\,k+1/2}}247 {\delta_{k+1/2}[\rho]}248 \\249 \end{split}250 \end{equation}251 252 %gm% rewrite this as the explanation is not very clear !!!253 %In practice, \eqref{Eq_ldfslp_iso} is of little help in evaluating the neutral surface slopes. Indeed, for an unsimplified equation of state, the density has a strong dependancy on pressure (here approximated as the depth), therefore applying \eqref{Eq_ldfslp_iso} using the $in situ$ density, $\rho$, computed at T-points leads to a flattening of slopes as the depth increases. This is due to the strong increase of the $in situ$ density with depth.254 255 %By definition, neutral surfaces are tangent to the local $in situ$ density \citep{McDougall1987}, therefore in \eqref{Eq_ldfslp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters).256 257 %In the $z$-coordinate, the derivative of the \eqref{Eq_ldfslp_iso} numerator is evaluated at the same depth \nocite{as what?} ($T$-level, which is the same as the $u$- and $v$-levels), so the $in situ$ density can be used for its evaluation.258 259 As the mixing is performed along neutral surfaces, the gradient of $\rho$ in260 \eqref{Eq_ldfslp_iso} has to be evaluated at the same local pressure (which,261 in decibars, is approximated by the depth in meters in the model). Therefore262 \eqref{Eq_ldfslp_iso} cannot be used as such, but further transformation is263 needed depending on the vertical coordinate used:264 265 \begin{description}266 267 \item[$z$-coordinate with full step : ] in \eqref{Eq_ldfslp_iso} the densities268 appearing in the $i$ and $j$ derivatives are taken at the same depth, thus269 the $in situ$ density can be used. This is not the case for the vertical270 derivatives: $\delta_{k+1/2}[\rho]$ is replaced by $-\rho N^2/g$, where $N^2$271 is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following272 \citet{McDougall1987} (see \S\ref{TRA_bn2}).273 274 \item[$z$-coordinate with partial step : ] this case is identical to the full step275 case except that at partial step level, the \emph{horizontal} density gradient276 is evaluated as described in \S\ref{TRA_zpshde}.277 278 \item[$s$- or hybrid $s$-$z$- coordinate : ] in the current release of \NEMO,279 iso-neutral mixing is only employed for $s$-coordinates if the280 Griffies scheme is used (\np{traldf\_grif}=true; see Appdx \ref{sec:triad}).281 In other words, iso-neutral mixing will only be accurately represented with a282 linear equation of state (\np{nn\_eos}=1 or 2). In the case of a "true" equation283 of state, the evaluation of $i$ and $j$ derivatives in \eqref{Eq_ldfslp_iso}284 will include a pressure dependent part, leading to the wrong evaluation of285 the neutral slopes.286 287 %gm%288 Note: The solution for $s$-coordinate passes trough the use of different289 (and better) expression for the constraint on iso-neutral fluxes. Following290 \citet{Griffies_Bk04}, instead of specifying directly that there is a zero neutral291 diffusive flux of locally referenced potential density, we stay in the $T$-$S$292 plane and consider the balance between the neutral direction diffusive fluxes293 of potential temperature and salinity:294 \begin{equation}295 \alpha \ \textbf{F}(T) = \beta \ \textbf{F}(S)296 \end{equation}297 %gm{ where vector F is ....}298 299 This constraint leads to the following definition for the slopes:300 301 \begin{equation} \label{Eq_ldfslp_iso2}302 \begin{split}303 r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac304 {\alpha_u \;\delta_{i+1/2}[T] - \beta_u \;\delta_{i+1/2}[S]}305 {\alpha_u \;\overline{\overline{\delta_{k+1/2}[T]}}^{\,i+1/2,\,k}306 -\beta_u \;\overline{\overline{\delta_{k+1/2}[S]}}^{\,i+1/2,\,k} }307 \\308 r_{2v} &= \frac{e_{3v}}{e_{2v}}\; \frac309 {\alpha_v \;\delta_{j+1/2}[T] - \beta_v \;\delta_{j+1/2}[S]}310 {\alpha_v \;\overline{\overline{\delta_{k+1/2}[T]}}^{\,j+1/2,\,k}311 -\beta_v \;\overline{\overline{\delta_{k+1/2}[S]}}^{\,j+1/2,\,k} }312 \\313 r_{1w} &= \frac{e_{3w}}{e_{1w}}\; \frac314 {\alpha_w \;\overline{\overline{\delta_{i+1/2}[T]}}^{\,i,\,k+1/2}315 -\beta_w \;\overline{\overline{\delta_{i+1/2}[S]}}^{\,i,\,k+1/2} }316 {\alpha_w \;\delta_{k+1/2}[T] - \beta_w \;\delta_{k+1/2}[S]}317 \\318 r_{2w} &= \frac{e_{3w}}{e_{2w}}\; \frac319 {\alpha_w \;\overline{\overline{\delta_{j+1/2}[T]}}^{\,j,\,k+1/2}320 -\beta_w \;\overline{\overline{\delta_{j+1/2}[S]}}^{\,j,\,k+1/2} }321 {\alpha_w \;\delta_{k+1/2}[T] - \beta_w \;\delta_{k+1/2}[S]}322 \\323 \end{split}324 \end{equation}325 where $\alpha$ and $\beta$, the thermal expansion and saline contraction326 coefficients introduced in \S\ref{TRA_bn2}, have to be evaluated at the three327 velocity points. In order to save computation time, they should be approximated328 by the mean of their values at $T$-points (for example in the case of $\alpha$:329 $\alpha_u=\overline{\alpha_T}^{i+1/2}$, $\alpha_v=\overline{\alpha_T}^{j+1/2}$330 and $\alpha_w=\overline{\alpha_T}^{k+1/2}$).331 332 Note that such a formulation could be also used in the $z$-coordinate and333 $z$-coordinate with partial steps cases.334 335 \end{description}336 337 This implementation is a rather old one. It is similar to the one338 proposed by Cox [1987], except for the background horizontal339 diffusion. Indeed, the Cox implementation of isopycnal diffusion in340 GFDL-type models requires a minimum background horizontal diffusion341 for numerical stability reasons. To overcome this problem, several342 techniques have been proposed in which the numerical schemes of the343 ocean model are modified \citep{Weaver_Eby_JPO97,344 Griffies_al_JPO98}. Griffies's scheme is now available in \NEMO if345 \np{traldf\_grif\_iso} is set true; see Appdx \ref{sec:triad}. Here,346 another strategy is presented \citep{Lazar_PhD97}: a local347 filtering of the iso-neutral slopes (made on 9 grid-points) prevents348 the development of grid point noise generated by the iso-neutral349 diffusion operator (Fig.~\ref{Fig_LDF_ZDF1}). This allows an350 iso-neutral diffusion scheme without additional background horizontal351 mixing. This technique can be viewed as a diffusion operator that acts352 along large-scale (2~$\Delta$x) \gmcomment{2deltax doesnt seem very353 large scale} iso-neutral surfaces. The diapycnal diffusion required354 for numerical stability is thus minimized and its net effect on the355 flow is quite small when compared to the effect of an horizontal356 background mixing.357 358 Nevertheless, this iso-neutral operator does not ensure that variance cannot increase,359 contrary to the \citet{Griffies_al_JPO98} operator which has that property.360 361 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>362 \begin{figure}[!ht] \begin{center}363 \includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_LDF_ZDF1.pdf}364 \caption { \label{Fig_LDF_ZDF1}365 averaging procedure for isopycnal slope computation.}366 \end{center} \end{figure}367 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>368 369 %There are three additional questions about the slope calculation.370 %First the expression for the rotation tensor has been obtain assuming the "small slope" approximation, so a bound has to be imposed on slopes.371 %Second, numerical stability issues also require a bound on slopes.372 %Third, the question of boundary condition specified on slopes...373 374 %from griffies: chapter 13.1....375 376 377 378 % In addition and also for numerical stability reasons \citep{Cox1987, Griffies_Bk04},379 % the slopes are bounded by $1/100$ everywhere. This limit is decreasing linearly380 % to zero fom $70$ meters depth and the surface (the fact that the eddies "feel" the381 % surface motivates this flattening of isopycnals near the surface).382 383 For numerical stability reasons \citep{Cox1987, Griffies_Bk04}, the slopes must also384 be bounded by $1/100$ everywhere. This constraint is applied in a piecewise linear385 fashion, increasing from zero at the surface to $1/100$ at $70$ metres and thereafter386 decreasing to zero at the bottom of the ocean. (the fact that the eddies "feel" the387 surface motivates this flattening of isopycnals near the surface).388 389 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>390 \begin{figure}[!ht] \begin{center}391 \includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_eiv_slp.pdf}392 \caption { \label{Fig_eiv_slp}393 Vertical profile of the slope used for lateral mixing in the mixed layer :394 \textit{(a)} in the real ocean the slope is the iso-neutral slope in the ocean interior,395 which has to be adjusted at the surface boundary (i.e. it must tend to zero at the396 surface since there is no mixing across the air-sea interface: wall boundary397 condition). Nevertheless, the profile between the surface zero value and the interior398 iso-neutral one is unknown, and especially the value at the base of the mixed layer ;399 \textit{(b)} profile of slope using a linear tapering of the slope near the surface and400 imposing a maximum slope of 1/100 ; \textit{(c)} profile of slope actually used in401 \NEMO: a linear decrease of the slope from zero at the surface to its ocean interior402 value computed just below the mixed layer. Note the huge change in the slope at the403 base of the mixed layer between \textit{(b)} and \textit{(c)}.}404 \end{center} \end{figure}405 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>406 407 \colorbox{yellow}{add here a discussion about the flattening of the slopes, vs tapering the coefficient.}408 409 \subsection{slopes for momentum iso-neutral mixing}410 411 The iso-neutral diffusion operator on momentum is the same as the one used on412 tracers but applied to each component of the velocity separately (see413 \eqref{Eq_dyn_ldf_iso} in section~\ref{DYN_ldf_iso}). The slopes between the414 surface along which the diffusion operator acts and the surface of computation415 ($z$- or $s$-surfaces) are defined at $T$-, $f$-, and \textit{uw}- points for the416 $u$-component, and $T$-, $f$- and \textit{vw}- points for the $v$-component.417 They are computed from the slopes used for tracer diffusion, $i.e.$418 \eqref{Eq_ldfslp_geo} and \eqref{Eq_ldfslp_iso} :419 420 \begin{equation} \label{Eq_ldfslp_dyn}421 \begin{aligned}422 &r_{1t}\ \ = \overline{r_{1u}}^{\,i} &&& r_{1f}\ \ &= \overline{r_{1u}}^{\,i+1/2} \\423 &r_{2f} \ \ = \overline{r_{2v}}^{\,j+1/2} &&& r_{2t}\ &= \overline{r_{2v}}^{\,j} \\424 &r_{1uw} = \overline{r_{1w}}^{\,i+1/2} &&\ \ \text{and} \ \ & r_{1vw}&= \overline{r_{1w}}^{\,j+1/2} \\425 &r_{2uw}= \overline{r_{2w}}^{\,j+1/2} &&& r_{2vw}&= \overline{r_{2w}}^{\,j+1/2}\\426 \end{aligned}427 \end{equation}428 429 The major issue remaining is in the specification of the boundary conditions.430 The same boundary conditions are chosen as those used for lateral431 diffusion along model level surfaces, i.e. using the shear computed along432 the model levels and with no additional friction at the ocean bottom (see433 {\S\ref{LBC_coast}).434 435 436 % ================================================================437 444 % Eddy Induced Mixing 438 445 % ================================================================ -
branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_Model_Basics.tex
r3294 r6040 247 247 sufficient to solve a linearized version of (\ref{Eq_PE_ssh}), which still allows 248 248 to 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. 249 250 250 251 The filtering of EGWs in models with a free surface is usually a matter of discretisation 251 of the temporal derivatives, using the time splitting method \citep{Killworth_al_JPO91, Zhang_Endoh_JGR92} 252 or the implicit scheme \citep{Dukowicz1994}. In \NEMO, we use a slightly different approach 253 developed by \citet{Roullet_Madec_JGR00}: the damping of EGWs is ensured by introducing an 254 additional 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} 260 where $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}$ 262 represents the collected contributions of the Coriolis, hydrostatic pressure gradient, 263 non-linear and viscous terms in \eqref{Eq_PE_dyn}. 264 265 The new force can be interpreted as a diffusion of vertically integrated volume flux divergence. 266 The time evolution of $D$ is thus governed by a balance of two terms, $-g$ \textbf{A} $\eta$ 267 and $g \, T_c \,$ \textbf{A} $D$, associated with a propagative regime and a diffusive regime 268 in 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 271 can 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 273 has to be computed implicitly. This is not surprising since the use of a large time step has a 274 necessarily numerical cost. Two gains arise in comparison with the previous formulations. 275 Firstly, the damping of EGWs can be quantified through the magnitude of the additional term. 276 Secondly, the numerical scheme does not need any tuning. Numerical stability is ensured as 277 soon as $T_c > \rdt$. 278 279 When the variations of free surface elevation are small compared to the thickness of the first 280 model layer, the free surface equation (\ref{Eq_PE_ssh}) can be linearized. As emphasized 281 by \citet{Roullet_Madec_JGR00} the linearization of (\ref{Eq_PE_ssh}) has consequences on the 282 conservation of salt in the model. With the nonlinear free surface equation, the time evolution 283 of 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} 288 where $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 290 integral over the free surface) vanishes when the nonlinear equation (\ref{Eq_PE_ssh}) 291 is satisfied, so that the salt is perfectly conserved. When the free surface equation is 292 linearized, \citet{Roullet_Madec_JGR00} show that the total salt content integrated in the fixed 293 volume $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 299 The 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 301 the globally averaged free surface is drifting. An increase in sea surface height \textit{$\eta $} 302 results in a decrease of the salinity in the fixed volume $D$. Even in that case though, 303 the 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 311 Although the total salt content is not exactly conserved with the linearized free surface, 312 its variations are driven by correlations of the time variation of surface salinity with the 313 sea surface height, which is a negligible term. This situation contrasts with the case of 314 the rigid lid approximation in which case freshwater forcing is represented by a virtual 315 salt 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 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 ligne 320 260 321 261 % ================================================================ … … 655 595 the surface pressure, is given by: 656 596 \begin{equation} \label{Eq_PE_spg} 657 p_s = \left\{ \begin{split} 658 \rho \,g \,\eta & \qquad \qquad \; \qquad \text{ standard free surface} \\ 659 \rho \,g \,\eta &+ \rho_o \,\mu \,\frac{\partial \eta }{\partial t} \qquad \text{ filtered free surface} 660 \end{split} 661 \right. 597 p_s = \rho \,g \,\eta 662 598 \end{equation} 663 599 with $\eta$ is solution of \eqref{Eq_PE_ssh} … … 773 709 \end{equation} 774 710 775 The equations solved by the ocean model \eqref{Eq_PE} in $s-$coordinate can be written as follows :711 The equations solved by the ocean model \eqref{Eq_PE} in $s-$coordinate can be written as follows (see Appendix~\ref{Apdx_A_momentum}): 776 712 777 713 \vspace{0.5cm} 778 * momentum equation:714 $\bullet$ Vector invariant form of the momentum equation : 779 715 \begin{multline} \label{Eq_PE_sco_u} 780 \frac{ 1}{e_3} \frac{\partial \left( e_3\,u \right)}{\partial t}=716 \frac{\partial u }{\partial t}= 781 717 + \left( {\zeta +f} \right)\,v 782 718 - \frac{1}{2\,e_1} \frac{\partial}{\partial i} \left( u^2+v^2 \right) … … 787 723 \end{multline} 788 724 \begin{multline} \label{Eq_PE_sco_v} 789 \frac{ 1}{e_3} \frac{\partial \left( e_3\,v \right)}{\partial t}=725 \frac{\partial v }{\partial t}= 790 726 - \left( {\zeta +f} \right)\,u 791 727 - \frac{1}{2\,e_2 }\frac{\partial }{\partial j}\left( u^2+v^2 \right) … … 795 731 + D_v^{\vect{U}} + F_v^{\vect{U}} \quad 796 732 \end{multline} 733 734 \vspace{0.5cm} 735 $\bullet$ Vector invariant form of the momentum equation : 736 \begin{multline} \label{Eq_PE_sco_u} 737 \frac{1}{e_3} \frac{\partial \left( e_3\,u \right) }{\partial t}= 738 + \left( { f + \frac{1}{e_1 \; e_2 } 739 \left( v \frac{\partial e_2}{\partial i} 740 -u \frac{\partial e_1}{\partial j} \right)} \right) \, v \\ 741 - \frac{1}{e_1 \; e_2 \; e_3 } \left( 742 \frac{\partial \left( {e_2 \, e_3 \, u\,u} \right)}{\partial i} 743 + \frac{\partial \left( {e_1 \, e_3 \, v\,u} \right)}{\partial j} \right) 744 - \frac{1}{e_3 }\frac{\partial \left( { \omega\,u} \right)}{\partial k} \\ 745 - \frac{1}{e_1} \frac{\partial}{\partial i} \left( \frac{p_s + p_h}{\rho _o} \right) 746 + g\frac{\rho }{\rho _o}\sigma _1 747 + D_u^{\vect{U}} + F_u^{\vect{U}} \quad 748 \end{multline} 749 \begin{multline} \label{Eq_PE_sco_v} 750 \frac{1}{e_3} \frac{\partial \left( e_3\,v \right) }{\partial t}= 751 - \left( { f + \frac{1}{e_1 \; e_2} 752 \left( v \frac{\partial e_2}{\partial i} 753 -u \frac{\partial e_1}{\partial j} \right)} \right) \, u \\ 754 - \frac{1}{e_1 \; e_2 \; e_3 } \left( 755 \frac{\partial \left( {e_2 \; e_3 \,u\,v} \right)}{\partial i} 756 + \frac{\partial \left( {e_1 \; e_3 \,v\,v} \right)}{\partial j} \right) 757 - \frac{1}{e_3 } \frac{\partial \left( { \omega\,v} \right)}{\partial k} \\ 758 - \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho _o} \right) 759 + g\frac{\rho }{\rho _o }\sigma _2 760 + D_v^{\vect{U}} + F_v^{\vect{U}} \quad 761 \end{multline} 762 797 763 where the relative vorticity, \textit{$\zeta $}, the surface pressure gradient, and the hydrostatic 798 764 pressure have the same expressions as in $z$-coordinates although they do not represent 799 765 exactly the same quantities. $\omega$ is provided by the continuity equation 800 766 (see Appendix~\ref{Apdx_A}): 801 802 767 \begin{equation} \label{Eq_PE_sco_continuity} 803 768 \frac{\partial e_3}{\partial t} + e_3 \; \chi + \frac{\partial \omega }{\partial s} = 0 … … 809 774 810 775 \vspace{0.5cm} 811 *tracer equations:776 $\bullet$ tracer equations: 812 777 \begin{multline} \label{Eq_PE_sco_t} 813 778 \frac{1}{e_3} \frac{\partial \left( e_3\,T \right) }{\partial t}= … … 1024 989 1025 990 The $\tilde{z}$-coordinate has been developed by \citet{Leclair_Madec_OM10s}. 1026 It is not available in the current version of \NEMO. 991 It is available in \NEMO since the version 3.4. Nevertheless, it is currently not robust enough 992 to be used in all possible configurations. Its use is therefore not recommended. 993 We 1027 994 1028 995 \newpage … … 1180 1147 ocean (see Appendix~\ref{Apdx_B}). 1181 1148 1149 For \textit{iso-level} diffusion, $r_1$ and $r_2 $ are zero. $\Re $ reduces to the identity 1150 in the horizontal direction, no rotation is applied. 1151 1182 1152 For \textit{geopotential} diffusion, $r_1$ and $r_2 $ are the slopes between the 1183 geopotential 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 1185 equal to $\sigma _1$ and $\sigma _2$, respectively (see \eqref{Eq_PE_sco_slope} ). 1153 geopotential and computational surfaces: they are equal to $\sigma _1$ and $\sigma _2$, 1154 respectively (see \eqref{Eq_PE_sco_slope} ). 1186 1155 1187 1156 For \textit{isoneutral} diffusion $r_1$ and $r_2$ are the slopes between the isoneutral … … 1235 1204 The lateral fourth order tracer diffusive operator is defined by: 1236 1205 \begin{equation} \label{Eq_PE_bilapT} 1237 D^{lT}=\Delta \left( {A^{lT}\;\Delta T}\right)1238 \qquad \text{where} \ D^{lT}=\Delta \left( {A^{lT}\;\Delta T} \right)1206 D^{lT}=\Delta \left( \;\Delta T \right) 1207 \qquad \text{where} \;\; \Delta \bullet = \nabla \left( {\sqrt{B^{lT}\,}\;\Re \;\nabla \bullet} \right) 1239 1208 \end{equation} 1240 1241 1209 It is the second order operator given by \eqref{Eq_PE_iso_tensor} applied twice with 1242 the eddy diffusion coefficient correctly placed.1210 the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. 1243 1211 1244 1212 … … 1262 1230 1263 1231 Such a formulation ensures a complete separation between the vorticity and 1264 horizontal divergence fields (see Appendix~\ref{Apdx_C}). Unfortunately, it is not1265 available for geopotential diffusion in $s-$coordinates and for isoneutral1266 diffusion in both $z$- and $s$-coordinates ($i.e.$ when a rotation is required).1267 In these two cases, the $u$ and $v-$fields are considered as independent scalar1268 fields, so that the diffusive operator is given by:1232 horizontal divergence fields (see Appendix~\ref{Apdx_C}). 1233 Unfortunately, it is only available in \textit{iso-level} direction. 1234 When a rotation is required ($i.e.$ geopotential diffusion in $s-$coordinates 1235 or isoneutral diffusion in both $z$- and $s$-coordinates), the $u$ and $v-$fields 1236 are considered as independent scalar fields, so that the diffusive operator is given by: 1269 1237 \begin{equation} \label{Eq_PE_lapU_iso} 1270 1238 \begin{split} 1271 D_u^{l{\rm {\bf U}}} &= \nabla .\left( { \Re \;\nabla u} \right) \\1272 D_v^{l{\rm {\bf U}}} &= \nabla .\left( { \Re \;\nabla v} \right)1239 D_u^{l{\rm {\bf U}}} &= \nabla .\left( {A^{lm} \;\Re \;\nabla u} \right) \\ 1240 D_v^{l{\rm {\bf U}}} &= \nabla .\left( {A^{lm} \;\Re \;\nabla v} \right) 1273 1241 \end{split} 1274 1242 \end{equation} … … 1282 1250 1283 1251 As for tracers, the fourth order momentum diffusive operator along $z$ or $s$-surfaces 1284 is a re-entering second order operator \eqref{Eq_PE_lapU} or \eqref{Eq_PE_lapU} 1285 with the eddy viscosity coefficient correctly placed: 1286 1287 geopotential diffusion in $z$-coordinate: 1288 \begin{equation} \label{Eq_PE_bilapU} 1289 \begin{split} 1290 {\rm {\bf D}}^{l{\rm {\bf U}}} &=\nabla _h \left\{ {\;\nabla _h {\rm {\bf 1291 .}}\left[ {A^{lm}\,\nabla _h \left( \chi \right)} \right]\;} 1292 \right\}\; \\ 1293 &+\nabla _h \times \left\{ {\;{\rm {\bf k}}\cdot \nabla \times 1294 \left[ {A^{lm}\,\nabla _h \times \left( {\zeta \;{\rm {\bf k}}} \right)} 1295 \right]\;} \right\} 1296 \end{split} 1297 \end{equation} 1298 1299 \gmcomment{ change the position of the coefficient, both here and in the code} 1300 1301 geopotential diffusion in $s$-coordinate: 1302 \begin{equation} \label{Eq_bilapU_iso} 1303 \left\{ \begin{aligned} 1304 D_u^{l{\rm {\bf U}}} =\Delta \left( {A^{lm}\;\Delta u} \right) \\ 1305 D_v^{l{\rm {\bf U}}} =\Delta \left( {A^{lm}\;\Delta v} \right) 1306 \end{aligned} \right. 1307 \quad \text{where} \quad 1308 \Delta \left( \bullet \right) = \nabla \cdot \left( \Re \nabla(\bullet) \right) 1309 \end{equation} 1310 1252 is a re-entering second order operator \eqref{Eq_PE_lapU} or \eqref{Eq_PE_lapU_iso} 1253 with the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. 1254 -
branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_Model_Basics_zstar.tex
r4147 r6040 1 1 % ================================================================ 2 % Chapter 1 �Model Basics2 % Chapter 1 ——— Model Basics 3 3 % ================================================================ 4 4 % ================================================================ -
branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_STP.tex
r4147 r6040 1 1 2 2 % ================================================================ 3 % Chapter 2 �Time Domain (step.F90)3 % Chapter 2 ——— Time Domain (step.F90) 4 4 % ================================================================ 5 5 \chapter{Time Domain (STP) } … … 21 21 22 22 Having defined the continuous equations in Chap.~\ref{PE}, we need now to choose 23 a time discretization. In the present chapter, we provide a general description of the \NEMO 23 a time discretization, a key feature of an ocean model as it exerts a strong influence 24 on the structure of the computer code ($i.e.$ on its flowchart). 25 In the present chapter, we provide a general description of the \NEMO 24 26 time stepping strategy and the consequences for the order in which the equations are 25 27 solved. … … 158 160 \end{equation} 159 161 162 %%gm 163 %%gm UPDATE the next paragraphs with time varying thickness ... 164 %%gm 165 160 166 This scheme is rather time consuming since it requires a matrix inversion, 161 167 but it becomes attractive since a value of 3 or more is needed for N in … … 188 194 189 195 % ------------------------------------------------------------------------------------------------------------- 190 % Hydrostatic Pressure gradient 191 % ------------------------------------------------------------------------------------------------------------- 192 \section{Hydrostatic Pressure Gradient --- semi-implicit scheme} 193 \label{STP_hpg_imp} 196 % Surface Pressure gradient 197 % ------------------------------------------------------------------------------------------------------------- 198 \section{Surface Pressure Gradient} 199 \label{STP_spg_ts} 200 201 ===>>>> TO BE written.... :-) 194 202 195 203 %\gmcomment{ … … 209 217 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 210 218 %} 211 212 The range of stability of the Leap-Frog scheme can be extended by a factor of two213 by introducing a semi-implicit computation of the hydrostatic pressure gradient term214 \citep{Brown_Campana_MWR78}. Instead of evaluating the pressure at $t$, a linear215 combination of values at $t-\rdt$, $t$ and $t+\rdt$ is used (see \S~\ref{DYN_hpg_imp}).216 This technique, controlled by the \np{nn\_dynhpg\_rst} namelist parameter, does not217 introduce a significant additional computational cost when tracers and thus density218 is time stepped before the dynamics. This time step ordering is used in \NEMO219 (Fig.\ref{Fig_TimeStep_flowchart}).220 221 222 This technique, used in several GCMs (\NEMO, POP or MOM for instance),223 makes the Leap-Frog scheme as efficient224 \footnote{The efficiency is defined as the maximum allowed Courant number of the time225 stepping scheme divided by the number of computations of the right-hand side per time step.}226 as the Forward-Backward scheme used in MOM \citep{Griffies_al_OS05} and more227 efficient than the LF-AM3 scheme (leapfrog time stepping combined with a third order228 Adams-Moulthon interpolation for the predictor phase) used in ROMS229 \citep{Shchepetkin_McWilliams_OM05}.230 231 In fact, this technique is efficient when the physical phenomenon that232 limits the time-step is internal gravity waves (IGWs). Indeed, it is233 equivalent to applying a time filter to the pressure gradient to eliminate high234 frequency IGWs. Obviously, the doubling of the time-step is achievable only235 if no other factors control the time-step, such as the stability limits associated236 with advection, diffusion or Coriolis terms. For example, it is inefficient in low resolution237 global ocean configurations, since inertial oscillations in the vicinity of the North Pole238 are the limiting factor for the time step. It is also often inefficient in very high239 resolution configurations where strong currents and small grid cells exert240 the strongest constraint on the time step.241 219 242 220 % ------------------------------------------------------------------------------------------------------------- -
branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_TRA.tex
r5102 r6040 1 1 % ================================================================ 2 % Chapter 1 �Ocean Tracers (TRA)2 % Chapter 1 ——— Ocean Tracers (TRA) 3 3 % ================================================================ 4 4 \chapter{Ocean Tracers (TRA)} … … 40 40 described in chapters \S\ref{SBC}, \S\ref{LDF} and \S\ref{ZDF}, respectively. 41 41 Note that \mdl{tranpc}, the non-penetrative convection module, although 42 (temporarily) located in the NEMO/OPA/TRA directory, is described with the 43 model vertical physics (ZDF). 44 %%% 45 \gmcomment{change the position of eosbn2 in the reference code} 46 %%% 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. 47 45 48 46 In the present chapter we also describe the diagnostic equations used to compute … … 50 48 freezing point with associated modules \mdl{eosbn2} and \mdl{phycst}). 51 49 52 The different options available to the user are managed by namelist logicals or 53 CPP keys. For each equation term \textit{ttt}, the namelist logicals are \textit{ln\_trattt\_xxx},50 The different options available to the user are managed by namelist logicals or CPP keys. 51 For each equation term \textit{TTT}, the namelist logicals are \textit{ln\_traTTT\_xxx}, 54 52 where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme. 55 The CPP key (when it exists) is \textbf{key\_tra ttt}. The equivalent code can be56 found in the \textit{tra ttt} or \textit{trattt\_xxx} module, in the NEMO/OPA/TRA directory.53 The CPP key (when it exists) is \textbf{key\_traTTT}. The equivalent code can be 54 found in the \textit{traTTT} or \textit{traTTT\_xxx} module, in the NEMO/OPA/TRA directory. 57 55 58 56 The user has the option of extracting each tendency term on the rhs of the tracer 59 equation for output (\ key{trdtra} is defined), as described in Chap.~\ref{MISC}.57 equation for output (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}~=~true), as described in Chap.~\ref{MISC}. 60 58 61 59 $\ $\newline % force a new ligne … … 82 80 implicitly requires the use of the continuity equation. Indeed, it is obtained 83 81 by using the following equality : $\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$ 84 which results from the use of the continuity equation, $\nabla \cdot \vect{U}=0$ or85 $ \partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ in constant volume or variable volume case, respectively.82 which results from the use of the continuity equation, $\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ 83 (which reduces to $\nabla \cdot \vect{U}=0$ in linear free surface, $i.e.$ \np{ln\_linssh}=true). 86 84 Therefore it is of paramount importance to design the discrete analogue of the 87 85 advection tendency so that it is consistent with the continuity equation in order to 88 86 enforce the conservation properties of the continuous equations. In other words, 89 by replacing $\tau$ by the number 1in (\ref{Eq_tra_adv}) we recover the discrete form of87 by setting $\tau = 1$ in (\ref{Eq_tra_adv}) we recover the discrete form of 90 88 the continuity equation which is used to calculate the vertical velocity. 91 89 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 113 111 boundary condition depends on the type of sea surface chosen: 114 112 \begin{description} 115 \item [linear free surface:] the first level thickness is constant in time:113 \item [linear free surface:] (\np{ln\_linssh}=true) the first level thickness is constant in time: 116 114 the vertical boundary condition is applied at the fixed surface $z=0$ 117 115 rather than on the moving surface $z=\eta$. There is a non-zero advective … … 119 117 $\left. {\tau _w } \right|_{k=1/2} =T_{k=1} $, $i.e.$ 120 118 the product of surface velocity (at $z=0$) by the first level tracer value. 121 \item [non-linear free surface:] (\ key{vvl} is defined)119 \item [non-linear free surface:] (\np{ln\_linssh}=false) 122 120 convergence/divergence in the first ocean level moves the free surface 123 121 up/down. There is no tracer advection through it so that the advective … … 125 123 \end{description} 126 124 In all cases, this boundary condition retains local conservation of tracer. 127 Global conservation is obtained in both rigid-lid and non-linear free surface128 cases, but not in the linear free surface case. Nevertheless, in the latter 129 case,it is achieved to a good approximation since the non-conservative125 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 130 128 term is the product of the time derivative of the tracer and the free surface 131 height, two quantities that are not correlated (see \S\ref{PE_free_surface}, 132 and also \citet{Roullet_Madec_JGR00,Griffies_al_MWR01,Campin2004}). 129 height, two quantities that are not correlated \citep{Roullet_Madec_JGR00,Griffies_al_MWR01,Campin2004}. 133 130 134 131 The velocity field that appears in (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_zco}) 135 is the centred (\textit{now}) \textit{eulerian} ocean velocity (see Chap.~\ref{DYN}). 136 When 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. 138 139 The choice of an advection scheme is made in the \textit{\ngn{nam\_traadv}} namelist, by 140 setting to \textit{true} one and only one of the logicals \textit{ln\_traadv\_xxx}. The 132 is the centred (\textit{now}) \textit{effective} ocean velocity, $i.e.$ the \textit{eulerian} velocity 133 (see Chap.~\ref{DYN}) plus the eddy induced velocity (\textit{eiv}) 134 and/or the mixed layer eddy induced velocity (\textit{eiv}) 135 when those parameterisations are used (see Chap.~\ref{LDF}). 136 137 The choice of an advection scheme is made in the \textit{\ngn{namtra\_adv}} namelist, by 138 setting to \textit{true} one of the logicals \textit{ln\_traadv\_xxx}. The 141 139 corresponding code can be found in the \textit{traadv\_xxx.F90} module, where 142 \textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme. Details 143 of the advection schemes are given below. The choice of an advection scheme 140 \textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme. 141 By default ($i.e.$ in the reference namelist, \ngn{namelist\_ref}), all the logicals 142 are set to \textit{false}. If the user does not select an advection scheme 143 in the configuration namelist (\ngn{namelist\_cfg}), the tracers will not be advected ! 144 145 Details of the advection schemes are given below. The choice of an advection scheme 144 146 is a complex matter which depends on the model physics, model resolution, 145 147 type of tracer, as well as the issue of numerical cost. 146 148 147 149 Note that 148 (1) cen2, cen4 and TVD schemes require an explicit diffusion 149 operator while the other schemes are diffusive enough so that they do not 150 require additional diffusion ; 151 (2) cen2, cen4, MUSCL2, and UBS are not \textit{positive} schemes 150 (1) CEN and FCT schemes require an explicit diffusion operator 151 while the other schemes are diffusive enough so that they do not necessarily require additional diffusion ; 152 (2) CEN and UBS are not \textit{positive} schemes 152 153 \footnote{negative values can appear in an initially strictly positive tracer field 153 154 which is advected} … … 163 164 164 165 % ------------------------------------------------------------------------------------------------------------- 166 % 2nd and 4th order centred schemes 167 % ------------------------------------------------------------------------------------------------------------- 168 \subsection [$2^{nd}$ and $4^{th}$ order centred schemes (CEN) (\np{ln\_traadv\_cen})] 169 {$2^{nd}$ and $4^{th}$ order centred schemes (CEN) (\np{ln\_traadv\_cen}=true)} 170 \label{TRA_adv_cen} 171 165 172 % 2nd order centred scheme 166 % ------------------------------------------------------------------------------------------------------------- 167 \subsection [$2^{nd}$ order centred scheme (cen2) (\np{ln\_traadv\_cen2})] 168 {$2^{nd}$ order centred scheme (cen2) (\np{ln\_traadv\_cen2}=true)} 169 \label{TRA_adv_cen2} 170 171 In the centred second order formulation, the tracer at velocity points is 173 174 In the $2^{nd}$ order centred formulation (CEN2), the tracer at velocity points is 172 175 evaluated as the mean of the two neighbouring $T$-point values. 173 176 For example, in the $i$-direction : … … 176 179 \end{equation} 177 180 178 The schemeis non diffusive ($i.e.$ it conserves the tracer variance, $\tau^2)$181 CEN2 is non diffusive ($i.e.$ it conserves the tracer variance, $\tau^2)$ 179 182 but dispersive ($i.e.$ it may create false extrema). It is therefore notoriously 180 183 noisy and must be used in conjunction with an explicit diffusion operator to 181 184 produce a sensible solution. The associated time-stepping is performed using 182 185 a leapfrog scheme in conjunction with an Asselin time-filter, so $T$ in 183 (\ref{Eq_tra_adv_cen2}) is the \textit{now} tracer value. The centered second 184 order advection is computed in the \mdl{traadv\_cen2} module. In this module, 185 it is advantageous to combine the \textit{cen2} scheme with an upstream scheme 186 in specific areas which require a strong diffusion in order to avoid the generation 187 of false extrema. These areas are the vicinity of large river mouths, some straits 188 with coarse resolution, and the vicinity of ice cover area ($i.e.$ when the ocean 189 temperature is close to the freezing point). 190 This combined scheme has been included for specific grid points in the ORCA2 191 and ORCA4 configurations only. This is an obsolescent feature as the recommended 192 advection scheme for the ORCA configuration is TVD (see \S\ref{TRA_adv_tvd}). 193 194 Note that using the cen2 scheme, the overall tracer advection is of second 186 (\ref{Eq_tra_adv_cen2}) is the \textit{now} tracer value. 187 CEN2 is computed in the \mdl{traadv\_cen} module. 188 189 Note that using the CEN2, the overall tracer advection is of second 195 190 order accuracy since both (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_cen2}) 196 191 have this order of accuracy. \gmcomment{Note also that ... blah, blah} 197 192 198 % -------------------------------------------------------------------------------------------------------------199 193 % 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 205 In the $4^{th}$ order formulation (to be implemented), tracer values are 206 evaluated at velocity points as a $4^{th}$ order interpolation, and thus depend on 207 the four neighbouring $T$-points. For example, in the $i$-direction: 194 195 In the $4^{th}$ order formulation (CEN4), tracer values are evaluated at velocity points as 196 a $4^{th}$ order interpolation, and thus depend on the four neighbouring $T$-points. 197 For example, in the $i$-direction: 208 198 \begin{equation} \label{Eq_tra_adv_cen4} 209 199 \tau _u^{cen4} … … 211 201 \end{equation} 212 202 213 Strictly speaking, the cen4 scheme is not a $4^{th}$ order advection scheme203 Strictly speaking, the CEN4 scheme is not a $4^{th}$ order advection scheme 214 204 but a $4^{th}$ order evaluation of advective fluxes, since the divergence of 215 advective fluxes \eqref{Eq_tra_adv} is kept at $2^{nd}$ order. The phrase ``$4^{th}$ 216 order scheme'' used in oceanographic literature is usually associated 217 with the scheme presented here. Introducing a \textit{true} $4^{th}$ order advection 218 scheme is feasible but, for consistency reasons, it requires changes in the 219 discretisation of the tracer advection together with changes in both the 220 continuity equation and the momentum advection terms. 205 advective fluxes \eqref{Eq_tra_adv} is kept at $2^{nd}$ order. 206 The expression \textit{$4^{th}$ order scheme} used in oceanographic literature 207 is usually associated with the scheme presented here. 208 Introducing a \textit{true} $4^{th}$ order advection scheme is feasible but, 209 for consistency reasons, it requires changes in the discretisation of the tracer 210 advection together with changes in the continuity equation, 211 and the momentum advection and pressure terms. 221 212 222 213 A direct consequence of the pseudo-fourth order nature of the scheme is that 223 it is not non-diffusive, i.e.the global variance of a tracer is not preserved using224 \textit{cen4}. Furthermore, it must be used in conjunction with an explicit225 diffusion operator to produce a sensible solution. The time-stepping is also214 it is not non-diffusive, $i.e.$ the global variance of a tracer is not preserved using 215 CEN4. Furthermore, it must be used in conjunction with an explicit 216 diffusion operator to produce a sensible solution. As in CEN2 case, the time-stepping is 226 217 performed using a leapfrog scheme in conjunction with an Asselin time-filter, 227 218 so $T$ in (\ref{Eq_tra_adv_cen4}) is the \textit{now} tracer. … … 235 226 236 227 % ------------------------------------------------------------------------------------------------------------- 237 % TVDscheme238 % ------------------------------------------------------------------------------------------------------------- 239 \subsection [ Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd})]240 { Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd}=true)}228 % FCT scheme 229 % ------------------------------------------------------------------------------------------------------------- 230 \subsection [$2^{nd}$ and $4^{th}$ Flux Corrected Transport schemes (FCT) (\np{ln\_traadv\_fct})] 231 {$2^{nd}$ and $4^{th}$ Flux Corrected Transport schemes (FCT) (\np{ln\_traadv\_fct}=true)} 241 232 \label{TRA_adv_tvd} 242 233 243 In the Total Variance Dissipation (TVD)formulation, the tracer at velocity234 In the Flux Corrected Transport formulation, the tracer at velocity 244 235 points is evaluated using a combination of an upstream and a centred scheme. 245 236 For example, in the $i$-direction : 246 \begin{equation} \label{Eq_tra_adv_ tvd}237 \begin{equation} \label{Eq_tra_adv_fct} 247 238 \begin{split} 248 239 \tau _u^{ups}&= \begin{cases} … … 251 242 \end{cases} \\ 252 243 \\ 253 \tau _u^{ tvd}&=\tau _u^{ups} +c_u \;\left( {\tau _u^{cen2} -\tau _u^{ups} } \right)244 \tau _u^{fct}&=\tau _u^{ups} +c_u \;\left( {\tau _u^{cen} -\tau _u^{ups} } \right) 254 245 \end{split} 255 246 \end{equation} … … 260 251 produces a local extremum in the tracer field. The resulting scheme is quite 261 252 expensive but \emph{positive}. It can be used on both active and passive tracers. 262 This scheme is tested and compared with MUSCL and the MPDATA scheme in 263 \citet{Levy_al_GRL01}; note that in this paper it is referred to as "FCT" (Flux corrected 264 transport) rather than TVD. The TVD scheme is implemented in the \mdl{traadv\_tvd} module. 253 This scheme is tested and compared with MUSCL and a MPDATA scheme in \citet{Levy_al_GRL01}. 254 The FCT scheme is implemented in the \mdl{traadv\_fct} module. 265 255 266 256 For stability reasons (see \S\ref{STP}), 267 $\tau _u^{cen 2}$ is evaluated in (\ref{Eq_tra_adv_tvd}) using the \textit{now} tracer while $\tau _u^{ups}$268 is evaluated using the \textit{before} tracer. In other words, the advective part of269 the scheme is time stepped with a leap-frog scheme while a forward scheme is270 used for the diffusive part.257 $\tau _u^{cen}$ is evaluated in (\ref{Eq_tra_adv_fct}) using the \textit{now} tracer 258 while $\tau _u^{ups}$ is evaluated using the \textit{before} tracer. In other words, 259 the advective part of the scheme is time stepped with a leap-frog scheme 260 while a forward scheme is used for the diffusive part. 271 261 272 262 % ------------------------------------------------------------------------------------------------------------- 273 263 % MUSCL scheme 274 264 % ------------------------------------------------------------------------------------------------------------- 275 \subsection[MUSCL scheme (\np{ln\_traadv\_mus cl})]276 {Monotone Upstream Scheme for Conservative Laws (MUSCL) (\np{ln\_traadv\_mus cl}=T)}277 \label{TRA_adv_mus cl}265 \subsection[MUSCL scheme (\np{ln\_traadv\_mus})] 266 {Monotone Upstream Scheme for Conservative Laws (MUSCL) (\np{ln\_traadv\_mus}=T)} 267 \label{TRA_adv_mus} 278 268 279 269 The Monotone Upstream Scheme for Conservative Laws (MUSCL) has been … … 281 271 is evaluated assuming a linear tracer variation between two $T$-points 282 272 (Fig.\ref{Fig_adv_scheme}). For example, in the $i$-direction : 283 \begin{equation} \label{Eq_tra_adv_mus cl}273 \begin{equation} \label{Eq_tra_adv_mus} 284 274 \tau _u^{mus} = \left\{ \begin{aligned} 285 275 &\tau _i &+ \frac{1}{2} \;\left( 1-\frac{u_{i+1/2} \;\rdt}{e_{1u}} \right) … … 296 286 297 287 For an ocean grid point adjacent to land and where the ocean velocity is 298 directed 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 301 the \textit{positive} character of the scheme. Only the former can be used 302 on both active and passive tracers. The two MUSCL schemes are implemented 303 in the \mdl{traadv\_tvd} and \mdl{traadv\_tvd2} modules. 288 directed toward land, an upstream flux is used. This choice ensure 289 the \textit{positive} character of the scheme. 304 290 305 291 % ------------------------------------------------------------------------------------------------------------- … … 310 296 \label{TRA_adv_ubs} 311 297 312 The UBS advection scheme is an upstream-biased third order scheme based on313 an upstream-biased parabolic interpolation. It is also known as the Cell314 Averaged QUICK scheme (Quadratic Upstream Interpolation for Convective298 The UBS advection scheme (also often called UP3) is an upstream-biased third order 299 scheme based on an upstream-biased parabolic interpolation. It is also known as 300 the Cell Averaged QUICK scheme (Quadratic Upstream Interpolation for Convective 315 301 Kinematics). For example, in the $i$-direction : 316 302 \begin{equation} \label{Eq_tra_adv_ubs} … … 324 310 325 311 This results in a dissipatively dominant (i.e. hyper-diffusive) truncation 326 error \citep{Shchepetkin_McWilliams_OM05}. The overall performance of the advection327 scheme is similar to that reported in \cite{Farrow1995}.312 error \citep{Shchepetkin_McWilliams_OM05}. The overall performance of 313 the advection scheme is similar to that reported in \cite{Farrow1995}. 328 314 It is a relatively good compromise between accuracy and smoothness. 329 315 It is not a \emph{positive} scheme, meaning that false extrema are permitted, 330 316 but the amplitude of such are significantly reduced over the centred second 331 or der method. Nevertheless it is not recommended that it should be applied332 to a passive tracer that requires positivity.317 or fourth order method. Nevertheless it is not recommended that it should be 318 applied to a passive tracer that requires positivity. 333 319 334 320 The intrinsic diffusion of UBS makes its use risky in the vertical direction 335 321 where the control of artificial diapycnal fluxes is of paramount importance. 336 Therefore the vertical flux is evaluated using the TVD scheme when337 \np{ln\_traadv\_ubs}=true.322 Therefore the vertical flux is evaluated using either a 2nd order FCT scheme 323 or a 4th order COMPACT scheme (\np{nn\_cen\_v}=2 or 4). 338 324 339 325 For stability reasons (see \S\ref{STP}), 340 the first term in \eqref{Eq_tra_adv_ubs} (which corresponds to a second order centred scheme)341 is evaluated using the \textit{now} tracer (centred in time) while the342 second term (which is the diffusive part of the scheme), is326 the first term in \eqref{Eq_tra_adv_ubs} (which corresponds to a second order 327 centred scheme) is evaluated using the \textit{now} tracer (centred in time) 328 while the second term (which is the diffusive part of the scheme), is 343 329 evaluated using the \textit{before} tracer (forward in time). 344 330 This choice is discussed by \citet{Webb_al_JAOT98} in the context of the … … 350 336 substitution in the \mdl{traadv\_ubs} module and obtain a QUICK scheme. 351 337 338 ??? 339 352 340 Four different options are possible for the vertical 353 341 component used in the UBS scheme. $\tau _w^{ubs}$ can be evaluated 354 342 using either \textit{(a)} a centred $2^{nd}$ order scheme, or \textit{(b)} 355 a TVDscheme, or \textit{(c)} an interpolation based on conservative343 a FCT scheme, or \textit{(c)} an interpolation based on conservative 356 344 parabolic splines following the \citet{Shchepetkin_McWilliams_OM05} 357 345 implementation of UBS in ROMS, or \textit{(d)} a UBS. The $3^{rd}$ case 358 346 has dispersion properties similar to an eighth-order accurate conventional scheme. 359 The current reference version uses method b) 347 The current reference version uses method (b). 348 349 ??? 360 350 361 351 Note that : … … 390 380 Thirdly, the diffusion term is in fact a biharmonic operator with an eddy 391 381 coefficient which is simply proportional to the velocity: 392 $A_u^{lm}= - \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note that NEMO v3.4still uses382 $A_u^{lm}= \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note the current version of NEMO still uses 393 383 \eqref{Eq_tra_adv_ubs}, not \eqref{Eq_traadv_ubs2}. 394 384 %%% … … 416 406 direction (as for the UBS case) should be implemented to restore this property. 417 407 418 419 % -------------------------------------------------------------------------------------------------------------420 % PPM scheme421 % -------------------------------------------------------------------------------------------------------------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 426 The Piecewise Parabolic Method (PPM) proposed by Colella and Woodward (1984)427 \sgacomment{reference?}428 is based on a quadradic piecewise construction. Like the QCK scheme, it is associated429 with the ULTIMATE QUICKEST limiter \citep{Leonard1991}. It has been implemented430 in \NEMO by G. Reffray (MERCATOR-ocean) but is not yet offered in the reference431 version 3.3.432 408 433 409 % ================================================================ … … 1167 1143 % Equation of State 1168 1144 % ------------------------------------------------------------------------------------------------------------- 1169 \subsection{Equation of State (\np{nn\_eos} = 0, 1 or 2)}1145 \subsection{Equation Of Seawater (\np{nn\_eos} = -1, 0, or 1)} 1170 1146 \label{TRA_eos} 1171 1147 1172 It is necessary to know the equation of state for the ocean very accurately 1173 to determine stability properties (especially the Brunt-Vais\"{a}l\"{a} frequency), 1174 particularly in the deep ocean. The ocean seawater volumic mass, $\rho$, 1175 abusively called density, is a non linear empirical function of \textit{in situ} 1176 temperature, salinity and pressure. The reference equation of state is that 1177 defined by the Joint Panel on Oceanographic Tables and Standards 1178 \citep{UNESCO1983}. It was the standard equation of state used in early 1179 releases of OPA. However, even though this computation is fully vectorised, 1180 it is quite time consuming ($15$ to $20${\%} of the total CPU time) since 1181 it requires the prior computation of the \textit{in situ} temperature from the 1182 model \textit{potential} temperature using the \citep{Bryden1973} polynomial 1183 for adiabatic lapse rate and a $4^th$ order Runge-Kutta integration scheme. 1184 Since OPA6, we have used the \citet{JackMcD1995} equation of state for 1185 seawater instead. It allows the computation of the \textit{in situ} ocean density 1186 directly as a function of \textit{potential} temperature relative to the surface 1187 (an \NEMO variable), the practical salinity (another \NEMO variable) and the 1188 pressure (assuming no pressure variation along geopotential surfaces, $i.e.$ 1189 the pressure in decibars is approximated by the depth in meters). 1190 Both the \citet{UNESCO1983} and \citet{JackMcD1995} equations of state 1191 have exactly the same except that the values of the various coefficients have 1192 been adjusted by \citet{JackMcD1995} in order to directly use the \textit{potential} 1193 temperature 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, 1195 while maintaining a quite accurate equation of state. 1196 1197 In the computer code, a \textit{true} density anomaly, $d_a= \rho / \rho_o - 1$, 1198 is computed, with $\rho_o$ a reference volumic mass. Called \textit{rau0} 1199 in the code, $\rho_o$ is defined in \mdl{phycst}, and a value of $1,035~Kg/m^3$. 1148 The Equation Of Seawater (EOS) is an empirical nonlinear thermodynamic relationship 1149 linking seawater density, $\rho$, to a number of state variables, 1150 most typically temperature, salinity and pressure. 1151 Because density gradients control the pressure gradient force through the hydrostatic balance, 1152 the equation of state provides a fundamental bridge between the distribution of active tracers 1153 and the fluid dynamics. Nonlinearities of the EOS are of major importance, in particular 1154 influencing the circulation through determination of the static stability below the mixed layer, 1155 thus controlling rates of exchange between the atmosphere and the ocean interior \citep{Roquet_JPO2015}. 1156 Therefore an accurate EOS based on either the 1980 equation of state (EOS-80, \cite{UNESCO1983}) 1157 or TEOS-10 \citep{TEOS10} standards should be used anytime a simulation of the real 1158 ocean circulation is attempted \citep{Roquet_JPO2015}. 1159 The use of TEOS-10 is highly recommended because 1160 \textit{(i)} it is the new official EOS, 1161 \textit{(ii)} it is more accurate, being based on an updated database of laboratory measurements, and 1162 \textit{(iii)} it uses Conservative Temperature and Absolute Salinity (instead of potential temperature 1163 and practical salinity for EOS-980, both variables being more suitable for use as model variables 1164 \citep{TEOS10, Graham_McDougall_JPO13}. 1165 EOS-80 is an obsolescent feature of the NEMO system, kept only for backward compatibility. 1166 For process studies, it is often convenient to use an approximation of the EOS. To that purposed, 1167 a simplified EOS (S-EOS) inspired by \citet{Vallis06} is also available. 1168 1169 In the computer code, a density anomaly, $d_a= \rho / \rho_o - 1$, 1170 is computed, with $\rho_o$ a reference density. Called \textit{rau0} 1171 in the code, $\rho_o$ is set in \mdl{phycst} to a value of $1,026~Kg/m^3$. 1200 1172 This is a sensible choice for the reference density used in a Boussinesq ocean 1201 1173 climate model, as, with the exception of only a small percentage of the ocean, 1202 density in the World Ocean varies by no more than 2$\%$ from $1,035~kg/m^3$ 1203 \citep{Gill1982}. 1204 1205 Options are defined through the \ngn{nameos} namelist variables. 1206 The default option (namelist parameter \np{nn\_eos}=0) is the \citet{JackMcD1995} 1207 equation of state. Its use is highly recommended. However, for process studies, 1208 it is often convenient to use a linear approximation of the density. 1174 density in the World Ocean varies by no more than 2$\%$ from that value \citep{Gill1982}. 1175 1176 Options are defined through the \ngn{nameos} namelist variables, and in particular \np{nn\_eos} 1177 which controls the EOS used (=-1 for TEOS10 ; =0 for EOS-80 ; =1 for S-EOS). 1178 \begin{description} 1179 1180 \item[\np{nn\_eos}$=-1$] the polyTEOS10-bsq equation of seawater \citep{Roquet_OM2015} is used. 1181 The accuracy of this approximation is comparable to the TEOS-10 rational function approximation, 1182 but it is optimized for a boussinesq fluid and the polynomial expressions have simpler 1183 and more computationally efficient expressions for their derived quantities 1184 which make them more adapted for use in ocean models. 1185 Note that a slightly higher precision polynomial form is now used replacement of the TEOS-10 1186 rational function approximation for hydrographic data analysis \citep{TEOS10}. 1187 A key point is that conservative state variables are used: 1188 Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: $\degres C$, notation: $\Theta$). 1189 The pressure in decibars is approximated by the depth in meters. 1190 With TEOS10, the specific heat capacity of sea water, $C_p$, is a constant. It is set to 1191 $C_p=3991.86795711963~J\,Kg^{-1}\,\degres K^{-1}$, according to \citet{TEOS10}. 1192 1193 Choosing polyTEOS10-bsq implies that the state variables used by the model are 1194 $\Theta$ and $S_A$. In particular, the initial state deined by the user have to be given as 1195 \textit{Conservative} Temperature and \textit{Absolute} Salinity. 1196 In addition, setting \np{ln\_useCT} to \textit{true} convert the Conservative SST to potential SST 1197 prior to either computing the air-sea and ice-sea fluxes (forced mode) 1198 or sending the SST field to the atmosphere (coupled mode). 1199 1200 \item[\np{nn\_eos}$=0$] the polyEOS80-bsq equation of seawater is used. 1201 It takes the same polynomial form as the polyTEOS10, but the coefficients have been optimized 1202 to accurately fit EOS80 (Roquet, personal comm.). The state variables used in both the EOS80 1203 and the ocean model are: 1204 the Practical Salinity ((unit: psu, notation: $S_p$)) and Potential Temperature (unit: $\degres C$, notation: $\theta$). 1205 The pressure in decibars is approximated by the depth in meters. 1206 With thsi EOS, the specific heat capacity of sea water, $C_p$, is a function of temperature, 1207 salinity and pressure \citep{UNESCO1983}. Nevertheless, a severe assumption is made in order to 1208 have a heat content ($C_p T_p$) which is conserved by the model: $C_p$ is set to a constant 1209 value, the TEOS10 value. 1210 1211 \item[\np{nn\_eos}$=1$] a simplified EOS (S-EOS) inspired by \citet{Vallis06} is chosen, 1212 the coefficients of which has been optimized to fit the behavior of TEOS10 (Roquet, personal comm.) 1213 (see also \citet{Roquet_JPO2015}). It provides a simplistic linear representation of both 1214 cabbeling and thermobaricity effects which is enough for a proper treatment of the EOS 1215 in theoretical studies \citep{Roquet_JPO2015}. 1209 1216 With such an equation of state there is no longer a distinction between 1210 \textit{in situ} and \textit{potential} density and both cabbeling and thermobaric 1211 effects are removed. 1212 Two linear formulations are available: a function of $T$ only (\np{nn\_eos}=1) 1213 and a function of both $T$ and $S$ (\np{nn\_eos}=2): 1214 \begin{equation} \label{Eq_tra_eos_linear} 1217 \textit{conservative} and \textit{potential} temperature, as well as between \textit{absolute} 1218 and \textit{practical} salinity. 1219 S-EOS takes the following expression: 1220 \begin{equation} \label{Eq_tra_S-EOS} 1215 1221 \begin{split} 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 1222 d_a(T,S,z) = ( & - a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * T_a \\ 1223 & + b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * S_a \\ 1224 & - \nu \; T_a \; S_a \; ) \; / \; \rho_o \\ 1225 with \ \ T_a = T-10 \; ; & \; S_a = S-35 \; ;\; \rho_o = 1026~Kg/m^3 1218 1226 \end{split} 1219 1227 \end{equation} 1220 where $\alpha$ and $\beta$ are the thermal and haline expansion 1221 coefficients, 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 1224 of $T$ only (\np{nn\_eos}=1), the salinity is a passive tracer and can be 1225 used as such. 1228 where the computer name of the coefficients as well as their standard value are given in \ref{Tab_SEOS}. 1229 In fact, when choosing S-EOS, various approximation of EOS can be specified simply by changing 1230 the associated coefficients. 1231 Setting to zero the two thermobaric coefficients ($\mu_1$, $\mu_2$) remove thermobaric effect from S-EOS. 1232 setting to zero the three cabbeling coefficients ($\lambda_1$, $\lambda_2$, $\nu$) remove cabbeling effect from S-EOS. 1233 Keeping non-zero value to $a_0$ and $b_0$ provide a linear EOS function of T and S. 1234 1235 \end{description} 1236 1237 1238 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1239 \begin{table}[!tb] 1240 \begin{center} \begin{tabular}{|p{26pt}|p{72pt}|p{56pt}|p{136pt}|} 1241 \hline 1242 coeff. & computer name & S-EOS & description \\ \hline 1243 $a_0$ & \np{nn\_a0} & 1.6550 $10^{-1}$ & linear thermal expansion coeff. \\ \hline 1244 $b_0$ & \np{nn\_b0} & 7.6554 $10^{-1}$ & linear haline expansion coeff. \\ \hline 1245 $\lambda_1$ & \np{nn\_lambda1}& 5.9520 $10^{-2}$ & cabbeling coeff. in $T^2$ \\ \hline 1246 $\lambda_2$ & \np{nn\_lambda2}& 5.4914 $10^{-4}$ & cabbeling coeff. in $S^2$ \\ \hline 1247 $\nu$ & \np{nn\_nu} & 2.4341 $10^{-3}$ & cabbeling coeff. in $T \, S$ \\ \hline 1248 $\mu_1$ & \np{nn\_mu1} & 1.4970 $10^{-4}$ & thermobaric coeff. in T \\ \hline 1249 $\mu_2$ & \np{nn\_mu2} & 1.1090 $10^{-5}$ & thermobaric coeff. in S \\ \hline 1250 \end{tabular} 1251 \caption{ \label{Tab_SEOS} 1252 Standard value of S-EOS coefficients. } 1253 \end{center} 1254 \end{table} 1255 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1256 1226 1257 1227 1258 % ------------------------------------------------------------------------------------------------------------- … … 1232 1263 1233 1264 An 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 1241 For \np{nn\_eos}=0 (\citet{JackMcD1995} equation of state), the \citet{McDougall1987} 1242 polynomial expression is used (with the pressure in decibar approximated by 1243 the depth in meters): 1265 frequency) is of paramount importance as determine the ocean stratification and 1266 is used in several ocean parameterisations (namely TKE, GLS, Richardson number dependent 1267 vertical diffusion, enhanced vertical diffusion, non-penetrative convection, tidal mixing 1268 parameterisation, iso-neutral diffusion). In particular, $N^2$ has to be computed at the local pressure 1269 (pressure in decibar being approximated by the depth in meters). The expression for $N^2$ 1270 is given by: 1244 1271 \begin{equation} \label{Eq_tra_bn2} 1245 N^2 = \frac{g}{e_{3w}} \; \beta \1246 \left( \alpha / \beta \ \delta_{k+1/2}[T] - \delta_{k+1/2}[S] \right)1247 \end{equation}1248 where $\alpha$ and $\beta$ are the thermal and haline expansion coefficients.1249 They are a function of $\overline{T}^{\,k+1/2},\widetilde{S}=\overline{S}^{\,k+1/2} - 35.$,1250 and $z_w$, with $T$ the \textit{potential} temperature and $\widetilde{S}$ a salinity anomaly.1251 Note that both $\alpha$ and $\beta$ depend on \textit{potential}1252 temperature and salinity which are averaged at $w$-points prior1253 to the computation instead of being computed at $T$-points and1254 then averaged to $w$-points.1255 1256 When 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}1259 1272 N^2 = \frac{g}{e_{3w}} \left( \beta \;\delta_{k+1/2}[S] - \alpha \;\delta_{k+1/2}[T] \right) 1260 1273 \end{equation} 1261 where $\alpha$ and $\beta $ are the constant coefficients used to 1262 defined 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 1271 The specific heat of sea water, $C_p$, is a function of temperature, salinity 1272 and pressure \citep{UNESCO1983}. It is only used in the model to convert 1273 surface heat fluxes into surface temperature increase and so the pressure 1274 dependence is neglected. The dependence on $T$ and $S$ is weak. 1275 For example, with $S=35~psu$, $C_p$ increases from $3989$ to $4002$ 1276 when $T$ varies from -2~\degres C to 31~\degres C. Therefore, $C_p$ has 1277 been chosen as a constant: $C_p=4.10^3~J\,Kg^{-1}\,\degres K^{-1}$. 1278 Its value is set in \mdl{phycst} module. 1279 1274 where $(T,S) = (\Theta, S_A)$ for TEOS10, $= (\theta, S_p)$ for TEOS-80, or $=(T,S)$ for S-EOS, 1275 and, $\alpha$ and $\beta$ are the thermal and haline expansion coefficients. 1276 The coefficients are a polynomial function of temperature, salinity and depth which expression 1277 depends on the chosen EOS. They are computed through \textit{eos\_rab}, a \textsc{Fortran} 1278 function that can be found in \mdl{eosbn2}. 1280 1279 1281 1280 % ------------------------------------------------------------------------------------------------------------- … … 1298 1297 sea water ($i.e.$ referenced to the surface $p=0$), thus the pressure dependent 1299 1298 terms in \eqref{Eq_tra_eos_fzp} (last term) have been dropped. The freezing 1300 point is computed through \textit{ tfreez}, a \textsc{Fortran} function that can be found1299 point is computed through \textit{eos\_fzp}, a \textsc{Fortran} function that can be found 1301 1300 in \mdl{eosbn2}. 1301 1302 1303 % ------------------------------------------------------------------------------------------------------------- 1304 % Potential Energy 1305 % ------------------------------------------------------------------------------------------------------------- 1306 %\subsection{Potential Energy anomalies} 1307 %\label{TRA_bn2} 1308 1309 % =====>>>>> TO BE written 1310 % 1311 1302 1312 1303 1313 % ================================================================ -
branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Introduction.tex
r4661 r6040 24 24 release 8.2, described in \citet{Madec1998}. This model has been used for a wide 25 25 range of applications, both regional or global, as a forced ocean model and as a 26 model coupled with the atmosphere. A complete list of references is found on the 27 \NEMO web site. 26 model coupled with the sea-ice and/or the atmosphere. 28 27 29 28 This manual is organised in as follows. Chapter~\ref{PE} presents the model basics, 30 29 $i.e.$ the equations and their assumptions, the vertical coordinates used, and the 31 30 subgrid scale physics. This part deals with the continuous equations of the model 32 (primitive equations, with potential temperature, salinity and an equation of state).31 (primitive equations, with temperature, salinity and an equation of seawater). 33 32 The equations are written in a curvilinear coordinate system, with a choice of vertical 34 coordinates ($z$ or $s$, with the rescaled height coordinate formulation \textit{z*}, or35 \textit{s*}). Momentum equations are formulated in the vector invariant form or in the36 flux form.Dimensional units in the meter, kilogram, second (MKS) international system33 coordinates ($z$, $s$, \textit{z*}, \textit{s*}, $\tilde{z}$, $\tilde{s}$, and a mixture of them). 34 Momentum equations are formulated in the vector invariant form or in the flux form. 35 Dimensional units in the meter, kilogram, second (MKS) international system 37 36 are used throughout. 38 37 … … 79 78 space and time variable coefficient \citet{Treguier1997}. The model has vertical harmonic 80 79 viscosity and diffusion with a space and time variable coefficient, with options to compute 81 the coefficients with \citet{Blanke1993}, \citet{ Large_al_RG94}, \citet{Pacanowski_Philander_JPO81},80 the coefficients with \citet{Blanke1993}, \citet{Pacanowski_Philander_JPO81}, 82 81 or \citet{Umlauf_Burchard_JMS03} mixing schemes. 83 82 \vspace{1cm} 84 83 85 84 %%gm To be put somewhere else .... 85 86 86 \noindent CPP keys and namelists are used for inputs to the code. \newline 87 87 … … 112 112 \vspace{1cm} 113 113 114 %%gm end 114 115 115 116 Model outputs management and specific online diagnostics are described in chapters~\ref{DIA}. … … 249 250 250 251 252 \vspace{1cm} 253 $\bullet$ The main modifications from NEMO/OPA v3.4 and v3.6 are :\\ 254 \begin{enumerate} 255 \item ... ; 256 \end{enumerate} 257 258 259 \vspace{1cm} 260 $\bullet$ The main modifications from NEMO/OPA v3.6 and v4.0 are :\\ 261 \begin{enumerate} 262 \item ... ; 263 \end{enumerate} 264 265
Note: See TracChangeset
for help on using the changeset viewer.