[707] | 1 | % ================================================================ |
---|
| 2 | % Chapter 1 Ñ Ocean Tracers (TRA) |
---|
| 3 | % ================================================================ |
---|
| 4 | \chapter{Ocean Tracers (TRA)} |
---|
| 5 | \label{TRA} |
---|
| 6 | \minitoc |
---|
| 7 | |
---|
| 8 | % missing/update |
---|
| 9 | % traqsr: need to coordinate with SBC module |
---|
| 10 | % trabbl : advective case to be discussed |
---|
| 11 | % diffusive case : add : only the bottom ocean cell is concerned |
---|
| 12 | % ==> addfigure on bbl |
---|
| 13 | |
---|
[817] | 14 | %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] | 15 | |
---|
[1224] | 16 | %\newpage |
---|
| 17 | \vspace{2.cm} |
---|
| 18 | %$\ $\newline % force a new ligne |
---|
[707] | 19 | |
---|
[817] | 20 | Using the representation described in Chap.~\ref{DOM}, several semi-discrete |
---|
| 21 | space forms of the tracer equations are available depending on the vertical |
---|
[707] | 22 | coordinate used and on the physics used. In all the equations presented |
---|
| 23 | here, the masking has been omitted for simplicity. One must be aware that |
---|
| 24 | all the quantities are masked fields and that each time a mean or difference |
---|
| 25 | operator is used, the resulting field is multiplied by a mask. |
---|
| 26 | |
---|
[817] | 27 | The two active tracers are potential temperature and salinity. Their prognostic |
---|
| 28 | equations can be summarized as follows: |
---|
[707] | 29 | \begin{equation*} |
---|
| 30 | \text{NXT} = \text{ADV}+\text{LDF}+\text{ZDF}+\text{SBC} |
---|
| 31 | \ (+\text{QSR})\ (+\text{BBC})\ (+\text{BBL})\ (+\text{DMP}) |
---|
| 32 | \end{equation*} |
---|
| 33 | |
---|
[817] | 34 | NXT stands for next, referring to the time-stepping. From left to right, the terms |
---|
| 35 | on the rhs of the tracer equations are the advection (ADV), the lateral diffusion |
---|
| 36 | (LDF), the vertical diffusion (ZDF), the contributions from the external forcings |
---|
| 37 | (SBC: Surface Boundary Condition, QSR: penetrative Solar Radiation, and BBC: |
---|
| 38 | Bottom Boundary Condition), the contribution from the bottom boundary Layer |
---|
| 39 | (BBL) parametrisation, and an internal damping (DMP) term. The terms QSR, |
---|
[1224] | 40 | BBC, BBL and DMP are optional. The external forcings and parameterisations |
---|
[817] | 41 | require complex inputs and complex calculations (e.g. bulk formulae, estimation |
---|
| 42 | of mixing coefficients) that are carried out in the SBC, LDF and ZDF modules and |
---|
| 43 | described in chapters \S\ref{SBC}, \S\ref{LDF} and \S\ref{ZDF}, respectively. |
---|
| 44 | Note that \mdl{tranpc}, the non-penetrative convection module, although |
---|
| 45 | (temporarily) located in the NEMO/OPA/TRA directory, is described with the |
---|
| 46 | model vertical physics (ZDF). |
---|
| 47 | %%% |
---|
| 48 | \gmcomment{change the position of eosbn2 in the reference code} |
---|
| 49 | %%% |
---|
[707] | 50 | |
---|
[817] | 51 | In the present chapter we also describe the diagnostic equations used to compute |
---|
| 52 | the sea-water properties (density, Brunt-Vais\"{a}l\"{a} frequency, specific heat and |
---|
| 53 | freezing point) although the associated modules ($i.e.$ \mdl{eosbn2}, \mdl{ocfzpt} |
---|
| 54 | and \mdl{phycst}) are (temporarily) located in the NEMO/OPA directory. |
---|
[707] | 55 | |
---|
[817] | 56 | The different options available to the user are managed by namelist logical or |
---|
[1224] | 57 | CPP keys. For each equation term \textit{ttt}, the namelist logicals are \textit{ln\_trattt\_xxx}, |
---|
[817] | 58 | where \textit{xxx} is a 3 or 4 letter acronym accounting for each optional scheme. |
---|
| 59 | The CPP key (when it exists) is \textbf{key\_trattt}. The corresponding code can be |
---|
| 60 | found in the \textit{trattt} or \textit{trattt\_xxx} module, in the NEMO/OPA/TRA directory. |
---|
[707] | 61 | |
---|
[817] | 62 | The user has the option of extracting each tendency term on the rhs of the tracer |
---|
| 63 | equation for output (\key{trdtra} is defined), as described in Chap.~\ref{MISC}. |
---|
[707] | 64 | |
---|
[1224] | 65 | $\ $\newline % force a new ligne |
---|
[707] | 66 | % ================================================================ |
---|
| 67 | % Tracer Advection |
---|
| 68 | % ================================================================ |
---|
[817] | 69 | \section [Tracer Advection (\textit{traadv})] |
---|
| 70 | {Tracer Advection (\mdl{traadv})} |
---|
[707] | 71 | \label{TRA_adv} |
---|
| 72 | %------------------------------------------nam_traadv----------------------------------------------------- |
---|
| 73 | \namdisplay{nam_traadv} |
---|
| 74 | %------------------------------------------------------------------------------------------------------------- |
---|
| 75 | |
---|
[817] | 76 | The advection tendency of a tracer in flux form is the divergence of the advective |
---|
[707] | 77 | fluxes. Its discrete expression is given by : |
---|
| 78 | \begin{equation} \label{Eq_tra_adv} |
---|
[1224] | 79 | ADV_\tau =-\frac{1}{b_T} \left( |
---|
| 80 | \;\delta _i \left[ e_{2u}\,e_{3u} \; u\; \tau _u \right] |
---|
| 81 | +\delta _j \left[ e_{1v}\,e_{3v} \; v\; \tau _v \right] \; \right) |
---|
| 82 | -\frac{1}{e_{3T}} \;\delta _k \left[ w\; \tau _w \right] |
---|
[707] | 83 | \end{equation} |
---|
[1224] | 84 | where $\tau$ is either T or S, and $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells. |
---|
| 85 | In pure $z$-coordinate (\key{zco} is defined), it reduces to : |
---|
[707] | 86 | \begin{equation} \label{Eq_tra_adv_zco} |
---|
| 87 | ADV_\tau =-\frac{1}{e_{1T} {\kern 1pt}e_{2T} {\kern 1pt}}\left( {\;\delta _i |
---|
| 88 | \left[ {e_{2u} {\kern 1pt}{\kern 1pt}\;u\;\tau _u } \right]+\delta _j \left[ |
---|
| 89 | {e_{1v} {\kern 1pt}v\;\tau _v } \right]\;} \right)-\frac{1}{\mathop |
---|
| 90 | e\nolimits_{3T} }\delta _k \left[ {w\;\tau _w } \right] |
---|
| 91 | \end{equation} |
---|
[1224] | 92 | since the vertical scale factors are functions of $k$ only, and thus |
---|
| 93 | $e_{3u} =e_{3v} =e_{3T} $. The flux form in \eqref{Eq_tra_adv} |
---|
| 94 | requires implicitly the use of the continuity equation. Indeed, it is obtained |
---|
| 95 | by using the following equality : $\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$ |
---|
| 96 | which results from the use of the continuity equation, $\nabla \cdot \vect{U}=0$ or |
---|
| 97 | $\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ in constant (default option) |
---|
| 98 | or variable (\key{vvl} defined) volume case, respectively. |
---|
| 99 | Therefore it is of paramount importance to design the discrete analogue of the |
---|
| 100 | advection tendency so that it is consistent with the continuity equation in order to |
---|
[817] | 101 | enforce the conservation properties of the continuous equations. In other words, |
---|
| 102 | by substituting $\tau$ by 1 in (\ref{Eq_tra_adv}) we recover the discrete form of |
---|
| 103 | the continuity equation which is used to calculate the vertical velocity. |
---|
[707] | 104 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 105 | \begin{figure}[!t] \label{Fig_adv_scheme} \begin{center} |
---|
[998] | 106 | \includegraphics[width=0.9\textwidth]{./TexFiles/Figures/Fig_adv_scheme.pdf} |
---|
[817] | 107 | \caption{Schematic representation of some ways used to evaluate the tracer value |
---|
| 108 | at $u$-point and the amount of tracer exchanged between two neighbouring grid |
---|
| 109 | points. Upsteam biased scheme (ups): the upstream value is used and the black |
---|
| 110 | area is exchanged. Piecewise parabolic method (ppm): a parabolic interpolation |
---|
| 111 | is used and the black and dark grey areas are exchanged. Monotonic upstream |
---|
| 112 | scheme for conservative laws (muscl): a parabolic interpolation is used and black, |
---|
| 113 | dark grey and grey areas are exchanged. Second order scheme (cen2): the mean |
---|
| 114 | value is used and black, dark grey, grey and light grey areas are exchanged. Note |
---|
| 115 | that this illustration does not include the flux limiter used in ppm and muscl schemes.} |
---|
[707] | 116 | \end{center} \end{figure} |
---|
| 117 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
[817] | 118 | |
---|
| 119 | The key difference between the advection schemes used in \NEMO is the choice |
---|
| 120 | made in space and time interpolation to define the value of the tracer at the |
---|
| 121 | velocity points (Fig.~\ref{Fig_adv_scheme}). |
---|
| 122 | |
---|
[707] | 123 | Along solid lateral and bottom boundaries a zero tracer flux is naturally |
---|
| 124 | specified, since the normal velocity is zero there. At the sea surface the |
---|
[817] | 125 | boundary condition depends on the type of sea surface chosen: |
---|
| 126 | \begin{description} |
---|
| 127 | \item [rigid-lid formulation:] $w=0$ at the surface, so the advective |
---|
| 128 | fluxes through the surface are zero. |
---|
| 129 | \item [linear free surface:] the first level thickness is constant in time: |
---|
| 130 | the vertical boundary condition is applied at the fixed surface $z=0$ |
---|
| 131 | rather than on the moving surface $z=\eta$. There is a non-zero advective |
---|
| 132 | flux which is set for all advection schemes as the product of surface |
---|
| 133 | velocity (at $z=0$) by the first level tracer value: |
---|
| 134 | $\left. {\tau _w } \right|_{k=1/2} =T_{k=1} $. |
---|
| 135 | \item [non-linear free surface:] (\key{vvl} is defined) |
---|
| 136 | convergence/divergence in the first ocean level moves the free surface |
---|
| 137 | up/down. There is no tracer advection through it so that the advective |
---|
| 138 | fluxes through the surface are also zero |
---|
| 139 | \end{description} |
---|
| 140 | In all cases, this boundary condition retains local conservation of tracer. |
---|
| 141 | Global conservation is obtained in both rigid-lid and non-linear free surface |
---|
| 142 | cases, but not in the linear free surface case. Nevertheless, in the latter |
---|
| 143 | case, it is achieved to a good approximation since the non-conservative |
---|
| 144 | term is the product of the time derivative of the tracer and the free surface |
---|
| 145 | height, two quantities that are not correlated (see \S\ref{PE_free_surface}, |
---|
| 146 | and also \citet{Roullet2000,Griffies2001,Campin2004}). |
---|
[707] | 147 | |
---|
[817] | 148 | The velocity field that appears in (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_zco}) |
---|
| 149 | is the centred (\textit{now}) \textit{eulerian} ocean velocity (see Chap.~\ref{DYN}). |
---|
| 150 | When advective bottom boundary layer (\textit{bbl}) and/or eddy induced velocity |
---|
| 151 | (\textit{eiv}) parameterisations are used it is the \textit{now} \textit{effective} |
---|
| 152 | velocity ($i.e.$ the sum of the eulerian, the bbl and/or the eiv velocities) which is used. |
---|
[707] | 153 | |
---|
| 154 | The choice of an advection scheme is made in the \np{nam\_traadv} namelist, by |
---|
| 155 | setting to \textit{true} one and only one of the logicals \textit{ln\_traadv\_xxx}. The |
---|
[817] | 156 | corresponding code can be found in the \textit{traadv\_xxx.F90} module, where |
---|
| 157 | \textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme. Details |
---|
[707] | 158 | of the advection schemes are given below. The choice of an advection scheme |
---|
| 159 | is a complex matter which depends on the model physics, model resolution, |
---|
| 160 | type of tracer, as well as the issue of numerical cost. |
---|
| 161 | |
---|
[817] | 162 | Note that |
---|
| 163 | (1) cen2, cen4 and TVD schemes require an explicit diffusion |
---|
[707] | 164 | operator while the other schemes are diffusive enough so that they do not |
---|
[817] | 165 | require additional diffusion ; |
---|
| 166 | (2) cen2, cen4, MUSCL2, and UBS are not \textit{positive} schemes |
---|
| 167 | \footnote{negative values can appear in an initially strictly positive tracer field |
---|
| 168 | which is advected} |
---|
| 169 | , implying that false extrema are permitted. Their use is not recommended on passive tracers ; |
---|
| 170 | (3) It is highly recommended that the same advection-diffusion scheme is |
---|
| 171 | used on both active and passive tracers. Indeed, if a source or sink of a |
---|
| 172 | passive tracer depends on an active one, the difference of treatment of |
---|
| 173 | active and passive tracers can create very nice-looking frontal structures |
---|
| 174 | that are pure numerical artefacts. |
---|
[707] | 175 | |
---|
| 176 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 177 | % 2nd order centred scheme |
---|
| 178 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 179 | \subsection [$2^{nd}$ order centred scheme (cen2) (\np{ln\_traadv\_cen2})] |
---|
| 180 | {$2^{nd}$ order centred scheme (cen2) (\np{ln\_traadv\_cen2}=.true.)} |
---|
[707] | 181 | \label{TRA_adv_cen2} |
---|
| 182 | |
---|
| 183 | In the centred second order formulation, the tracer at velocity points is |
---|
[817] | 184 | evaluated as the mean of the two neighbouring $T$-point values. |
---|
| 185 | For example, in the $i$-direction : |
---|
[707] | 186 | \begin{equation} \label{Eq_tra_adv_cen2} |
---|
| 187 | \tau _u^{cen2} =\overline T ^{i+1/2} |
---|
| 188 | \end{equation} |
---|
| 189 | |
---|
[817] | 190 | The scheme is non diffusive ($i.e.$ it conserves the tracer variance, $\tau^2)$ |
---|
| 191 | but dispersive ($i.e.$ it may create false extrema). It is therefore notoriously |
---|
| 192 | noisy and must be used in conjunction with an explicit diffusion operator to |
---|
| 193 | produce a sensible solution. The associated time-stepping is performed using |
---|
| 194 | a leapfrog scheme in conjunction with an Asselin time-filter, so $T$ in |
---|
[1224] | 195 | (\ref{Eq_tra_adv_cen2}) is the \textit{now} tracer value. The centered second |
---|
| 196 | order advection is computed in the \mdl{traadv\_cen2} module. In this module, |
---|
| 197 | it is also proposed to combine the \textit{cen2} scheme with an upstream scheme |
---|
| 198 | in specific areas which requires a strong diffusion in order to avoid the generation |
---|
| 199 | of false extrema. These areas are the vicinity of large river mouths, some straits |
---|
| 200 | with coarse resolution, and the vicinity of ice cover area ($i.e.$ when the ocean |
---|
| 201 | temperature is close to the freezing point). |
---|
[707] | 202 | |
---|
[817] | 203 | Note that using the cen2 scheme, the overall tracer advection is of second |
---|
| 204 | order accuracy since both (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_cen2}) |
---|
[1224] | 205 | have this order of accuracy. Note also that |
---|
[707] | 206 | |
---|
| 207 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 208 | % 4nd order centred scheme |
---|
| 209 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 210 | \subsection [$4^{nd}$ order centred scheme (cen4) (\np{ln\_traadv\_cen4})] |
---|
[994] | 211 | {$4^{nd}$ order centred scheme (cen4) (\np{ln\_traadv\_cen4}=.true.)} |
---|
[707] | 212 | \label{TRA_adv_cen4} |
---|
| 213 | |
---|
[817] | 214 | In the $4^{th}$ order formulation (to be implemented), tracer values are |
---|
| 215 | evaluated at velocity points as a $4^{th}$ order interpolation, and thus uses |
---|
| 216 | the four neighbouring $T$-points. For example, in the $i$-direction: |
---|
[707] | 217 | \begin{equation} \label{Eq_tra_adv_cen4} |
---|
| 218 | \tau _u^{cen4} |
---|
| 219 | =\overline{ T - \frac{1}{6}\,\delta _i \left[ \delta_{i+1/2}[T] \,\right] }^{\,i+1/2} |
---|
| 220 | \end{equation} |
---|
| 221 | |
---|
| 222 | Strictly speaking, the cen4 scheme is not a $4^{th}$ order advection scheme |
---|
[817] | 223 | but a $4^{th}$ order evaluation of advective fluxes, since the divergence of |
---|
| 224 | advective fluxes \eqref{Eq_tra_adv} is kept at $2^{nd}$ order. The phrase ``$4^{th}$ |
---|
| 225 | order scheme'' used in oceanographic literature is usually associated |
---|
[707] | 226 | with the scheme presented here. Introducing a \textit{true} $4^{th}$ order advection |
---|
| 227 | scheme is feasible but, for consistency reasons, it requires changes in the |
---|
| 228 | discretisation of the tracer advection together with changes in both the |
---|
[817] | 229 | continuity equation and the momentum advection terms. |
---|
[707] | 230 | |
---|
| 231 | A direct consequence of the pseudo-fourth order nature of the scheme is that |
---|
[1224] | 232 | it is not non-diffusive, i.e. the global variance of a tracer is not preserved using |
---|
| 233 | \textit{cen4}. Furthermore, it must be used in conjunction with an explicit |
---|
| 234 | diffusion operator to produce a sensible solution. The time-stepping is also |
---|
| 235 | performed using a leapfrog scheme in conjunction with an Asselin time-filter, |
---|
| 236 | so $T$ in (\ref{Eq_tra_adv_cen4}) is the \textit{now} tracer. |
---|
[707] | 237 | |
---|
[817] | 238 | At a $T$-grid cell adjacent to a boundary (coastline, bottom and surface), an |
---|
[707] | 239 | additional hypothesis must be made to evaluate $\tau _u^{cen4}$. This |
---|
| 240 | hypothesis usually reduces the order of the scheme. Here we choose to set |
---|
| 241 | the gradient of $T$ across the boundary to zero. Alternative conditions can be |
---|
[817] | 242 | specified, such as a reduction to a second order scheme for these near boundary |
---|
| 243 | grid points. |
---|
[707] | 244 | |
---|
| 245 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 246 | % TVD scheme |
---|
| 247 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 248 | \subsection [Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd})] |
---|
| 249 | {Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd}=.true.)} |
---|
[707] | 250 | \label{TRA_adv_tvd} |
---|
| 251 | |
---|
| 252 | In the Total Variance Dissipation (TVD) formulation, the tracer at velocity |
---|
[1224] | 253 | points is evaluated using a combination of an upstream and a centred scheme. |
---|
| 254 | For example, in the $i$-direction : |
---|
[707] | 255 | \begin{equation} \label{Eq_tra_adv_tvd} |
---|
| 256 | \begin{split} |
---|
| 257 | \tau _u^{ups}&= \begin{cases} |
---|
| 258 | T_{i+1} & \text{if $\ u_{i+1/2} < 0$} \hfill \\ |
---|
| 259 | T_i & \text{if $\ u_{i+1/2} \geq 0$} \hfill \\ |
---|
| 260 | \end{cases} \\ |
---|
| 261 | \\ |
---|
| 262 | \tau _u^{tvd}&=\tau _u^{ups} +c_u \;\left( {\tau _u^{cen2} -\tau _u^{ups} } \right) |
---|
| 263 | \end{split} |
---|
| 264 | \end{equation} |
---|
[1224] | 265 | where $c_u$ is a flux limiter function taking values between 0 and 1. |
---|
| 266 | There exist many ways to define $c_u$, each correcponding to a different |
---|
| 267 | total variance decreasing scheme. The one chosen in \NEMO is described in |
---|
[996] | 268 | \citet{Zalesak1979}. $c_u$ only departs from $1$ when the advective term |
---|
[817] | 269 | produces a local extremum in the tracer field. The resulting scheme is quite |
---|
| 270 | expensive but \emph{positive}. It can be used on both active and passive tracers. |
---|
| 271 | This scheme is tested and compared with MUSCL and the MPDATA scheme in |
---|
| 272 | \citet{Levy2001}; note that in this paper it is referred to as "FCT" (Flux corrected |
---|
[1224] | 273 | transport) rather than TVD. The TVD scheme is computed in the \mdl{traadv\_tvd} module. |
---|
[707] | 274 | |
---|
[817] | 275 | For stability reasons (see \S\ref{DOM_nxt}), in (\ref{Eq_tra_adv_tvd}) |
---|
| 276 | $\tau _u^{cen2}$ is evaluated using the \textit{now} tracer while $\tau _u^{ups}$ |
---|
| 277 | is evaluated using the \textit{before} tracer. In other words, the advective part of |
---|
| 278 | the scheme is time stepped with a leap-frog scheme while a forward scheme is |
---|
| 279 | used for the diffusive part. |
---|
[707] | 280 | |
---|
| 281 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 282 | % MUSCL scheme |
---|
| 283 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 284 | \subsection[MUSCL scheme (\np{ln\_traadv\_muscl})] |
---|
[707] | 285 | {Monotone Upstream Scheme for Conservative Laws (MUSCL) (\np{ln\_traadv\_muscl}=T)} |
---|
| 286 | \label{TRA_adv_muscl} |
---|
| 287 | |
---|
| 288 | The Monotone Upstream Scheme for Conservative Laws (MUSCL) has been |
---|
[817] | 289 | implemented by \citet{Levy2001}. In its formulation, the tracer at velocity points |
---|
| 290 | is evaluated assuming a linear tracer variation between two $T$-points |
---|
| 291 | (Fig.\ref{Fig_adv_scheme}). For example, in the $i$-direction : |
---|
[707] | 292 | \begin{equation} \label{Eq_tra_adv_muscl} |
---|
| 293 | \tau _u^{mus} = \left\{ \begin{aligned} |
---|
| 294 | &\tau _i &+ \frac{1}{2} \;\left( 1-\frac{u_{i+1/2} \;\Delta t}{e_{1u}} \right) |
---|
| 295 | &\ \widetilde{\partial _i \tau} & \quad \text{if }\;u_{i+1/2} \geqslant 0 \\ |
---|
| 296 | &\tau _{i+1/2} &+\frac{1}{2}\;\left( 1+\frac{u_{i+1/2} \;\Delta t}{e_{1u} } \right) |
---|
| 297 | &\ \widetilde{\partial_{i+1/2} \tau } & \text{if }\;u_{i+1/2} <0 |
---|
| 298 | \end{aligned} \right. |
---|
| 299 | \end{equation} |
---|
[817] | 300 | where $\widetilde{\partial _i \tau}$ is the slope of the tracer on which a limitation |
---|
| 301 | is imposed to ensure the \textit{positive} character of the scheme. |
---|
[707] | 302 | |
---|
[817] | 303 | The time stepping is performed using a forward scheme, that is the \textit{before} |
---|
| 304 | tracer field is used to evaluate $\tau _u^{mus}$. |
---|
[707] | 305 | |
---|
[817] | 306 | For an ocean grid point adjacent to land and where the ocean velocity is |
---|
| 307 | directed toward land, two choices are available: an upstream flux |
---|
| 308 | (\np{ln\_traadv\_muscl}=.true.) or a second order flux |
---|
| 309 | (\np{ln\_traadv\_muscl2}=.true.). Note that the latter choice does not ensure |
---|
| 310 | the \textit{positive} character of the scheme. Only the former can be used |
---|
[1224] | 311 | on both active and passive tracers. The two MUSCL schemes are computed |
---|
| 312 | in the \mdl{traadv\_tvd} and \mdl{traadv\_tvd2} modules. |
---|
[707] | 313 | |
---|
| 314 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 315 | % UBS scheme |
---|
| 316 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 317 | \subsection [Upstream-Biased Scheme (UBS) (\np{ln\_traadv\_ubs})] |
---|
| 318 | {Upstream-Biased Scheme (UBS) (\np{ln\_traadv\_ubs}=.true.)} |
---|
[707] | 319 | \label{TRA_adv_ubs} |
---|
| 320 | |
---|
[817] | 321 | The UBS advection scheme is an upstream-biased third order scheme based on |
---|
| 322 | an upstream-biased parabolic interpolation. It is also known as the Cell |
---|
[707] | 323 | Averaged QUICK scheme (Quadratic Upstream Interpolation for Convective |
---|
| 324 | Kinematics). For example, in the $i$-direction : |
---|
| 325 | \begin{equation} \label{Eq_tra_adv_ubs} |
---|
| 326 | \tau _u^{ubs} =\overline T ^{i+1/2}-\;\frac{1}{6} \left\{ |
---|
| 327 | \begin{aligned} |
---|
| 328 | &\tau"_i & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\ |
---|
| 329 | &\tau"_{i+1} & \quad \text{if }\ u_{i+1/2} < 0 |
---|
| 330 | \end{aligned} \right. |
---|
| 331 | \end{equation} |
---|
| 332 | where $\tau "_i =\delta _i \left[ {\delta _{i+1/2} \left[ \tau \right]} \right]$. |
---|
| 333 | |
---|
| 334 | This results in a dissipatively dominant (i.e. hyper-diffusive) truncation |
---|
[1224] | 335 | error \citep{Sacha2005}. The overall performance of the advection |
---|
| 336 | scheme is similar to that reported in \cite{Farrow1995}. |
---|
| 337 | It is a relatively good compromise between accuracy and smoothness. |
---|
| 338 | It is not a \emph{positive} scheme, meaning that false extrema are permitted, |
---|
| 339 | but the amplitude of such are significantly reduced over the centred second |
---|
| 340 | order method. Nevertheless it is not recommended that it should be applied |
---|
| 341 | to a passive tracer that requires positivity. |
---|
[707] | 342 | |
---|
| 343 | The intrinsic diffusion of UBS makes its use risky in the vertical direction |
---|
| 344 | where the control of artificial diapycnal fluxes is of paramount importance. |
---|
[1224] | 345 | Therefore the vertical flux is evaluated using the TVD scheme when |
---|
| 346 | \np{ln\_traadv\_ubs}=.true.. |
---|
[707] | 347 | |
---|
[817] | 348 | For stability reasons (see \S\ref{DOM_nxt}), in \eqref{Eq_tra_adv_ubs}, |
---|
| 349 | the first term (which corresponds to a second order centred scheme) |
---|
| 350 | is evaluated using the \textit{now} tracer (centred in time) while the |
---|
| 351 | second term (which is the diffusive part of the scheme), is |
---|
| 352 | evaluated using the \textit{before} tracer (forward in time). |
---|
[1224] | 353 | This choice is discussed by \citet{Webb1998} in the context of the |
---|
| 354 | QUICK advection scheme. UBS and QUICK schemes only differ |
---|
| 355 | by one coefficient. Replacing 1/6 with 1/8 in \eqref{Eq_tra_adv_ubs} |
---|
| 356 | leads to the QUICK advection scheme \citep{Webb1998}. |
---|
| 357 | This option is not available through a namelist parameter, since the |
---|
| 358 | 1/6 coefficient is hard coded. Nevertheless it is quite easy to make the |
---|
| 359 | substitution in the \mdl{traadv\_ubs} module and obtain a QUICK scheme. |
---|
[707] | 360 | |
---|
[817] | 361 | Note that : |
---|
| 362 | |
---|
[1224] | 363 | (1) When a high vertical resolution $O(1m)$ is used, the model stability can |
---|
[707] | 364 | be controlled by vertical advection (not vertical diffusion which is usually |
---|
| 365 | solved using an implicit scheme). Computer time can be saved by using a |
---|
[1224] | 366 | time-splitting technique on vertical advection. Such a technique has been |
---|
| 367 | implemented and validated in ORCA05 with 301 levels. It is not available |
---|
| 368 | in the current reference version. |
---|
[707] | 369 | |
---|
[1224] | 370 | (2) In a forthcoming release four options will be available for the vertical |
---|
[817] | 371 | component used in the UBS scheme. $\tau _w^{ubs}$ will be evaluated |
---|
[1224] | 372 | using either \textit{(a)} a centred $2^{nd}$ order scheme, or \textit{(b)} |
---|
[817] | 373 | a TVD scheme, or \textit{(c)} an interpolation based on conservative |
---|
| 374 | parabolic splines following the \citet{Sacha2005} implementation of UBS |
---|
| 375 | in ROMS, or \textit{(d)} a UBS. The $3^{rd}$ case has dispersion properties |
---|
| 376 | similar to an eighth-order accurate conventional scheme. |
---|
| 377 | |
---|
[1224] | 378 | (3) It is straightforward to rewrite \eqref{Eq_tra_adv_ubs} as follows: |
---|
| 379 | \begin{equation} \label{Eq_traadv_ubs2} |
---|
| 380 | \tau _u^{ubs} = \tau _u^{cen4} + \frac{1}{12} \left\{ |
---|
| 381 | \begin{aligned} |
---|
| 382 | & + \tau"_i & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\ |
---|
| 383 | & - \tau"_{i+1} & \quad \text{if }\ u_{i+1/2} < 0 |
---|
| 384 | \end{aligned} \right. |
---|
[707] | 385 | \end{equation} |
---|
| 386 | or equivalently |
---|
[1224] | 387 | \begin{equation} \label{Eq_traadv_ubs2b} |
---|
[707] | 388 | u_{i+1/2} \ \tau _u^{ubs} |
---|
| 389 | =u_{i+1/2} \ \overline{ T - \frac{1}{6}\,\delta _i\left[ \delta_{i+1/2}[T] \,\right] }^{\,i+1/2} |
---|
| 390 | - \frac{1}{2} |u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i] |
---|
| 391 | \end{equation} |
---|
[1224] | 392 | |
---|
| 393 | \eqref{Eq_traadv_ubs2} has several advantages. Firstly, it clearly reveals |
---|
[817] | 394 | that the UBS scheme is based on the fourth order scheme to which an |
---|
| 395 | upstream-biased diffusion term is added. Secondly, this emphasises that the |
---|
| 396 | $4^{th}$ order part (as well as the $2^{nd}$ order part as stated above) has |
---|
| 397 | to be evaluated at the \emph{now} time step using \eqref{Eq_tra_adv_ubs}. |
---|
| 398 | Thirdly, the diffusion term is in fact a biharmonic operator with an eddy |
---|
| 399 | coefficient which is simply proportional to the velocity: |
---|
| 400 | $A_u^{lm}= - \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note that NEMO v2.3 still uses |
---|
[1224] | 401 | \eqref{Eq_tra_adv_ubs}, not \eqref{Eq_traadv_ubs2}. This should be |
---|
[817] | 402 | changed in forthcoming release. |
---|
| 403 | %%% |
---|
| 404 | \gmcomment{the change in UBS scheme has to be done} |
---|
| 405 | %%% |
---|
[707] | 406 | |
---|
| 407 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 408 | % QCK scheme |
---|
| 409 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 410 | \subsection [QUICKEST scheme (QCK) (\np{ln\_traadv\_qck})] |
---|
| 411 | {QUICKEST scheme (QCK) (\np{ln\_traadv\_qck}=.true.)} |
---|
[707] | 412 | \label{TRA_adv_qck} |
---|
| 413 | |
---|
| 414 | The Quadratic Upstream Interpolation for Convective Kinematics with |
---|
[817] | 415 | Estimated Streaming Terms (QUICKEST) scheme proposed by \citet{Leonard1979} |
---|
| 416 | is the third order Godunov scheme. It is associated with the ULTIMATE QUICKEST |
---|
[707] | 417 | limiter \citep{Leonard1991}. It has been implemented in NEMO by G. Reffray |
---|
[1224] | 418 | (MERCATOR-ocean) and can be found in the \mdl{traadv\_qck} module. |
---|
| 419 | The resulting scheme is quite expensive but \emph{positive}. |
---|
| 420 | It can be used on both active and passive tracers. |
---|
| 421 | Nevertheless, the intrinsic diffusion of QCK makes its use risky in the vertical |
---|
| 422 | direction where the control of artificial diapycnal fluxes is of paramount importance. |
---|
| 423 | Therefore the vertical flux is evaluated using the CEN2 scheme. |
---|
| 424 | This no more ensure the positivity of the scheme. The use of TVD in the vertical |
---|
| 425 | direction as for the UBS case should be implemented to maintain the property. |
---|
[707] | 426 | |
---|
| 427 | |
---|
| 428 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 429 | % PPM scheme |
---|
| 430 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 431 | \subsection [Piecewise Parabolic Method (PPM) (\np{ln\_traadv\_ppm})] |
---|
| 432 | {Piecewise Parabolic Method (PPM) (\np{ln\_traadv\_ppm}=.true.)} |
---|
[707] | 433 | \label{TRA_adv_ppm} |
---|
| 434 | |
---|
| 435 | The Piecewise Parabolic Method (PPM) proposed by Colella and Woodward (1984) |
---|
[817] | 436 | is based on a quadradic piecewise rebuilding. Like the QCK scheme, it is associated |
---|
| 437 | with the ULTIMATE QUICKEST limiter \citep{Leonard1991}. It has been implemented |
---|
| 438 | in \NEMO by G. Reffray (MERCATOR-ocean) but is not yet offered in the reference |
---|
[1224] | 439 | version 3.0. |
---|
[707] | 440 | |
---|
| 441 | % ================================================================ |
---|
| 442 | % Tracer Lateral Diffusion |
---|
| 443 | % ================================================================ |
---|
[817] | 444 | \section [Tracer Lateral Diffusion (\textit{traldf})] |
---|
| 445 | {Tracer Lateral Diffusion (\mdl{traldf})} |
---|
[707] | 446 | \label{TRA_ldf} |
---|
| 447 | %-----------------------------------------nam_traldf------------------------------------------------------ |
---|
| 448 | \namdisplay{nam_traldf} |
---|
| 449 | %------------------------------------------------------------------------------------------------------------- |
---|
| 450 | |
---|
[817] | 451 | The options available for lateral diffusion are a laplacian (rotated or not) |
---|
| 452 | or a biharmonic operator, the latter being more scale-selective (more |
---|
[707] | 453 | diffusive at small scales). The specification of eddy diffusivity |
---|
[817] | 454 | coefficients (either constant or variable in space and time) as well as the |
---|
| 455 | computation of the slope along which the operators act, are performed in the |
---|
[1224] | 456 | \mdl{ldftra} and \mdl{ldfslp} modules, respectively. This is described in Chap.~\ref{LDF}. |
---|
| 457 | The lateral diffusion of tracers is evaluated using a forward scheme, |
---|
[817] | 458 | $i.e.$ the tracers appearing in its expression are the \textit{before} tracers in time, |
---|
| 459 | except for the pure vertical component that appears when a rotation tensor |
---|
| 460 | is used. This latter term is solved implicitly together with the |
---|
| 461 | vertical diffusion term (see \S\ref{DOM_nxt}). |
---|
[707] | 462 | |
---|
| 463 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 464 | % Iso-level laplacian operator |
---|
| 465 | % ------------------------------------------------------------------------------------------------------------- |
---|
[1224] | 466 | \subsection [Iso-level laplacian operator (lap) (\np{ln\_traldf\_lap})] |
---|
| 467 | {Iso-level laplacian operator (lap) (\np{ln\_traldf\_lap}=.true.) } |
---|
[707] | 468 | \label{TRA_ldf_lap} |
---|
| 469 | |
---|
[1224] | 470 | A laplacian diffusion operator ($i.e.$ a harmonic operator) acting along the model |
---|
[817] | 471 | surfaces is given by: |
---|
[707] | 472 | \begin{equation} \label{Eq_tra_ldf_lap} |
---|
[1224] | 473 | D_T^{lT} =\frac{1}{b_T} \left( \; |
---|
| 474 | \delta _{i}\left[ A_u^{lT} \; \frac{e_{2u}\,e_{3u}}{e_{1u}} \;\delta _{i+1/2} [T] \right] |
---|
| 475 | + \delta _{j}\left[ A_v^{lT} \; \frac{e_{1v}\,e_{3v}}{e_{2v}} \;\delta _{j+1/2} [T] \right] \;\right) |
---|
[707] | 476 | \end{equation} |
---|
[1224] | 477 | where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells. |
---|
| 478 | It can be found in the \mdl{traadv\_lap} module. |
---|
[707] | 479 | |
---|
[1224] | 480 | This lateral operator is computed in \mdl{traldf\_lap}. It is a \emph{horizontal} |
---|
| 481 | operator ($i.e.$ acting along geopotential surfaces) in the $z$-coordinate with |
---|
| 482 | or without partial step, but is simply an iso-level operator in the $s$-coordinate. |
---|
[817] | 483 | It is thus used when, in addition to \np{ln\_traldf\_lap}=.true., we have |
---|
[1224] | 484 | \np{ln\_traldf\_level}=.true., or \np{ln\_traldf\_hor}=\np{ln\_zco}=.true.. |
---|
| 485 | In both cases, it significantly contributes to diapycnal mixing. |
---|
| 486 | It is therefore not recommended. |
---|
[707] | 487 | |
---|
[817] | 488 | Note that |
---|
[1224] | 489 | (a) In pure $z$-coordinate (\key{zco} is defined), $e_{3u}=e_{3v}=e_{3T}$, |
---|
| 490 | so that the vertical scale factors disappear from (\ref{Eq_tra_ldf_lap}) ; |
---|
| 491 | (b) In partial step $z$-coordinate (\np{ln\_zps}=.true.), tracers in horizontally |
---|
[817] | 492 | adjacent cells are located at different depths in the vicinity of the bottom. |
---|
| 493 | In this case, horizontal derivatives in (\ref{Eq_tra_ldf_lap}) at the bottom level |
---|
| 494 | require a specific treatment. They are calculated in the \mdl{zpshde} module, |
---|
| 495 | described in \S\ref{TRA_zpshde}. |
---|
[707] | 496 | |
---|
| 497 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 498 | % Rotated laplacian operator |
---|
| 499 | % ------------------------------------------------------------------------------------------------------------- |
---|
[1224] | 500 | \subsection [Rotated laplacian operator (iso) (\np{ln\_traldf\_lap})] |
---|
| 501 | {Rotated laplacian operator (iso) (\np{ln\_traldf\_lap}=.true.)} |
---|
[707] | 502 | \label{TRA_ldf_iso} |
---|
| 503 | |
---|
[817] | 504 | The general form of the second order lateral tracer subgrid scale physics |
---|
| 505 | (\ref{Eq_PE_zdf}) takes the following semi-discrete space form in $z$- and |
---|
| 506 | $s$-coordinates: |
---|
[707] | 507 | \begin{equation} \label{Eq_tra_ldf_iso} |
---|
| 508 | \begin{split} |
---|
[1224] | 509 | D_T^{lT} = \frac{1}{b_T} & \left\{ \,\;\delta_i \left[ A_u^{lT} \left( |
---|
| 510 | \frac{e_{2u}\;e_{3u}}{e_{1u}} \,\delta_{i+1/2}[T] |
---|
| 511 | - e_{2u}\;r_{1u} \,\overline{\overline{ \delta_{k+1/2}[T] }}^{\,i+1/2,k} |
---|
| 512 | \right) \right] \right. \\ |
---|
| 513 | & +\delta_j \left[ A_v^{lT} \left( |
---|
| 514 | \frac{e_{1v}\,e_{3v}}{e_{2v}} \,\delta_{j+1/2} [T] |
---|
| 515 | - e_{1v}\,r_{2v} \,\overline{\overline{ \delta_{k+1/2} [T] }}^{\,j+1/2,k} |
---|
| 516 | \right) \right] \\ |
---|
| 517 | & +\delta_k \left[ A_w^{lT} \left( |
---|
| 518 | -\;e_{2w}\,r_{1w} \,\overline{\overline{ \delta_{i+1/2} [T] }}^{\,i,k+1/2} |
---|
| 519 | \right. \right. \\ |
---|
[707] | 520 | & \qquad \qquad \quad |
---|
[1224] | 521 | - e_{1w}\,r_{2w} \,\overline{\overline{ \delta_{j+1/2} [T] }}^{\,j,k+1/2} \\ |
---|
| 522 | & \left. {\left. { \qquad \qquad \ \ \ \left. { |
---|
| 523 | +\;\frac{e_{1w}\,e_{2w}}{e_{3w}} \,\left( r_{1w}^2 + r_{2w}^2 \right) |
---|
| 524 | \,\delta_{k+1/2} [T] } \right) } \right] \quad } \right\} |
---|
[707] | 525 | \end{split} |
---|
| 526 | \end{equation} |
---|
[1224] | 527 | where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells, |
---|
| 528 | $r_1$ and $r_2$ are the slopes between the surface of computation |
---|
[817] | 529 | ($z$- or $s$-surfaces) and the surface along which the diffusion operator |
---|
| 530 | acts ($i.e.$ horizontal or iso-neutral surfaces). It is thus used when, |
---|
[1224] | 531 | in addition to \np{ln\_traldf\_lap}= .true., we have \np{ln\_traldf\_iso}=.true., |
---|
[817] | 532 | or both \np{ln\_traldf\_hor}=.true. and \np{ln\_zco}=.true.. The way these |
---|
| 533 | slopes are evaluated is given in \S\ref{LDF_slp}. At the surface, bottom |
---|
| 534 | and lateral boundaries, the turbulent fluxes of heat and salt are set to zero |
---|
| 535 | using the mask technique (see \S\ref{LBC_coast}). |
---|
[707] | 536 | |
---|
[817] | 537 | The operator in \eqref{Eq_tra_ldf_iso} involves both lateral and vertical |
---|
| 538 | derivatives. For numerical stability, the vertical second derivative must |
---|
| 539 | be solved using the same implicit time scheme as that used in the vertical |
---|
| 540 | physics (see \S\ref{TRA_zdf}). For computer efficiency reasons, this term |
---|
[1224] | 541 | is not computed in the \mdl{traldf\_iso} module, but in the \mdl{trazdf} module |
---|
[817] | 542 | where, if iso-neutral mixing is used, the vertical mixing coefficient is simply |
---|
| 543 | increased by $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$. |
---|
[707] | 544 | |
---|
[817] | 545 | This formulation conserves the tracer but does not ensure the decrease |
---|
| 546 | of the tracer variance. Nevertheless the treatment performed on the slopes |
---|
| 547 | (see \S\ref{LDF}) allows the model to run safely without any additional |
---|
| 548 | background horizontal diffusion \citep{Guily2001}. An alternative scheme |
---|
[1224] | 549 | developed by \cite{Griffies1998} which preserves both tracer and its variance |
---|
| 550 | is currently been tested in \NEMO. It should be available in a forthcoming |
---|
| 551 | release. |
---|
[707] | 552 | |
---|
[817] | 553 | Note that in the partial step $z$-coordinate (\np{ln\_zps}=.true.), the horizontal |
---|
| 554 | derivatives at the bottom level in \eqref{Eq_tra_ldf_iso} require a specific |
---|
| 555 | treatment. They are calculated in module zpshde, described in \S\ref{TRA_zpshde}. |
---|
[707] | 556 | |
---|
| 557 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 558 | % Iso-level bilaplacian operator |
---|
| 559 | % ------------------------------------------------------------------------------------------------------------- |
---|
[1224] | 560 | \subsection [Iso-level bilaplacian operator (bilap) (\np{ln\_traldf\_bilap})] |
---|
| 561 | {Iso-level bilaplacian operator (bilap) (\np{ln\_traldf\_bilap}=.true.)} |
---|
[707] | 562 | \label{TRA_ldf_bilap} |
---|
| 563 | |
---|
[817] | 564 | The lateral fourth order bilaplacian operator on tracers is obtained by |
---|
| 565 | applying (\ref{Eq_tra_ldf_lap}) twice. It requires an additional assumption |
---|
| 566 | on boundary conditions: the first and third derivative terms normal to the |
---|
[1224] | 567 | coast are set to zero. It is used when, in addition to \np{ln\_traldf\_bilap}=.true., |
---|
| 568 | we have \np{ln\_traldf\_level}=.true., or both \np{ln\_traldf\_hor}=.true. and |
---|
[817] | 569 | \np{ln\_zco}=.false.. In both cases, it can contribute diapycnal mixing, |
---|
| 570 | although less than in the laplacian case. It is therefore not recommended. |
---|
[707] | 571 | |
---|
[817] | 572 | Note that in the code, the bilaplacian routine does not call the laplacian |
---|
[1224] | 573 | routine twice but is rather a separate routine that can be found in the |
---|
| 574 | \mdl{traldf\_bilap} module. This is due to the fact that we introduce the |
---|
| 575 | eddy diffusivity coefficient, A, in the operator as: |
---|
| 576 | $\nabla \cdot \nabla \left( {A\nabla \cdot \nabla T} \right)$, |
---|
| 577 | instead of |
---|
| 578 | $-\nabla \cdot a\nabla \left( {\nabla \cdot a\nabla T} \right)$ |
---|
| 579 | where $a=\sqrt{|A|}$ and $A<0$. This was a mistake: both formulations |
---|
| 580 | ensure the total variance decrease, but the former requires a larger |
---|
| 581 | number of code-lines. It will be corrected in a forthcoming release. |
---|
[707] | 582 | |
---|
| 583 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 584 | % Rotated bilaplacian operator |
---|
| 585 | % ------------------------------------------------------------------------------------------------------------- |
---|
[1224] | 586 | \subsection [Rotated bilaplacian operator (bilapg) (\np{ln\_traldf\_bilap})] |
---|
| 587 | {Rotated bilaplacian operator (bilapg) (\np{ln\_traldf\_bilap}=.true.)} |
---|
[707] | 588 | \label{TRA_ldf_bilapg} |
---|
| 589 | |
---|
[817] | 590 | The lateral fourth order operator formulation on tracers is obtained by |
---|
| 591 | applying (\ref{Eq_tra_ldf_iso}) twice. It requires an additional assumption |
---|
| 592 | on boundary conditions: first and third derivative terms normal to the |
---|
[1224] | 593 | coast, the bottom and the surface are set to zero. It can be found in the |
---|
| 594 | \mdl{traldf\_bilapg}. |
---|
[707] | 595 | |
---|
[1224] | 596 | It is used when, in addition to \np{ln\_traldf\_bilap}=.true., we have |
---|
| 597 | \np{ln\_traldf\_iso}= .true, or both \np{ln\_traldf\_hor}=.true. and \np{ln\_zco}=.true.. |
---|
[817] | 598 | Nevertheless, this rotated bilaplacian operator has never been seriously |
---|
| 599 | tested. No warranties that it is neither free of bugs or correctly formulated. |
---|
| 600 | Moreover, the stability range of such an operator will be probably quite |
---|
| 601 | narrow, requiring a significantly smaller time-step than the one used on |
---|
| 602 | unrotated operator. |
---|
[707] | 603 | |
---|
| 604 | % ================================================================ |
---|
| 605 | % Tracer Vertical Diffusion |
---|
| 606 | % ================================================================ |
---|
[817] | 607 | \section [Tracer Vertical Diffusion (\textit{trazdf})] |
---|
| 608 | {Tracer Vertical Diffusion (\mdl{trazdf})} |
---|
[707] | 609 | \label{TRA_zdf} |
---|
| 610 | %--------------------------------------------namzdf--------------------------------------------------------- |
---|
| 611 | \namdisplay{namzdf} |
---|
| 612 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 613 | |
---|
[817] | 614 | The formulation of the vertical subgrid scale tracer physics is the same |
---|
| 615 | for all the vertical coordinates, and is based on a laplacian operator. |
---|
| 616 | The vertical diffusion operator given by (\ref{Eq_PE_zdf}) takes the |
---|
| 617 | following semi-discrete space form: |
---|
[707] | 618 | \begin{equation} \label{Eq_tra_zdf} |
---|
| 619 | \begin{split} |
---|
[1224] | 620 | D^{vT}_T &= \frac{1}{e_{3T}} \; \delta_k \left[ \;\frac{A^{vT}_w}{e_{3w}} \delta_{k+1/2}[T] \;\right] |
---|
[707] | 621 | \\ |
---|
[1224] | 622 | D^{vS}_T &= \frac{1}{e_{3T}} \; \delta_k \left[ \;\frac{A^{vS}_w}{e_{3w}} \delta_{k+1/2}[S] \;\right] |
---|
[707] | 623 | \end{split} |
---|
| 624 | \end{equation} |
---|
| 625 | where $A_w^{vT}$ and $A_w^{vS}$ are the vertical eddy diffusivity |
---|
[1224] | 626 | coefficients on temperature and salinity, respectively. Generally, |
---|
[817] | 627 | $A_w^{vT}=A_w^{vS}$ except when double diffusive mixing is |
---|
[1224] | 628 | parameterised ($i.e.$ \key{zdfddm} is defined). The way these coefficients |
---|
[817] | 629 | are evaluated is given in \S\ref{ZDF} (ZDF). Furthermore, when |
---|
| 630 | iso-neutral mixing is used, both mixing coefficients are increased |
---|
| 631 | by $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$ |
---|
| 632 | to account for the vertical second derivative of \eqref{Eq_tra_ldf_iso}. |
---|
[707] | 633 | |
---|
[817] | 634 | At the surface and bottom boundaries, the turbulent fluxes of |
---|
[1224] | 635 | heat and salt must be specified. At the surface they are prescribed |
---|
| 636 | from the surface forcing and added in a dedicated routine (see \S\ref{TRA_sbc}), |
---|
[817] | 637 | whilst at the bottom they are set to zero for heat and salt unless |
---|
| 638 | a geothermal flux forcing is prescribed as a bottom boundary |
---|
[1224] | 639 | condition (see \S\ref{TRA_bbc}). |
---|
[707] | 640 | |
---|
| 641 | The large eddy coefficient found in the mixed layer together with high |
---|
[817] | 642 | vertical resolution implies that in the case of explicit time stepping |
---|
| 643 | (\np{ln\_zdfexp}=.true.) there would be too restrictive a constraint on |
---|
| 644 | the time step. Therefore, the default implicit time stepping is preferred |
---|
| 645 | for the vertical diffusion since it overcomes the stability constraint. |
---|
| 646 | A forward time differencing scheme (\np{ln\_zdfexp}=.true.) using a time |
---|
| 647 | splitting technique (\np{n\_zdfexp} $> 1$) is provided as an alternative. |
---|
| 648 | Namelist variables \np{ln\_zdfexp} and \np{n\_zdfexp} apply to both |
---|
| 649 | tracers and dynamics. |
---|
[707] | 650 | |
---|
| 651 | % ================================================================ |
---|
| 652 | % External Forcing |
---|
| 653 | % ================================================================ |
---|
| 654 | \section{External Forcing} |
---|
| 655 | \label{TRA_sbc_qsr_bbc} |
---|
| 656 | |
---|
| 657 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 658 | % surface boundary condition |
---|
| 659 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 660 | \subsection [Surface boundary condition (\textit{trasbc})] |
---|
| 661 | {Surface boundary condition (\mdl{trasbc})} |
---|
[707] | 662 | \label{TRA_sbc} |
---|
| 663 | |
---|
| 664 | The surface boundary condition for tracers is implemented in a separate |
---|
| 665 | module (\mdl{trasbc}) instead of entering as a boundary condition on the vertical |
---|
| 666 | diffusion operator (as in the case of momentum). This has been found to |
---|
| 667 | enhance readability of the code. The two formulations are completely |
---|
| 668 | equivalent; the forcing terms in trasbc are the surface fluxes divided by |
---|
| 669 | the thickness of the top model layer. Following \citet{Roullet2000} the |
---|
[817] | 670 | forcing on an ocean tracer, $c$, can be split into two parts: $F_{ext}$, the |
---|
[707] | 671 | flux of tracer crossing the sea surface and not linked with the water |
---|
[817] | 672 | exchange with the atmosphere, $F_{wf}^d$, and $F_{wf}^i$ the forcing |
---|
| 673 | on the tracer concentration associated with this water exchange. |
---|
| 674 | The latter forcing has two components: a direct effect of change |
---|
| 675 | in concentration associated with the tracer carried by the water flux, |
---|
| 676 | and an indirect concentration/dilution effect : |
---|
[707] | 677 | \begin{equation*} |
---|
| 678 | \begin{split} |
---|
| 679 | F^C &= F_{ext} + F_{wf}^d +F_{wf}^i \\ |
---|
| 680 | &= F_{ext} - \left( c_E \, E - c_p \,P - c_R \,R \right) +c\left( E-P-R \right) |
---|
| 681 | \end{split} |
---|
| 682 | \end{equation*} |
---|
| 683 | |
---|
[817] | 684 | \gmcomment{add here a description of the variable names used in the above equation} |
---|
| 685 | |
---|
| 686 | Two cases must be distinguished, the nonlinear free surface case |
---|
| 687 | (\key{vvl} is defined) and the linear free surface case. The first case |
---|
| 688 | is simpler, because the indirect concentration/dilution effect is naturally |
---|
| 689 | taken into account by letting the vertical scale factors vary in time. |
---|
| 690 | The salinity of water exchanged at the surface is assumed to be zero, |
---|
| 691 | so there is no salt flux at the free surface, except in the presence of sea ice. |
---|
| 692 | The heat flux at the free surface is the sum of $F_{ext}$, the direct |
---|
| 693 | heating/cooling (by the total non-penetrative heat flux) and $F_{wf}^e$ |
---|
| 694 | the heat carried by the water exchanged through the surface (evaporation, |
---|
| 695 | precipitation, runoff). The temperature of precipitation is not well known. |
---|
| 696 | In the model we assume that this water has the same temperature as |
---|
| 697 | the sea surface temperature. The resulting forcing terms for temperature |
---|
| 698 | T and salinity S are: |
---|
[707] | 699 | \begin{equation} \label{Eq_tra_forcing} |
---|
| 700 | \begin{aligned} |
---|
| 701 | &F^T =\frac{Q_{ns} }{\rho _o \;C_p \,e_{3T} }-\frac{\text{EMP}\;\left. T |
---|
| 702 | \right|_{k=1} }{e_{3T} } & \\ |
---|
| 703 | \\ |
---|
| 704 | & F^S =\frac{\text{EMP}_S\;\left. S \right|_{k=1} }{e_{3T} } & |
---|
| 705 | \end{aligned} |
---|
| 706 | \end{equation} |
---|
[817] | 707 | where EMP is the freshwater budget (evaporation minus precipitation |
---|
| 708 | minus river runoff) which forces the ocean volume, $Q_{ns}$ is the |
---|
| 709 | non-penetrative part of the net surface heat flux (difference between |
---|
| 710 | the total surface heat flux and the fraction of the short wave flux that |
---|
| 711 | penetrates into the water column), the product EMP$_S\;.\left. S \right|_{k=1}$ |
---|
| 712 | is the ice-ocean salt flux, and $\left. S\right|_{k=1}$ is the sea surface |
---|
| 713 | salinity (\textit{SSS}). The total salt content is conserved in this formulation |
---|
| 714 | (except for the effect of the Asselin filter). |
---|
[707] | 715 | |
---|
[1224] | 716 | %AMT note: the ice-ocean flux had been forgotten in the first release of the key_vvl option, has this been corrected in the code? ===> gm : NO to be added at NOCS |
---|
[707] | 717 | |
---|
[817] | 718 | In the second case (linear free surface), the vertical scale factors are |
---|
| 719 | fixed in time so that the concentration/dilution effect must be added in |
---|
| 720 | the \mdl{trasbc} module. Because of the hypothesis made for the |
---|
| 721 | temperature of precipitation and runoffs, $F_{wf}^e +F_{wf}^i =0$ |
---|
| 722 | for temperature. The resulting forcing term for temperature is: |
---|
[707] | 723 | \begin{equation} \label{Eq_tra_forcing_q} |
---|
| 724 | F^T=\frac{Q_{ns} }{\rho _o \;C_p \,e_{3T} } |
---|
| 725 | \end{equation} |
---|
| 726 | |
---|
[817] | 727 | The salinity forcing is still given by \eqref{Eq_tra_forcing} but the |
---|
| 728 | definition of EMP$_S$ is different: it is the total surface freshwater |
---|
| 729 | budget (evaporation minus precipitation minus river runoff plus |
---|
| 730 | the rate of change of the sea ice thickness). The total salt content |
---|
| 731 | is not exactly conserved (\citet{Roullet2000}. |
---|
| 732 | See also \S\ref{PE_free_surface}). |
---|
[707] | 733 | |
---|
| 734 | In the case of the rigid lid approximation, the surface salinity forcing $F^s$ |
---|
[817] | 735 | is also expressed by \eqref{Eq_tra_forcing}, but now the global integral of |
---|
| 736 | the product of EMP and S, is not compensated by the advection of fluid |
---|
| 737 | through the top level: this is because in the rigid lid case \textit{w(k=1) = 0} |
---|
| 738 | (in contrast to the linear free surface case). As a result, even if the budget |
---|
| 739 | of \textit{EMP} is zero on average over the whole ocean domain, the |
---|
| 740 | associated salt flux is not, since sea-surface salinity and \textit{EMP} are |
---|
| 741 | intrinsically correlated (high \textit{SSS} are found where evaporation is |
---|
| 742 | strong whilst low \textit{SSS} is usually associated with high precipitation |
---|
| 743 | or river runoff). |
---|
[707] | 744 | |
---|
[817] | 745 | The $Q_{ns} $ and \textit{EMP} fields are defined and updated in the |
---|
| 746 | \mdl{sbcmod} module (see \S\ref{SBC}). |
---|
[707] | 747 | |
---|
| 748 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 749 | % Solar Radiation Penetration |
---|
| 750 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 751 | \subsection [Solar Radiation Penetration (\textit{traqsr})] |
---|
| 752 | {Solar Radiation Penetration (\mdl{traqsr})} |
---|
[707] | 753 | \label{TRA_qsr} |
---|
[817] | 754 | %--------------------------------------------namqsr-------------------------------------------------------- |
---|
[707] | 755 | \namdisplay{namqsr} |
---|
| 756 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 757 | |
---|
[817] | 758 | When the penetrative solar radiation option is used (\np{ln\_flxqsr}=.true.), |
---|
| 759 | the solar radiation penetrates the top few meters of the ocean, otherwise |
---|
| 760 | all the heat flux is absorbed in the first ocean level (\np{ln\_flxqsr}=.false.). |
---|
| 761 | Thus, in the former case a term is added to the time evolution equation of |
---|
| 762 | temperature \eqref{Eq_PE_tra_T} whilst the surface boundary condition is |
---|
| 763 | modified to take into account only the non-penetrative part of the surface |
---|
| 764 | heat flux: |
---|
[707] | 765 | \begin{equation} \label{Eq_PE_qsr} |
---|
| 766 | \begin{split} |
---|
| 767 | \frac{\partial T}{\partial t} &= {\ldots} + \frac{1}{\rho_o\, C_p \,e_3} \; \frac{\partial I}{\partial k} \\ |
---|
| 768 | Q_{ns} &= Q_\text{Total} - Q_{sr} |
---|
| 769 | \end{split} |
---|
| 770 | \end{equation} |
---|
| 771 | |
---|
[1224] | 772 | where $I$ is the downward irradiance. The additional term in \eqref{Eq_PE_qsr} |
---|
| 773 | is discretized as follows: |
---|
[707] | 774 | \begin{equation} \label{Eq_tra_qsr} |
---|
| 775 | \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] |
---|
| 776 | \end{equation} |
---|
| 777 | |
---|
[817] | 778 | A formulation involving two extinction coefficients is assumed for the |
---|
| 779 | downward irradiance $I$ \citep{Paulson1977}: |
---|
[707] | 780 | \begin{equation} \label{Eq_traqsr_iradiance} |
---|
| 781 | I(z) = Q_{sr} \left[Re^{-z / \xi_1} + \left( 1-R\right) e^{-z / \xi_2} \right] |
---|
| 782 | \end{equation} |
---|
| 783 | where $Q_{sr}$ is the penetrative part of the surface heat flux, |
---|
[817] | 784 | $\xi_1$ and $\xi_2$ are two extinction length scales and $R$ |
---|
| 785 | determines the relative contribution of the two terms. |
---|
| 786 | The default values used correspond to a Type I water in Jerlov's [1968] |
---|
| 787 | % |
---|
| 788 | \gmcomment : Jerlov reference to be added |
---|
| 789 | % |
---|
[1224] | 790 | classification: $\xi_1 = 0.35~m$, $\xi_2 = 23~m$ and $R = 0.58$ |
---|
[817] | 791 | (corresponding to \np{xsi1}, \np{xsi2} and \np{rabs} namelist parameters, |
---|
| 792 | respectively). $I$ is masked (no flux through the ocean bottom), |
---|
| 793 | so all the solar radiation that reaches the last ocean level is absorbed |
---|
| 794 | in that level. The trend in \eqref{Eq_tra_qsr} associated with the |
---|
| 795 | penetration of the solar radiation is added to the temperature trend, |
---|
| 796 | and the surface heat flux is modified in routine \mdl{traqsr}. |
---|
| 797 | Note that in the $z$-coordinate, the depth of $T-$levels depends |
---|
| 798 | on the single variable $k$. A one dimensional array of the coefficients |
---|
| 799 | $gdsr(k) = Re^{-z_w (k)/\xi_1} + (1-R)e^{-z_w (k)/\xi_2}$ can then |
---|
| 800 | be computed once and saved in memory. Moreover \textit{nksr}, |
---|
| 801 | the level at which $gdrs$ becomes negligible (less than the |
---|
| 802 | computer precision) is computed once, and the trend associated |
---|
| 803 | with the penetration of the solar radiation is only added until that level. |
---|
| 804 | Finally, note that when the ocean is shallow (< 200~m), part of the |
---|
| 805 | solar radiation can reach the ocean floor. In this case, we have |
---|
| 806 | chosen that all remaining radiation is absorbed in the last ocean |
---|
| 807 | level ($i.e.$ $I_w$ is masked). |
---|
[707] | 808 | |
---|
[817] | 809 | When coupling with a biological model (for example PISCES or LOBSTER), |
---|
| 810 | it is possible to calculate the light attenuation using information from |
---|
| 811 | the biology model. Without biological model, it is still possible to introduce |
---|
| 812 | a horizontal variation of the light attenuation by using the observed ocean |
---|
| 813 | surface color. At the time of writing, the latter has not been implemented |
---|
| 814 | in the reference version. |
---|
| 815 | % |
---|
| 816 | \gmcomment{ {yellow}{case 4 bands and bio-coupling to add !!!} } |
---|
| 817 | % |
---|
[707] | 818 | |
---|
| 819 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 820 | % Bottom Boundary Condition |
---|
| 821 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 822 | \subsection [Bottom Boundary Condition (\textit{trabbc} - \key{bbc})] |
---|
| 823 | {Bottom Boundary Condition (\mdl{trabbc} - \key{bbc})} |
---|
[707] | 824 | \label{TRA_bbc} |
---|
| 825 | %--------------------------------------------nambbc-------------------------------------------------------- |
---|
| 826 | \namdisplay{nambbc} |
---|
| 827 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 828 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 829 | \begin{figure}[!t] \label{Fig_geothermal} \begin{center} |
---|
[998] | 830 | \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_TRA_geoth.pdf} |
---|
[1224] | 831 | \caption{Geothermal Heat flux (in $mW.m^{-2}$) used by \cite{Emile-Geay_Madec_OSD08}. |
---|
| 832 | It is inferred from the age of the sea floor and the formulae of \citet{Stein1992}.} |
---|
[707] | 833 | \end{center} \end{figure} |
---|
| 834 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 835 | |
---|
[817] | 836 | Usually it is assumed that there is no exchange of heat or salt through |
---|
| 837 | the ocean bottom, $i.e.$ a no flux boundary condition is applied on active |
---|
| 838 | tracers at the bottom. This is the default option in \NEMO, and it is |
---|
[1224] | 839 | implemented using the masking technique. However, there is a |
---|
[817] | 840 | non-zero heat flux across the seafloor that is associated with solid |
---|
| 841 | earth cooling. This flux is weak compared to surface fluxes (a mean |
---|
| 842 | global value of $\sim0.1\;W/m^2$ \citep{Stein1992}), but it is |
---|
[1224] | 843 | systematically positive and acts on the densest water masses. |
---|
| 844 | Taking this flux into account in a global ocean model increases |
---|
| 845 | the deepest overturning cell ($i.e.$ the one associated with the Antarctic |
---|
| 846 | Bottom Water) by a few Sverdrups \citep{Emile-Geay_Madec_OSD08}. |
---|
[707] | 847 | |
---|
[817] | 848 | The presence or not of geothermal heating is controlled by the namelist |
---|
| 849 | parameter \np{ngeo\_flux}. If this parameter is set to 1, a constant |
---|
| 850 | geothermal heating is introduced whose value is given by the |
---|
| 851 | \np{ngeo\_flux\_const}, which is also a namelist parameter. If it is set to 2, |
---|
| 852 | a spatially varying geothermal heat flux is introduced which is provided |
---|
| 853 | in the geothermal\_heating.nc NetCDF file (Fig.\ref{Fig_geothermal}). |
---|
[707] | 854 | |
---|
| 855 | % ================================================================ |
---|
| 856 | % Bottom Boundary Layer |
---|
| 857 | % ================================================================ |
---|
[817] | 858 | \section [Bottom Boundary Layer (\textit{trabbl}, \textit{trabbl\_adv} )] |
---|
| 859 | {Bottom Boundary Layer (\mdl{trabbl}, \mdl{trabbl\_adv})} |
---|
[707] | 860 | \label{TRA_bbl} |
---|
| 861 | %--------------------------------------------nambbl--------------------------------------------------------- |
---|
| 862 | \namdisplay{nambbl} |
---|
| 863 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 864 | |
---|
[817] | 865 | In a $z$-coordinate configuration, the bottom topography is represented by a |
---|
[707] | 866 | series of discrete steps. This is not adequate to represent gravity driven |
---|
| 867 | downslope flows. Such flows arise downstream of sills such as the Strait of |
---|
| 868 | Gibraltar, Bab El Mandeb, or Denmark Strait, where dense water formed in |
---|
| 869 | marginal seas flows into a basin filled with less dense water. The amount of |
---|
[817] | 870 | entrainment that occurs in these gravity plumes is critical in determining the |
---|
| 871 | density and volume flux of the densest waters of the ocean, such as |
---|
| 872 | Antarctic Bottom Water, or North Atlantic Deep Water. $z$-coordinate |
---|
| 873 | models tend to overestimate the entrainment, because the gravity flow is |
---|
[707] | 874 | mixed down vertically by convection as it goes ``downstairs'' following the |
---|
| 875 | step topography, sometimes over a thickness much larger than the thickness |
---|
[817] | 876 | of the observed gravity plume. A similar problem occurs in the $s$-coordinate when |
---|
[707] | 877 | the thickness of the bottom level varies in large proportions downstream of |
---|
[817] | 878 | a sill \citep{Willebrand2001}, and the thickness of the plume is not resolved. |
---|
[707] | 879 | |
---|
[1224] | 880 | The idea of the bottom boundary layer (BBL) parameterisation first introduced by |
---|
[707] | 881 | \citet{BeckDos1998} is to allow a direct communication between |
---|
[817] | 882 | two adjacent bottom cells at different levels, whenever the densest water is |
---|
| 883 | located above the less dense water. The communication can be by a diffusive |
---|
| 884 | (diffusive BBL), advective fluxes (advective BBL), or both. In the current |
---|
| 885 | implementation of the BBL, only the tracers are modified, not the velocities. |
---|
| 886 | Furthermore, it only connects ocean bottom cells, and therefore does not include |
---|
| 887 | the improvment proposed by \citet{Campin_Goosse_Tel99}. |
---|
[707] | 888 | |
---|
| 889 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 890 | % Diffusive BBL |
---|
| 891 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 892 | \subsection{Diffusive Bottom Boundary layer (\key{bbl\_diff})} |
---|
[707] | 893 | \label{TRA_bbl_diff} |
---|
| 894 | |
---|
[817] | 895 | When applying sigma-diffusion (\key{trabbl} is defined), the diffusive flux between |
---|
| 896 | two adjacent cells living at the ocean bottom is given by |
---|
| 897 | \begin{equation} \label{Eq_tra_bbl_diff} |
---|
| 898 | {\rm {\bf F}}_\sigma=A_l^\sigma \; \nabla_\sigma T |
---|
| 899 | \end{equation} |
---|
| 900 | with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, |
---|
| 901 | and $A_l^\sigma $ the lateral diffusivity in the BBL. Following \citet{BeckDos1998}, |
---|
| 902 | the latter is prescribed with a spatial dependence, $e.g.$ in the conditional form |
---|
| 903 | \begin{equation} \label{Eq_tra_bbl_coef} |
---|
[707] | 904 | A_l^\sigma (i,j,t)=\left\{ {\begin{array}{l} |
---|
[817] | 905 | A_{bbl} \quad \quad \mbox{if} \quad \nabla_\sigma \rho \cdot \nabla H<0 \\ |
---|
[707] | 906 | \\ |
---|
| 907 | 0\quad \quad \;\,\mbox{otherwise} \\ |
---|
| 908 | \end{array}} \right. |
---|
| 909 | \end{equation} |
---|
[817] | 910 | where $A_{bbl}$ is the BBL diffusivity coefficient, given by the namelist |
---|
| 911 | parameter \np{atrbbl}. $A_{bbl}$ is usually set to a value much larger |
---|
| 912 | than the one used on lateral mixing in open ocean. |
---|
| 913 | Note that in practice, \eqref{Eq_tra_bbl_coef} constraint is applied |
---|
| 914 | separately in the two horizontal directions, and the density gradient in |
---|
| 915 | \eqref{Eq_tra_bbl_coef} is evaluated at $\overline{H}^i$ ($\overline{H}^j$) |
---|
| 916 | using the along bottom mean temperature and salinity. |
---|
[707] | 917 | |
---|
| 918 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 919 | % Advective BBL |
---|
| 920 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 921 | \subsection {Advective Bottom Boundary Layer (\key{bbl\_adv})} |
---|
[707] | 922 | \label{TRA_bbl_adv} |
---|
| 923 | |
---|
| 924 | |
---|
[817] | 925 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 926 | \begin{figure}[!t] \label{Fig_bbl} \begin{center} |
---|
[1224] | 927 | \includegraphics[width=0.8\textwidth]{./TexFiles/Figures/Fig_BBL_adv.pdf} |
---|
[817] | 928 | \caption{Advective Bottom Boundary Layer.} |
---|
| 929 | \end{center} \end{figure} |
---|
| 930 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
[707] | 931 | |
---|
[817] | 932 | %%%gmcomment : this section has to be really written |
---|
| 933 | |
---|
| 934 | The advective BBL is in fact not only an advective one but include a diffusive |
---|
| 935 | component as we chose an upstream scheme to perform the advection within |
---|
| 936 | the BBL. The associated diffusion only act in the stream direction and is |
---|
| 937 | proportional to the velocity. |
---|
| 938 | |
---|
| 939 | When applying sigma-advection (\key{trabbl\_adv} defined), the advective |
---|
| 940 | flux between two adjacent cells living at the ocean bottom is given by |
---|
| 941 | \begin{equation} \label{Eq_tra_bbl_fadv} |
---|
| 942 | {\rm {\bf F}}_\sigma={\rm {\bf U}}_h^\sigma \; \overline{T}^\sigma |
---|
| 943 | \end{equation} |
---|
| 944 | with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, |
---|
| 945 | and $A_l^\sigma$ the lateral diffusivity in the BBL. Following \citet{BeckDos1998}, |
---|
| 946 | the latter is prescribed with a spatial dependence, $e.g.$ in the conditional form |
---|
| 947 | \begin{equation} \label{Eq_tra_bbl_Aadv} |
---|
| 948 | A_l^\sigma (i,j,t)=\left\{ {\begin{array}{l} |
---|
| 949 | A_{bbl} \quad \quad \mbox{if} \quad \nabla_\sigma \rho \cdot \nabla H<0 |
---|
| 950 | \quad \quad \mbox{and} \quad {\rm {\bf U}}_h \cdot \nabla H<0 \\ |
---|
| 951 | \\ |
---|
| 952 | 0\quad \quad \;\,\mbox{otherwise} \\ |
---|
| 953 | \end{array}} \right. |
---|
| 954 | \end{equation} |
---|
| 955 | |
---|
[707] | 956 | % ================================================================ |
---|
| 957 | % Tracer damping |
---|
| 958 | % ================================================================ |
---|
[817] | 959 | \section [Tracer damping (\textit{tradmp})] |
---|
| 960 | {Tracer damping (\mdl{tradmp})} |
---|
[707] | 961 | \label{TRA_dmp} |
---|
[1225] | 962 | %--------------------------------------------namtdp----------------------------------------------------- |
---|
| 963 | \namdisplay{namtdp} |
---|
[707] | 964 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 965 | |
---|
[817] | 966 | In some applications it can be useful to add a Newtonian damping term |
---|
| 967 | into the temperature and salinity equations: |
---|
[707] | 968 | \begin{equation} \label{Eq_tra_dmp} |
---|
| 969 | \begin{split} |
---|
| 970 | \frac{\partial T}{\partial t}=\;\cdots \;-\gamma \,\left( {T-T_o } \right) \\ |
---|
| 971 | \\ |
---|
| 972 | \frac{\partial S}{\partial t}=\;\cdots \;-\gamma \,\left( {S-S_o } \right) |
---|
| 973 | \end{split} |
---|
| 974 | \end{equation} |
---|
[817] | 975 | where $\gamma$ is the inverse of a time scale, and $T_o$ and $S_o$ |
---|
| 976 | are given temperature and salinity fields (usually a climatology). |
---|
| 977 | The restoring term is added when \key{tradmp} is defined. |
---|
| 978 | It also requires that both \key{temdta} and \key{saldta} are defined |
---|
| 979 | ($i.e.$ that $T_o$ and $S_o$ are read). The restoring coefficient |
---|
| 980 | $S_o$ is a three-dimensional array initialized by the user in routine |
---|
| 981 | \rou{dtacof} also located in module \mdl{tradmp}. |
---|
[707] | 982 | |
---|
[817] | 983 | The two main cases in which \eqref{Eq_tra_dmp} is used are \textit{(a)} |
---|
| 984 | the specification of the boundary conditions along artificial walls of a |
---|
| 985 | limited domain basin and \textit{(b)} the computation of the velocity |
---|
| 986 | field associated with a given $T$-$S$ field (for example to build the |
---|
| 987 | initial state of a prognostic simulation, or to use the resulting velocity |
---|
| 988 | field for a passive tracer study). The first case applies to regional |
---|
| 989 | models that have artificial walls instead of open boundaries. |
---|
| 990 | In the vicinity of these walls, $S_o$ takes large values (equivalent to |
---|
| 991 | a time scale of a few days) whereas it is zero in the interior of the |
---|
| 992 | model domain. The second case corresponds to the use of the robust |
---|
| 993 | diagnostic method \citep{Sarmiento1982}. It allows us to find the velocity |
---|
| 994 | field consistent with the model dynamics whilst having a $T$-$S$ field |
---|
| 995 | close to a given climatological field ($T_o -S_o$). The time scale |
---|
| 996 | associated with $S_o$ is generally not a constant but spatially varying |
---|
| 997 | in order to respect other properties. For example, it is usually set to zero |
---|
| 998 | in the mixed layer (defined either on a density or $S_o$ criterion) |
---|
| 999 | \citep{Madec1996} and in the equatorial region |
---|
| 1000 | \citep{Reverdin1991, Fujio1991, MartiTh1992} since these two regions |
---|
| 1001 | have a short time scale of adjustment; while smaller $S_o$ are used |
---|
| 1002 | in the deep ocean where the typical time scale is long \citep{Sarmiento1982}. |
---|
| 1003 | In addition the time scale is reduced (even to zero) along the western |
---|
| 1004 | boundary to allow the model to reconstruct its own western boundary |
---|
| 1005 | structure in equilibrium with its physics. The choice of a |
---|
| 1006 | Newtonian damping acting in the mixed layer or not is controlled by |
---|
| 1007 | namelist parameter \np{nmldmp}. |
---|
[707] | 1008 | |
---|
[817] | 1009 | The robust diagnostic method is very efficient in preventing temperature |
---|
| 1010 | drift in intermediate waters but it produces artificial sources of heat and salt |
---|
| 1011 | within the ocean. It also has undesirable effects on the ocean convection. |
---|
| 1012 | It tends to prevent deep convection and subsequent deep-water formation, |
---|
| 1013 | by stabilising the water column too much. |
---|
[707] | 1014 | |
---|
[817] | 1015 | An example of the computation of $S_o$ for robust diagnostic experiments |
---|
| 1016 | with the ORCA2 model is provided in the \mdl{tradmp} module |
---|
| 1017 | (subroutines \rou{dtacof} and \rou{cofdis} which compute the coefficient |
---|
| 1018 | and the distance to the bathymetry, respectively). These routines are |
---|
| 1019 | provided as examples and can be customised by the user. |
---|
[707] | 1020 | |
---|
| 1021 | % ================================================================ |
---|
| 1022 | % Tracer time evolution |
---|
| 1023 | % ================================================================ |
---|
[817] | 1024 | \section [Tracer time evolution (\textit{tranxt})] |
---|
| 1025 | {Tracer time evolution (\mdl{tranxt})} |
---|
[707] | 1026 | \label{TRA_nxt} |
---|
| 1027 | %--------------------------------------------namdom----------------------------------------------------- |
---|
| 1028 | \namdisplay{namdom} |
---|
| 1029 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 1030 | |
---|
[817] | 1031 | The general framework for tracer time stepping is a leap-frog scheme, |
---|
| 1032 | $i.e.$ a three level centred time scheme associated with a Asselin time |
---|
| 1033 | filter (cf. \S\ref{DOM_nxt}): |
---|
[707] | 1034 | \begin{equation} \label{Eq_tra_nxt} |
---|
| 1035 | \begin{split} |
---|
| 1036 | T^{t+\Delta t} &= T^{t-\Delta t} + 2 \, \Delta t \ \text{RHS}_T^t \\ |
---|
| 1037 | \\ |
---|
| 1038 | T_f^t \;\ \quad &= T^t \;\quad +\gamma \,\left[ {T_f^{t-\Delta t} -2T^t+T^{t+\Delta t}} \right] |
---|
| 1039 | \end{split} |
---|
| 1040 | \end{equation} |
---|
[817] | 1041 | where $\text{RHS}_T$ is the right hand side of the temperature equation, |
---|
| 1042 | the subscript $f$ denotes filtered values and $\gamma$ is the Asselin |
---|
| 1043 | coefficient. $\gamma$ is initialized as \np{atfp} (\textbf{namelist} parameter). |
---|
| 1044 | Its default value is \np{atfp=0.1}. |
---|
[707] | 1045 | |
---|
[817] | 1046 | When the vertical mixing is solved implicitly, the update of the \textit{next} tracer |
---|
| 1047 | fields is done in module \mdl{trazdf}. In this case only the swapping of arrays |
---|
| 1048 | and the Asselin filtering is done in the \mdl{tranxt} module. |
---|
[707] | 1049 | |
---|
[817] | 1050 | In order to prepare for the computation of the \textit{next} time step, |
---|
| 1051 | a swap of tracer arrays is performed: $T^{t-\Delta t} = T^t$ and $T^t = T_f$. |
---|
[707] | 1052 | |
---|
| 1053 | % ================================================================ |
---|
| 1054 | % Equation of State (eosbn2) |
---|
| 1055 | % ================================================================ |
---|
[817] | 1056 | \section [Equation of State (\textit{eosbn2}) ] |
---|
| 1057 | {Equation of State (\mdl{eosbn2}) } |
---|
[707] | 1058 | \label{TRA_eosbn2} |
---|
| 1059 | %--------------------------------------------nameos----------------------------------------------------- |
---|
| 1060 | \namdisplay{nameos} |
---|
| 1061 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 1062 | |
---|
| 1063 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 1064 | % Equation of State |
---|
| 1065 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 1066 | \subsection{Equation of State (\np{neos} = 0, 1 or 2)} |
---|
| 1067 | \label{TRA_eos} |
---|
| 1068 | |
---|
[817] | 1069 | It is necessary to know the equation of state for the ocean very accurately |
---|
| 1070 | to determine stability properties (especially the Brunt-Vais\"{a}l\"{a} frequency), |
---|
| 1071 | particularly in the deep ocean. The ocean density is a non linear empirical |
---|
| 1072 | function of \textit{in situ }temperature, salinity and pressure. The reference |
---|
| 1073 | equation of state is that defined by the Joint Panel on Oceanographic Tables |
---|
| 1074 | and Standards \citep{UNESCO1983}. It was the standard equation of state |
---|
| 1075 | used in early releases of OPA. However, even though this computation is |
---|
| 1076 | fully vectorised, it is quite time consuming ($15$ to $20${\%} of the total |
---|
| 1077 | CPU time) since it requires the prior computation of the \textit{in situ} |
---|
| 1078 | temperature from the model \textit{potential} temperature using the |
---|
| 1079 | \citep{Bryden1973} polynomial for adiabatic lapse rate and a $4^th$ order |
---|
| 1080 | Runge-Kutta integration scheme. Since OPA6, we have used the |
---|
| 1081 | \citet{JackMcD1995} equation of state for seawater instead. It allows the |
---|
| 1082 | computation of the \textit{in situ} ocean density directly as a function of |
---|
| 1083 | \textit{potential} temperature relative to the surface (an \NEMO variable), |
---|
| 1084 | the practical salinity (another \NEMO variable) and the pressure (assuming no |
---|
| 1085 | pressure variation along geopotential surfaces, i.e. the pressure in decibars is |
---|
[1224] | 1086 | approximated by the depth in meters). Both the \citet{UNESCO1983} and |
---|
| 1087 | \citet{JackMcD1995} equations of state have exactly the same except that |
---|
[817] | 1088 | the values of the various coefficients have been adjusted by \citet{JackMcD1995} |
---|
| 1089 | in order to directly use the \textit{potential} temperature instead of the |
---|
| 1090 | \textit{in situ} one. This reduces the CPU time of the in situ density computation |
---|
| 1091 | to about $3${\%} of the total CPU time, while maintaining a quite accurate |
---|
| 1092 | equation of state. |
---|
[707] | 1093 | |
---|
[817] | 1094 | In the computer code, a \textit{true} density $d$ is computed, $i.e.$ the ratio |
---|
| 1095 | of seawater volumic mass to $\rho_o$, a reference volumic mass (\textit{rau0} |
---|
| 1096 | defined in \mdl{phycst}, usually $rau0= 1,020~Kg/m^3$). The default option |
---|
| 1097 | (namelist prameter \np{neos}=0) is the \citet{JackMcD1995} equation of state. |
---|
| 1098 | Its use is highly recommended. However, for process studies, it is often |
---|
| 1099 | convenient to use a linear approximation of the density$^{\ast}$ |
---|
| 1100 | \footnote{$^{\ast }$ With the linear equation of state there is no longer |
---|
| 1101 | a distinction between \textit{in situ} and \textit{potential} density. Cabling |
---|
| 1102 | and thermobaric effects are also removed.}. |
---|
| 1103 | Two linear formulations are available: a function of $T$ only (\np{neos}=1) |
---|
| 1104 | and a function of both $T$ and $S$ (\np{neos}=2): |
---|
[707] | 1105 | \begin{equation} \label{Eq_tra_eos_linear} |
---|
| 1106 | \begin{aligned} |
---|
| 1107 | d(T) &= {\rho (T)} / {\rho _0 } &&= 1.028 - \alpha \;T \\ |
---|
| 1108 | d(T,S) &= {\rho (T,S)} &&= \ \ \ \beta \;S - \alpha \;T |
---|
| 1109 | \end{aligned} |
---|
| 1110 | \end{equation} |
---|
[817] | 1111 | where $\alpha$ and $\beta$ are the thermal and haline expansion |
---|
| 1112 | coefficients, and $\rho_o$, the reference volumic mass, $rau0$. |
---|
| 1113 | ($\alpha$ and $\beta$ can be modified through the \np{ralpha} and |
---|
| 1114 | \np{rbeta} namelist parameters). Note that when $d$ is a function |
---|
| 1115 | of $T$ only (\np{neos}=1), the salinity is a passive tracer and can be |
---|
| 1116 | used as such. |
---|
[707] | 1117 | |
---|
| 1118 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 1119 | % Brunt-Vais\"{a}l\"{a} Frequency |
---|
| 1120 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 1121 | \subsection{Brunt-Vais\"{a}l\"{a} Frequency (\np{neos} = 0, 1 or 2)} |
---|
| 1122 | \label{TRA_bn2} |
---|
| 1123 | |
---|
[817] | 1124 | An accurate computation of the ocean stability (i.e. of $N$, the brunt-Vais\"{a}l\"{a} |
---|
| 1125 | frequency) is of paramount importance as it is used in several ocean |
---|
| 1126 | parameterisations (namely TKE, KPP, Richardson number dependent |
---|
| 1127 | vertical diffusion, enhanced vertical diffusion, non-penetrative convection, |
---|
| 1128 | iso-neutral diffusion). In particular, one must be aware that $N^2$ has to |
---|
| 1129 | be computed with an \textit{in situ} reference. The expression for $N^2$ |
---|
| 1130 | depends on the type of equation of state used (\np{neos} namelist parameter). |
---|
[707] | 1131 | |
---|
| 1132 | For \np{neos}=0 (\citet{JackMcD1995} equation of state), the \citet{McDougall1987} |
---|
[817] | 1133 | polynomial expression is used (with the pressure in decibar approximated by |
---|
| 1134 | the depth in meters): |
---|
| 1135 | \begin{equation} \label{Eq_tra_bn2} |
---|
| 1136 | N^2 = \frac{g}{e_{3w}} \; \beta \ |
---|
| 1137 | \left( \alpha / \beta \ \delta_{k+1/2}[T] - \delta_{k+1/2}[S] \right) |
---|
| 1138 | \end{equation} |
---|
| 1139 | where $\alpha$ ($\beta$) is the thermal (haline) expansion coefficient. |
---|
| 1140 | They are a function of |
---|
| 1141 | $\overline{T}^{\,k+1/2},\widetilde{S}=\overline{S}^{\,k+1/2} - 35.$, |
---|
| 1142 | and $z_w$, with $T$ the \textit{potential} temperature and |
---|
| 1143 | $\widetilde{S}$ a salinity anomaly. |
---|
| 1144 | Note that both $\alpha$ and $\beta$ depend on \textit{potential} |
---|
| 1145 | temperature and salinity which are averaged at $w$-points prior |
---|
| 1146 | to the computation instead of being computed at $T$-points and |
---|
| 1147 | then averaged to $w$-points. |
---|
[707] | 1148 | |
---|
[817] | 1149 | When a linear equation of state is used (\np{neos}=1 or 2, |
---|
| 1150 | \eqref{Eq_tra_bn2} reduces to: |
---|
[707] | 1151 | \begin{equation} \label{Eq_tra_bn2_linear} |
---|
| 1152 | N^2 = \frac{g}{e_{3w}} \left( \beta \;\delta_{k+1/2}[S] - \alpha \;\delta_{k+1/2}[T] \right) |
---|
| 1153 | \end{equation} |
---|
[817] | 1154 | where $\alpha$ and $\beta $ are the constant coefficients used to |
---|
| 1155 | defined the linear equation of state \eqref{Eq_tra_eos_linear}. |
---|
[707] | 1156 | |
---|
| 1157 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 1158 | % Specific Heat |
---|
| 1159 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 1160 | \subsection [Specific Heat (\textit{phycst})] |
---|
| 1161 | {Specific Heat (\mdl{phycst})} |
---|
[707] | 1162 | \label{TRA_adv_ldf} |
---|
| 1163 | |
---|
[817] | 1164 | The specific heat of sea water, $C_p$, is a function of temperature, salinity |
---|
| 1165 | and pressure \citep{UNESCO1983}. It is only used in the model to convert |
---|
| 1166 | surface heat fluxes into surface temperature increase and so the pressure |
---|
| 1167 | dependence is neglected. The dependence on $T$ and $S$ is weak. |
---|
| 1168 | For example, with $S=35~psu$, $C_p$ increases from $3989$ to $4002$ |
---|
| 1169 | when $T$ varies from -2~\degres C to 31~\degres C. Therefore, $C_p$ has |
---|
| 1170 | been chosen as a constant: $C_p=4.10^3~J\,Kg^{-1}\,\degres K^{-1}$. |
---|
| 1171 | Its value is set in \mdl{phycst} module. |
---|
[707] | 1172 | |
---|
[817] | 1173 | %%% |
---|
| 1174 | \gmcomment{ STEVEN: consistency, no other computer variable names are |
---|
| 1175 | supplied, so why this one} |
---|
| 1176 | %%% |
---|
| 1177 | |
---|
[707] | 1178 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 1179 | % Freezing Point of Seawater |
---|
| 1180 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 1181 | \subsection [Freezing Point of Seawater (\textit{ocfzpt})] |
---|
| 1182 | {Freezing Point of Seawater (\mdl{ocfzpt})} |
---|
[707] | 1183 | \label{TRA_fzp} |
---|
| 1184 | |
---|
| 1185 | The freezing point of seawater is a function of salinity and pressure \citep{UNESCO1983}: |
---|
| 1186 | \begin{equation} \label{Eq_tra_eos_fzp} |
---|
| 1187 | \begin{split} |
---|
[1224] | 1188 | T_f (S,p) = \left( -0.0575 + 1.710523 \;10^{-3} \, \sqrt{S} |
---|
[707] | 1189 | - 2.154996 \;10^{-4} \,S \right) \ S \\ |
---|
[1224] | 1190 | - 7.53\,10^{-3} \ \ p |
---|
[707] | 1191 | \end{split} |
---|
| 1192 | \end{equation} |
---|
| 1193 | |
---|
[817] | 1194 | \eqref{Eq_tra_eos_fzp} is only used to compute the potential freezing point of |
---|
| 1195 | sea water ($i.e.$ referenced to the surface $p=0$), thus the pressure dependent |
---|
| 1196 | terms in \eqref{Eq_tra_eos_fzp} (last term) have been dropped. The \textit{before} |
---|
| 1197 | and \textit{now} surface freezing point is introduced in the code as $fzptb$ and |
---|
| 1198 | $fzptn$ 2D arrays together with a \textit{now} mask (\textit{freezn}) which takes |
---|
| 1199 | the value 0 or 1 depending on whether the ocean temperature is above or at the |
---|
| 1200 | freezing point. Caution: do not confuse \textit{freezn} with the fraction of lead |
---|
| 1201 | (\textit{frld}) defined in LIM. |
---|
[707] | 1202 | |
---|
[817] | 1203 | %%% |
---|
| 1204 | \gmcomment{STEVEN: consistency, not many computer variable names are supplied, so why these ===> gm I agree this should evolve both here and in the code itself} |
---|
| 1205 | %%% |
---|
| 1206 | |
---|
[707] | 1207 | % ================================================================ |
---|
| 1208 | % Horizontal Derivative in zps-coordinate |
---|
| 1209 | % ================================================================ |
---|
[817] | 1210 | \section [Horizontal Derivative in \textit{zps}-coordinate (\textit{zpshde})] |
---|
| 1211 | {Horizontal Derivative in \textit{zps}-coordinate (\mdl{zpshde})} |
---|
[707] | 1212 | \label{TRA_zpshde} |
---|
| 1213 | |
---|
[817] | 1214 | \gmcomment{STEVEN: to be consistent with earlier discussion of differencing and averaging operators, I've changed "derivative" to "difference" and "mean" to "average"} |
---|
| 1215 | |
---|
| 1216 | With partial bottom cells (\np{ln\_zps}=.true.), in general, tracers in horizontally |
---|
| 1217 | adjacent cells live at different depths. Horizontal gradients of tracers are needed |
---|
| 1218 | for horizontal diffusion (\mdl{traldf} module) and for the hydrostatic pressure |
---|
| 1219 | gradient (\mdl{dynhpg} module) to be active. |
---|
| 1220 | \gmcomment{STEVEN from gm : question: not sure of what -to be active- means} |
---|
| 1221 | Before taking horizontal gradients between the tracers next to the bottom, a linear |
---|
| 1222 | interpolation in the vertical is used to approximate the deeper tracer as if it actually |
---|
| 1223 | lived at the depth of the shallower tracer point (Fig.~\ref{Fig_Partial_step_scheme}). |
---|
| 1224 | For example, for temperature in the $i$-direction the needed interpolated |
---|
| 1225 | temperature, $\widetilde{T}$, is: |
---|
| 1226 | |
---|
[707] | 1227 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 1228 | \begin{figure}[!p] \label{Fig_Partial_step_scheme} \begin{center} |
---|
[998] | 1229 | \includegraphics[width=0.9\textwidth]{./TexFiles/Figures/Partial_step_scheme.pdf} |
---|
[817] | 1230 | \caption{ Discretisation of the horizontal difference and average of tracers in the $z$-partial step coordinate (\np{ln\_zps}=.true.) in the case $( e3w_k^{i+1} - e3w_k^i )>0$. A linear interpolation is used to estimate $\widetilde{T}_k^{i+1}$, the tracer value at the depth of the shallower tracer point of the two adjacent bottom $T$-points. The horizontal difference is then given by: $\delta _{i+1/2} T_k= \widetilde{T}_k^{\,i+1} -T_k^{\,i}$ and the average by: $\overline{T}_k^{\,i+1/2}= ( \widetilde{T}_k^{\,i+1/2} - T_k^{\,i} ) / 2$. } |
---|
[707] | 1231 | \end{center} \end{figure} |
---|
| 1232 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 1233 | \begin{equation*} |
---|
| 1234 | \widetilde{T}= \left\{ \begin{aligned} |
---|
| 1235 | &T^{\,i+1} -\frac{ \left( e_{3w}^{i+1} -e_{3w}^i \right)}{ e_{3w}^{i+1} }\;\delta _k T^{i+1} |
---|
| 1236 | && \quad\text{if $\ e_{3w}^{i+1} \geq e_{3w}^i$ } \\ |
---|
| 1237 | \\ |
---|
| 1238 | &T^{\,i} \ \ \ \,+\frac{ \left( e_{3w}^{i+1} -e_{3w}^i \right) }{e_{3w}^i }\;\delta _k T^{i+1} |
---|
| 1239 | && \quad\text{if $\ e_{3w}^{i+1} < e_{3w}^i$ } |
---|
| 1240 | \end{aligned} \right. |
---|
| 1241 | \end{equation*} |
---|
[817] | 1242 | and the resulting forms for the horizontal difference and the horizontal average |
---|
| 1243 | value of $T$ at a $U$-point are: |
---|
[707] | 1244 | \begin{equation} \label{Eq_zps_hde} |
---|
| 1245 | \begin{aligned} |
---|
| 1246 | \delta _{i+1/2} T= \begin{cases} |
---|
| 1247 | \ \ \ \widetilde {T}\quad\ -T^i & \ \ \quad\quad\text{if $\ e_{3w}^{i+1} \geq e_{3w}^i$ } \\ |
---|
| 1248 | \\ |
---|
| 1249 | \ \ \ T^{\,i+1}-\widetilde{T} & \ \ \quad\quad\text{if $\ e_{3w}^{i+1} < e_{3w}^i$ } |
---|
| 1250 | \end{cases} \\ |
---|
| 1251 | \\ |
---|
| 1252 | \overline {T}^{\,i+1/2} \ = \begin{cases} |
---|
| 1253 | ( \widetilde {T}\ \ \;\,-T^{\,i}) / 2 & \;\ \ \quad\text{if $\ e_{3w}^{i+1} \geq e_{3w}^i$ } \\ |
---|
| 1254 | \\ |
---|
| 1255 | ( T^{\,i+1}-\widetilde{T} ) / 2 & \;\ \ \quad\text{if $\ e_{3w}^{i+1} < e_{3w}^i$ } |
---|
| 1256 | \end{cases} |
---|
| 1257 | \end{aligned} |
---|
| 1258 | \end{equation} |
---|
| 1259 | |
---|
[817] | 1260 | The computation of horizontal derivative of tracers as well as of density is |
---|
| 1261 | performed once for all at each time step in \mdl{zpshde} module and stored |
---|
| 1262 | in shared arrays to be used when needed. It has to be emphasized that the |
---|
| 1263 | procedure used to compute the interpolated density, $\widetilde{\rho}$, is not |
---|
| 1264 | the same as that used for $T$ and $S$. Instead of forming a linear approximation |
---|
| 1265 | of density, we compute $\widetilde{\rho }$ from the interpolated values of $T$ |
---|
| 1266 | and $S$, and the pressure at a $u$-point (in the equation of state pressure is |
---|
| 1267 | approximated by depth, see \S\ref{TRA_eos} ) : |
---|
[707] | 1268 | \begin{equation} \label{Eq_zps_hde_rho} |
---|
| 1269 | \widetilde{\rho } = \rho ( {\widetilde{T},\widetilde {S},z_u }) |
---|
| 1270 | \quad \text{where }\ z_u = \min \left( {z_T^{i+1} ,z_T^i } \right) |
---|
| 1271 | \end{equation} |
---|
| 1272 | |
---|
[817] | 1273 | This is a much better approximation as the variation of $\rho$ with depth (and |
---|
| 1274 | thus pressure) is highly non-linear with a true equation of state and thus is badly |
---|
| 1275 | approximated with a linear interpolation. This approximation is used to compute |
---|
| 1276 | both the horizontal pressure gradient (\S\ref{DYN_hpg}) and the slopes of neutral |
---|
| 1277 | surfaces (\S\ref{LDF_slp}) |
---|
[707] | 1278 | |
---|
[817] | 1279 | Note that in almost all the advection schemes presented in this Chapter, both |
---|
| 1280 | averaging and differencing operators appear. Yet \eqref{Eq_zps_hde} has not |
---|
| 1281 | been used in these schemes: in contrast to diffusion and pressure gradient |
---|
| 1282 | computations, no correction for partial steps is applied for advection. The main |
---|
| 1283 | motivation is to preserve the domain averaged mean variance of the advected |
---|
| 1284 | field when using the $2^{nd}$ order centred scheme. Sensitivity of the advection |
---|
| 1285 | schemes to the way horizontal averages are performed in the vicinity of partial |
---|
| 1286 | cells should be further investigated in the near future. |
---|
| 1287 | %%% |
---|
| 1288 | \gmcomment{gm : this last remark has to be done} |
---|
| 1289 | %%% |
---|