- Timestamp:
- 2010-04-12T16:59:59+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_MISC.tex
r1225 r1831 2 2 % Chapter Ñ Miscellaneous Topics 3 3 % ================================================================ 4 \chapter{Miscellaneous Topics (xxx)}4 \chapter{Miscellaneous Topics} 5 5 \label{MISC} 6 6 \minitoc 7 7 8 \newpage 9 $\ $\newline % force a new ligne 10 8 11 % ================================================================ 9 12 % Representation of Unresolved Straits … … 12 15 \label{MISC_strait} 13 16 14 % In climate modeling, it is often necessary to allow water masses that are separated by land to exchange properties. This situation arises in models when the grid mesh is too coarse to resolve narrow passageways that in reality provide crucial connections between water masses. For example, coarse grid spacing typically closes off the Mediterranean from the Atlantic at the Straits of Gibraltar. In this case, it is important for climate models to include the effects of salty water entering the Atlantic from the Mediterranean. Likewise, it is important for the Mediterranean to replenish its supply of water from the Atlantic to balance the net evaporation occurring over the Mediterranean region. 15 % Same problem occurs as soon as important straits are not properly resolved by the mesh. For example it arises in ORCA 1/4\r{} in several straits of the Indonesian archipelago (Ombai, Lombok...). 16 % We describe here the three methods that can be used in \NEMO to handle such improperly resolved straits. 17 %The problem is resolved here by allowing the mixing of tracers and mass/volume between non-adjacent water columns at nominated regions within the model. Momentum is not mixed. The scheme conserves total tracer content, and total volume (the latter in $z*$- or $s*$-coordinate), and maintains compatibility between the tracer and mass/volume budgets. 18 19 % Note: such changes are so specific to a given configuration that no attempt has been made to set them in a generic way. Nevertheless, an example of how they can be set up is given in the ORCA2 configuration (\key{ORCA_R2}). 17 In climate modeling, it often occurs that a crucial connections between water masses 18 is broken as the grid mesh is too coarse to resolve narrow straits. For example, coarse 19 grid spacing typically closes off the Mediterranean from the Atlantic at the Strait of 20 Gibraltar. In this case, it is important for climate models to include the effects of salty 21 water entering the Atlantic from the Mediterranean. Likewise, it is important for the 22 Mediterranean to replenish its supply of water from the Atlantic to balance the net 23 evaporation occurring over the Mediterranean region. This problem occurs even in 24 eddy permitting simulations. For example, in ORCA 1/4\r{} several straits of the Indonesian 25 archipelago (Ombai, Lombok...) are much narrow than even a single ocean grid-point. 26 27 We describe briefly here the three methods that can be used in \NEMO to handle 28 such improperly resolved straits. The first two consist of opening the strait by hand 29 while ensuring that the mass exchanges through the strait are not too large by 30 either artificially reducing the surface of the strait grid-cells or, locally increasing 31 the lateral friction. In the third one, the strait is closed but exchanges of mass, 32 heat and salt across the land are allowed. 33 Note that such modifications are so specific to a given configuration that no attempt 34 has been made to set them in a generic way. However, examples of how 35 they can be set up is given in the ORCA 2\r{} and 0.5\r{} configurations (search for 36 \key{ORCA\_R2} or \key{ORCA\_R05} in the code). 20 37 21 38 % ------------------------------------------------------------------------------------------------------------- … … 25 42 \label{MISC_strait_hand} 26 43 27 $\bullet$ reduced scale factor, also called partially open face (Hallberg, personnal communication 2006) 28 %also called partially open cell or partial closed cells 29 30 $\bullet$ increase of the viscous boundary layer thickness by local increase of the fmask value at the coast 31 32 \colorbox{yellow}{Add a short description of scale factor changes staff and fmask increase} 33 44 $\bullet$ reduced scale factor in the cross-strait direction to a value in better agreement 45 with the true mean width of the strait. (Fig.~\ref{Fig_MISC_strait_hand}). 46 This technique is sometime called "partially open face" or "partially closed cells". 47 The key issue here is only to reduce the faces of $T$-cell ($i.e.$ change the value 48 of the horizontal scale factors at $u$- or $v$-point) but not the volume of the $T$-cell. 49 Indeed, reducing the volume of strait $T$-cell can easily produce a numerical 50 instability at that grid point that would require a reduction of the model time step. 51 The changes associated with strait management are done in \mdl{domhgr}, 52 just after the definition or reading of the horizontal scale factors. 53 54 $\bullet$ increase of the viscous boundary layer thickness by local increase of the 55 fmask value at the coast (Fig.~\ref{Fig_MISC_strait_hand}). This is done in 56 \mdl{dommsk} together with the setting of the coastal value of fmask 57 (see Section \ref{LBC_coast}) 34 58 35 59 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 51 75 % Cross Land Advection 52 76 % ------------------------------------------------------------------------------------------------------------- 53 \subsection{Cross Land Advection (\ textit{tracla} module)}77 \subsection{Cross Land Advection (\mdl{tracla})} 54 78 \label{MISC_strait_cla} 55 56 79 %--------------------------------------------namcla-------------------------------------------------------- 57 80 \namdisplay{namcla} … … 60 83 \colorbox{yellow}{Add a short description of CLA staff here or in lateral boundary condition chapter?} 61 84 85 %The problem is resolved here by allowing the mixing of tracers and mass/volume between non-adjacent water columns at nominated regions within the model. Momentum is not mixed. The scheme conserves total tracer content, and total volume (the latter in $z*$- or $s*$-coordinate), and maintains compatibility between the tracer and mass/volume budgets. 86 62 87 % ================================================================ 63 88 % Closed seas 64 89 % ================================================================ 65 \section{Closed seas }90 \section{Closed seas (\mdl{closea})} 66 91 \label{MISC_closea} 92 93 \colorbox{yellow}{Add here a short description of the way closed seas are managed} 67 94 68 95 … … 81 108 model configuration while having a small computer memory requirement. 82 109 It can also be used to easily test specific physics in a sub-domain (for example, 83 see \citep{Madec 1996} for a test of the coupling used in the global ocean110 see \citep{Madec_al_JPO96} for a test of the coupling used in the global ocean 84 111 version of OPA between sea-ice and ocean model over the Arctic or Antarctic 85 112 ocean, using a sub-domain). In the standard model, this option does not … … 120 147 \label{MISC_1D} 121 148 122 The 1D model option simulates a stand alone water column within the 3D \NEMO 123 system. It can be applied to the ocean alone or to the ocean-ice system and can 124 include passive tracers or a biogeochemical model. It is set up by defining the 125 \key{cfg\_1d} CPP key. This 1D model is a very useful tool \textit{(a)} to learn 126 about the physics and numerical treatment of vertical mixing processes ; \textit{(b)} 127 to investigate suitable parameterisations of unresolved turbulence (wind steering, 128 langmuir circulation, skin layers) ; \textit{(c)} to compare the behaviour of different 129 vertical mixing schemes ; \textit{(d)} to perform sensitivity studies on the vertical 130 diffusion at a particular point of the ocean global domain ; \textit{(d)} to produce 131 extra diagnostics, without the large memory requirement of the full 3D model. 132 133 The methodology is based on the use of the zoom functionality (see 134 \S\ref{MISC_zoom}), with some extra routines. There is no need to define a new 135 mesh, bathymetry, initial state or forcing, since the 1D model will use those of the 136 configuration it is a zoom of. 149 The 1D model option simulates a stand alone water column within the 3D \NEMO system. 150 It can be applied to the ocean alone or to the ocean-ice system and can include passive tracers 151 or a biogeochemical model. It is set up by defining the \key{cfg\_1d} CPP key. 152 The 1D model is a very useful tool 153 \textit{(a)} to learn about the physics and numerical treatment of vertical mixing processes ; 154 \textit{(b)} to investigate suitable parameterisations of unresolved turbulence (wind steering, 155 langmuir circulation, skin layers) ; 156 \textit{(c)} to compare the behaviour of different vertical mixing schemes ; 157 \textit{(d)} to perform sensitivity studies on the vertical diffusion at a particular point of an ocean domain ; 158 \textit{(d)} to produce extra diagnostics, without the large memory requirement of the full 3D model. 159 160 The methodology is based on the use of the zoom functionality over the smallest possible 161 domain : a 3 x 3 domain centred on the grid point of interest (see \S\ref{MISC_zoom}), 162 with some extra routines. There is no need to define a new mesh, bathymetry, 163 initial state or forcing, since the 1D model will use those of the configuration it is a zoom of. 164 The chosen grid point is set in par\_oce.F90 module by setting the jpizoom and jpjzoom 165 parameters to the indices of the location of the chosen grid point. 137 166 138 167 % ================================================================ … … 145 174 %-------------------------------------------------------------------------------------------------------------- 146 175 147 % Steven update just bellow... not sure it is what I vant to say 148 %Searching an equilibrium state with an ocean model requires repeated very long time 149 %integrations (each of a few thousand years for a global model). Due to the size of 150 Searching an equilibrium state with an ocean model requires very long time 151 integration (a few thousand years for a global model). Due to the size of 152 the time step required for numerical stability (less than a 153 few hours), this usually requires a large elapsed time. In order to overcome 154 this problem, \citet{Bryan1984} introduces a technique that is intended to accelerate 176 Searching an equilibrium state with an global ocean model requires a very long time 177 integration period (a few thousand years for a global model). Due to the size of 178 the time step required for numerical stability (less than a few hours), 179 this usually requires a large elapsed time. In order to overcome this problem, 180 \citet{Bryan1984} introduces a technique that is intended to accelerate 155 181 the spin up to equilibrium. It uses a larger time step in 156 the t hermodynamic evolution equations than in the dynamicevolution182 the tracer evolution equations than in the momentum evolution 157 183 equations. It does not affect the equilibrium solution but modifies the 158 184 trajectory to reach it. 159 185 160 186 The acceleration of convergence option is used when \np{nn\_acc}=1. In that case, 161 $\Delta t=rdt$ is the time step of dynamics while $\widetilde{\Delta t}=rdttra$ is the 162 tracer time-step. Their values are set from the \np{rdt} and \np{rdttra} namelist 163 parameters. The set of prognostic equations to solve becomes: 187 $\Delta t=rn\_rdt$ is the time step of dynamics while $\widetilde{\Delta t}=rdttra$ is the 188 tracer time-step. the former is set from the \np{rn\_rdt} namelist parameter while the latter 189 is computed using a hyperbolic tangent profile and the following namelist parameters : 190 \np{rn\_rdtmin}, \np{rn\_rdtmax} and \np{rn\_rdth}. Those three parameters correspond 191 to the surface value the deep ocean value and the depth at which the transition occurs, respectively. 192 The set of prognostic equations to solve becomes: 164 193 \begin{equation} \label{Eq_acc} 165 194 \begin{split} … … 183 212 $\widetilde{ \Delta t}$ implies that the heat and salt contents are no longer 184 213 conserved due to the vertical coupling of the ocean level through both 185 advection and diffusion. 186 % first mention of depth varying tilda-delta-t! and namelist parameter associated to that are to be described 187 214 advection and diffusion. Therefore \np{rn\_rdtmin} = \np{rn\_rdtmax} should be 215 a more clever choice. 188 216 189 217 % ================================================================ … … 202 230 $\bullet$ Vector and memory optimisation: 203 231 204 \key{vectopt\_loop} enables internal loop collapse, a very efficient way to increase 205 the length of vector calculations and thus speed up the model on vector computers. 232 \key{vectopt\_loop} enables the internal loops to collapse. This is very 233 a very efficient way to increase the length of vector calculations and thus 234 to speed up the model on vector computers. 206 235 207 236 % Add here also one word on NPROMA technique that has been found useless, since compiler have made significant progress during the last decade. … … 209 238 % Add also one word on NEC specific optimisation (Novercheck option for example) 210 239 211 \key{vectopt\_memory} has been introduced in order to reduce the memory 212 requirement of the model. This is obviously done at the cost of increasing the 213 CPU time requirement, since it suppress intermediate computations which would 214 have been saved in in-core memory. This feature has not been intensively used. 215 In fact, currently it is only implemented for the TKE physics, in which, when 216 \key{vectopt\_memory} is defined, the coefficients used for horizontal smoothing 217 of $A_v^T$ and $A_v^m$ are no longer computed once and for all. This reduces 218 the memory requirement by three 2D arrays. 240 \key{vectopt\_memory} is an obsolescent option. It has been introduced in order 241 to reduce the memory requirement of the model at a time when in-core memory 242 were rather limited. This is obviously done at the cost of increasing the CPU 243 time requirement, since it suppress intermediate computations which would have 244 been saved in in-core memory. Currently it is only used in the old implementation 245 of the TKE physics (\key{tke\_old}) where, when \key{vectopt\_memory} 246 is defined, the coefficients used for horizontal smoothing of $A_v^T$ and $A_v^m$ 247 are no longer computed once and for all. This reduces the memory requirement by three 248 3D arrays. This option will disappear in the next \NEMO release. 219 249 220 250 … … 231 261 management. When defined, this key forces the activation of all options and 232 262 CPP keys. For example, all tracer and momentum advection schemes are called! 233 There is therefore no physical meaning associated with the model results.263 Therefore the model results have no physical meaning. 234 264 However, this option forces both the compiler and the model to run through 235 265 all the \textsc{Fortran} lines of the model. This allows the user to check for obvious 236 266 compilation or execution errors with all CPP options, and errors in namelist options. 237 267 238 3- \key{esopa} (to be rename key\_nemo) : which is another option for model 239 management. When defined, this key forces the activation of all options and CPP keys. 240 For example, all tracer and momentum advection schemes are called! There is 241 therefore no physical meaning associated with the model results. However, this option 242 forces both the compiler and the model to run through all the Fortran lines of the model. 243 This allows the user to check for obvious compilation or execution errors with all CPP 244 options, and errors in namelist options. 245 246 4- last digit comparison (\np{nbit\_cmp}). In an MPP simulation, the computation of 268 4- last digit comparison (\np{nn\_bit\_cmp}). In an MPP simulation, the computation of 247 269 a sum over the whole domain is performed as the summation over all processors of 248 270 each of their sums over their interior domains. This double sum never gives exactly … … 251 273 mono-processor and multi-processor runs give exactly the same results. 252 274 253 $\bullet$ Benchmark (\np{n bench}). This option defines a benchmark run based on275 $\bullet$ Benchmark (\np{nn\_bench}). This option defines a benchmark run based on 254 276 a GYRE configuration in which the resolution remains the same whatever the domain 255 277 size. This allows a very large model domain to be used, just by changing the domain … … 261 283 % Elliptic solvers (SOL) 262 284 % ================================================================ 263 \section{Elliptic solvers (SOL directory)}285 \section{Elliptic solvers (SOL)} 264 286 \label{MISC_sol} 265 287 %--------------------------------------------namdom------------------------------------------------------- … … 267 289 %-------------------------------------------------------------------------------------------------------------- 268 290 269 The computation of the surface pressure gradient with a rigid lid assumption270 requires the computation of $\partial_t \psi$, the time evolution of the 271 barotropic streamfunction. $\partial_t \psi$ is the solution of an elliptic272 equation \eqref{Eq_PE_psi} for which three solvers are available: a273 Successive-Over-Relaxation scheme (SOR), a preconditioned conjugate gradient 274 scheme(PCG) \citep{Madec1988, Madec1990} and a Finite Elements Tearing and275 Interconnecting scheme (FETI) \citep{Guyon_al_EP99, Guyon_al_CalPar99}. 291 When the filtered sea surface height option is used, the surface pressure gradient is 292 computed in \mdl{dynspg\_flt}. The force added in the momentum equation is solved implicitely. 293 It is thus solution of an elliptic equation \eqref{Eq_PE_flt} for which two solvers are available: 294 a Successive-Over-Relaxation scheme (SOR) and a preconditioned conjugate gradient 295 scheme(PCG) \citep{Madec_al_OM88, Madec_PhD90}. The solver is selected trough the 296 the value of \np{nn\_solv} (namelist parameter). 297 276 298 The PCG is a very efficient method for solving elliptic equations on vector computers. 277 299 It is a fast and rather easy method to use; which are attractive features for a large … … 280 302 a search for an optimal parameter as in the SOR method. However, the SOR has 281 303 been retained because it is a linear solver, which is a very useful property when 282 using the adjoint model of \NEMO. The FETI solver is a very efficient method on 283 massively parallel computers. However, it has never been used since OPA 8.0. 284 The current version in \NEMO should not even successfully go through the 285 compilation phase. 286 287 The surface pressure gradient is computed in \mdl{dynspg}. The solver used 288 (PCG or SOR) depending on the value of \np{nsolv} (namelist parameter). 289 At each time step the time derivative of the barotropic streamfunction is 290 the solution of \eqref{Eq_psi_total}. Introducing the following coefficients: 291 292 \begin{equation} 304 using the adjoint model of \NEMO. 305 306 At each time step, the time derivative of the sea surface height at time step $t+1$ 307 (or equivalently the divergence of the \textit{after} barotropic transport) that appears 308 in the filtering forced is the solution of the elliptic equation obtained from the horizontal 309 divergence of the vertical summation of \eqref{Eq_PE_flt}. 310 Introducing the following coefficients: 311 \begin{equation} \label{Eq_sol_matrix} 293 312 \begin{aligned} 294 & C_{i,j}^{NS} &&= \frac{e_{2v}(i,j)}{\left( H_v (i,j) e_{1v} (i,j)\right)}\\295 & C_{i,j}^{EW} &&= \frac{e_{1u}(i,j)}{\left( H_u (i,j) e_{2u} (i,j)\right)}\\296 & B_{i,j} &&= \delta_i \left( e_{2v}M_v \right) - \delta_j \left( e_{1u}M_u \right)\\313 &c_{i,j}^{NS} &&= {2 \Delta t }^2 \; \frac{H_v (i,j) \; e_{1v} (i,j)}{e_{2v}(i,j)} \\ 314 &c_{i,j}^{EW} &&= {2 \Delta t }^2 \; \frac{H_u (i,j) \; e_{2u} (i,j)}{e_{1u}(i,j)} \\ 315 &b_{i,j} &&= \delta_i \left[ e_{2u}M_u \right] - \delta_j \left[ e_{1v}M_v \right]\ , \\ 297 316 \end{aligned} 298 317 \end{equation} 299 300 318 the five-point finite difference equation \eqref{Eq_psi_total} can be rewritten as: 301 302 \begin{multline} \label{Eq_solmat} 303 C_{i+1,j}^{NS} { \left( \frac{\partial \psi}{\partial t } \right) }_{i+1,j} + \; 304 C_{i,j+1}^{EW}{ \left( \frac{\partial \psi}{\partial t } \right) }_{i,j+1} + \; 305 C_{i,j} ^{NS} { \left( \frac{\partial \psi}{\partial t } \right) }_{i-1,j} + \; 306 C_{i,j} ^{EW}{ \left( \frac{\partial \psi}{\partial t } \right) }_{i,j-1}\\ 307 -\left(C_{i+1,j}^{NS} + C_{i,j+1}^{EW} + C_{i,j}^{NS} + C_{i,j}^{EW} \right) { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j} = B_{i,j} 308 \end{multline} 319 \begin{equation} \label{Eq_solmat} 320 \begin{split} 321 c_{i+1,j}^{NS} D_{i+1,j} + \; c_{i,j+1}^{EW} D_{i,j+1} 322 + c_{i,j} ^{NS} D_{i-1,j} + \; c_{i,j} ^{EW} D_{i,j-1} & \\ 323 - \left(c_{i+1,j}^{NS} + c_{i,j+1}^{EW} + c_{i,j}^{NS} + c_{i,j}^{EW} \right) D_{i,j} &= b_{i,j} 324 \end{split} 325 \end{equation} 309 326 \eqref{Eq_solmat} is a linear symmetric system of equations. All the elements of 310 327 the corresponding matrix \textbf{A} vanish except those of five diagonals. With … … 315 332 \eqref{Eq_solmat}, is a vector. 316 333 334 Note that in the linear free surface case, the depth that appears in \eqref{Eq_sol_matrix} 335 does not vary with time, and thus the matrix can be computed once for all. In non-linear free surface 336 (\key{vvl} defined) the matrix have to be updated at each time step. 337 317 338 % ------------------------------------------------------------------------------------------------------------- 318 339 % Successive Over Relaxation 319 340 % ------------------------------------------------------------------------------------------------------------- 320 \subsection{Successive Over Relaxation \np{nsolv}=2}341 \subsection{Successive Over Relaxation (\np{nn\_solv}=2, \mdl{solsor})} 321 342 \label{MISC_solsor} 322 343 323 Let us introduce the four cardinal coefficients: $A_{i,j}^S = C_{i,j}^{NS}/D_{i,j}\,$,324 $A_{i,j}^W = C_{i,j}^{EW}/D_{i,j}\,$, $A_{i,j}^E = C_{i,j+1}^{EW}/D_{i,j}\,$ and $A_{i,j}^N = C_{i+1,j}^{NS}/D_{i,j}\,$, and define 325 $\tilde B_{i,j} = B_{i,j}/D_{i,j}\,$, where $D_{i,j} = C_{i,j}^{NS}+ C_{i+1,j}^{NS} + C_{i,j}^{EW} + C_{i,j+1}^{EW} $ (i.e. the diagonal of 326 \textbf{A}). (VII.5.1) can be rewritten as: 327 328 \begin{multline} \label{Eq_solmat_p} 329 A_{i,j}^{N} { \left( \frac{\partial \psi}{\partial t } \right) }_{i+1,j } 330 +\,A_{i,j}^{E} { \left( \frac{\partial \psi}{\partial t } \right) }_{i ,j+1} \\ 331 +\,A_{i,j}^{S} { \left( \frac{\partial \psi}{\partial t } \right) }_{i-1 ,j } 332 +\,A_{i,j}^{W} { \left( \frac{\partial \psi}{\partial t } \right) }_{i ,j-1} 333 - { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j} = \tilde B_{i,j}334 \end{ multline}335 336 The SOR method used is an iterative method. Its algorithm can be summarised337 as follows [see \citet{Haltiner1980} for a further discussion]:344 Let us introduce the four cardinal coefficients: 345 \begin{align*} 346 a_{i,j}^S &= c_{i,j }^{NS}/d_{i,j} &\qquad a_{i,j}^W &= c_{i,j}^{EW}/d_{i,j} \\ 347 a_{i,j}^E &= c_{i,j+1}^{EW}/d_{i,j} &\qquad a_{i,j}^N &= c_{i+1,j}^{NS}/d_{i,j} 348 \end{align*} 349 where $d_{i,j} = c_{i,j}^{NS}+ c_{i+1,j}^{NS} + c_{i,j}^{EW} + c_{i,j+1}^{EW}$ 350 (i.e. the diagonal of the matrix). \eqref{Eq_solmat} can be rewritten as: 351 \begin{equation} \label{Eq_solmat_p} 352 \begin{split} 353 a_{i,j}^{N} D_{i+1,j} +\,a_{i,j}^{E} D_{i,j+1} +\, a_{i,j}^{S} D_{i-1,j} +\,a_{i,j}^{W} D_{i,j-1} - D_{i,j} = \tilde{b}_{i,j} 354 \end{split} 355 \end{equation} 356 with $\tilde b_{i,j} = b_{i,j}/d_{i,j}$. \eqref{Eq_solmat_p} is the equation actually solved 357 with the SOR method. This method used is an iterative one. Its algorithm can be 358 summarised as follows (see \citet{Haltiner1980} for a further discussion): 338 359 339 360 initialisation (evaluate a first guess from previous time step computations) 340 361 \begin{equation} 341 \left( \frac{\partial \psi}{\partial t } \right)_{i,j}^0 = 2{\left( \frac{\partial \psi}{\partial t } \right)_{i,j}^t} - {\left( \frac{\partial \psi}{\partial t } \right)_{i,j}^{t-1}}362 D_{i,j}^0 = 2 \, D_{i,j}^t - D_{i,j}^{t-1} 342 363 \end{equation} 343 364 iteration $n$, from $n=0$ until convergence, do : 344 \begin{multline} 345 R_{i,j}^n 346 = A_{i,j}^{N} { \left( \frac{\partial \psi}{\partial t } \right) }_{i+1,j}^n 347 +\,A_{i,j}^{E} { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j+1} ^n \\ 348 +\,A_{i,j}^{S} { \left( \frac{\partial \psi}{\partial t } \right) }_{i-1,j} ^{n+1} 349 +\,A_{i,j}^{W} { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j-1} ^{n+1} 350 - { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j}^n - \tilde B_{i,j} 351 \end{multline} 352 \begin{equation} 353 { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j} ^{n+1} 354 = { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j} ^{n} 355 + \omega \;R_{i,j}^n 365 \begin{equation} \label{Eq_sor_algo} 366 \begin{split} 367 R_{i,j}^n = &a_{i,j}^{N} D_{i+1,j}^n +\,a_{i,j}^{E} D_{i,j+1} ^n 368 +\, a_{i,j}^{S} D_{i-1,j} ^{n+1}+\,a_{i,j}^{W} D_{i,j-1} ^{n+1} 369 - D_{i,j}^n - \tilde{b}_{i,j} \\ 370 D_{i,j} ^{n+1} = &D_{i,j} ^{n} + \omega \;R_{i,j}^n 371 \end{split} 356 372 \end{equation} 357 373 where \textit{$\omega $ }satisfies $1\leq \omega \leq 2$. An optimal value exists for … … 359 375 adjusted empirically for each model domain (except for a uniform grid where an 360 376 analytical expression for \textit{$\omega$} can be found \citep{Richtmyer1967}). 361 The value of $\omega$ is set using \ textit{sor}, a \textbf{namelist} parameter.377 The value of $\omega$ is set using \np{rn\_sor}, a \textbf{namelist} parameter. 362 378 The convergence test is of the form: 363 379 \begin{equation} 364 \delta = \frac{\sum\limits_{i,j}{R_{i,j}^n}{R_{i,j}^n}}{\sum\limits_{i,j}{ \tilde B_{i,j}^n}{\tilde B_{i,j}^n}} \leq \epsilon 380 \delta = \frac{\sum\limits_{i,j}{R_{i,j}^n}{R_{i,j}^n}} 381 {\sum\limits_{i,j}{ \tilde{b}_{i,j}^n}{\tilde{b}_{i,j}^n}} \leq \epsilon 365 382 \end{equation} 366 383 where $\epsilon$ is the absolute precision that is required. It is recommended 367 that a value smaller or equal to $10^{-3}$ is used for $\epsilon$ since larger 368 values may lead to numerically induced basin scale barotropic oscillations. In fact, 369 for an eddy resolving configuration or in a filtered free surface case, a value three 370 orders of magnitude smaller than this should be used. The precision is specified by 371 setting \textit{eps} (\textbf{namelist parameter}). In addition, two other tests are 372 used to halt the iterative algorithm. They involve the number of iterations and the 373 modulus of the right hand side. If the former exceeds a specified value, \textit{nmax} 374 (\textbf{namelist parameter}), or the latter is greater than $10^{15}$, the whole 375 model computation is stopped and the last computed time step fields are saved 376 in the standard output file. In both cases, this usually indicates that there is something 377 wrong in the model configuration (an error in the mesh, the initial state, the input forcing, 384 that a value smaller or equal to $10^{-6}$ is used for $\epsilon$ since larger 385 values may lead to numerically induced basin scale barotropic oscillations. 386 The precision is specified by setting \np{rn\_eps} (\textbf{namelist} parameter). 387 In addition, two other tests are used to halt the iterative algorithm. They involve 388 the number of iterations and the modulus of the right hand side. If the former 389 exceeds a specified value, \np{nn\_max} (\textbf{namelist} parameter), 390 or the latter is greater than $10^{15}$, the whole model computation is stopped 391 and the last computed time step fields are saved in a abort.nc NetCDF file. 392 In both cases, this usually indicates that there is something wrong in the model 393 configuration (an error in the mesh, the initial state, the input forcing, 378 394 or the magnitude of the time step or of the mixing coefficients). A typical value of 379 $n max$ is a few hundred when $\epsilon = 10^{-6}$, increasing to a few395 $nn\_max$ is a few hundred when $\epsilon = 10^{-6}$, increasing to a few 380 396 thousand when $\epsilon = 10^{-12}$. 381 397 The vectorization of the SOR algorithm is not straightforward. The scheme 382 398 contains two linear recurrences on $i$ and $j$. This inhibits the vectorisation. 383 Therefore it hasbeen rewritten as:384 385 \begin{ multline}399 \eqref{Eq_sor_algo} can be been rewritten as: 400 \begin{equation} 401 \begin{split} 386 402 R_{i,j}^n 387 = A_{i,j}^{N} { \left( \frac{\partial \psi}{\partial t } \right) }_{i+1,j}^n 388 +\,A_{i,j}^{E} { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j+1} ^n \\ 389 +\,A_{i,j}^{S} { \left( \frac{\partial \psi}{\partial t } \right) }_{i-1,j} ^{n} 390 +\,A_{i,j}^{W} { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j-1} ^{n} 391 - { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j}^n - \tilde B_{i,j} 392 \end{multline} 393 394 \begin{equation} 395 R_{i,j}^n = R_{i,j}^n - \omega \;A_{i,j}^{S}\; R_{i,j-1}^n 396 \end{equation} 397 398 \begin{equation} 399 R_{i,j}^n = R_{i,j}^n - \omega \;A_{i,j}^{W}\; R_{i-1,j}^n 400 \end{equation} 401 402 The SOR method is very flexible and can be used under a wide 403 range of conditions, including irregular boundaries, interior boundary 404 points, etc. Proofs of convergence, etc. may be found in the standard 405 numerical methods texts for partial differential equations. 403 = &a_{i,j}^{N} D_{i+1,j}^n +\,a_{i,j}^{E} D_{i,j+1} ^n 404 +\,a_{i,j}^{S} D_{i-1,j} ^{n}+\,_{i,j}^{W} D_{i,j-1} ^{n} - D_{i,j}^n - \tilde{b}_{i,j} \\ 405 R_{i,j}^n = &R_{i,j}^n - \omega \;a_{i,j}^{S}\; R_{i,j-1}^n \\ 406 R_{i,j}^n = &R_{i,j}^n - \omega \;a_{i,j}^{W}\; R_{i-1,j}^n 407 \end{split} 408 \end{equation} 409 This technique slightly increases the number of iteration required to reach the convergence, 410 but this is largely compensated by the gain obtained by the suppression of the recurrences. 411 412 Another technique have been chosen, the so-called red-black SOR. It consist in solving successively 413 \eqref{Eq_sor_algo} for odd and even grid points. It also slightly reduced the convergence rate 414 but allows the vectorisation. In addition, and this is the reason why it has been chosen, it is able to handle the north fold boundary condition used in ORCA configuration ($i.e.$ tri-polar global ocean mesh). 415 416 The SOR method is very flexible and can be used under a wide range of conditions, 417 including irregular boundaries, interior boundary points, etc. Proofs of convergence, etc. 418 may be found in the standard numerical methods texts for partial differential equations. 406 419 407 420 % ------------------------------------------------------------------------------------------------------------- 408 421 % Preconditioned Conjugate Gradient 409 422 % ------------------------------------------------------------------------------------------------------------- 410 \subsection{Preconditioned Conjugate Gradient }423 \subsection{Preconditioned Conjugate Gradient (\np{nn\_solv}=1, \mdl{solpcg}) } 411 424 \label{MISC_solpcg} 412 413 \begin{flushright}414 (\textit{nbsfs=1}, \textbf{namelist parameter})415 \end{flushright}416 425 417 426 \textbf{A} is a definite positive symmetric matrix, thus solving the linear … … 448 457 449 458 initialisation : 450 451 459 \begin{equation*} 452 \textbf{x}^0 = \left( \frac{\partial \psi}{\partial t } \right)_{i,j}^0 = 2{\left( \frac{\partial \psi}{\partial t } \right)_{i,j}^t} - {\left( \frac{\partial \psi}{\partial t } \right)_{i,j}^{t-1}} 453 , \text{the initial guess } 460 \begin{split} 461 \textbf{x}^0 &= D_{i,j}^0 = 2 D_{i,j}^t - D_{i,j}^{t-1} \quad, \text{the initial guess } \\ 462 \textbf{r}^0 &= \textbf{d}^0 = \textbf{b} - \textbf{A x}^0 \\ 463 \gamma_0 &= \langle{ \textbf{r}^0 , \textbf{r}^0} \rangle 464 \end{split} 454 465 \end{equation*} 455 466 456 \begin{equation*}457 \textbf{r}^0 = \textbf{d}^0 = \textbf{b} - \textbf{A x}^0458 \end{equation*}459 \begin{equation*}460 \gamma_0 = \langle{ \textbf{r}^0 , \textbf{r}^0} \rangle461 \end{equation*}462 467 iteration $n,$ from $n=0$ until convergence, do : 463 464 468 \begin{equation} 465 469 \begin{split} 466 470 \text{z}^n& = \textbf{A d}^n \\ 467 \alpha_n &= \gamma_n \langle{ \textbf{z}^n , \textbf{d}^n} \rangle \\471 \alpha_n &= \gamma_n / \langle{ \textbf{z}^n , \textbf{d}^n} \rangle \\ 468 472 \textbf{x}^{n+1} &= \textbf{x}^n + \alpha_n \,\textbf{d}^n \\ 469 473 \textbf{r}^{n+1} &= \textbf{r}^n - \alpha_n \,\textbf{z}^n \\ … … 476 480 477 481 The convergence test is: 478 479 482 \begin{equation} 480 483 \delta = \gamma_{n}\; / \langle{ \textbf{b} , \textbf{b}} \rangle \leq \epsilon 481 484 \end{equation} 482 485 where $\epsilon $ is the absolute precision that is required. As for the SOR algorithm, 483 the whole model computation is stopped when the number of iterations, $nmax$, or486 the whole model computation is stopped when the number of iterations, \np{nn\_max}, or 484 487 the modulus of the right hand side of the convergence equation exceeds a 485 488 specified value (see \S\ref{MISC_solsor} for a further discussion). The required 486 489 precision and the maximum number of iterations allowed are specified by setting 487 $eps$ and $nmax$ (\textbf{namelist parameters}).490 \np{rn\_eps} and \np{nn\_max} (\textbf{namelist} parameters). 488 491 489 492 It can be demonstrated that the above algorithm is optimal, provides the exact 490 493 solution in a number of iterations equal to the size of the matrix, and that 491 494 the convergence rate is faster as the matrix is closer to the identity matrix, 492 $i.e.$ its eigenvalues are closer to 1. Therefore, it is more efficient to solve a493 better conditioned system which has the same solution. For that purpose, we494 introduce a preconditioning matrix \textbf{Q} which is an approximation of495 \textbf{A} but much easier to invert than \textbf{A}, and solve the system:495 $i.e.$ its eigenvalues are closer to 1. Therefore, it is more efficient to solve 496 a better conditioned system which has the same solution. For that purpose, 497 we introduce a preconditioning matrix \textbf{Q} which is an approximation 498 of \textbf{A} but much easier to invert than \textbf{A}, and solve the system: 496 499 \begin{equation} \label{Eq_pmat} 497 500 \textbf{Q}^{-1} \textbf{A x} = \textbf{Q}^{-1} \textbf{b} … … 502 505 ${\langle{ \textbf{a} , \textbf{b}} \rangle}_Q = \langle{ \textbf{a} , \textbf{Q b}} \rangle$, and 503 506 if $\textbf{\~{b}} = \textbf{Q}^{-1}\;\textbf{b}$ and $\textbf{\~{A}} = \textbf{Q}^{-1}\;\textbf{A}$ 504 are substituted to \textbf{b} and \textbf{A} \citep{Madec 1988}.507 are substituted to \textbf{b} and \textbf{A} \citep{Madec_al_OM88}. 505 508 In \NEMO, \textbf{Q} is chosen as the diagonal of \textbf{ A}, i.e. the simplest form for 506 509 \textbf{Q} so that it can be easily inverted. In this case, the discrete formulation of … … 508 511 right hand side are computed independently from the solver used. 509 512 510 % -------------------------------------------------------------------------------------------------------------511 % FETI512 % -------------------------------------------------------------------------------------------------------------513 \subsection{FETI}514 \label{MISC_solfeti}515 516 FETI is a powerful solver that was developed by Marc Guyon \citep{Guyon_al_EP99,517 Guyon_al_CalPar99}. It has been converted from Fortan 77 to 90, but never518 successfully tested after that.519 520 Its main advantage is to save a lot of CPU time when compared to the SOR or PCG521 algorithms. However, its main drawback is that the solution is dependent on the522 domain decomposition chosen. Using a different number of processors, the solution523 is the same at the precision required, but not the same at the computer precision.524 This make it hard to debug.525 526 % -------------------------------------------------------------------------------------------------------------527 % Boundary Conditions --- Islands528 % -------------------------------------------------------------------------------------------------------------529 \subsection{Boundary Conditions --- Islands (\key{islands})}530 \label{MISC_solisl}531 532 The boundary condition used for both recommended solvers is that the time533 derivative of the barotropic streamfunction is zero along all the coastlines.534 When islands are present in the model domain, additional computations must535 be performed to determine the barotropic streamfunction with the correct boundary536 conditions. This is detailed below.537 538 The model does not have specialised code for islands. They must instead be539 identified to the solvers by the user via bathymetry information, i.e. the $mbathy$540 array should equal $-1$ over the first island, $-2$ over the second, ... , $-N$ over541 the $N^{th}$ island. The model determines the position of island grid-points and542 defines a closed contour around each island which is used to compute the circulation543 around each one. The closed contour is formed from the ocean grid-points which544 are the closest to the island.545 546 First, the island barotropic streamfunctions $\psi_n$ are computed using the SOR547 or PCG scheme (they are solutions of \eqref{Eq_solmat} with the right-hand side548 equal to zero and with $\psi_n = 1$ along the island $n$ and $\psi_n = 0$ along549 the other coastlines) (Note that specifying $1$ as boundary condition on an island550 for $\psi$ is equivalent to defining a specific right hand side for \eqref{Eq_solmat}).551 The requested precision of this computation can be very high since it is only552 performed once. The absolute precision, $epsisl$, and the maximum number of553 iterations allowed, $nmisl$, are the \textbf{namelist parameters} used for this554 computation. Their typical values are $epsisl = 10^{-10}$ and $nmisl = 4000$.555 Then the island matrix A is computed from \eqref{Eq_PE_psi_isl_matrix} and556 inverted. At each time step, $\psi_0$, the solution of \eqref{Eq_solmat} with $\psi = 0$557 along all coastlines, is computed using either SOR or PCG. (It should558 be noted that the first guess of this computation remains as usual except that559 $\partial_t \psi_0$ is used, instead of $\partial_t \psi$. Indeed we are560 computing $\partial_t \psi_0$ which is usually very different from $\partial_t \psi$.)561 Then, it is easy to find the time evolution of the barotropic streamfunction on each562 island and to deduce $\partial_t \psi$, and to use it to compute the surface pressure563 gradient. Note that the value of the barotropic streamfunction itself is also computed564 as the time integration of $\partial_t \psi$ for further diagnostics.565 566 513 % ================================================================ 567 514 % Diagnostics 568 515 % ================================================================ 569 \section{Diagnostics }516 \section{Diagnostics (DIA, IOM)} 570 517 \label{MISC_diag} 571 518 … … 577 524 578 525 %to be updated with Seb documentation on the IO 579 \amtcomment{ 526 580 527 The model outputs are of three types: the restart file, the output listing, 581 528 and the output file(s). The restart file is used internally by the code when … … 592 539 the $ocean.output$ file. The information is printed from within the code on the 593 540 logical unit $numout$. To locate these prints, use the UNIX command 594 "$grep -i numout^*$" in the source code directory. 595 596 In the standard configuration, the user will find the model results in two 597 output files containing values for every time-step where output is demanded: 598 a \textbf{VO} file containing all the three dimensional fields in logical unit 599 $numwvo$, and a \textbf{SO} file containing all the two dimensional (horizontal) 600 fields in logical unit $numwso$. These outputs are defined in the $diawri.F$ 601 routine. The default and user-selectable diagnostics are described in {\S}III-12-c. 602 603 The default output (for all output files apart from the restart file) is 32 bits binary 604 IEEE format, compatible with the Vairmer software (see the Climate Modelling 605 and global Change team WEB server at CERFACS: http://www.cerfacs.fr). 606 The model's reference directory also contains a visualisation tool based on 607 \textbf{NCAR Graphics} (http://ngwww.ucar.edu). If a user has access to the 608 NCAR software, he or she can copy the \textbf{LODMODEL/UTILS/OPADRA} 609 directory from the reference and, following the \textbf{README}, create 610 graphical outputs from the model's results. 611 } 541 "\textit{grep -i numout}" in the source code directory. 542 543 In the standard configuration, the user will find the model results in 544 NetCDF files containing mean values (or instantaneous values if 545 \key{diainstant} is defined) for every time-step where output is demanded. 546 These outputs are defined in the \mdl{diawri} module. 547 When defining \key{dimgout}, the output are written in DIMG format, 548 an IEEE output format. 549 550 Since version 3.2, an I/O server has been added which provides more 551 flexibility in the choice of the fields to be output as well as how the 552 writing work is distributed over the processors in massively parallel 553 computing. It is activated when \key{dimgout} is defined. 612 554 613 555 % ------------------------------------------------------------------------------------------------------------- … … 635 577 % On-line Floats trajectories 636 578 % ------------------------------------------------------------------------------------------------------------- 637 \subsection{On-line Floats trajectories }579 \subsection{On-line Floats trajectories (FLO)} 638 580 \label{FLO} 639 581 %--------------------------------------------namflo------------------------------------------------------- … … 669 611 - the meridional heat and salt transports and their decomposition (\mdl{diamfl}) 670 612 671 - the surface pressure (\mdl{diaspr}) 613 In addition, a series of diagnostics has been added in the \mdl{diaar5}. 614 They corresponds to outputs that are required for AR5 simulations 615 (see Section \ref{MISC_steric} below for one of them). 616 Activating those outputs requires to define the \key{diaar5} CPP key. 617 618 619 620 % ================================================================ 621 % Steric effect in sea surface height 622 % ================================================================ 623 \section{Steric effect in sea surface height} 624 \label{MISC_steric} 625 626 627 Changes in steric sea level are caused when changes in the density of the water 628 column imply an expansion or contraction of the column. It is essentially produced 629 through surface heating/cooling and to a lesser extent through non-linear effects of 630 the equation of state (cabbeling, thermobaricity...). 631 Non-Boussinesq models contain all ocean effects within the ocean acting 632 on the sea level. In particular, they include the steric effect. In contrast, 633 Boussinesq models, such as \NEMO, conserve volume, rather than mass, 634 and so do not properly represent expansion or contraction. The steric effect is 635 therefore not explicitely represented. 636 This approximation does not represent a serious error with respect to the flow field 637 calculated by the model \citep{Greatbatch_JGR94}, but extra attention is required 638 when investigating sea level, as steric changes are an important 639 contribution to local changes in sea level on seasonal and climatic time scales. 640 This is especially true for investigation into sea level rise due to global warming. 641 642 Fortunately, the steric contribution to the sea level consists of a spatially uniform 643 component that can be diagnosed by considering the mass budget of the world 644 ocean \citep{Greatbatch_JGR94}. 645 In order to better understand how global mean sea level evolves and thus how 646 the steric sea level can be diagnosed, we compare, in the following, the 647 non-Boussinesq and Boussinesq cases. 648 649 Let denote 650 $\mathcal{M}$ the total mass of liquid seawater ($\mathcal{M}=\int_D \rho dv$), 651 $\mathcal{V}$ the total volume of seawater ($\mathcal{V}=\int_D dv$), 652 $\mathcal{A}$ the total surface of the ocean ($\mathcal{A}=\int_S ds$), 653 $\bar{\rho}$ the global mean seawater (\textit{in situ}) density ($\bar{\rho}= 1/\mathcal{V} \int_D \rho \,dv$), and 654 $\bar{\eta}$ the global mean sea level ($\bar{\eta}=1/\mathcal{A}\int_S \eta \,ds$). 655 656 A non-Boussinesq fluid conserves mass. It satisfies the following relations: 657 \begin{equation} \label{Eq_MV_nBq} 658 \begin{split} 659 \mathcal{M} &= \mathcal{V} \;\bar{\rho} \\ 660 \mathcal{V} &= \mathcal{A} \;\bar{\eta} 661 \end{split} 662 \end{equation} 663 Temporal changes in total mass is obtained from the density conservation equation : 664 \begin{equation} \label{Eq_Co_nBq} 665 \frac{1}{e_3} \partial_t ( e_3\,\rho) + \nabla( \rho \, \textbf{U} ) = \left. \frac{\textit{emp}}{e_3}\right|_\textit{surface} 666 \end{equation} 667 where $\rho$ is the \textit{in situ} density, and \textit{emp} the surface mass 668 exchanges with the other media of the Earth system (atmosphere, sea-ice, land). 669 Its global averaged leads to the total mass change 670 \begin{equation} \label{Eq_Mass_nBq} 671 \partial_t \mathcal{M} = \mathcal{A} \;\overline{\textit{emp}} 672 \end{equation} 673 where $\overline{\textit{emp}}=\int_S \textit{emp}\,ds$ is the net mass flux 674 through the ocean surface. 675 Bringing \eqref{Eq_Mass_nBq} and the time derivative of \eqref{Eq_MV_nBq} 676 together leads to the evolution equation of the mean sea level 677 \begin{equation} \label{Eq_ssh_nBq} 678 \partial_t \bar{\eta} = \frac{\overline{\textit{emp}}}{ \bar{\rho}} 679 - \frac{\mathcal{V}}{\mathcal{A}} \;\frac{\partial_t \bar{\rho} }{\bar{\rho}} 680 \end{equation} 681 The first term in equation \eqref{Eq_ssh_nBq} alters sea level by adding or 682 subtracting mass from the ocean. 683 The second term arises from temporal changes in the global mean 684 density; $i.e.$ from steric effects. 685 686 In a Boussinesq fluid, $\rho$ is replaced by $\rho_o$ in all the equation except when $\rho$ 687 appears multiplied by the gravity ($i.e.$ in the hydrostatic balance of the primitive Equations). 688 In particular, the mass conservation equation, \eqref{Eq_Co_nBq}, degenerates into 689 the incompressibility equation: 690 \begin{equation} \label{Eq_Co_Bq} 691 \frac{1}{e_3} \partial_t ( e_3 ) + \nabla( \textbf{U} ) = \left. \frac{\textit{emp}}{\rho_o \,e_3}\right|_ \textit{surface} 692 \end{equation} 693 and the global average of this equation now gives the temporal change of the total volume, 694 \begin{equation} \label{Eq_V_Bq} 695 \partial_t \mathcal{V} = \mathcal{A} \;\frac{\overline{\textit{emp}}}{\rho_o} 696 \end{equation} 697 Only the volume is conserved, not mass, or, more precisely, the mass which is conserved is the 698 Boussinesq mass, $\mathcal{M}_o = \rho_o \mathcal{V}$. The total volume (or equivalently 699 the global mean sea level) is altered only by net volume fluxes across the ocean surface, 700 not by changes in mean mass of the ocean: the steric effect is missing in a Boussinesq fluid. 701 702 Nevertheless, following \citep{Greatbatch_JGR94}, the steric effect on the volume can be 703 diagnosed by considering the mass budget of the ocean. 704 The apparent changes in $\mathcal{M}$, mass of the ocean, which are not induced by surface 705 mass flux must be compensated by a spatially uniform change in the mean sea level due to 706 expansion/contraction of the ocean \citep{Greatbatch_JGR94}. In others words, the Boussinesq 707 mass, $\mathcal{M}_o$, can be related to $\mathcal{M}$, the total mass of the ocean seen 708 by the Boussinesq model, via the steric contribution to the sea level, $\eta_s$, a spatially 709 uniform variable, as follows: 710 \begin{equation} \label{Eq_M_Bq} 711 \mathcal{M}_o = \mathcal{M} + \rho_o \,\eta_s \,\mathcal{A} 712 \end{equation} 713 Any change in $\mathcal{M}$ which cannot be explained by the net mass flux through 714 the ocean surface is converted into a mean change in sea level. Introducing the total density 715 anomaly, $\mathcal{D}= \int_D d_a \,dv$, where $d_a= (\rho -\rho_o ) / \rho_o$ 716 is the density anomaly used in \NEMO (cf. \S\ref{TRA_eos}) in \eqref{Eq_M_Bq} 717 leads to a very simple form for the steric height: 718 \begin{equation} \label{Eq_steric_Bq} 719 \eta_s = - \frac{1}{\mathcal{A}} \mathcal{D} 720 \end{equation} 721 722 The above formulation of the steric height of a Boussinesq ocean requires four remarks. 723 First, one can be tempted to define $\rho_o$ as the initial value of $\mathcal{M}/\mathcal{V}$, 724 $i.e.$ set $\mathcal{D}_{t=0}=0$, so that the initial steric height is zero. We do not 725 recommend that. Indeed, in this case $\rho_o$ depends on the initial state of the ocean. 726 Since $\rho_o$ has a direct effect on the dynamics of the ocean (it appears in the pressure 727 gradient term of the momentum equation) it is definitively not a good idea when 728 inter-comparing experiments. 729 We better recommend to fixe once for all $\rho_o$ to $1035\;Kg\,m^{-3}$. This value is a 730 sensible choice for the reference density used in a Boussinesq ocean climate model since, 731 with the exception of only a small percentage of the ocean, density in the World Ocean 732 varies by no more than 2$\%$ from this value (\cite{Gill1982}, page 47). 733 734 Second, we have assumed here that the total ocean surface, $\mathcal{A}$, does not 735 change when the sea level is changing as it is the case in all global ocean GCMs 736 (wetting and drying of grid point is not allowed). 737 738 Third, the discretisation of \eqref{Eq_steric_Bq} depends on the type of free surface 739 which is considered. In the non linear free surface case, $i.e.$ \key{vvl} defined, it is 740 given by 741 \begin{equation} \label{Eq_discrete_steric_Bq} 742 \eta_s = - \frac{ \sum_{i,\,j,\,k} d_a\; e_{1t} e_{2t} e_{3t} } 743 { \sum_{i,\,j,\,k} e_{1t} e_{2t} e_{3t} } 744 \end{equation} 745 whereas in the linear free surface, the volume above the \textit{z=0} surface must be explicitly taken 746 into account to better approximate the total ocean mass and thus the steric sea level: 747 \begin{equation} \label{Eq_discrete_steric_Bq} 748 \eta_s = - \frac{ \sum_{i,\,j,\,k} d_a\; e_{1t}e_{2t}e_{3t} + \sum_{i,\,j} d_a\; e_{1t}e_{2t} \eta } 749 {\sum_{i,\,j,\,k} e_{1t}e_{2t}e_{3t} + \sum_{i,\,j} e_{1t}e_{2t} \eta } 750 \end{equation} 751 752 The fourth and last remark concerns the effective sea level and the presence of sea-ice. 753 In the real ocean, sea ice (and snow above it) depresses the liquid seawater through 754 its mass loading. This depression is a result of the mass of sea ice/snow system acting 755 on the liquid ocean. There is, however, no dynamical effect associated with these depressions 756 in the liquid ocean sea level, so that there are no associated ocean currents. Hence, the 757 dynamically relevant sea level is the effective sea level, $i.e.$ the sea level as if sea ice 758 (and snow) were converted to liquid seawater \citep{Campin_al_OM08}. However, 759 in the current version of \NEMO the sea-ice is levitating above the ocean without 760 mass exchanges between ice and ocean. Therefore the model effective sea level 761 is always given by $\eta + \eta_s$, whether or not there is sea ice present. 762 763 In AR5 outputs, the thermosteric sea level is demanded. It is steric sea level due to 764 changes in ocean density arising just from changes in temperature. It is given by: 765 \begin{equation} \label{Eq_thermosteric_Bq} 766 \eta_s = - \frac{1}{\mathcal{A}} \int_D d_a(T,S_o,p_o) \,dv 767 \end{equation} 768 where $S_o$ and $p_o$ are the initial salinity and pressure, respectively. 769 770 Both steric and thermosteric sea level are computed in \mdl{diaar5} which needs 771 the \key{diaar5} defined to be called. 772
Note: See TracChangeset
for help on using the changeset viewer.
