New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 1223 – NEMO

Changeset 1223


Ignore:
Timestamp:
2008-11-26T13:12:16+01:00 (15 years ago)
Author:
gm
Message:

minor corrections in the Appendix from Steven see ticket #282

Location:
trunk/DOC/TexFiles/Chapters
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/DOC/TexFiles/Chapters/Annex_A.tex

    r996 r1223  
    88 
    99 
    10 In order to establish the set of Primitive Equation in curvilinear $s$-coordinates ($i.e.$  
    11 orthogonal curvilinear coordinate in the horizontal and $s$-coordinate in the vertical), we  
    12 start from the set of equation established in \S\ref{PE_zco_Eq} for the special case  
    13 $k = z$ and thus $e_3 = 1$, and we introduce an arbitrary vertical coordinate  
    14 $s = s(i,j,z,t)$. Let us define a new vertical scale factor by $e_3 = \partial z / \partial s$  
    15 (which now depends on $(i,j,z,t)$) and the horizontal slope of $s$-surfaces by : 
     10In order to establish the set of Primitive Equation in curvilinear $s$-coordinates  
     11($i.e.$ an orthogonal curvilinear coordinate in the horizontal and $s$-coordinate  
     12in the vertical), we start from the set of equations established in \S\ref{PE_zco_Eq}  
     13for the special case $k = z$ and thus $e_3 = 1$, and we introduce an arbitrary  
     14vertical coordinate $s = s(i,j,z,t)$. Let us define a new vertical scale factor by  
     15$e_3 = \partial z / \partial s$ (which now depends on $(i,j,z,t)$) and the horizontal  
     16slope of $s$-surfaces by : 
    1617\begin{equation} \label{Apdx_A_s_slope} 
    1718\sigma _1 =\frac{1}{e_1 }\;\left. {\frac{\partial z}{\partial i}} \right|_s  
    … …  
    2021\end{equation} 
    2122 
    22 The chain rule to establish the model equations in the curvilinear $s$-coordinate system  
    23 is: 
     23The chain rule to establish the model equations in the curvilinear $s$-coordinate  
     24system is: 
    2425\begin{equation} \label{Apdx_A_s_chain_rule} 
    2526\begin{aligned} 
    … …  
    4142\end{equation} 
    4243 
    43 In particular applying the time derivative chain rule to $z$ provide the expression of $w_s$,  the vertical velocity of the $s-$surfaces: 
     44In particular applying the time derivative chain rule to $z$ provides the  
     45expression for $w_s$,  the vertical velocity of the $s-$surfaces: 
    4446\begin{equation} \label{Apdx_A_w_in_s} 
    4547w_s   =  \left.   \frac{\partial z }{\partial t}   \right|_s  
    … …  
    5456\label{Apdx_B_continuity} 
    5557 
    56 Using (\ref{Apdx_A_s_chain_rule}) and the fact that the horizontal scale factors $e_1$ and $e_2$ do not depend on the vertical coordinate, the divergence of the velocity relative to the ($i$,$j$,$z$) coordinate system is transformed as follows: 
     58Using (\ref{Apdx_A_s_chain_rule}) and the fact that the horizontal scale factors  
     59$e_1$ and $e_2$ do not depend on the vertical coordinate, the divergence of  
     60the velocity relative to the ($i$,$j$,$z$) coordinate system is transformed as follows: 
    5761 
    5862\begin{align*} 
    … …  
    109113 \end{align*}  
    110114  
    111 Here, $w$ is the vertical velocity relative to the $z-$coordinate system. Introducing the dia-surface velocity component, $\omega $, defined as the velocity relative to the moving $s$-surfaces and normal to them: 
     115Here, $w$ is the vertical velocity relative to the $z-$coordinate system.  
     116Introducing the dia-surface velocity component, $\omega $, defined as  
     117the velocity relative to the moving $s$-surfaces and normal to them: 
    112118\begin{equation} \label{Apdx_A_w_s} 
    113119\omega  = w - w_s - \sigma _1 \,u - \sigma _2 \,v    \\ 
    114120\end{equation} 
    115 with $w_s$ given by \eqref{Apdx_A_w_in_s}, we obtain the expression of the divergence of the velocity in the curvilinear $s$-coordinate system: 
     121with $w_s$ given by \eqref{Apdx_A_w_in_s}, we obtain the expression for  
     122the divergence of the velocity in the curvilinear $s$-coordinate system: 
    116123\begin{align*} \label{Apdx_A_A4} 
    117124\nabla \cdot {\rm {\bf U}}  
    … …  
    142149\end{align*} 
    143150 
    144 As a result, the continuity equation \eqref{Eq_PE_continuity} in $s$-coordinates becomes: 
     151As a result, the continuity equation \eqref{Eq_PE_continuity} in the  
     152$s$-coordinates becomes: 
    145153\begin{equation} \label{Apdx_A_A5} 
    146154\frac{1}{e_3 } \frac{\partial e_3}{\partial t}  
    … …  
    157165\label{Apdx_B_momentum} 
    158166 
    159 Let us consider \eqref{Eq_PE_dyn_vect}, the first component of the  
    160 momentum equation in the vector invariant form (similar manipulations can be performed on the second one). Its non linear term can be transformed  
    161 as follows: 
     167Let us consider \eqref{Eq_PE_dyn_vect}, the first component of the momentum  
     168equation in the vector invariant form (similar manipulations can be performed  
     169on the second component). Its non-linear term can be transformed as follows: 
    162170 
    163171\begin{align*} 
    … …  
    205213\end{align*} 
    206214 
    207 Therefore, the non-linear terms of the momentum equation have the same form  
    208 in $z-$ and $s-$coordinates but with the addition of the time derivative of the velocity:  
     215Therefore, the non-linear terms of the momentum equation have the same  
     216form in $z-$ and $s-$coordinates but with the addition of the time derivative  
     217of the velocity:  
    209218\begin{multline}  \label{Apdx_A_momentum_NL} 
    210219+\left. \zeta \right|_z v-\frac{1}{2e_1 }\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_z  
    … …  
    224233\end{equation} 
    225234 
    226 An additional term appears in (\ref{Apdx_A_grad_p}) which accounts for the tilt of model  
    227 levels. 
    228  
    229 Introducing \eqref{Apdx_A_momentum_NL} and \eqref{Apdx_A_grad_p} in \eqref{Eq_PE_dyn_vect} and regrouping the time derivative terms in the left hand side, and performing the same manipulation on the second component, we obtain the vector invariant form of momentum equation in $s-$coordinate : 
     235An additional term appears in (\ref{Apdx_A_grad_p}) which accounts for the  
     236tilt of model levels. 
     237 
     238Introducing \eqref{Apdx_A_momentum_NL} and \eqref{Apdx_A_grad_p} in \eqref{Eq_PE_dyn_vect} and regrouping the time derivative terms in the left  
     239hand side, and performing the same manipulation on the second component,  
     240we obtain the vector invariant form of the momentum equations in the  
     241$s-$coordinate : 
    230242\begin{subequations} \label{Apdx_A_dyn_vect} 
    231243\begin{multline} \label{Apdx_A_PE_dyn_vect_u} 
    … …  
    249261\end{subequations} 
    250262 
    251 It has the same form as in $z-$coordinate but the vertical scale factor that has appeared inside the time derivative. The form of the vertical physics and forcing terms remain unchanged. The form of the lateral physics is discussed in appendix~\ref{Apdx_B}.   
     263It has the same form as in the $z-$coordinate but for the vertical scale factor  
     264that has appeared inside the time derivative. The form of the vertical physics  
     265and forcing terms remains unchanged. The form of the lateral physics is  
     266discussed in appendix~\ref{Apdx_B}.   
    252267 
    253268% ================================================================ 
    … …  
    257272\label{Apdx_B_tracer} 
    258273 
    259 The tracer equation is obtained using the same calculation as for the  
    260 continuity equation and then regrouping the time derivative terms in the left hand side : 
     274The tracer equation is obtained using the same calculation as for the continuity  
     275equation and then regrouping the time derivative terms in the left hand side : 
    261276 
    262277\begin{multline} \label{Apdx_A_tracer} 
    … …  
    269284 
    270285 
    271 The expression of the advection term is a straight consequence of (A.4), the  
    272 expression of the 3D divergence in $s$-coordinates established above.  
    273  
     286The expression for the advection term is a straight consequence of (A.4), the  
     287expression of the 3D divergence in the $s$-coordinates established above.  
     288 

  • trunk/DOC/TexFiles/Chapters/Annex_B.tex

    r994 r1223  
    1313 
    1414 
    15 In $z$-coordinate, the horizontal/vertical second order tracer diffusive operator is given by: 
     15In the $z$-coordinate, the horizontal/vertical second order tracer diffusion operator  
     16is given by: 
    1617\begin{multline} \label{Apdx_B1} 
    1718 D^T = \frac{1}{e_1 \, e_2}      \left[  
    … …  
    2223\end{multline} 
    2324 
    24 In $s$-coordinate, we defined the slopes of $s$-surfaces, $\sigma_1$ and $\sigma_2$ by (A.1) and the vertical/horizontal ratio of diffusive coefficient by $\epsilon = A^{vT} / A^{lT}$. The diffusive operator is given by: 
     25In the $s$-coordinate, we defined the slopes of $s$-surfaces, $\sigma_1$ and  
     26$\sigma_2$ by (!!!A.1!!!) and the vertical/horizontal ratio of diffusion coefficient  
     27by $\epsilon = A^{vT} / A^{lT}$. The diffusion operator is given by: 
    2528 
    2629\begin{equation} \label{Apdx_B2} 
    … …  
    3437\end{array} }} \right) 
    3538\end{equation} 
    36 or in expended form: 
     39or in expanded form: 
    3740\begin{multline} \label{Apdx_B3} 
    3841D^T=\frac{1}{e_1\,e_2\,e_3 }\;\left[ {\quad \; \; e_2\,e_3\,A^{lT} \;\left.  
    … …  
    4548\end{multline} 
    4649 
    47  
    48 Equation \eqref{Apdx_B2} (or equivalently \eqref{Apdx_B3}) is obtained from \eqref{Apdx_B1} without any additional assumption. Indeed, for the special case $k=z$ and thus $e_3 =1$, we introduce an arbitrary vertical coordinate $s = s (i,j,z)$ as in Appendix~\ref{Apdx_A} and use \eqref{Apdx_A_s_slope} and \eqref{Apdx_A_s_chain_rule}. Since no cross horizontal derivate $\partial _i \partial _j $ appears in \eqref{Apdx_B1}, the ($i$,$z$) and ($j$,$z$) planes are independent. The demonstration can then be done for the ($i$,$z$)~$\to$~($j$,$s$) transformation without any loss of generality: 
     50Equation \eqref{Apdx_B2} (or equivalently \eqref{Apdx_B3}) is obtained from  
     51\eqref{Apdx_B1} without any additional assumption. Indeed, for the special  
     52case $k=z$ and thus $e_3 =1$, we introduce an arbitrary vertical coordinate  
     53$s = s (i,j,z)$ as in Appendix~\ref{Apdx_A} and use \eqref{Apdx_A_s_slope}  
     54and \eqref{Apdx_A_s_chain_rule}. Since no cross horizontal derivative  
     55$\partial _i \partial _j $ appears in \eqref{Apdx_B1}, the ($i$,$z$) and ($j$,$z$)  
     56planes are independent. The derivation can then be demonstrated for the  
     57($i$,$z$)~$\to$~($j$,$s$) transformation without any loss of generality: 
    4958 
    5059\begin{equation*} 
    … …  
    8796\end{multline*} 
    8897 
    89 Since the horizontal scale factor do not depend on the vertical coordinate, the last term of the first line and the first term of the last line cancel, while the second line reduces to a single vertical derivative, so it becomes: 
     98Since the horizontal scale factors do not depend on the vertical coordinate,  
     99the last term of the first line and the first term of the last line cancel, while  
     100the second line reduces to a single vertical derivative, so it becomes: 
    90101 
    91102\begin{multline*} 
    … …  
    94105\end{multline*} 
    95106 
    96 in other words, the horizontal Laplacian operator in the ($i$,$s$)plane takes the following expression : 
     107in other words, the horizontal Laplacian operator in the ($i$,$s$) plane takes the following form : 
    97108 
    98109\begin{equation*} 
    … …  
    123134 
    124135 
    125 The iso/diapycnal diffusive tensor $\textbf {A}_{\textbf I}$ expressed in the  ($i$,$j$,$k$) curvilinear coordinate system in which the equations of the ocean circulation model are formulated, takes the following expression \citep{Redi_JPO82}: 
     136The iso/diapycnal diffusive tensor $\textbf {A}_{\textbf I}$ expressed in the ($i$,$j$,$k$)  
     137curvilinear coordinate system in which the equations of the ocean circulation model are  
     138formulated, takes the following form \citep{Redi_JPO82}: 
    126139 
    127140\begin{equation*} 
    … …  
    141154\end{equation*} 
    142155 
    143 In practice, the isopycnal slopes are generally less than $10^{-2}$ in the ocean, so $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{Cox1987}: 
     156In practice, the isopycnal slopes are generally less than $10^{-2}$ in the ocean, so  
     157$\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{Cox1987}: 
    144158\begin{equation*} 
    145159{\textbf{A}_{\textbf{I}}} \approx A^{lT} 
    … …  
    151165\end{equation*} 
    152166 
    153 The resulting isopycnal operator conserves the quantity and dissipates its square. The demonstration of the first property is trivial as \eqref{Apdx_B2} is the divergence of fluxes. Let us demonstrate the second one: 
     167The resulting isopycnal operator conserves the quantity and dissipates its square.  
     168The demonstration of the first property is trivial as \eqref{Apdx_B2} is the divergence  
     169of fluxes. Let us demonstrate the second one: 
    154170\begin{equation*} 
    155171\iiint\limits_D T\;\nabla .\left( {\textbf{A}}_{\textbf{I}} \nabla T \right)\,dv = -\iiint\limits_D \nabla T\;.\left( {\textbf{A}}_{\textbf{I}} \nabla T \right)\,dv 
    … …  
    169185the property becomes obvious.  
    170186 
    171 The resulting diffusive operator in $z$-coordinates has the following  
    172 expression : 
     187The resulting diffusion operator in $z$-coordinate has the following form : 
    173188\begin{multline*} \label{Apdx_B_ldfiso} 
    174189 D^T=\frac{1}{e_1 e_2 }\left\{ {\;\frac{\partial }{\partial i}\left[ {A_h \left( {\frac{e_2 }{e_1 }\frac{\partial T}{\partial i}-a_1 \frac{e_2}{e_3}\frac{\partial T}{\partial k}} \right)} \right]\;} \right.\;\; \\  
    … …  
    177192\end{multline*} 
    178193 
    179 It has to be emphasised that the simplification introduced leads to a decoupling between ($i$,$z$) and ($j$,$z$) planes. The operator has therefore the same expression as \eqref{Apdx_B3}, the diffusive operator obtained for geopotential diffusion in $s$-coordinate.  
     194It has to be emphasised that the simplification introduced, leads to a decoupling  
     195between ($i$,$z$) and ($j$,$z$) planes. The operator has therefore the same  
     196expression as \eqref{Apdx_B3}, the diffusion operator obtained for geopotential  
     197diffusion in the $s$-coordinate.  
    180198 
    181199% ================================================================ 
    … …  
    185203\label{Apdx_B_3} 
    186204 
    187 The second order momentum diffusive operator (Laplacian) in $z$-coordinate is found by applying \eqref{Eq_PE_lap_vector}, the expression of the Laplacian of a vector,  to the horizontal velocity vector :  
     205The second order momentum diffusion operator (Laplacian) in the $z$-coordinate  
     206is found by applying \eqref{Eq_PE_lap_vector}, the expression for the Laplacian  
     207of a vector,  to the horizontal velocity vector :  
    188208\begin{align*} 
    189209\Delta {\textbf{U}}_h  
    … …  
    221241\end{array} }} \right) 
    222242\end{align*} 
    223 Using \eqref{Eq_PE_div}, the definition of the horizontal divergence, the third componant of the second vector is obviously zero and thus : 
     243Using \eqref{Eq_PE_div}, the definition of the horizontal divergence, the third  
     244componant of the second vector is obviously zero and thus : 
    224245\begin{equation*} 
    225246\Delta {\textbf{U}}_h = \nabla _h \left( \chi \right) - \nabla _h \times \left( \zeta \right) + \frac {1}{e_3 } \frac {\partial }{\partial k} \left( {\frac {1}{e_3 } \frac{\partial {\textbf{ U}}_h }{\partial k}} \right) 
    226247\end{equation*} 
    227248 
    228 Note that this operator ensures a full separation between the vorticity and horizontal divergence fields (see Appendix~\ref{Apdx_C}). It is only equal to a Laplacian applied on each component in Cartesian coordinate, not on the sphere. 
     249Note that this operator ensures a full separation between the vorticity and horizontal  
     250divergence fields (see Appendix~\ref{Apdx_C}). It is only equal to a Laplacian  
     251applied to each component in Cartesian coordinates, not on the sphere. 
    229252 
    230253The horizontal/vertical second order (Laplacian type) operator used to diffuse  
    231 horizontal momentum in $z$-coordinate takes therefore the following expression : 
     254horizontal momentum in the $z$-coordinate therefore takes the following form : 
    232255\begin{equation} \label{Apdx_B_Lap_U} 
    233256 {\textbf{D}}^{\textbf{U}} =  
    … …  
    237260            \frac{\partial {\rm {\bf U}}_h }{\partial k}} \right) \\  
    238261\end{equation} 
    239 that is in expanded form: 
     262that is, in expanded form: 
    240263\begin{align*} 
    241264D^{\textbf{U}}_u  
    … …  
    249272\end{align*} 
    250273 
    251 Note Bene: introducing a rotation in \eqref{Apdx_B_Lap_U} does not lead to any useful expression for the iso/diapycnal Laplacian operator in $z$-coordinate. Similarly, we did not found an expression of practical use for the geopotential horizontal/vertical Laplacian operator in $s$-coordinate. Generally, \eqref{Apdx_B_Lap_U} is used in both $z$- and $s$-coordinate system, that is a Laplacian diffusion is applied on momentum along the coordinate directions. 
     274Note Bene: introducing a rotation in \eqref{Apdx_B_Lap_U} does not lead to a  
     275useful expression for the iso/diapycnal Laplacian operator in the $z$-coordinate.  
     276Similarly, we did not found an expression of practical use for the geopotential  
     277horizontal/vertical Laplacian operator in the $s$-coordinate. Generally,  
     278\eqref{Apdx_B_Lap_U} is used in both $z$- and $s$-coordinate systems, that is  
     279a Laplacian diffusion is applied on momentum along the coordinate directions. 

  • trunk/DOC/TexFiles/Chapters/Annex_C.tex

    r999 r1223  
    1717\label{Apdx_C.1} 
    1818 
    19  
    20 First, the boundary condition on the vertical velocity (no flux through the surface and the bottom) is established for the discrete set of momentum equations. Then, it is shown that the non linear terms of the momentum equation are written such that the potential enstrophy of a horizontally non divergent flow is preserved while all the other non-diffusive terms preserve the kinetic energy: the energy is also preserved in practice. In addition, an option is also offer for the vorticity term discretization which provides  
    21 a total kinetic energy conserving discretization for that term.  
    22  
    23 Nota Bene: these properties are established here in the rigid-lid case and for the 2nd order centered scheme. A forthcoming update will be their generalisation to the free surface case 
    24 and higher order scheme. 
     19First, the boundary condition on the vertical velocity (no flux through the surface  
     20and the bottom) is established for the discrete set of momentum equations.  
     21Then, it is shown that the non-linear terms of the momentum equation are written  
     22such that the potential enstrophy of a horizontally non-divergent flow is preserved  
     23while all the other non-diffusive terms preserve the kinetic energy; in practice the  
     24energy is also preserved. In addition, an option is also offered for the vorticity term  
     25discretization which provides a total kinetic energy conserving discretization for  
     26that term.  
     27 
     28Nota Bene: these properties are established here in the rigid-lid case and for the  
     292nd order centered scheme. A forthcoming update will be their generalisation to  
     30the free surface case and higher order scheme. 
    2531 
    2632% ------------------------------------------------------------------------------------------------------------- 
    … …  
    3137 
    3238 
    33 The discrete set of momentum equations used in rigid lid approximation  
     39The discrete set of momentum equations used in the rigid-lid approximation  
    3440automatically satisfies the surface and bottom boundary conditions  
    3541(no flux through the surface and the bottom: $w_{surface} =w_{bottom} =~0$).  
    3642Indeed, taking the discrete horizontal divergence of the vertical sum of the  
    37 horizontal momentum equations (Eqs. (II.2.1) and (II.2.2)~) wheighted by the  
     43horizontal momentum equations (!!!Eqs. (II.2.1) and (II.2.2)!!!) weighted by the  
    3844vertical scale factors, it becomes: 
    3945\begin{flalign*} 
    … …  
    8389 
    8490 
    85 The surface boundary condition associated with the rigid lid approximation ($w_{surface} =0)$ is imposed in the computation of the vertical velocity (II.2.5). Therefore, it turns out to be: 
     91The surface boundary condition associated with the rigid lid approximation  
     92($w_{surface} =0)$ is imposed in the computation of the vertical velocity (!!! II.2.5!!!!).  
     93Therefore, it turns out to be: 
    8694\begin{equation*} 
    8795\frac{\partial } {\partial t}w_{bottom} \equiv 0 
    8896\end{equation*} 
    89 As the bottom velocity is initially set to zero, it remains zero all the time. Symmetrically, if $w_{bottom} =0$ is used in the computation of the vertical velocity (upward integral of the horizontal divergence), the same computation leads to $w_{surface} =0$ as soon as the surface vertical velocity is initially set to zero. 
     97As the bottom velocity is initially set to zero, it remains zero all the time.  
     98Symmetrically, if $w_{bottom} =0$ is used in the computation of the vertical  
     99velocity (upward integral of the horizontal divergence), the same computation  
     100leads to $w_{surface} =0$ as soon as the surface vertical velocity is initially  
     101set to zero. 
    90102 
    91103% ------------------------------------------------------------------------------------------------------------- 
    … …  
    101113\label{Apdx_C_vor}  
    102114 
    103 Potential vorticity is located at $f$-points and defined as: $\zeta / e_{3f}$. The standard discrete formulation of the relative vorticity term obviously conserves potential vorticity (ENS scheme). It also conserves the potential enstrophy for a horizontally non-divergent flow (i.e. $\chi $=0) but not the total kinetic energy. Indeed, using the symmetry or skew symmetry properties of the operators (Eqs \eqref{DOM_mi_adj} and \eqref{DOM_di_adj}), it can be shown that: 
     115Potential vorticity is located at $f$-points and defined as: $\zeta / e_{3f}$.  
     116The standard discrete formulation of the relative vorticity term obviously  
     117conserves potential vorticity (ENS scheme). It also conserves the potential  
     118enstrophy for a horizontally non-divergent flow (i.e. $\chi $=0) but not the  
     119total kinetic energy. Indeed, using the symmetry or skew symmetry properties  
     120of the operators (Eqs \eqref{DOM_mi_adj} and \eqref{DOM_di_adj}), it can  
     121be shown that: 
    104122\begin{equation} \label{Apdx_C_1.1} 
    105123\int_D {\zeta / e_3\,\;{\textbf{k}}\cdot \frac{1} {e_3} \nabla \times \left( {\zeta \;{\textbf{k}}\times {\textbf{U}}_h } \right)\;dv} \equiv 0 
    … …  
    164182\end{align*} 
    165183 
    166 Note that the demonstration is done here for the relative potential  
    167 vorticity but it still hold for the planetary ($f/e_3$) and the total  
    168 potential vorticity $((\zeta +f) /e_3 )$. Another formulation of  
    169 the two components of the vorticity term is optionally offered (ENE scheme) : 
     184Note that the derivation is demonstrated here for the relative potential  
     185vorticity but it applies also to the planetary ($f/e_3$) and the total  
     186potential vorticity $((\zeta +f) /e_3 )$. Another formulation of the two  
     187components of the vorticity term is optionally offered (ENE scheme) : 
    170188\begin{equation*} 
    171189   - \zeta \;{\textbf{k}}\times {\textbf {U}}_h  
    … …  
    185203\end{equation*} 
    186204 
    187 This formulation does not conserve the enstrophy but the total kinetic  
    188 energy. It is also possible to mix the two formulations in order to conserve  
    189 enstrophy on the relative vorticity term and energy on the Coriolis term. 
     205This formulation does not conserve the enstrophy but it does conserve the  
     206total kinetic energy. It is also possible to mix the two formulations in order  
     207to conserve enstrophy on the relative vorticity term and energy on the  
     208Coriolis term. 
    190209\begin{flalign*} 
    191210&\int\limits_D - \textbf{U}_h \cdot   \left(  \zeta \;\textbf{k} \times \textbf{U}_h  \right)  \;  dv &&  \\ 
    … …  
    219238 = - \int_D \textbf{U}_h \cdot w \frac{\partial \textbf{U}_h} {\partial k}\;dv 
    220239\end{equation*} 
    221 Indeed, using successively \eqref{DOM_di_adj} ($i.e.$ the skew symmetry property of the $\delta$ operator) and the incompressibility, then again \eqref{DOM_di_adj}, then  
    222 the commutativity of operators $\overline {\,\cdot \,}$ and $\delta$, and finally \eqref{DOM_mi_adj} ($i.e.$ the  symmetry property of the $\overline {\,\cdot \,}$ operator) applied in the horizontal and vertical direction, it becomes: 
     240Indeed, using successively \eqref{DOM_di_adj} ($i.e.$ the skew symmetry  
     241property of the $\delta$ operator) and the incompressibility, then  
     242\eqref{DOM_di_adj} again, then the commutativity of operators  
     243$\overline {\,\cdot \,}$ and $\delta$, and finally \eqref{DOM_mi_adj}  
     244($i.e.$ the symmetry property of the $\overline {\,\cdot \,}$ operator)  
     245applied in the horizontal and vertical directions, it becomes: 
    223246\begin{flalign*} 
    224247&\int_D \textbf{U}_h \cdot \nabla_h \left( \frac{1}{2}\;{\textbf{U}_h}^2 \right)\;dv   &&&\\ 
    … …  
    252275   \biggl\{  \overline {e_{1T}\,e_{2T} \,w}^{\,i+1/2}\;2     
    253276         \overline {u}^{\,k+1/2}\; \delta_{k+1/2}         \left[ u \right]     %&&&  \\ 
    254     + \overline {e_{1T}\,e_{2T} \,w}^{\,j+1/2}\;2  \overline {v}^{\,k+1/2}\; \delta_{k+1/2}        \left[ v \right]  \; 
     277    + \overline {e_{1T}\,e_{2T} \,w}^{\,j+1/2}\;2  \overline {v}^{\,k+1/2}\; \delta_{k+1/2}  \left[ v \right]  \; 
    255278   \biggr\}  
    256279   &&\displaybreak[0] \\  
    … …  
    272295\end{flalign*} 
    273296 
    274 The main point here is that the satisfaction of this property links the choice of the discrete formulation of vertical advection and of horizontal gradient of KE. Choosing one imposes the other. For example KE can also be discretized as $1/2\,({\overline u^{\,i}}^2 + {\overline v^{\,j}}^2)$. This leads to the following expression for the vertical advection: 
     297The main point here is that the satisfaction of this property links the choice of  
     298the discrete formulation of the vertical advection and of the horizontal gradient  
     299of KE. Choosing one imposes the other. For example KE can also be discretized  
     300as $1/2\,({\overline u^{\,i}}^2 + {\overline v^{\,j}}^2)$. This leads to the following  
     301expression for the vertical advection: 
    275302\begin{equation*} 
    276303\frac{1} {e_3 }\; w\; \frac{\partial \textbf{U}_h } {\partial k} 
    … …  
    284311\end{array}} } \right) 
    285312\end{equation*} 
    286 a formulation that requires a additional horizontal mean compare to the one used in NEMO. Nine velocity points have to be used instead of 3. This is the reason why it has not been choosen. 
     313a formulation that requires an additional horizontal mean in contrast with  
     314the one used in NEMO. Nine velocity points have to be used instead of 3.  
     315This is the reason why it has not been chosen. 
    287316 
    288317% ------------------------------------------------------------------------------------------------------------- 
    … …  
    298327\label{Apdx_C.1.3.1}  
    299328 
    300 In flux from the vorticity term reduces to a Coriolis term in which the Coriolis parameter has been modified to account for the ``metric'' term. This altered Coriolis parameter is discretised at F-point. It is given by:  
     329In flux from the vorticity term reduces to a Coriolis term in which the Coriolis  
     330parameter has been modified to account for the ``metric'' term. This altered  
     331Coriolis parameter is discretised at an f-point. It is given by:  
    301332\begin{equation*} 
    302333f+\frac{1} {e_1 e_2 }  
    … …  
    308339\end{equation*} 
    309340 
    310 The ENE scheme is then applied to obtain the vorticity term in flux form. It therefore conserves the total KE. The demonstration is same as for the vorticity term in vector invariant form (\S\ref{Apdx_C_vor}). 
     341The ENE scheme is then applied to obtain the vorticity term in flux form.  
     342It therefore conserves the total KE. The derivation is the same as for the  
     343vorticity term in the vector invariant form (\S\ref{Apdx_C_vor}). 
    311344 
    312345% ------------------------------------------------------------------------------------------------------------- 
    … …  
    316349\label{Apdx_C.1.3.2}  
    317350 
    318 The flux form operator of the momentum advection is evaluated using a centered second order finite difference scheme. Because of the flux form, the discrete operator does not contribute to the global budget of linear momentum. Because of the centered second order scheme, it conserves the horizontal kinetic energy, that is : 
     351The flux form operator of the momentum advection is evaluated using a  
     352centered second order finite difference scheme. Because of the flux form,  
     353the discrete operator does not contribute to the global budget of linear  
     354momentum. Because of the centered second order scheme, it conserves  
     355the horizontal kinetic energy, that is : 
    319356 
    320357\begin{equation} \label{Apdx_C_I.3.10} 
    … …  
    326363\end{equation} 
    327364 
    328 Let us demonstrate this property for the first term of the scalar product (i.e. considering just the the terms associated with the i-component of the advection): 
     365Let us demonstrate this property for the first term of the scalar product  
     366($i.e.$ considering just the the terms associated with the i-component of  
     367the advection): 
    329368\begin{flalign*} 
    330369&\int_D u \cdot \nabla \cdot \left(   \textbf{U}\,u   \right) \; dv    &&&\\ 
    … …  
    382421\end{flalign*} 
    383422 
    384 When the UBS scheme is used to evaluate the flux form momentum advection, the discrete operator does not contribute to the global budget of linear momentum (flux form). The horizontal kinetic energy is not conserved, but forced to decrease (the scheme is diffusive).  
     423When the UBS scheme is used to evaluate the flux form momentum advection,  
     424the discrete operator does not contribute to the global budget of linear momentum  
     425(flux form). The horizontal kinetic energy is not conserved, but forced to decay  
     426($i.e.$ the scheme is diffusive).  
    385427 
    386428% ------------------------------------------------------------------------------------------------------------- 
    … …  
    391433 
    392434 
    393 A pressure gradient has no contribution to the evolution of the vorticity as the curl of a gradient is zero. In $z$-coordinate, this properties is satisfied locally on a C-grid with 2nd order finite differences (property \eqref{Eq_DOM_curl_grad}). When the equation of state is linear ($i.e.$ when an advective-diffusive equation for density can be derived from those of temperature and salinity) the change of KE due to the work of pressure forces is balanced by the change of potential energy due to buoyancy forces:  
     435A pressure gradient has no contribution to the evolution of the vorticity as the  
     436curl of a gradient is zero. In the $z$-coordinate, this property is satisfied locally  
     437on a C-grid with 2nd order finite differences (property \eqref{Eq_DOM_curl_grad}).  
     438When the equation of state is linear ($i.e.$ when an advection-diffusion equation  
     439for density can be derived from those of temperature and salinity) the change of  
     440KE due to the work of pressure forces is balanced by the change of potential  
     441energy due to buoyancy forces:  
    394442\begin{equation*} 
    395443\int_D - \frac{1} {\rho_o} \left. \nabla p^h \right|_z \cdot \textbf{U}_h \;dv  
    … …  
    397445\end{equation*} 
    398446 
    399 This property can be satisfied in discrete sense for both $z$- and $s$-coordinates. Indeed, defining the depth of a $T$-point, $z_T$ defined as the sum of the vertical scale factors at $w$-points starting from the surface, the workof pressure forces can be written as: 
     447This property can be satisfied in a discrete sense for both $z$- and $s$-coordinates.  
     448Indeed, defining the depth of a $T$-point, $z_T$, as the sum of the vertical scale  
     449factors at $w$-points starting from the surface, the work of pressure forces can be  
     450written as: 
    400451\begin{flalign*} 
    401452&\int_D - \frac{1} {\rho_o} \left. \nabla p^h \right|_z \cdot \textbf{U}_h \;dv   &&& \\ 
    … …  
    410461\end{flalign*} 
    411462 
    412 Using  \eqref{DOM_di_adj}, $i.e.$ the skew symmetry property of the $\delta$ operator, \eqref{Eq_wzv}, the continuity equation), and \eqref{Eq_dynhpg_sco}, the hydrostatic  
    413 equation in $s$-coordinate, it turns out to be:  
     463Using  \eqref{DOM_di_adj}, $i.e.$ the skew symmetry property of the $\delta$  
     464operator, \eqref{Eq_wzv}, the continuity equation), and \eqref{Eq_dynhpg_sco},  
     465the hydrostatic equation in the $s$-coordinate, it becomes:  
    414466\begin{flalign*}  
    415467\equiv& \frac{1} {\rho_o} \sum\limits_{i,j,k}    \biggl\{  
    … …  
    467519\end{flalign*} 
    468520 
    469 Note that this property strongly constraints the discrete expression of both  
    470 the depth of $T-$points and of the term added to the pressure gradient in  
    471 $s$-coordinate. Nevertheless, it is almost never satisfied as a linear equation of state  
    472 is rarely used. 
     521Note that this property strongly constrains the discrete expression of both  
     522the depth of $T-$points and of the term added to the pressure gradient in the 
     523$s$-coordinate. Nevertheless, it is almost never satisfied since a linear equation  
     524of state is rarely used. 
    473525 
    474526% ------------------------------------------------------------------------------------------------------------- 
    … …  
    479531 
    480532 
    481 The surface pressure gradient has no contribution to the evolution of the vorticity. This property is trivially satisfied locally as the equation verified by $\psi $ has been derived from the discrete formulation of the momentum equation and of the curl. But it has to be noticed that since the elliptic equation verified by $\psi $ is solved numerically by an iterative solver (preconditioned conjugate gradient or successive over relaxation), the  
    482 property is only satisfied at the precision required on the solver used. 
    483  
    484 With the rigid-lid approximation, the change of KE due to the work of surface pressure forces is exactly zero. This is satisfied in discrete form, at the precision required on the elliptic solver used to solve this equation. This can be demonstrated as follows: 
     533The surface pressure gradient has no contribution to the evolution of the vorticity.  
     534This property is trivially satisfied locally since the equation verified by $\psi$ has  
     535been derived from the discrete formulation of the momentum equation and of the curl.  
     536But it has to be noted that since the elliptic equation satisfied by $\psi$ is solved  
     537numerically by an iterative solver (preconditioned conjugate gradient or successive  
     538over relaxation), the property is only satisfied at the precision requested for the  
     539solver used. 
     540 
     541With the rigid-lid approximation, the change of KE due to the work of surface  
     542pressure forces is exactly zero. This is satisfied in discrete form, at the precision  
     543requested for the elliptic solver used to solve this equation. This can be  
     544demonstrated as follows: 
    485545\begin{flalign*} 
    486546\int\limits_D  - \frac{1} {\rho_o} \nabla_h \left( p_s \right) \cdot \textbf{U}_h \;dv%   &&& \\ 
    … …  
    503563   && \\  
    504564% 
    505 \intertext{using the relation between \textit{$\psi $} and the vertically sum of the velocity, it becomes:} 
     565\intertext{using the relation between \textit{$\psi $} and the vertical sum of the velocity, it becomes:} 
    506566% 
    507567&\equiv \sum\limits_{i,j}  
    … …  
    537597\end{flalign*} 
    538598 
    539 The last equality is obtained using \eqref{Eq_dynspg_rl}, the discrete barotropic streamfunction time evolution equation. By the way, this shows that \eqref{Eq_dynspg_rl} is the only way do compute the streamfunction, otherwise the surface pressure forces will work. Nevertheless, since the elliptic equation verified by $\psi $ is solved numerically by an iterative solver, the property is only satisfied at the precision required on the solver. 
     599The last equality is obtained using \eqref{Eq_dynspg_rl}, the discrete barotropic  
     600streamfunction time evolution equation. By the way, this shows that  
     601\eqref{Eq_dynspg_rl} is the only way to compute the streamfunction,  
     602otherwise the surface pressure forces will do work. Nevertheless, since  
     603the elliptic equation satisfied by $\psi $ is solved numerically by an iterative  
     604solver, the property is only satisfied at the precision requested for the solver. 
    540605 
    541606% ================================================================ 
    … …  
    546611 
    547612 
    548 All the numerical schemes used in NEMO are written such that the tracer content is conserved by the internal dynamics and physics (equations in flux form). For advection, only the CEN2 scheme ($i.e.$ 2nd order finite different scheme) conserves the global variance of tracer. Nevertheless the other schemes ensure that the global variance decreases ($i.e.$ they are at least slightly diffusive). For diffusion, all the schemes ensure the decrease of the total tracer variance, but the iso-neutral operator. There is generally no strict conservation of mass, as the equation of state is non linear with respect to $T$ and $S$. In practice, the mass is conserved with a very good accuracy.  
     613All the numerical schemes used in NEMO are written such that the tracer content  
     614is conserved by the internal dynamics and physics (equations in flux form).  
     615For advection, only the CEN2 scheme ($i.e.$ $2^{nd}$ order finite different scheme)  
     616conserves the global variance of tracer. Nevertheless the other schemes ensure  
     617that the global variance decreases ($i.e.$ they are at least slightly diffusive).  
     618For diffusion, all the schemes ensure the decrease of the total tracer variance,  
     619except the iso-neutral operator. There is generally no strict conservation of mass,  
     620as the equation of state is non linear with respect to $T$ and $S$. In practice,  
     621the mass is conserved to a very high accuracy.  
    549622% ------------------------------------------------------------------------------------------------------------- 
    550623%       Advection Term 
    … …  
    553626\label{Apdx_C.2.1} 
    554627 
    555 Whatever the advection scheme considered it conserves of the tracer content as all the scheme are written in flux form. Let $\tau$ be the tracer interpolated at velocity point (whatever the interpolation is). The conservation of the tracer content is obtained as follows:  
     628Whatever the advection scheme considered it conserves of the tracer content as all  
     629the scheme are written in flux form. Let $\tau$ be the tracer interpolated at velocity point  
     630(whatever the interpolation is). The conservation of the tracer content is obtained as follows:  
    556631\begin{flalign*} 
    557632&\int_D \nabla \cdot \left( T \textbf{U} \right)\;dv    &&&\\ 
    … …  
    571646\end{flalign*} 
    572647 
    573 The conservation of the variance of tracer can be achieved only with the CEN2 scheme. It can be demonstarted as follows: 
     648The conservation of the variance of tracer can be achieved only with the CEN2 scheme.  
     649It can be demonstarted as follows: 
    574650\begin{flalign*} 
    575651&\int\limits_D T\;\nabla \cdot \left( T\; \textbf{U} \right)\;dv &&&\\ 
    … …  
    614690 
    615691The discrete formulation of the horizontal diffusion of momentum ensures the  
    616 conservation of potential vorticity and horizontal divergence and the  
     692conservation of potential vorticity and the horizontal divergence, and the  
    617693dissipation of the square of these quantities (i.e. enstrophy and the  
    618694variance of the horizontal divergence) as well as the dissipation of the  
    … …  
    623699horizontal divergence) and \textit{vice versa}.  
    624700 
    625 These properties of the horizontal diffusive operator are a direct  
    626 consequence of properties \eqref{Eq_DOM_curl_grad} and \eqref{Eq_DOM_div_curl}. When the vertical curl of the horizontal diffusion of momentum (discrete sense) is taken, the term associated to the horizontal gradient of the divergence is zero locally.  
     701These properties of the horizontal diffusion operator are a direct consequence  
     702of properties \eqref{Eq_DOM_curl_grad} and \eqref{Eq_DOM_div_curl}.  
     703When the vertical curl of the horizontal diffusion of momentum (discrete sense)  
     704is taken, the term associated with the horizontal gradient of the divergence is  
     705locally zero.  
    627706 
    628707% ------------------------------------------------------------------------------------------------------------- 
    … …  
    789868 
    790869When the horizontal divergence of the horizontal diffusion of momentum  
    791 (discrete sense) is taken, the term associated to the vertical curl of the  
    792 vorticity is zero locally, due to (II.1.8). The resulting term conserves the  
     870(discrete sense) is taken, the term associated with the vertical curl of the  
     871vorticity is zero locally, due to (!!! II.1.8  !!!!!). The resulting term conserves the  
    793872$\chi$ and dissipates $\chi^2$ when the eddy coefficients are  
    794873horizontally uniform. 
    … …  
    873952 
    874953 
    875 As for the lateral momentum physics, the continuous form of the vertical diffusion of momentum satisfies the several integral constraints. The first two are associated to the conservation of momentum and the dissipation of horizontal kinetic energy: 
     954As for the lateral momentum physics, the continuous form of the vertical diffusion  
     955of momentum satisfies several integral constraints. The first two are associated  
     956with the conservation of momentum and the dissipation of horizontal kinetic energy: 
    876957\begin{align*} 
    877958 \int\limits_D   
    … …  
    9231004   &&&\\  
    9241005% 
    925 \intertext{as the horizontal scale factor do not depend on $k$, it follows:} 
     1006\intertext{since the horizontal scale factor does not depend on $k$, it follows:} 
    9261007% 
    9271008&\equiv - \sum\limits_{i,j,k}  
    … …  
    9821063\end{flalign*} 
    9831064If the vertical diffusion coefficient is uniform over the whole domain, the  
    984 enstrophy is dissipated, i.e. 
     1065enstrophy is dissipated, $i.e.$ 
    9851066\begin{flalign*} 
    9861067\int\limits_D \zeta \, \textbf{k} \cdot \nabla \times  
    … …  
    10531134   &&\\  
    10541135\end{flalign*} 
    1055 Using the fact that the vertical diffusive coefficients are uniform and that in $z$-coordinates, the vertical scale factors do not depends on $i$ and $j$ so that: $e_{3f} =e_{3u} =e_{3v} =e_{3T} $ and $e_{3w} =e_{3uw} =e_{3vw} $, it follows: 
     1136Using the fact that the vertical diffusion coefficients are uniform, and that in  
     1137$z$-coordinate, the vertical scale factors do not depend on $i$ and $j$ so  
     1138that: $e_{3f} =e_{3u} =e_{3v} =e_{3T} $ and $e_{3w} =e_{3uw} =e_{3vw} $,  
     1139it follows: 
    10561140\begin{flalign*} 
    10571141\equiv A^{\,vm} \sum\limits_{i,j,k} \zeta \;\delta_k  
    … …  
    10901174   &&&\\ 
    10911175\end{flalign*} 
    1092 and the square of the horizontal divergence decreases (i.e. the horizontal divergence is dissipated) if vertical diffusion coefficient is uniform over the whole domain: 
     1176and the square of the horizontal divergence decreases ($i.e.$ the horizontal  
     1177divergence is dissipated) if the vertical diffusion coefficient is uniform over the  
     1178whole domain: 
    10931179 
    10941180\begin{flalign*} 
    … …  
    11031189   &&&\\ 
    11041190\end{flalign*} 
    1105 This property is only satisfied in $z$-coordinates: 
    1106  
     1191This property is only satisfied in the $z$-coordinate: 
    11071192\begin{flalign*} 
    11081193\int\limits_D \chi \;\nabla \cdot  
    … …  
    12091294\label{Apdx_C.5} 
    12101295 
    1211  
    1212  
    1213 The numerical schemes used for tracer subgridscale physics are written such that the heat and salt contents are conserved (equations in flux form, second order centered finite differences). As a form flux is used to compute the temperature and salinity, the quadratic form of these quantities (i.e. their variance) globally tends to diminish. As for the advection term, there is generally no strict conservation of mass even if, in practice, the mass is conserved with a very good accuracy.  
     1296The numerical schemes used for tracer subgridscale physics are written such  
     1297that the heat and salt contents are conserved (equations in flux form, second  
     1298order centered finite differences). Since a flux form is used to compute the  
     1299temperature and salinity, the quadratic form of these quantities (i.e. their variance)  
     1300globally tends to diminish. As for the advection term, there is generally no strict  
     1301conservation of mass, even if in practice the mass is conserved to a very high  
     1302accuracy.  
    12141303 
    12151304% ------------------------------------------------------------------------------------------------------------- 
    … …  
    12451334\end{flalign*} 
    12461335 
    1247 In fact, this property is simply resulting from the flux form of the operator.  
     1336In fact, this property simply results from the flux form of the operator.  
    12481337 
    12491338% ------------------------------------------------------------------------------------------------------------- 
    … …  
    12531342\label{Apdx_C.5.2} 
    12541343 
    1255 constraint of dissipation of tracer variance: 
     1344constraint on the dissipation of tracer variance: 
    12561345\begin{flalign*} 
    12571346\int\limits_D T\;\nabla & \cdot \left( A\;\nabla T \right)\;dv &&&\\  

  • trunk/DOC/TexFiles/Chapters/Annex_D.tex

    r994 r1223  
    77 
    88 
    9 A "model life" is more than ten years. Its software, composed by a few  
    10 hundred modules, is used by many people who are scientists or students and  
    11 do not necessary know very well all computer aspects. Moreover, a well  
    12 thought-out programme is easy to read and understand, less difficult to  
    13 modify, produces fewer bugs and is easier to maintain. Therefore, it is  
    14 essential that the model development follows some rules : 
     9A "model life" is more than ten years. Its software, composed of a few  
     10hundred modules, is used by many people who are scientists or students  
     11and do not necessarily know every aspect of computing very well.  
     12Moreover, a well thought-out program is easier to read and understand,  
     13less difficult to modify, produces fewer bugs and is easier to maintain.  
     14Therefore, it is essential that the model development follows some rules : 
    1515 
    1616- well planned and designed 
    … …  
    2727 
    2828To satisfy part of these aims, \NEMO is written with a coding standard which  
    29 is close to the ECMWF rules, named DOCTOR \citep{Gibson_TR86}. These rules present some advantages like : 
     29is close to the ECMWF rules, named DOCTOR \citep{Gibson_TR86}.  
     30These rules present some advantages like : 
    3031 
    3132- to provide a well presented program 
    … …  
    5152- the detail of input and output interfaces 
    5253 
    53 - the external routines and functions used (if exist) 
     54- the external routines and functions used (if they exist) 
    5455 
    55 - references (if exist) 
     56- references (if they exist) 
    5657 
    57 - the author name (s), the date of creation and of updates. 
     58- the author name(s), the date of creation and any updates. 
    5859 
    59 - Each program is splitted into several well separated sections and  
    60 sub-sections with a underlined title and specific labelled statements. 
     60- Each program is split into several well separated sections and  
     61sub-sections with an underlined title and specific labelled statements. 
    6162 
    6263- A program has not more than 200 to 300 lines. 
    … …  
    6869\label{Apdx_D_coding} 
    6970 
    70 - Use of the universal language \textsc{Fortran} 5 ANSI 77, with non  
    71 standard extensions, NAMELIST and \textsc{Fortran} 90 (matrix resolution  
    72 algorithm), and some well-identified particular statements or functions for  
    73 weak and massive parallelism, and vectorization  
     71- Use of the universal language \textsc{Fortran} 90, and try to avoid obsolescent 
     72features like statement functions, do not use GO TO and EQUIVALENCE statements. 
    7473 
    75 - A comment line begins with a uppercase character C at the first column. A  
    76 space line must have a C. For the on-line documentation, comments are  
    77 classified into three levels : 
     74- A continuation line begins with the character {\$} indented by three spaces  
     75compared to the previous line. 
    7876 
    79 Overview, triggered by CCC in columns 1 to 3. Only the title and the purpose  
    80 of the program are identified like that. This overview documentation can be  
    81 extracted by the UNIX function : grep -e '\^{}CCC' * 
     77- Always use a three spaces indentation in DO loop, CASE, or IF-ELSEIF-ELSE-ENDIF  
     78statements. 
    8279 
    83 External, triggered by CC in columns 1 to 2, and which correspond to  
    84 headlines of each programme, extracted by : grep -e '\^{}CC' * 
    85  
    86 Internal which are all the comments, extracted by : grep -e '\^{}C' * 
    87  
    88 - Statements GO TO, EQUIVALENCE are forbidden. 
    89  
    90 - A section is numbered with labels which are in agreement with the  
    91 paragraph label and increase from the begin to the end of routine. Labels of  
    92 a hundred ( 200, 201.. 220..) are reserved to a unique section. The  
    93 \textsc{Fortran} 90 extension syntax DO/ENDDO is used except for multitasked  
    94 do-loop. In this case labels 1000, 2000, ... are used. The FORMAT statement  
    95 are labelled with numbers in the range 9000 to 9999. 
    96  
    97 - A continuation line begins with the character {\$} in column 6. 
    98  
    99 - All statements begin in column 7 with the following gaps : 
    100  
    101 2 spaces toward the right in a DO loop. 
    102  
    103 4 spaces toward the right in IF, ELSEIF, ELSE and ENDIF statements, with  
    104 only 2 spaces for ELSE and ELSEIF lines. All IF statement must be followed  
    105 by a ELSE statement. 
    106  
    107 Some spaces in the continuation line for alignment. 
    108  
    109 - Use of different labels for each DO loop statement. 
    110  
    111 - STOP must be well documented with the name of the subroutine or a number. 
     80- use call to ctl\_stop routine instead of just a STOP. 
    11281 
    11382% ================================================================ 
    … …  
    11887 
    11988The purpose of the naming conventions is to use prefix letters to classify  
    120 model variables. These conventions allow to know easily the variable type  
    121 and to identify them rapidly: 
    122  
     89model variables. These conventions allow the variable type to be easily known 
     90and rapidly identified: 
    12391 
    12492%--------------------------------------------------TABLE-------------------------------------------------- 
    … …  
    12896\hline  Type \par / Status &   integer&   real&   logical &   character&   double \par precision&   complex \\   
    12997\hline 
    130 public &  
    131 \textbf{m n} \par \textit{but not } \par \textbf{nam}&  
    132 \textbf{a b e f g h o} \textbf{q} \textit{to} \textbf{x} \par but not \par \textbf{sf}&  
    133 \textbf{l} \par \textit{but not} \par \textbf{lp ld ll}&  
    134 \textbf{c} \par \textit{but not} \par \textbf{cp cd cl } \par \textbf{com cim}&  
    135 \textbf{d} \par \textit{but not} \par \textbf{dp dd dl}&  
     98public  \par or  \par module variable&  
     99\textbf{m n} \par \textit{but not } \par \textbf{nn\_}&  
     100\textbf{a b e f g h o} \textbf{q} \textit{to} \textbf{x} \par but not \par \textbf{fs rn\_}&  
     101\textbf{l} \par \textit{but not} \par \textbf{lp ld ll ln\_}&  
     102\textbf{c} \par \textit{but not} \par \textbf{cp cd cl cn\_}&  
     103\textbf{d} \par \textit{but not} \par \textbf{dp dd dl dn\_}&  
    136104\textbf{y} \par \textit{but not} \par \textbf{yp yd yl} \\ 
    137105\hline 
    … …  
    168136\textbf{yp} \\ 
    169137\hline 
    170 statement \par function&  
     138 
     139namelist& 
     140\textbf{nn\_}&  
     141\textbf{rn\_}&  
     142\textbf{ln\_}&  
     143\textbf{cn\_}&  
     144\textbf{dn\_}&  
     145\\ 
     146\hline 
     147CPP \par macro&  
    171148\textbf{kf}&  
    172149\textbf{sf} \par &  

  • trunk/DOC/TexFiles/Chapters/Annex_E.tex

    r996 r1223  
    2929- \frac{1}{2}\, |U|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i] 
    3030\end{equation} 
    31 where $U_{i+1/2} = e_{1u}\,e_{3u}\,u_{i+1/2}$ and $\tau "_i =\delta _i \left[ {\delta _{i+1/2} \left[ \tau \right]} \right]$. By choosing this expression for $\tau "$ we consider a fourth order approximation of $\partial_i^2$ with a constant i-grid spacing ($\Delta i=1$).  
    32  
    33 Alternative choice: introduce the scale factors:  $\tau "_i =\frac{e_{1T}}{e_{2T}\,e_{3T}}\delta _i \left[ \frac{e_{2u} e_{3u} }{e_{1u} }\delta _{i+1/2}[\tau] \right]$. 
     31where $U_{i+1/2} = e_{1u}\,e_{3u}\,u_{i+1/2}$ and  
     32$\tau "_i =\delta _i \left[ {\delta _{i+1/2} \left[ \tau \right]} \right]$.  
     33By choosing this expression for $\tau "$ we consider a fourth order approximation  
     34of $\partial_i^2$ with a constant i-grid spacing ($\Delta i=1$).  
     35 
     36Alternative choice: introduce the scale factors:   
     37$\tau "_i =\frac{e_{1T}}{e_{2T}\,e_{3T}}\delta _i \left[ \frac{e_{2u} e_{3u} }{e_{1u} }\delta _{i+1/2}[\tau] \right]$. 
    3438 
    3539 
    … …  
    4852scheme when \np{ln\_traadv\_ubs}=T. 
    4953 
    50 For stability reasons, in \eqref{Eq_tra_adv_ubs}, the first term which corresponds to a  
    51 second order centred scheme is evaluated using the \textit{now} velocity (centred in  
    52 time) while the second term which is the diffusive part of the scheme, is  
    53 evaluated using the \textit{before} velocity (forward in time. This is discussed by \citet{Webb1998} in the context of the Quick advection scheme. UBS and QUICK  
     54For stability reasons, in \eqref{Eq_tra_adv_ubs}, the first term which corresponds  
     55to a second order centred scheme is evaluated using the \textit{now} velocity  
     56(centred in time) while the second term which is the diffusive part of the scheme,  
     57is evaluated using the \textit{before} velocity (forward in time. This is discussed  
     58by \citet{Webb1998} in the context of the Quick advection scheme. UBS and QUICK  
    5459schemes only differ by one coefficient. Substituting 1/6 with 1/8 in  
    55 (\ref{Eq_tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb1998}. This  
    56 option is not available through a namelist parameter, since the 1/6  
     60(\ref{Eq_tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb1998}.  
     61This option is not available through a namelist parameter, since the 1/6  
    5762coefficient is hard coded. Nevertheless it is quite easy to make the  
    5863substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme 
    … …  
    6772NB 2 : In a forthcoming release four options will be proposed for the  
    6873vertical component used in the UBS scheme. $\tau _w^{ubs}$ will be  
    69 evaluated using either \textit{(a)} a centred $2^{nd}$ order scheme , or  \textit{(b)} a TVD  
    70 scheme, or  \textit{(c)} an interpolation based on conservative parabolic splines  
    71 following \citet{Sacha2005} implementation of UBS in ROMS,  
     74evaluated using either \textit{(a)} a centred $2^{nd}$ order scheme ,  
     75or  \textit{(b)} a TVD scheme, or  \textit{(c)} an interpolation based on conservative  
     76parabolic splines following \citet{Sacha2005} implementation of UBS in ROMS,  
    7277or  \textit{(d)} an UBS. The $3^{rd}$ case has dispersion properties similar to an  
    7378eight-order accurate conventional scheme. 
    … …  
    8893\end{split} 
    8994\end{equation} 
    90 \eqref{Eq_tra_adv_ubs2} has several advantages. First it clearly evidence that the UBS scheme is based on the fourth order scheme to which is added an upstream biased diffusive term. Second, this emphasises that the $4^{th}$ order part have to be evaluated at \emph{now} time step, not only the $2^{th}$ order part as stated above using \eqref{Eq_tra_adv_ubs}. Third, the diffusive term is in fact a biharmonic operator with a eddy coefficient with is simply proportional to the velocity. 
    91  
     95\eqref{Eq_tra_adv_ubs2} has several advantages. First it clearly evidence that  
     96the UBS scheme is based on the fourth order scheme to which is added an  
     97upstream biased diffusive term. Second, this emphasises that the $4^{th}$ order  
     98part have to be evaluated at \emph{now} time step, not only the $2^{th}$ order  
     99part as stated above using \eqref{Eq_tra_adv_ubs}. Third, the diffusive term is  
     100in fact a biharmonic operator with a eddy coefficient with is simply proportional  
     101to the velocity. 
    92102 
    93103laplacian diffusion: 
    … …  
    179189\end{align} 
    180190\end{subequations} 
    181 As for space operator, the adjoint of the derivation and averaging time operators are $\delta_t^*=\delta_{t+\Delta t/2}$ and $\overline{\cdot}^{\,t\,*}= \overline{\cdot}^{\,t+\Delta/2}$, respectively.  
     191As for space operator, the adjoint of the derivation and averaging time operators are  
     192$\delta_t^*=\delta_{t+\Delta t/2}$ and $\overline{\cdot}^{\,t\,*}= \overline{\cdot}^{\,t+\Delta/2}$ 
     193, respectively.  
    182194 
    183195The Leap-frog time stepping given by \eqref{Eq_DOM_nxt} can be defined as: 
    … …  
    187199      =         \frac{q^{t+\Delta t}-q^{t-\Delta t}}{2\Delta t} 
    188200\end{equation}  
    189 Note that \eqref{LF} shows that the leapfrog time step is $\Delta t$, not $2\Delta t$ as it can be found sometime in literature.  
    190 The leap-Frog time stepping is a second order centered scheme. As such it respects the quadratic invariant in integral forms, $i.e.$ the following continuous property, 
     201Note that \eqref{LF} shows that the leapfrog time step is $\Delta t$, not $2\Delta t$  
     202as it can be found sometime in literature.  
     203The leap-Frog time stepping is a second order centered scheme. As such it respects  
     204the quadratic invariant in integral forms, $i.e.$ the following continuous property, 
    191205\begin{equation} \label{Energy} 
    192206\int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt}  
    … …  
    205219      \equiv \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right) 
    206220\end{split} \end{equation} 
    207 NB here pb of boundary condition when applying the adjoin! In space, setting to 0 the quantity in land area is sufficient to get rid of the boundary condition (equivalently of the boundary value of the integration by part). In time this boundary condition is not physical and \textbf{add something here!!!} 
    208  
    209  
    210  
     221NB here pb of boundary condition when applying the adjoin! In space, setting to 0  
     222the quantity in land area is sufficient to get rid of the boundary condition  
     223(equivalently of the boundary value of the integration by part). In time this boundary  
     224condition is not physical and \textbf{add something here!!!} 
     225 
     226 
     227 

  • trunk/DOC/TexFiles/Chapters/Chap_Conservation.tex

    r996 r1223  
    33% Invariant of the Equations 
    44% ================================================================ 
    5 \chapter{Annex --- Invariants of the Primitive Equations} 
     5\chapter{Invariants of the Primitive Equations} 
    66\label{Invariant} 
    77\minitoc 
Note: See TracChangeset for help on using the changeset viewer.