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