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

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_DOM.tex

    r1224 r1831  
    3131as well as the DOM (DOMain) directory.  
    3232 
     33$\ $\newline    % force a new ligne 
     34 
    3335% ================================================================ 
    3436% Fundamentals of the Discretisation 
     
    4648\begin{figure}[!tb] \label{Fig_cell}  \begin{center} 
    4749\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_cell.pdf} 
    48 \caption{Arrangement of variables. $T$ indicates scalar points where temperature,  
     50\caption{Arrangement of variables. $t$ indicates scalar points where temperature,  
    4951salinity, density, pressure and horizontal divergence are defined. ($u$,$v$,$w$)  
    5052indicates vector points, and $f$ indicates vorticity points where both relative and  
     
    5759Special attention has been given to the homogeneity of the solution in the three  
    5860space directions. The arrangement of variables is the same in all directions.  
    59 It consists of cells centred on scalar points ($T$, $S$, $p$, $\rho$) with vector  
     61It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector  
    6062points $(u, v, w)$ defined in the centre of each face of the cells (Fig. \ref{Fig_cell}).  
    6163This is the generalisation to three dimensions of the well-known ``C'' grid in  
     
    120122Similar operators are defined with respect to $i+1/2$, $j$, $j+1/2$, $k$, and  
    121123$k+1/2$. Following \eqref{Eq_PE_grad} and \eqref{Eq_PE_lap}, the gradient of a  
    122 variable $q$ defined at a $T$-point has its three components defined at $u$-, $v$-  
    123 and $w$-points while its Laplacien is defined at $T$-point. These operators have  
     124variable $q$ defined at a $t$-point has its three components defined at $u$-, $v$-  
     125and $w$-points while its Laplacien is defined at $t$-point. These operators have  
    124126the following discrete forms in the curvilinear $s$-coordinate system: 
    125127\begin{equation} \label{Eq_DOM_grad} 
    126 \nabla q\equiv    \frac{1}{e_{1u} }\delta _{i+1/2} \left[ q \right]\;\,{\rm {\bf i}} 
    127          +  \frac{1}{e_{2v} }\delta _{j+1/2} \left[ q \right]\;\,{\rm {\bf j}} 
    128          +  \frac{1}{e_{3w} }\delta _{k+1/2} \left[ q \right]\;\,{\rm {\bf k}} 
     128\nabla q\equiv    \frac{1}{e_{1u} } \delta _{i+1/2 } [q] \;\,{\rm {\bf i}} 
     129      +  \frac{1}{e_{2v} } \delta _{j+1/2 } [q] \;\,{\rm {\bf j}} 
     130      +  \frac{1}{e_{3w}} \delta _{k+1/2} [q] \;\,{\rm {\bf k}} 
    129131\end{equation} 
    130132\begin{multline} \label{Eq_DOM_lap} 
    131 \Delta q\equiv \frac{1}{e_{1T} {\kern 1pt}e_{2T} {\kern 1pt}e_{3T} }\;\left(  
    132 {\delta _i \left[ {\frac{e_{2u} e_{3u} }{e_{1u} }\;\delta _{i+1/2}  
    133 \left[ q \right]} \right] 
    134 +\delta _j \left[ {\frac{e_{1v} e_{3v} }{e_{2v}  
    135 }\;\delta _{j+1/2} \left[ q \right]} \right]\;} \right)     \\ 
    136 +\frac{1}{e_{3T} }\delta _k \left[ {\frac{1}{e_{3w} }\;\delta _{k+1/2}  
    137 \left[ q \right]} \right] 
     133\Delta q\equiv \frac{1}{e_{1t}\,e_{2t}\,e_{3t} }  
     134       \;\left(          \delta_i  \left[ \frac{e_{2u}\,e_{3u}} {e_{1u}} \;\delta_{i+1/2} [q] \right] 
     135+                        \delta_j  \left[ \frac{e_{1v}\,e_{3v}}  {e_{2v}} \;\delta_{j+1/2} [q] \right] \;  \right)      \\ 
     136+\frac{1}{e_{3t}} \delta_k \left[ \frac{1}{e_{3w} }                     \;\delta_{k+1/2} [q] \right] 
    138137\end{multline} 
    139138 
    140 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  
    141 components defined at $vw$-, $uw$, and $f$-points, and its divergence defined  
    142 at $T$-points: 
    143 \begin{equation} \label{Eq_DOM_curl} 
    144 \begin{split} 
    145  \nabla \times {\rm {\bf A}}\equiv \frac{1}{e_{2v} {\kern 1pt}e_{3vw}  
    146 }{\kern 1pt}\,\;\left( {\delta _{j+1/2} \left[ {e_{3w} a_3 } \right]-\delta  
    147 _{k+1/2} \left[ {e_{2v} a_2 } \right]} \right)  &\;\;{\rm {\bf i}} \\  
    148  +\frac{1}{e_{2u} {\kern 1pt}e_{3uw} }\;\left( {\delta _{k+1/2} \left[ {e_{1u} a_1 }  
    149 \right]-\delta _{i+1/2} \left[ {e_{3w} a_3 } \right]} \right)  &\;\;{\rm{\bf j}} \\  
    150  +\frac{e_{3f} }{e_{1f} {\kern 1pt}e_{2f} }\,{\kern 1pt}\;\left( {\delta  
    151 _{i+1/2} \left[ {e_{2v} a_2 } \right]-\delta _{j+12} \left[ {e_{1u} a_1 } \right]}  
    152 \right)  &\;\;{\rm {\bf k}} 
    153  \end{split} 
    154 \end{equation} 
     139Following \eqref{Eq_PE_curl} and \eqref{Eq_PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$  
     140defined at vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$,  
     141and $f$-points, and its divergence defined at $t$-points: 
     142\begin{eqnarray}  \label{Eq_DOM_curl} 
     143 \nabla \times {\rm {\bf A}}\equiv & 
     144      \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} \\  
     145 +& \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} \\  
     146 +& \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} 
     147 \end{eqnarray} 
    155148\begin{equation} \label{Eq_DOM_div} 
    156 \nabla \cdot {\rm {\bf A}}=\frac{1}{e_{1T} e_{2T} e_{3T} }\left( {\delta  
    157 _i \left[ {e_{2u} e_{3u} a_1 } \right]+\delta _j \left[ {e_{1v} e_{3v} a_2 }  
    158 \right]} \right)+\frac{1}{e_{3T} }\delta _k \left[ {a_3 } \right] 
     149\nabla \cdot \rm{\bf A}=\frac{1}{e_{1t}\,e_{2t}\,e_{3t}}\left( \delta_i \left[e_{2u}\,e_{3u}\,a_1 \right] 
     150                                                                                         +\delta_j \left[e_{1v}\,e_{3v}\,a_2 \right] \right)+\frac{1}{e_{3t} }\delta_k \left[a_3 \right] 
    159151\end{equation} 
    160152 
     
    164156horizontal location of a grid point. For example \eqref{Eq_DOM_div} reduces to:  
    165157\begin{equation*} 
    166 \nabla \cdot {\rm {\bf A}}=\frac{1}{e_{1T} e_{2T} }\left( {\delta  
    167 _i \left[ {e_{2u} a_1 } \right]+\delta _j \left[ {e_{1v}  a_2 }  
    168 \right]} \right)+\frac{1}{e_{3T} }\delta _k \left[ {a_3 } \right] 
     158\nabla \cdot \rm{\bf A}=\frac{1}{e_{1t}\,e_{2t}} \left( \delta_i \left[e_{2u}\,a_1 \right]  
     159                                                                              +\delta_j \left[e_{1v}\, a_2 \right]  \right) 
     160                                                     +\frac{1}{e_{3t}} \delta_k \left[             a_3 \right] 
    169161\end{equation*} 
    170162 
     
    172164for a quantity $q$ which is a masked field (i.e. equal to zero inside solid area): 
    173165\begin{equation} \label{DOM_bar} 
    174 \bar q   = \frac{1}{H}\int_{k^b}^{k^o} {q\;e_{3q} \,dk}  
     166\bar q   =         \frac{1}{H}    \int_{k^b}^{k^o} {q\;e_{3q} \,dk}  
    175167      \equiv \frac{1}{H_q }\sum\limits_k {q\;e_{3q} } 
    176168\end{equation} 
     
    189181 
    190182It is straightforward to demonstrate that these properties are verified locally in  
    191 discrete form as soon as the scalar $q$ is taken at $T$-points and the vector  
     183discrete form as soon as the scalar $q$ is taken at $t$-points and the vector  
    192184\textbf{A} has its components defined at vector points $(u,v,w)$. 
    193185 
     
    230222The array representation used in the \textsc{Fortran} code requires an integer  
    231223indexing while the analytical definition of the mesh (see \S\ref{DOM_cell}) is  
    232 associated with the use of integer values for $T$-points and both integer and  
     224associated with the use of integer values for $t$-points and both integer and  
    233225integer and a half values for all the other points. Therefore a specific integer  
    234 indexing must be defined for points other than $T$-points ($i.e.$ velocity and  
     226indexing must be defined for points other than $t$-points ($i.e.$ velocity and  
    235227vorticity grid-points). Furthermore, the direction of the vertical indexing has  
    236228been changed so that the surface level is at $k=1$. 
     
    242234\label{DOM_Num_Index_hor} 
    243235 
    244 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  
    245 and the eastward $u$-point (northward $v$-point) have the same index  
    246 (see the dashed area in Fig.\ref{Fig_index_hor}). A $T$-point and its  
    247 nearest northeast $f$-point have the same $i$-and $j$-indices. 
     236The indexing in the horizontal plane has been chosen as shown in Fig.\ref{Fig_index_hor}.  
     237For an increasing $i$ index ($j$ index), the $t$-point and the eastward $u$-point  
     238(northward $v$-point) have the same index (see the dashed area in Fig.\ref{Fig_index_hor}).  
     239A $t$-point and its nearest northeast $f$-point have the same $i$-and $j$-indices. 
    248240 
    249241% ----------------------------------- 
     
    257249to the indexing used in the semi-discrete equations and given in \S\ref{DOM_cell}.  
    258250The sea surface corresponds to the $w$-level $k=1$ which is the same index  
    259 as $T$-level just below (Fig.\ref{Fig_index_vert}). The last $w$-level ($k=jpk$)  
     251as $t$-level just below (Fig.\ref{Fig_index_vert}). The last $w$-level ($k=jpk$)  
    260252either corresponds to the ocean floor or is inside the bathymetry while the last  
    261 $T$-level is always inside the bathymetry (Fig.\ref{Fig_index_vert}). Note that  
    262 for an increasing $k$ index, a $w$-point and the $T$-point just below have the  
     253$t$-level is always inside the bathymetry (Fig.\ref{Fig_index_vert}). Note that  
     254for an increasing $k$ index, a $w$-point and the $t$-point just below have the  
    263255same $k$ index, in opposition to what is done in the horizontal plane where  
    264 it is the $T$-point and the nearest velocity points in the direction of the horizontal  
     256it is the $t$-point and the nearest velocity points in the direction of the horizontal  
    265257axis that have the same $i$ or $j$ index (compare the dashed area in  
    266258Fig.\ref{Fig_index_hor} and \ref{Fig_index_vert}). Since the scale factors are  
     
    309301\S\ref{LBC_mpp}). 
    310302 
     303 
     304$\ $\newline    % force a new ligne 
     305 
    311306% ================================================================ 
    312307% Domain: Horizontal Grid (mesh)  
     
    341336factors in the horizontal plane as follows: 
    342337\begin{flalign*} 
    343 \lambda_T &\equiv \text{glamt}= \lambda(i)     & \varphi_T &\equiv \text{gphit} = \varphi(j)\\ 
     338\lambda_t &\equiv \text{glamt}= \lambda(i)     & \varphi_t &\equiv \text{gphit} = \varphi(j)\\ 
    344339\lambda_u &\equiv \text{glamu}= \lambda(i+1/2)& \varphi_u &\equiv \text{gphiu}= \varphi(j)\\ 
    345340\lambda_v &\equiv \text{glamv}= \lambda(i)       & \varphi_v &\equiv \text{gphiv} = \varphi(j+1/2)\\ 
     
    347342\end{flalign*} 
    348343\begin{flalign*} 
    349 e_{1T} &\equiv \text{e1t} = r_a |\lambda'(i)    \; \cos\varphi(j)  |& 
    350 e_{2T} &\equiv \text{e2t} = r_a |\varphi'(j)|  \\ 
     344e_{1t} &\equiv \text{e1t} = r_a |\lambda'(i)    \; \cos\varphi(j)  |& 
     345e_{2t} &\equiv \text{e2t} = r_a |\varphi'(j)|  \\ 
    351346e_{1u} &\equiv \text{e1t} = r_a |\lambda'(i+1/2)   \; \cos\varphi(j)  |& 
    352347e_{2u} &\equiv \text{e2t} = r_a |\varphi'(j)|\\ 
     
    359354considered and $r_a$ is the earth radius (defined in \mdl{phycst} along with  
    360355all universal constants). Note that the horizontal position of and scale factors  
    361 at $w$-points are exactly equal to those of $T$-points, thus no specific arrays  
     356at $w$-points are exactly equal to those of $t$-points, thus no specific arrays  
    362357are defined at $w$-points.  
    363358 
    364359Note that the definition of the scale factors ($i.e.$ as the analytical first derivative  
    365360of the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) is  
    366 specific to the \NEMO model \citep{Marti1992}. As an example, $e_{1T}$ is defined  
    367 locally at a $T$-point, whereas many other models on a C grid choose to define  
     361specific to the \NEMO model \citep{Marti_al_JGR92}. As an example, $e_{1t}$ is defined  
     362locally at a $t$-point, whereas many other models on a C grid choose to define  
    368363such a scale factor as the distance between the $U$-points on each side of the  
    369 $T$-point. Relying on an analytical transformation has two advantages: firstly, there  
     364$t$-point. Relying on an analytical transformation has two advantages: firstly, there  
    370365is no ambiguity in the scale factors appearing in the discrete equations, since they  
    371366are first introduced in the continuous equations; secondly, analytical transformations  
     
    380375in the vertical, and (b) analytically derived grid-point position and scale factors. For  
    381376both grids here,  the same $w$-point depth has been chosen but in (a) the  
    382 $T$-points are set half way between $w$-points while in (b) they are defined from  
     377$t$-points are set half way between $w$-points while in (b) they are defined from  
    383378an analytical function: $z(k)=5\,(i-1/2)^3 - 45\,(i-1/2)^2 + 140\,(i-1/2) - 150$.  
    384379Note the resulting difference between the value of the grid-size $\Delta_k$ and  
     
    422417and \pp{ppgphi0}). Note that for the Mercator grid the user need only provide  
    423418an approximate starting latitude: the real latitude will be recalculated analytically,  
    424 in order to ensure that the equator corresponds to line passing through $T$-  
     419in order to ensure that the equator corresponds to line passing through $t$-  
    425420and $u$-points.   
    426421 
     
    431426and given by the parameters \pp{ppe1\_m} and \pp{ppe2\_m} respectively.  
    432427The zonal grid coordinate (\textit{glam} arrays) is in kilometers, starting at zero  
    433 with the first $T$-point. The meridional coordinate (gphi. arrays) is in kilometers,  
    434 and the second $T$-point corresponds to coordinate $gphit=0$. The input  
     428with the first $t$-point. The meridional coordinate (gphi. arrays) is in kilometers,  
     429and the second $t$-point corresponds to coordinate $gphit=0$. The input  
    435430parameter \pp{ppglam0} is ignored. \pp{ppgphi0} is used to set the reference  
    436431latitude for computation of the Coriolis parameter. In the case of the beta plane,  
     
    460455the output grid written when $\np{nmsh} \not=0$ is no more equal to the input grid. 
    461456 
     457$\ $\newline    % force a new ligne 
     458 
    462459% ================================================================ 
    463460% Domain: Vertical Grid (domzgr) 
     
    467464\label{DOM_zgr} 
    468465%-----------------------------------------nam_zgr & namdom------------------------------------------- 
    469 \namdisplay{nam_zgr}  
     466\namdisplay{namzgr}  
    470467\namdisplay{namdom}  
    471468%------------------------------------------------------------------------------------------------------------- 
     
    593590 
    594591The reference coordinate transformation $z_0 (k)$ defines the arrays $gdept_0$  
    595 and $gdepw_0$ for $T$- and $w$-points, respectively. As indicated on  
     592and $gdepw_0$ for $t$- and $w$-points, respectively. As indicated on  
    596593Fig.\ref{Fig_index_vert} \jp{jpk} is the number of $w$-levels. $gdepw_0(1)$ is the  
    597 ocean surface. There are at most \jp{jpk}-1 $T$-points inside the ocean, the  
    598 additional $T$-point at $jk=jpk$ is below the sea floor and is not used.  
    599 The vertical location of $w$- and $T$-levels is defined from the analytic expression  
     594ocean surface. There are at most \jp{jpk}-1 $t$-points inside the ocean, the  
     595additional $t$-point at $jk=jpk$ is below the sea floor and is not used.  
     596The vertical location of $w$- and $t$-levels is defined from the analytic expression  
    600597of the depth $z_0(k)$ whose analytical derivative with respect to $k$ provides the  
    601598vertical scale factors. The user must provide the analytical expression of both  
     
    663660\begin{center} \begin{tabular}{c||r|r|r|r} 
    664661\hline 
    665 \textbf{LEVEL}& \textbf{GDEPT}& \textbf{GDEPW}& \textbf{E3T }& \textbf{E3W  } \\ \hline 
     662\textbf{LEVEL}& \textbf{gdept}& \textbf{gdepw}& \textbf{e3t }& \textbf{e3w  } \\ \hline 
    6666631  &  \textbf{  5.00}   &       0.00 & \textbf{ 10.00} &  10.00 \\   \hline 
    6676642  &  \textbf{15.00} &    10.00 &   \textbf{ 10.00} &  10.00 \\   \hline 
     
    728725Two variables in the namdom namelist are used to define the partial step  
    729726vertical grid. The mimimum water thickness (in meters) allowed for a cell  
    730 partially filled with bathymetry at level jk is the minimum of \np{e3zpsmin}  
    731 (thickness in meters, usually $20~m$) or $e_{3t}(jk)*\np{e3zps\_rat}$ (a fraction,  
     727partially filled with bathymetry at level jk is the minimum of \np{rn\_e3zps\_min}  
     728(thickness in meters, usually $20~m$) or $e_{3t}(jk)*\np{rn\_e3zps\_rat}$ (a fraction,  
    732729usually 10\%, of the default thickness $e_{3t}(jk)$). 
    733730 
     
    738735% ------------------------------------------------------------------------------------------------------------- 
    739736\subsection   [$s$-coordinate (\np{ln\_sco})] 
    740          {$s$-coordinate (\np{ln\_sco}=true)} 
     737           {$s$-coordinate (\np{ln\_sco}=true)} 
    741738\label{DOM_sco} 
    742739%------------------------------------------nam_zgr_sco--------------------------------------------------- 
    743 \namdisplay{nam_zgr_sco}  
     740\namdisplay{namzgr_sco}  
    744741%-------------------------------------------------------------------------------------------------------------- 
    745742In $s$-coordinate (\key{sco} is defined), the depth and thickness of the model  
     
    752749\end{split} 
    753750\end{equation} 
    754 where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $T$-point  
     751where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point  
    755752location in the horizontal and $z_0(k)$ is a function which varies from $0$ at the sea  
    756753surface to $1$ at the ocean bottom. The depth field $h$ is not necessary the ocean  
     
    760757sharp bathymetric gradients. 
    761758 
    762 A new flexible stretching function, modified from \citet{Song1994} is provided as an example: 
     759A new flexible stretching function, modified from \citet{Song_Haidvogel_JCP94} is provided as an example: 
    763760\begin{equation} \label{DOM_sco_function} 
    764761\begin{split} 
     
    801798Whatever the vertical coordinate used, the model offers the possibility of  
    802799representing the bottom topography with steps that follow the face of the  
    803 model cells (step like topography) \citep{Madec1996}. The distribution of  
     800model cells (step like topography) \citep{Madec_al_JPO96}. The distribution of  
    804801the steps in the horizontal is defined in a 2D integer array, mbathy, which  
    805802gives the number of ocean levels ($i.e.$ those that are not masked) at each  
    806 $T$-point. mbathy is computed from the meter bathymetry using the definiton of  
    807 gdept as the number of $T$-points which gdept $\leq$ bathy. Note that in version  
     803$t$-point. mbathy is computed from the meter bathymetry using the definiton of  
     804gdept as the number of $t$-points which gdept $\leq$ bathy. Note that in version  
    808805NEMO v2.3, the user still has to provide the "level" bathymetry in a NetCDF 
    809806file when using the full step option (\np{ln\_zco}), rather than the bathymetry  
     
    817814domain (\np{ln\_dynspg\_rl}=.true. and \key{island} is defined), the \textit{mbathy}  
    818815array must be provided and takes values from $-N$ to \jp{jpk}-1. It provides the  
    819 following information: $mbathy(i,j) = -n, \ n \in \left] 0,N \right]$, $T$-points are  
    820 land points on the $n^{th}$ island ; $mbathy(i,j) =0$, $T$-points are land points  
    821 on the main land (continent) ; $mbathy(i,j) =k$, the first $k$ $T$- and $w$-points  
     816following information: $mbathy(i,j) = -n, \ n \in \left] 0,N \right]$, $t$-points are  
     817land points on the $n^{th}$ island ; $mbathy(i,j) =0$, $t$-points are land points  
     818on the main land (continent) ; $mbathy(i,j) =k$, the first $k$ $t$- and $w$-points  
    822819are ocean points, the others are points below the ocean floor.  
    823820 
     
    849846%%% 
    850847 
     848$\ $\newline    % force a new ligne 
     849 
    851850% ================================================================ 
    852851% Time Discretisation 
     
    860859   x^{t+\Delta t} = x^{t-\Delta t} + 2 \, \Delta t \  \text{RHS}_x^{t-\Delta t,t,t+\Delta t} 
    861860\end{equation}  
    862 where $x$ stands for $u$, $v$, $T$ or $S$; RHS is the Right-Hand-Side of the  
     861where $x$ stands for $u$, $v$, $t$ or $S$; RHS is the Right-Hand-Side of the  
    863862corresponding time evolution equation; $\Delta t$ is the time step; and the  
    864863superscripts indicate the time at which a quantity is evaluated. Each term of the  
     
    995994 
    996995This is diffusive in time and conditionally stable. For example, the  
    997 conditions for stability of second and fourth order horizontal diffusion schemes are \citep{Griffies2004}: 
     996conditions for stability of second and fourth order horizontal diffusion schemes are \citep{Griffies_Bk04}: 
    998997\begin{equation} \label{Eq_DOM_nxt_euler_stability} 
    999998A^h < \left\{ 
     
    10401039approximation of the temperature equation is: 
    10411040\begin{equation} \label{Eq_DOM_nxt_imp_zdf} 
    1042 \frac{T(k)^{t+1}-T(k)^{t-1}}{2\;\Delta t}\equiv \text{RHS}+\frac{1}{e_{3T} }\delta  
     1041\frac{T(k)^{t+1}-T(k)^{t-1}}{2\;\Delta t}\equiv \text{RHS}+\frac{1}{e_{3t} }\delta  
    10431042_k \left[ {\frac{A_w^{vT} }{e_{3w} }\delta _{k+1/2} \left[ {T^{t+1}} \right]}  
    10441043\right] 
     
    11091108}        %% end add 
    11101109 
     1110 
     1111 
     1112 
     1113 
     1114Implicit time stepping in case of variable volume thickness. 
     1115 
     1116Tracer case (NB for momentum in vector invariant form take care!) 
     1117 
     1118\begin{flalign*} 
     1119&\frac{\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1}}{2\Delta t} 
     1120\equiv \text{RHS}+ \delta _k \left[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta _{k+1/2} \left[ {T^{t+1}} \right]}  
     1121\right]      \\ 
     1122&\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1} 
     1123\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]}  
     1124\right]      \\ 
     1125&\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1} 
     1126\equiv 2\Delta t \ \text{RHS} 
     1127+ 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} ] 
     1128                          - \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} [ T_k       ^{t+1} - T_{k-1}^{t+1} ]  \right\}     \\ 
     1129&\\ 
     1130&\left( e_{3t}\,T \right)_k^{t+1} 
     1131-  {2\Delta t} \           \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}                  T_{k+1}^{t+1}  
     1132+ {2\Delta t} \ \left\{  \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}  
     1133                            +  \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2}     \right\}   T_{k    }^{t+1} 
     1134-  {2\Delta t} \           \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2}                  T_{k-1}^{t+1}      \\ 
     1135&\equiv \left( e_{3t}\,T \right)_k^{t-1} + {2\Delta t} \ \text{RHS}    \\ 
     1136% 
     1137\end{flalign*} 
     1138 
     1139\begin{flalign*} 
     1140\allowdisplaybreaks 
     1141\intertext{ Tracer case } 
     1142% 
     1143&  \qquad \qquad  \quad   -  {2\Delta t}                  \ \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}    
     1144                                                                                                      \qquad \qquad \qquad  \qquad  T_{k+1}^{t+1}   \\ 
     1145&+ {2\Delta t} \ \biggl\{  (e_{3t})_{k   }^{t+1}  \bigg. +    \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}   
     1146                                                                               +   \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} \bigg. \biggr\}  \ \ \ T_{k   }^{t+1}  &&\\ 
     1147& \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}      
     1148\ \equiv \ \left( e_{3t}\,T \right)_k^{t-1} + {2\Delta t} \ \text{RHS}  \\ 
     1149% 
     1150\end{flalign*} 
     1151\begin{flalign*} 
     1152\allowdisplaybreaks 
     1153\intertext{ Tracer content case } 
     1154% 
     1155& -  {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}   &\\ 
     1156& + {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}}  
     1157                                    + & \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}  &\\ 
     1158& -  {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}    
     1159\equiv \left( e_{3t}\,T \right)_k^{t-1} + {2\Delta t} \ \text{RHS}  & 
     1160\end{flalign*} 
     1161 
Note: See TracChangeset for help on using the changeset viewer.