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