- Timestamp:
- 2016-07-19T10:38:35+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5549_BDY_ZEROGRAD/DOC/TexFiles/Chapters/Chap_TRA.tex
r5102 r6808 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 141 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 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 Several tracer advection scheme are proposed, namely 138 a $2^{nd}$ or $4^{th}$ order centred schemes (CEN), 139 a $2^{nd}$ or $4^{th}$ order Flux Corrected Transport scheme (FCT), 140 a Monotone Upstream Scheme for Conservative Laws scheme (MUSCL), 141 a $3^{rd}$ Upstream Biased Scheme (UBS, also often called UP3), and 142 a Quadratic Upstream Interpolation for Convective Kinematics with 143 Estimated Streaming Terms scheme (QUICKEST). 144 The choice is made in the \textit{\ngn{namtra\_adv}} namelist, by 145 setting to \textit{true} one of the logicals \textit{ln\_traadv\_xxx}. 146 The corresponding code can be found in the \textit{traadv\_xxx.F90} module, 147 where \textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme. 148 By default ($i.e.$ in the reference namelist, \ngn{namelist\_ref}), all the logicals 149 are set to \textit{false}. If the user does not select an advection scheme 150 in the configuration namelist (\ngn{namelist\_cfg}), the tracers will \textit{not} be advected ! 151 152 Details of the advection schemes are given below. The choosing an advection scheme 144 153 is a complex matter which depends on the model physics, model resolution, 145 type of tracer, as well as the issue of numerical cost. 146 147 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 154 type of tracer, as well as the issue of numerical cost. In particular, we note that 155 (1) CEN and FCT schemes require an explicit diffusion operator 156 while the other schemes are diffusive enough so that they do not necessarily need additional diffusion ; 157 (2) CEN and UBS are not \textit{positive} schemes 152 158 \footnote{negative values can appear in an initially strictly positive tracer field 153 159 which is advected} … … 163 169 164 170 % ------------------------------------------------------------------------------------------------------------- 171 % 2nd and 4th order centred schemes 172 % ------------------------------------------------------------------------------------------------------------- 173 \subsection [centred schemes (CEN) (\np{ln\_traadv\_cen})] 174 {centred schemes (CEN) (\np{ln\_traadv\_cen}=true)} 175 \label{TRA_adv_cen} 176 165 177 % 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 172 evaluated as the mean of the two neighbouring $T$-point values. 178 179 The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen}~=~\textit{true}. 180 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) 181 and vertical direction by setting \np{nn\_cen\_h} and \np{nn\_cen\_v} to $2$ or $4$. 182 CEN implementation can be found in the \mdl{traadv\_cen} module. 183 184 In the $2^{nd}$ order centred formulation (CEN2), the tracer at velocity points 185 is evaluated as the mean of the two neighbouring $T$-point values. 173 186 For example, in the $i$-direction : 174 187 \begin{equation} \label{Eq_tra_adv_cen2} … … 176 189 \end{equation} 177 190 178 The schemeis non diffusive ($i.e.$ it conserves the tracer variance, $\tau^2)$191 CEN2 is non diffusive ($i.e.$ it conserves the tracer variance, $\tau^2)$ 179 192 but dispersive ($i.e.$ it may create false extrema). It is therefore notoriously 180 193 noisy and must be used in conjunction with an explicit diffusion operator to 181 194 produce a sensible solution. The associated time-stepping is performed using 182 195 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 196 (\ref{Eq_tra_adv_cen2}) is the \textit{now} tracer value. 197 198 Note that using the CEN2, the overall tracer advection is of second 195 199 order accuracy since both (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_cen2}) 196 have this order of accuracy. \gmcomment{Note also that ... blah, blah} 197 198 % ------------------------------------------------------------------------------------------------------------- 200 have this order of accuracy. 201 199 202 % 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: 203 204 In the $4^{th}$ order formulation (CEN4), tracer values are evaluated at u- and v-points as 205 a $4^{th}$ order interpolation, and thus depend on the four neighbouring $T$-points. 206 For example, in the $i$-direction: 208 207 \begin{equation} \label{Eq_tra_adv_cen4} 209 208 \tau _u^{cen4} 210 209 =\overline{ T - \frac{1}{6}\,\delta _i \left[ \delta_{i+1/2}[T] \,\right] }^{\,i+1/2} 211 210 \end{equation} 212 213 Strictly speaking, the cen4 scheme is not a $4^{th}$ order advection scheme 211 In the vertical direction (\np{nn\_cen\_v}=$4$), a $4^{th}$ COMPACT interpolation 212 has been prefered \citep{Demange_PhD2014}. 213 In the COMPACT scheme, both the field and its derivative are interpolated, 214 which leads, after a matrix inversion, spectral characteristics 215 similar to schemes of higher order \citep{Lele_JCP1992}. 216 217 218 Strictly speaking, the CEN4 scheme is not a $4^{th}$ order advection scheme 214 219 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. 220 advective fluxes \eqref{Eq_tra_adv} is kept at $2^{nd}$ order. 221 The expression \textit{$4^{th}$ order scheme} used in oceanographic literature 222 is usually associated with the scheme presented here. 223 Introducing a \textit{true} $4^{th}$ order advection scheme is feasible but, 224 for consistency reasons, it requires changes in the discretisation of the tracer 225 advection together with changes in the continuity equation, 226 and the momentum advection and pressure terms. 221 227 222 228 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 also226 performed using a leapfrog scheme in conjunction with an Asselin time-filter,227 so $T$ in (\ref{Eq_tra_adv_cen4}) is the \textit{now} tracer.228 229 At a $T$-grid cell adjacent to a boundary (coastline, bottom and surface), an230 a dditional hypothesis must be made to evaluate $\tau _u^{cen4}$. This231 hypothesis usually reduces the order of the scheme. Here we choose to set232 the gradient of $T$ across the boundary to zero. Alternative conditions can be233 specified, such as a reduction to a second order scheme for these near boundary234 grid points.235 236 % ------------------------------------------------------------------------------------------------------------- 237 % TVDscheme238 % ------------------------------------------------------------------------------------------------------------- 239 \subsection [ Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd})]240 { Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd}=true)}229 it is not non-diffusive, $i.e.$ the global variance of a tracer is not preserved using CEN4. 230 Furthermore, it must be used in conjunction with an explicit diffusion operator 231 to produce a sensible solution. 232 As in CEN2 case, the time-stepping is performed using a leapfrog scheme in conjunction 233 with an Asselin time-filter, so $T$ in (\ref{Eq_tra_adv_cen4}) is the \textit{now} tracer. 234 235 At a $T$-grid cell adjacent to a boundary (coastline, bottom and surface), 236 an additional hypothesis must be made to evaluate $\tau _u^{cen4}$. 237 This hypothesis usually reduces the order of the scheme. 238 Here we choose to set the gradient of $T$ across the boundary to zero. 239 Alternative conditions can be specified, such as a reduction to a second order scheme 240 for these near boundary grid points. 241 242 % ------------------------------------------------------------------------------------------------------------- 243 % FCT scheme 244 % ------------------------------------------------------------------------------------------------------------- 245 \subsection [Flux Corrected Transport schemes (FCT) (\np{ln\_traadv\_fct})] 246 {Flux Corrected Transport schemes (FCT) (\np{ln\_traadv\_fct}=true)} 241 247 \label{TRA_adv_tvd} 242 248 243 In the Total Variance Dissipation (TVD) formulation, the tracer at velocity 244 points is evaluated using a combination of an upstream and a centred scheme. 245 For example, in the $i$-direction : 246 \begin{equation} \label{Eq_tra_adv_tvd} 249 The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct}~=~\textit{true}. 250 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) 251 and vertical direction by setting \np{nn\_fct\_h} and \np{nn\_fct\_v} to $2$ or $4$. 252 FCT implementation can be found in the \mdl{traadv\_fct} module. 253 254 In FCT formulation, the tracer at velocity points is evaluated using a combination of 255 an upstream and a centred scheme. For example, in the $i$-direction : 256 \begin{equation} \label{Eq_tra_adv_fct} 247 257 \begin{split} 248 258 \tau _u^{ups}&= \begin{cases} … … 251 261 \end{cases} \\ 252 262 \\ 253 \tau _u^{ tvd}&=\tau _u^{ups} +c_u \;\left( {\tau _u^{cen2} -\tau _u^{ups} } \right)263 \tau _u^{fct}&=\tau _u^{ups} +c_u \;\left( {\tau _u^{cen} -\tau _u^{ups} } \right) 254 264 \end{split} 255 265 \end{equation} 256 266 where $c_u$ is a flux limiter function taking values between 0 and 1. 267 The FCT order is the one of the centred scheme used ($i.e.$ it depends on the setting of 268 \np{nn\_fct\_h} and \np{nn\_fct\_v}. 257 269 There exist many ways to define $c_u$, each corresponding to a different 258 total variance decreasing scheme. The one chosen in \NEMO is described in 259 \citet{Zalesak_JCP79}. $c_u$ only departs from $1$ when the advective term 260 produces a local extremum in the tracer field. The resulting scheme is quite 261 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. 265 266 For stability reasons (see \S\ref{STP}), 267 $\tau _u^{cen2}$ 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 of 269 the scheme is time stepped with a leap-frog scheme while a forward scheme is 270 used for the diffusive part. 270 FCT scheme. The one chosen in \NEMO is described in \citet{Zalesak_JCP79}. 271 $c_u$ only departs from $1$ when the advective term produces a local extremum in the tracer field. 272 The resulting scheme is quite expensive but \emph{positive}. 273 It can be used on both active and passive tracers. 274 A comparison of FCT-2 with MUSCL and a MPDATA scheme can be found in \citet{Levy_al_GRL01}. 275 276 An additional option has been added controlled by \np{nn\_fct\_zts}. By setting this integer to 277 a value larger than zero, a $2^{nd}$ order FCT scheme is used on both horizontal and vertical direction, 278 but on the latter, a split-explicit time stepping is used, with a number of sub-timestep equals 279 to \np{nn\_fct\_zts}. This option can be useful when the size of the timestep is limited 280 by vertical advection \citep{Lemarie_OM2015)}. Note that in this case, a similar split-explicit 281 time stepping should be used on vertical advection of momentum to insure a better stability 282 (see \S\ref{DYN_zad}). 283 284 For stability reasons (see \S\ref{STP}), $\tau _u^{cen}$ is evaluated in (\ref{Eq_tra_adv_fct}) 285 using the \textit{now} tracer while $\tau _u^{ups}$ is evaluated using the \textit{before} tracer. In other words, 286 the advective part of the scheme is time stepped with a leap-frog scheme 287 while a forward scheme is used for the diffusive part. 271 288 272 289 % ------------------------------------------------------------------------------------------------------------- 273 290 % MUSCL scheme 274 291 % ------------------------------------------------------------------------------------------------------------- 275 \subsection[MUSCL scheme (\np{ln\_traadv\_muscl})] 276 {Monotone Upstream Scheme for Conservative Laws (MUSCL) (\np{ln\_traadv\_muscl}=T)} 277 \label{TRA_adv_muscl} 278 279 The Monotone Upstream Scheme for Conservative Laws (MUSCL) has been 280 implemented by \citet{Levy_al_GRL01}. In its formulation, the tracer at velocity points 292 \subsection[MUSCL scheme (\np{ln\_traadv\_mus})] 293 {Monotone Upstream Scheme for Conservative Laws (MUSCL) (\np{ln\_traadv\_mus}=T)} 294 \label{TRA_adv_mus} 295 296 The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln\_traadv\_mus}~=~\textit{true}. 297 MUSCL implementation can be found in the \mdl{traadv\_mus} module. 298 299 MUSCL has been first implemented in \NEMO by \citet{Levy_al_GRL01}. In its formulation, the tracer at velocity points 281 300 is evaluated assuming a linear tracer variation between two $T$-points 282 301 (Fig.\ref{Fig_adv_scheme}). For example, in the $i$-direction : 283 \begin{equation} \label{Eq_tra_adv_mus cl}302 \begin{equation} \label{Eq_tra_adv_mus} 284 303 \tau _u^{mus} = \left\{ \begin{aligned} 285 304 &\tau _i &+ \frac{1}{2} \;\left( 1-\frac{u_{i+1/2} \;\rdt}{e_{1u}} \right) … … 296 315 297 316 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. 317 directed toward land, an upstream flux is used. This choice ensure 318 the \textit{positive} character of the scheme. 319 In addition, fluxes round a grid-point where a runoff is applied can optionally be 320 computed using upstream fluxes (\np{ln\_mus\_ups}~=~\textit{true}). 304 321 305 322 % ------------------------------------------------------------------------------------------------------------- … … 310 327 \label{TRA_adv_ubs} 311 328 312 The UBS advection scheme is an upstream-biased third order scheme based on 313 an upstream-biased parabolic interpolation. It is also known as the Cell 314 Averaged QUICK scheme (Quadratic Upstream Interpolation for Convective 315 Kinematics). For example, in the $i$-direction : 329 The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs}~=~\textit{true}. 330 UBS implementation can be found in the \mdl{traadv\_mus} module. 331 332 The UBS scheme, often called UP3, is also known as the Cell Averaged QUICK scheme 333 (Quadratic Upstream Interpolation for Convective Kinematics). It is an upstream-biased 334 third order scheme based on an upstream-biased parabolic interpolation. 335 For example, in the $i$-direction : 316 336 \begin{equation} \label{Eq_tra_adv_ubs} 317 337 \tau _u^{ubs} =\overline T ^{i+1/2}-\;\frac{1}{6} \left\{ … … 324 344 325 345 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}.346 error \citep{Shchepetkin_McWilliams_OM05}. The overall performance of 347 the advection scheme is similar to that reported in \cite{Farrow1995}. 328 348 It is a relatively good compromise between accuracy and smoothness. 329 It is not a \emph{positive} scheme, meaning that false extrema are permitted,349 Nevertheless the scheme is not \emph{positive}, meaning that false extrema are permitted, 330 350 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.351 or fourth order method. therefore it is not recommended that it should be 352 applied to a passive tracer that requires positivity. 333 353 334 354 The intrinsic diffusion of UBS makes its use risky in the vertical direction 335 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.355 where the control of artificial diapycnal fluxes is of paramount importance \citep{Shchepetkin_McWilliams_OM05, Demange_PhD2014}. 356 Therefore the vertical flux is evaluated using either a $2^nd$ order FCT scheme 357 or a $4^th$ order COMPACT scheme (\np{nn\_cen\_v}=2 or 4). 338 358 339 359 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), is360 the first term in \eqref{Eq_tra_adv_ubs} (which corresponds to a second order 361 centred scheme) is evaluated using the \textit{now} tracer (centred in time) 362 while the second term (which is the diffusive part of the scheme), is 343 363 evaluated using the \textit{before} tracer (forward in time). 344 364 This choice is discussed by \citet{Webb_al_JAOT98} in the context of the … … 350 370 substitution in the \mdl{traadv\_ubs} module and obtain a QUICK scheme. 351 371 352 Four different options are possible for the vertical 353 component used in the UBS scheme. $\tau _w^{ubs}$ can be evaluated 354 using either \textit{(a)} a centred $2^{nd}$ order scheme, or \textit{(b)} 355 a TVD scheme, or \textit{(c)} an interpolation based on conservative 356 parabolic splines following the \citet{Shchepetkin_McWilliams_OM05} 357 implementation of UBS in ROMS, or \textit{(d)} a UBS. The $3^{rd}$ case 358 has dispersion properties similar to an eighth-order accurate conventional scheme. 359 The current reference version uses method b) 360 361 Note that : 362 363 (1) When a high vertical resolution $O(1m)$ is used, the model stability can 364 be controlled by vertical advection (not vertical diffusion which is usually 365 solved using an implicit scheme). Computer time can be saved by using a 366 time-splitting technique on vertical advection. Such a technique has been 367 implemented and validated in ORCA05 with 301 levels. It is not available 368 in the current reference version. 369 370 (2) It is straightforward to rewrite \eqref{Eq_tra_adv_ubs} as follows: 372 Note that it is straightforward to rewrite \eqref{Eq_tra_adv_ubs} as follows: 371 373 \begin{equation} \label{Eq_traadv_ubs2} 372 374 \tau _u^{ubs} = \tau _u^{cen4} + \frac{1}{12} \left\{ … … 390 392 Thirdly, the diffusion term is in fact a biharmonic operator with an eddy 391 393 coefficient which is simply proportional to the velocity: 392 $A_u^{lm}= - \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note that NEMO v3.4 still uses 393 \eqref{Eq_tra_adv_ubs}, not \eqref{Eq_traadv_ubs2}. 394 %%% 395 \gmcomment{the change in UBS scheme has to be done} 396 %%% 394 $A_u^{lm}= \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note the current version of NEMO uses 395 the computationally more efficient formulation \eqref{Eq_tra_adv_ubs}. 397 396 398 397 % ------------------------------------------------------------------------------------------------------------- … … 405 404 The Quadratic Upstream Interpolation for Convective Kinematics with 406 405 Estimated Streaming Terms (QUICKEST) scheme proposed by \citet{Leonard1979} 407 is the third order Godunov scheme. It is associated with the ULTIMATE QUICKEST 406 is used when \np{ln\_traadv\_qck}~=~\textit{true}. 407 QUICKEST implementation can be found in the \mdl{traadv\_mus} module. 408 409 QUICKEST is the third order Godunov scheme which is associated with the ULTIMATE QUICKEST 408 410 limiter \citep{Leonard1991}. It has been implemented in NEMO by G. Reffray 409 411 (MERCATOR-ocean) and can be found in the \mdl{traadv\_qck} module. … … 416 418 direction (as for the UBS case) should be implemented to restore this property. 417 419 418 419 % ------------------------------------------------------------------------------------------------------------- 420 % PPM scheme 421 % ------------------------------------------------------------------------------------------------------------- 422 \subsection [Piecewise Parabolic Method (PPM) (\np{ln\_traadv\_ppm})] 423 {Piecewise Parabolic Method (PPM) (\np{ln\_traadv\_ppm}=true)} 424 \label{TRA_adv_ppm} 425 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 associated 429 with the ULTIMATE QUICKEST limiter \citep{Leonard1991}. It has been implemented 430 in \NEMO by G. Reffray (MERCATOR-ocean) but is not yet offered in the reference 431 version 3.3. 420 %%%gmcomment : Cross term are missing in the current implementation.... 421 432 422 433 423 % ================================================================ … … 1167 1157 % Equation of State 1168 1158 % ------------------------------------------------------------------------------------------------------------- 1169 \subsection{Equation of State (\np{nn\_eos} = 0, 1 or 2)}1159 \subsection{Equation Of Seawater (\np{nn\_eos} = -1, 0, or 1)} 1170 1160 \label{TRA_eos} 1171 1161 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$. 1162 The Equation Of Seawater (EOS) is an empirical nonlinear thermodynamic relationship 1163 linking seawater density, $\rho$, to a number of state variables, 1164 most typically temperature, salinity and pressure. 1165 Because density gradients control the pressure gradient force through the hydrostatic balance, 1166 the equation of state provides a fundamental bridge between the distribution of active tracers 1167 and the fluid dynamics. Nonlinearities of the EOS are of major importance, in particular 1168 influencing the circulation through determination of the static stability below the mixed layer, 1169 thus controlling rates of exchange between the atmosphere and the ocean interior \citep{Roquet_JPO2015}. 1170 Therefore an accurate EOS based on either the 1980 equation of state (EOS-80, \cite{UNESCO1983}) 1171 or TEOS-10 \citep{TEOS10} standards should be used anytime a simulation of the real 1172 ocean circulation is attempted \citep{Roquet_JPO2015}. 1173 The use of TEOS-10 is highly recommended because 1174 \textit{(i)} it is the new official EOS, 1175 \textit{(ii)} it is more accurate, being based on an updated database of laboratory measurements, and 1176 \textit{(iii)} it uses Conservative Temperature and Absolute Salinity (instead of potential temperature 1177 and practical salinity for EOS-980, both variables being more suitable for use as model variables 1178 \citep{TEOS10, Graham_McDougall_JPO13}. 1179 EOS-80 is an obsolescent feature of the NEMO system, kept only for backward compatibility. 1180 For process studies, it is often convenient to use an approximation of the EOS. To that purposed, 1181 a simplified EOS (S-EOS) inspired by \citet{Vallis06} is also available. 1182 1183 In the computer code, a density anomaly, $d_a= \rho / \rho_o - 1$, 1184 is computed, with $\rho_o$ a reference density. Called \textit{rau0} 1185 in the code, $\rho_o$ is set in \mdl{phycst} to a value of $1,026~Kg/m^3$. 1200 1186 This is a sensible choice for the reference density used in a Boussinesq ocean 1201 1187 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. 1188 density in the World Ocean varies by no more than 2$\%$ from that value \citep{Gill1982}. 1189 1190 Options are defined through the \ngn{nameos} namelist variables, and in particular \np{nn\_eos} 1191 which controls the EOS used (=-1 for TEOS10 ; =0 for EOS-80 ; =1 for S-EOS). 1192 \begin{description} 1193 1194 \item[\np{nn\_eos}$=-1$] the polyTEOS10-bsq equation of seawater \citep{Roquet_OM2015} is used. 1195 The accuracy of this approximation is comparable to the TEOS-10 rational function approximation, 1196 but it is optimized for a boussinesq fluid and the polynomial expressions have simpler 1197 and more computationally efficient expressions for their derived quantities 1198 which make them more adapted for use in ocean models. 1199 Note that a slightly higher precision polynomial form is now used replacement of the TEOS-10 1200 rational function approximation for hydrographic data analysis \citep{TEOS10}. 1201 A key point is that conservative state variables are used: 1202 Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: $\degres C$, notation: $\Theta$). 1203 The pressure in decibars is approximated by the depth in meters. 1204 With TEOS10, the specific heat capacity of sea water, $C_p$, is a constant. It is set to 1205 $C_p=3991.86795711963~J\,Kg^{-1}\,\degres K^{-1}$, according to \citet{TEOS10}. 1206 1207 Choosing polyTEOS10-bsq implies that the state variables used by the model are 1208 $\Theta$ and $S_A$. In particular, the initial state deined by the user have to be given as 1209 \textit{Conservative} Temperature and \textit{Absolute} Salinity. 1210 In addition, setting \np{ln\_useCT} to \textit{true} convert the Conservative SST to potential SST 1211 prior to either computing the air-sea and ice-sea fluxes (forced mode) 1212 or sending the SST field to the atmosphere (coupled mode). 1213 1214 \item[\np{nn\_eos}$=0$] the polyEOS80-bsq equation of seawater is used. 1215 It takes the same polynomial form as the polyTEOS10, but the coefficients have been optimized 1216 to accurately fit EOS80 (Roquet, personal comm.). The state variables used in both the EOS80 1217 and the ocean model are: 1218 the Practical Salinity ((unit: psu, notation: $S_p$)) and Potential Temperature (unit: $\degres C$, notation: $\theta$). 1219 The pressure in decibars is approximated by the depth in meters. 1220 With thsi EOS, the specific heat capacity of sea water, $C_p$, is a function of temperature, 1221 salinity and pressure \citep{UNESCO1983}. Nevertheless, a severe assumption is made in order to 1222 have a heat content ($C_p T_p$) which is conserved by the model: $C_p$ is set to a constant 1223 value, the TEOS10 value. 1224 1225 \item[\np{nn\_eos}$=1$] a simplified EOS (S-EOS) inspired by \citet{Vallis06} is chosen, 1226 the coefficients of which has been optimized to fit the behavior of TEOS10 (Roquet, personal comm.) 1227 (see also \citet{Roquet_JPO2015}). It provides a simplistic linear representation of both 1228 cabbeling and thermobaricity effects which is enough for a proper treatment of the EOS 1229 in theoretical studies \citep{Roquet_JPO2015}. 1209 1230 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} 1231 \textit{conservative} and \textit{potential} temperature, as well as between \textit{absolute} 1232 and \textit{practical} salinity. 1233 S-EOS takes the following expression: 1234 \begin{equation} \label{Eq_tra_S-EOS} 1215 1235 \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 1236 d_a(T,S,z) = ( & - a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * T_a \\ 1237 & + b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * S_a \\ 1238 & - \nu \; T_a \; S_a \; ) \; / \; \rho_o \\ 1239 with \ \ T_a = T-10 \; ; & \; S_a = S-35 \; ;\; \rho_o = 1026~Kg/m^3 1218 1240 \end{split} 1219 1241 \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. 1242 where the computer name of the coefficients as well as their standard value are given in \ref{Tab_SEOS}. 1243 In fact, when choosing S-EOS, various approximation of EOS can be specified simply by changing 1244 the associated coefficients. 1245 Setting to zero the two thermobaric coefficients ($\mu_1$, $\mu_2$) remove thermobaric effect from S-EOS. 1246 setting to zero the three cabbeling coefficients ($\lambda_1$, $\lambda_2$, $\nu$) remove cabbeling effect from S-EOS. 1247 Keeping non-zero value to $a_0$ and $b_0$ provide a linear EOS function of T and S. 1248 1249 \end{description} 1250 1251 1252 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1253 \begin{table}[!tb] 1254 \begin{center} \begin{tabular}{|p{26pt}|p{72pt}|p{56pt}|p{136pt}|} 1255 \hline 1256 coeff. & computer name & S-EOS & description \\ \hline 1257 $a_0$ & \np{rn\_a0} & 1.6550 $10^{-1}$ & linear thermal expansion coeff. \\ \hline 1258 $b_0$ & \np{rn\_b0} & 7.6554 $10^{-1}$ & linear haline expansion coeff. \\ \hline 1259 $\lambda_1$ & \np{rn\_lambda1}& 5.9520 $10^{-2}$ & cabbeling coeff. in $T^2$ \\ \hline 1260 $\lambda_2$ & \np{rn\_lambda2}& 5.4914 $10^{-4}$ & cabbeling coeff. in $S^2$ \\ \hline 1261 $\nu$ & \np{rn\_nu} & 2.4341 $10^{-3}$ & cabbeling coeff. in $T \, S$ \\ \hline 1262 $\mu_1$ & \np{rn\_mu1} & 1.4970 $10^{-4}$ & thermobaric coeff. in T \\ \hline 1263 $\mu_2$ & \np{rn\_mu2} & 1.1090 $10^{-5}$ & thermobaric coeff. in S \\ \hline 1264 \end{tabular} 1265 \caption{ \label{Tab_SEOS} 1266 Standard value of S-EOS coefficients. } 1267 \end{center} 1268 \end{table} 1269 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1270 1226 1271 1227 1272 % ------------------------------------------------------------------------------------------------------------- … … 1232 1277 1233 1278 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): 1279 frequency) is of paramount importance as determine the ocean stratification and 1280 is used in several ocean parameterisations (namely TKE, GLS, Richardson number dependent 1281 vertical diffusion, enhanced vertical diffusion, non-penetrative convection, tidal mixing 1282 parameterisation, iso-neutral diffusion). In particular, $N^2$ has to be computed at the local pressure 1283 (pressure in decibar being approximated by the depth in meters). The expression for $N^2$ 1284 is given by: 1244 1285 \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 1286 N^2 = \frac{g}{e_{3w}} \left( \beta \;\delta_{k+1/2}[S] - \alpha \;\delta_{k+1/2}[T] \right) 1260 1287 \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 1288 where $(T,S) = (\Theta, S_A)$ for TEOS10, $= (\theta, S_p)$ for TEOS-80, or $=(T,S)$ for S-EOS, 1289 and, $\alpha$ and $\beta$ are the thermal and haline expansion coefficients. 1290 The coefficients are a polynomial function of temperature, salinity and depth which expression 1291 depends on the chosen EOS. They are computed through \textit{eos\_rab}, a \textsc{Fortran} 1292 function that can be found in \mdl{eosbn2}. 1280 1293 1281 1294 % ------------------------------------------------------------------------------------------------------------- … … 1298 1311 sea water ($i.e.$ referenced to the surface $p=0$), thus the pressure dependent 1299 1312 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 found1313 point is computed through \textit{eos\_fzp}, a \textsc{Fortran} function that can be found 1301 1314 in \mdl{eosbn2}. 1315 1316 1317 % ------------------------------------------------------------------------------------------------------------- 1318 % Potential Energy 1319 % ------------------------------------------------------------------------------------------------------------- 1320 %\subsection{Potential Energy anomalies} 1321 %\label{TRA_bn2} 1322 1323 % =====>>>>> TO BE written 1324 % 1325 1302 1326 1303 1327 % ================================================================
Note: See TracChangeset
for help on using the changeset viewer.