- Timestamp:
- 2010-10-15T16:42:00+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_Model_Basics.tex
r1224 r2282 7 7 \minitoc 8 8 9 10 \newpage 11 $\ $\newline % force a new ligne 9 12 10 13 % ================================================================ … … 164 167 165 168 169 \newpage 170 $\ $\newline % force a new ligne 171 166 172 % ================================================================ 167 173 % The Horizontal Pressure Gradient … … 196 202 the sea surface is the surface $z=0$. This well known approximation increases the surface 197 203 wave speed to infinity and modifies certain other longwave dynamics ($e.g.$ barotropic 198 Rossby or planetary waves). In the present release of \NEMO, both strategies are still available. 199 They are further described in the next two sub-sections. 204 Rossby or planetary waves). The rigid-lid hypothesis is an obsolescent feature in modern 205 OGCMs. It has been available until the release 3.1 of \NEMO, and it has been removed 206 in release 3.2 and followings. Only the free surface formulation is now described in the 207 this document (see the next sub-section). 200 208 201 209 % ------------------------------------------------------------------------------------------------------------- … … 221 229 short with respect to the other processes described by the primitive equations. 222 230 223 T hreechoices can be made regarding the implementation of the free surface in the model,231 Two choices can be made regarding the implementation of the free surface in the model, 224 232 depending on the physical processes of interest. 225 233 … … 238 246 of heat and salt contents is sufficient for the problem solved, then it is 239 247 sufficient to solve a linearized version of (\ref{Eq_PE_ssh}), which still allows 240 to take into account freshwater fluxes applied at the ocean surface \citep{Roullet2000}. 241 242 $\bullet$ For process studies not involving external waves nor surface freshwater 243 fluxes, it is possible to use the rigid lid approximation see (next 244 section). The ocean surface is then considered as a fixed surface, so that all 245 external waves are removed from the system. 248 to take into account freshwater fluxes applied at the ocean surface \citep{Roullet_Madec_JGR00}. 246 249 247 250 The filtering of EGWs in models with a free surface is usually a matter of discretisation 248 of the temporal derivatives, using the time splitting method \citep{Killworth 1991, Zhang1992}251 of the temporal derivatives, using the time splitting method \citep{Killworth_al_JPO91, Zhang_Endoh_JGR92} 249 252 or the implicit scheme \citep{Dukowicz1994}. In \NEMO, we use a slightly different approach 250 developed by \citet{Roullet 2000}: the damping of EGWs is ensured by introducing an253 developed by \citet{Roullet_Madec_JGR00}: the damping of EGWs is ensured by introducing an 251 254 additional force in the momentum equation. \eqref{Eq_PE_dyn} becomes: 252 255 \begin{equation} \label{Eq_PE_flt} … … 266 269 $i.e.$ they are stationary and damped. The diffusion regime applies to the modes shorter than 267 270 $T_c$. For longer ones, the diffusion term vanishes. Hence, the temporally unresolved EGWs 268 can be damped by choosing $T_c > \ Delta t$. \citet{Roullet2000} demonstrate that271 can be damped by choosing $T_c > \rdt$. \citet{Roullet_Madec_JGR00} demonstrate that 269 272 (\ref{Eq_PE_flt}) can be integrated with a leap frog scheme except the additional term which 270 273 has to be computed implicitly. This is not surprising since the use of a large time step has a … … 272 275 Firstly, the damping of EGWs can be quantified through the magnitude of the additional term. 273 276 Secondly, the numerical scheme does not need any tuning. Numerical stability is ensured as 274 soon as $T_c > \ Deltat$.277 soon as $T_c > \rdt$. 275 278 276 279 When the variations of free surface elevation are small compared to the thickness of the first 277 280 model layer, the free surface equation (\ref{Eq_PE_ssh}) can be linearized. As emphasized 278 by \citet{Roullet 2000} the linearization of (\ref{Eq_PE_ssh}) has consequences on the281 by \citet{Roullet_Madec_JGR00} the linearization of (\ref{Eq_PE_ssh}) has consequences on the 279 282 conservation of salt in the model. With the nonlinear free surface equation, the time evolution 280 283 of the total salt content is … … 287 290 integral over the free surface) vanishes when the nonlinear equation (\ref{Eq_PE_ssh}) 288 291 is satisfied, so that the salt is perfectly conserved. When the free surface equation is 289 linearized, \citet{Roullet 2000} show that the total salt content integrated in the fixed292 linearized, \citet{Roullet_Madec_JGR00} show that the total salt content integrated in the fixed 290 293 volume $D$ (bounded by the surface $z=0$) is no longer conserved: 291 294 \begin{equation} \label{Eq_PE_salt_content_linear} … … 295 298 296 299 The right hand side of (\ref{Eq_PE_salt_content_linear}) is small in equilibrium solutions 297 \citep{Roullet 2000}. It can be significant when the freshwater forcing is not balanced and300 \citep{Roullet_Madec_JGR00}. It can be significant when the freshwater forcing is not balanced and 298 301 the globally averaged free surface is drifting. An increase in sea surface height \textit{$\eta $} 299 302 results in a decrease of the salinity in the fixed volume $D$. Even in that case though, … … 309 312 its variations are driven by correlations of the time variation of surface salinity with the 310 313 sea surface height, which is a negligible term. This situation contrasts with the case of 311 the rigid lid approximation (following section) in which case freshwater forcing is 312 represented by a virtual salt flux, leading to a spurious source of salt at the ocean 313 surface \citep{Roullet2000}. 314 315 % ------------------------------------------------------------------------------------------------------------- 316 % Rigid-Lid Formulation 317 % ------------------------------------------------------------------------------------------------------------- 318 \subsection{Rigid-Lid formulation} 319 \label{PE_rigid_lid} 320 321 With the rigid lid approximation, we assume that the ocean surface ($z=0$) is a rigid lid 322 on which a pressure $p_s$ is exerted. This implies that the vertical velocity at the surface 323 is equal to zero. From the continuity equation \eqref{Eq_PE_continuity} and the kinematic 324 condition at the bottom \eqref{Eq_PE_w_bbc} (no flux across the bottom), it can be shown 325 that the vertically integrated flow $H{\rm {\bf \bar {U}}}_h$ is non-divergent (where the 326 overbar indicates a vertical average over the whole water column, i.e. from $z=-H$, 327 the ocean bottom, to $z=0$ , the rigid-lid). Thus, $\rm {\bf \bar {U}}_h$ can be derived 328 from a volume transport streamfunction $\psi$: 329 \begin{equation} \label{Eq_PE_u_psi} 330 \overline{\vect{U}}_h =\frac{1}{H}\left( \vect{k} \times \nabla \psi \right) 331 \end{equation} 332 333 As $p_s$ does not depend on depth, its horizontal gradient is obtained by forming the 334 vertical average of \eqref{Eq_PE_dyn} and using \eqref{Eq_PE_u_psi}: 335 336 \begin{equation} \label{Eq_PE_u_barotrope} 337 \frac{1}{\rho _o }\nabla _h p_s 338 =\overline{\vect{M}} -\frac{\partial \overline{\vect{U}} _h }{\partial t} 339 =\overline{\vect{M}} 340 -\frac{1}{H} \left[ \vect{k} \times \nabla \left( \frac{\partial \psi}{\partial t} \right) \right] 341 \end{equation} 342 343 Here ${\rm {\bf M}} = \left( M_u,M_v \right)$ represents the collected contributions of the 344 Coriolis, hydrostatic pressure gradient, nonlinear and viscous terms in \eqref{Eq_PE_dyn}. 345 The time derivative of $\psi $ is the solution of an elliptic equation which is obtained from 346 the vertical component of the curl of (\ref{Eq_PE_u_barotrope}): 347 \begin{equation} \label{Eq_PE_psi} 348 \left[ {\nabla \times \left[ {\frac{1}{H} \vect{\bf k} 349 \times \nabla \left( {\frac{\partial \psi }{\partial t}} \right)} \right]} \; \right]_z 350 =\left[ {\nabla \times \overline{\vect{M}} } \right]_z 351 \end{equation} 352 353 Using the proper boundary conditions, \eqref{Eq_PE_psi} can be solved to find $\partial_t \psi$ 354 and thus using \eqref{Eq_PE_u_barotrope} the horizontal surface pressure gradient. 355 It should be noted that $p_s$ can be computed by taking the divergence of 356 \eqref{Eq_PE_u_barotrope} and solving the resulting elliptic equation. Thus the surface 357 pressure is a diagnostic quantity that can be recovered for analysis purposes. 358 359 A difficulty lies in the determination of the boundary condition on $\partial_t \psi$. 360 The boundary condition on velocity is that there is no flow normal to a solid wall, 361 $i.e.$ the coastlines are streamlines. Therefore \eqref{Eq_PE_psi} is solved with 362 the following Dirichlet boundary condition: $\partial_t \psi$ is constant along each 363 coastline of the same continent or of the same island. When all the coastlines are 364 connected (there are no islands), the constant value of $\partial_t \psi$ along the 365 coast can be arbitrarily chosen to be zero. When islands are present in the domain, 366 the value of the barotropic streamfunction will generally be different for each island 367 and for the continent, and will vary with respect to time. So the boundary condition is: 368 $\psi=0$ along the continent and $\psi=\mu_n$ along island $n$ ($1 \leq n \leq Q$), 369 where $Q$ is the number of islands present in the domain and $\mu_n$ is a time 370 dependent variable. A time evolution equation of the unknown $\mu_n$ can be found 371 by evaluating the circulation of the time derivative of the vertical average (barotropic) 372 velocity field along a closed contour around each island. Since the circulation of a 373 gradient field along a closed contour is zero, from \eqref{Eq_PE_u_barotrope} we have: 374 \begin{equation} \label{Eq_PE_isl_circulation} 375 \oint_n {\frac{1}{H}\left[ {{\rm {\bf k}}\times \nabla \left( 376 {\frac{\partial \psi }{\partial t}} \right)} \right] \cdot {\rm {\bf d}}\ell } 377 = \oint_n {\overline {\rm {\bf M}} \cdot {\rm {\bf d}}\ell } 378 \qquad 1 \leq n \leq Q 379 \end{equation} 380 381 Since (\ref{Eq_PE_psi}) is linear, its solution \textit{$\psi $} can be decomposed 382 as follows: 383 \begin{equation} \label{Eq_PE_psi_isl} 384 \psi =\psi _o +\sum\limits_{n=1}^{n=Q} {\mu _n \psi _n } 385 \end{equation} 386 where $\psi _o$ is the solution of \eqref{Eq_PE_psi} with $\psi _o=0$ long all 387 the coastlines, and where $\psi _n$ is the solution of \eqref{Eq_PE_psi} with 388 the right-hand side equal to $0$, and with $\psi _n =1$ long the island $n$, 389 $\psi _n =0$ along the other boundaries. The function $\psi _n$ is thus 390 independent of time. Introducing \eqref{Eq_PE_psi_isl} into 391 \eqref{Eq_PE_isl_circulation} yields: 392 \begin{multline} \label{Eq_PE_psi_isl_circulation} 393 \left[ {\oint_n {\frac{1}{H} \left[ {{\rm {\bf k}}\times \nabla \psi _m } \right]\cdot {\rm {\bf d}}\ell } } \right]_{1\leq m\leqslant Q \atop 1\leq n\leqslant Q } 394 \left( {\frac{\partial \mu _n }{\partial t}} 395 \right)_{1\leqslant n\leqslant Q} \\ 396 =\left( {\oint_n {\left[ {\overline {\rm 397 {\bf M}} -\frac{1}{H}\left[ {{\rm {\bf k}}\times \nabla \left( 398 {\frac{\partial \psi _o }{\partial t}} \right)} \right]} \right]\cdot {\rm 399 {\bf d}}\ell } } \right)_{1\leqslant n\leqslant Q} 400 \end{multline} 401 which can be rewritten as: 402 \begin{equation} \label{Eq_PE_psi_isl_matrix} 403 {\rm {\bf A}}\;\left( {\frac{\partial \mu _n }{\partial t}} 404 \right)_{1\leqslant n\leqslant Q} ={\rm {\bf B}} 405 \end{equation} 406 where \textbf{A} is a $Q \times Q$ matrix and \textbf{B} is a time dependent vector. 407 As \textbf{A} is independent of time, it can be calculated and inverted once. The time 408 derivative of the streamfunction when islands are present is thus given by: 409 \begin{equation} \label{Eq_PE_psi_isl_dt} 410 \frac{\partial \psi }{\partial t}=\frac{\partial \psi _o }{\partial 411 t}+\sum\limits_{n=1}^{n=Q} {{\rm {\bf A}}^{-1}{\rm {\bf B}}\;\psi _n } 412 \end{equation} 413 414 314 the rigid lid approximation in which case freshwater forcing is represented by a virtual 315 salt flux, leading to a spurious source of salt at the ocean surface 316 \citep{Huang_JPO93, Roullet_Madec_JGR00}. 317 318 \newpage 319 $\ $\newline % force a new ligne 415 320 416 321 % ================================================================ … … 435 340 cannot be easily treated in a global model without filtering. A solution consists of introducing 436 341 an appropriate coordinate transformation that shifts the singular point onto land 437 \citep{Madec Imb1996, Murray1996}. As a consequence, it is important to solve the primitive342 \citep{Madec_Imbard_CD96, Murray_JCP96}. As a consequence, it is important to solve the primitive 438 343 equations in various curvilinear coordinate systems. An efficient way of introducing an 439 344 appropriate coordinate transform can be found when using a tensorial formalism. … … 444 349 of the conservation laws of fluid dynamics. 445 350 446 Let (\text bf{i},\textbf{j},\textbf{k}) be a set of orthogonal curvilinear coordinates on the sphere351 Let (\textit{i},\textit{j},\textit{k}) be a set of orthogonal curvilinear coordinates on the sphere 447 352 associated with the positively oriented orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) 448 353 linked to the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are … … 682 587 term and can be viewed as a modification of the Coriolis parameter: 683 588 \begin{equation} \label{Eq_PE_cor+metric} 684 f \to f + \frac{1}{e_1 \; e_2} \left(v \frac{\partial e_2}{\partial i}685 -u \frac{\partial e_1}{\partial j}\right)589 f \to f + \frac{1}{e_1\;e_2} \left( v \frac{\partial e_2}{\partial i} 590 -u \frac{\partial e_1}{\partial j} \right) 686 591 \end{equation} 687 592 … … 690 595 the Coriolis parameter $f \to f+(u/a) \tan \varphi$. 691 596 692 To sum up, the equations solved by the ocean model can be written in the following tensorial formalism: 597 598 $\ $\newline % force a new ligne 599 600 To sum up, the curvilinear $z$-coordinate equations solved by the ocean model can be 601 written in the following tensorial formalism: 693 602 694 603 \vspace{+10pt} 695 $\bullet$ \textit{momentum equations} : 696 697 vector invariant form : 604 $\bullet$ \textbf{Vector invariant form of the momentum equations} : 605 698 606 \begin{subequations} \label{Eq_PE_dyn_vect} 699 \begin{ multline} \label{Eq_PE_dyn_vect_u}700 \frac{\partial u}{\partial t} =701 702 - \frac{1}{2\,e_1}\frac{\partial}{\partial i} \left( u^2+v^2 \right)703 - \frac{1}{e_3} w \frac{\partial u}{\partial k}\\704 - \frac{1}{e_1}\frac{\partial}{\partial i} \left( \frac{p_s+p_h }{\rho _o} \right)705 + D_u^{\vect{U}} + F_u^{\vect{U}} 706 \ end{multline}707 \ begin{multline} \label{Eq_PE_dyn_vect_v}708 \frac{\partial v}{\partial t}= 709 - \left( {\zeta +f} \right)\,u710 - \frac{1}{2\,e_2 }\frac{\partial }{\partial j}\left( u^2+v^2 \right) - \frac{1}{e_3 }w\frac{\partial v}{\partial k}\\711 - \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho _o} \right)712 + D_v^{\vect{U}} + F_v^{\vect{U}}713 \end{ multline}607 \begin{equation} \label{Eq_PE_dyn_vect_u} \begin{split} 608 \frac{\partial u}{\partial t} 609 = + \left( {\zeta +f} \right)\,v 610 - \frac{1}{2\,e_1} \frac{\partial}{\partial i} \left( u^2+v^2 \right) 611 - \frac{1}{e_3 } w \frac{\partial u}{\partial k} & \\ 612 - \frac{1}{e_1 } \frac{\partial}{\partial i} \left( \frac{p_s+p_h }{\rho _o} \right) 613 &+ D_u^{\vect{U}} + F_u^{\vect{U}} \\ 614 \\ 615 \frac{\partial v}{\partial t} = 616 - \left( {\zeta +f} \right)\,u 617 - \frac{1}{2\,e_2 } \frac{\partial }{\partial j}\left( u^2+v^2 \right) 618 - \frac{1}{e_3 } w \frac{\partial v}{\partial k} & \\ 619 - \frac{1}{e_2 } \frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho _o} \right) 620 &+ D_v^{\vect{U}} + F_v^{\vect{U}} 621 \end{split} \end{equation} 714 622 \end{subequations} 715 623 716 flux form: 624 625 \vspace{+10pt} 626 $\bullet$ \textbf{flux form of the momentum equations} : 717 627 \begin{subequations} \label{Eq_PE_dyn_flux} 718 628 \begin{multline} \label{Eq_PE_dyn_flux_u} … … 741 651 \end{multline} 742 652 \end{subequations} 743 where $\zeta$ is given by \eqref{Eq_PE_curl_Uh} and the surface pressure gradient formulation 744 depends on the one of the free surface: 745 746 $*$ free surface formulation 747 \begin{equation}\label{Eq_PE_dyn_sco} 748 \frac{1}{\rho _o }\nabla _h p_s =\left( {{\begin{array}{*{20}c} 749 {\frac{g}{\;e_1 }\frac{\partial \eta }{\partial i}} \hfill \\ 750 {\frac{g}{\;e_2 }\frac{\partial \eta }{\partial j}} \hfill \\ 751 \end{array} }} \right) 752 \qquad \text{where $\eta$ is solution of \eqref{Eq_PE_ssh} } 753 \end{equation} 754 755 $*$ rigid-lid approximation 756 \begin{equation}\label{Eq_PE_dyn_zco} 757 \frac{1}{\rho _o }\nabla _h p_s =\left( {{\begin{array}{*{20}c} 758 {\overline M _u +\frac{1}{H\;e_2 }\frac{\partial }{\partial j}\left( 759 {\frac{\partial \psi }{\partial t}} \right)} \\ 760 {\overline M _v -\frac{1}{H\;e_1 }\frac{\partial }{\partial i}\left( 761 {\frac{\partial \psi }{\partial t}} \right)} \\ 762 \end{array} }} \right) 763 \end{equation} 764 where ${\vect{M}}= \left( M_u,M_v \right)$ represents the collected contributions of nonlinear, 765 viscosity and hydrostatic pressure gradient terms in \eqref{Eq_PE_dyn_vect} and the overbar 766 indicates a vertical average over the whole water column ($i.e.$ from $z=-H$, the ocean bottom, 767 to $z=0$, the rigid-lid), and where the time derivative of $\psi$ is the solution of an elliptic equation: 768 \begin{multline} \label{Eq_psi_total} 769 \frac{\partial }{\partial i}\left[ {\frac{e_2 }{H\,e_1}\frac{\partial}{\partial i} 770 \left( {\frac{\partial \psi }{\partial t}} \right)} \right] 771 +\frac{\partial }{\partial j}\left[ {\frac{e_1 }{H\,e_2}\frac{\partial }{\partial j} 772 \left( {\frac{\partial \psi }{\partial t}} \right)} \right] 773 = \\ 774 \frac{\partial }{\partial i}\left( {e_2 \overline M _v } \right) 775 - \frac{\partial }{\partial j}\left( {e_1 \overline M _u } \right) 776 \end{multline} 653 where $\zeta$, the relative vorticity, is given by \eqref{Eq_PE_curl_Uh} and $p_s $, 654 the surface pressure, is given by: 655 \begin{equation} \label{Eq_PE_spg} 656 p_s = \left\{ \begin{split} 657 \rho \,g \,\eta & \qquad \qquad \; \qquad \text{ standard free surface} \\ 658 \rho \,g \,\eta &+ \rho_o \,\mu \,\frac{\partial \eta }{\partial t} \qquad \text{ filtered free surface} 659 \end{split} 660 \right. 661 \end{equation} 662 with $\eta$ is solution of \eqref{Eq_PE_ssh} 777 663 778 664 The vertical velocity and the hydrostatic pressure are diagnosed from the following equations: … … 783 669 \frac{\partial p_h }{\partial k}=-\rho \;g\;e_3 784 670 \end{equation} 785 786 671 where the divergence of the horizontal velocity, $\chi$ is given by \eqref{Eq_PE_div_Uh}. 787 672 … … 809 694 in Chapter~\ref{SBC}. 810 695 696 811 697 \newpage 698 $\ $\newline % force a new ligne 812 699 % ================================================================ 813 % Curvilinear z*-coordinate System700 % Curvilinear generalised vertical coordinate System 814 701 % ================================================================ 815 \section{Curvilinear \textit{z*}-coordinate System} 816 817 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 818 \begin{figure}[!b] \label{Fig_z_zstar} \begin{center} 819 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_z_zstar.pdf} 820 \caption{(a) $z$-coordinate in linear free-surface case ; (b) $z-$coordinate in non-linear 821 free surface case (c) re-scaled height coordinate (become popular as the \textit{z*-}coordinate 822 \citep{Adcroft_Campin_OM04} ).} 823 \end{center} \end{figure} 824 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 825 826 827 In that case, the free surface equation is nonlinear, and the variations of volume are fully 828 taken into account. These coordinates systems is presented in a report \citep{Levier2007} 829 available on the \NEMO web site. 830 831 \gmcomment{ 832 The \textit{z*} coordinate approach is an unapproximated, non-linear free surface implementation 833 which allows one to deal with large amplitude free-surface 834 variations relative to the vertical resolution \citep{Adcroft_Campin_OM04}. In 835 the \textit{z*} formulation, the variation of the column thickness due to sea-surface 836 undulations is not concentrated in the surface level, as in the $z$-coordinate formulation, 837 but is equally distributed over the full water column. Thus vertical 838 levels naturally follow sea-surface variations, with a linear attenuation with 839 depth, as illustrated by figure fig.1c . Note that with a flat bottom, such as in 840 fig.1c, the bottom-following $z$ coordinate and \textit{z*} are equivalent. 841 The definition and modified oceanic equations for the rescaled vertical coordinate 842 \textit{z*}, including the treatment of fresh-water flux at the surface, are 843 detailed in Adcroft and Campin (2004). The major points are summarized 844 here. The position ( \textit{z*}) and vertical discretization (\textit{z*}) are expressed as: 845 846 $H + \textit{z*} = (H + z) / r$ and $\delta \textit{z*} = \delta z / r$ with $r = \frac{H+\eta} {H}$ 847 848 Since the vertical displacement of the free surface is incorporated in the vertical 849 coordinate \textit{z*}, the upper and lower boundaries are at fixed \textit{z*} position, 850 $\textit{z*} = 0$ and $\textit{z*} = ?H$ respectively. Also the divergence of the flow field 851 is no longer zero as shown by the continuity equation: 852 853 $\frac{\partial r}{\partial t} = \nabla_{\textit{z*}} \cdot \left( r \; \rm{\bf U}_h \right) 854 \left( r \; w\textit{*} \right) = 0 $ 855 856 } 857 858 859 \newpage 860 % ================================================================ 861 % Curvilinear s-coordinate System 862 % ================================================================ 863 \section{Curvilinear \textit{s}-coordinate System} 864 \label{PE_sco} 865 866 % ------------------------------------------------------------------------------------------------------------- 867 % Introduction 868 % ------------------------------------------------------------------------------------------------------------- 869 \subsection{Introduction} 870 871 Several important aspects of the ocean circulation are influenced by bottom topography. 872 Of course, the most important is that bottom topography determines deep ocean sub-basins, 873 barriers, sills and channels that strongly constrain the path of water masses, but more subtle 874 effects exist. For example, the topographic $\beta$-effect is usually larger than the planetary 875 one along continental slopes. Topographic Rossby waves can be excited and can interact 876 with the mean current. In the $z-$coordinate system presented in the previous section 877 (\S\ref{PE_zco}), $z-$surfaces are geopotential surfaces. The bottom topography is 878 discretised by steps. This often leads to a misrepresentation of a gradually sloping bottom 879 and to large localized depth gradients associated with large localized vertical velocities. 880 The response to such a velocity field often leads to numerical dispersion effects. 881 One solution to strongly reduce this error is to use a partial step representation of bottom 882 topography instead of a full step one \cite{Pacanowski_Gnanadesikan_MWR98}. 883 Another solution is to introduce a terrain-following coordinate system (hereafter $s-$coordinate) 884 885 The $s$-coordinate avoids the discretisation error in the depth field since the layers of 886 computation are gradually adjusted with depth to the ocean bottom. Relatively small 887 topographic features as well as gentle, large-scale slopes of the sea floor in the deep 888 ocean, which would be ignored in typical $z$-model applications with the largest grid 889 spacing at greatest depths, can easily be represented (with relatively low vertical resolution). 890 A terrain-following model (hereafter $s-$model) also facilitates the modelling of the 891 boundary layer flows over a large depth range, which in the framework of the $z$-model 892 would require high vertical resolution over the whole depth range. Moreover, with a 893 $s$-coordinate it is possible, at least in principle, to have the bottom and the sea surface 894 as the only boundaries of the domain (nomore lateral boundary condition to specify). 895 Nevertheless, a $s$-coordinate also has its drawbacks. Perfectly adapted to a 896 homogeneous ocean, it has strong limitations as soon as stratification is introduced. 897 The main two problems come from the truncation error in the horizontal pressure 898 gradient and a possibly increased diapycnal diffusion. The horizontal pressure force 899 in $s$-coordinate consists of two terms (see Appendix~\ref{Apdx_A}), 900 901 \begin{equation} \label{Eq_PE_p_sco} 902 \left. {\nabla p} \right|_z =\left. {\nabla p} \right|_s -\frac{\partial 903 p}{\partial s}\left. {\nabla z} \right|_s 904 \end{equation} 905 906 The second term in \eqref{Eq_PE_p_sco} depends on the tilt of the coordinate surface 907 and introduces a truncation error that is not present in a $z$-model. In the special case 908 of a $\sigma-$coordinate (i.e. a depth-normalised coordinate system $\sigma = z/H$), 909 \citet{Haney1991} and \citet{Beckmann1993} have given estimates of the magnitude 910 of this truncation error. It depends on topographic slope, stratification, horizontal and 911 vertical resolution, the equation of state, and the finite difference scheme. This error 912 limits the possible topographic slopes that a model can handle at a given horizontal 913 and vertical resolution. This is a severe restriction for large-scale applications using 914 realistic bottom topography. The large-scale slopes require high horizontal resolution, 915 and the computational cost becomes prohibitive. This problem can be at least partially 916 overcome by mixing $s$-coordinate and step-like representation of bottom topography \citep{Gerdes1993a,Gerdes1993b,Madec1996}. However, the definition of the model 917 domain vertical coordinate becomes then a non-trivial thing for a realistic bottom 918 topography: a envelope topography is defined in $s$-coordinate on which a full or 919 partial step bottom topography is then applied in order to adjust the model depth to 920 the observed one (see \S\ref{DOM_zgr}. 921 922 For numerical reasons a minimum of diffusion is required along the coordinate surfaces 923 of any finite difference model. It causes spurious diapycnal mixing when coordinate 924 surfaces do not coincide with isoneutral surfaces. This is the case for a $z$-model as 925 well as for a $s$-model. However, density varies more strongly on $s-$surfaces than 926 on horizontal surfaces in regions of large topographic slopes, implying larger diapycnal 927 diffusion in a $s$-model than in a $z$-model. Whereas such a diapycnal diffusion in a 928 $z$-model tends to weaken horizontal density (pressure) gradients and thus the horizontal 929 circulation, it usually reinforces these gradients in a $s$-model, creating spurious circulation. 930 For example, imagine an isolated bump of topography in an ocean at rest with a horizontally 931 uniform stratification. Spurious diffusion along $s$-surfaces will induce a bump of isoneutral 932 surfaces over the topography, and thus will generate there a baroclinic eddy. In contrast, 933 the ocean will stay at rest in a $z$-model. As for the truncation error, the problem can be reduced by introducing the terrain-following coordinate below the strongly stratified portion of the water column 934 ($i.e.$ the main thermocline) \citep{Madec1996}. An alternate solution consists of rotating 935 the lateral diffusive tensor to geopotential or to isoneutral surfaces (see \S\ref{PE_ldf}. 936 Unfortunately, the slope of isoneutral surfaces relative to the $s$-surfaces can very large, 937 strongly exceeding the stability limit of such a operator when it is discretized (see Chapter~\ref{LDF}). 938 939 The $s-$coordinates introduced here \citep{Lott1990,Madec1996} differ mainly in two 940 aspects from similar models: it allows a representation of bottom topography with mixed 941 full or partial step-like/terrain following topography ; It also offers a completely general 942 transformation, $s=s(i,j,z)$ for the vertical coordinate. 702 \section{Curvilinear generalised vertical coordinate System} 703 \label{PE_gco} 704 705 %\gmcomment{ 706 The ocean domain presents a huge diversity of situation in the vertical. First the ocean surface is a time dependent surface (moving surface). Second the ocean floor depends on the geographical position, varying from more than 6,000 meters in abyssal trenches to zero at the coast. Last but not least, the ocean stratification exerts a strong barrier to vertical motions and mixing. 707 Therefore, in order to represent the ocean with respect to the first point a space and time dependent vertical coordinate that follows the variation of the sea surface height $e.g.$ an $z$*-coordinate; for the second point, a space variation to fit the change of bottom topography $e.g.$ a terrain-following or $\sigma$-coordinate; and for the third point, one will be tempted to use a space and time dependent coordinate that follows the isopycnal surfaces, $e.g.$ an isopycnic coordinate. 708 709 In order to satisfy two or more constrains one can even be tempted to mixed these coordinate systems, as in HYCOM (mixture of $z$-coordinate at the surface, isopycnic coordinate in the ocean interior and $\sigma$ at the ocean bottom) \citep{Chassignet_al_JPO03} or OPA (mixture of $z$-coordinate in vicinity the surface and steep topography areas and $\sigma$-coordinate elsewhere) \citep{Madec_al_JPO96} among others. 710 711 In fact one is totally free to choose any space and time vertical coordinate by introducing an arbitrary vertical coordinate : 712 \begin{equation} \label{Eq_s} 713 s=s(i,j,k,t) 714 \end{equation} 715 with the restriction that the above equation gives a single-valued monotonic relationship between $s$ and $k$, when $i$, $j$ and $t$ are held fixed. \eqref{Eq_s} is a transformation from the $(i,j,k,t)$ coordinate system with independent variables into the $(i,j,s,t)$ generalised coordinate system with $s$ depending on the other three variables through \eqref{Eq_s}. 716 This so-called \textit{generalised vertical coordinate} \citep{Kasahara_MWR74} is in fact an Arbitrary Lagrangian--Eulerian (ALE) coordinate. Indeed, choosing an expression for $s$ is an arbitrary choice that determines which part of the vertical velocity (defined from a fixed referential) will cross the levels (Eulerian part) and which part will be used to move them (Lagrangian part). 717 The coordinate is also sometime referenced as an adaptive coordinate \citep{Hofmeister_al_OM09}, since the coordinate system is adapted in the course of the simulation. Its most often used implementation is via an ALE algorithm, in which a pure lagrangian step is followed by regridding and remapping steps, the later step implicitly embedding the vertical advection \citep{Hirt_al_JCP74, Chassignet_al_JPO03, White_al_JCP09}. Here we follow the \citep{Kasahara_MWR74} strategy : a regridding step (an update of the vertical coordinate) followed by an eulerian step with an explicit computation of vertical advection relative to the moving s-surfaces. 718 719 A key point here is that the $s$-coordinate depends on $(i,j)$ ==> horizontal pressure gradient... 720 721 the generalized vertical coordinates used in ocean modelling are not orthogonal, which contrasts with many other applications in mathematical physics. Hence, it is useful to keep in mind the following properties that may seem odd on initial encounter. 722 723 the horizontal velocity in ocean models measures motions in the horizontal plane, perpendicular to the local gravitational field. That is, horizontal velocity is mathematically the same regardless the vertical coordinate, be it geopotential, isopycnal, pressure, or terrain following. The key motivation for maintaining the same horizontal velocity component is that the hydrostatic and geostrophic balances are dominant in the large-scale ocean. Use of an alternative quasi-horizontal velocity, for example one oriented parallel to the generalized surface, would lead to unacceptable numerical errors. Correspondingly, the vertical direction is anti-parallel to the gravitational force in all of the coordinate systems. We do not choose the alternative of a quasi-vertical direction oriented normal to the surface of a constant generalized vertical coordinate. 724 725 It is the method used to measure transport across the generalized vertical coordinate surfaces which differs between the vertical coordinate choices. That is, computation of the dia-surface velocity component represents the fundamental distinction between the various coordinates. In some models, such as geopotential, pressure, 726 and terrain following, this transport is typically diagnosed from volume or mass conservation. In other models, such as isopycnal layered models, this transport is prescribed based on assumptions about the physical processes producing a flux across the layer interfaces. 727 728 729 In this section we first establish the PE in the generalised vertical $s$-coordinate, then we discuss the particular cases available in \NEMO, namely $z$, $z$*, $s$, and $\tilde z$. 730 %} 943 731 944 732 % ------------------------------------------------------------------------------------------------------------- … … 1023 811 Add a few works on z and zps and s and underlies the differences between all of them 1024 812 \colorbox{yellow}{ $< = =$ end update} } 813 814 815 816 % ------------------------------------------------------------------------------------------------------------- 817 % Curvilinear z*-coordinate System 818 % ------------------------------------------------------------------------------------------------------------- 819 \subsection{Curvilinear \textit{z*}--coordinate System} 820 \label{PE_zco_star} 821 822 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 823 \begin{figure}[!b] \label{Fig_z_zstar} \begin{center} 824 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_z_zstar.pdf} 825 \caption{(a) $z$-coordinate in linear free-surface case ; (b) $z-$coordinate in non-linear 826 free surface case (c) re-scaled height coordinate (become popular as the \textit{z*-}coordinate 827 \citep{Adcroft_Campin_OM04} ).} 828 \end{center} \end{figure} 829 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 830 831 832 In that case, the free surface equation is nonlinear, and the variations of volume are fully 833 taken into account. These coordinates systems is presented in a report \citep{Levier2007} 834 available on the \NEMO web site. 835 836 %\gmcomment{ 837 The \textit{z*} coordinate approach is an unapproximated, non-linear free surface implementation 838 which allows one to deal with large amplitude free-surface 839 variations relative to the vertical resolution \citep{Adcroft_Campin_OM04}. In 840 the \textit{z*} formulation, the variation of the column thickness due to sea-surface 841 undulations is not concentrated in the surface level, as in the $z$-coordinate formulation, 842 but is equally distributed over the full water column. Thus vertical 843 levels naturally follow sea-surface variations, with a linear attenuation with 844 depth, as illustrated by figure fig.1c . Note that with a flat bottom, such as in 845 fig.1c, the bottom-following $z$ coordinate and \textit{z*} are equivalent. 846 The definition and modified oceanic equations for the rescaled vertical coordinate 847 \textit{z*}, including the treatment of fresh-water flux at the surface, are 848 detailed in Adcroft and Campin (2004). The major points are summarized 849 here. The position ( \textit{z*}) and vertical discretization (\textit{z*}) are expressed as: 850 \begin{equation} \label{Eq_z-star} 851 H + \textit{z*} = (H + z) / r \quad \text{and} \ \delta \textit{z*} = \delta z / r \quad \text{with} \ r = \frac{H+\eta} {H} 852 \end{equation} 853 Since the vertical displacement of the free surface is incorporated in the vertical 854 coordinate \textit{z*}, the upper and lower boundaries are at fixed \textit{z*} position, 855 $\textit{z*} = 0$ and $\textit{z*} = -H$ respectively. Also the divergence of the flow field 856 is no longer zero as shown by the continuity equation: 857 \begin{equation*} 858 \frac{\partial r}{\partial t} = \nabla_{\textit{z*}} \cdot \left( r \; \rm{\bf U}_h \right) 859 \left( r \; w\textit{*} \right) = 0 860 \end{equation*} 861 %} 862 863 864 % from MOM4p1 documentation 865 866 To overcome problems with vanishing surface and/or bottom cells, we consider the 867 zstar coordinate 868 \begin{equation} \label{PE_} 869 z^\star = H \left( \frac{z-\eta}{H+\eta} \right) 870 \end{equation} 871 872 This coordinate is closely related to the "eta" coordinate used in many atmospheric 873 models (see Black (1994) for a review of eta coordinate atmospheric models). It 874 was originally used in ocean models by Stacey et al. (1995) for studies of tides 875 next to shelves, and it has been recently promoted by Adcroft and Campin (2004) 876 for global climate modelling. 877 878 The surfaces of constant $z^\star$ are quasi-horizontal. Indeed, the $z^\star$ coordinate reduces to $z$ when $\eta$ is zero. In general, when noting the large differences between 879 undulations of the bottom topography versus undulations in the surface height, it 880 is clear that surfaces constant $z^\star$ are very similar to the depth surfaces. These properties greatly reduce difficulties of computing the horizontal pressure gradient relative to terrain following sigma models discussed in \S\ref{PE_sco}. 881 Additionally, since $z^\star$ when $\eta = 0$, no flow is spontaneously generated in an 882 unforced ocean starting from rest, regardless the bottom topography. This behaviour is in contrast to the case with "s"-models, where pressure gradient errors in 883 the presence of nontrivial topographic variations can generate nontrivial spontaneous flow from a resting state, depending on the sophistication of the pressure 884 gradient solver. The quasi-horizontal nature of the coordinate surfaces also facilitates the implementation of neutral physics parameterizations in $z^\star$ models using 885 the same techniques as in $z$-models (see Chapters 13-16 of \cite{Griffies_Bk04}) for a 886 discussion of neutral physics in $z$-models, as well as Section \S\ref{LDF_slp} 887 in this document for treatment in \NEMO). 888 889 The range over which $z^\star$ varies is time independent $-H \leq z^\star \leq 0$. Hence, all 890 cells remain nonvanishing, so long as the surface height maintains $\eta > ?H$. This 891 is a minor constraint relative to that encountered on the surface height when using 892 $s = z$ or $s = z - \eta$. 893 894 Because $z^\star$ has a time independent range, all grid cells have static increments 895 ds, and the sum of the ver tical increments yields the time independent ocean 896 depth %·k ds = H. 897 The $z^\star$ coordinate is therefore invisible to undulations of the 898 free surface, since it moves along with the free surface. This proper ty means that 899 no spurious ver tical transpor t is induced across surfaces of constant $z^\star$ by the 900 motion of external gravity waves. Such spurious transpor t can be a problem in 901 z-models, especially those with tidal forcing. Quite generally, the time independent 902 range for the $z^\star$ coordinate is a very convenient proper ty that allows for a nearly 903 arbitrary ver tical resolution even in the presence of large amplitude fluctuations of 904 the surface height, again so long as $\eta > -H$. 905 906 %end MOM doc %%% 907 908 909 910 \newpage 911 % ------------------------------------------------------------------------------------------------------------- 912 % Terrain following coordinate System 913 % ------------------------------------------------------------------------------------------------------------- 914 \subsection{Curvilinear Terrain-following \textit{s}--coordinate} 915 \label{PE_sco} 916 917 % ------------------------------------------------------------------------------------------------------------- 918 % Introduction 919 % ------------------------------------------------------------------------------------------------------------- 920 \subsubsection{Introduction} 921 922 Several important aspects of the ocean circulation are influenced by bottom topography. 923 Of course, the most important is that bottom topography determines deep ocean sub-basins, 924 barriers, sills and channels that strongly constrain the path of water masses, but more subtle 925 effects exist. For example, the topographic $\beta$-effect is usually larger than the planetary 926 one along continental slopes. Topographic Rossby waves can be excited and can interact 927 with the mean current. In the $z-$coordinate system presented in the previous section 928 (\S\ref{PE_zco}), $z-$surfaces are geopotential surfaces. The bottom topography is 929 discretised by steps. This often leads to a misrepresentation of a gradually sloping bottom 930 and to large localized depth gradients associated with large localized vertical velocities. 931 The response to such a velocity field often leads to numerical dispersion effects. 932 One solution to strongly reduce this error is to use a partial step representation of bottom 933 topography instead of a full step one \cite{Pacanowski_Gnanadesikan_MWR98}. 934 Another solution is to introduce a terrain-following coordinate system (hereafter $s-$coordinate) 935 936 The $s$-coordinate avoids the discretisation error in the depth field since the layers of 937 computation are gradually adjusted with depth to the ocean bottom. Relatively small 938 topographic features as well as gentle, large-scale slopes of the sea floor in the deep 939 ocean, which would be ignored in typical $z$-model applications with the largest grid 940 spacing at greatest depths, can easily be represented (with relatively low vertical resolution). 941 A terrain-following model (hereafter $s-$model) also facilitates the modelling of the 942 boundary layer flows over a large depth range, which in the framework of the $z$-model 943 would require high vertical resolution over the whole depth range. Moreover, with a 944 $s$-coordinate it is possible, at least in principle, to have the bottom and the sea surface 945 as the only boundaries of the domain (nomore lateral boundary condition to specify). 946 Nevertheless, a $s$-coordinate also has its drawbacks. Perfectly adapted to a 947 homogeneous ocean, it has strong limitations as soon as stratification is introduced. 948 The main two problems come from the truncation error in the horizontal pressure 949 gradient and a possibly increased diapycnal diffusion. The horizontal pressure force 950 in $s$-coordinate consists of two terms (see Appendix~\ref{Apdx_A}), 951 952 \begin{equation} \label{Eq_PE_p_sco} 953 \left. {\nabla p} \right|_z =\left. {\nabla p} \right|_s -\frac{\partial 954 p}{\partial s}\left. {\nabla z} \right|_s 955 \end{equation} 956 957 The second term in \eqref{Eq_PE_p_sco} depends on the tilt of the coordinate surface 958 and introduces a truncation error that is not present in a $z$-model. In the special case 959 of a $\sigma-$coordinate (i.e. a depth-normalised coordinate system $\sigma = z/H$), 960 \citet{Haney1991} and \citet{Beckmann1993} have given estimates of the magnitude 961 of this truncation error. It depends on topographic slope, stratification, horizontal and 962 vertical resolution, the equation of state, and the finite difference scheme. This error 963 limits the possible topographic slopes that a model can handle at a given horizontal 964 and vertical resolution. This is a severe restriction for large-scale applications using 965 realistic bottom topography. The large-scale slopes require high horizontal resolution, 966 and the computational cost becomes prohibitive. This problem can be at least partially 967 overcome by mixing $s$-coordinate and step-like representation of bottom topography \citep{Gerdes1993a,Gerdes1993b,Madec_al_JPO96}. However, the definition of the model 968 domain vertical coordinate becomes then a non-trivial thing for a realistic bottom 969 topography: a envelope topography is defined in $s$-coordinate on which a full or 970 partial step bottom topography is then applied in order to adjust the model depth to 971 the observed one (see \S\ref{DOM_zgr}. 972 973 For numerical reasons a minimum of diffusion is required along the coordinate surfaces 974 of any finite difference model. It causes spurious diapycnal mixing when coordinate 975 surfaces do not coincide with isoneutral surfaces. This is the case for a $z$-model as 976 well as for a $s$-model. However, density varies more strongly on $s-$surfaces than 977 on horizontal surfaces in regions of large topographic slopes, implying larger diapycnal 978 diffusion in a $s$-model than in a $z$-model. Whereas such a diapycnal diffusion in a 979 $z$-model tends to weaken horizontal density (pressure) gradients and thus the horizontal 980 circulation, it usually reinforces these gradients in a $s$-model, creating spurious circulation. 981 For example, imagine an isolated bump of topography in an ocean at rest with a horizontally 982 uniform stratification. Spurious diffusion along $s$-surfaces will induce a bump of isoneutral 983 surfaces over the topography, and thus will generate there a baroclinic eddy. In contrast, 984 the ocean will stay at rest in a $z$-model. As for the truncation error, the problem can be reduced by introducing the terrain-following coordinate below the strongly stratified portion of the water column 985 ($i.e.$ the main thermocline) \citep{Madec_al_JPO96}. An alternate solution consists of rotating 986 the lateral diffusive tensor to geopotential or to isoneutral surfaces (see \S\ref{PE_ldf}. 987 Unfortunately, the slope of isoneutral surfaces relative to the $s$-surfaces can very large, 988 strongly exceeding the stability limit of such a operator when it is discretized (see Chapter~\ref{LDF}). 989 990 The $s-$coordinates introduced here \citep{Lott_al_OM90,Madec_al_JPO96} differ mainly in two 991 aspects from similar models: it allows a representation of bottom topography with mixed 992 full or partial step-like/terrain following topography ; It also offers a completely general 993 transformation, $s=s(i,j,z)$ for the vertical coordinate. 994 995 996 \newpage 997 % ------------------------------------------------------------------------------------------------------------- 998 % Curvilinear z-tilde coordinate System 999 % ------------------------------------------------------------------------------------------------------------- 1000 \subsection{Curvilinear $\tilde{z}$--coordinate} 1001 \label{PE_zco_tilde} 1002 1003 1025 1004 1026 1005 \newpage
Note: See TracChangeset
for help on using the changeset viewer.
