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