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