# Changeset 1831 for branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_DOM.tex

Ignore:
Timestamp:
2010-04-12T16:59:59+02:00 (11 years ago)
Message:

cover, namelist, rigid-lid, e3t, appendices, see ticket: #658

File:
1 edited

### Legend:

Unmodified
 r1224 as well as the DOM (DOMain) directory. $\$\newline    % force a new ligne % ================================================================ % Fundamentals of the Discretisation \begin{figure}[!tb] \label{Fig_cell}  \begin{center} \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_cell.pdf} \caption{Arrangement of variables. $T$ indicates scalar points where temperature, \caption{Arrangement of variables. $t$ indicates scalar points where temperature, salinity, density, pressure and horizontal divergence are defined. ($u$,$v$,$w$) indicates vector points, and $f$ indicates vorticity points where both relative and Special attention has been given to the homogeneity of the solution in the three space directions. The arrangement of variables is the same in all directions. It consists of cells centred on scalar points ($T$, $S$, $p$, $\rho$) with vector It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in the centre of each face of the cells (Fig. \ref{Fig_cell}). This is the generalisation to three dimensions of the well-known C'' grid in Similar operators are defined with respect to $i+1/2$, $j$, $j+1/2$, $k$, and $k+1/2$. Following \eqref{Eq_PE_grad} and \eqref{Eq_PE_lap}, the gradient of a variable $q$ defined at a $T$-point has its three components defined at $u$-, $v$- and $w$-points while its Laplacien is defined at $T$-point. These operators have variable $q$ defined at a $t$-point has its three components defined at $u$-, $v$- and $w$-points while its Laplacien is defined at $t$-point. These operators have the following discrete forms in the curvilinear $s$-coordinate system: \label{Eq_DOM_grad} \nabla q\equiv    \frac{1}{e_{1u} }\delta _{i+1/2} \left[ q \right]\;\,{\rm {\bf i}} +  \frac{1}{e_{2v} }\delta _{j+1/2} \left[ q \right]\;\,{\rm {\bf j}} +  \frac{1}{e_{3w} }\delta _{k+1/2} \left[ q \right]\;\,{\rm {\bf k}} \nabla q\equiv    \frac{1}{e_{1u} } \delta _{i+1/2 } [q] \;\,{\rm {\bf i}} +  \frac{1}{e_{2v} } \delta _{j+1/2 } [q] \;\,{\rm {\bf j}} +  \frac{1}{e_{3w}} \delta _{k+1/2} [q] \;\,{\rm {\bf k}} \begin{multline} \label{Eq_DOM_lap} \Delta q\equiv \frac{1}{e_{1T} {\kern 1pt}e_{2T} {\kern 1pt}e_{3T} }\;\left( {\delta _i \left[ {\frac{e_{2u} e_{3u} }{e_{1u} }\;\delta _{i+1/2} \left[ q \right]} \right] +\delta _j \left[ {\frac{e_{1v} e_{3v} }{e_{2v} }\;\delta _{j+1/2} \left[ q \right]} \right]\;} \right)     \\ +\frac{1}{e_{3T} }\delta _k \left[ {\frac{1}{e_{3w} }\;\delta _{k+1/2} \left[ q \right]} \right] \Delta q\equiv \frac{1}{e_{1t}\,e_{2t}\,e_{3t} } \;\left(          \delta_i  \left[ \frac{e_{2u}\,e_{3u}} {e_{1u}} \;\delta_{i+1/2} [q] \right] +                        \delta_j  \left[ \frac{e_{1v}\,e_{3v}}  {e_{2v}} \;\delta_{j+1/2} [q] \right] \;  \right)      \\ +\frac{1}{e_{3t}} \delta_k \left[ \frac{1}{e_{3w} }                     \;\delta_{k+1/2} [q] \right] \end{multline} Following \eqref{Eq_PE_curl} and \eqref{Eq_PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$ defined at vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and its divergence defined at $T$-points: \label{Eq_DOM_curl} \begin{split} \nabla \times {\rm {\bf A}}\equiv \frac{1}{e_{2v} {\kern 1pt}e_{3vw} }{\kern 1pt}\,\;\left( {\delta _{j+1/2} \left[ {e_{3w} a_3 } \right]-\delta _{k+1/2} \left[ {e_{2v} a_2 } \right]} \right)  &\;\;{\rm {\bf i}} \\ +\frac{1}{e_{2u} {\kern 1pt}e_{3uw} }\;\left( {\delta _{k+1/2} \left[ {e_{1u} a_1 } \right]-\delta _{i+1/2} \left[ {e_{3w} a_3 } \right]} \right)  &\;\;{\rm{\bf j}} \\ +\frac{e_{3f} }{e_{1f} {\kern 1pt}e_{2f} }\,{\kern 1pt}\;\left( {\delta _{i+1/2} \left[ {e_{2v} a_2 } \right]-\delta _{j+12} \left[ {e_{1u} a_1 } \right]} \right)  &\;\;{\rm {\bf k}} \end{split} Following \eqref{Eq_PE_curl} and \eqref{Eq_PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$ defined at vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and its divergence defined at $t$-points: \begin{eqnarray}  \label{Eq_DOM_curl} \nabla \times {\rm {\bf A}}\equiv & \frac{1}{e_{2v}\,e_{3vw} } \ \left( \delta_{j +1/2} \left[e_{3w}\,a_3 \right] -\delta_{k+1/2} \left[e_{2v} \,a_2 \right] \right)  &\ \rm{\bf i} \\ +& \frac{1}{e_{2u}\,e_{3uw}} \ \left( \delta_{k+1/2} \left[e_{1u}\,a_1  \right] -\delta_{i +1/2} \left[e_{3w}\,a_3 \right] \right)  &\ \rm{\bf j} \\ +& \frac{1}{e_{1f} \,e_{2f}    } \ \left( \delta_{i +1/2} \left[e_{2v}\,a_2  \right] -\delta_{j +1/2} \left[e_{1u}\,a_1 \right] \right)  &\ \rm{\bf k} \end{eqnarray} \label{Eq_DOM_div} \nabla \cdot {\rm {\bf A}}=\frac{1}{e_{1T} e_{2T} e_{3T} }\left( {\delta _i \left[ {e_{2u} e_{3u} a_1 } \right]+\delta _j \left[ {e_{1v} e_{3v} a_2 } \right]} \right)+\frac{1}{e_{3T} }\delta _k \left[ {a_3 } \right] \nabla \cdot \rm{\bf A}=\frac{1}{e_{1t}\,e_{2t}\,e_{3t}}\left( \delta_i \left[e_{2u}\,e_{3u}\,a_1 \right] +\delta_j \left[e_{1v}\,e_{3v}\,a_2 \right] \right)+\frac{1}{e_{3t} }\delta_k \left[a_3 \right] horizontal location of a grid point. For example \eqref{Eq_DOM_div} reduces to: \begin{equation*} \nabla \cdot {\rm {\bf A}}=\frac{1}{e_{1T} e_{2T} }\left( {\delta _i \left[ {e_{2u} a_1 } \right]+\delta _j \left[ {e_{1v}  a_2 } \right]} \right)+\frac{1}{e_{3T} }\delta _k \left[ {a_3 } \right] \nabla \cdot \rm{\bf A}=\frac{1}{e_{1t}\,e_{2t}} \left( \delta_i \left[e_{2u}\,a_1 \right] +\delta_j \left[e_{1v}\, a_2 \right]  \right) +\frac{1}{e_{3t}} \delta_k \left[             a_3 \right] \end{equation*} for a quantity $q$ which is a masked field (i.e. equal to zero inside solid area): \label{DOM_bar} \bar q   = \frac{1}{H}\int_{k^b}^{k^o} {q\;e_{3q} \,dk} \bar q   =         \frac{1}{H}    \int_{k^b}^{k^o} {q\;e_{3q} \,dk} \equiv \frac{1}{H_q }\sum\limits_k {q\;e_{3q} } It is straightforward to demonstrate that these properties are verified locally in discrete form as soon as the scalar $q$ is taken at $T$-points and the vector discrete form as soon as the scalar $q$ is taken at $t$-points and the vector \textbf{A} has its components defined at vector points $(u,v,w)$. The array representation used in the \textsc{Fortran} code requires an integer indexing while the analytical definition of the mesh (see \S\ref{DOM_cell}) is associated with the use of integer values for $T$-points and both integer and associated with the use of integer values for $t$-points and both integer and integer and a half values for all the other points. Therefore a specific integer indexing must be defined for points other than $T$-points ($i.e.$ velocity and indexing must be defined for points other than $t$-points ($i.e.$ velocity and vorticity grid-points). Furthermore, the direction of the vertical indexing has been changed so that the surface level is at $k=1$. \label{DOM_Num_Index_hor} The indexing in the horizontal plane has been chosen as shown in Fig.\ref{Fig_index_hor}. For an increasing $i$ index ($j$ index), the $T$-point and the eastward $u$-point (northward $v$-point) have the same index (see the dashed area in Fig.\ref{Fig_index_hor}). A $T$-point and its nearest northeast $f$-point have the same $i$-and $j$-indices. The indexing in the horizontal plane has been chosen as shown in Fig.\ref{Fig_index_hor}. For an increasing $i$ index ($j$ index), the $t$-point and the eastward $u$-point (northward $v$-point) have the same index (see the dashed area in Fig.\ref{Fig_index_hor}). A $t$-point and its nearest northeast $f$-point have the same $i$-and $j$-indices. % ----------------------------------- to the indexing used in the semi-discrete equations and given in \S\ref{DOM_cell}. The sea surface corresponds to the $w$-level $k=1$ which is the same index as $T$-level just below (Fig.\ref{Fig_index_vert}). The last $w$-level ($k=jpk$) as $t$-level just below (Fig.\ref{Fig_index_vert}). The last $w$-level ($k=jpk$) either corresponds to the ocean floor or is inside the bathymetry while the last $T$-level is always inside the bathymetry (Fig.\ref{Fig_index_vert}). Note that for an increasing $k$ index, a $w$-point and the $T$-point just below have the $t$-level is always inside the bathymetry (Fig.\ref{Fig_index_vert}). Note that for an increasing $k$ index, a $w$-point and the $t$-point just below have the same $k$ index, in opposition to what is done in the horizontal plane where it is the $T$-point and the nearest velocity points in the direction of the horizontal it is the $t$-point and the nearest velocity points in the direction of the horizontal axis that have the same $i$ or $j$ index (compare the dashed area in Fig.\ref{Fig_index_hor} and \ref{Fig_index_vert}). Since the scale factors are \S\ref{LBC_mpp}). $\$\newline    % force a new ligne % ================================================================ % Domain: Horizontal Grid (mesh) factors in the horizontal plane as follows: \begin{flalign*} \lambda_T &\equiv \text{glamt}= \lambda(i)     & \varphi_T &\equiv \text{gphit} = \varphi(j)\\ \lambda_t &\equiv \text{glamt}= \lambda(i)     & \varphi_t &\equiv \text{gphit} = \varphi(j)\\ \lambda_u &\equiv \text{glamu}= \lambda(i+1/2)& \varphi_u &\equiv \text{gphiu}= \varphi(j)\\ \lambda_v &\equiv \text{glamv}= \lambda(i)       & \varphi_v &\equiv \text{gphiv} = \varphi(j+1/2)\\ \end{flalign*} \begin{flalign*} e_{1T} &\equiv \text{e1t} = r_a |\lambda'(i)    \; \cos\varphi(j)  |& e_{2T} &\equiv \text{e2t} = r_a |\varphi'(j)|  \\ e_{1t} &\equiv \text{e1t} = r_a |\lambda'(i)    \; \cos\varphi(j)  |& e_{2t} &\equiv \text{e2t} = r_a |\varphi'(j)|  \\ e_{1u} &\equiv \text{e1t} = r_a |\lambda'(i+1/2)   \; \cos\varphi(j)  |& e_{2u} &\equiv \text{e2t} = r_a |\varphi'(j)|\\ considered and $r_a$ is the earth radius (defined in \mdl{phycst} along with all universal constants). Note that the horizontal position of and scale factors at $w$-points are exactly equal to those of $T$-points, thus no specific arrays at $w$-points are exactly equal to those of $t$-points, thus no specific arrays are defined at $w$-points. Note that the definition of the scale factors ($i.e.$ as the analytical first derivative of the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) is specific to the \NEMO model \citep{Marti1992}. As an example, $e_{1T}$ is defined locally at a $T$-point, whereas many other models on a C grid choose to define specific to the \NEMO model \citep{Marti_al_JGR92}. As an example, $e_{1t}$ is defined locally at a $t$-point, whereas many other models on a C grid choose to define such a scale factor as the distance between the $U$-points on each side of the $T$-point. Relying on an analytical transformation has two advantages: firstly, there $t$-point. Relying on an analytical transformation has two advantages: firstly, there is no ambiguity in the scale factors appearing in the discrete equations, since they are first introduced in the continuous equations; secondly, analytical transformations in the vertical, and (b) analytically derived grid-point position and scale factors. For both grids here,  the same $w$-point depth has been chosen but in (a) the $T$-points are set half way between $w$-points while in (b) they are defined from $t$-points are set half way between $w$-points while in (b) they are defined from an analytical function: $z(k)=5\,(i-1/2)^3 - 45\,(i-1/2)^2 + 140\,(i-1/2) - 150$. Note the resulting difference between the value of the grid-size $\Delta_k$ and and \pp{ppgphi0}). Note that for the Mercator grid the user need only provide an approximate starting latitude: the real latitude will be recalculated analytically, in order to ensure that the equator corresponds to line passing through $T$- in order to ensure that the equator corresponds to line passing through $t$- and $u$-points. and given by the parameters \pp{ppe1\_m} and \pp{ppe2\_m} respectively. The zonal grid coordinate (\textit{glam} arrays) is in kilometers, starting at zero with the first $T$-point. The meridional coordinate (gphi. arrays) is in kilometers, and the second $T$-point corresponds to coordinate $gphit=0$. The input with the first $t$-point. The meridional coordinate (gphi. arrays) is in kilometers, and the second $t$-point corresponds to coordinate $gphit=0$. The input parameter \pp{ppglam0} is ignored. \pp{ppgphi0} is used to set the reference latitude for computation of the Coriolis parameter. In the case of the beta plane, the output grid written when $\np{nmsh} \not=0$ is no more equal to the input grid. $\$\newline    % force a new ligne % ================================================================ % Domain: Vertical Grid (domzgr) \label{DOM_zgr} %-----------------------------------------nam_zgr & namdom------------------------------------------- \namdisplay{nam_zgr} \namdisplay{namzgr} \namdisplay{namdom} %------------------------------------------------------------------------------------------------------------- The reference coordinate transformation $z_0 (k)$ defines the arrays $gdept_0$ and $gdepw_0$ for $T$- and $w$-points, respectively. As indicated on and $gdepw_0$ for $t$- and $w$-points, respectively. As indicated on Fig.\ref{Fig_index_vert} \jp{jpk} is the number of $w$-levels. $gdepw_0(1)$ is the ocean surface. There are at most \jp{jpk}-1 $T$-points inside the ocean, the additional $T$-point at $jk=jpk$ is below the sea floor and is not used. The vertical location of $w$- and $T$-levels is defined from the analytic expression ocean surface. There are at most \jp{jpk}-1 $t$-points inside the ocean, the additional $t$-point at $jk=jpk$ is below the sea floor and is not used. The vertical location of $w$- and $t$-levels is defined from the analytic expression of the depth $z_0(k)$ whose analytical derivative with respect to $k$ provides the vertical scale factors. The user must provide the analytical expression of both \begin{center} \begin{tabular}{c||r|r|r|r} \hline \textbf{LEVEL}& \textbf{GDEPT}& \textbf{GDEPW}& \textbf{E3T }& \textbf{E3W  } \\ \hline \textbf{LEVEL}& \textbf{gdept}& \textbf{gdepw}& \textbf{e3t }& \textbf{e3w  } \\ \hline 1  &  \textbf{  5.00}   &       0.00 & \textbf{ 10.00} &  10.00 \\   \hline 2  &  \textbf{15.00} &    10.00 &   \textbf{ 10.00} &  10.00 \\   \hline Two variables in the namdom namelist are used to define the partial step vertical grid. The mimimum water thickness (in meters) allowed for a cell partially filled with bathymetry at level jk is the minimum of \np{e3zpsmin} (thickness in meters, usually $20~m$) or $e_{3t}(jk)*\np{e3zps\_rat}$ (a fraction, partially filled with bathymetry at level jk is the minimum of \np{rn\_e3zps\_min} (thickness in meters, usually $20~m$) or $e_{3t}(jk)*\np{rn\_e3zps\_rat}$ (a fraction, usually 10\%, of the default thickness $e_{3t}(jk)$). % ------------------------------------------------------------------------------------------------------------- \subsection   [$s$-coordinate (\np{ln\_sco})] {$s$-coordinate (\np{ln\_sco}=true)} {$s$-coordinate (\np{ln\_sco}=true)} \label{DOM_sco} %------------------------------------------nam_zgr_sco--------------------------------------------------- \namdisplay{nam_zgr_sco} \namdisplay{namzgr_sco} %-------------------------------------------------------------------------------------------------------------- In $s$-coordinate (\key{sco} is defined), the depth and thickness of the model \end{split} where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $T$-point where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point location in the horizontal and $z_0(k)$ is a function which varies from $0$ at the sea surface to $1$ at the ocean bottom. The depth field $h$ is not necessary the ocean sharp bathymetric gradients. A new flexible stretching function, modified from \citet{Song1994} is provided as an example: A new flexible stretching function, modified from \citet{Song_Haidvogel_JCP94} is provided as an example: \label{DOM_sco_function} \begin{split} Whatever the vertical coordinate used, the model offers the possibility of representing the bottom topography with steps that follow the face of the model cells (step like topography) \citep{Madec1996}. The distribution of model cells (step like topography) \citep{Madec_al_JPO96}. The distribution of the steps in the horizontal is defined in a 2D integer array, mbathy, which gives the number of ocean levels ($i.e.$ those that are not masked) at each $T$-point. mbathy is computed from the meter bathymetry using the definiton of gdept as the number of $T$-points which gdept $\leq$ bathy. Note that in version $t$-point. mbathy is computed from the meter bathymetry using the definiton of gdept as the number of $t$-points which gdept $\leq$ bathy. Note that in version NEMO v2.3, the user still has to provide the "level" bathymetry in a NetCDF file when using the full step option (\np{ln\_zco}), rather than the bathymetry domain (\np{ln\_dynspg\_rl}=.true. and \key{island} is defined), the \textit{mbathy} array must be provided and takes values from $-N$ to \jp{jpk}-1. It provides the following information: $mbathy(i,j) = -n, \ n \in \left] 0,N \right]$, $T$-points are land points on the $n^{th}$ island ; $mbathy(i,j) =0$, $T$-points are land points on the main land (continent) ; $mbathy(i,j) =k$, the first $k$ $T$- and $w$-points following information: $mbathy(i,j) = -n, \ n \in \left] 0,N \right]$, $t$-points are land points on the $n^{th}$ island ; $mbathy(i,j) =0$, $t$-points are land points on the main land (continent) ; $mbathy(i,j) =k$, the first $k$ $t$- and $w$-points are ocean points, the others are points below the ocean floor. %%% $\$\newline    % force a new ligne % ================================================================ % Time Discretisation x^{t+\Delta t} = x^{t-\Delta t} + 2 \, \Delta t \  \text{RHS}_x^{t-\Delta t,t,t+\Delta t} where $x$ stands for $u$, $v$, $T$ or $S$; RHS is the Right-Hand-Side of the where $x$ stands for $u$, $v$, $t$ or $S$; RHS is the Right-Hand-Side of the corresponding time evolution equation; $\Delta t$ is the time step; and the superscripts indicate the time at which a quantity is evaluated. Each term of the This is diffusive in time and conditionally stable. For example, the conditions for stability of second and fourth order horizontal diffusion schemes are \citep{Griffies2004}: conditions for stability of second and fourth order horizontal diffusion schemes are \citep{Griffies_Bk04}: \label{Eq_DOM_nxt_euler_stability} A^h < \left\{ approximation of the temperature equation is: \label{Eq_DOM_nxt_imp_zdf} \frac{T(k)^{t+1}-T(k)^{t-1}}{2\;\Delta t}\equiv \text{RHS}+\frac{1}{e_{3T} }\delta \frac{T(k)^{t+1}-T(k)^{t-1}}{2\;\Delta t}\equiv \text{RHS}+\frac{1}{e_{3t} }\delta _k \left[ {\frac{A_w^{vT} }{e_{3w} }\delta _{k+1/2} \left[ {T^{t+1}} \right]} \right] }        %% end add Implicit time stepping in case of variable volume thickness. Tracer case (NB for momentum in vector invariant form take care!) \begin{flalign*} &\frac{\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1}}{2\Delta t} \equiv \text{RHS}+ \delta _k \left[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta _{k+1/2} \left[ {T^{t+1}} \right]} \right]      \\ &\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1} \equiv {2\Delta t} \ \text{RHS}+ {2\Delta t} \ \delta _k \left[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta _{k+1/2} \left[ {T^{t+1}} \right]} \right]      \\ &\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1} \equiv 2\Delta t \ \text{RHS} + 2\Delta t \ \left\{ \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2} [ T_{k+1}^{t+1} - T_k      ^{t+1} ] - \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} [ T_k       ^{t+1} - T_{k-1}^{t+1} ]  \right\}     \\ &\\ &\left( e_{3t}\,T \right)_k^{t+1} -  {2\Delta t} \           \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}                  T_{k+1}^{t+1} + {2\Delta t} \ \left\{  \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2} +  \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2}     \right\}   T_{k    }^{t+1} -  {2\Delta t} \           \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2}                  T_{k-1}^{t+1}      \\ &\equiv \left( e_{3t}\,T \right)_k^{t-1} + {2\Delta t} \ \text{RHS}    \\ % \end{flalign*} \begin{flalign*} \allowdisplaybreaks \intertext{ Tracer case } % &  \qquad \qquad  \quad   -  {2\Delta t}                  \ \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2} \qquad \qquad \qquad  \qquad  T_{k+1}^{t+1}   \\ &+ {2\Delta t} \ \biggl\{  (e_{3t})_{k   }^{t+1}  \bigg. +    \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2} +   \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} \bigg. \biggr\}  \ \ \ T_{k   }^{t+1}  &&\\ & \qquad \qquad  \qquad \qquad \qquad \quad \ \ -  {2\Delta t} \                          \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2}                          \quad \ \ T_{k-1}^{t+1} \ \equiv \ \left( e_{3t}\,T \right)_k^{t-1} + {2\Delta t} \ \text{RHS}  \\ % \end{flalign*} \begin{flalign*} \allowdisplaybreaks \intertext{ Tracer content case } % & -  {2\Delta t} \              & \frac{(A_w^{vt})_{k+1/2}} {(e_{3w})_{k+1/2}^{t+1}\;(e_{3t})_{k+1}^{t+1}}  && \  \left( e_{3t}\,T \right)_{k+1}^{t+1}   &\\ & + {2\Delta t} \ \left[ 1  \right.+ & \frac{(A_w^{vt})_{k+1/2}} {(e_{3w})_{k+1/2}^{t+1}\;(e_{3t})_k^{t+1}} + & \frac{(A_w^{vt})_{k -1/2}} {(e_{3w})_{k -1/2}^{t+1}\;(e_{3t})_k^{t+1}}  \left.  \right]  & \left( e_{3t}\,T \right)_{k   }^{t+1}  &\\ & -  {2\Delta t} \               & \frac{(A_w^{vt})_{k -1/2}} {(e_{3w})_{k -1/2}^{t+1}\;(e_{3t})_{k-1}^{t+1}}     &\  \left( e_{3t}\,T \right)_{k-1}^{t+1} \equiv \left( e_{3t}\,T \right)_k^{t-1} + {2\Delta t} \ \text{RHS}  & \end{flalign*}