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