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