Changeset 10354 for NEMO/trunk/doc/latex/NEMO/subfiles/chap_TRA.tex
- Timestamp:
- 2018-11-21T17:59:55+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/doc/latex/NEMO/subfiles/chap_TRA.tex
r10146 r10354 17 17 %$\ $\newline % force a new ligne 18 18 19 Using the representation described in \autoref{chap:DOM}, several semi-discrete 20 space forms of the tracer equations are available depending on the vertical 21 coordinate used and on the physics used. In all the equations presented 22 here, the masking has been omitted for simplicity. One must be aware that 23 all the quantities are masked fields and that each time a mean or difference 24 operator is used, the resulting field is multiplied by a mask. 25 26 The two active tracers are potential temperature and salinity. Their prognostic 27 equations can be summarized as follows: 19 Using the representation described in \autoref{chap:DOM}, 20 several semi-discrete space forms of the tracer equations are available depending on 21 the vertical coordinate used and on the physics used. 22 In all the equations presented here, the masking has been omitted for simplicity. 23 One must be aware that all the quantities are masked fields and 24 that each time a mean or difference operator is used, 25 the resulting field is multiplied by a mask. 26 27 The two active tracers are potential temperature and salinity. 28 Their prognostic equations can be summarized as follows: 28 29 \begin{equation*} 29 30 \text{NXT} = \text{ADV}+\text{LDF}+\text{ZDF}+\text{SBC} … … 31 32 \end{equation*} 32 33 33 NXT stands for next, referring to the time-stepping. From left to right, the terms 34 on the rhs of the tracer equations are the advection (ADV), the lateral diffusion 35 (LDF), the vertical diffusion (ZDF), the contributions from the external forcings 36 (SBC: Surface Boundary Condition, QSR: penetrative Solar Radiation, and BBC: 37 Bottom Boundary Condition), the contribution from the bottom boundary Layer 38 (BBL) parametrisation, and an internal damping (DMP) term. The terms QSR, 39 BBC, BBL and DMP are optional. The external forcings and parameterisations 40 require complex inputs and complex calculations ($e.g.$ bulk formulae, estimation 41 of mixing coefficients) that are carried out in the SBC, LDF and ZDF modules and 42 described in \autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively. 43 Note that \mdl{tranpc}, the non-penetrative convection module, although 44 located in the NEMO/OPA/TRA directory as it directly modifies the tracer fields, 45 is described with the model vertical physics (ZDF) together with other available 46 parameterization of convection. 47 48 In the present chapter we also describe the diagnostic equations used to compute 49 the sea-water properties (density, Brunt-V\"{a}is\"{a}l\"{a} frequency, specific heat and 50 freezing point with associated modules \mdl{eosbn2} and \mdl{phycst}). 51 52 The different options available to the user are managed by namelist logicals or CPP keys. 53 For each equation term \textit{TTT}, the namelist logicals are \textit{ln\_traTTT\_xxx}, 54 where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme. 55 The CPP key (when it exists) is \key{traTTT}. The equivalent code can be 56 found in the \textit{traTTT} or \textit{traTTT\_xxx} module, in the NEMO/OPA/TRA directory. 57 58 The user has the option of extracting each tendency term on the RHS of the tracer 59 equation for output (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}\forcode{ = .true.}), as described in \autoref{chap:DIA}. 34 NXT stands for next, referring to the time-stepping. 35 From left to right, the terms on the rhs of the tracer equations are the advection (ADV), 36 the lateral diffusion (LDF), the vertical diffusion (ZDF), the contributions from the external forcings 37 (SBC: Surface Boundary Condition, QSR: penetrative Solar Radiation, and BBC: Bottom Boundary Condition), 38 the contribution from the bottom boundary Layer (BBL) parametrisation, and an internal damping (DMP) term. 39 The terms QSR, BBC, BBL and DMP are optional. 40 The external forcings and parameterisations require complex inputs and complex calculations 41 ($e.g.$ bulk formulae, estimation of mixing coefficients) that are carried out in the SBC, LDF and ZDF modules and 42 described in \autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively. 43 Note that \mdl{tranpc}, the non-penetrative convection module, although located in the NEMO/OPA/TRA directory as 44 it directly modifies the tracer fields, is described with the model vertical physics (ZDF) together with 45 other available parameterization of convection. 46 47 In the present chapter we also describe the diagnostic equations used to compute the sea-water properties 48 (density, Brunt-V\"{a}is\"{a}l\"{a} frequency, specific heat and freezing point with 49 associated modules \mdl{eosbn2} and \mdl{phycst}). 50 51 The different options available to the user are managed by namelist logicals or CPP keys. 52 For each equation term \textit{TTT}, the namelist logicals are \textit{ln\_traTTT\_xxx}, 53 where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme. 54 The CPP key (when it exists) is \key{traTTT}. 55 The equivalent code can be found in the \textit{traTTT} or \textit{traTTT\_xxx} module, 56 in the NEMO/OPA/TRA directory. 57 58 The user has the option of extracting each tendency term on the RHS of the tracer equation for output 59 (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}\forcode{ = .true.}), as described in \autoref{chap:DIA}. 60 60 61 61 $\ $\newline % force a new ligne … … 70 70 %------------------------------------------------------------------------------------------------------------- 71 71 72 When considered ($i.e.$ when \np{ln\_traadv\_NONE} is not set to \forcode{.true.}), 73 the advection tendency of a tracer is expressed in flux form, 74 $i.e.$ as the divergence of the advective fluxes. Its discrete expression is given by : 72 When considered ($i.e.$ when \np{ln\_traadv\_NONE} is not set to \forcode{.true.}), 73 the advection tendency of a tracer is expressed in flux form, 74 $i.e.$ as the divergence of the advective fluxes. 75 Its discrete expression is given by : 75 76 \begin{equation} \label{eq:tra_adv} 76 77 ADV_\tau =-\frac{1}{b_t} \left( … … 79 80 -\frac{1}{e_{3t}} \;\delta _k \left[ w\; \tau _w \right] 80 81 \end{equation} 81 where $\tau$ is either T or S, and $b_t= e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$-cells. 82 The flux form in \autoref{eq:tra_adv} 83 implicitly requires the use of the continuity equation. Indeed, it is obtained 84 by using the following equality : $\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$ 85 which results from the use of the continuity equation, $\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ 86 (which reduces to $\nabla \cdot \vect{U}=0$ in linear free surface, $i.e.$ \np{ln\_linssh}\forcode{ = .true.}). 87 Therefore it is of paramount importance to design the discrete analogue of the 88 advection tendency so that it is consistent with the continuity equation in order to 89 enforce the conservation properties of the continuous equations. In other words, 90 by setting $\tau = 1$ in (\autoref{eq:tra_adv}) we recover the discrete form of 82 where $\tau$ is either T or S, and $b_t= e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$-cells. 83 The flux form in \autoref{eq:tra_adv} implicitly requires the use of the continuity equation. 84 Indeed, it is obtained by using the following equality: 85 $\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$ which 86 results from the use of the continuity equation, $\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ 87 (which reduces to $\nabla \cdot \vect{U}=0$ in linear free surface, $i.e.$ \np{ln\_linssh}\forcode{ = .true.}). 88 Therefore it is of paramount importance to design the discrete analogue of the advection tendency so that 89 it is consistent with the continuity equation in order to enforce the conservation properties of 90 the continuous equations. 91 In other words, by setting $\tau = 1$ in (\autoref{eq:tra_adv}) we recover the discrete form of 91 92 the continuity equation which is used to calculate the vertical velocity. 92 93 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 93 \begin{figure}[!t] \begin{center} 94 \includegraphics[width=0.9\textwidth]{Fig_adv_scheme} 95 \caption{ \protect\label{fig:adv_scheme} 96 Schematic representation of some ways used to evaluate the tracer value 97 at $u$-point and the amount of tracer exchanged between two neighbouring grid 98 points. Upsteam biased scheme (ups): the upstream value is used and the black 99 area is exchanged. Piecewise parabolic method (ppm): a parabolic interpolation 100 is used and the black and dark grey areas are exchanged. Monotonic upstream 101 scheme for conservative laws (muscl): a parabolic interpolation is used and black, 102 dark grey and grey areas are exchanged. Second order scheme (cen2): the mean 103 value is used and black, dark grey, grey and light grey areas are exchanged. Note 104 that this illustration does not include the flux limiter used in ppm and muscl schemes.} 105 \end{center} \end{figure} 94 \begin{figure}[!t] 95 \begin{center} 96 \includegraphics[width=0.9\textwidth]{Fig_adv_scheme} 97 \caption{ \protect\label{fig:adv_scheme} 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 } 110 \end{center} 111 \end{figure} 106 112 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 107 113 108 The key difference between the advection schemes available in \NEMO is the choice 109 made in space and time interpolation to define the value of the tracer at the 110 velocity points (\autoref{fig:adv_scheme}). 111 112 Along solid lateral and bottom boundaries a zero tracer flux is automatically 113 s pecified, since the normal velocity is zero there. At the sea surface the114 boundary condition depends on the type of sea surface chosen: 114 The key difference between the advection schemes available in \NEMO is the choice made in space and 115 time interpolation to define the value of the tracer at the velocity points 116 (\autoref{fig:adv_scheme}). 117 118 Along solid lateral and bottom boundaries a zero tracer flux is automatically specified, 119 since the normal velocity is zero there. 120 At the sea surface the boundary condition depends on the type of sea surface chosen: 115 121 \begin{description} 116 \item [linear free surface:] (\np{ln\_linssh}\forcode{ = .true.}) the first level thickness is constant in time: 117 the vertical boundary condition is applied at the fixed surface $z=0$ 118 rather than on the moving surface $z=\eta$. There is a non-zero advective 119 flux which is set for all advection schemes as 120 $\left. {\tau _w } \right|_{k=1/2} =T_{k=1} $, $i.e.$ 121 the product of surface velocity (at $z=0$) by the first level tracer value. 122 \item [non-linear free surface:] (\np{ln\_linssh}\forcode{ = .false.}) 123 convergence/divergence in the first ocean level moves the free surface 124 up/down. There is no tracer advection through it so that the advective 125 fluxes through the surface are also zero 122 \item[linear free surface:] 123 (\np{ln\_linssh}\forcode{ = .true.}) 124 the first level thickness is constant in time: 125 the vertical boundary condition is applied at the fixed surface $z=0$ rather than on the moving surface $z=\eta$. 126 There is a non-zero advective flux which is set for all advection schemes as 127 $\left. {\tau _w } \right|_{k=1/2} =T_{k=1} $, 128 $i.e.$ the product of surface velocity (at $z=0$) by the first level tracer value. 129 \item[non-linear free surface:] 130 (\np{ln\_linssh}\forcode{ = .false.}) 131 convergence/divergence in the first ocean level moves the free surface up/down. 132 There is no tracer advection through it so that the advective fluxes through the surface are also zero. 126 133 \end{description} 127 In all cases, this boundary condition retains local conservation of tracer. 128 Global conservation is obtained in non-linear free surface case, 129 but \textit{not} in the linear free surface case. Nevertheless, in the latter case, 130 it is achieved to a good approximation since the non-conservative 131 term is the product of the time derivative of the tracer and the free surface 132 height, two quantities that are not correlated \citep{Roullet_Madec_JGR00,Griffies_al_MWR01,Campin2004}. 133 134 The velocity field that appears in (\autoref{eq:tra_adv}) and (\autoref{eq:tra_adv_zco}) 135 is the centred (\textit{now}) \textit{effective} ocean velocity, $i.e.$ the \textit{eulerian} velocity 136 (see \autoref{chap:DYN}) plus the eddy induced velocity (\textit{eiv}) 137 and/or the mixed layer eddy induced velocity (\textit{eiv}) 138 when those parameterisations are used (see \autoref{chap:LDF}). 139 140 Several tracer advection scheme are proposed, namely 141 a $2^{nd}$ or $4^{th}$ order centred schemes (CEN), 134 In all cases, this boundary condition retains local conservation of tracer. 135 Global conservation is obtained in non-linear free surface case, but \textit{not} in the linear free surface case. 136 Nevertheless, in the latter case, it is achieved to a good approximation since 137 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}) 141 is the centred (\textit{now}) \textit{effective} ocean velocity, 142 $i.e.$ the \textit{eulerian} velocity (see \autoref{chap:DYN}) plus 143 the eddy induced velocity (\textit{eiv}) and/or 144 the mixed layer eddy induced velocity (\textit{eiv}) when 145 those parameterisations are used (see \autoref{chap:LDF}). 146 147 Several tracer advection scheme are proposed, namely a $2^{nd}$ or $4^{th}$ order centred schemes (CEN), 142 148 a $2^{nd}$ or $4^{th}$ order Flux Corrected Transport scheme (FCT), 143 149 a Monotone Upstream Scheme for Conservative Laws scheme (MUSCL), 144 a $3^{rd}$ Upstream Biased Scheme (UBS, also often called UP3), and145 a Quadratic Upstream Interpolation for Convective Kinematics with150 a $3^{rd}$ Upstream Biased Scheme (UBS, also often called UP3), 151 and a Quadratic Upstream Interpolation for Convective Kinematics with 146 152 Estimated Streaming Terms scheme (QUICKEST). 147 The choice is made in the \textit{\ngn{namtra\_adv}} namelist, by148 setting to \forcode{.true.} one of the logicals \textit{ln\_traadv\_xxx}. 149 The corresponding code can be found in the \mdl{traadv\_xxx} module, 150 where \textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme. 151 By default ($i.e.$ in the reference namelist, \ngn{namelist\_ref}), all the logicals 152 are set to \forcode{.false.}. If the user does not select an advection scheme 153 in the configuration namelist (\ngn{namelist\_cfg}), the tracers will \textit{not} be advected!154 155 Details of the advection schemes are given below. The choosing an advection scheme156 is a complex matter which depends on the model physics, model resolution, 153 The choice is made in the \textit{\ngn{namtra\_adv}} namelist, 154 by setting to \forcode{.true.} one of the logicals \textit{ln\_traadv\_xxx}. 155 The corresponding code can be found in the \mdl{traadv\_xxx} module, 156 where \textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme. 157 By default ($i.e.$ in the reference namelist, \ngn{namelist\_ref}), all the logicals are set to \forcode{.false.}. 158 If the user does not select an advection scheme in the configuration namelist (\ngn{namelist\_cfg}), 159 the tracers will \textit{not} be advected! 160 161 Details of the advection schemes are given below. 162 The choosing an advection scheme is a complex matter which depends on the model physics, model resolution, 157 163 type of tracer, as well as the issue of numerical cost. In particular, we note that 158 (1) CEN and FCT schemes require an explicit diffusion operator 159 while the other schemes are diffusive enough so that they do not necessarily need additional diffusion ; 164 (1) CEN and FCT schemes require an explicit diffusion operator while the other schemes are diffusive enough so that 165 they do not necessarily need additional diffusion; 160 166 (2) CEN and UBS are not \textit{positive} schemes 161 \footnote{negative values can appear in an initially strictly positive tracer field 162 which is advected} 163 , implying that false extrema are permitted. Their use is not recommended on passive tracers ; 164 (3) It is recommended that the same advection-diffusion scheme is 165 used on both active and passive tracers. Indeed, if a source or sink of a 166 passive tracer depends on an active one, the difference of treatment of 167 a ctive and passive tracers can create very nice-looking frontal structures168 that are pure numerical artefacts. Nevertheless, most of our users set a different 169 t reatment on passive and active tracers, that's the reason why this possibility170 is offered. We strongly suggest them to perform a sensitivity experiment 171 using a same treatment toassess the robustness of their results.167 \footnote{negative values can appear in an initially strictly positive tracer field which is advected}, 168 implying that false extrema are permitted. 169 Their use is not recommended on passive tracers; 170 (3) It is recommended that the same advection-diffusion scheme is used on both active and passive tracers. 171 Indeed, if a source or sink of a passive tracer depends on an active one, 172 the difference of treatment of active and passive tracers can create very nice-looking frontal structures that 173 are pure numerical artefacts. 174 Nevertheless, most of our users set a different treatment on passive and active tracers, 175 that's the reason why this possibility is offered. 176 We strongly suggest them to perform a sensitivity experiment using a same treatment to 177 assess the robustness of their results. 172 178 173 179 % ------------------------------------------------------------------------------------------------------------- … … 179 185 % 2nd order centred scheme 180 186 181 The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen}\forcode{ = .true.}. 182 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) 183 and vertical direction by setting \np{nn\_cen\_h} and \np{nn\_cen\_v} to $2$ or $4$. 187 The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen}\forcode{ = .true.}. 188 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) and vertical direction by 189 setting \np{nn\_cen\_h} and \np{nn\_cen\_v} to $2$ or $4$. 184 190 CEN implementation can be found in the \mdl{traadv\_cen} module. 185 191 186 In the $2^{nd}$ order centred formulation (CEN2), the tracer at velocity points 187 is evaluated as the mean of the two neighbouring $T$-point values. 192 In the $2^{nd}$ order centred formulation (CEN2), the tracer at velocity points is evaluated as the mean of 193 the two neighbouring $T$-point values. 188 194 For example, in the $i$-direction : 189 195 \begin{equation} \label{eq:tra_adv_cen2} … … 191 197 \end{equation} 192 198 193 CEN2 is non diffusive ($i.e.$ it conserves the tracer variance, $\tau^2)$ 194 but dispersive ($i.e.$ it may create false extrema). It is therefore notoriously 195 noisy and must be used in conjunction with an explicit diffusion operator to 196 produce a sensible solution. The associated time-stepping is performed using 197 a leapfrog scheme in conjunction with an Asselin time-filter, so $T$ in 198 (\autoref{eq:tra_adv_cen2}) is the \textit{now} tracer value. 199 200 Note that using the CEN2, the overall tracer advection is of second 201 order accuracy since both (\autoref{eq:tra_adv}) and (\autoref{eq:tra_adv_cen2}) 202 have this order of accuracy. 199 CEN2 is non diffusive ($i.e.$ it conserves the tracer variance, $\tau^2)$ but dispersive 200 ($i.e.$ it may create false extrema). 201 It is therefore notoriously noisy and must be used in conjunction with an explicit diffusion operator to 202 produce a sensible solution. 203 The associated time-stepping is performed using a leapfrog scheme in conjunction with an Asselin time-filter, 204 so $T$ in (\autoref{eq:tra_adv_cen2}) is the \textit{now} tracer value. 205 206 Note that using the CEN2, the overall tracer advection is of second order accuracy since 207 both (\autoref{eq:tra_adv}) and (\autoref{eq:tra_adv_cen2}) have this order of accuracy. 203 208 204 209 % 4nd order centred scheme 205 210 206 In the $4^{th}$ order formulation (CEN4), tracer values are evaluated at u- and v-points as 207 a $4^{th}$ order interpolation, and thus depend on the four neighbouring $T$-points. 211 In the $4^{th}$ order formulation (CEN4), tracer values are evaluated at u- and v-points as 212 a $4^{th}$ order interpolation, and thus depend on the four neighbouring $T$-points. 208 213 For example, in the $i$-direction: 209 214 \begin{equation} \label{eq:tra_adv_cen4} … … 211 216 =\overline{ T - \frac{1}{6}\,\delta _i \left[ \delta_{i+1/2}[T] \,\right] }^{\,i+1/2} 212 217 \end{equation} 213 In the vertical direction (\np{nn\_cen\_v}\forcode{ = 4}), a $4^{th}$ COMPACT interpolation 214 has been prefered \citep{Demange_PhD2014}. 215 In the COMPACT scheme, both the field and its derivative are interpolated, 216 which leads, after a matrix inversion, spectral characteristics 217 similar to schemes of higher order \citep{Lele_JCP1992}. 218 In the vertical direction (\np{nn\_cen\_v}\forcode{ = 4}), 219 a $4^{th}$ COMPACT interpolation has been prefered \citep{Demange_PhD2014}. 220 In the COMPACT scheme, both the field and its derivative are interpolated, which leads, after a matrix inversion, 221 spectral characteristics similar to schemes of higher order \citep{Lele_JCP1992}. 218 222 219 223 220 Strictly speaking, the CEN4 scheme is not a $4^{th}$ order advection scheme 221 but a $4^{th}$ order evaluation of advective fluxes, since the divergence of 222 advective fluxes \autoref{eq:tra_adv} is kept at $2^{nd}$ order. 223 The expression \textit{$4^{th}$ order scheme} used in oceanographic literature 224 is usually associated with the scheme presented here. 225 Introducing a \forcode{.true.} $4^{th}$ order advection scheme is feasible but, 226 for consistency reasons, it requires changes in the discretisation of the tracer 227 advection together with changes in the continuity equation, 228 and the momentum advection and pressure terms. 229 230 A direct consequence of the pseudo-fourth order nature of the scheme is that 231 it is not non-diffusive, $i.e.$ the global variance of a tracer is not preserved using CEN4. 232 Furthermore, it must be used in conjunction with an explicit diffusion operator 233 to produce a sensible solution. 234 As in CEN2 case, the time-stepping is performed using a leapfrog scheme in conjunction 235 with an Asselin time-filter, so $T$ in (\autoref{eq:tra_adv_cen4}) is the \textit{now} tracer. 236 237 At a $T$-grid cell adjacent to a boundary (coastline, bottom and surface), 238 an additional hypothesis must be made to evaluate $\tau _u^{cen4}$. 239 This hypothesis usually reduces the order of the scheme. 240 Here we choose to set the gradient of $T$ across the boundary to zero. 241 Alternative conditions can be specified, such as a reduction to a second order scheme 242 for these near boundary grid points. 224 Strictly speaking, the CEN4 scheme is not a $4^{th}$ order advection scheme but 225 a $4^{th}$ order evaluation of advective fluxes, 226 since the divergence of advective fluxes \autoref{eq:tra_adv} is kept at $2^{nd}$ order. 227 The expression \textit{$4^{th}$ order scheme} used in oceanographic literature is usually associated with 228 the scheme presented here. 229 Introducing a \forcode{.true.} $4^{th}$ order advection scheme is feasible but, for consistency reasons, 230 it requires changes in the discretisation of the tracer advection together with changes in the continuity equation, 231 and the momentum advection and pressure terms. 232 233 A direct consequence of the pseudo-fourth order nature of the scheme is that it is not non-diffusive, 234 $i.e.$ the global variance of a tracer is not preserved using CEN4. 235 Furthermore, it must be used in conjunction with an explicit diffusion operator to produce a sensible solution. 236 As in CEN2 case, the time-stepping is performed using a leapfrog scheme in conjunction with an Asselin time-filter, 237 so $T$ in (\autoref{eq:tra_adv_cen4}) is the \textit{now} tracer. 238 239 At a $T$-grid cell adjacent to a boundary (coastline, bottom and surface), 240 an additional hypothesis must be made to evaluate $\tau _u^{cen4}$. 241 This hypothesis usually reduces the order of the scheme. 242 Here we choose to set the gradient of $T$ across the boundary to zero. 243 Alternative conditions can be specified, such as a reduction to a second order scheme for 244 these near boundary grid points. 243 245 244 246 % ------------------------------------------------------------------------------------------------------------- … … 248 250 \label{subsec:TRA_adv_tvd} 249 251 250 The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct}\forcode{ = .true.}. 251 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) 252 and vertical direction bysetting \np{nn\_fct\_h} and \np{nn\_fct\_v} to $2$ or $4$.252 The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct}\forcode{ = .true.}. 253 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) and vertical direction by 254 setting \np{nn\_fct\_h} and \np{nn\_fct\_v} to $2$ or $4$. 253 255 FCT implementation can be found in the \mdl{traadv\_fct} module. 254 256 255 In FCT formulation, the tracer at velocity points is evaluated using a combination of 256 an upstream and a centred scheme. For example, in the $i$-direction : 257 In FCT formulation, the tracer at velocity points is evaluated using a combination of an upstream and 258 a centred scheme. 259 For example, in the $i$-direction : 257 260 \begin{equation} \label{eq:tra_adv_fct} 258 261 \begin{split} … … 265 268 \end{split} 266 269 \end{equation} 267 where $c_u$ is a flux limiter function taking values between 0 and 1. 268 The FCT order is the one of the centred scheme used ($i.e.$ it depends on the setting of 269 \np{nn\_fct\_h} and \np{nn\_fct\_v}. 270 There exist many ways to define $c_u$, each corresponding to a different 271 FCT scheme. The one chosen in \NEMO is described in \citet{Zalesak_JCP79}. 272 $c_u$ only departs from $1$ when the advective term produces a local extremum in the tracer field. 273 The resulting scheme is quite expensive but \emph{positive}. 274 It can be used on both active and passive tracers. 275 A comparison of FCT-2 with MUSCL and a MPDATA scheme can be found in \citet{Levy_al_GRL01}. 276 277 An additional option has been added controlled by \np{nn\_fct\_zts}. By setting this integer to 278 a value larger than zero, a $2^{nd}$ order FCT scheme is used on both horizontal and vertical direction, 279 but on the latter, a split-explicit time stepping is used, with a number of sub-timestep equals 280 to \np{nn\_fct\_zts}. This option can be useful when the size of the timestep is limited 281 by vertical advection \citep{Lemarie_OM2015}. Note that in this case, a similar split-explicit 282 time stepping should be used on vertical advection of momentum to insure a better stability 283 (see \autoref{subsec:DYN_zad}). 284 285 For stability reasons (see \autoref{chap:STP}), $\tau _u^{cen}$ is evaluated in (\autoref{eq:tra_adv_fct}) 286 using the \textit{now} tracer while $\tau _u^{ups}$ is evaluated using the \textit{before} tracer. In other words, 287 the advective part of the scheme is time stepped with a leap-frog scheme 270 where $c_u$ is a flux limiter function taking values between 0 and 1. 271 The FCT order is the one of the centred scheme used 272 ($i.e.$ it depends on the setting of \np{nn\_fct\_h} and \np{nn\_fct\_v}). 273 There exist many ways to define $c_u$, each corresponding to a different FCT scheme. 274 The one chosen in \NEMO is described in \citet{Zalesak_JCP79}. 275 $c_u$ only departs from $1$ when the advective term produces a local extremum in the tracer field. 276 The resulting scheme is quite expensive but \emph{positive}. 277 It can be used on both active and passive tracers. 278 A comparison of FCT-2 with MUSCL and a MPDATA scheme can be found in \citet{Levy_al_GRL01}. 279 280 An additional option has been added controlled by \np{nn\_fct\_zts}. 281 By setting this integer to a value larger than zero, 282 a $2^{nd}$ order FCT scheme is used on both horizontal and vertical direction, but on the latter, 283 a split-explicit time stepping is used, with a number of sub-timestep equals to \np{nn\_fct\_zts}. 284 This option can be useful when the size of the timestep is limited by vertical advection \citep{Lemarie_OM2015}. 285 Note that in this case, a similar split-explicit time stepping should be used on vertical advection of momentum to 286 insure a better stability (see \autoref{subsec:DYN_zad}). 287 288 For stability reasons (see \autoref{chap:STP}), 289 $\tau _u^{cen}$ is evaluated in (\autoref{eq:tra_adv_fct}) using the \textit{now} tracer while 290 $\tau _u^{ups}$ is evaluated using the \textit{before} tracer. 291 In other words, the advective part of the scheme is time stepped with a leap-frog scheme 288 292 while a forward scheme is used for the diffusive part. 289 293 … … 294 298 \label{subsec:TRA_adv_mus} 295 299 296 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.}. 297 301 MUSCL implementation can be found in the \mdl{traadv\_mus} module. 298 302 299 MUSCL has been first implemented in \NEMO by \citet{Levy_al_GRL01}. In its formulation, the tracer at velocity points 300 is evaluated assuming a linear tracer variation between two $T$-points 301 (\autoref{fig:adv_scheme}). For example, in the $i$-direction : 303 MUSCL has been first implemented in \NEMO by \citet{Levy_al_GRL01}. 304 In its formulation, the tracer at velocity points is evaluated assuming a linear tracer variation between 305 two $T$-points (\autoref{fig:adv_scheme}). 306 For example, in the $i$-direction : 302 307 \begin{equation} \label{eq:tra_adv_mus} 303 308 \tau _u^{mus} = \left\{ \begin{aligned} … … 308 313 \end{aligned} \right. 309 314 \end{equation} 310 where $\widetilde{\partial _i \tau}$ is the slope of the tracer on which a limitation 311 is imposed toensure the \textit{positive} character of the scheme.312 313 The time stepping is performed using a forward scheme, that is the \textit{before}314 t racer field is used to evaluate $\tau _u^{mus}$.315 316 For an ocean grid point adjacent to land and where the ocean velocity is 317 directed toward land, an upstream flux is used. This choice ensure 318 the \textit{positive} character of the scheme. 319 In addition, fluxes round a grid-point where a runoff is applied can optionally be 320 computed using upstream fluxes(\np{ln\_mus\_ups}\forcode{ = .true.}).315 where $\widetilde{\partial _i \tau}$ is the slope of the tracer on which a limitation is imposed to 316 ensure the \textit{positive} character of the scheme. 317 318 The time stepping is performed using a forward scheme, 319 that is the \textit{before} tracer field is used to evaluate $\tau _u^{mus}$. 320 321 For an ocean grid point adjacent to land and where the ocean velocity is directed toward land, 322 an upstream flux is used. 323 This choice ensure the \textit{positive} character of the scheme. 324 In addition, fluxes round a grid-point where a runoff is applied can optionally be computed using upstream fluxes 325 (\np{ln\_mus\_ups}\forcode{ = .true.}). 321 326 322 327 % ------------------------------------------------------------------------------------------------------------- … … 326 331 \label{subsec:TRA_adv_ubs} 327 332 328 The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs}\forcode{ = .true.}. 333 The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs}\forcode{ = .true.}. 329 334 UBS implementation can be found in the \mdl{traadv\_mus} module. 330 335 331 The UBS scheme, often called UP3, is also known as the Cell Averaged QUICK scheme 332 (Quadratic Upstream Interpolation for Convective Kinematics). It is an upstream-biased333 third order scheme based on an upstream-biased parabolic interpolation. 334 For example, in the $i$-direction 336 The UBS scheme, often called UP3, is also known as the Cell Averaged QUICK scheme 337 (Quadratic Upstream Interpolation for Convective Kinematics). 338 It is an upstream-biased third order scheme based on an upstream-biased parabolic interpolation. 339 For example, in the $i$-direction: 335 340 \begin{equation} \label{eq:tra_adv_ubs} 336 341 \tau _u^{ubs} =\overline T ^{i+1/2}-\;\frac{1}{6} \left\{ … … 342 347 where $\tau "_i =\delta _i \left[ {\delta _{i+1/2} \left[ \tau \right]} \right]$. 343 348 344 This results in a dissipatively dominant (i.e. hyper-diffusive) truncation 345 error \citep{Shchepetkin_McWilliams_OM05}. The overall performance of 346 the advection scheme is similar to that reported in \cite{Farrow1995}. 347 It is a relatively good compromise between accuracy and smoothness. 348 Nevertheless the scheme is not \emph{positive}, meaning that false extrema are permitted, 349 but the amplitude of such are significantly reduced over the centred second 350 or fourth order method. therefore it is not recommended that it should be 351 applied to a passive tracer that requires positivity. 352 353 The intrinsic diffusion of UBS makes its use risky in the vertical direction 354 where the control of artificial diapycnal fluxes is of paramount importance \citep{Shchepetkin_McWilliams_OM05, Demange_PhD2014}. 355 Therefore the vertical flux is evaluated using either a $2^nd$ order FCT scheme 356 or a $4^th$ order COMPACT scheme (\np{nn\_cen\_v}\forcode{ = 2 or 4}). 357 358 For stability reasons (see \autoref{chap:STP}), 359 the first term in \autoref{eq:tra_adv_ubs} (which corresponds to a second order 360 centred scheme) is evaluated using the \textit{now} tracer (centred in time) 361 while the second term (which is the diffusive part of the scheme), is 362 evaluated using the \textit{before} tracer (forward in time). 363 This choice is discussed by \citet{Webb_al_JAOT98} in the context of the 364 QUICK advection scheme. UBS and QUICK schemes only differ 365 by one coefficient. Replacing 1/6 with 1/8 in \autoref{eq:tra_adv_ubs} 366 leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. 367 This option is not available through a namelist parameter, since the 368 1/6 coefficient is hard coded. Nevertheless it is quite easy to make the 369 substitution in the \mdl{traadv\_ubs} module and obtain a QUICK scheme. 349 This results in a dissipatively dominant (i.e. hyper-diffusive) truncation error 350 \citep{Shchepetkin_McWilliams_OM05}. 351 The overall performance of the advection scheme is similar to that reported in \cite{Farrow1995}. 352 It is a relatively good compromise between accuracy and smoothness. 353 Nevertheless the scheme is not \emph{positive}, meaning that false extrema are permitted, 354 but the amplitude of such are significantly reduced over the centred second or fourth order method. 355 Therefore it is not recommended that it should be applied to a passive tracer that requires positivity. 356 357 The intrinsic diffusion of UBS makes its use risky in the vertical direction where 358 the control of artificial diapycnal fluxes is of paramount importance 359 \citep{Shchepetkin_McWilliams_OM05, Demange_PhD2014}. 360 Therefore the vertical flux is evaluated using either a $2^nd$ order FCT scheme or a $4^th$ order COMPACT scheme 361 (\np{nn\_cen\_v}\forcode{ = 2 or 4}). 362 363 For stability reasons (see \autoref{chap:STP}), the first term in \autoref{eq:tra_adv_ubs} 364 (which corresponds to a second order centred scheme) 365 is evaluated using the \textit{now} tracer (centred in time) while the second term 366 (which is the diffusive part of the scheme), 367 is evaluated using the \textit{before} tracer (forward in time). 368 This choice is discussed by \citet{Webb_al_JAOT98} in the context of the QUICK advection scheme. 369 UBS and QUICK schemes only differ by one coefficient. 370 Replacing 1/6 with 1/8 in \autoref{eq:tra_adv_ubs} leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. 371 This option is not available through a namelist parameter, since the 1/6 coefficient is hard coded. 372 Nevertheless it is quite easy to make the substitution in the \mdl{traadv\_ubs} module and obtain a QUICK scheme. 370 373 371 374 Note that it is straightforward to rewrite \autoref{eq:tra_adv_ubs} as follows: … … 384 387 \end{equation} 385 388 386 \autoref{eq:traadv_ubs2} has several advantages. Firstly, it clearly reveals387 that the UBS scheme is based on the fourth order scheme to which an 388 upstream-biased diffusion term is added. Secondly, this emphasises that the 389 $4^{th}$ order part (as well as the $2^{nd}$ order part as stated above) has 390 to be evaluated at the \emph{now} time step using \autoref{eq:tra_adv_ubs}. 391 Thirdly, the diffusion term is in fact a biharmonic operator with an eddy 392 coefficient whichis simply proportional to the velocity:393 $A_u^{lm}= \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note the current version of NEMO uses 394 the computationally more efficient formulation \autoref{eq:tra_adv_ubs}.389 \autoref{eq:traadv_ubs2} has several advantages. 390 Firstly, it clearly reveals that the UBS scheme is based on the fourth order scheme to which 391 an upstream-biased diffusion term is added. 392 Secondly, this emphasises that the $4^{th}$ order part (as well as the $2^{nd}$ order part as stated above) has to 393 be evaluated at the \emph{now} time step using \autoref{eq:tra_adv_ubs}. 394 Thirdly, the diffusion term is in fact a biharmonic operator with an eddy coefficient which 395 is simply proportional to the velocity: 396 $A_u^{lm}= \frac{1}{12}\,{e_{1u}}^3\,|u|$. 397 Note the current version of NEMO uses the computationally more efficient formulation \autoref{eq:tra_adv_ubs}. 395 398 396 399 % ------------------------------------------------------------------------------------------------------------- … … 400 403 \label{subsec:TRA_adv_qck} 401 404 402 The Quadratic Upstream Interpolation for Convective Kinematics with 403 Estimated Streaming Terms (QUICKEST) scheme proposed by \citet{Leonard1979} 404 is used when \np{ln\_traadv\_qck}\forcode{ = .true.}. 405 The Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms (QUICKEST) scheme 406 proposed by \citet{Leonard1979} is used when \np{ln\_traadv\_qck}\forcode{ = .true.}. 405 407 QUICKEST implementation can be found in the \mdl{traadv\_qck} module. 406 408 407 QUICKEST is the third order Godunov scheme which is associated with the ULTIMATE QUICKEST 408 limiter \citep{Leonard1991}. It has been implemented in NEMO by G. Reffray 409 (MERCATOR-ocean) and can be found in the \mdl{traadv\_qck} module. 410 The resulting scheme is quite expensive but \emph{positive}. 411 It can be used on both active and passive tracers. 412 However, the intrinsic diffusion of QCK makes its use risky in the vertical 413 direction where the control of artificial diapycnal fluxes is of paramount importance. 414 Therefore the vertical flux is evaluated using the CEN2 scheme. 415 This no longer guarantees the positivity of the scheme. 416 The use of FCT in the vertical direction (as for the UBS case) should be implemented 417 to restore this property. 409 QUICKEST is the third order Godunov scheme which is associated with the ULTIMATE QUICKEST limiter 410 \citep{Leonard1991}. 411 It has been implemented in NEMO by G. Reffray (MERCATOR-ocean) and can be found in the \mdl{traadv\_qck} module. 412 The resulting scheme is quite expensive but \emph{positive}. 413 It can be used on both active and passive tracers. 414 However, the intrinsic diffusion of QCK makes its use risky in the vertical direction where 415 the control of artificial diapycnal fluxes is of paramount importance. 416 Therefore the vertical flux is evaluated using the CEN2 scheme. 417 This no longer guarantees the positivity of the scheme. 418 The use of FCT in the vertical direction (as for the UBS case) should be implemented to restore this property. 418 419 419 420 %%%gmcomment : Cross term are missing in the current implementation.... … … 432 433 Options are defined through the \ngn{namtra\_ldf} namelist variables. 433 434 They are regrouped in four items, allowing to specify 434 $(i)$ the type of operator used (none, laplacian, bilaplacian), 435 $(ii)$ the direction along which the operator acts (iso-level, horizontal, iso-neutral), 436 $(iii)$ some specific options related to the rotated operators ($i.e.$ non-iso-level operator), and 435 $(i)$ the type of operator used (none, laplacian, bilaplacian), 436 $(ii)$ the direction along which the operator acts (iso-level, horizontal, iso-neutral), 437 $(iii)$ some specific options related to the rotated operators ($i.e.$ non-iso-level operator), and 437 438 $(iv)$ the specification of eddy diffusivity coefficient (either constant or variable in space and time). 438 Item $(iv)$ will be described in \autoref{chap:LDF} . 439 The direction along which the operators act is defined through the slope between this direction and the iso-level surfaces. 440 The slope is computed in the \mdl{ldfslp} module and will also be described in \autoref{chap:LDF}. 441 442 The lateral diffusion of tracers is evaluated using a forward scheme, 443 $i.e.$ the tracers appearing in its expression are the \textit{before} tracers in time, 444 except for the pure vertical component that appears when a rotation tensor is used. 445 This latter component is solved implicitly together with the vertical diffusion term (see \autoref{chap:STP}). 446 When \np{ln\_traldf\_msc}\forcode{ = .true.}, a Method of Stabilizing Correction is used in which 439 Item $(iv)$ will be described in \autoref{chap:LDF}. 440 The direction along which the operators act is defined through the slope between 441 this direction and the iso-level surfaces. 442 The slope is computed in the \mdl{ldfslp} module and will also be described in \autoref{chap:LDF}. 443 444 The lateral diffusion of tracers is evaluated using a forward scheme, 445 $i.e.$ the tracers appearing in its expression are the \textit{before} tracers in time, 446 except for the pure vertical component that appears when a rotation tensor is used. 447 This latter component is solved implicitly together with the vertical diffusion term (see \autoref{chap:STP}). 448 When \np{ln\_traldf\_msc}\forcode{ = .true.}, a Method of Stabilizing Correction is used in which 447 449 the pure vertical component is split into an explicit and an implicit part \citep{Lemarie_OM2012}. 448 450 … … 456 458 Three operator options are proposed and, one and only one of them must be selected: 457 459 \begin{description} 458 \item [\np{ln\_traldf\_NONE}\forcode{ = .true.}]: no operator selected, the lateral diffusive tendency will not be 459 applied to the tracer equation. This option can be used when the selected advection scheme 460 is diffusive enough (MUSCL scheme for example). 461 \item [\np{ln\_traldf\_lap}\forcode{ = .true.}]: a laplacian operator is selected. This harmonic operator 462 takes the following expression: $\mathpzc{L}(T)=\nabla \cdot A_{ht}\;\nabla T $, 463 where the gradient operates along the selected direction (see \autoref{subsec:TRA_ldf_dir}), 464 and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see \autoref{chap:LDF}). 465 \item [\np{ln\_traldf\_blp}\forcode{ = .true.}]: a bilaplacian operator is selected. This biharmonic operator 466 takes the following expression: 467 $\mathpzc{B}=- \mathpzc{L}\left(\mathpzc{L}(T) \right) = -\nabla \cdot b\nabla \left( {\nabla \cdot b\nabla T} \right)$ 468 where the gradient operats along the selected direction, 469 and $b^2=B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$ (see \autoref{chap:LDF}). 470 In the code, the bilaplacian operator is obtained by calling the laplacian twice. 460 \item[\np{ln\_traldf\_NONE}\forcode{ = .true.}:] 461 no operator selected, the lateral diffusive tendency will not be applied to the tracer equation. 462 This option can be used when the selected advection scheme is diffusive enough (MUSCL scheme for example). 463 \item[\np{ln\_traldf\_lap}\forcode{ = .true.}:] 464 a laplacian operator is selected. 465 This harmonic operator takes the following expression: $\mathpzc{L}(T)=\nabla \cdot A_{ht}\;\nabla T $, 466 where the gradient operates along the selected direction (see \autoref{subsec:TRA_ldf_dir}), 467 and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see \autoref{chap:LDF}). 468 \item[\np{ln\_traldf\_blp}\forcode{ = .true.}]: 469 a bilaplacian operator is selected. 470 This biharmonic operator takes the following expression: 471 $\mathpzc{B}=- \mathpzc{L}\left(\mathpzc{L}(T) \right) = -\nabla \cdot b\nabla \left( {\nabla \cdot b\nabla T} \right)$ 472 where the gradient operats along the selected direction, 473 and $b^2=B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$ (see \autoref{chap:LDF}). 474 In the code, the bilaplacian operator is obtained by calling the laplacian twice. 471 475 \end{description} 472 476 473 Both laplacian and bilaplacian operators ensure the total tracer variance decrease. 474 Their primary role is to provide strong dissipation at the smallest scale supported 475 by the grid while minimizing the impact on the larger scale features. 476 The main difference between the two operators is the scale selectiveness. 477 The bilaplacian damping time ($i.e.$ its spin down time) scales like $\lambda^{-4}$ 478 for disturbances of wavelength $\lambda$ (so that short waves damped more rapidelly than long ones), 477 Both laplacian and bilaplacian operators ensure the total tracer variance decrease. 478 Their primary role is to provide strong dissipation at the smallest scale supported by the grid while 479 minimizing the impact on the larger scale features. 480 The main difference between the two operators is the scale selectiveness. 481 The bilaplacian damping time ($i.e.$ its spin down time) scales like $\lambda^{-4}$ for 482 disturbances of wavelength $\lambda$ (so that short waves damped more rapidelly than long ones), 479 483 whereas the laplacian damping time scales only like $\lambda^{-2}$. 480 484 … … 487 491 \label{subsec:TRA_ldf_dir} 488 492 489 The choice of a direction of action determines the form of operator used. 490 The operator is a simple (re-entrant) laplacian acting in the (\textbf{i},\textbf{j}) plane 491 when iso-level option is used (\np{ln\_traldf\_lev}\forcode{ = .true.}) 492 or when a horizontal ($i.e.$ geopotential) operator is demanded in \textit{z}-coordinate 493 (\np{ln\_traldf\_hor} and \np{ln\_zco} equal \forcode{.true.}). 493 The choice of a direction of action determines the form of operator used. 494 The operator is a simple (re-entrant) laplacian acting in the (\textbf{i},\textbf{j}) plane when 495 iso-level option is used (\np{ln\_traldf\_lev}\forcode{ = .true.}) or 496 when a horizontal ($i.e.$ geopotential) operator is demanded in \textit{z}-coordinate 497 (\np{ln\_traldf\_hor} and \np{ln\_zco} equal \forcode{.true.}). 494 498 The associated code can be found in the \mdl{traldf\_lap\_blp} module. 495 The operator is a rotated (re-entrant) laplacian when the direction along which it acts 496 does not coincide with the iso-level surfaces, 497 that is when standard or triad iso-neutral option is used (\np{ln\_traldf\_iso} or 498 \np{ln\_traldf\_triad} equals \forcode{.true.}, see \mdl{traldf\_iso} or \mdl{traldf\_triad} module, resp.), 499 or when a horizontal ($i.e.$ geopotential) operator is demanded in \textit{s}-coordinate 499 The operator is a rotated (re-entrant) laplacian when 500 the direction along which it acts does not coincide with the iso-level surfaces, 501 that is when standard or triad iso-neutral option is used 502 (\np{ln\_traldf\_iso} or \np{ln\_traldf\_triad} equals \forcode{.true.}, 503 see \mdl{traldf\_iso} or \mdl{traldf\_triad} module, resp.), or 504 when a horizontal ($i.e.$ geopotential) operator is demanded in \textit{s}-coordinate 500 505 (\np{ln\_traldf\_hor} and \np{ln\_sco} equal \forcode{.true.}) 501 \footnote{In this case, the standard iso-neutral operator will be automatically selected}. 502 In that case, a rotation is applied to the gradient(s) that appears in the operator 503 so thatdiffusive fluxes acts on the three spatial direction.504 505 The resulting discret form of the three operators (one iso-level and two rotated one) 506 is given inthe next two sub-sections.506 \footnote{In this case, the standard iso-neutral operator will be automatically selected}. 507 In that case, a rotation is applied to the gradient(s) that appears in the operator so that 508 diffusive fluxes acts on the three spatial direction. 509 510 The resulting discret form of the three operators (one iso-level and two rotated one) is given in 511 the next two sub-sections. 507 512 508 513 … … 519 524 + \delta _{j}\left[ A_v^{lT} \; \frac{e_{1v}\,e_{3v}}{e_{2v}} \;\delta _{j+1/2} [T] \right] \;\right) 520 525 \end{equation} 521 where $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$-cells 522 and where zero diffusive fluxes is assumed across solid boundaries, 523 first (and third in bilaplacian case) horizontal tracer derivative are masked. 524 It is implemented in the \rou{traldf\_lap} subroutine found in the \mdl{traldf\_lap} module. 525 The module also contains \rou{traldf\_blp}, the subroutine calling twice \rou{traldf\_lap} 526 in order tocompute the iso-level bilaplacian operator.527 528 It is a \emph{horizontal} operator ($i.e.$ acting along geopotential surfaces) in the $z$-coordinate529 with or without partial steps, but is simply an iso-level operator in the $s$-coordinate. 530 It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp}\forcode{ = .true.}, 531 we have \np{ln\_traldf\_lev}\forcode{ = .true.} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}\forcode{ = .true.}. 532 In both cases, it significantly contributes to diapycnal mixing. 526 where $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$-cells and 527 where zero diffusive fluxes is assumed across solid boundaries, 528 first (and third in bilaplacian case) horizontal tracer derivative are masked. 529 It is implemented in the \rou{traldf\_lap} subroutine found in the \mdl{traldf\_lap} module. 530 The module also contains \rou{traldf\_blp}, the subroutine calling twice \rou{traldf\_lap} in order to 531 compute the iso-level bilaplacian operator. 532 533 It is a \emph{horizontal} operator ($i.e.$ acting along geopotential surfaces) in 534 the $z$-coordinate with or without partial steps, but is simply an iso-level operator in the $s$-coordinate. 535 It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp}\forcode{ = .true.}, 536 we have \np{ln\_traldf\_lev}\forcode{ = .true.} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}\forcode{ = .true.}. 537 In both cases, it significantly contributes to diapycnal mixing. 533 538 It is therefore never recommended, even when using it in the bilaplacian case. 534 539 535 Note that in the partial step $z$-coordinate (\np{ln\_zps}\forcode{ = .true.}), tracers in horizontally 536 adjacent cells are located at different depths in the vicinity of the bottom. 537 In this case, horizontal derivatives in (\autoref{eq:tra_ldf_lap}) at the bottom level 538 require a specific treatment. They are calculated in the \mdl{zpshde} module, 539 described in \autoref{sec:TRA_zpshde}. 540 Note that in the partial step $z$-coordinate (\np{ln\_zps}\forcode{ = .true.}), 541 tracers in horizontally adjacent cells are located at different depths in the vicinity of the bottom. 542 In this case, horizontal derivatives in (\autoref{eq:tra_ldf_lap}) at the bottom level require a specific treatment. 543 They are calculated in the \mdl{zpshde} module, described in \autoref{sec:TRA_zpshde}. 540 544 541 545 … … 550 554 \subsubsection{Standard rotated (bi-)laplacian operator (\protect\mdl{traldf\_iso})} 551 555 \label{subsec:TRA_ldf_iso} 552 The general form of the second order lateral tracer subgrid scale physics 553 (\autoref{eq:PE_zdf})takes the following semi-discrete space form in $z$- and $s$-coordinates:556 The general form of the second order lateral tracer subgrid scale physics (\autoref{eq:PE_zdf}) 557 takes the following semi-discrete space form in $z$- and $s$-coordinates: 554 558 \begin{equation} \label{eq:tra_ldf_iso} 555 559 \begin{split} … … 569 573 & \left. {\left. { \qquad \qquad \ \ \ \left. { 570 574 +\;\frac{e_{1w}\,e_{2w}}{e_{3w}} \,\left( r_{1w}^2 + r_{2w}^2 \right) 571 \,\delta_{k+1/2} [T] } \right) } \right] \quad } \right\} 572 \end{split} 573 \end{equation} 574 where $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$-cells, 575 $r_1$ and $r_2$ are the slopes between the surface of computation 576 ($z$- or $s$-surfaces) and the surface along which the diffusion operator 577 acts ($i.e.$ horizontal or iso-neutral surfaces). It is thus used when, 578 in addition to \np{ln\_traldf\_lap}\forcode{ = .true.}, we have \np{ln\_traldf\_iso}\forcode{ = .true.}, 579 or both \np{ln\_traldf\_hor}\forcode{ = .true.} and \np{ln\_zco}\forcode{ = .true.}. The way these 580 slopes are evaluated is given in \autoref{sec:LDF_slp}. At the surface, bottom 581 and lateral boundaries, the turbulent fluxes of heat and salt are set to zero 582 using the mask technique (see \autoref{sec:LBC_coast}). 583 584 The operator in \autoref{eq:tra_ldf_iso} involves both lateral and vertical 585 derivatives. For numerical stability, the vertical second derivative must 586 be solved using the same implicit time scheme as that used in the vertical 587 physics (see \autoref{sec:TRA_zdf}). For computer efficiency reasons, this term 588 is not computed in the \mdl{traldf\_iso} module, but in the \mdl{trazdf} module 589 where, if iso-neutral mixing is used, the vertical mixing coefficient is simply 590 increased by $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$. 591 592 This formulation conserves the tracer but does not ensure the decrease 593 of the tracer variance. Nevertheless the treatment performed on the slopes 594 (see \autoref{chap:LDF}) allows the model to run safely without any additional 595 background horizontal diffusion \citep{Guilyardi_al_CD01}. 596 597 Note that in the partial step $z$-coordinate (\np{ln\_zps}\forcode{ = .true.}), the horizontal derivatives 598 at the bottom level in \autoref{eq:tra_ldf_iso} require a specific treatment. 575 \,\delta_{k+1/2} [T] } \right) } \right] \quad } \right\} 576 \end{split} 577 \end{equation} 578 where $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$-cells, 579 $r_1$ and $r_2$ are the slopes between the surface of computation ($z$- or $s$-surfaces) and 580 the surface along which the diffusion operator acts ($i.e.$ horizontal or iso-neutral surfaces). 581 It is thus used when, in addition to \np{ln\_traldf\_lap}\forcode{ = .true.}, 582 we have \np{ln\_traldf\_iso}\forcode{ = .true.}, 583 or both \np{ln\_traldf\_hor}\forcode{ = .true.} and \np{ln\_zco}\forcode{ = .true.}. 584 The way these slopes are evaluated is given in \autoref{sec:LDF_slp}. 585 At the surface, bottom and lateral boundaries, the turbulent fluxes of heat and salt are set to zero using 586 the mask technique (see \autoref{sec:LBC_coast}). 587 588 The operator in \autoref{eq:tra_ldf_iso} involves both lateral and vertical derivatives. 589 For numerical stability, the vertical second derivative must be solved using the same implicit time scheme as that 590 used in the vertical physics (see \autoref{sec:TRA_zdf}). 591 For computer efficiency reasons, this term is not computed in the \mdl{traldf\_iso} module, 592 but in the \mdl{trazdf} module where, if iso-neutral mixing is used, 593 the vertical mixing coefficient is simply increased by 594 $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$. 595 596 This formulation conserves the tracer but does not ensure the decrease of the tracer variance. 597 Nevertheless the treatment performed on the slopes (see \autoref{chap:LDF}) allows the model to run safely without 598 any additional background horizontal diffusion \citep{Guilyardi_al_CD01}. 599 600 Note that in the partial step $z$-coordinate (\np{ln\_zps}\forcode{ = .true.}), 601 the horizontal derivatives at the bottom level in \autoref{eq:tra_ldf_iso} require a specific treatment. 599 602 They are calculated in module zpshde, described in \autoref{sec:TRA_zpshde}. 600 603 … … 604 607 \label{subsec:TRA_ldf_triad} 605 608 606 If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}\forcode{ = .true.} ; see \autoref{apdx:triad}) 607 608 An alternative scheme developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases 609 is also available in \NEMO (\np{ln\_traldf\_grif}\forcode{ = .true.}). A complete description of 610 the algorithm is given in \autoref{apdx:triad}. 611 612 The lateral fourth order bilaplacian operator on tracers is obtained by 613 applying (\autoref{eq:tra_ldf_lap}) twice. The operator requires an additional assumption 614 on boundary conditions: both first and third derivative terms normal to the 615 coast are set to zero. 616 617 The lateral fourth order operator formulation on tracers is obtained by 618 applying (\autoref{eq:tra_ldf_iso}) twice. It requires an additional assumption 619 on boundary conditions: first and third derivative terms normal to the 620 coast, normal to the bottom and normal to the surface are set to zero. 609 If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}\forcode{ = .true.}; see \autoref{apdx:triad}) 610 611 An alternative scheme developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases 612 is also available in \NEMO (\np{ln\_traldf\_grif}\forcode{ = .true.}). 613 A complete description of the algorithm is given in \autoref{apdx:triad}. 614 615 The lateral fourth order bilaplacian operator on tracers is obtained by applying (\autoref{eq:tra_ldf_lap}) twice. 616 The operator requires an additional assumption on boundary conditions: 617 both first and third derivative terms normal to the coast are set to zero. 618 619 The lateral fourth order operator formulation on tracers is obtained by applying (\autoref{eq:tra_ldf_iso}) twice. 620 It requires an additional assumption on boundary conditions: 621 first and third derivative terms normal to the coast, 622 normal to the bottom and normal to the surface are set to zero. 621 623 622 624 %&& Option for the rotated operators … … 631 633 \np{ln\_triad\_iso} = pure horizontal mixing in ML (triad only) 632 634 633 \np{rn\_sw\_triad} =1 switching triad ; =0 all 4 triads used (triad only) 635 \np{rn\_sw\_triad} =1 switching triad; 636 =0 all 4 triads used (triad only) 634 637 635 638 \np{ln\_botmix\_triad} = lateral mixing on bottom (triad only) … … 646 649 647 650 Options are defined through the \ngn{namzdf} namelist variables. 648 The formulation of the vertical subgrid scale tracer physics is the same 649 for all the vertical coordinates, and is based on a laplacian operator. 650 The vertical diffusion operator given by (\autoref{eq:PE_zdf}) takes the 651 following semi-discrete space form: 651 The formulation of the vertical subgrid scale tracer physics is the same for all the vertical coordinates, 652 and is based on a laplacian operator. 653 The vertical diffusion operator given by (\autoref{eq:PE_zdf}) takes the following semi-discrete space form: 652 654 \begin{equation} \label{eq:tra_zdf} 653 655 \begin{split} … … 657 659 \end{split} 658 660 \end{equation} 659 where $A_w^{vT}$ and $A_w^{vS}$ are the vertical eddy diffusivity 660 coefficients on temperature and salinity, respectively. Generally, 661 $A_w^{vT}=A_w^{vS}$ except when double diffusive mixing is 662 parameterised ($i.e.$ \key{zdfddm} is defined). The way these coefficients 663 are evaluated is given in \autoref{chap:ZDF} (ZDF). Furthermore, when 664 iso-neutral mixing is used, both mixing coefficients are increased 665 by $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$ 666 to account for the vertical second derivative of \autoref{eq:tra_ldf_iso}. 667 668 At the surface and bottom boundaries, the turbulent fluxes of 669 heat and salt must be specified. At the surface they are prescribed 670 from the surface forcing and added in a dedicated routine (see \autoref{subsec:TRA_sbc}), 671 whilst at the bottom they are set to zero for heat and salt unless 672 a geothermal flux forcing is prescribed as a bottom boundary 673 condition (see \autoref{subsec:TRA_bbc}). 674 675 The large eddy coefficient found in the mixed layer together with high 676 vertical resolution implies that in the case of explicit time stepping 677 (\np{ln\_zdfexp}\forcode{ = .true.}) there would be too restrictive a constraint on 678 the time step. Therefore, the default implicit time stepping is preferred 679 for the vertical diffusion since it overcomes the stability constraint. 680 A forward time differencing scheme (\np{ln\_zdfexp}\forcode{ = .true.}) using a time 681 splitting technique (\np{nn\_zdfexp} $> 1$) is provided as an alternative. 682 Namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both 683 tracers and dynamics. 661 where $A_w^{vT}$ and $A_w^{vS}$ are the vertical eddy diffusivity coefficients on temperature and salinity, 662 respectively. 663 Generally, $A_w^{vT}=A_w^{vS}$ except when double diffusive mixing is parameterised ($i.e.$ \key{zdfddm} is defined). 664 The way these coefficients are evaluated is given in \autoref{chap:ZDF} (ZDF). 665 Furthermore, when iso-neutral mixing is used, both mixing coefficients are increased by 666 $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$ to account for 667 the vertical second derivative of \autoref{eq:tra_ldf_iso}. 668 669 At the surface and bottom boundaries, the turbulent fluxes of heat and salt must be specified. 670 At the surface they are prescribed from the surface forcing and added in a dedicated routine 671 (see \autoref{subsec:TRA_sbc}), whilst at the bottom they are set to zero for heat and salt unless 672 a geothermal flux forcing is prescribed as a bottom boundary condition (see \autoref{subsec:TRA_bbc}). 673 674 The large eddy coefficient found in the mixed layer together with high vertical resolution implies that 675 in the case of explicit time stepping (\np{ln\_zdfexp}\forcode{ = .true.}) 676 there would be too restrictive a constraint on the time step. 677 Therefore, the default implicit time stepping is preferred for the vertical diffusion since 678 it overcomes the stability constraint. 679 A forward time differencing scheme (\np{ln\_zdfexp}\forcode{ = .true.}) using 680 a time splitting technique (\np{nn\_zdfexp} $> 1$) is provided as an alternative. 681 Namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics. 684 682 685 683 % ================================================================ … … 695 693 \label{subsec:TRA_sbc} 696 694 697 The surface boundary condition for tracers is implemented in a separate 698 module (\mdl{trasbc}) instead of entering as a boundary condition on the vertical 699 diffusion operator (as in the case of momentum). This has been found to 700 enhance readability of the code. The two formulations are completely 701 equivalent; the forcing terms in trasbc are the surface fluxes divided by702 the thickness of the top model layer. 703 704 Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components 705 ($i.e.$ atmosphere, sea-ice, land), the change in the heat and salt content of the surface layer 706 of the ocean is due both to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$) 707 and to the heat and salt content of the mass exchange. They are both included directly in $Q_{ns}$, 708 the surface heat flux,and $F_{salt}$, the surface salt flux (see \autoref{chap:SBC} for further details).695 The surface boundary condition for tracers is implemented in a separate module (\mdl{trasbc}) instead of 696 entering as a boundary condition on the vertical diffusion operator (as in the case of momentum). 697 This has been found to enhance readability of the code. 698 The two formulations are completely equivalent; 699 the forcing terms in trasbc are the surface fluxes divided by the thickness of the top model layer. 700 701 Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components 702 ($i.e.$ atmosphere, sea-ice, land), the change in the heat and salt content of the surface layer of the ocean is due 703 both to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$) and 704 to the heat and salt content of the mass exchange. 705 They are both included directly in $Q_{ns}$, the surface heat flux, 706 and $F_{salt}$, the surface salt flux (see \autoref{chap:SBC} for further details). 709 707 By doing this, the forcing formulation is the same for any tracer (including temperature and salinity). 710 708 711 The surface module (\mdl{sbcmod}, see \autoref{chap:SBC}) provides the following 712 forcing fields (used on tracers): 713 714 $\bullet$ $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface 709 The surface module (\mdl{sbcmod}, see \autoref{chap:SBC}) provides the following forcing fields (used on tracers): 710 711 $\bullet$ $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface 715 712 (i.e. the difference between the total surface heat flux and the fraction of the short wave flux that 716 penetrates into the water column, see \autoref{subsec:TRA_qsr}) plus the heat content associated with717 of the mass exchange with the atmosphere and lands.713 penetrates into the water column, see \autoref{subsec:TRA_qsr}) 714 plus the heat content associated with of the mass exchange with the atmosphere and lands. 718 715 719 716 $\bullet$ $\textit{sfx}$, the salt flux resulting from ice-ocean mass exchange (freezing, melting, ridging...) 720 717 721 $\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation) 722 andpossibly with the sea-ice and ice-shelves.723 724 $\bullet$ \textit{rnf}, the mass flux associated with runoff 718 $\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation) and 719 possibly with the sea-ice and ice-shelves. 720 721 $\bullet$ \textit{rnf}, the mass flux associated with runoff 725 722 (see \autoref{sec:SBC_rnf} for further detail of how it acts on temperature and salinity tendencies) 726 723 727 $\bullet$ \textit{fwfisf}, the mass flux associated with ice shelf melt, 724 $\bullet$ \textit{fwfisf}, the mass flux associated with ice shelf melt, 728 725 (see \autoref{sec:SBC_isf} for further details on how the ice shelf melt is computed and applied). 729 726 … … 735 732 \end{aligned} 736 733 \end{equation} 737 where $\overline{x }^t$ means that $x$ is averaged over two consecutive time steps 738 ($t-\rdt/2$ and $t+\rdt/2$). Such time averaging prevents the 739 divergence of odd and even time step (see \autoref{chap:STP}). 740 741 In the linear free surface case (\np{ln\_linssh}\forcode{ = .true.}), 742 an additional term has to be added on both temperature and salinity. 743 On temperature, this term remove the heat content associated with mass exchange 744 that has been added to $Q_{ns}$. On salinity, this term mimics the concentration/dilution effect that 745 would have resulted from a change in the volume of the first level. 734 where $\overline{x }^t$ means that $x$ is averaged over two consecutive time steps ($t-\rdt/2$ and $t+\rdt/2$). 735 Such time averaging prevents the divergence of odd and even time step (see \autoref{chap:STP}). 736 737 In the linear free surface case (\np{ln\_linssh}\forcode{ = .true.}), an additional term has to be added on 738 both temperature and salinity. 739 On temperature, this term remove the heat content associated with mass exchange that has been added to $Q_{ns}$. 740 On salinity, this term mimics the concentration/dilution effect that would have resulted from a change in 741 the volume of the first level. 746 742 The resulting surface boundary condition is applied as follows: 747 743 \begin{equation} \label{eq:tra_sbc_lin} … … 754 750 \end{aligned} 755 751 \end{equation} 756 Note that an exact conservation of heat and salt content is only achieved with non-linear free surface. 757 In the linear free surface case, there is a small imbalance. The imbalance is larger758 than the imbalance associated with the Asselin time filter \citep{Leclair_Madec_OM09}. 752 Note that an exact conservation of heat and salt content is only achieved with non-linear free surface. 753 In the linear free surface case, there is a small imbalance. 754 The imbalance is larger than the imbalance associated with the Asselin time filter \citep{Leclair_Madec_OM09}. 759 755 This is the reason why the modified filter is not applied in the linear free surface case (see \autoref{chap:STP}). 760 756 … … 770 766 771 767 Options are defined through the \ngn{namtra\_qsr} namelist variables. 772 When the penetrative solar radiation option is used (\np{ln\_flxqsr}\forcode{ = .true.}), 773 the solar radiation penetrates the top few tens of meters of the ocean. If it is not used 774 (\np{ln\_flxqsr}\forcode{ = .false.}) all the heat flux is absorbed in the first ocean level. 775 Thus, in the former case a term is added to the time evolution equation of 776 temperature \autoref{eq:PE_tra_T} and the surface boundary condition is 777 modified to take into account only the non-penetrative part of the surface 768 When the penetrative solar radiation option is used (\np{ln\_flxqsr}\forcode{ = .true.}), 769 the solar radiation penetrates the top few tens of meters of the ocean. 770 If it is not used (\np{ln\_flxqsr}\forcode{ = .false.}) all the heat flux is absorbed in the first ocean level. 771 Thus, in the former case a term is added to the time evolution equation of temperature \autoref{eq:PE_tra_T} and 772 the surface boundary condition is modified to take into account only the non-penetrative part of the surface 778 773 heat flux: 779 774 \begin{equation} \label{eq:PE_qsr} … … 783 778 \end{split} 784 779 \end{equation} 785 where $Q_{sr}$ is the penetrative part of the surface heat flux ($i.e.$ the shortwave radiation) 786 and $I$ is the downward irradiance ($\left. I \right|_{z=\eta}=Q_{sr}$). 780 where $Q_{sr}$ is the penetrative part of the surface heat flux ($i.e.$ the shortwave radiation) and 781 $I$ is the downward irradiance ($\left. I \right|_{z=\eta}=Q_{sr}$). 787 782 The additional term in \autoref{eq:PE_qsr} is discretized as follows: 788 783 \begin{equation} \label{eq:tra_qsr} … … 790 785 \end{equation} 791 786 792 The shortwave radiation, $Q_{sr}$, consists of energy distributed across a wide spectral range. 793 The ocean is strongly absorbing for wavelengths longer than 700~nm and these 794 wavelengths contribute to heating the upper few tens of centimetres. The fraction of $Q_{sr}$ 795 that resides in these almost non-penetrative wavebands, $R$, is $\sim 58\%$ (specified 796 through namelist parameter \np{rn\_abs}). It is assumed to penetrate the ocean 797 with a decreasing exponential profile, with an e-folding depth scale, $\xi_0$, 787 The shortwave radiation, $Q_{sr}$, consists of energy distributed across a wide spectral range. 788 The ocean is strongly absorbing for wavelengths longer than 700~nm and these wavelengths contribute to 789 heating the upper few tens of centimetres. 790 The fraction of $Q_{sr}$ that resides in these almost non-penetrative wavebands, $R$, is $\sim 58\%$ 791 (specified through namelist parameter \np{rn\_abs}). 792 It is assumed to penetrate the ocean with a decreasing exponential profile, with an e-folding depth scale, $\xi_0$, 798 793 of a few tens of centimetres (typically $\xi_0=0.35~m$ set as \np{rn\_si0} in the \ngn{namtra\_qsr} namelist). 799 For shorter wavelengths (400-700~nm), the ocean is more transparent, and solar energy 800 propagates to larger depths where it contributes to 801 local heating. 802 The way this second part of the solar energy penetrates into the ocean depends on 803 which formulation is chosen. In the simple 2-waveband light penetration scheme (\np{ln\_qsr\_2bd}\forcode{ = .true.}) 804 a chlorophyll-independent monochromatic formulation is chosen for the shorter wavelengths, 794 For shorter wavelengths (400-700~nm), the ocean is more transparent, and solar energy propagates to 795 larger depths where it contributes to local heating. 796 The way this second part of the solar energy penetrates into the ocean depends on which formulation is chosen. 797 In the simple 2-waveband light penetration scheme (\np{ln\_qsr\_2bd}\forcode{ = .true.}) 798 a chlorophyll-independent monochromatic formulation is chosen for the shorter wavelengths, 805 799 leading to the following expression \citep{Paulson1977}: 806 800 \begin{equation} \label{eq:traqsr_iradiance} 807 801 I(z) = Q_{sr} \left[Re^{-z / \xi_0} + \left( 1-R\right) e^{-z / \xi_1} \right] 808 802 \end{equation} 809 where $\xi_1$ is the second extinction length scale associated with the shorter wavelengths. 810 It is usually chosen to be 23~m by setting the \np{rn\_si0} namelist parameter. 811 The set of default values ($\xi_0$, $\xi_1$, $R$) corresponds to a Type I water in 812 Jerlov's (1968) classification (oligotrophic waters). 813 814 Such assumptions have been shown to provide a very crude and simplistic 815 representation of observed light penetration profiles (\cite{Morel_JGR88}, see also 816 \autoref{fig:traqsr_irradiance}). Light absorption in the ocean depends on 817 particle concentration and is spectrally selective. \cite{Morel_JGR88} has shown 818 that an accurate representation of light penetration can be provided by a 61 waveband 819 formulation. Unfortunately, such a model is very computationally expensive. 820 Thus, \cite{Lengaigne_al_CD07} have constructed a simplified version of this 821 formulation in which visible light is split into three wavebands: blue (400-500 nm), 822 green (500-600 nm) and red (600-700nm). For each wave-band, the chlorophyll-dependent 823 attenuation coefficient is fitted to the coefficients computed from the full spectral model 824 of \cite{Morel_JGR88} (as modified by \cite{Morel_Maritorena_JGR01}), assuming 825 the same power-law relationship. As shown in \autoref{fig:traqsr_irradiance}, 826 this formulation, called RGB (Red-Green-Blue), reproduces quite closely 827 the light penetration profiles predicted by the full spectal model, but with much greater 828 computational efficiency. The 2-bands formulation does not reproduce the full model very well. 829 830 The RGB formulation is used when \np{ln\_qsr\_rgb}\forcode{ = .true.}. The RGB attenuation coefficients 831 ($i.e.$ the inverses of the extinction length scales) are tabulated over 61 nonuniform 832 chlorophyll classes ranging from 0.01 to 10 g.Chl/L (see the routine \rou{trc\_oce\_rgb} 833 in \mdl{trc\_oce} module). Four types of chlorophyll can be chosen in the RGB formulation: 803 where $\xi_1$ is the second extinction length scale associated with the shorter wavelengths. 804 It is usually chosen to be 23~m by setting the \np{rn\_si0} namelist parameter. 805 The set of default values ($\xi_0$, $\xi_1$, $R$) corresponds to a Type I water in Jerlov's (1968) classification 806 (oligotrophic waters). 807 808 Such assumptions have been shown to provide a very crude and simplistic representation of 809 observed light penetration profiles (\cite{Morel_JGR88}, see also \autoref{fig:traqsr_irradiance}). 810 Light absorption in the ocean depends on particle concentration and is spectrally selective. 811 \cite{Morel_JGR88} has shown that an accurate representation of light penetration can be provided by 812 a 61 waveband formulation. 813 Unfortunately, such a model is very computationally expensive. 814 Thus, \cite{Lengaigne_al_CD07} have constructed a simplified version of this formulation in which 815 visible light is split into three wavebands: blue (400-500 nm), green (500-600 nm) and red (600-700nm). 816 For each wave-band, the chlorophyll-dependent attenuation coefficient is fitted to the coefficients computed from 817 the full spectral model of \cite{Morel_JGR88} (as modified by \cite{Morel_Maritorena_JGR01}), 818 assuming the same power-law relationship. 819 As shown in \autoref{fig:traqsr_irradiance}, this formulation, called RGB (Red-Green-Blue), 820 reproduces quite closely the light penetration profiles predicted by the full spectal model, 821 but with much greater computational efficiency. 822 The 2-bands formulation does not reproduce the full model very well. 823 824 The RGB formulation is used when \np{ln\_qsr\_rgb}\forcode{ = .true.}. 825 The RGB attenuation coefficients ($i.e.$ the inverses of the extinction length scales) are tabulated over 826 61 nonuniform chlorophyll classes ranging from 0.01 to 10 g.Chl/L 827 (see the routine \rou{trc\_oce\_rgb} in \mdl{trc\_oce} module). 828 Four types of chlorophyll can be chosen in the RGB formulation: 834 829 \begin{description} 835 \item[\np{nn\_chdta}\forcode{ = 0}] 836 a constant 0.05 g.Chl/L value everywhere ;837 \item[\np{nn\_chdta}\forcode{ = 1}] 838 an observed time varying chlorophyll deduced from satellite surface ocean color measurement 839 spread uniformly in the vertical direction ; 840 \item[\np{nn\_chdta}\forcode{ = 2}] 841 same as previous case except that a vertical profile of chlorophyl is used. 842 Following \cite{Morel_Berthon_LO89}, the profile is computed from the local surface chlorophyll value;843 \item[\np{ln\_qsr\_bio}\forcode{ = .true.}] 844 simulated time varying chlorophyll by TOP biogeochemical model. 845 In this case, the RGB formulation is used to calculate both the phytoplankton 846 light limitation inPISCES or LOBSTER and the oceanic heating rate.830 \item[\np{nn\_chdta}\forcode{ = 0}] 831 a constant 0.05 g.Chl/L value everywhere ; 832 \item[\np{nn\_chdta}\forcode{ = 1}] 833 an observed time varying chlorophyll deduced from satellite surface ocean color measurement spread uniformly in 834 the vertical direction; 835 \item[\np{nn\_chdta}\forcode{ = 2}] 836 same as previous case except that a vertical profile of chlorophyl is used. 837 Following \cite{Morel_Berthon_LO89}, the profile is computed from the local surface chlorophyll value; 838 \item[\np{ln\_qsr\_bio}\forcode{ = .true.}] 839 simulated time varying chlorophyll by TOP biogeochemical model. 840 In this case, the RGB formulation is used to calculate both the phytoplankton light limitation in 841 PISCES or LOBSTER and the oceanic heating rate. 847 842 \end{description} 848 The trend in \autoref{eq:tra_qsr} associated with the penetration of the solar radiation 849 is added to the temperature trend, and the surface heat flux is modified in routine \mdl{traqsr}. 850 851 When the $z$-coordinate is preferred to the $s$-coordinate, the depth of $w-$levels does852 not significantly vary with location. The level at which the light has been totally 853 absorbed ($i.e.$ it is less than the computer precision) is computed once, 854 and the trend associated with the penetration of the solar radiation is only added down to that level. 855 Finally, note that when the ocean is shallow ($<$ 200~m), part of the 856 solar radiation can reach the ocean floor. In this case, we have 857 chosen that all remaining radiation is absorbed in the last ocean 858 level($i.e.$ $I$ is masked).843 The trend in \autoref{eq:tra_qsr} associated with the penetration of the solar radiation is added to 844 the temperature trend, and the surface heat flux is modified in routine \mdl{traqsr}. 845 846 When the $z$-coordinate is preferred to the $s$-coordinate, 847 the depth of $w-$levels does not significantly vary with location. 848 The level at which the light has been totally absorbed 849 ($i.e.$ it is less than the computer precision) is computed once, 850 and the trend associated with the penetration of the solar radiation is only added down to that level. 851 Finally, note that when the ocean is shallow ($<$ 200~m), part of the solar radiation can reach the ocean floor. 852 In this case, we have chosen that all remaining radiation is absorbed in the last ocean level 853 ($i.e.$ $I$ is masked). 859 854 860 855 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 861 \begin{figure}[!t] \begin{center} 862 \includegraphics[width=1.0\textwidth]{Fig_TRA_Irradiance} 863 \caption{ \protect\label{fig:traqsr_irradiance} 864 Penetration profile of the downward solar irradiance calculated by four models. 865 Two waveband chlorophyll-independent formulation (blue), a chlorophyll-dependent 866 monochromatic formulation (green), 4 waveband RGB formulation (red), 867 61 waveband Morel (1988) formulation (black) for a chlorophyll concentration of 868 (a) Chl=0.05 mg/m$^3$ and (b) Chl=0.5 mg/m$^3$. From \citet{Lengaigne_al_CD07}.} 869 \end{center} \end{figure} 856 \begin{figure}[!t] 857 \begin{center} 858 \includegraphics[width=1.0\textwidth]{Fig_TRA_Irradiance} 859 \caption{ \protect\label{fig:traqsr_irradiance} 860 Penetration profile of the downward solar irradiance calculated by four models. 861 Two waveband chlorophyll-independent formulation (blue), 862 a chlorophyll-dependent monochromatic formulation (green), 863 4 waveband RGB formulation (red), 864 61 waveband Morel (1988) formulation (black) for a chlorophyll concentration of 865 (a) Chl=0.05 mg/m$^3$ and (b) Chl=0.5 mg/m$^3$. 866 From \citet{Lengaigne_al_CD07}. 867 } 868 \end{center} 869 \end{figure} 870 870 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 871 871 … … 880 880 %-------------------------------------------------------------------------------------------------------------- 881 881 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 882 \begin{figure}[!t] \begin{center} 883 \includegraphics[width=1.0\textwidth]{Fig_TRA_geoth} 884 \caption{ \protect\label{fig:geothermal} 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_Nat92}.} 887 \end{center} \end{figure} 882 \begin{figure}[!t] 883 \begin{center} 884 \includegraphics[width=1.0\textwidth]{Fig_TRA_geoth} 885 \caption{ \protect\label{fig:geothermal} 886 Geothermal Heat flux (in $mW.m^{-2}$) used by \cite{Emile-Geay_Madec_OS09}. 887 It is inferred from the age of the sea floor and the formulae of \citet{Stein_Stein_Nat92}. 888 } 889 \end{center} 890 \end{figure} 888 891 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 889 892 890 Usually it is assumed that there is no exchange of heat or salt through 891 the ocean bottom, $i.e.$ a no flux boundary condition is applied on active 892 tracers at the bottom. This is the default option in \NEMO, and it is 893 implemented using the masking technique. However, there is a 894 non-zero heat flux across the seafloor that is associated with solid 895 earth cooling. This flux is weak compared to surface fluxes (a mean 896 global value of $\sim0.1\;W/m^2$ \citep{Stein_Stein_Nat92}), but it warms 897 systematically the ocean and acts on the densest water masses. 898 Taking this flux into account in a global ocean model increases 899 the deepest overturning cell ($i.e.$ the one associated with the Antarctic 900 Bottom Water) by a few Sverdrups \citep{Emile-Geay_Madec_OS09}. 893 Usually it is assumed that there is no exchange of heat or salt through the ocean bottom, 894 $i.e.$ a no flux boundary condition is applied on active tracers at the bottom. 895 This is the default option in \NEMO, and it is implemented using the masking technique. 896 However, there is a non-zero heat flux across the seafloor that is associated with solid earth cooling. 897 This flux is weak compared to surface fluxes (a mean global value of $\sim0.1\;W/m^2$ \citep{Stein_Stein_Nat92}), 898 but it warms systematically the ocean and acts on the densest water masses. 899 Taking this flux into account in a global ocean model increases the deepest overturning cell 900 ($i.e.$ the one associated with the Antarctic Bottom Water) by a few Sverdrups \citep{Emile-Geay_Madec_OS09}. 901 901 902 902 Options are defined through the \ngn{namtra\_bbc} namelist variables. 903 The presence of geothermal heating is controlled by setting the namelist 904 parameter \np{ln\_trabbc} to true. Then, when \np{nn\_geoflx} is set to 1, 905 a constant geothermal heating is introduced whose value is given by the 906 \np{nn\_geoflx\_cst}, which is also a namelist parameter. 907 When \np{nn\_geoflx} is set to 2, a spatially varying geothermal heat flux is 908 introduced which is provided in the \ifile{geothermal\_heating} NetCDF file 909 (\autoref{fig:geothermal}) \citep{Emile-Geay_Madec_OS09}. 903 The presence of geothermal heating is controlled by setting the namelist parameter \np{ln\_trabbc} to true. 904 Then, when \np{nn\_geoflx} is set to 1, a constant geothermal heating is introduced whose value is given by 905 the \np{nn\_geoflx\_cst}, which is also a namelist parameter. 906 When \np{nn\_geoflx} is set to 2, a spatially varying geothermal heat flux is introduced which is provided in 907 the \ifile{geothermal\_heating} NetCDF file (\autoref{fig:geothermal}) \citep{Emile-Geay_Madec_OS09}. 910 908 911 909 % ================================================================ … … 920 918 921 919 Options are defined through the \ngn{nambbl} namelist variables. 922 In a $z$-coordinate configuration, the bottom topography is represented by a 923 series of discrete steps. This is not adequate to represent gravity driven 924 downslope flows. Such flows arise either downstream of sills such as the Strait of 925 Gibraltar or Denmark Strait, where dense water formed in marginal seas flows 926 into a basin filled with less dense water, or along the continental slope when dense 927 water masses are formed on a continental shelf. The amount of entrainment 928 that occurs in these gravity plumes is critical in determining the density 929 and volume flux of the densest waters of the ocean, such as Antarctic Bottom Water, 930 or North Atlantic Deep Water. $z$-coordinate models tend to overestimate the 931 entrainment, because the gravity flow is mixed vertically by convection 932 as it goes ''downstairs'' following the step topography, sometimes over a thickness 933 much larger than the thickness of the observed gravity plume. A similar problem 934 occurs in the $s$-coordinate when the thickness of the bottom level varies rapidly 935 downstream of a sill \citep{Willebrand_al_PO01}, and the thickness 936 of the plume is not resolved. 937 938 The idea of the bottom boundary layer (BBL) parameterisation, first introduced by 939 \citet{Beckmann_Doscher1997}, is to allow a direct communication between 940 two adjacent bottom cells at different levels, whenever the densest water is 941 located above the less dense water. The communication can be by a diffusive flux 942 (diffusive BBL), an advective flux (advective BBL), or both. In the current 943 implementation of the BBL, only the tracers are modified, not the velocities. 944 Furthermore, it only connects ocean bottom cells, and therefore does not include 945 all the improvements introduced by \citet{Campin_Goosse_Tel99}. 920 In a $z$-coordinate configuration, the bottom topography is represented by a series of discrete steps. 921 This is not adequate to represent gravity driven downslope flows. 922 Such flows arise either downstream of sills such as the Strait of Gibraltar or Denmark Strait, 923 where dense water formed in marginal seas flows into a basin filled with less dense water, 924 or along the continental slope when dense water masses are formed on a continental shelf. 925 The amount of entrainment that occurs in these gravity plumes is critical in determining the density and 926 volume flux of the densest waters of the ocean, such as Antarctic Bottom Water, or North Atlantic Deep Water. 927 $z$-coordinate models tend to overestimate the entrainment, 928 because the gravity flow is mixed vertically by convection as it goes ''downstairs'' following the step topography, 929 sometimes over a thickness much larger than the thickness of the observed gravity plume. 930 A similar problem occurs in the $s$-coordinate when the thickness of the bottom level varies rapidly downstream of 931 a sill \citep{Willebrand_al_PO01}, and the thickness of the plume is not resolved. 932 933 The idea of the bottom boundary layer (BBL) parameterisation, first introduced by \citet{Beckmann_Doscher1997}, 934 is to allow a direct communication between two adjacent bottom cells at different levels, 935 whenever the densest water is located above the less dense water. 936 The communication can be by a diffusive flux (diffusive BBL), an advective flux (advective BBL), or both. 937 In the current implementation of the BBL, only the tracers are modified, not the velocities. 938 Furthermore, it only connects ocean bottom cells, and therefore does not include all the improvements introduced by 939 \citet{Campin_Goosse_Tel99}. 946 940 947 941 % ------------------------------------------------------------------------------------------------------------- … … 951 945 \label{subsec:TRA_bbl_diff} 952 946 953 When applying sigma-diffusion (\key{trabbl} defined and \np{nn\_bbl\_ldf} set to 1), 947 When applying sigma-diffusion (\key{trabbl} defined and \np{nn\_bbl\_ldf} set to 1), 954 948 the diffusive flux between two adjacent cells at the ocean floor is given by 955 949 \begin{equation} \label{eq:tra_bbl_diff} 956 950 {\rm {\bf F}}_\sigma=A_l^\sigma \; \nabla_\sigma T 957 951 \end{equation} 958 with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, 959 and $A_l^\sigma$ the lateral diffusivity in the BBL. Following \citet{Beckmann_Doscher1997}, 960 the latter is prescribed with a spatial dependence, $i.e.$ in the conditional form 952 with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, 953 and $A_l^\sigma$ the lateral diffusivity in the BBL. 954 Following \citet{Beckmann_Doscher1997}, the latter is prescribed with a spatial dependence, 955 $i.e.$ in the conditional form 961 956 \begin{equation} \label{eq:tra_bbl_coef} 962 957 A_l^\sigma (i,j,t)=\left\{ {\begin{array}{l} … … 966 961 \end{array}} \right. 967 962 \end{equation} 968 where $A_{bbl}$ is the BBL diffusivity coefficient, given by the namelist 969 parameter \np{rn\_ahtbbl} and usually set to a value much larger 970 than the one used for lateral mixing in the open ocean. The constraint in \autoref{eq:tra_bbl_coef} 971 implies that sigma-like diffusion only occurs when the density above the sea floor, at the top of 972 the slope, is larger than in the deeper ocean (see green arrow in \autoref{fig:bbl}). 973 In practice, this constraint is applied separately in the two horizontal directions, 963 where $A_{bbl}$ is the BBL diffusivity coefficient, given by the namelist parameter \np{rn\_ahtbbl} and 964 usually set to a value much larger than the one used for lateral mixing in the open ocean. 965 The constraint in \autoref{eq:tra_bbl_coef} implies that sigma-like diffusion only occurs when 966 the density above the sea floor, at the top of the slope, is larger than in the deeper ocean 967 (see green arrow in \autoref{fig:bbl}). 968 In practice, this constraint is applied separately in the two horizontal directions, 974 969 and the density gradient in \autoref{eq:tra_bbl_coef} is evaluated with the log gradient formulation: 975 970 \begin{equation} \label{eq:tra_bbl_Drho} 976 971 \nabla_\sigma \rho / \rho = \alpha \,\nabla_\sigma T + \beta \,\nabla_\sigma S 977 972 \end{equation} 978 where $\rho$, $\alpha$ and $\beta$ are functions of $\overline{T}^\sigma$, 979 $\overline{S}^\sigma$ and $\overline{H}^\sigma$, the along bottom mean temperature, 980 salinity and depth, respectively. 973 where $\rho$, $\alpha$ and $\beta$ are functions of $\overline{T}^\sigma$, 974 $\overline{S}^\sigma$ and $\overline{H}^\sigma$, the along bottom mean temperature, salinity and depth, respectively. 981 975 982 976 % ------------------------------------------------------------------------------------------------------------- … … 990 984 991 985 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 992 \begin{figure}[!t] \begin{center} 993 \includegraphics[width=0.7\textwidth]{Fig_BBL_adv} 994 \caption{ \protect\label{fig:bbl} 995 Advective/diffusive Bottom Boundary Layer. The BBL parameterisation is 996 activated when $\rho^i_{kup}$ is larger than $\rho^{i+1}_{kdnw}$. 997 Red arrows indicate the additional overturning circulation due to the advective BBL. 998 The transport of the downslope flow is defined either as the transport of the bottom 999 ocean cell (black arrow), or as a function of the along slope density gradient. 1000 The green arrow indicates the diffusive BBL flux directly connecting $kup$ and $kdwn$ 1001 ocean bottom cells. 1002 connection} 1003 \end{center} \end{figure} 986 \begin{figure}[!t] 987 \begin{center} 988 \includegraphics[width=0.7\textwidth]{Fig_BBL_adv} 989 \caption{ \protect\label{fig:bbl} 990 Advective/diffusive Bottom Boundary Layer. 991 The BBL parameterisation is activated when $\rho^i_{kup}$ is larger than $\rho^{i+1}_{kdnw}$. 992 Red arrows indicate the additional overturning circulation due to the advective BBL. 993 The transport of the downslope flow is defined either as the transport of the bottom ocean cell (black arrow), 994 or as a function of the along slope density gradient. 995 The green arrow indicates the diffusive BBL flux directly connecting $kup$ and $kdwn$ ocean bottom cells. 996 } 997 \end{center} 998 \end{figure} 1004 999 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1005 1000 … … 1011 1006 %%%gmcomment : this section has to be really written 1012 1007 1013 When applying an advective BBL (\np{nn\_bbl\_adv}\forcode{ = 1..2}), an overturning 1014 c irculation is added which connects two adjacent bottom grid-points only if dense1015 water overlies less dense water on the slope. The density difference causes dense1016 water to move down the slope. 1017 1018 \np{nn\_bbl\_adv}\forcode{ = 1} : the downslope velocity is chosen to be the Eulerian 1019 ocean velocity just above the topographic step (see black arrow in \autoref{fig:bbl}) 1020 \citep{Beckmann_Doscher1997}. It is a \textit{conditional advection}, that is, advection 1021 i s allowed only if dense water overlies less dense water on the slope ($i.e.$1022 $\nabla_\sigma \rho \cdot \nabla H<0$) and if the velocity is directed towards 1023 greater depth ($i.e.$ $\vect{U} \cdot \nabla H>0$). 1024 1025 \np{nn\_bbl\_adv}\forcode{ = 2} :the downslope velocity is chosen to be proportional to $\Delta \rho$,1008 When applying an advective BBL (\np{nn\_bbl\_adv}\forcode{ = 1..2}), an overturning circulation is added which 1009 connects two adjacent bottom grid-points only if dense water overlies less dense water on the slope. 1010 The density difference causes dense water to move down the slope. 1011 1012 \np{nn\_bbl\_adv}\forcode{ = 1}: 1013 the downslope velocity is chosen to be the Eulerian ocean velocity just above the topographic step 1014 (see black arrow in \autoref{fig:bbl}) \citep{Beckmann_Doscher1997}. 1015 It is a \textit{conditional advection}, that is, advection is allowed only 1016 if dense water overlies less dense water on the slope ($i.e.$ $\nabla_\sigma \rho \cdot \nabla H<0$) and 1017 if the velocity is directed towards greater depth ($i.e.$ $\vect{U} \cdot \nabla H>0$). 1018 1019 \np{nn\_bbl\_adv}\forcode{ = 2}: 1020 the downslope velocity is chosen to be proportional to $\Delta \rho$, 1026 1021 the density difference between the higher cell and lower cell densities \citep{Campin_Goosse_Tel99}. 1027 The advection is allowed only if dense water overlies less dense water on the slope ($i.e.$1028 $\nabla_\sigma \rho \cdot \nabla H<0$). For example, the resulting transport of the 1029 downslope flow, here in the $i$-direction (\autoref{fig:bbl}), is simply given by the 1030 following expression:1022 The advection is allowed only if dense water overlies less dense water on the slope 1023 ($i.e.$ $\nabla_\sigma \rho \cdot \nabla H<0$). 1024 For example, the resulting transport of the downslope flow, here in the $i$-direction (\autoref{fig:bbl}), 1025 is simply given by the following expression: 1031 1026 \begin{equation} \label{eq:bbl_Utr} 1032 1027 u^{tr}_{bbl} = \gamma \, g \frac{\Delta \rho}{\rho_o} e_{1u} \; min \left( {e_{3u}}_{kup},{e_{3u}}_{kdwn} \right) 1033 1028 \end{equation} 1034 where $\gamma$, expressed in seconds, is the coefficient of proportionality 1035 provided as \np{rn\_gambbl}, a namelist parameter, and \textit{kup} and \textit{kdwn} 1036 are the vertical index of the higher and lower cells, respectively. 1037 The parameter $\gamma$ should take a different value for each bathymetric 1038 step, but for simplicity, and because no direct estimation of this parameter is 1039 available, a uniform value has been assumed. The possible values for $\gamma$ 1040 range between 1 and $10~s$ \citep{Campin_Goosse_Tel99}. 1041 1042 Scalar properties are advected by this additional transport $( u^{tr}_{bbl}, v^{tr}_{bbl} )$ 1043 using the upwind scheme. Such a diffusive advective scheme has been chosen 1044 to mimic the entrainment between the downslope plume and the surrounding 1045 water at intermediate depths. The entrainment is replaced by the vertical mixing 1046 implicit in the advection scheme. Let us consider as an example the 1047 case displayed in \autoref{fig:bbl} where the density at level $(i,kup)$ is 1048 larger than the one at level $(i,kdwn)$. The advective BBL scheme 1049 modifies the tracer time tendency of the ocean cells near the 1050 topographic step by the downslope flow \autoref{eq:bbl_dw}, 1051 the horizontal \autoref{eq:bbl_hor} and the upward \autoref{eq:bbl_up} 1052 return flows as follows: 1029 where $\gamma$, expressed in seconds, is the coefficient of proportionality provided as \np{rn\_gambbl}, 1030 a namelist parameter, and \textit{kup} and \textit{kdwn} are the vertical index of the higher and lower cells, 1031 respectively. 1032 The parameter $\gamma$ should take a different value for each bathymetric step, but for simplicity, 1033 and because no direct estimation of this parameter is available, a uniform value has been assumed. 1034 The possible values for $\gamma$ range between 1 and $10~s$ \citep{Campin_Goosse_Tel99}. 1035 1036 Scalar properties are advected by this additional transport $( u^{tr}_{bbl}, v^{tr}_{bbl} )$ using the upwind scheme. 1037 Such a diffusive advective scheme has been chosen to mimic the entrainment between the downslope plume and 1038 the surrounding water at intermediate depths. 1039 The entrainment is replaced by the vertical mixing implicit in the advection scheme. 1040 Let us consider as an example the case displayed in \autoref{fig:bbl} where 1041 the density at level $(i,kup)$ is larger than the one at level $(i,kdwn)$. 1042 The advective BBL scheme modifies the tracer time tendency of the ocean cells near the topographic step by 1043 the downslope flow \autoref{eq:bbl_dw}, the horizontal \autoref{eq:bbl_hor} and 1044 the upward \autoref{eq:bbl_up} return flows as follows: 1053 1045 \begin{align} 1054 1046 \partial_t T^{do}_{kdw} &\equiv \partial_t T^{do}_{kdw} … … 1065 1057 where $b_t$ is the $T$-cell volume. 1066 1058 1067 Note that the BBL transport, $( u^{tr}_{bbl}, v^{tr}_{bbl} )$, is available in 1068 the model outputs. It has to be used to compute the effective velocity 1069 as well as the effective overturning circulation. 1059 Note that the BBL transport, $( u^{tr}_{bbl}, v^{tr}_{bbl} )$, is available in the model outputs. 1060 It has to be used to compute the effective velocity as well as the effective overturning circulation. 1070 1061 1071 1062 % ================================================================ … … 1079 1070 %-------------------------------------------------------------------------------------------------------------- 1080 1071 1081 In some applications it can be useful to add a Newtonian damping term 1082 into the temperature and salinity equations: 1072 In some applications it can be useful to add a Newtonian damping term into the temperature and salinity equations: 1083 1073 \begin{equation} \label{eq:tra_dmp} 1084 1074 \begin{split} 1085 1075 \frac{\partial T}{\partial t}=\;\cdots \;-\gamma \,\left( {T-T_o } \right) \\ 1086 \frac{\partial S}{\partial t}=\;\cdots \;-\gamma \,\left( {S-S_o } \right) 1087 1088 1089 where $\gamma$ is the inverse of a time scale, and $T_o$ and $S_o$ 1090 are given temperature and salinity fields (usually a climatology). 1076 \frac{\partial S}{\partial t}=\;\cdots \;-\gamma \,\left( {S-S_o } \right) 1077 \end{split} 1078 \end{equation} 1079 where $\gamma$ is the inverse of a time scale, and $T_o$ and $S_o$ are given temperature and salinity fields 1080 (usually a climatology). 1091 1081 Options are defined through the \ngn{namtra\_dmp} namelist variables. 1092 The restoring term is added when the namelist parameter \np{ln\_tradmp} is set to true. 1093 It also requires that both \np{ln\_tsd\_init} and \np{ln\_tsd\_tradmp} are set to true 1094 in \textit{namtsd} namelist as well as \np{sn\_tem} and \np{sn\_sal} structures are 1095 correctly set ($i.e.$ that $T_o$ and $S_o$ are provided in input files and read 1096 using \mdl{fldread}, see \autoref{subsec:SBC_fldread}). 1097 The restoring coefficient $\gamma$ is a three-dimensional array read in during the \rou{tra\_dmp\_init} routine. The file name is specified by the namelist variable \np{cn\_resto}. The DMP\_TOOLS tool is provided to allow users to generate the netcdf file. 1098 1099 The two main cases in which \autoref{eq:tra_dmp} is used are \textit{(a)} 1100 the specification of the boundary conditions along artificial walls of a 1101 limited domain basin and \textit{(b)} the computation of the velocity 1102 field associated with a given $T$-$S$ field (for example to build the 1103 initial state of a prognostic simulation, or to use the resulting velocity 1104 field for a passive tracer study). The first case applies to regional 1105 models that have artificial walls instead of open boundaries. 1106 In the vicinity of these walls, $\gamma$ takes large values (equivalent to 1107 a time scale of a few days) whereas it is zero in the interior of the 1108 model domain. The second case corresponds to the use of the robust 1109 diagnostic method \citep{Sarmiento1982}. It allows us to find the velocity 1110 field consistent with the model dynamics whilst having a $T$, $S$ field 1111 close to a given climatological field ($T_o$, $S_o$). 1112 1113 The robust diagnostic method is very efficient in preventing temperature 1114 drift in intermediate waters but it produces artificial sources of heat and salt 1115 within the ocean. It also has undesirable effects on the ocean convection. 1116 It tends to prevent deep convection and subsequent deep-water formation, 1117 by stabilising the water column too much. 1118 1119 The namelist parameter \np{nn\_zdmp} sets whether the damping should be applied in the whole water column or only below the mixed layer (defined either on a density or $S_o$ criterion). It is common to set the damping to zero in the mixed layer as the adjustment time scale is short here \citep{Madec_al_JPO96}. 1082 The restoring term is added when the namelist parameter \np{ln\_tradmp} is set to true. 1083 It also requires that both \np{ln\_tsd\_init} and \np{ln\_tsd\_tradmp} are set to true in 1084 \textit{namtsd} namelist as well as \np{sn\_tem} and \np{sn\_sal} structures are correctly set 1085 ($i.e.$ that $T_o$ and $S_o$ are provided in input files and read using \mdl{fldread}, 1086 see \autoref{subsec:SBC_fldread}). 1087 The restoring coefficient $\gamma$ is a three-dimensional array read in during the \rou{tra\_dmp\_init} routine. 1088 The file name is specified by the namelist variable \np{cn\_resto}. 1089 The DMP\_TOOLS tool is provided to allow users to generate the netcdf file. 1090 1091 The two main cases in which \autoref{eq:tra_dmp} is used are 1092 \textit{(a)} the specification of the boundary conditions along artificial walls of a limited domain basin and 1093 \textit{(b)} the computation of the velocity field associated with a given $T$-$S$ field 1094 (for example to build the initial state of a prognostic simulation, 1095 or to use the resulting velocity field for a passive tracer study). 1096 The first case applies to regional models that have artificial walls instead of open boundaries. 1097 In the vicinity of these walls, $\gamma$ takes large values (equivalent to a time scale of a few days) whereas 1098 it is zero in the interior of the model domain. 1099 The second case corresponds to the use of the robust diagnostic method \citep{Sarmiento1982}. 1100 It allows us to find the velocity field consistent with the model dynamics whilst 1101 having a $T$, $S$ field close to a given climatological field ($T_o$, $S_o$). 1102 1103 The robust diagnostic method is very efficient in preventing temperature drift in intermediate waters but 1104 it produces artificial sources of heat and salt within the ocean. 1105 It also has undesirable effects on the ocean convection. 1106 It tends to prevent deep convection and subsequent deep-water formation, by stabilising the water column too much. 1107 1108 The namelist parameter \np{nn\_zdmp} sets whether the damping should be applied in the whole water column or 1109 only below the mixed layer (defined either on a density or $S_o$ criterion). 1110 It is common to set the damping to zero in the mixed layer as the adjustment time scale is short here 1111 \citep{Madec_al_JPO96}. 1120 1112 1121 1113 \subsection{Generating \ifile{resto} using DMP\_TOOLS} 1122 1114 1123 DMP\_TOOLS can be used to generate a netcdf file containing the restoration coefficient $\gamma$. 1124 Note that in order to maintain bit comparison with previous NEMO versions DMP\_TOOLS must be compiled 1125 and run on the same machine as the NEMO model. A \ifile{mesh\_mask} file for the model configuration is required as an input. 1126 This can be generated by carrying out a short model run with the namelist parameter \np{nn\_msh} set to 1. 1127 The namelist parameter \np{ln\_tradmp} will also need to be set to .false. for this to work. 1128 The \nl{nam\_dmp\_create} namelist in the DMP\_TOOLS directory is used to specify options for the restoration coefficient. 1115 DMP\_TOOLS can be used to generate a netcdf file containing the restoration coefficient $\gamma$. 1116 Note that in order to maintain bit comparison with previous NEMO versions DMP\_TOOLS must be compiled and 1117 run on the same machine as the NEMO model. 1118 A \ifile{mesh\_mask} file for the model configuration is required as an input. 1119 This can be generated by carrying out a short model run with the namelist parameter \np{nn\_msh} set to 1. 1120 The namelist parameter \np{ln\_tradmp} will also need to be set to .false. for this to work. 1121 The \nl{nam\_dmp\_create} namelist in the DMP\_TOOLS directory is used to specify options for 1122 the restoration coefficient. 1129 1123 1130 1124 %--------------------------------------------nam_dmp_create------------------------------------------------- … … 1132 1126 %------------------------------------------------------------------------------------------------------- 1133 1127 1134 \np{cp\_cfg}, \np{cp\_cpz}, \np{jp\_cfg} and \np{jperio} specify the model configuration being used and should be the same as specified in \nl{namcfg}. The variable \nl{lzoom} is used to specify that the damping is being used as in case \textit{a} above to provide boundary conditions to a zoom configuration. In the case of the arctic or antarctic zoom configurations this includes some specific treatment. Otherwise damping is applied to the 6 grid points along the ocean boundaries. The open boundaries are specified by the variables \np{lzoom\_n}, \np{lzoom\_e}, \np{lzoom\_s}, \np{lzoom\_w} in the \nl{nam\_zoom\_dmp} name list. 1135 1136 The remaining switch namelist variables determine the spatial variation of the restoration coefficient in non-zoom configurations. 1137 \np{ln\_full\_field} specifies that newtonian damping should be applied to the whole model domain. 1138 \np{ln\_med\_red\_seas} specifies grid specific restoration coefficients in the Mediterranean Sea 1139 for the ORCA4, ORCA2 and ORCA05 configurations. 1140 If \np{ln\_old\_31\_lev\_code} is set then the depth variation of the coeffients will be specified as 1141 a function of the model number. This option is included to allow backwards compatability of the ORCA2 reference 1142 configurations with previous model versions. 1143 \np{ln\_coast} specifies that the restoration coefficient should be reduced near to coastlines. 1144 This option only has an effect if \np{ln\_full\_field} is true. 1145 \np{ln\_zero\_top\_layer} specifies that the restoration coefficient should be zero in the surface layer. 1146 Finally \np{ln\_custom} specifies that the custom module will be called. 1147 This module is contained in the file \mdl{custom} and can be edited by users. For example damping could be applied in a specific region. 1148 1149 The restoration coefficient can be set to zero in equatorial regions by specifying a positive value of \np{nn\_hdmp}. 1128 \np{cp\_cfg}, \np{cp\_cpz}, \np{jp\_cfg} and \np{jperio} specify the model configuration being used and 1129 should be the same as specified in \nl{namcfg}. 1130 The variable \nl{lzoom} is used to specify that the damping is being used as in case \textit{a} above to 1131 provide boundary conditions to a zoom configuration. 1132 In the case of the arctic or antarctic zoom configurations this includes some specific treatment. 1133 Otherwise damping is applied to the 6 grid points along the ocean boundaries. 1134 The open boundaries are specified by the variables \np{lzoom\_n}, \np{lzoom\_e}, \np{lzoom\_s}, \np{lzoom\_w} in 1135 the \nl{nam\_zoom\_dmp} name list. 1136 1137 The remaining switch namelist variables determine the spatial variation of the restoration coefficient in 1138 non-zoom configurations. 1139 \np{ln\_full\_field} specifies that newtonian damping should be applied to the whole model domain. 1140 \np{ln\_med\_red\_seas} specifies grid specific restoration coefficients in the Mediterranean Sea for 1141 the ORCA4, ORCA2 and ORCA05 configurations. 1142 If \np{ln\_old\_31\_lev\_code} is set then the depth variation of the coeffients will be specified as 1143 a function of the model number. 1144 This option is included to allow backwards compatability of the ORCA2 reference configurations with 1145 previous model versions. 1146 \np{ln\_coast} specifies that the restoration coefficient should be reduced near to coastlines. 1147 This option only has an effect if \np{ln\_full\_field} is true. 1148 \np{ln\_zero\_top\_layer} specifies that the restoration coefficient should be zero in the surface layer. 1149 Finally \np{ln\_custom} specifies that the custom module will be called. 1150 This module is contained in the file \mdl{custom} and can be edited by users. 1151 For example damping could be applied in a specific region. 1152 1153 The restoration coefficient can be set to zero in equatorial regions by 1154 specifying a positive value of \np{nn\_hdmp}. 1150 1155 Equatorward of this latitude the restoration coefficient will be zero with a smooth transition to 1151 1156 the full values of a 10\deg latitud band. 1152 This is often used because of the short adjustment time scale in the equatorial region 1153 \citep{Reverdin1991, Fujio1991, Marti_PhD92}. The time scale associated with the damping depends on the depth as a 1154 hyperbolic tangent, with \np{rn\_surf} as surface value, \np{rn\_bot} as bottom value and a transition depth of \np{rn\_dep}. 1157 This is often used because of the short adjustment time scale in the equatorial region 1158 \citep{Reverdin1991, Fujio1991, Marti_PhD92}. 1159 The time scale associated with the damping depends on the depth as a hyperbolic tangent, 1160 with \np{rn\_surf} as surface value, \np{rn\_bot} as bottom value and a transition depth of \np{rn\_dep}. 1155 1161 1156 1162 % ================================================================ … … 1165 1171 1166 1172 Options are defined through the \ngn{namdom} namelist variables. 1167 The general framework for tracer time stepping is a modified leap-frog scheme 1168 \citep{Leclair_Madec_OM09}, $i.e.$ a three level centred time scheme associated 1169 with a Asselin time filter (cf. \autoref{sec:STP_mLF}): 1173 The general framework for tracer time stepping is a modified leap-frog scheme \citep{Leclair_Madec_OM09}, 1174 $i.e.$ a three level centred time scheme associated with a Asselin time filter (cf. \autoref{sec:STP_mLF}): 1170 1175 \begin{equation} \label{eq:tra_nxt} 1171 1176 \begin{aligned} … … 1177 1182 \end{aligned} 1178 1183 \end{equation} 1179 where RHS is the right hand side of the temperature equation, 1180 the subscript $f$ denotes filtered values, $\gamma$ is the Asselin coefficient, 1181 and $S$ is the total forcing applied on $T$ ($i.e.$ fluxes plus content in mass exchanges). 1182 $\gamma$ is initialized as \np{rn\_atfp} (\textbf{namelist} parameter). 1183 Its default value is \np{rn\_atfp}\forcode{ = 10.e-3}. Note that the forcing correction term in the filter1184 is not applied in linear free surface (\jp{lk\_vvl}\forcode{ = .false.}) (see \autoref{subsec:TRA_sbc}. 1185 Not also that in constant volume case, the time stepping is performed on $T$, 1186 not on its content, $e_{3t}T$.1187 1188 When the vertical mixing is solved implicitly, the update of the \textit{next} tracer1189 fields is done in module \mdl{trazdf}. In this case only the swapping of arrays 1190 and the Asselin filtering is done in the \mdl{tranxt} module.1191 1192 In order to prepare for the computation of the \textit{next} time step, 1193 a swap of tracer arrays is performed:$T^{t-\rdt} = T^t$ and $T^t = T_f$.1184 where RHS is the right hand side of the temperature equation, the subscript $f$ denotes filtered values, 1185 $\gamma$ is the Asselin coefficient, and $S$ is the total forcing applied on $T$ 1186 ($i.e.$ fluxes plus content in mass exchanges). 1187 $\gamma$ is initialized as \np{rn\_atfp} (\textbf{namelist} parameter). 1188 Its default value is \np{rn\_atfp}\forcode{ = 10.e-3}. 1189 Note that the forcing correction term in the filter is not applied in linear free surface 1190 (\jp{lk\_vvl}\forcode{ = .false.}) (see \autoref{subsec:TRA_sbc}. 1191 Not also that in constant volume case, the time stepping is performed on $T$, not on its content, $e_{3t}T$. 1192 1193 When the vertical mixing is solved implicitly, 1194 the update of the \textit{next} tracer fields is done in module \mdl{trazdf}. 1195 In this case only the swapping of arrays and the Asselin filtering is done in the \mdl{tranxt} module. 1196 1197 In order to prepare for the computation of the \textit{next} time step, a swap of tracer arrays is performed: 1198 $T^{t-\rdt} = T^t$ and $T^t = T_f$. 1194 1199 1195 1200 % ================================================================ … … 1209 1214 \label{subsec:TRA_eos} 1210 1215 1211 The Equation Of Seawater (EOS) is an empirical nonlinear thermodynamic relationship 1212 linking seawater density, $\rho$, to a number of state variables, 1213 most typically temperature, salinity and pressure. 1214 Because density gradients control the pressure gradient force through the hydrostatic balance, 1215 the equation of state provides a fundamental bridge between the distribution of active tracers1216 and the fluid dynamics. Nonlinearities of the EOS are of major importance, in particular 1217 influencing the circulation through determination of the static stability below the mixed layer, 1218 thus controlling rates of exchange between the atmosphere and the ocean interior \citep{Roquet_JPO2015}.1219 Therefore an accurate EOS based on either the 1980 equation of state (EOS-80, \cite{UNESCO1983}) 1220 or TEOS-10 \citep{TEOS10} standards should be used anytime a simulation of the real 1221 ocean circulation is attempted \citep{Roquet_JPO2015}. 1222 The use of TEOS-10 is highly recommended because 1223 \textit{(i)} it is the new official EOS,1224 \textit{(ii)} it is more accurate, being based on an updated database of laboratory measurements, and1225 \textit{(iii)} it uses Conservative Temperature and Absolute Salinity (instead of potential temperature 1226 and practical salinity for EOS-980, both variables being more suitable for use as model variables 1227 \citep{TEOS10, Graham_McDougall_JPO13}. 1216 The Equation Of Seawater (EOS) is an empirical nonlinear thermodynamic relationship linking seawater density, 1217 $\rho$, to a number of state variables, most typically temperature, salinity and pressure. 1218 Because density gradients control the pressure gradient force through the hydrostatic balance, 1219 the equation of state provides a fundamental bridge between the distribution of active tracers and 1220 the fluid dynamics. 1221 Nonlinearities of the EOS are of major importance, in particular influencing the circulation through 1222 determination of the static stability below the mixed layer, 1223 thus controlling rates of exchange between the atmosphere and the ocean interior \citep{Roquet_JPO2015}. 1224 Therefore an accurate EOS based on either the 1980 equation of state (EOS-80, \cite{UNESCO1983}) or 1225 TEOS-10 \citep{TEOS10} standards should be used anytime a simulation of the real ocean circulation is attempted 1226 \citep{Roquet_JPO2015}. 1227 The use of TEOS-10 is highly recommended because 1228 \textit{(i)} it is the new official EOS, 1229 \textit{(ii)} it is more accurate, being based on an updated database of laboratory measurements, and 1230 \textit{(iii)} it uses Conservative Temperature and Absolute Salinity (instead of potential temperature and 1231 practical salinity for EOS-980, both variables being more suitable for use as model variables 1232 \citep{TEOS10, Graham_McDougall_JPO13}. 1228 1233 EOS-80 is an obsolescent feature of the NEMO system, kept only for backward compatibility. 1229 For process studies, it is often convenient to use an approximation of the EOS. To that purposed, 1230 a simplified EOS (S-EOS) inspired by \citet{Vallis06} is also available. 1231 1232 In the computer code, a density anomaly, $d_a= \rho / \rho_o - 1$, 1233 is computed, with $\rho_o$ a reference density. Called \textit{rau0} 1234 in the code, $\rho_o$ is set in \mdl{phycst} to a value of $1,026~Kg/m^3$. 1235 This is a sensible choice for the reference density used in a Boussinesq ocean 1236 climate model, as, with the exception of only a small percentage of the ocean, 1234 For process studies, it is often convenient to use an approximation of the EOS. 1235 To that purposed, a simplified EOS (S-EOS) inspired by \citet{Vallis06} is also available. 1236 1237 In the computer code, a density anomaly, $d_a= \rho / \rho_o - 1$, is computed, with $\rho_o$ a reference density. 1238 Called \textit{rau0} in the code, $\rho_o$ is set in \mdl{phycst} to a value of $1,026~Kg/m^3$. 1239 This is a sensible choice for the reference density used in a Boussinesq ocean climate model, as, 1240 with the exception of only a small percentage of the ocean, 1237 1241 density in the World Ocean varies by no more than 2$\%$ from that value \citep{Gill1982}. 1238 1242 1239 Options are defined through the \ngn{nameos} namelist variables, and in particular \np{nn\_eos} 1240 whichcontrols the EOS used (\forcode{= -1} for TEOS10 ; \forcode{= 0} for EOS-80 ; \forcode{= 1} for S-EOS).1243 Options are defined through the \ngn{nameos} namelist variables, and in particular \np{nn\_eos} which 1244 controls the EOS used (\forcode{= -1} for TEOS10 ; \forcode{= 0} for EOS-80 ; \forcode{= 1} for S-EOS). 1241 1245 \begin{description} 1242 1243 \item[\np{nn\_eos}\forcode{ = -1}] the polyTEOS10-bsq equation of seawater \citep{Roquet_OM2015} is used. 1244 The accuracy of this approximation is comparable to the TEOS-10 rational function approximation, 1245 but it is optimized for a boussinesq fluid and the polynomial expressions have simpler 1246 and more computationally efficient expressions for their derived quantities 1247 which make them more adapted for use in ocean models. 1248 Note that a slightly higher precision polynomial form is now used replacement of the TEOS-10 1249 rational function approximation for hydrographic data analysis \citep{TEOS10}. 1250 A key point is that conservative state variables are used: 1251 Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: \degC, notation: $\Theta$). 1252 The pressure in decibars is approximated by the depth in meters. 1253 With TEOS10, the specific heat capacity of sea water, $C_p$, is a constant. It is set to 1254 $C_p=3991.86795711963~J\,Kg^{-1}\,^{\circ}K^{-1}$, according to \citet{TEOS10}. 1255 1256 Choosing polyTEOS10-bsq implies that the state variables used by the model are 1257 $\Theta$ and $S_A$. In particular, the initial state deined by the user have to be given as 1258 \textit{Conservative} Temperature and \textit{Absolute} Salinity. 1259 In addition, setting \np{ln\_useCT} to \forcode{.true.} convert the Conservative SST to potential SST 1260 prior to either computing the air-sea and ice-sea fluxes (forced mode) 1261 or sending the SST field to the atmosphere (coupled mode). 1262 1263 \item[\np{nn\_eos}\forcode{ = 0}] the polyEOS80-bsq equation of seawater is used. 1264 It takes the same polynomial form as the polyTEOS10, but the coefficients have been optimized 1265 to accurately fit EOS80 (Roquet, personal comm.). The state variables used in both the EOS80 1266 and the ocean model are: 1267 the Practical Salinity ((unit: psu, notation: $S_p$)) and Potential Temperature (unit: $^{\circ}C$, notation: $\theta$). 1268 The pressure in decibars is approximated by the depth in meters. 1269 With thsi EOS, the specific heat capacity of sea water, $C_p$, is a function of temperature, 1270 salinity and pressure \citep{UNESCO1983}. Nevertheless, a severe assumption is made in order to 1271 have a heat content ($C_p T_p$) which is conserved by the model: $C_p$ is set to a constant 1272 value, the TEOS10 value. 1246 \item[\np{nn\_eos}\forcode{ = -1}] 1247 the polyTEOS10-bsq equation of seawater \citep{Roquet_OM2015} is used. 1248 The accuracy of this approximation is comparable to the TEOS-10 rational function approximation, 1249 but it is optimized for a boussinesq fluid and the polynomial expressions have simpler and 1250 more computationally efficient expressions for their derived quantities which make them more adapted for 1251 use in ocean models. 1252 Note that a slightly higher precision polynomial form is now used replacement of 1253 the TEOS-10 rational function approximation for hydrographic data analysis \citep{TEOS10}. 1254 A key point is that conservative state variables are used: 1255 Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: \degC, notation: $\Theta$). 1256 The pressure in decibars is approximated by the depth in meters. 1257 With TEOS10, the specific heat capacity of sea water, $C_p$, is a constant. 1258 It is set to $C_p=3991.86795711963~J\,Kg^{-1}\,^{\circ}K^{-1}$, according to \citet{TEOS10}. 1259 1260 Choosing polyTEOS10-bsq implies that the state variables used by the model are $\Theta$ and $S_A$. 1261 In particular, the initial state deined by the user have to be given as \textit{Conservative} Temperature and 1262 \textit{Absolute} Salinity. 1263 In addition, setting \np{ln\_useCT} to \forcode{.true.} convert the Conservative SST to potential SST prior to 1264 either computing the air-sea and ice-sea fluxes (forced mode) or 1265 sending the SST field to the atmosphere (coupled mode). 1266 1267 \item[\np{nn\_eos}\forcode{ = 0}] 1268 the polyEOS80-bsq equation of seawater is used. 1269 It takes the same polynomial form as the polyTEOS10, but the coefficients have been optimized to 1270 accurately fit EOS80 (Roquet, personal comm.). 1271 The state variables used in both the EOS80 and the ocean model are: 1272 the Practical Salinity ((unit: psu, notation: $S_p$)) and 1273 Potential Temperature (unit: $^{\circ}C$, notation: $\theta$). 1274 The pressure in decibars is approximated by the depth in meters. 1275 With thsi EOS, the specific heat capacity of sea water, $C_p$, is a function of temperature, salinity and 1276 pressure \citep{UNESCO1983}. 1277 Nevertheless, a severe assumption is made in order to have a heat content ($C_p T_p$) which 1278 is conserved by the model: $C_p$ is set to a constant value, the TEOS10 value. 1273 1279 1274 \item[\np{nn\_eos}\forcode{ = 1}] a simplified EOS (S-EOS) inspired by \citet{Vallis06} is chosen, 1275 the coefficients of which has been optimized to fit the behavior of TEOS10 (Roquet, personal comm.) 1276 (see also \citet{Roquet_JPO2015}). It provides a simplistic linear representation of both 1277 cabbeling and thermobaricity effects which is enough for a proper treatment of the EOS 1278 in theoretical studies \citep{Roquet_JPO2015}. 1279 With such an equation of state there is no longer a distinction between 1280 \textit{conservative} and \textit{potential} temperature, as well as between \textit{absolute} 1281 and \textit{practical} salinity. 1282 S-EOS takes the following expression: 1283 \begin{equation} \label{eq:tra_S-EOS} 1284 \begin{split} 1285 d_a(T,S,z) = ( & - a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * T_a \\ 1286 & + b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * S_a \\ 1287 & - \nu \; T_a \; S_a \; ) \; / \; \rho_o \\ 1288 with \ \ T_a = T-10 \; ; & \; S_a = S-35 \; ;\; \rho_o = 1026~Kg/m^3 1289 \end{split} 1290 \end{equation} 1291 where the computer name of the coefficients as well as their standard value are given in \autoref{tab:SEOS}. 1292 In fact, when choosing S-EOS, various approximation of EOS can be specified simply by changing 1293 the associated coefficients. 1294 Setting to zero the two thermobaric coefficients ($\mu_1$, $\mu_2$) remove thermobaric effect from S-EOS. 1295 setting to zero the three cabbeling coefficients ($\lambda_1$, $\lambda_2$, $\nu$) remove cabbeling effect from S-EOS. 1296 Keeping non-zero value to $a_0$ and $b_0$ provide a linear EOS function of T and S. 1297 1280 \item[\np{nn\_eos}\forcode{ = 1}] 1281 a simplified EOS (S-EOS) inspired by \citet{Vallis06} is chosen, 1282 the coefficients of which has been optimized to fit the behavior of TEOS10 1283 (Roquet, personal comm.) (see also \citet{Roquet_JPO2015}). 1284 It provides a simplistic linear representation of both cabbeling and thermobaricity effects which 1285 is enough for a proper treatment of the EOS in theoretical studies \citep{Roquet_JPO2015}. 1286 With such an equation of state there is no longer a distinction between 1287 \textit{conservative} and \textit{potential} temperature, 1288 as well as between \textit{absolute} and \textit{practical} salinity. 1289 S-EOS takes the following expression: 1290 \begin{equation} \label{eq:tra_S-EOS} 1291 \begin{split} 1292 d_a(T,S,z) = ( & - a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * T_a \\ 1293 & + b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * S_a \\ 1294 & - \nu \; T_a \; S_a \; ) \; / \; \rho_o \\ 1295 with \ \ T_a = T-10 \; ; & \; S_a = S-35 \; ;\; \rho_o = 1026~Kg/m^3 1296 \end{split} 1297 \end{equation} 1298 where the computer name of the coefficients as well as their standard value are given in \autoref{tab:SEOS}. 1299 In fact, when choosing S-EOS, various approximation of EOS can be specified simply by changing the associated coefficients. 1300 Setting to zero the two thermobaric coefficients ($\mu_1$, $\mu_2$) remove thermobaric effect from S-EOS. 1301 setting to zero the three cabbeling coefficients ($\lambda_1$, $\lambda_2$, $\nu$) remove cabbeling effect from S-EOS. 1302 Keeping non-zero value to $a_0$ and $b_0$ provide a linear EOS function of T and S. 1298 1303 \end{description} 1299 1304 … … 1313 1318 \end{tabular} 1314 1319 \caption{ \protect\label{tab:SEOS} 1315 Standard value of S-EOS coefficients. } 1320 Standard value of S-EOS coefficients. 1321 } 1316 1322 \end{center} 1317 1323 \end{table} … … 1325 1331 \label{subsec:TRA_bn2} 1326 1332 1327 An accurate computation of the ocean stability (i.e. of $N$, the brunt-V\"{a}is\"{a}l\"{a} 1328 frequency) is of paramount importance as determine the ocean stratification and 1329 is used in several ocean parameterisations (namely TKE, GLS, Richardson number dependent 1330 vertical diffusion, enhanced vertical diffusion, non-penetrative convection, tidal mixing 1331 parameterisation, iso-neutral diffusion). In particular, $N^2$ has to be computed at the local pressure 1332 (pressure in decibar being approximated by the depth in meters). The expression for $N^2$ 1333 is given by:1333 An accurate computation of the ocean stability (i.e. of $N$, the brunt-V\"{a}is\"{a}l\"{a} frequency) is of 1334 paramount importance as determine the ocean stratification and is used in several ocean parameterisations 1335 (namely TKE, GLS, Richardson number dependent vertical diffusion, enhanced vertical diffusion, 1336 non-penetrative convection, tidal mixing parameterisation, iso-neutral diffusion). 1337 In particular, $N^2$ has to be computed at the local pressure 1338 (pressure in decibar being approximated by the depth in meters). 1339 The expression for $N^2$ is given by: 1334 1340 \begin{equation} \label{eq:tra_bn2} 1335 1341 N^2 = \frac{g}{e_{3w}} \left( \beta \;\delta_{k+1/2}[S] - \alpha \;\delta_{k+1/2}[T] \right) 1336 1342 \end{equation} 1337 where $(T,S) = (\Theta, S_A)$ for TEOS10, $= (\theta, S_p)$ for TEOS-80, or $=(T,S)$ for S-EOS, 1338 and, $\alpha$ and $\beta$ are the thermal and haline expansion coefficients. 1339 The coefficients are a polynomial function of temperature, salinity and depth which expression1340 depends on the chosen EOS. They are computed through \textit{eos\_rab}, a \textsc{Fortran} 1341 function that can be found in \mdl{eosbn2}.1343 where $(T,S) = (\Theta, S_A)$ for TEOS10, $= (\theta, S_p)$ for TEOS-80, or $=(T,S)$ for S-EOS, 1344 and, $\alpha$ and $\beta$ are the thermal and haline expansion coefficients. 1345 The coefficients are a polynomial function of temperature, salinity and depth which 1346 expression depends on the chosen EOS. 1347 They are computed through \textit{eos\_rab}, a \textsc{Fortran} function that can be found in \mdl{eosbn2}. 1342 1348 1343 1349 % ------------------------------------------------------------------------------------------------------------- … … 1356 1362 \end{equation} 1357 1363 1358 \autoref{eq:tra_eos_fzp} is only used to compute the potential freezing point of 1359 sea water ($i.e.$ referenced to the surface $p=0$), thus the pressure dependent 1360 t erms in \autoref{eq:tra_eos_fzp} (last term) have been dropped. The freezing1361 point is computed through \textit{eos\_fzp}, a \textsc{Fortran} function that can be found 1362 in \mdl{eosbn2}.1364 \autoref{eq:tra_eos_fzp} is only used to compute the potential freezing point of sea water 1365 ($i.e.$ referenced to the surface $p=0$), 1366 thus the pressure dependent terms in \autoref{eq:tra_eos_fzp} (last term) have been dropped. 1367 The freezing point is computed through \textit{eos\_fzp}, 1368 a \textsc{Fortran} function that can be found in \mdl{eosbn2}. 1363 1369 1364 1370 … … 1380 1386 1381 1387 \gmcomment{STEVEN: to be consistent with earlier discussion of differencing and averaging operators, 1382 I've changed "derivative" to "difference" and "mean" to "average"} 1383 1384 With partial cells (\np{ln\_zps}\forcode{ = .true.}) at bottom and top (\np{ln\_isfcav}\forcode{ = .true.}), in general, 1385 tracers in horizontally adjacent cells live at different depths. 1386 Horizontal gradients of tracers are needed for horizontal diffusion (\mdl{traldf} module) 1387 and the hydrostatic pressure gradient calculations (\mdl{dynhpg} module). 1388 The partial cell properties at the top (\np{ln\_isfcav}\forcode{ = .true.}) are computed in the same way as for the bottom. 1388 I've changed "derivative" to "difference" and "mean" to "average"} 1389 1390 With partial cells (\np{ln\_zps}\forcode{ = .true.}) at bottom and top (\np{ln\_isfcav}\forcode{ = .true.}), 1391 in general, tracers in horizontally adjacent cells live at different depths. 1392 Horizontal gradients of tracers are needed for horizontal diffusion (\mdl{traldf} module) and 1393 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 as 1395 for the bottom. 1389 1396 So, only the bottom interpolation is explained below. 1390 1397 1391 Before taking horizontal gradients between the tracers next to the bottom, a linear 1392 interpolation in the vertical is used to approximate the deeper tracer as if it actually 1393 lived at the depth of the shallower tracer point (\autoref{fig:Partial_step_scheme}). 1394 For example, for temperature in the $i$-direction the needed interpolated 1395 temperature, $\widetilde{T}$, is: 1398 Before taking horizontal gradients between the tracers next to the bottom, 1399 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}). 1401 For example, for temperature in the $i$-direction the needed interpolated temperature, $\widetilde{T}$, is: 1396 1402 1397 1403 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1398 \begin{figure}[!p] \begin{center} 1399 \includegraphics[width=0.9\textwidth]{Partial_step_scheme} 1400 \caption{ \protect\label{fig:Partial_step_scheme} 1401 Discretisation of the horizontal difference and average of tracers in the $z$-partial 1402 step coordinate (\protect\np{ln\_zps}\forcode{ = .true.}) in the case $( e3w_k^{i+1} - e3w_k^i )>0$. 1403 A linear interpolation is used to estimate $\widetilde{T}_k^{i+1}$, the tracer value 1404 at the depth of the shallower tracer point of the two adjacent bottom $T$-points. 1405 The horizontal difference is then given by: $\delta _{i+1/2} T_k= \widetilde{T}_k^{\,i+1} -T_k^{\,i}$ 1406 and the average by: $\overline{T}_k^{\,i+1/2}= ( \widetilde{T}_k^{\,i+1/2} - T_k^{\,i} ) / 2$. } 1407 \end{center} \end{figure} 1404 \begin{figure}[!p] 1405 \begin{center} 1406 \includegraphics[width=0.9\textwidth]{Fig_partial_step_scheme} 1407 \caption{ \protect\label{fig:Partial_step_scheme} 1408 Discretisation of the horizontal difference and average of tracers in the $z$-partial step coordinate 1409 (\protect\np{ln\_zps}\forcode{ = .true.}) in the case $( e3w_k^{i+1} - e3w_k^i )>0$. 1410 A linear interpolation is used to estimate $\widetilde{T}_k^{i+1}$, 1411 the tracer value at the depth of the shallower tracer point of the two adjacent bottom $T$-points. 1412 The horizontal difference is then given by: $\delta _{i+1/2} T_k= \widetilde{T}_k^{\,i+1} -T_k^{\,i}$ and 1413 the average by: $\overline{T}_k^{\,i+1/2}= ( \widetilde{T}_k^{\,i+1/2} - T_k^{\,i} ) / 2$. 1414 } 1415 \end{center} 1416 \end{figure} 1408 1417 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1409 1418 \begin{equation*} … … 1416 1425 \end{aligned} \right. 1417 1426 \end{equation*} 1418 and the resulting forms for the horizontal difference and the horizontal average 1419 value of $T$ at a $U$-point are: 1427 and the resulting forms for the horizontal difference and the horizontal average value of $T$ at a $U$-point are: 1420 1428 \begin{equation} \label{eq:zps_hde} 1421 1429 \begin{aligned} … … 1434 1442 \end{equation} 1435 1443 1436 The computation of horizontal derivative of tracers as well as of density is 1437 performed once for all at each time step in \mdl{zpshde} module and stored 1438 in shared arrays to be used when needed. It has to be emphasized that the 1439 procedure used to compute the interpolated density, $\widetilde{\rho}$, is not 1440 the same as that used for $T$ and $S$. Instead of forming a linear approximation 1441 of density, we compute $\widetilde{\rho }$ from the interpolated values of $T$ 1442 and $S$, and the pressure at a $u$-point (in the equation of state pressure is 1443 approximated by depth, see \autoref{subsec:TRA_eos} ) : 1444 The computation of horizontal derivative of tracers as well as of density is performed once for all at 1445 each time step in \mdl{zpshde} module and stored in shared arrays to be used when needed. 1446 It has to be emphasized that the procedure used to compute the interpolated density, $\widetilde{\rho}$, 1447 is not the same as that used for $T$ and $S$. 1448 Instead of forming a linear approximation of density, we compute $\widetilde{\rho }$ from the interpolated values of 1449 $T$ and $S$, and the pressure at a $u$-point 1450 (in the equation of state pressure is approximated by depth, see \autoref{subsec:TRA_eos} ): 1444 1451 \begin{equation} \label{eq:zps_hde_rho} 1445 1452 \widetilde{\rho } = \rho ( {\widetilde{T},\widetilde {S},z_u }) … … 1447 1454 \end{equation} 1448 1455 1449 This is a much better approximation as the variation of $\rho$ with depth (and 1450 thus pressure) is highly non-linear with a true equation of state and thus is badly 1451 approximated with a linear interpolation. This approximation is used to compute 1452 both the horizontal pressure gradient (\autoref{sec:DYN_hpg}) and the slopes of neutral 1453 surfaces (\autoref{sec:LDF_slp}) 1454 1455 Note that in almost all the advection schemes presented in this Chapter, both 1456 averaging and differencing operators appear. Yet \autoref{eq:zps_hde} has not 1457 been used in these schemes: in contrast to diffusion and pressure gradient 1458 computations, no correction for partial steps is applied for advection. The main 1459 motivation is to preserve the domain averaged mean variance of the advected 1460 field when using the $2^{nd}$ order centred scheme. Sensitivity of the advection 1461 schemes to the way horizontal averages are performed in the vicinity of partial 1462 cells should be further investigated in the near future.1456 This is a much better approximation as the variation of $\rho$ with depth (and thus pressure) 1457 is highly non-linear with a true equation of state and thus is badly approximated with a linear interpolation. 1458 This approximation is used to compute both the horizontal pressure gradient (\autoref{sec:DYN_hpg}) and 1459 the slopes of neutral surfaces (\autoref{sec:LDF_slp}). 1460 1461 Note that in almost all the advection schemes presented in this Chapter, 1462 both averaging and differencing operators appear. 1463 Yet \autoref{eq:zps_hde} has not been used in these schemes: 1464 in contrast to diffusion and pressure gradient computations, 1465 no correction for partial steps is applied for advection. 1466 The main motivation is to preserve the domain averaged mean variance of the advected field when 1467 using the $2^{nd}$ order centred scheme. 1468 Sensitivity of the advection schemes to the way horizontal averages are performed in the vicinity of 1469 partial cells should be further investigated in the near future. 1463 1470 %%% 1464 1471 \gmcomment{gm : this last remark has to be done}
Note: See TracChangeset
for help on using the changeset viewer.