Changeset 11564 for NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/doc/latex/NEMO/subfiles/chap_TRA.tex
- Timestamp:
- 2019-09-18T16:11:52+02:00 (5 years ago)
- Location:
- NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/doc
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/doc
- Property svn:ignore deleted
-
NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/doc/latex
- Property svn:ignore
-
old new 1 *.aux 2 *.bbl 3 *.blg 4 *.dvi 5 *.fdb* 6 *.fls 7 *.idx 8 *.ilg 9 *.ind 10 *.log 11 *.maf 12 *.mtc* 13 *.out 14 *.pdf 15 *.toc 16 _minted-* 1 figures
-
- Property svn:ignore
-
NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/doc/latex/NEMO
- Property svn:ignore deleted
-
NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/doc/latex/NEMO/subfiles/chap_TRA.tex
r10544 r11564 8 8 \label{chap:TRA} 9 9 10 \ minitoc11 12 % missing/update 10 \chaptertoc 11 12 % missing/update 13 13 % traqsr: need to coordinate with SBC module 14 14 15 %STEVEN : is the use of the word "positive" to describe a scheme enough, or should it be "positive definite"? I added a comment to this effect on some instances of this below 15 %STEVEN : is the use of the word "positive" to describe a scheme enough, or should it be "positive definite"? 16 %I added a comment to this effect on some instances of this below 16 17 17 18 Using the representation described in \autoref{chap:DOM}, several semi -discrete space forms of … … 35 36 The terms QSR, BBC, BBL and DMP are optional. 36 37 The external forcings and parameterisations require complex inputs and complex calculations 37 (\eg bulk formulae, estimation of mixing coefficients) that are carried out in the SBC,38 (\eg\ bulk formulae, estimation of mixing coefficients) that are carried out in the SBC, 38 39 LDF and ZDF modules and described in \autoref{chap:SBC}, \autoref{chap:LDF} and 39 40 \autoref{chap:ZDF}, respectively. … … 47 48 associated modules \mdl{eosbn2} and \mdl{phycst}). 48 49 49 The different options available to the user are managed by namelist logicals or CPP keys.50 The different options available to the user are managed by namelist logicals. 50 51 For each equation term \textit{TTT}, the namelist logicals are \textit{ln\_traTTT\_xxx}, 51 52 where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme. 52 The CPP key (when it exists) is \key{traTTT}.53 53 The equivalent code can be found in the \textit{traTTT} or \textit{traTTT\_xxx} module, 54 54 in the \path{./src/OCE/TRA} directory. 55 55 56 56 The user has the option of extracting each tendency term on the RHS of the tracer equation for output 57 (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl} ~\forcode{=.true.}), as described in \autoref{chap:DIA}.57 (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}\forcode{=.true.}), as described in \autoref{chap:DIA}. 58 58 59 59 % ================================================================ 60 60 % Tracer Advection 61 61 % ================================================================ 62 \section{Tracer advection (\protect\mdl{traadv})} 62 \section[Tracer advection (\textit{traadv.F90})] 63 {Tracer advection (\protect\mdl{traadv})} 63 64 \label{sec:TRA_adv} 64 65 %------------------------------------------namtra_adv----------------------------------------------------- 65 66 66 \nlst{namtra_adv} 67 \begin{listing} 68 \nlst{namtra_adv} 69 \caption{\texttt{namtra\_adv}} 70 \label{lst:namtra_adv} 71 \end{listing} 67 72 %------------------------------------------------------------------------------------------------------------- 68 73 69 When considered (\ie when \np{ln\_traadv\_NONE} is not set to \forcode{.true.}),74 When considered (\ie\ when \np{ln\_traadv\_OFF} is not set to \forcode{.true.}), 70 75 the advection tendency of a tracer is expressed in flux form, 71 \ie as the divergence of the advective fluxes.76 \ie\ as the divergence of the advective fluxes. 72 77 Its discrete expression is given by : 73 78 \begin{equation} 74 \label{eq: tra_adv}79 \label{eq:TRA_adv} 75 80 ADV_\tau = - \frac{1}{b_t} \Big( \delta_i [ e_{2u} \, e_{3u} \; u \; \tau_u] 76 81 + \delta_j [ e_{1v} \, e_{3v} \; v \; \tau_v] \Big) … … 78 83 \end{equation} 79 84 where $\tau$ is either T or S, and $b_t = e_{1t} \, e_{2t} \, e_{3t}$ is the volume of $T$-cells. 80 The flux form in \autoref{eq: tra_adv} implicitly requires the use of the continuity equation.85 The flux form in \autoref{eq:TRA_adv} implicitly requires the use of the continuity equation. 81 86 Indeed, it is obtained by using the following equality: $\nabla \cdot (\vect U \, T) = \vect U \cdot \nabla T$ which 82 87 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, \ie \np{ln\_linssh}~\forcode{=.true.}).88 (which reduces to $\nabla \cdot \vect U = 0$ in linear free surface, \ie\ \np{ln\_linssh}\forcode{=.true.}). 84 89 Therefore it is of paramount importance to design the discrete analogue of the advection tendency so that 85 90 it is consistent with the continuity equation in order to enforce the conservation properties of 86 91 the continuous equations. 87 In other words, by setting $\tau = 1$ in (\autoref{eq: tra_adv}) we recover the discrete form of92 In other words, by setting $\tau = 1$ in (\autoref{eq:TRA_adv}) we recover the discrete form of 88 93 the continuity equation which is used to calculate the vertical velocity. 89 94 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 90 95 \begin{figure}[!t] 91 \begin{center} 92 \includegraphics[]{Fig_adv_scheme} 93 \caption{ 94 \protect\label{fig:adv_scheme} 95 Schematic representation of some ways used to evaluate the tracer value at $u$-point and 96 the amount of tracer exchanged between two neighbouring grid points. 97 Upsteam biased scheme (ups): 98 the upstream value is used and the black area is exchanged. 99 Piecewise parabolic method (ppm): 100 a parabolic interpolation is used and the black and dark grey areas are exchanged. 101 Monotonic upstream scheme for conservative laws (muscl): 102 a parabolic interpolation is used and black, dark grey and grey areas are exchanged. 103 Second order scheme (cen2): 104 the mean value is used and black, dark grey, grey and light grey areas are exchanged. 105 Note that this illustration does not include the flux limiter used in ppm and muscl schemes. 106 } 107 \end{center} 96 \centering 97 \includegraphics[width=\textwidth]{Fig_adv_scheme} 98 \caption[Ways to evaluate the tracer value and the amount of tracer exchanged]{ 99 Schematic representation of some ways used to evaluate the tracer value at $u$-point and 100 the amount of tracer exchanged between two neighbouring grid points. 101 Upsteam biased scheme (ups): 102 the upstream value is used and the black area is exchanged. 103 Piecewise parabolic method (ppm): 104 a parabolic interpolation is used and the black and dark grey areas are exchanged. 105 Monotonic upstream scheme for conservative laws (muscl): 106 a parabolic interpolation is used and black, dark grey and grey areas are exchanged. 107 Second order scheme (cen2): 108 the mean value is used and black, dark grey, grey and light grey areas are exchanged. 109 Note that this illustration does not include the flux limiter used in ppm and muscl schemes.} 110 \label{fig:TRA_adv_scheme} 108 111 \end{figure} 109 112 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 110 113 111 The key difference between the advection schemes available in \NEMO is the choice made in space and114 The key difference between the advection schemes available in \NEMO\ is the choice made in space and 112 115 time interpolation to define the value of the tracer at the velocity points 113 (\autoref{fig: adv_scheme}).116 (\autoref{fig:TRA_adv_scheme}). 114 117 115 118 Along solid lateral and bottom boundaries a zero tracer flux is automatically specified, … … 119 122 \begin{description} 120 123 \item[linear free surface:] 121 (\np{ln\_linssh} ~\forcode{=.true.})124 (\np{ln\_linssh}\forcode{=.true.}) 122 125 the first level thickness is constant in time: 123 126 the vertical boundary condition is applied at the fixed surface $z = 0$ rather than on 124 127 the moving surface $z = \eta$. 125 128 There is a non-zero advective flux which is set for all advection schemes as 126 $\tau_w|_{k = 1/2} = T_{k = 1}$, \ie the product of surface velocity (at $z = 0$) by129 $\tau_w|_{k = 1/2} = T_{k = 1}$, \ie\ the product of surface velocity (at $z = 0$) by 127 130 the first level tracer value. 128 131 \item[non-linear free surface:] 129 (\np{ln\_linssh} ~\forcode{=.false.})132 (\np{ln\_linssh}\forcode{=.false.}) 130 133 convergence/divergence in the first ocean level moves the free surface up/down. 131 134 There is no tracer advection through it so that the advective fluxes through the surface are also zero. … … 136 139 Nevertheless, in the latter case, it is achieved to a good approximation since 137 140 the non-conservative term is the product of the time derivative of the tracer and the free surface height, 138 two quantities that are not correlated \citep{ Roullet_Madec_JGR00, Griffies_al_MWR01, Campin2004}.139 140 The velocity field that appears in (\autoref{eq: tra_adv} and \autoref{eq:tra_adv_zco})is141 the centred (\textit{now}) \textit{effective} ocean velocity, \ie the \textit{eulerian} velocity141 two quantities that are not correlated \citep{roullet.madec_JGR00, griffies.pacanowski.ea_MWR01, campin.adcroft.ea_OM04}. 142 143 The velocity field that appears in (\autoref{eq:TRA_adv} is 144 the centred (\textit{now}) \textit{effective} ocean velocity, \ie\ the \textit{eulerian} velocity 142 145 (see \autoref{chap:DYN}) plus the eddy induced velocity (\textit{eiv}) and/or 143 146 the mixed layer eddy induced velocity (\textit{eiv}) when those parameterisations are used … … 148 151 Conservative Laws scheme (MUSCL), a $3^{rd}$ Upstream Biased Scheme (UBS, also often called UP3), 149 152 and a Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms scheme (QUICKEST). 150 The choice is made in the \n gn{namtra\_adv} namelist, by setting to \forcode{.true.} one of153 The choice is made in the \nam{tra\_adv} namelist, by setting to \forcode{.true.} one of 151 154 the logicals \textit{ln\_traadv\_xxx}. 152 155 The corresponding code can be found in the \textit{traadv\_xxx.F90} module, where 153 156 \textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme. 154 By default (\ie in the reference namelist, \textit{namelist\_ref}), all the logicals are set to \forcode{.false.}.157 By default (\ie\ in the reference namelist, \textit{namelist\_ref}), all the logicals are set to \forcode{.false.}. 155 158 If the user does not select an advection scheme in the configuration namelist (\textit{namelist\_cfg}), 156 159 the tracers will \textit{not} be advected! … … 183 186 % 2nd and 4th order centred schemes 184 187 % ------------------------------------------------------------------------------------------------------------- 185 \subsection{CEN: Centred scheme (\protect\np{ln\_traadv\_cen}~\forcode{= .true.})} 188 \subsection[CEN: Centred scheme (\forcode{ln_traadv_cen=.true.})] 189 {CEN: Centred scheme (\protect\np{ln\_traadv\_cen}\forcode{=.true.})} 186 190 \label{subsec:TRA_adv_cen} 187 191 188 % 2nd order centred scheme 189 190 The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen} ~\forcode{=.true.}.192 % 2nd order centred scheme 193 194 The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen}\forcode{=.true.}. 191 195 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) and vertical direction by 192 196 setting \np{nn\_cen\_h} and \np{nn\_cen\_v} to $2$ or $4$. … … 197 201 For example, in the $i$-direction : 198 202 \begin{equation} 199 \label{eq: tra_adv_cen2}203 \label{eq:TRA_adv_cen2} 200 204 \tau_u^{cen2} = \overline T ^{i + 1/2} 201 205 \end{equation} 202 206 203 CEN2 is non diffusive (\ie it conserves the tracer variance, $\tau^2$) but dispersive204 (\ie it may create false extrema).207 CEN2 is non diffusive (\ie\ it conserves the tracer variance, $\tau^2$) but dispersive 208 (\ie\ it may create false extrema). 205 209 It is therefore notoriously noisy and must be used in conjunction with an explicit diffusion operator to 206 210 produce a sensible solution. 207 211 The associated time-stepping is performed using a leapfrog scheme in conjunction with an Asselin time-filter, 208 so $T$ in (\autoref{eq: tra_adv_cen2}) is the \textit{now} tracer value.212 so $T$ in (\autoref{eq:TRA_adv_cen2}) is the \textit{now} tracer value. 209 213 210 214 Note that using the CEN2, the overall tracer advection is of second order accuracy since 211 both (\autoref{eq: tra_adv}) and (\autoref{eq:tra_adv_cen2}) have this order of accuracy.212 213 % 4nd order centred scheme 215 both (\autoref{eq:TRA_adv}) and (\autoref{eq:TRA_adv_cen2}) have this order of accuracy. 216 217 % 4nd order centred scheme 214 218 215 219 In the $4^{th}$ order formulation (CEN4), tracer values are evaluated at u- and v-points as … … 217 221 For example, in the $i$-direction: 218 222 \begin{equation} 219 \label{eq: tra_adv_cen4}223 \label{eq:TRA_adv_cen4} 220 224 \tau_u^{cen4} = \overline{T - \frac{1}{6} \, \delta_i \Big[ \delta_{i + 1/2}[T] \, \Big]}^{\,i + 1/2} 221 225 \end{equation} 222 In the vertical direction (\np{nn\_cen\_v} ~\forcode{=4}),223 a $4^{th}$ COMPACT interpolation has been prefered \citep{ Demange_PhD2014}.226 In the vertical direction (\np{nn\_cen\_v}\forcode{=4}), 227 a $4^{th}$ COMPACT interpolation has been prefered \citep{demange_phd14}. 224 228 In the COMPACT scheme, both the field and its derivative are interpolated, which leads, after a matrix inversion, 225 spectral characteristics similar to schemes of higher order \citep{ Lele_JCP1992}.229 spectral characteristics similar to schemes of higher order \citep{lele_JCP92}. 226 230 227 231 Strictly speaking, the CEN4 scheme is not a $4^{th}$ order advection scheme but 228 232 a $4^{th}$ order evaluation of advective fluxes, 229 since the divergence of advective fluxes \autoref{eq: tra_adv} is kept at $2^{nd}$ order.233 since the divergence of advective fluxes \autoref{eq:TRA_adv} is kept at $2^{nd}$ order. 230 234 The expression \textit{$4^{th}$ order scheme} used in oceanographic literature is usually associated with 231 235 the scheme presented here. … … 235 239 236 240 A direct consequence of the pseudo-fourth order nature of the scheme is that it is not non-diffusive, 237 \ie the global variance of a tracer is not preserved using CEN4.241 \ie\ the global variance of a tracer is not preserved using CEN4. 238 242 Furthermore, it must be used in conjunction with an explicit diffusion operator to produce a sensible solution. 239 243 As in CEN2 case, the time-stepping is performed using a leapfrog scheme in conjunction with an Asselin time-filter, 240 so $T$ in (\autoref{eq: tra_adv_cen4}) is the \textit{now} tracer.244 so $T$ in (\autoref{eq:TRA_adv_cen4}) is the \textit{now} tracer. 241 245 242 246 At a $T$-grid cell adjacent to a boundary (coastline, bottom and surface), … … 248 252 249 253 % ------------------------------------------------------------------------------------------------------------- 250 % FCT scheme 251 % ------------------------------------------------------------------------------------------------------------- 252 \subsection{FCT: Flux Corrected Transport scheme (\protect\np{ln\_traadv\_fct}~\forcode{= .true.})} 254 % FCT scheme 255 % ------------------------------------------------------------------------------------------------------------- 256 \subsection[FCT: Flux Corrected Transport scheme (\forcode{ln_traadv_fct=.true.})] 257 {FCT: Flux Corrected Transport scheme (\protect\np{ln\_traadv\_fct}\forcode{=.true.})} 253 258 \label{subsec:TRA_adv_tvd} 254 259 255 The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct} ~\forcode{=.true.}.260 The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct}\forcode{=.true.}. 256 261 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) and vertical direction by 257 262 setting \np{nn\_fct\_h} and \np{nn\_fct\_v} to $2$ or $4$. … … 262 267 For example, in the $i$-direction : 263 268 \begin{equation} 264 \label{eq: tra_adv_fct}269 \label{eq:TRA_adv_fct} 265 270 \begin{split} 266 271 \tau_u^{ups} &= … … 275 280 where $c_u$ is a flux limiter function taking values between 0 and 1. 276 281 The FCT order is the one of the centred scheme used 277 (\ie it depends on the setting of \np{nn\_fct\_h} and \np{nn\_fct\_v}).282 (\ie\ it depends on the setting of \np{nn\_fct\_h} and \np{nn\_fct\_v}). 278 283 There exist many ways to define $c_u$, each corresponding to a different FCT scheme. 279 The one chosen in \NEMO is described in \citet{Zalesak_JCP79}.284 The one chosen in \NEMO\ is described in \citet{zalesak_JCP79}. 280 285 $c_u$ only departs from $1$ when the advective term produces a local extremum in the tracer field. 281 286 The resulting scheme is quite expensive but \textit{positive}. 282 287 It can be used on both active and passive tracers. 283 A comparison of FCT-2 with MUSCL and a MPDATA scheme can be found in \citet{Levy_al_GRL01}. 284 285 An additional option has been added controlled by \np{nn\_fct\_zts}. 286 By setting this integer to a value larger than zero, 287 a $2^{nd}$ order FCT scheme is used on both horizontal and vertical direction, but on the latter, 288 a split-explicit time stepping is used, with a number of sub-timestep equals to \np{nn\_fct\_zts}. 289 This option can be useful when the size of the timestep is limited by vertical advection \citep{Lemarie_OM2015}. 290 Note that in this case, a similar split-explicit time stepping should be used on vertical advection of momentum to 291 insure a better stability (see \autoref{subsec:DYN_zad}). 292 293 For stability reasons (see \autoref{chap:STP}), 294 $\tau_u^{cen}$ is evaluated in (\autoref{eq:tra_adv_fct}) using the \textit{now} tracer while 288 A comparison of FCT-2 with MUSCL and a MPDATA scheme can be found in \citet{levy.estublier.ea_GRL01}. 289 290 291 For stability reasons (see \autoref{chap:TD}), 292 $\tau_u^{cen}$ is evaluated in (\autoref{eq:TRA_adv_fct}) using the \textit{now} tracer while 295 293 $\tau_u^{ups}$ is evaluated using the \textit{before} tracer. 296 294 In other words, the advective part of the scheme is time stepped with a leap-frog scheme … … 298 296 299 297 % ------------------------------------------------------------------------------------------------------------- 300 % MUSCL scheme 301 % ------------------------------------------------------------------------------------------------------------- 302 \subsection{MUSCL: Monotone Upstream Scheme for Conservative Laws (\protect\np{ln\_traadv\_mus}~\forcode{= .true.})} 298 % MUSCL scheme 299 % ------------------------------------------------------------------------------------------------------------- 300 \subsection[MUSCL: Monotone Upstream Scheme for Conservative Laws (\forcode{ln_traadv_mus=.true.})] 301 {MUSCL: Monotone Upstream Scheme for Conservative Laws (\protect\np{ln\_traadv\_mus}\forcode{=.true.})} 303 302 \label{subsec:TRA_adv_mus} 304 303 305 The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln\_traadv\_mus} ~\forcode{=.true.}.304 The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln\_traadv\_mus}\forcode{=.true.}. 306 305 MUSCL implementation can be found in the \mdl{traadv\_mus} module. 307 306 308 MUSCL has been first implemented in \NEMO by \citet{Levy_al_GRL01}.307 MUSCL has been first implemented in \NEMO\ by \citet{levy.estublier.ea_GRL01}. 309 308 In its formulation, the tracer at velocity points is evaluated assuming a linear tracer variation between 310 two $T$-points (\autoref{fig: adv_scheme}).309 two $T$-points (\autoref{fig:TRA_adv_scheme}). 311 310 For example, in the $i$-direction : 312 311 \begin{equation} 313 % \label{eq: tra_adv_mus}312 % \label{eq:TRA_adv_mus} 314 313 \tau_u^{mus} = \lt\{ 315 314 \begin{split} … … 331 330 This choice ensure the \textit{positive} character of the scheme. 332 331 In addition, fluxes round a grid-point where a runoff is applied can optionally be computed using upstream fluxes 333 (\np{ln\_mus\_ups}~\forcode{= .true.}). 334 335 % ------------------------------------------------------------------------------------------------------------- 336 % UBS scheme 337 % ------------------------------------------------------------------------------------------------------------- 338 \subsection{UBS a.k.a. UP3: Upstream-Biased Scheme (\protect\np{ln\_traadv\_ubs}~\forcode{= .true.})} 332 (\np{ln\_mus\_ups}\forcode{=.true.}). 333 334 % ------------------------------------------------------------------------------------------------------------- 335 % UBS scheme 336 % ------------------------------------------------------------------------------------------------------------- 337 \subsection[UBS a.k.a. UP3: Upstream-Biased Scheme (\forcode{ln_traadv_ubs=.true.})] 338 {UBS a.k.a. UP3: Upstream-Biased Scheme (\protect\np{ln\_traadv\_ubs}\forcode{=.true.})} 339 339 \label{subsec:TRA_adv_ubs} 340 340 341 The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs} ~\forcode{=.true.}.341 The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs}\forcode{=.true.}. 342 342 UBS implementation can be found in the \mdl{traadv\_mus} module. 343 343 … … 347 347 For example, in the $i$-direction: 348 348 \begin{equation} 349 \label{eq: tra_adv_ubs}349 \label{eq:TRA_adv_ubs} 350 350 \tau_u^{ubs} = \overline T ^{i + 1/2} - \frac{1}{6} 351 351 \begin{cases} … … 358 358 359 359 This results in a dissipatively dominant (i.e. hyper-diffusive) truncation error 360 \citep{ Shchepetkin_McWilliams_OM05}.361 The overall performance of the advection scheme is similar to that reported in \cite{ Farrow1995}.360 \citep{shchepetkin.mcwilliams_OM05}. 361 The overall performance of the advection scheme is similar to that reported in \cite{farrow.stevens_JPO95}. 362 362 It is a relatively good compromise between accuracy and smoothness. 363 363 Nevertheless the scheme is not \textit{positive}, meaning that false extrema are permitted, … … 367 367 The intrinsic diffusion of UBS makes its use risky in the vertical direction where 368 368 the control of artificial diapycnal fluxes is of paramount importance 369 \citep{ Shchepetkin_McWilliams_OM05, Demange_PhD2014}.369 \citep{shchepetkin.mcwilliams_OM05, demange_phd14}. 370 370 Therefore the vertical flux is evaluated using either a $2^nd$ order FCT scheme or a $4^th$ order COMPACT scheme 371 (\np{nn\_ cen\_v}~\forcode{=2 or 4}).372 373 For stability reasons (see \autoref{chap: STP}), the first term in \autoref{eq:tra_adv_ubs}371 (\np{nn\_ubs\_v}\forcode{=2 or 4}). 372 373 For stability reasons (see \autoref{chap:TD}), the first term in \autoref{eq:TRA_adv_ubs} 374 374 (which corresponds to a second order centred scheme) 375 375 is evaluated using the \textit{now} tracer (centred in time) while the second term 376 376 (which is the diffusive part of the scheme), 377 377 is evaluated using the \textit{before} tracer (forward in time). 378 This choice is discussed by \citet{ Webb_al_JAOT98} in the context of the QUICK advection scheme.378 This choice is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the QUICK advection scheme. 379 379 UBS and QUICK schemes only differ by one coefficient. 380 Replacing 1/6 with 1/8 in \autoref{eq: tra_adv_ubs} leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.380 Replacing 1/6 with 1/8 in \autoref{eq:TRA_adv_ubs} leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}. 381 381 This option is not available through a namelist parameter, since the 1/6 coefficient is hard coded. 382 382 Nevertheless it is quite easy to make the substitution in the \mdl{traadv\_ubs} module and obtain a QUICK scheme. 383 383 384 Note that it is straightforward to rewrite \autoref{eq: tra_adv_ubs} as follows:384 Note that it is straightforward to rewrite \autoref{eq:TRA_adv_ubs} as follows: 385 385 \begin{gather} 386 \label{eq: traadv_ubs2}386 \label{eq:TRA_adv_ubs2} 387 387 \tau_u^{ubs} = \tau_u^{cen4} + \frac{1}{12} 388 388 \begin{cases} … … 391 391 \end{cases} 392 392 \intertext{or equivalently} 393 % \label{eq: traadv_ubs2b}393 % \label{eq:TRA_adv_ubs2b} 394 394 u_{i + 1/2} \ \tau_u^{ubs} = u_{i + 1/2} \, \overline{T - \frac{1}{6} \, \delta_i \Big[ \delta_{i + 1/2}[T] \Big]}^{\,i + 1/2} 395 395 - \frac{1}{2} |u|_{i + 1/2} \, \frac{1}{6} \, \delta_{i + 1/2} [\tau"_i] \nonumber 396 396 \end{gather} 397 397 398 \autoref{eq: traadv_ubs2} has several advantages.398 \autoref{eq:TRA_adv_ubs2} has several advantages. 399 399 Firstly, it clearly reveals that the UBS scheme is based on the fourth order scheme to which 400 400 an upstream-biased diffusion term is added. 401 401 Secondly, this emphasises that the $4^{th}$ order part (as well as the $2^{nd}$ order part as stated above) has to 402 be evaluated at the \textit{now} time step using \autoref{eq: tra_adv_ubs}.402 be evaluated at the \textit{now} time step using \autoref{eq:TRA_adv_ubs}. 403 403 Thirdly, the diffusion term is in fact a biharmonic operator with an eddy coefficient which 404 404 is simply proportional to the velocity: $A_u^{lm} = \frac{1}{12} \, {e_{1u}}^3 \, |u|$. 405 Note the current version of NEMO uses the computationally more efficient formulation \autoref{eq:tra_adv_ubs}. 406 407 % ------------------------------------------------------------------------------------------------------------- 408 % QCK scheme 409 % ------------------------------------------------------------------------------------------------------------- 410 \subsection{QCK: QuiCKest scheme (\protect\np{ln\_traadv\_qck}~\forcode{= .true.})} 405 Note the current version of \NEMO\ uses the computationally more efficient formulation \autoref{eq:TRA_adv_ubs}. 406 407 % ------------------------------------------------------------------------------------------------------------- 408 % QCK scheme 409 % ------------------------------------------------------------------------------------------------------------- 410 \subsection[QCK: QuiCKest scheme (\forcode{ln_traadv_qck=.true.})] 411 {QCK: QuiCKest scheme (\protect\np{ln\_traadv\_qck}\forcode{=.true.})} 411 412 \label{subsec:TRA_adv_qck} 412 413 413 414 The Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms (QUICKEST) scheme 414 proposed by \citet{ Leonard1979} is used when \np{ln\_traadv\_qck}~\forcode{=.true.}.415 proposed by \citet{leonard_CMAME79} is used when \np{ln\_traadv\_qck}\forcode{=.true.}. 415 416 QUICKEST implementation can be found in the \mdl{traadv\_qck} module. 416 417 417 418 QUICKEST is the third order Godunov scheme which is associated with the ULTIMATE QUICKEST limiter 418 \citep{ Leonard1991}.419 It has been implemented in NEMOby G. Reffray (MERCATOR-ocean) and can be found in the \mdl{traadv\_qck} module.419 \citep{leonard_CMAME91}. 420 It has been implemented in \NEMO\ by G. Reffray (MERCATOR-ocean) and can be found in the \mdl{traadv\_qck} module. 420 421 The resulting scheme is quite expensive but \textit{positive}. 421 422 It can be used on both active and passive tracers. … … 431 432 % Tracer Lateral Diffusion 432 433 % ================================================================ 433 \section{Tracer lateral diffusion (\protect\mdl{traldf})} 434 \section[Tracer lateral diffusion (\textit{traldf.F90})] 435 {Tracer lateral diffusion (\protect\mdl{traldf})} 434 436 \label{sec:TRA_ldf} 435 437 %-----------------------------------------nam_traldf------------------------------------------------------ 436 438 437 \nlst{namtra_ldf} 439 \begin{listing} 440 \nlst{namtra_ldf} 441 \caption{\texttt{namtra\_ldf}} 442 \label{lst:namtra_ldf} 443 \end{listing} 438 444 %------------------------------------------------------------------------------------------------------------- 439 440 Options are defined through the \n gn{namtra\_ldf} namelist variables.441 They are regrouped in four items, allowing to specify 445 446 Options are defined through the \nam{tra\_ldf} namelist variables. 447 They are regrouped in four items, allowing to specify 442 448 $(i)$ the type of operator used (none, laplacian, bilaplacian), 443 449 $(ii)$ the direction along which the operator acts (iso-level, horizontal, iso-neutral), 444 $(iii)$ some specific options related to the rotated operators (\ie non-iso-level operator), and450 $(iii)$ some specific options related to the rotated operators (\ie\ non-iso-level operator), and 445 451 $(iv)$ the specification of eddy diffusivity coefficient (either constant or variable in space and time). 446 452 Item $(iv)$ will be described in \autoref{chap:LDF}. … … 450 456 451 457 The lateral diffusion of tracers is evaluated using a forward scheme, 452 \ie the tracers appearing in its expression are the \textit{before} tracers in time,458 \ie\ the tracers appearing in its expression are the \textit{before} tracers in time, 453 459 except for the pure vertical component that appears when a rotation tensor is used. 454 This latter component is solved implicitly together with the vertical diffusion term (see \autoref{chap: STP}).455 When \np{ln\_traldf\_msc} ~\forcode{=.true.}, a Method of Stabilizing Correction is used in which456 the pure vertical component is split into an explicit and an implicit part \citep{ Lemarie_OM2012}.460 This latter component is solved implicitly together with the vertical diffusion term (see \autoref{chap:TD}). 461 When \np{ln\_traldf\_msc}\forcode{=.true.}, a Method of Stabilizing Correction is used in which 462 the pure vertical component is split into an explicit and an implicit part \citep{lemarie.debreu.ea_OM12}. 457 463 458 464 % ------------------------------------------------------------------------------------------------------------- 459 465 % Type of operator 460 466 % ------------------------------------------------------------------------------------------------------------- 461 \subsection[Type of operator (\protect\np{ln\_traldf}\{\_NONE,\_lap,\_blp\}\})]{Type of operator (\protect\np{ln\_traldf\_NONE}, \protect\np{ln\_traldf\_lap}, or \protect\np{ln\_traldf\_blp}) } 467 \subsection[Type of operator (\texttt{ln\_traldf}\{\texttt{\_OFF,\_lap,\_blp}\})] 468 {Type of operator (\protect\np{ln\_traldf\_OFF}, \protect\np{ln\_traldf\_lap}, or \protect\np{ln\_traldf\_blp}) } 462 469 \label{subsec:TRA_ldf_op} 463 470 … … 465 472 466 473 \begin{description} 467 \item[\np{ln\_traldf\_ NONE}~\forcode{=.true.}:]474 \item[\np{ln\_traldf\_OFF}\forcode{=.true.}:] 468 475 no operator selected, the lateral diffusive tendency will not be applied to the tracer equation. 469 476 This option can be used when the selected advection scheme is diffusive enough (MUSCL scheme for example). 470 \item[\np{ln\_traldf\_lap} ~\forcode{=.true.}:]477 \item[\np{ln\_traldf\_lap}\forcode{=.true.}:] 471 478 a laplacian operator is selected. 472 This harmonic operator takes the following expression: $\math pzc{L}(T) = \nabla \cdot A_{ht} \; \nabla T $,479 This harmonic operator takes the following expression: $\mathcal{L}(T) = \nabla \cdot A_{ht} \; \nabla T $, 473 480 where the gradient operates along the selected direction (see \autoref{subsec:TRA_ldf_dir}), 474 481 and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see \autoref{chap:LDF}). 475 \item[\np{ln\_traldf\_blp} ~\forcode{=.true.}]:482 \item[\np{ln\_traldf\_blp}\forcode{=.true.}]: 476 483 a bilaplacian operator is selected. 477 484 This biharmonic operator takes the following expression: 478 $\math pzc{B} = - \mathpzc{L}(\mathpzc{L}(T)) = - \nabla \cdot b \nabla (\nabla \cdot b \nabla T)$485 $\mathcal{B} = - \mathcal{L}(\mathcal{L}(T)) = - \nabla \cdot b \nabla (\nabla \cdot b \nabla T)$ 479 486 where the gradient operats along the selected direction, 480 487 and $b^2 = B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$ (see \autoref{chap:LDF}). … … 486 493 minimizing the impact on the larger scale features. 487 494 The main difference between the two operators is the scale selectiveness. 488 The bilaplacian damping time (\ie its spin down time) scales like $\lambda^{-4}$ for495 The bilaplacian damping time (\ie\ its spin down time) scales like $\lambda^{-4}$ for 489 496 disturbances of wavelength $\lambda$ (so that short waves damped more rapidelly than long ones), 490 497 whereas the laplacian damping time scales only like $\lambda^{-2}$. … … 493 500 % Direction of action 494 501 % ------------------------------------------------------------------------------------------------------------- 495 \subsection[Action direction (\protect\np{ln\_traldf}\{\_lev,\_hor,\_iso,\_triad\})]{Direction of action (\protect\np{ln\_traldf\_lev}, \protect\np{ln\_traldf\_hor}, \protect\np{ln\_traldf\_iso}, or \protect\np{ln\_traldf\_triad}) } 502 \subsection[Action direction (\texttt{ln\_traldf}\{\texttt{\_lev,\_hor,\_iso,\_triad}\})] 503 {Direction of action (\protect\np{ln\_traldf\_lev}, \protect\np{ln\_traldf\_hor}, \protect\np{ln\_traldf\_iso}, or \protect\np{ln\_traldf\_triad}) } 496 504 \label{subsec:TRA_ldf_dir} 497 505 498 506 The choice of a direction of action determines the form of operator used. 499 507 The operator is a simple (re-entrant) laplacian acting in the (\textbf{i},\textbf{j}) plane when 500 iso-level option is used (\np{ln\_traldf\_lev} ~\forcode{=.true.}) or501 when a horizontal (\ie geopotential) operator is demanded in \textit{z}-coordinate502 (\np{ln\_traldf\_hor} and \np{ln\_zco} equal \forcode{.true.}).508 iso-level option is used (\np{ln\_traldf\_lev}\forcode{=.true.}) or 509 when a horizontal (\ie\ geopotential) operator is demanded in \textit{z}-coordinate 510 (\np{ln\_traldf\_hor} and \np{ln\_zco}\forcode{=.true.}). 503 511 The associated code can be found in the \mdl{traldf\_lap\_blp} module. 504 512 The operator is a rotated (re-entrant) laplacian when 505 513 the direction along which it acts does not coincide with the iso-level surfaces, 506 514 that is when standard or triad iso-neutral option is used 507 (\np{ln\_traldf\_iso} or \np{ln\_traldf\_triad} equals\forcode{.true.},515 (\np{ln\_traldf\_iso} or \np{ln\_traldf\_triad} = \forcode{.true.}, 508 516 see \mdl{traldf\_iso} or \mdl{traldf\_triad} module, resp.), or 509 when a horizontal (\ie geopotential) operator is demanded in \textit{s}-coordinate510 (\np{ln\_traldf\_hor} and \np{ln\_sco} equal\forcode{.true.})517 when a horizontal (\ie\ geopotential) operator is demanded in \textit{s}-coordinate 518 (\np{ln\_traldf\_hor} and \np{ln\_sco} = \forcode{.true.}) 511 519 \footnote{In this case, the standard iso-neutral operator will be automatically selected}. 512 520 In that case, a rotation is applied to the gradient(s) that appears in the operator so that … … 519 527 % iso-level operator 520 528 % ------------------------------------------------------------------------------------------------------------- 521 \subsection{Iso-level (bi -)laplacian operator ( \protect\np{ln\_traldf\_iso}) } 529 \subsection[Iso-level (bi-)laplacian operator (\texttt{ln\_traldf\_iso})] 530 {Iso-level (bi-)laplacian operator ( \protect\np{ln\_traldf\_iso})} 522 531 \label{subsec:TRA_ldf_lev} 523 532 524 The laplacian diffusion operator acting along the model (\textit{i,j})-surfaces is given by: 525 \begin{equation} 526 \label{eq: tra_ldf_lap}533 The laplacian diffusion operator acting along the model (\textit{i,j})-surfaces is given by: 534 \begin{equation} 535 \label{eq:TRA_ldf_lap} 527 536 D_t^{lT} = \frac{1}{b_t} \Bigg( \delta_{i} \lt[ A_u^{lT} \; \frac{e_{2u} \, e_{3u}}{e_{1u}} \; \delta_{i + 1/2} [T] \rt] 528 537 + \delta_{j} \lt[ A_v^{lT} \; \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [T] \rt] \Bigg) … … 531 540 where zero diffusive fluxes is assumed across solid boundaries, 532 541 first (and third in bilaplacian case) horizontal tracer derivative are masked. 533 It is implemented in the \rou{tra ldf\_lap} subroutine found in the \mdl{traldf\_lap} module.534 The module also contains \rou{tra ldf\_blp}, the subroutine calling twice \rou{traldf\_lap} in order to542 It is implemented in the \rou{tra\_ldf\_lap} subroutine found in the \mdl{traldf\_lap\_blp} module. 543 The module also contains \rou{tra\_ldf\_blp}, the subroutine calling twice \rou{tra\_ldf\_lap} in order to 535 544 compute the iso-level bilaplacian operator. 536 545 537 546 It is a \textit{horizontal} operator (\ie acting along geopotential surfaces) in 538 547 the $z$-coordinate with or without partial steps, but is simply an iso-level operator in the $s$-coordinate. 539 It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp} ~\forcode{=.true.},540 we have \np{ln\_traldf\_lev} ~\forcode{= .true.} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}~\forcode{=.true.}.548 It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp}\forcode{=.true.}, 549 we have \np{ln\_traldf\_lev}\forcode{=.true.} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}\forcode{=.true.}. 541 550 In both cases, it significantly contributes to diapycnal mixing. 542 551 It is therefore never recommended, even when using it in the bilaplacian case. 543 552 544 Note that in the partial step $z$-coordinate (\np{ln\_zps} ~\forcode{=.true.}),553 Note that in the partial step $z$-coordinate (\np{ln\_zps}\forcode{=.true.}), 545 554 tracers in horizontally adjacent cells are located at different depths in the vicinity of the bottom. 546 In this case, horizontal derivatives in (\autoref{eq: tra_ldf_lap}) at the bottom level require a specific treatment.555 In this case, horizontal derivatives in (\autoref{eq:TRA_ldf_lap}) at the bottom level require a specific treatment. 547 556 They are calculated in the \mdl{zpshde} module, described in \autoref{sec:TRA_zpshde}. 548 557 … … 550 559 % Rotated laplacian operator 551 560 % ------------------------------------------------------------------------------------------------------------- 552 \subsection{Standard and triad (bi 561 \subsection{Standard and triad (bi-)laplacian operator} 553 562 \label{subsec:TRA_ldf_iso_triad} 554 563 555 %&& Standard rotated (bi 564 %&& Standard rotated (bi-)laplacian operator 556 565 %&& ---------------------------------------------- 557 \subsubsection{Standard rotated (bi -)laplacian operator (\protect\mdl{traldf\_iso})} 566 \subsubsection[Standard rotated (bi-)laplacian operator (\textit{traldf\_iso.F90})] 567 {Standard rotated (bi-)laplacian operator (\protect\mdl{traldf\_iso})} 558 568 \label{subsec:TRA_ldf_iso} 559 The general form of the second order lateral tracer subgrid scale physics (\autoref{eq: PE_zdf})569 The general form of the second order lateral tracer subgrid scale physics (\autoref{eq:MB_zdf}) 560 570 takes the following semi -discrete space form in $z$- and $s$-coordinates: 561 571 \begin{equation} 562 \label{eq: tra_ldf_iso}572 \label{eq:TRA_ldf_iso} 563 573 \begin{split} 564 574 D_T^{lT} = \frac{1}{b_t} \Bigg[ \quad &\delta_i A_u^{lT} \lt( \frac{e_{2u} e_{3u}}{e_{1u}} \, \delta_{i + 1/2} [T] … … 573 583 where $b_t = e_{1t} \, e_{2t} \, e_{3t}$ is the volume of $T$-cells, 574 584 $r_1$ and $r_2$ are the slopes between the surface of computation ($z$- or $s$-surfaces) and 575 the surface along which the diffusion operator acts (\ie horizontal or iso-neutral surfaces).576 It is thus used when, in addition to \np{ln\_traldf\_lap} ~\forcode{=.true.},577 we have \np{ln\_traldf\_iso} ~\forcode{=.true.},578 or both \np{ln\_traldf\_hor} ~\forcode{= .true.} and \np{ln\_zco}~\forcode{=.true.}.585 the surface along which the diffusion operator acts (\ie\ horizontal or iso-neutral surfaces). 586 It is thus used when, in addition to \np{ln\_traldf\_lap}\forcode{=.true.}, 587 we have \np{ln\_traldf\_iso}\forcode{=.true.}, 588 or both \np{ln\_traldf\_hor}\forcode{=.true.} and \np{ln\_zco}\forcode{=.true.}. 579 589 The way these slopes are evaluated is given in \autoref{sec:LDF_slp}. 580 590 At the surface, bottom and lateral boundaries, the turbulent fluxes of heat and salt are set to zero using 581 591 the mask technique (see \autoref{sec:LBC_coast}). 582 592 583 The operator in \autoref{eq: tra_ldf_iso} involves both lateral and vertical derivatives.593 The operator in \autoref{eq:TRA_ldf_iso} involves both lateral and vertical derivatives. 584 594 For numerical stability, the vertical second derivative must be solved using the same implicit time scheme as that 585 595 used in the vertical physics (see \autoref{sec:TRA_zdf}). … … 590 600 This formulation conserves the tracer but does not ensure the decrease of the tracer variance. 591 601 Nevertheless the treatment performed on the slopes (see \autoref{chap:LDF}) allows the model to run safely without 592 any additional background horizontal diffusion \citep{ Guilyardi_al_CD01}.593 594 Note that in the partial step $z$-coordinate (\np{ln\_zps} ~\forcode{=.true.}),595 the horizontal derivatives at the bottom level in \autoref{eq: tra_ldf_iso} require a specific treatment.602 any additional background horizontal diffusion \citep{guilyardi.madec.ea_CD01}. 603 604 Note that in the partial step $z$-coordinate (\np{ln\_zps}\forcode{=.true.}), 605 the horizontal derivatives at the bottom level in \autoref{eq:TRA_ldf_iso} require a specific treatment. 596 606 They are calculated in module zpshde, described in \autoref{sec:TRA_zpshde}. 597 607 598 %&& Triad rotated (bi 608 %&& Triad rotated (bi-)laplacian operator 599 609 %&& ------------------------------------------- 600 \subsubsection{Triad rotated (bi -)laplacian operator (\protect\np{ln\_traldf\_triad})} 610 \subsubsection[Triad rotated (bi-)laplacian operator (\textit{ln\_traldf\_triad})] 611 {Triad rotated (bi-)laplacian operator (\protect\np{ln\_traldf\_triad})} 601 612 \label{subsec:TRA_ldf_triad} 602 613 603 If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}~\forcode{= .true.}; see \autoref{apdx:triad}) 604 605 An alternative scheme developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases 606 is also available in \NEMO (\np{ln\_traldf\_grif}~\forcode{= .true.}). 607 A complete description of the algorithm is given in \autoref{apdx:triad}. 608 609 The lateral fourth order bilaplacian operator on tracers is obtained by applying (\autoref{eq:tra_ldf_lap}) twice. 614 An alternative scheme developed by \cite{griffies.gnanadesikan.ea_JPO98} which ensures tracer variance decreases 615 is also available in \NEMO\ (\np{ln\_traldf\_triad}\forcode{=.true.}). 616 A complete description of the algorithm is given in \autoref{apdx:TRIADS}. 617 618 The lateral fourth order bilaplacian operator on tracers is obtained by applying (\autoref{eq:TRA_ldf_lap}) twice. 610 619 The operator requires an additional assumption on boundary conditions: 611 620 both first and third derivative terms normal to the coast are set to zero. 612 621 613 The lateral fourth order operator formulation on tracers is obtained by applying (\autoref{eq: tra_ldf_iso}) twice.622 The lateral fourth order operator formulation on tracers is obtained by applying (\autoref{eq:TRA_ldf_iso}) twice. 614 623 It requires an additional assumption on boundary conditions: 615 624 first and third derivative terms normal to the coast, … … 625 634 \item \np{rn\_slpmax} = slope limit (both operators) 626 635 \item \np{ln\_triad\_iso} = pure horizontal mixing in ML (triad only) 627 \item \np{rn\_sw\_triad} $= 1$ switching triad; $= 0$ all 4 triads used (triad only) 636 \item \np{rn\_sw\_triad} $= 1$ switching triad; $= 0$ all 4 triads used (triad only) 628 637 \item \np{ln\_botmix\_triad} = lateral mixing on bottom (triad only) 629 638 \end{itemize} … … 632 641 % Tracer Vertical Diffusion 633 642 % ================================================================ 634 \section{Tracer vertical diffusion (\protect\mdl{trazdf})} 643 \section[Tracer vertical diffusion (\textit{trazdf.F90})] 644 {Tracer vertical diffusion (\protect\mdl{trazdf})} 635 645 \label{sec:TRA_zdf} 636 646 %--------------------------------------------namzdf--------------------------------------------------------- 637 647 638 \nlst{namzdf}639 648 %-------------------------------------------------------------------------------------------------------------- 640 649 641 Options are defined through the \n gn{namzdf} namelist variables.650 Options are defined through the \nam{zdf} namelist variables. 642 651 The formulation of the vertical subgrid scale tracer physics is the same for all the vertical coordinates, 643 652 and is based on a laplacian operator. 644 The vertical diffusion operator given by (\autoref{eq: PE_zdf}) takes the following semi -discrete space form:653 The vertical diffusion operator given by (\autoref{eq:MB_zdf}) takes the following semi -discrete space form: 645 654 \begin{gather*} 646 % \label{eq: tra_zdf}655 % \label{eq:TRA_zdf} 647 656 D^{vT}_T = \frac{1}{e_{3t}} \, \delta_k \lt[ \, \frac{A^{vT}_w}{e_{3w}} \delta_{k + 1/2}[T] \, \rt] \\ 648 657 D^{vS}_T = \frac{1}{e_{3t}} \; \delta_k \lt[ \, \frac{A^{vS}_w}{e_{3w}} \delta_{k + 1/2}[S] \, \rt] … … 651 660 respectively. 652 661 Generally, $A_w^{vT} = A_w^{vS}$ except when double diffusive mixing is parameterised 653 (\ie \key{zdfddm} is defined).662 (\ie\ \np{ln\_zdfddm}\forcode{=.true.},). 654 663 The way these coefficients are evaluated is given in \autoref{chap:ZDF} (ZDF). 655 664 Furthermore, when iso-neutral mixing is used, both mixing coefficients are increased by 656 665 $\frac{e_{1w} e_{2w}}{e_{3w} }({r_{1w}^2 + r_{2w}^2})$ to account for the vertical second derivative of 657 \autoref{eq: tra_ldf_iso}.666 \autoref{eq:TRA_ldf_iso}. 658 667 659 668 At the surface and bottom boundaries, the turbulent fluxes of heat and salt must be specified. … … 663 672 664 673 The large eddy coefficient found in the mixed layer together with high vertical resolution implies that 665 in the case of explicit time stepping (\np{ln\_zdfexp}~\forcode{= .true.}) 666 there would be too restrictive a constraint on the time step. 667 Therefore, the default implicit time stepping is preferred for the vertical diffusion since 674 there would be too restrictive constraint on the time step if we use explicit time stepping. 675 Therefore an implicit time stepping is preferred for the vertical diffusion since 668 676 it overcomes the stability constraint. 669 A forward time differencing scheme (\np{ln\_zdfexp}~\forcode{= .true.}) using670 a time splitting technique (\np{nn\_zdfexp} $> 1$) is provided as an alternative.671 Namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics.672 677 673 678 % ================================================================ … … 680 685 % surface boundary condition 681 686 % ------------------------------------------------------------------------------------------------------------- 682 \subsection{Surface boundary condition (\protect\mdl{trasbc})} 687 \subsection[Surface boundary condition (\textit{trasbc.F90})] 688 {Surface boundary condition (\protect\mdl{trasbc})} 683 689 \label{subsec:TRA_sbc} 684 690 … … 690 696 691 697 Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components 692 (\ie atmosphere, sea-ice, land), the change in the heat and salt content of the surface layer of the ocean is due698 (\ie\ atmosphere, sea-ice, land), the change in the heat and salt content of the surface layer of the ocean is due 693 699 both to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$) and 694 700 to the heat and salt content of the mass exchange. … … 702 708 \item 703 709 $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface 704 (\ie the difference between the total surface heat flux and the fraction of the short wave flux that710 (\ie\ the difference between the total surface heat flux and the fraction of the short wave flux that 705 711 penetrates into the water column, see \autoref{subsec:TRA_qsr}) 706 712 plus the heat content associated with of the mass exchange with the atmosphere and lands. … … 720 726 The surface boundary condition on temperature and salinity is applied as follows: 721 727 \begin{equation} 722 \label{eq: tra_sbc}728 \label{eq:TRA_sbc} 723 729 \begin{alignedat}{2} 724 730 F^T &= \frac{1}{C_p} &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} &\overline{Q_{ns} }^t \\ … … 728 734 where $\overline x^t$ means that $x$ is averaged over two consecutive time steps 729 735 ($t - \rdt / 2$ and $t + \rdt / 2$). 730 Such time averaging prevents the divergence of odd and even time step (see \autoref{chap: STP}).731 732 In the linear free surface case (\np{ln\_linssh} ~\forcode{=.true.}), an additional term has to be added on736 Such time averaging prevents the divergence of odd and even time step (see \autoref{chap:TD}). 737 738 In the linear free surface case (\np{ln\_linssh}\forcode{=.true.}), an additional term has to be added on 733 739 both temperature and salinity. 734 740 On temperature, this term remove the heat content associated with mass exchange that has been added to $Q_{ns}$. … … 737 743 The resulting surface boundary condition is applied as follows: 738 744 \begin{equation} 739 \label{eq: tra_sbc_lin}745 \label{eq:TRA_sbc_lin} 740 746 \begin{alignedat}{2} 741 747 F^T &= \frac{1}{C_p} &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} … … 744 750 &\overline{(\textit{sfx} - \textit{emp} \lt. S \rt|_{k = 1})}^t 745 751 \end{alignedat} 746 \end{equation} 752 \end{equation} 747 753 Note that an exact conservation of heat and salt content is only achieved with non-linear free surface. 748 754 In the linear free surface case, there is a small imbalance. 749 The imbalance is larger than the imbalance associated with the Asselin time filter \citep{Leclair_Madec_OM09}. 750 This is the reason why the modified filter is not applied in the linear free surface case (see \autoref{chap:STP}). 751 752 % ------------------------------------------------------------------------------------------------------------- 753 % Solar Radiation Penetration 754 % ------------------------------------------------------------------------------------------------------------- 755 \subsection{Solar radiation penetration (\protect\mdl{traqsr})} 755 The imbalance is larger than the imbalance associated with the Asselin time filter \citep{leclair.madec_OM09}. 756 This is the reason why the modified filter is not applied in the linear free surface case (see \autoref{chap:TD}). 757 758 % ------------------------------------------------------------------------------------------------------------- 759 % Solar Radiation Penetration 760 % ------------------------------------------------------------------------------------------------------------- 761 \subsection[Solar radiation penetration (\textit{traqsr.F90})] 762 {Solar radiation penetration (\protect\mdl{traqsr})} 756 763 \label{subsec:TRA_qsr} 757 764 %--------------------------------------------namqsr-------------------------------------------------------- 758 765 759 \nlst{namtra_qsr} 766 \begin{listing} 767 \nlst{namtra_qsr} 768 \caption{\texttt{namtra\_qsr}} 769 \label{lst:namtra_qsr} 770 \end{listing} 760 771 %-------------------------------------------------------------------------------------------------------------- 761 772 762 Options are defined through the \n gn{namtra\_qsr} namelist variables.763 When the penetrative solar radiation option is used (\np{ln\_ flxqsr}~\forcode{=.true.}),773 Options are defined through the \nam{tra\_qsr} namelist variables. 774 When the penetrative solar radiation option is used (\np{ln\_traqsr}\forcode{=.true.}), 764 775 the solar radiation penetrates the top few tens of meters of the ocean. 765 If it is not used (\np{ln\_ flxqsr}~\forcode{=.false.}) all the heat flux is absorbed in the first ocean level.766 Thus, in the former case a term is added to the time evolution equation of temperature \autoref{eq: PE_tra_T} and767 the surface boundary condition is modified to take into account only the non-penetrative part of the surface 776 If it is not used (\np{ln\_traqsr}\forcode{=.false.}) all the heat flux is absorbed in the first ocean level. 777 Thus, in the former case a term is added to the time evolution equation of temperature \autoref{eq:MB_PE_tra_T} and 778 the surface boundary condition is modified to take into account only the non-penetrative part of the surface 768 779 heat flux: 769 780 \begin{equation} 770 \label{eq: PE_qsr}781 \label{eq:TRA_PE_qsr} 771 782 \begin{gathered} 772 783 \pd[T]{t} = \ldots + \frac{1}{\rho_o \, C_p \, e_3} \; \pd[I]{k} \\ … … 774 785 \end{gathered} 775 786 \end{equation} 776 where $Q_{sr}$ is the penetrative part of the surface heat flux (\ie the shortwave radiation) and787 where $Q_{sr}$ is the penetrative part of the surface heat flux (\ie\ the shortwave radiation) and 777 788 $I$ is the downward irradiance ($\lt. I \rt|_{z = \eta} = Q_{sr}$). 778 The additional term in \autoref{eq: PE_qsr} is discretized as follows:779 \begin{equation} 780 \label{eq: tra_qsr}789 The additional term in \autoref{eq:TRA_PE_qsr} is discretized as follows: 790 \begin{equation} 791 \label{eq:TRA_qsr} 781 792 \frac{1}{\rho_o \, C_p \, e_3} \, \pd[I]{k} \equiv \frac{1}{\rho_o \, C_p \, e_{3t}} \delta_k [I_w] 782 793 \end{equation} … … 788 799 (specified through namelist parameter \np{rn\_abs}). 789 800 It is assumed to penetrate the ocean with a decreasing exponential profile, with an e-folding depth scale, $\xi_0$, 790 of a few tens of centimetres (typically $\xi_0 = 0.35~m$ set as \np{rn\_si0} in the \n gn{namtra\_qsr} namelist).801 of a few tens of centimetres (typically $\xi_0 = 0.35~m$ set as \np{rn\_si0} in the \nam{tra\_qsr} namelist). 791 802 For shorter wavelengths (400-700~nm), the ocean is more transparent, and solar energy propagates to 792 803 larger depths where it contributes to local heating. 793 804 The way this second part of the solar energy penetrates into the ocean depends on which formulation is chosen. 794 In the simple 2-waveband light penetration scheme (\np{ln\_qsr\_2bd} ~\forcode{=.true.})805 In the simple 2-waveband light penetration scheme (\np{ln\_qsr\_2bd}\forcode{=.true.}) 795 806 a chlorophyll-independent monochromatic formulation is chosen for the shorter wavelengths, 796 leading to the following expression \citep{ Paulson1977}:807 leading to the following expression \citep{paulson.simpson_JPO77}: 797 808 \[ 798 % \label{eq: traqsr_iradiance}809 % \label{eq:TRA_qsr_iradiance} 799 810 I(z) = Q_{sr} \lt[ Re^{- z / \xi_0} + (1 - R) e^{- z / \xi_1} \rt] 800 811 \] … … 805 816 806 817 Such assumptions have been shown to provide a very crude and simplistic representation of 807 observed light penetration profiles (\cite{ Morel_JGR88}, see also \autoref{fig:traqsr_irradiance}).818 observed light penetration profiles (\cite{morel_JGR88}, see also \autoref{fig:TRA_qsr_irradiance}). 808 819 Light absorption in the ocean depends on particle concentration and is spectrally selective. 809 \cite{ Morel_JGR88} has shown that an accurate representation of light penetration can be provided by820 \cite{morel_JGR88} has shown that an accurate representation of light penetration can be provided by 810 821 a 61 waveband formulation. 811 822 Unfortunately, such a model is very computationally expensive. 812 Thus, \cite{ Lengaigne_al_CD07} have constructed a simplified version of this formulation in which823 Thus, \cite{lengaigne.menkes.ea_CD07} have constructed a simplified version of this formulation in which 813 824 visible light is split into three wavebands: blue (400-500 nm), green (500-600 nm) and red (600-700nm). 814 825 For each wave-band, the chlorophyll-dependent attenuation coefficient is fitted to the coefficients computed from 815 the full spectral model of \cite{ Morel_JGR88} (as modified by \cite{Morel_Maritorena_JGR01}),826 the full spectral model of \cite{morel_JGR88} (as modified by \cite{morel.maritorena_JGR01}), 816 827 assuming the same power-law relationship. 817 As shown in \autoref{fig: traqsr_irradiance}, this formulation, called RGB (Red-Green-Blue),828 As shown in \autoref{fig:TRA_qsr_irradiance}, this formulation, called RGB (Red-Green-Blue), 818 829 reproduces quite closely the light penetration profiles predicted by the full spectal model, 819 830 but with much greater computational efficiency. 820 831 The 2-bands formulation does not reproduce the full model very well. 821 832 822 The RGB formulation is used when \np{ln\_qsr\_rgb} ~\forcode{=.true.}.823 The RGB attenuation coefficients (\ie the inverses of the extinction length scales) are tabulated over833 The RGB formulation is used when \np{ln\_qsr\_rgb}\forcode{=.true.}. 834 The RGB attenuation coefficients (\ie\ the inverses of the extinction length scales) are tabulated over 824 835 61 nonuniform chlorophyll classes ranging from 0.01 to 10 g.Chl/L 825 836 (see the routine \rou{trc\_oce\_rgb} in \mdl{trc\_oce} module). … … 827 838 828 839 \begin{description} 829 \item[\np{nn\_ch dta}~\forcode{=0}]830 a constant 0.05 g.Chl/L value everywhere ; 831 \item[\np{nn\_ch dta}~\forcode{=1}]840 \item[\np{nn\_chldta}\forcode{=0}] 841 a constant 0.05 g.Chl/L value everywhere ; 842 \item[\np{nn\_chldta}\forcode{=1}] 832 843 an observed time varying chlorophyll deduced from satellite surface ocean color measurement spread uniformly in 833 844 the vertical direction; 834 \item[\np{nn\_ch dta}~\forcode{=2}]845 \item[\np{nn\_chldta}\forcode{=2}] 835 846 same as previous case except that a vertical profile of chlorophyl is used. 836 Following \cite{ Morel_Berthon_LO89}, the profile is computed from the local surface chlorophyll value;837 \item[\np{ln\_qsr\_bio} ~\forcode{=.true.}]847 Following \cite{morel.berthon_LO89}, the profile is computed from the local surface chlorophyll value; 848 \item[\np{ln\_qsr\_bio}\forcode{=.true.}] 838 849 simulated time varying chlorophyll by TOP biogeochemical model. 839 850 In this case, the RGB formulation is used to calculate both the phytoplankton light limitation in 840 PISCES or LOBSTERand the oceanic heating rate.841 \end{description} 842 843 The trend in \autoref{eq: tra_qsr} associated with the penetration of the solar radiation is added to851 PISCES and the oceanic heating rate. 852 \end{description} 853 854 The trend in \autoref{eq:TRA_qsr} associated with the penetration of the solar radiation is added to 844 855 the temperature trend, and the surface heat flux is modified in routine \mdl{traqsr}. 845 856 … … 847 858 the depth of $w-$levels does not significantly vary with location. 848 859 The level at which the light has been totally absorbed 849 (\ie it is less than the computer precision) is computed once,860 (\ie\ it is less than the computer precision) is computed once, 850 861 and the trend associated with the penetration of the solar radiation is only added down to that level. 851 862 Finally, note that when the ocean is shallow ($<$ 200~m), part of the solar radiation can reach the ocean floor. 852 863 In this case, we have chosen that all remaining radiation is absorbed in the last ocean level 853 (\ie $I$ is masked).864 (\ie\ $I$ is masked). 854 865 855 866 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 856 867 \begin{figure}[!t] 857 \begin{center} 858 \includegraphics[]{Fig_TRA_Irradiance} 859 \caption{ 860 \protect\label{fig:traqsr_irradiance} 861 Penetration profile of the downward solar irradiance calculated by four models. 862 Two waveband chlorophyll-independent formulation (blue), 863 a chlorophyll-dependent monochromatic formulation (green), 864 4 waveband RGB formulation (red), 865 61 waveband Morel (1988) formulation (black) for a chlorophyll concentration of 866 (a) Chl=0.05 mg/m$^3$ and (b) Chl=0.5 mg/m$^3$. 867 From \citet{Lengaigne_al_CD07}. 868 } 869 \end{center} 868 \centering 869 \includegraphics[width=\textwidth]{Fig_TRA_Irradiance} 870 \caption[Penetration profile of the downward solar irradiance calculated by four models]{ 871 Penetration profile of the downward solar irradiance calculated by four models. 872 Two waveband chlorophyll-independent formulation (blue), 873 a chlorophyll-dependent monochromatic formulation (green), 874 4 waveband RGB formulation (red), 875 61 waveband Morel (1988) formulation (black) for a chlorophyll concentration of 876 (a) Chl=0.05 mg/m$^3$ and (b) Chl=0.5 mg/m$^3$. 877 From \citet{lengaigne.menkes.ea_CD07}.} 878 \label{fig:TRA_qsr_irradiance} 870 879 \end{figure} 871 880 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 874 883 % Bottom Boundary Condition 875 884 % ------------------------------------------------------------------------------------------------------------- 876 \subsection{Bottom boundary condition (\protect\mdl{trabbc})} 885 \subsection[Bottom boundary condition (\textit{trabbc.F90}) - \forcode{ln_trabbc=.true.})] 886 {Bottom boundary condition (\protect\mdl{trabbc})} 877 887 \label{subsec:TRA_bbc} 878 888 %--------------------------------------------nambbc-------------------------------------------------------- 879 889 880 \nlst{nambbc} 890 \begin{listing} 891 \nlst{nambbc} 892 \caption{\texttt{nambbc}} 893 \label{lst:nambbc} 894 \end{listing} 881 895 %-------------------------------------------------------------------------------------------------------------- 882 896 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 883 897 \begin{figure}[!t] 884 \begin{center} 885 \includegraphics[]{Fig_TRA_geoth} 886 \caption{ 887 \protect\label{fig:geothermal} 888 Geothermal Heat flux (in $mW.m^{-2}$) used by \cite{Emile-Geay_Madec_OS09}. 889 It is inferred from the age of the sea floor and the formulae of \citet{Stein_Stein_Nat92}. 890 } 891 \end{center} 898 \centering 899 \includegraphics[width=\textwidth]{Fig_TRA_geoth} 900 \caption[Geothermal heat flux]{ 901 Geothermal Heat flux (in $mW.m^{-2}$) used by \cite{emile-geay.madec_OS09}. 902 It is inferred from the age of the sea floor and the formulae of \citet{stein.stein_N92}.} 903 \label{fig:TRA_geothermal} 892 904 \end{figure} 893 905 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 894 906 895 907 Usually it is assumed that there is no exchange of heat or salt through the ocean bottom, 896 \ie a no flux boundary condition is applied on active tracers at the bottom.908 \ie\ a no flux boundary condition is applied on active tracers at the bottom. 897 909 This is the default option in \NEMO, and it is implemented using the masking technique. 898 910 However, there is a non-zero heat flux across the seafloor that is associated with solid earth cooling. 899 This flux is weak compared to surface fluxes (a mean global value of $\sim 0.1 \, W/m^2$ \citep{ Stein_Stein_Nat92}),911 This flux is weak compared to surface fluxes (a mean global value of $\sim 0.1 \, W/m^2$ \citep{stein.stein_N92}), 900 912 but it warms systematically the ocean and acts on the densest water masses. 901 913 Taking this flux into account in a global ocean model increases the deepest overturning cell 902 (\ie the one associated with the Antarctic Bottom Water) by a few Sverdrups \citep{Emile-Geay_Madec_OS09}.903 904 Options are defined through the \ngn{namtra\_bbc} namelist variables.914 (\ie\ the one associated with the Antarctic Bottom Water) by a few Sverdrups \citep{emile-geay.madec_OS09}. 915 916 Options are defined through the \nam{bbc} namelist variables. 905 917 The presence of geothermal heating is controlled by setting the namelist parameter \np{ln\_trabbc} to true. 906 918 Then, when \np{nn\_geoflx} is set to 1, a constant geothermal heating is introduced whose value is given by 907 the \np{ nn\_geoflx\_cst}, which is also a namelist parameter.919 the \np{rn\_geoflx\_cst}, which is also a namelist parameter. 908 920 When \np{nn\_geoflx} is set to 2, a spatially varying geothermal heat flux is introduced which is provided in 909 the \ifile{geothermal\_heating} NetCDF file (\autoref{fig: geothermal}) \citep{Emile-Geay_Madec_OS09}.921 the \ifile{geothermal\_heating} NetCDF file (\autoref{fig:TRA_geothermal}) \citep{emile-geay.madec_OS09}. 910 922 911 923 % ================================================================ 912 924 % Bottom Boundary Layer 913 925 % ================================================================ 914 \section{Bottom boundary layer (\protect\mdl{trabbl} - \protect\key{trabbl})} 926 \section[Bottom boundary layer (\textit{trabbl.F90} - \forcode{ln_trabbl=.true.})] 927 {Bottom boundary layer (\protect\mdl{trabbl} - \protect\np{ln\_trabbl}\forcode{=.true.})} 915 928 \label{sec:TRA_bbl} 916 929 %--------------------------------------------nambbl--------------------------------------------------------- 917 930 918 \nlst{nambbl} 931 \begin{listing} 932 \nlst{nambbl} 933 \caption{\texttt{nambbl}} 934 \label{lst:nambbl} 935 \end{listing} 919 936 %-------------------------------------------------------------------------------------------------------------- 920 937 921 Options are defined through the \n gn{nambbl} namelist variables.938 Options are defined through the \nam{bbl} namelist variables. 922 939 In a $z$-coordinate configuration, the bottom topography is represented by a series of discrete steps. 923 940 This is not adequate to represent gravity driven downslope flows. … … 931 948 sometimes over a thickness much larger than the thickness of the observed gravity plume. 932 949 A similar problem occurs in the $s$-coordinate when the thickness of the bottom level varies rapidly downstream of 933 a sill \citep{ Willebrand_al_PO01}, and the thickness of the plume is not resolved.934 935 The idea of the bottom boundary layer (BBL) parameterisation, first introduced by \citet{ Beckmann_Doscher1997},950 a sill \citep{willebrand.barnier.ea_PO01}, and the thickness of the plume is not resolved. 951 952 The idea of the bottom boundary layer (BBL) parameterisation, first introduced by \citet{beckmann.doscher_JPO97}, 936 953 is to allow a direct communication between two adjacent bottom cells at different levels, 937 954 whenever the densest water is located above the less dense water. … … 939 956 In the current implementation of the BBL, only the tracers are modified, not the velocities. 940 957 Furthermore, it only connects ocean bottom cells, and therefore does not include all the improvements introduced by 941 \citet{ Campin_Goosse_Tel99}.958 \citet{campin.goosse_T99}. 942 959 943 960 % ------------------------------------------------------------------------------------------------------------- 944 961 % Diffusive BBL 945 962 % ------------------------------------------------------------------------------------------------------------- 946 \subsection{Diffusive bottom boundary layer (\protect\np{nn\_bbl\_ldf}~\forcode{= 1})} 963 \subsection[Diffusive bottom boundary layer (\forcode{nn_bbl_ldf=1})] 964 {Diffusive bottom boundary layer (\protect\np{nn\_bbl\_ldf}\forcode{=1})} 947 965 \label{subsec:TRA_bbl_diff} 948 966 949 When applying sigma-diffusion (\ key{trabbl} definedand \np{nn\_bbl\_ldf} set to 1),950 the diffusive flux between two adjacent cells at the ocean floor is given by 967 When applying sigma-diffusion (\np{ln\_trabbl}\forcode{=.true.} and \np{nn\_bbl\_ldf} set to 1), 968 the diffusive flux between two adjacent cells at the ocean floor is given by 951 969 \[ 952 % \label{eq: tra_bbl_diff}970 % \label{eq:TRA_bbl_diff} 953 971 \vect F_\sigma = A_l^\sigma \, \nabla_\sigma T 954 972 \] 955 973 with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, and 956 974 $A_l^\sigma$ the lateral diffusivity in the BBL. 957 Following \citet{ Beckmann_Doscher1997}, the latter is prescribed with a spatial dependence,958 \ie in the conditional form959 \begin{equation} 960 \label{eq: tra_bbl_coef}975 Following \citet{beckmann.doscher_JPO97}, the latter is prescribed with a spatial dependence, 976 \ie\ in the conditional form 977 \begin{equation} 978 \label{eq:TRA_bbl_coef} 961 979 A_l^\sigma (i,j,t) = 962 980 \begin{cases} … … 968 986 where $A_{bbl}$ is the BBL diffusivity coefficient, given by the namelist parameter \np{rn\_ahtbbl} and 969 987 usually set to a value much larger than the one used for lateral mixing in the open ocean. 970 The constraint in \autoref{eq: tra_bbl_coef} implies that sigma-like diffusion only occurs when988 The constraint in \autoref{eq:TRA_bbl_coef} implies that sigma-like diffusion only occurs when 971 989 the density above the sea floor, at the top of the slope, is larger than in the deeper ocean 972 (see green arrow in \autoref{fig: bbl}).990 (see green arrow in \autoref{fig:TRA_bbl}). 973 991 In practice, this constraint is applied separately in the two horizontal directions, 974 and the density gradient in \autoref{eq: tra_bbl_coef} is evaluated with the log gradient formulation:992 and the density gradient in \autoref{eq:TRA_bbl_coef} is evaluated with the log gradient formulation: 975 993 \[ 976 % \label{eq: tra_bbl_Drho}994 % \label{eq:TRA_bbl_Drho} 977 995 \nabla_\sigma \rho / \rho = \alpha \, \nabla_\sigma T + \beta \, \nabla_\sigma S 978 996 \] … … 983 1001 % Advective BBL 984 1002 % ------------------------------------------------------------------------------------------------------------- 985 \subsection{Advective bottom boundary layer (\protect\np{nn\_bbl\_adv}~\forcode{= 1..2})} 1003 \subsection[Advective bottom boundary layer (\forcode{nn_bbl_adv=[12]})] 1004 {Advective bottom boundary layer (\protect\np{nn\_bbl\_adv}\forcode{=[12]})} 986 1005 \label{subsec:TRA_bbl_adv} 987 1006 … … 993 1012 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 994 1013 \begin{figure}[!t] 995 \ begin{center}996 \includegraphics[]{Fig_BBL_adv}997 \caption{998 \protect\label{fig:bbl}999 Advective/diffusive Bottom Boundary Layer.1000 The BBL parameterisation is activated when $\rho^i_{kup}$ is larger than $\rho^{i + 1}_{kdnw}$.1001 Red arrows indicate the additional overturning circulation due to the advective BBL.1002 The transport of the downslope flow is defined eitheras the transport of the bottom ocean cell (black arrow),1003 1004 The green arrow indicates the diffusive BBL flux directly connecting $kup$ and $kdwn$ ocean bottom cells.1005 }1006 \ end{center}1014 \centering 1015 \includegraphics[width=\textwidth]{Fig_BBL_adv} 1016 \caption[Advective/diffusive bottom boundary layer]{ 1017 Advective/diffusive Bottom Boundary Layer. 1018 The BBL parameterisation is activated when $\rho^i_{kup}$ is larger than $\rho^{i + 1}_{kdnw}$. 1019 Red arrows indicate the additional overturning circulation due to the advective BBL. 1020 The transport of the downslope flow is defined either 1021 as the transport of the bottom ocean cell (black arrow), 1022 or as a function of the along slope density gradient. 1023 The green arrow indicates the diffusive BBL flux directly connecting 1024 $kup$ and $kdwn$ ocean bottom cells.} 1025 \label{fig:TRA_bbl} 1007 1026 \end{figure} 1008 1027 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 1014 1033 %%%gmcomment : this section has to be really written 1015 1034 1016 When applying an advective BBL (\np{nn\_bbl\_adv} ~\forcode{=1..2}), an overturning circulation is added which1035 When applying an advective BBL (\np{nn\_bbl\_adv}\forcode{=1..2}), an overturning circulation is added which 1017 1036 connects two adjacent bottom grid-points only if dense water overlies less dense water on the slope. 1018 1037 The density difference causes dense water to move down the slope. 1019 1038 1020 \np{nn\_bbl\_adv} ~\forcode{=1}:1039 \np{nn\_bbl\_adv}\forcode{=1}: 1021 1040 the downslope velocity is chosen to be the Eulerian ocean velocity just above the topographic step 1022 (see black arrow in \autoref{fig: bbl}) \citep{Beckmann_Doscher1997}.1041 (see black arrow in \autoref{fig:TRA_bbl}) \citep{beckmann.doscher_JPO97}. 1023 1042 It is a \textit{conditional advection}, that is, advection is allowed only 1024 if dense water overlies less dense water on the slope (\ie $\nabla_\sigma \rho \cdot \nabla H < 0$) and1025 if the velocity is directed towards greater depth (\ie $\vect U \cdot \nabla H > 0$).1026 1027 \np{nn\_bbl\_adv} ~\forcode{=2}:1043 if dense water overlies less dense water on the slope (\ie\ $\nabla_\sigma \rho \cdot \nabla H < 0$) and 1044 if the velocity is directed towards greater depth (\ie\ $\vect U \cdot \nabla H > 0$). 1045 1046 \np{nn\_bbl\_adv}\forcode{=2}: 1028 1047 the downslope velocity is chosen to be proportional to $\Delta \rho$, 1029 the density difference between the higher cell and lower cell densities \citep{ Campin_Goosse_Tel99}.1048 the density difference between the higher cell and lower cell densities \citep{campin.goosse_T99}. 1030 1049 The advection is allowed only if dense water overlies less dense water on the slope 1031 (\ie $\nabla_\sigma \rho \cdot \nabla H < 0$).1032 For example, the resulting transport of the downslope flow, here in the $i$-direction (\autoref{fig: bbl}),1050 (\ie\ $\nabla_\sigma \rho \cdot \nabla H < 0$). 1051 For example, the resulting transport of the downslope flow, here in the $i$-direction (\autoref{fig:TRA_bbl}), 1033 1052 is simply given by the following expression: 1034 1053 \[ 1035 % \label{eq: bbl_Utr}1054 % \label{eq:TRA_bbl_Utr} 1036 1055 u^{tr}_{bbl} = \gamma g \frac{\Delta \rho}{\rho_o} e_{1u} \, min ({e_{3u}}_{kup},{e_{3u}}_{kdwn}) 1037 1056 \] … … 1041 1060 The parameter $\gamma$ should take a different value for each bathymetric step, but for simplicity, 1042 1061 and because no direct estimation of this parameter is available, a uniform value has been assumed. 1043 The possible values for $\gamma$ range between 1 and $10~s$ \citep{ Campin_Goosse_Tel99}.1062 The possible values for $\gamma$ range between 1 and $10~s$ \citep{campin.goosse_T99}. 1044 1063 1045 1064 Scalar properties are advected by this additional transport $(u^{tr}_{bbl},v^{tr}_{bbl})$ using the upwind scheme. … … 1047 1066 the surrounding water at intermediate depths. 1048 1067 The entrainment is replaced by the vertical mixing implicit in the advection scheme. 1049 Let us consider as an example the case displayed in \autoref{fig: bbl} where1068 Let us consider as an example the case displayed in \autoref{fig:TRA_bbl} where 1050 1069 the density at level $(i,kup)$ is larger than the one at level $(i,kdwn)$. 1051 1070 The advective BBL scheme modifies the tracer time tendency of the ocean cells near the topographic step by 1052 the downslope flow \autoref{eq: bbl_dw}, the horizontal \autoref{eq:bbl_hor} and1053 the upward \autoref{eq: bbl_up} return flows as follows:1071 the downslope flow \autoref{eq:TRA_bbl_dw}, the horizontal \autoref{eq:TRA_bbl_hor} and 1072 the upward \autoref{eq:TRA_bbl_up} return flows as follows: 1054 1073 \begin{alignat}{3} 1055 \label{eq: bbl_dw}1074 \label{eq:TRA_bbl_dw} 1056 1075 \partial_t T^{do}_{kdw} &\equiv \partial_t T^{do}_{kdw} 1057 1076 &&+ \frac{u^{tr}_{bbl}}{{b_t}^{do}_{kdw}} &&\lt( T^{sh}_{kup} - T^{do}_{kdw} \rt) \\ 1058 \label{eq: bbl_hor}1077 \label{eq:TRA_bbl_hor} 1059 1078 \partial_t T^{sh}_{kup} &\equiv \partial_t T^{sh}_{kup} 1060 1079 &&+ \frac{u^{tr}_{bbl}}{{b_t}^{sh}_{kup}} &&\lt( T^{do}_{kup} - T^{sh}_{kup} \rt) \\ … … 1062 1081 \intertext{and for $k =kdw-1,\;..., \; kup$ :} 1063 1082 % 1064 \label{eq: bbl_up}1083 \label{eq:TRA_bbl_up} 1065 1084 \partial_t T^{do}_{k} &\equiv \partial_t S^{do}_{k} 1066 1085 &&+ \frac{u^{tr}_{bbl}}{{b_t}^{do}_{k}} &&\lt( T^{do}_{k +1} - T^{sh}_{k} \rt) … … 1074 1093 % Tracer damping 1075 1094 % ================================================================ 1076 \section{Tracer damping (\protect\mdl{tradmp})} 1095 \section[Tracer damping (\textit{tradmp.F90})] 1096 {Tracer damping (\protect\mdl{tradmp})} 1077 1097 \label{sec:TRA_dmp} 1078 1098 %--------------------------------------------namtra_dmp------------------------------------------------- 1079 1099 1080 \nlst{namtra_dmp} 1100 \begin{listing} 1101 \nlst{namtra_dmp} 1102 \caption{\texttt{namtra\_dmp}} 1103 \label{lst:namtra_dmp} 1104 \end{listing} 1081 1105 %-------------------------------------------------------------------------------------------------------------- 1082 1106 1083 1107 In some applications it can be useful to add a Newtonian damping term into the temperature and salinity equations: 1084 1108 \begin{equation} 1085 \label{eq: tra_dmp}1109 \label{eq:TRA_dmp} 1086 1110 \begin{gathered} 1087 1111 \pd[T]{t} = \cdots - \gamma (T - T_o) \\ 1088 1112 \pd[S]{t} = \cdots - \gamma (S - S_o) 1089 1113 \end{gathered} 1090 \end{equation} 1114 \end{equation} 1091 1115 where $\gamma$ is the inverse of a time scale, and $T_o$ and $S_o$ are given temperature and salinity fields 1092 1116 (usually a climatology). 1093 Options are defined through the \n gn{namtra\_dmp} namelist variables.1117 Options are defined through the \nam{tra\_dmp} namelist variables. 1094 1118 The restoring term is added when the namelist parameter \np{ln\_tradmp} is set to true. 1095 It also requires that both \np{ln\_tsd\_init} and \np{ln\_tsd\_ tradmp} are set to true in1096 \n gn{namtsd} namelist as well as \np{sn\_tem} and \np{sn\_sal} structures are correctly set1097 (\ie that $T_o$ and $S_o$ are provided in input files and read using \mdl{fldread},1119 It also requires that both \np{ln\_tsd\_init} and \np{ln\_tsd\_dmp} are set to true in 1120 \nam{tsd} namelist as well as \np{sn\_tem} and \np{sn\_sal} structures are correctly set 1121 (\ie\ that $T_o$ and $S_o$ are provided in input files and read using \mdl{fldread}, 1098 1122 see \autoref{subsec:SBC_fldread}). 1099 1123 The restoring coefficient $\gamma$ is a three-dimensional array read in during the \rou{tra\_dmp\_init} routine. … … 1101 1125 The DMP\_TOOLS tool is provided to allow users to generate the netcdf file. 1102 1126 1103 The two main cases in which \autoref{eq: tra_dmp} is used are1127 The two main cases in which \autoref{eq:TRA_dmp} is used are 1104 1128 \textit{(a)} the specification of the boundary conditions along artificial walls of a limited domain basin and 1105 1129 \textit{(b)} the computation of the velocity field associated with a given $T$-$S$ field … … 1109 1133 In the vicinity of these walls, $\gamma$ takes large values (equivalent to a time scale of a few days) whereas 1110 1134 it is zero in the interior of the model domain. 1111 The second case corresponds to the use of the robust diagnostic method \citep{ Sarmiento1982}.1135 The second case corresponds to the use of the robust diagnostic method \citep{sarmiento.bryan_JGR82}. 1112 1136 It allows us to find the velocity field consistent with the model dynamics whilst 1113 1137 having a $T$, $S$ field close to a given climatological field ($T_o$, $S_o$). … … 1121 1145 only below the mixed layer (defined either on a density or $S_o$ criterion). 1122 1146 It is common to set the damping to zero in the mixed layer as the adjustment time scale is short here 1123 \citep{ Madec_al_JPO96}.1147 \citep{madec.delecluse.ea_JPO96}. 1124 1148 1125 1149 For generating \ifile{resto}, see the documentation for the DMP tool provided with the source code under … … 1129 1153 % Tracer time evolution 1130 1154 % ================================================================ 1131 \section{Tracer time evolution (\protect\mdl{tranxt})} 1155 \section[Tracer time evolution (\textit{tranxt.F90})] 1156 {Tracer time evolution (\protect\mdl{tranxt})} 1132 1157 \label{sec:TRA_nxt} 1133 1158 %--------------------------------------------namdom----------------------------------------------------- 1134 1135 \nlst{namdom}1136 1159 %-------------------------------------------------------------------------------------------------------------- 1137 1160 1138 Options are defined through the \n gn{namdom} namelist variables.1139 The general framework for tracer time stepping is a modified leap-frog scheme \citep{ Leclair_Madec_OM09},1140 \ie a three level centred time scheme associated with a Asselin time filter (cf. \autoref{sec:STP_mLF}):1141 \begin{equation} 1142 \label{eq: tra_nxt}1161 Options are defined through the \nam{dom} namelist variables. 1162 The general framework for tracer time stepping is a modified leap-frog scheme \citep{leclair.madec_OM09}, 1163 \ie\ a three level centred time scheme associated with a Asselin time filter (cf. \autoref{sec:TD_mLF}): 1164 \begin{equation} 1165 \label{eq:TRA_nxt} 1143 1166 \begin{alignedat}{3} 1144 1167 &(e_{3t}T)^{t + \rdt} &&= (e_{3t}T)_f^{t - \rdt} &&+ 2 \, \rdt \,e_{3t}^t \ \text{RHS}^t \\ 1145 1168 &(e_{3t}T)_f^t &&= (e_{3t}T)^t &&+ \, \gamma \, \lt[ (e_{3t}T)_f^{t - \rdt} - 2(e_{3t}T)^t + (e_{3t}T)^{t + \rdt} \rt] \\ 1146 & && &&- \, \gamma \, \rdt \, \lt[ Q^{t + \rdt/2} - Q^{t - \rdt/2} \rt] 1169 & && &&- \, \gamma \, \rdt \, \lt[ Q^{t + \rdt/2} - Q^{t - \rdt/2} \rt] 1147 1170 \end{alignedat} 1148 \end{equation} 1171 \end{equation} 1149 1172 where RHS is the right hand side of the temperature equation, the subscript $f$ denotes filtered values, 1150 1173 $\gamma$ is the Asselin coefficient, and $S$ is the total forcing applied on $T$ 1151 (\ie fluxes plus content in mass exchanges).1174 (\ie\ fluxes plus content in mass exchanges). 1152 1175 $\gamma$ is initialized as \np{rn\_atfp} (\textbf{namelist} parameter). 1153 Its default value is \np{rn\_atfp} ~\forcode{=10.e-3}.1176 Its default value is \np{rn\_atfp}\forcode{=10.e-3}. 1154 1177 Note that the forcing correction term in the filter is not applied in linear free surface 1155 (\jp{l k\_vvl}~\forcode{= .false.}) (see \autoref{subsec:TRA_sbc}).1178 (\jp{ln\_linssh}\forcode{=.true.}) (see \autoref{subsec:TRA_sbc}). 1156 1179 Not also that in constant volume case, the time stepping is performed on $T$, not on its content, $e_{3t}T$. 1157 1180 … … 1164 1187 1165 1188 % ================================================================ 1166 % Equation of State (eosbn2) 1167 % ================================================================ 1168 \section{Equation of state (\protect\mdl{eosbn2}) } 1189 % Equation of State (eosbn2) 1190 % ================================================================ 1191 \section[Equation of state (\textit{eosbn2.F90})] 1192 {Equation of state (\protect\mdl{eosbn2})} 1169 1193 \label{sec:TRA_eosbn2} 1170 1194 %--------------------------------------------nameos----------------------------------------------------- 1171 1195 1172 \nlst{nameos} 1196 \begin{listing} 1197 \nlst{nameos} 1198 \caption{\texttt{nameos}} 1199 \label{lst:nameos} 1200 \end{listing} 1173 1201 %-------------------------------------------------------------------------------------------------------------- 1174 1202 … … 1176 1204 % Equation of State 1177 1205 % ------------------------------------------------------------------------------------------------------------- 1178 \subsection{Equation of seawater (\protect\np{nn\_eos}~\forcode{= -1..1})} 1206 \subsection[Equation of seawater (\texttt{ln}\{\texttt{\_teso10,\_eos80,\_seos}\})] 1207 {Equation of seawater (\protect\np{ln\_teos10}, \protect\np{ln\_teos80}, or \protect\np{ln\_seos}) } 1179 1208 \label{subsec:TRA_eos} 1209 1180 1210 1181 1211 The Equation Of Seawater (EOS) is an empirical nonlinear thermodynamic relationship linking seawater density, … … 1186 1216 Nonlinearities of the EOS are of major importance, in particular influencing the circulation through 1187 1217 determination of the static stability below the mixed layer, 1188 thus controlling rates of exchange between the atmosphere and the ocean interior \citep{ Roquet_JPO2015}.1189 Therefore an accurate EOS based on either the 1980 equation of state (EOS-80, \cite{ UNESCO1983}) or1190 TEOS-10 \citep{ TEOS10} standards should be used anytime a simulation of the real ocean circulation is attempted1191 \citep{ Roquet_JPO2015}.1218 thus controlling rates of exchange between the atmosphere and the ocean interior \citep{roquet.madec.ea_JPO15}. 1219 Therefore an accurate EOS based on either the 1980 equation of state (EOS-80, \cite{fofonoff.millard_bk83}) or 1220 TEOS-10 \citep{ioc.iapso_bk10} standards should be used anytime a simulation of the real ocean circulation is attempted 1221 \citep{roquet.madec.ea_JPO15}. 1192 1222 The use of TEOS-10 is highly recommended because 1193 1223 \textit{(i)} it is the new official EOS, 1194 1224 \textit{(ii)} it is more accurate, being based on an updated database of laboratory measurements, and 1195 1225 \textit{(iii)} it uses Conservative Temperature and Absolute Salinity (instead of potential temperature and 1196 practical salinity for EOS- 980, both variables being more suitable for use as model variables1197 \citep{ TEOS10, Graham_McDougall_JPO13}.1198 EOS-80 is an obsolescent feature of the NEMOsystem, kept only for backward compatibility.1226 practical salinity for EOS-80, both variables being more suitable for use as model variables 1227 \citep{ioc.iapso_bk10, graham.mcdougall_JPO13}. 1228 EOS-80 is an obsolescent feature of the \NEMO\ system, kept only for backward compatibility. 1199 1229 For process studies, it is often convenient to use an approximation of the EOS. 1200 To that purposed, a simplified EOS (S-EOS) inspired by \citet{ Vallis06} is also available.1230 To that purposed, a simplified EOS (S-EOS) inspired by \citet{vallis_bk06} is also available. 1201 1231 1202 1232 In the computer code, a density anomaly, $d_a = \rho / \rho_o - 1$, is computed, with $\rho_o$ a reference density. … … 1204 1234 This is a sensible choice for the reference density used in a Boussinesq ocean climate model, as, 1205 1235 with the exception of only a small percentage of the ocean, 1206 density in the World Ocean varies by no more than 2$\%$ from that value \citep{Gill1982}. 1207 1208 Options are defined through the \ngn{nameos} namelist variables, and in particular \np{nn\_eos} which 1209 controls the EOS used (\forcode{= -1} for TEOS10 ; \forcode{= 0} for EOS-80 ; \forcode{= 1} for S-EOS). 1236 density in the World Ocean varies by no more than 2$\%$ from that value \citep{gill_bk82}. 1237 1238 Options which control the EOS used are defined through the \nam{eos} namelist variables. 1210 1239 1211 1240 \begin{description} 1212 \item[\np{ nn\_eos}~\forcode{= -1}]1213 the polyTEOS10-bsq equation of seawater \citep{ Roquet_OM2015} is used.1241 \item[\np{ln\_teos10}\forcode{=.true.}] 1242 the polyTEOS10-bsq equation of seawater \citep{roquet.madec.ea_OM15} is used. 1214 1243 The accuracy of this approximation is comparable to the TEOS-10 rational function approximation, 1215 1244 but it is optimized for a boussinesq fluid and the polynomial expressions have simpler and … … 1217 1246 use in ocean models. 1218 1247 Note that a slightly higher precision polynomial form is now used replacement of 1219 the TEOS-10 rational function approximation for hydrographic data analysis \citep{ TEOS10}.1248 the TEOS-10 rational function approximation for hydrographic data analysis \citep{ioc.iapso_bk10}. 1220 1249 A key point is that conservative state variables are used: 1221 1250 Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: \deg{C}, notation: $\Theta$). 1222 1251 The pressure in decibars is approximated by the depth in meters. 1223 1252 With TEOS10, the specific heat capacity of sea water, $C_p$, is a constant. 1224 It is set to $C_p = 3991.86795711963~J\,Kg^{-1}\,^{\circ}K^{-1}$, according to \citet{ TEOS10}.1253 It is set to $C_p = 3991.86795711963~J\,Kg^{-1}\,^{\circ}K^{-1}$, according to \citet{ioc.iapso_bk10}. 1225 1254 Choosing polyTEOS10-bsq implies that the state variables used by the model are $\Theta$ and $S_A$. 1226 1255 In particular, the initial state deined by the user have to be given as \textit{Conservative} Temperature and 1227 1256 \textit{Absolute} Salinity. 1228 In addition, setting \np{ln\_useCT} to \forcode{.true.} convert the Conservative SSTto potential SST prior to1257 In addition, when using TEOS10, the Conservative SST is converted to potential SST prior to 1229 1258 either computing the air-sea and ice-sea fluxes (forced mode) or 1230 1259 sending the SST field to the atmosphere (coupled mode). 1231 \item[\np{ nn\_eos}~\forcode{= 0}]1260 \item[\np{ln\_eos80}\forcode{=.true.}] 1232 1261 the polyEOS80-bsq equation of seawater is used. 1233 1262 It takes the same polynomial form as the polyTEOS10, but the coefficients have been optimized to … … 1238 1267 The pressure in decibars is approximated by the depth in meters. 1239 1268 With thsi EOS, the specific heat capacity of sea water, $C_p$, is a function of temperature, salinity and 1240 pressure \citep{ UNESCO1983}.1269 pressure \citep{fofonoff.millard_bk83}. 1241 1270 Nevertheless, a severe assumption is made in order to have a heat content ($C_p T_p$) which 1242 1271 is conserved by the model: $C_p$ is set to a constant value, the TEOS10 value. 1243 \item[\np{ nn\_eos}~\forcode{= 1}]1244 a simplified EOS (S-EOS) inspired by \citet{ Vallis06} is chosen,1272 \item[\np{ln\_seos}\forcode{=.true.}] 1273 a simplified EOS (S-EOS) inspired by \citet{vallis_bk06} is chosen, 1245 1274 the coefficients of which has been optimized to fit the behavior of TEOS10 1246 (Roquet, personal comm.) (see also \citet{ Roquet_JPO2015}).1275 (Roquet, personal comm.) (see also \citet{roquet.madec.ea_JPO15}). 1247 1276 It provides a simplistic linear representation of both cabbeling and thermobaricity effects which 1248 is enough for a proper treatment of the EOS in theoretical studies \citep{ Roquet_JPO2015}.1277 is enough for a proper treatment of the EOS in theoretical studies \citep{roquet.madec.ea_JPO15}. 1249 1278 With such an equation of state there is no longer a distinction between 1250 1279 \textit{conservative} and \textit{potential} temperature, 1251 1280 as well as between \textit{absolute} and \textit{practical} salinity. 1252 1281 S-EOS takes the following expression: 1282 1253 1283 \begin{gather*} 1254 % \label{eq: tra_S-EOS}1284 % \label{eq:TRA_S-EOS} 1255 1285 \begin{alignedat}{2} 1256 1286 &d_a(T,S,z) = \frac{1}{\rho_o} \big[ &- a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * &T_a \big. \\ 1257 & &+ b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * &S_a \\ 1287 & &+ b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * &S_a \\ 1258 1288 & \big. &- \nu \; T_a &S_a \big] \\ 1259 1289 \end{alignedat} … … 1261 1291 \text{with~} T_a = T - 10 \, ; \, S_a = S - 35 \, ; \, \rho_o = 1026~Kg/m^3 1262 1292 \end{gather*} 1263 where the computer name of the coefficients as well as their standard value are given in \autoref{tab: SEOS}.1293 where the computer name of the coefficients as well as their standard value are given in \autoref{tab:TRA_SEOS}. 1264 1294 In fact, when choosing S-EOS, various approximation of EOS can be specified simply by 1265 1295 changing the associated coefficients. … … 1272 1302 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1273 1303 \begin{table}[!tb] 1274 \begin{center} 1275 \begin{tabular}{|l|l|l|l|} 1276 \hline 1277 coeff. & computer name & S-EOS & description \\ 1278 \hline 1279 $a_0$ & \np{rn\_a0} & $1.6550~10^{-1}$ & linear thermal expansion coeff. \\ 1280 \hline 1281 $b_0$ & \np{rn\_b0} & $7.6554~10^{-1}$ & linear haline expansion coeff. \\ 1282 \hline 1283 $\lambda_1$ & \np{rn\_lambda1}& $5.9520~10^{-2}$ & cabbeling coeff. in $T^2$ \\ 1284 \hline 1285 $\lambda_2$ & \np{rn\_lambda2}& $5.4914~10^{-4}$ & cabbeling coeff. in $S^2$ \\ 1286 \hline 1287 $\nu$ & \np{rn\_nu} & $2.4341~10^{-3}$ & cabbeling coeff. in $T \, S$ \\ 1288 \hline 1289 $\mu_1$ & \np{rn\_mu1} & $1.4970~10^{-4}$ & thermobaric coeff. in T \\ 1290 \hline 1291 $\mu_2$ & \np{rn\_mu2} & $1.1090~10^{-5}$ & thermobaric coeff. in S \\ 1292 \hline 1293 \end{tabular} 1294 \caption{ 1295 \protect\label{tab:SEOS} 1296 Standard value of S-EOS coefficients. 1297 } 1298 \end{center} 1304 \centering 1305 \begin{tabular}{|l|l|l|l|} 1306 \hline 1307 coeff. & computer name & S-EOS & description \\ 1308 \hline 1309 $a_0$ & \np{rn\_a0} & $1.6550~10^{-1}$ & linear thermal expansion coeff. \\ 1310 \hline 1311 $b_0$ & \np{rn\_b0} & $7.6554~10^{-1}$ & linear haline expansion coeff. \\ 1312 \hline 1313 $\lambda_1$ & \np{rn\_lambda1}& $5.9520~10^{-2}$ & cabbeling coeff. in $T^2$ \\ 1314 \hline 1315 $\lambda_2$ & \np{rn\_lambda2}& $5.4914~10^{-4}$ & cabbeling coeff. in $S^2$ \\ 1316 \hline 1317 $\nu$ & \np{rn\_nu} & $2.4341~10^{-3}$ & cabbeling coeff. in $T \, S$ \\ 1318 \hline 1319 $\mu_1$ & \np{rn\_mu1} & $1.4970~10^{-4}$ & thermobaric coeff. in T \\ 1320 \hline 1321 $\mu_2$ & \np{rn\_mu2} & $1.1090~10^{-5}$ & thermobaric coeff. in S \\ 1322 \hline 1323 \end{tabular} 1324 \caption{Standard value of S-EOS coefficients} 1325 \label{tab:TRA_SEOS} 1299 1326 \end{table} 1300 1327 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 1303 1330 % Brunt-V\"{a}is\"{a}l\"{a} Frequency 1304 1331 % ------------------------------------------------------------------------------------------------------------- 1305 \subsection{Brunt-V\"{a}is\"{a}l\"{a} frequency (\protect\np{nn\_eos}~\forcode{= 0..2})} 1332 \subsection[Brunt-V\"{a}is\"{a}l\"{a} frequency] 1333 {Brunt-V\"{a}is\"{a}l\"{a} frequency} 1306 1334 \label{subsec:TRA_bn2} 1307 1335 … … 1312 1340 In particular, $N^2$ has to be computed at the local pressure 1313 1341 (pressure in decibar being approximated by the depth in meters). 1314 The expression for $N^2$ is given by: 1342 The expression for $N^2$ is given by: 1315 1343 \[ 1316 % \label{eq: tra_bn2}1344 % \label{eq:TRA_bn2} 1317 1345 N^2 = \frac{g}{e_{3w}} \lt( \beta \; \delta_{k + 1/2}[S] - \alpha \; \delta_{k + 1/2}[T] \rt) 1318 1346 \] … … 1321 1349 The coefficients are a polynomial function of temperature, salinity and depth which expression depends on 1322 1350 the chosen EOS. 1323 They are computed through \textit{eos\_rab}, a \fortran function that can be found in \mdl{eosbn2}.1351 They are computed through \textit{eos\_rab}, a \fortran\ function that can be found in \mdl{eosbn2}. 1324 1352 1325 1353 % ------------------------------------------------------------------------------------------------------------- … … 1329 1357 \label{subsec:TRA_fzp} 1330 1358 1331 The freezing point of seawater is a function of salinity and pressure \citep{ UNESCO1983}:1332 \begin{equation} 1333 \label{eq: tra_eos_fzp}1359 The freezing point of seawater is a function of salinity and pressure \citep{fofonoff.millard_bk83}: 1360 \begin{equation} 1361 \label{eq:TRA_eos_fzp} 1334 1362 \begin{split} 1335 1363 &T_f (S,p) = \lt( a + b \, \sqrt{S} + c \, S \rt) \, S + d \, p \\ 1336 &\text{where~} a = -0.0575, \, b = 1.710523~10^{-3}, \, c = -2.154996~10^{-4} \\ 1364 &\text{where~} a = -0.0575, \, b = 1.710523~10^{-3}, \, c = -2.154996~10^{-4} \\ 1337 1365 &\text{and~} d = -7.53~10^{-3} 1338 1366 \end{split} 1339 1367 \end{equation} 1340 1368 1341 \autoref{eq: tra_eos_fzp} is only used to compute the potential freezing point of sea water1342 (\ie referenced to the surface $p = 0$),1343 thus the pressure dependent terms in \autoref{eq: tra_eos_fzp} (last term) have been dropped.1369 \autoref{eq:TRA_eos_fzp} is only used to compute the potential freezing point of sea water 1370 (\ie\ referenced to the surface $p = 0$), 1371 thus the pressure dependent terms in \autoref{eq:TRA_eos_fzp} (last term) have been dropped. 1344 1372 The freezing point is computed through \textit{eos\_fzp}, 1345 a \fortran function that can be found in \mdl{eosbn2}.1346 1347 % ------------------------------------------------------------------------------------------------------------- 1348 % Potential Energy 1373 a \fortran\ function that can be found in \mdl{eosbn2}. 1374 1375 % ------------------------------------------------------------------------------------------------------------- 1376 % Potential Energy 1349 1377 % ------------------------------------------------------------------------------------------------------------- 1350 1378 %\subsection{Potential Energy anomalies} … … 1355 1383 1356 1384 % ================================================================ 1357 % Horizontal Derivative in zps-coordinate 1358 % ================================================================ 1359 \section{Horizontal derivative in \textit{zps}-coordinate (\protect\mdl{zpshde})} 1385 % Horizontal Derivative in zps-coordinate 1386 % ================================================================ 1387 \section[Horizontal derivative in \textit{zps}-coordinate (\textit{zpshde.F90})] 1388 {Horizontal derivative in \textit{zps}-coordinate (\protect\mdl{zpshde})} 1360 1389 \label{sec:TRA_zpshde} 1361 1390 1362 \gmcomment{STEVEN: to be consistent with earlier discussion of differencing and averaging operators, 1391 \gmcomment{STEVEN: to be consistent with earlier discussion of differencing and averaging operators, 1363 1392 I've changed "derivative" to "difference" and "mean" to "average"} 1364 1393 1365 With partial cells (\np{ln\_zps} ~\forcode{= .true.}) at bottom and top (\np{ln\_isfcav}~\forcode{=.true.}),1394 With partial cells (\np{ln\_zps}\forcode{=.true.}) at bottom and top (\np{ln\_isfcav}\forcode{=.true.}), 1366 1395 in general, tracers in horizontally adjacent cells live at different depths. 1367 1396 Horizontal gradients of tracers are needed for horizontal diffusion (\mdl{traldf} module) and 1368 1397 the hydrostatic pressure gradient calculations (\mdl{dynhpg} module). 1369 The partial cell properties at the top (\np{ln\_isfcav} ~\forcode{=.true.}) are computed in the same way as1398 The partial cell properties at the top (\np{ln\_isfcav}\forcode{=.true.}) are computed in the same way as 1370 1399 for the bottom. 1371 1400 So, only the bottom interpolation is explained below. … … 1373 1402 Before taking horizontal gradients between the tracers next to the bottom, 1374 1403 a linear interpolation in the vertical is used to approximate the deeper tracer as if 1375 it actually lived at the depth of the shallower tracer point (\autoref{fig: Partial_step_scheme}).1404 it actually lived at the depth of the shallower tracer point (\autoref{fig:TRA_Partial_step_scheme}). 1376 1405 For example, for temperature in the $i$-direction the needed interpolated temperature, $\widetilde T$, is: 1377 1406 1378 1407 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1379 1408 \begin{figure}[!p] 1380 \begin{center} 1381 \includegraphics[]{Fig_partial_step_scheme} 1382 \caption{ 1383 \protect\label{fig:Partial_step_scheme} 1384 Discretisation of the horizontal difference and average of tracers in the $z$-partial step coordinate 1385 (\protect\np{ln\_zps}~\forcode{= .true.}) in the case $(e3w_k^{i + 1} - e3w_k^i) > 0$. 1386 A linear interpolation is used to estimate $\widetilde T_k^{i + 1}$, 1387 the tracer value at the depth of the shallower tracer point of the two adjacent bottom $T$-points. 1388 The horizontal difference is then given by: $\delta_{i + 1/2} T_k = \widetilde T_k^{\, i + 1} -T_k^{\, i}$ and 1389 the average by: $\overline T_k^{\, i + 1/2} = (\widetilde T_k^{\, i + 1/2} - T_k^{\, i}) / 2$. 1390 } 1391 \end{center} 1409 \centering 1410 \includegraphics[width=\textwidth]{Fig_partial_step_scheme} 1411 \caption[Discretisation of the horizontal difference and average of tracers in 1412 the $z$-partial step coordinate]{ 1413 Discretisation of the horizontal difference and average of tracers in 1414 the $z$-partial step coordinate (\protect\np{ln\_zps}\forcode{=.true.}) in 1415 the case $(e3w_k^{i + 1} - e3w_k^i) > 0$. 1416 A linear interpolation is used to estimate $\widetilde T_k^{i + 1}$, 1417 the tracer value at the depth of the shallower tracer point of 1418 the two adjacent bottom $T$-points. 1419 The horizontal difference is then given by: 1420 $\delta_{i + 1/2} T_k = \widetilde T_k^{\, i + 1} -T_k^{\, i}$ and 1421 the average by: 1422 $\overline T_k^{\, i + 1/2} = (\widetilde T_k^{\, i + 1/2} - T_k^{\, i}) / 2$.} 1423 \label{fig:TRA_Partial_step_scheme} 1392 1424 \end{figure} 1393 1425 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 1402 1434 \rt. 1403 1435 \] 1404 and the resulting forms for the horizontal difference and the horizontal average value of $T$ at a $U$-point are: 1405 \begin{equation} 1406 \label{eq: zps_hde}1436 and the resulting forms for the horizontal difference and the horizontal average value of $T$ at a $U$-point are: 1437 \begin{equation} 1438 \label{eq:TRA_zps_hde} 1407 1439 \begin{split} 1408 1440 \delta_{i + 1/2} T &= … … 1428 1460 Instead of forming a linear approximation of density, we compute $\widetilde \rho$ from the interpolated values of 1429 1461 $T$ and $S$, and the pressure at a $u$-point 1430 (in the equation of state pressure is approximated by depth, see \autoref{subsec:TRA_eos}): 1462 (in the equation of state pressure is approximated by depth, see \autoref{subsec:TRA_eos}): 1431 1463 \[ 1432 % \label{eq: zps_hde_rho}1464 % \label{eq:TRA_zps_hde_rho} 1433 1465 \widetilde \rho = \rho (\widetilde T,\widetilde S,z_u) \quad \text{where~} z_u = \min \lt( z_T^{i + 1},z_T^i \rt) 1434 1466 \] … … 1441 1473 Note that in almost all the advection schemes presented in this Chapter, 1442 1474 both averaging and differencing operators appear. 1443 Yet \autoref{eq: zps_hde} has not been used in these schemes:1475 Yet \autoref{eq:TRA_zps_hde} has not been used in these schemes: 1444 1476 in contrast to diffusion and pressure gradient computations, 1445 1477 no correction for partial steps is applied for advection.
Note: See TracChangeset
for help on using the changeset viewer.