Changeset 11597 for NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex
- Timestamp:
- 2019-09-25T20:20:19+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex
r11596 r11597 53 53 MISC correspond to "extracting tendency terms" or "vorticity balance"?} 54 54 55 %% ================================================================================================= 55 56 \section{Sea surface height and diagnostic variables ($\eta$, $\zeta$, $\chi$, $w$)} 56 57 \label{sec:DYN_divcur_wzv} 57 58 58 %-------------------------------------------------------------------------------------------------------------- 59 % Horizontal divergence and relative vorticity 60 %-------------------------------------------------------------------------------------------------------------- 59 %% ================================================================================================= 61 60 \subsection[Horizontal divergence and relative vorticity (\textit{divcur.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{divcur})} 62 61 \label{subsec:DYN_divcur} … … 92 91 the nonlinear advection and of the vertical velocity respectively. 93 92 94 %-------------------------------------------------------------------------------------------------------------- 95 % Sea Surface Height evolution 96 %-------------------------------------------------------------------------------------------------------------- 93 %% ================================================================================================= 97 94 \subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} 98 95 \label{subsec:DYN_sshwzv} … … 152 149 (see \autoref{subsec:DOM_Num_Index_vertical}). 153 150 151 152 %% ================================================================================================= 154 153 \section{Coriolis and advection: vector invariant form} 155 154 \label{sec:DYN_adv_cor_vect} 156 %-----------------------------------------nam_dynadv----------------------------------------------------157 155 158 156 \begin{listing} … … 161 159 \label{lst:namdyn_adv} 162 160 \end{listing} 163 %-------------------------------------------------------------------------------------------------------------164 161 165 162 The vector invariant form of the momentum equations is the one most often used in … … 172 169 \autoref{chap:LBC}. 173 170 171 %% ================================================================================================= 174 172 \subsection[Vorticity term (\textit{dynvor.F90})]{Vorticity term (\protect\mdl{dynvor})} 175 173 \label{subsec:DYN_vor} 176 %------------------------------------------nam_dynvor----------------------------------------------------177 174 178 175 \begin{listing} … … 181 178 \label{lst:namdyn_vor} 182 179 \end{listing} 183 %-------------------------------------------------------------------------------------------------------------184 180 185 181 Options are defined through the \nam{dyn_vor}{dyn\_vor} namelist variables. … … 195 191 The vorticity terms are all computed in dedicated routines that can be found in the \mdl{dynvor} module. 196 192 197 %-------------------------------------------------------------198 193 % enstrophy conserving scheme 199 % -------------------------------------------------------------194 %% ================================================================================================= 200 195 \subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens})]{Enstrophy conserving scheme (\protect\np{ln_dynvor_ens}{ln\_dynvor\_ens})} 201 196 \label{subsec:DYN_vor_ens} … … 218 213 \end{equation} 219 214 220 %-------------------------------------------------------------221 215 % energy conserving scheme 222 % -------------------------------------------------------------216 %% ================================================================================================= 223 217 \subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene})]{Energy conserving scheme (\protect\np{ln_dynvor_ene}{ln\_dynvor\_ene})} 224 218 \label{subsec:DYN_vor_ene} … … 238 232 \end{equation} 239 233 240 %-------------------------------------------------------------241 234 % mix energy/enstrophy conserving scheme 242 % -------------------------------------------------------------235 %% ================================================================================================= 243 236 \subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix})]{Mixed energy/enstrophy conserving scheme (\protect\np{ln_dynvor_mix}{ln\_dynvor\_mix})} 244 237 \label{subsec:DYN_vor_mix} … … 263 256 \] 264 257 265 %-------------------------------------------------------------266 258 % energy and enstrophy conserving scheme 267 % -------------------------------------------------------------259 %% ================================================================================================= 268 260 \subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een})]{Energy and enstrophy conserving scheme (\protect\np{ln_dynvor_een}{ln\_dynvor\_een})} 269 261 \label{subsec:DYN_vor_een} … … 353 345 leading to a larger topostrophy of the flow \citep{barnier.madec.ea_OD06, penduff.le-sommer.ea_OS07}. 354 346 355 %-------------------------------------------------------------------------------------------------------------- 356 % Kinetic Energy Gradient term 357 %-------------------------------------------------------------------------------------------------------------- 347 %% ================================================================================================= 358 348 \subsection[Kinetic energy gradient term (\textit{dynkeg.F90})]{Kinetic energy gradient term (\protect\mdl{dynkeg})} 359 349 \label{subsec:DYN_keg} … … 373 363 \] 374 364 375 %-------------------------------------------------------------------------------------------------------------- 376 % Vertical advection term 377 %-------------------------------------------------------------------------------------------------------------- 365 %% ================================================================================================= 378 366 \subsection[Vertical advection term (\textit{dynzad.F90})]{Vertical advection term (\protect\mdl{dynzad})} 379 367 \label{subsec:DYN_zad} … … 400 388 an option which is only available with a TVD scheme (see \np{ln_traadv_tvd_zts}{ln\_traadv\_tvd\_zts} in \autoref{subsec:TRA_adv_tvd}). 401 389 390 391 %% ================================================================================================= 402 392 \section{Coriolis and advection: flux form} 403 393 \label{sec:DYN_adv_cor_flux} 404 %------------------------------------------nam_dynadv----------------------------------------------------405 406 %-------------------------------------------------------------------------------------------------------------407 394 408 395 Options are defined through the \nam{dyn_adv}{dyn\_adv} namelist variables. … … 413 400 no slip or partial slip boundary conditions are applied following \autoref{chap:LBC}. 414 401 415 %-------------------------------------------------------------------------------------------------------------- 416 % Coriolis plus curvature metric terms 417 %-------------------------------------------------------------------------------------------------------------- 402 %% ================================================================================================= 418 403 \subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})]{Coriolis plus curvature metric terms (\protect\mdl{dynvor})} 419 404 \label{subsec:DYN_cor_flux} … … 434 419 This term is evaluated using a leapfrog scheme, \ie\ the velocity is centred in time (\textit{now} velocity). 435 420 436 %-------------------------------------------------------------------------------------------------------------- 437 % Flux form Advection term 438 %-------------------------------------------------------------------------------------------------------------- 421 %% ================================================================================================= 439 422 \subsection[Flux form advection term (\textit{dynadv.F90})]{Flux form advection term (\protect\mdl{dynadv})} 440 423 \label{subsec:DYN_adv_flux} … … 467 450 and $uw$-points for $u$ and at the $f$-, $T$- and $vw$-points for $v$. 468 451 469 %-------------------------------------------------------------470 452 % 2nd order centred scheme 471 % -------------------------------------------------------------453 %% ================================================================================================= 472 454 \subsubsection[CEN2: $2^{nd}$ order centred scheme (\forcode{ln_dynadv_cen2})]{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln_dynadv_cen2}{ln\_dynadv\_cen2})} 473 455 \label{subsec:DYN_adv_cen2} … … 490 472 so $u$ and $v$ are the \emph{now} velocities. 491 473 492 %-------------------------------------------------------------493 474 % UBS scheme 494 % -------------------------------------------------------------475 %% ================================================================================================= 495 476 \subsubsection[UBS: Upstream Biased Scheme (\forcode{ln_dynadv_ubs})]{UBS: Upstream Biased Scheme (\protect\np{ln_dynadv_ubs}{ln\_dynadv\_ubs})} 496 477 \label{subsec:DYN_adv_ubs} … … 543 524 %%% 544 525 526 %% ================================================================================================= 545 527 \section[Hydrostatic pressure gradient (\textit{dynhpg.F90})]{Hydrostatic pressure gradient (\protect\mdl{dynhpg})} 546 528 \label{sec:DYN_hpg} 547 %------------------------------------------nam_dynhpg---------------------------------------------------548 529 549 530 \begin{listing} … … 552 533 \label{lst:namdyn_hpg} 553 534 \end{listing} 554 %-------------------------------------------------------------------------------------------------------------555 535 556 536 Options are defined through the \nam{dyn_hpg}{dyn\_hpg} namelist variables. … … 566 546 At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied. 567 547 568 %-------------------------------------------------------------------------------------------------------------- 569 % z-coordinate with full step 570 %-------------------------------------------------------------------------------------------------------------- 548 %% ================================================================================================= 571 549 \subsection[Full step $Z$-coordinate (\forcode{ln_dynhpg_zco})]{Full step $Z$-coordinate (\protect\np{ln_dynhpg_zco}{ln\_dynhpg\_zco})} 572 550 \label{subsec:DYN_hpg_zco} … … 611 589 \autoref{eq:DYN_hpg_zco} through the space and time variations of the vertical scale factor $e_{3w}$. 612 590 613 %-------------------------------------------------------------------------------------------------------------- 614 % z-coordinate with partial step 615 %-------------------------------------------------------------------------------------------------------------- 591 %% ================================================================================================= 616 592 \subsection[Partial step $Z$-coordinate (\forcode{ln_dynhpg_zps})]{Partial step $Z$-coordinate (\protect\np{ln_dynhpg_zps}{ln\_dynhpg\_zps})} 617 593 \label{subsec:DYN_hpg_zps} … … 632 608 module \mdl{zpsdhe} located in the TRA directory and described in \autoref{sec:TRA_zpshde}. 633 609 634 %-------------------------------------------------------------------------------------------------------------- 635 % s- and s-z-coordinates 636 %-------------------------------------------------------------------------------------------------------------- 610 %% ================================================================================================= 637 611 \subsection{$S$- and $Z$-$S$-coordinates} 638 612 \label{subsec:DYN_hpg_sco} … … 681 655 This method can provide a more accurate calculation of the horizontal pressure gradient than the standard scheme. 682 656 657 %% ================================================================================================= 683 658 \subsection{Ice shelf cavity} 684 659 \label{subsec:DYN_hpg_isf} … … 697 672 \autoref{subsec:DYN_hpg_sco}. 698 673 699 %-------------------------------------------------------------------------------------------------------------- 700 % Time-scheme 701 %-------------------------------------------------------------------------------------------------------------- 674 %% ================================================================================================= 702 675 \subsection[Time-scheme (\forcode{ln_dynhpg_imp})]{Time-scheme (\protect\np{ln_dynhpg_imp}{ln\_dynhpg\_imp})} 703 676 \label{subsec:DYN_hpg_imp} … … 758 731 This option is controlled by \np{nn_dynhpg_rst}{nn\_dynhpg\_rst}, a namelist parameter. 759 732 733 %% ================================================================================================= 760 734 \section[Surface pressure gradient (\textit{dynspg.F90})]{Surface pressure gradient (\protect\mdl{dynspg})} 761 735 \label{sec:DYN_spg} 762 %-----------------------------------------nam_dynspg----------------------------------------------------763 736 764 737 \begin{listing} … … 767 740 \label{lst:namdyn_spg} 768 741 \end{listing} 769 %------------------------------------------------------------------------------------------------------------770 742 771 743 Options are defined through the \nam{dyn_spg}{dyn\_spg} namelist variables. … … 795 767 As a consequence the update of the $next$ velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. 796 768 797 %-------------------------------------------------------------------------------------------------------------- 798 % Explicit free surface formulation 799 %-------------------------------------------------------------------------------------------------------------- 769 %% ================================================================================================= 800 770 \subsection[Explicit free surface (\forcode{ln_dynspg_exp})]{Explicit free surface (\protect\np{ln_dynspg_exp}{ln\_dynspg\_exp})} 801 771 \label{subsec:DYN_spg_exp} … … 821 791 Thus, nothing is done in the \mdl{dynspg\_exp} module. 822 792 823 %-------------------------------------------------------------------------------------------------------------- 824 % Split-explict free surface formulation 825 %-------------------------------------------------------------------------------------------------------------- 793 %% ================================================================================================= 826 794 \subsection[Split-explicit free surface (\forcode{ln_dynspg_ts})]{Split-explicit free surface (\protect\np{ln_dynspg_ts}{ln\_dynspg\_ts})} 827 795 \label{subsec:DYN_spg_ts} 828 %------------------------------------------namsplit-----------------------------------------------------------829 796 % 830 797 %\nlst{namsplit} 831 %-------------------------------------------------------------------------------------------------------------832 798 833 799 The split-explicit free surface formulation used in \NEMO\ (\np{ln_dynspg_ts}{ln\_dynspg\_ts} set to true), … … 1065 1031 %>>>>>=============== 1066 1032 1067 %-------------------------------------------------------------------------------------------------------------- 1068 % Filtered free surface formulation 1069 %-------------------------------------------------------------------------------------------------------------- 1033 %% ================================================================================================= 1070 1034 \subsection{Filtered free surface (\forcode{dynspg_flt?})} 1071 1035 \label{subsec:DYN_spg_fltp} … … 1093 1057 It is computed once and for all and applies to all ocean time steps. 1094 1058 1059 %% ================================================================================================= 1095 1060 \section[Lateral diffusion term and operators (\textit{dynldf.F90})]{Lateral diffusion term and operators (\protect\mdl{dynldf})} 1096 1061 \label{sec:DYN_ldf} 1097 %------------------------------------------nam_dynldf----------------------------------------------------1098 1062 1099 1063 \begin{listing} … … 1102 1066 \label{lst:namdyn_ldf} 1103 1067 \end{listing} 1104 %-------------------------------------------------------------------------------------------------------------1105 1068 1106 1069 Options are defined through the \nam{dyn_ldf}{dyn\_ldf} namelist variables. … … 1130 1093 } 1131 1094 1132 % Vertical diffusion term 1133 % External Forcing 1134 % Wetting and drying 1135 % Time evolution term 1095 %% ================================================================================================= 1096 \subsection[Iso-level laplacian (\forcode{ln_dynldf_lap})]{Iso-level laplacian operator (\protect\np{ln_dynldf_lap}{ln\_dynldf\_lap})} 1097 \label{subsec:DYN_ldf_lap} 1098 1099 For lateral iso-level diffusion, the discrete operator is: 1100 \begin{equation} 1101 \label{eq:DYN_ldf_lap} 1102 \left\{ 1103 \begin{aligned} 1104 D_u^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm} 1105 \;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta_j \left[ 1106 {A_f^{lm} \;e_{3f} \zeta } \right] \\ \\ 1107 D_v^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm} 1108 \;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta_i \left[ 1109 {A_f^{lm} \;e_{3f} \zeta } \right] 1110 \end{aligned} 1111 \right. 1112 \end{equation} 1113 1114 As explained in \autoref{subsec:MB_ldf}, 1115 this formulation (as the gradient of a divergence and curl of the vorticity) preserves symmetry and 1116 ensures a complete separation between the vorticity and divergence parts of the momentum diffusion. 1117 1118 %% ================================================================================================= 1119 \subsection[Rotated laplacian (\forcode{ln_dynldf_iso})]{Rotated laplacian operator (\protect\np{ln_dynldf_iso}{ln\_dynldf\_iso})} 1120 \label{subsec:DYN_ldf_iso} 1121 1122 A rotation of the lateral momentum diffusion operator is needed in several cases: 1123 for iso-neutral diffusion in the $z$-coordinate (\np[=.true.]{ln_dynldf_iso}{ln\_dynldf\_iso}) and 1124 for either iso-neutral (\np[=.true.]{ln_dynldf_iso}{ln\_dynldf\_iso}) or 1125 geopotential (\np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}) diffusion in the $s$-coordinate. 1126 In the partial step case, coordinates are horizontal except at the deepest level and 1127 no rotation is performed when \np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}. 1128 The diffusion operator is defined simply as the divergence of down gradient momentum fluxes on 1129 each momentum component. 1130 It must be emphasized that this formulation ignores constraints on the stress tensor such as symmetry. 1131 The resulting discrete representation is: 1132 \begin{equation} 1133 \label{eq:DYN_ldf_iso} 1134 \begin{split} 1135 D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\ 1136 & \left\{\quad {\delta_{i+1/2} \left[ {A_T^{lm} \left( 1137 {\frac{e_{2t} \; e_{3t} }{e_{1t} } \,\delta_{i}[u] 1138 -e_{2t} \; r_{1t} \,\overline{\overline {\delta_{k+1/2}[u]}}^{\,i,\,k}} 1139 \right)} \right]} \right. \\ 1140 & \qquad +\ \delta_j \left[ {A_f^{lm} \left( {\frac{e_{1f}\,e_{3f} }{e_{2f} 1141 }\,\delta_{j+1/2} [u] - e_{1f}\, r_{2f} 1142 \,\overline{\overline {\delta_{k+1/2} [u]}} ^{\,j+1/2,\,k}} 1143 \right)} \right] \\ 1144 &\qquad +\ \delta_k \left[ {A_{uw}^{lm} \left( {-e_{2u} \, r_{1uw} \,\overline{\overline 1145 {\delta_{i+1/2} [u]}}^{\,i+1/2,\,k+1/2} } 1146 \right.} \right. \\ 1147 & \ \qquad \qquad \qquad \quad\ 1148 - e_{1u} \, r_{2uw} \,\overline{\overline {\delta_{j+1/2} [u]}} ^{\,j,\,k+1/2} \\ 1149 & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\ 1150 +\frac{e_{1u}\, e_{2u} }{e_{3uw} }\,\left( {r_{1uw}^2+r_{2uw}^2} 1151 \right)\,\delta_{k+1/2} [u]} \right)} \right]\;\;\;} \right\} \\ \\ 1152 D_v^{l\textbf{V}} &= \frac{1}{e_{1v} \, e_{2v} \, e_{3v} } \\ 1153 & \left\{\quad {\delta_{i+1/2} \left[ {A_f^{lm} \left( 1154 {\frac{e_{2f} \; e_{3f} }{e_{1f} } \,\delta_{i+1/2}[v] 1155 -e_{2f} \; r_{1f} \,\overline{\overline {\delta_{k+1/2}[v]}}^{\,i+1/2,\,k}} 1156 \right)} \right]} \right. \\ 1157 & \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1t}\,e_{3t} }{e_{2t} 1158 }\,\delta_{j} [v] - e_{1t}\, r_{2t} 1159 \,\overline{\overline {\delta_{k+1/2} [v]}} ^{\,j,\,k}} 1160 \right)} \right] \\ 1161 & \qquad +\ \delta_k \left[ {A_{vw}^{lm} \left( {-e_{2v} \, r_{1vw} \,\overline{\overline 1162 {\delta_{i+1/2} [v]}}^{\,i+1/2,\,k+1/2} }\right.} \right. \\ 1163 & \ \qquad \qquad \qquad \quad\ 1164 - e_{1v} \, r_{2vw} \,\overline{\overline {\delta_{j+1/2} [v]}} ^{\,j+1/2,\,k+1/2} \\ 1165 & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\ 1166 +\frac{e_{1v}\, e_{2v} }{e_{3vw} }\,\left( {r_{1vw}^2+r_{2vw}^2} 1167 \right)\,\delta_{k+1/2} [v]} \right)} \right]\;\;\;} \right\} 1168 \end{split} 1169 \end{equation} 1170 where $r_1$ and $r_2$ are the slopes between the surface along which the diffusion operator acts and 1171 the surface of computation ($z$- or $s$-surfaces). 1172 The way these slopes are evaluated is given in the lateral physics chapter (\autoref{chap:LDF}). 1173 1174 %% ================================================================================================= 1175 \subsection[Iso-level bilaplacian (\forcode{ln_dynldf_bilap})]{Iso-level bilaplacian operator (\protect\np{ln_dynldf_bilap}{ln\_dynldf\_bilap})} 1176 \label{subsec:DYN_ldf_bilap} 1177 1178 The lateral fourth order operator formulation on momentum is obtained by applying \autoref{eq:DYN_ldf_lap} twice. 1179 It requires an additional assumption on boundary conditions: 1180 the first derivative term normal to the coast depends on the free or no-slip lateral boundary conditions chosen, 1181 while the third derivative terms normal to the coast are set to zero (see \autoref{chap:LBC}). 1182 %%% 1183 \gmcomment{add a remark on the the change in the position of the coefficient} 1184 %%% 1185 1186 %% ================================================================================================= 1187 \section[Vertical diffusion term (\textit{dynzdf.F90})]{Vertical diffusion term (\protect\mdl{dynzdf})} 1188 \label{sec:DYN_zdf} 1189 1190 Options are defined through the \nam{zdf}{zdf} namelist variables. 1191 The large vertical diffusion coefficient found in the surface mixed layer together with high vertical resolution implies that in the case of explicit time stepping there would be too restrictive a constraint on the time step. 1192 Two time stepping schemes can be used for the vertical diffusion term: 1193 $(a)$ a forward time differencing scheme 1194 (\np[=.true.]{ln_zdfexp}{ln\_zdfexp}) using a time splitting technique (\np{nn_zdfexp}{nn\_zdfexp} $>$ 1) or 1195 $(b)$ a backward (or implicit) time differencing scheme (\np[=.false.]{ln_zdfexp}{ln\_zdfexp}) 1196 (see \autoref{chap:TD}). 1197 Note that namelist variables \np{ln_zdfexp}{ln\_zdfexp} and \np{nn_zdfexp}{nn\_zdfexp} apply to both tracers and dynamics. 1198 1199 The formulation of the vertical subgrid scale physics is the same whatever the vertical coordinate is. 1200 The vertical diffusion operators given by \autoref{eq:MB_zdf} take the following semi-discrete space form: 1201 \[ 1202 % \label{eq:DYN_zdf} 1203 \left\{ 1204 \begin{aligned} 1205 D_u^{vm} &\equiv \frac{1}{e_{3u}} \ \delta_k \left[ \frac{A_{uw}^{vm} }{e_{3uw} } 1206 \ \delta_{k+1/2} [\,u\,] \right] \\ 1207 \\ 1208 D_v^{vm} &\equiv \frac{1}{e_{3v}} \ \delta_k \left[ \frac{A_{vw}^{vm} }{e_{3vw} } 1209 \ \delta_{k+1/2} [\,v\,] \right] 1210 \end{aligned} 1211 \right. 1212 \] 1213 where $A_{uw}^{vm} $ and $A_{vw}^{vm} $ are the vertical eddy viscosity and diffusivity coefficients. 1214 The way these coefficients are evaluated depends on the vertical physics used (see \autoref{chap:ZDF}). 1215 1216 The surface boundary condition on momentum is the stress exerted by the wind. 1217 At the surface, the momentum fluxes are prescribed as the boundary condition on 1218 the vertical turbulent momentum fluxes, 1219 \begin{equation} 1220 \label{eq:DYN_zdf_sbc} 1221 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1} 1222 = \frac{1}{\rho_o} \binom{\tau_u}{\tau_v } 1223 \end{equation} 1224 where $\left( \tau_u ,\tau_v \right)$ are the two components of the wind stress vector in 1225 the (\textbf{i},\textbf{j}) coordinate system. 1226 The high mixing coefficients in the surface mixed layer ensure that the surface wind stress is distributed in 1227 the vertical over the mixed layer depth. 1228 If the vertical mixing coefficient is small (when no mixed layer scheme is used) 1229 the surface stress enters only the top model level, as a body force. 1230 The surface wind stress is calculated in the surface module routines (SBC, see \autoref{chap:SBC}). 1231 1232 The turbulent flux of momentum at the bottom of the ocean is specified through a bottom friction parameterisation 1233 (see \autoref{sec:ZDF_drg}) 1234 1235 %% ================================================================================================= 1236 \section{External forcings} 1237 \label{sec:DYN_forcing} 1238 1239 Besides the surface and bottom stresses (see the above section) 1240 which are introduced as boundary conditions on the vertical mixing, 1241 three other forcings may enter the dynamical equations by affecting the surface pressure gradient. 1242 1243 (1) When \np[=.true.]{ln_apr_dyn}{ln\_apr\_dyn} (see \autoref{sec:SBC_apr}), 1244 the atmospheric pressure is taken into account when computing the surface pressure gradient. 1245 1246 (2) When \np[=.true.]{ln_tide_pot}{ln\_tide\_pot} and \np[=.true.]{ln_tide}{ln\_tide} (see \autoref{sec:SBC_tide}), 1247 the tidal potential is taken into account when computing the surface pressure gradient. 1248 1249 (3) When \np[=2]{nn_ice_embd}{nn\_ice\_embd} and LIM or CICE is used 1250 (\ie\ when the sea-ice is embedded in the ocean), 1251 the snow-ice mass is taken into account when computing the surface pressure gradient. 1252 1253 1254 \gmcomment{ missing : the lateral boundary condition !!! another external forcing 1255 } 1256 1257 %% ================================================================================================= 1258 \section{Wetting and drying } 1259 \label{sec:DYN_wetdry} 1260 1261 There are two main options for wetting and drying code (wd): 1262 (a) an iterative limiter (il) and (b) a directional limiter (dl). 1263 The directional limiter is based on the scheme developed by \cite{warner.defne.ea_CG13} for RO 1264 MS 1265 which was in turn based on ideas developed for POM by \cite{oey_OM06}. The iterative 1266 limiter is a new scheme. The iterative limiter is activated by setting $\mathrm{ln\_wd\_il} = \mathrm{.true.}$ 1267 and $\mathrm{ln\_wd\_dl} = \mathrm{.false.}$. The directional limiter is activated 1268 by setting $\mathrm{ln\_wd\_dl} = \mathrm{.true.}$ and $\mathrm{ln\_wd\_il} = \mathrm{.false.}$. 1269 1270 \begin{listing} 1271 \nlst{namwad} 1272 \caption{\forcode{&namwad}} 1273 \label{lst:namwad} 1274 \end{listing} 1275 1276 The following terminology is used. The depth of the topography (positive downwards) 1277 at each $(i,j)$ point is the quantity stored in array $\mathrm{ht\_wd}$ in the \NEMO\ code. 1278 The height of the free surface (positive upwards) is denoted by $ \mathrm{ssh}$. Given the sign 1279 conventions used, the water depth, $h$, is the height of the free surface plus the depth of the 1280 topography (i.e. $\mathrm{ssh} + \mathrm{ht\_wd}$). 1281 1282 Both wd schemes take all points in the domain below a land elevation of $\mathrm{rn\_wdld}$ to be 1283 covered by water. They require the topography specified with a model 1284 configuration to have negative depths at points where the land is higher than the 1285 topography's reference sea-level. The vertical grid in \NEMO\ is normally computed relative to an 1286 initial state with zero sea surface height elevation. 1287 The user can choose to compute the vertical grid and heights in the model relative to 1288 a non-zero reference height for the free surface. This choice affects the calculation of the metrics and depths 1289 (i.e. the $\mathrm{e3t\_0, ht\_0}$ etc. arrays). 1290 1291 Points where the water depth is less than $\mathrm{rn\_wdmin1}$ are interpreted as ``dry''. 1292 $\mathrm{rn\_wdmin1}$ is usually chosen to be of order $0.05$m but extreme topographies 1293 with very steep slopes require larger values for normal choices of time-step. Surface fluxes 1294 are also switched off for dry cells to prevent freezing, boiling etc. of very thin water layers. 1295 The fluxes are tappered down using a $\mathrm{tanh}$ weighting function 1296 to no flux as the dry limit $\mathrm{rn\_wdmin1}$ is approached. Even wet cells can be very shallow. 1297 The depth at which to start tapering is controlled by the user by setting $\mathrm{rn\_wd\_sbcdep}$. 1298 The fraction $(<1)$ of sufrace fluxes to use at this depth is set by $\mathrm{rn\_wd\_sbcfra}$. 1299 1300 Both versions of the code have been tested in six test cases provided in the WAD\_TEST\_CASES configuration 1301 and in ``realistic'' configurations covering parts of the north-west European shelf. 1302 All these configurations have used pure sigma coordinates. It is expected that 1303 the wetting and drying code will work in domains with more general s-coordinates provided 1304 the coordinates are pure sigma in the region where wetting and drying actually occurs. 1305 1306 The next sub-section descrbies the directional limiter and the following sub-section the iterative limiter. 1307 The final sub-section covers some additional considerations that are relevant to both schemes. 1308 1309 1310 % Iterative limiters 1311 %% ================================================================================================= 1312 \subsection[Directional limiter (\textit{wet\_dry.F90})]{Directional limiter (\mdl{wet\_dry})} 1313 \label{subsec:DYN_wd_directional_limiter} 1314 1315 The principal idea of the directional limiter is that 1316 water should not be allowed to flow out of a dry tracer cell (i.e. one whose water depth is less than \np{rn_wdmin1}{rn\_wdmin1}). 1317 1318 All the changes associated with this option are made to the barotropic solver for the non-linear 1319 free surface code within dynspg\_ts. 1320 On each barotropic sub-step the scheme determines the direction of the flow across each face of all the tracer cells 1321 and sets the flux across the face to zero when the flux is from a dry tracer cell. This prevents cells 1322 whose depth is rn\_wdmin1 or less from drying out further. The scheme does not force $h$ (the water depth) at tracer cells 1323 to be at least the minimum depth and hence is able to conserve mass / volume. 1324 1325 The flux across each $u$-face of a tracer cell is multiplied by a factor zuwdmask (an array which depends on ji and jj). 1326 If the user sets \np[=.false.]{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp} then zuwdmask is 1 when the 1327 flux is from a cell with water depth greater than \np{rn_wdmin1}{rn\_wdmin1} and 0 otherwise. If the user sets 1328 \np[=.true.]{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp} the flux across the face is ramped down as the water depth decreases 1329 from 2 * \np{rn_wdmin1}{rn\_wdmin1} to \np{rn_wdmin1}{rn\_wdmin1}. The use of this ramp reduced grid-scale noise in idealised test cases. 1330 1331 At the point where the flux across a $u$-face is multiplied by zuwdmask , we have chosen 1332 also to multiply the corresponding velocity on the ``now'' step at that face by zuwdmask. We could have 1333 chosen not to do that and to allow fairly large velocities to occur in these ``dry'' cells. 1334 The rationale for setting the velocity to zero is that it is the momentum equations that are being solved 1335 and the total momentum of the upstream cell (treating it as a finite volume) should be considered 1336 to be its depth times its velocity. This depth is considered to be zero at ``dry'' $u$-points consistent with its 1337 treatment in the calculation of the flux of mass across the cell face. 1338 1339 1340 \cite{warner.defne.ea_CG13} state that in their scheme the velocity masks at the cell faces for the baroclinic 1341 timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic sub-steps is respectively less than 1342 or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer 1343 fields (tracers independent of $x$, $y$ and $z$). Our scheme conserves constant tracers because 1344 the velocities used at the tracer cell faces on the baroclinic timesteps are carefully calculated by dynspg\_ts 1345 to equal their mean value during the barotropic steps. If the user sets \np[=.true.]{ln_wd_dl_bc}{ln\_wd\_dl\_bc}, the 1346 baroclinic velocities are also multiplied by a suitably weighted average of zuwdmask. 1347 1348 % Iterative limiters 1349 1350 %% ================================================================================================= 1351 \subsection[Iterative limiter (\textit{wet\_dry.F90})]{Iterative limiter (\mdl{wet\_dry})} 1352 \label{subsec:DYN_wd_iterative_limiter} 1353 1354 %% ================================================================================================= 1355 \subsubsection[Iterative flux limiter (\textit{wet\_dry.F90})]{Iterative flux limiter (\mdl{wet\_dry})} 1356 \label{subsec:DYN_wd_il_spg_limiter} 1357 1358 The iterative limiter modifies the fluxes across the faces of cells that are either already ``dry'' 1359 or may become dry within the next time-step using an iterative method. 1360 1361 The flux limiter for the barotropic flow (devised by Hedong Liu) can be understood as follows: 1362 1363 The continuity equation for the total water depth in a column 1364 \begin{equation} 1365 \label{eq:DYN_wd_continuity} 1366 \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 . 1367 \end{equation} 1368 can be written in discrete form as 1369 1370 \begin{align} 1371 \label{eq:DYN_wd_continuity_2} 1372 \frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) ) 1373 &= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j} + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\ 1374 &= \mathrm{zzflx}_{i,j} . 1375 \end{align} 1376 1377 In the above $h$ is the depth of the water in the column at point $(i,j)$, 1378 $\mathrm{flxu}_{i+1,j}$ is the flux out of the ``eastern'' face of the cell and 1379 $\mathrm{flxv}_{i,j+1}$ the flux out of the ``northern'' face of the cell; $t_{n+1}$ is 1380 the new timestep, $t_e$ is the old timestep (either $t_b$ or $t_n$) and $ \Delta t = 1381 t_{n+1} - t_e$; $e_1 e_2$ is the area of the tracer cells centred at $(i,j)$ and 1382 $\mathrm{zzflx}$ is the sum of the fluxes through all the faces. 1383 1384 The flux limiter splits the flux $\mathrm{zzflx}$ into fluxes that are out of the cell 1385 (zzflxp) and fluxes that are into the cell (zzflxn). Clearly 1386 1387 \begin{equation} 1388 \label{eq:DYN_wd_zzflx_p_n_1} 1389 \mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} . 1390 \end{equation} 1391 1392 The flux limiter iteratively adjusts the fluxes $\mathrm{flxu}$ and $\mathrm{flxv}$ until 1393 none of the cells will ``dry out''. To be precise the fluxes are limited until none of the 1394 cells has water depth less than $\mathrm{rn\_wdmin1}$ on step $n+1$. 1395 1396 Let the fluxes on the $m$th iteration step be denoted by $\mathrm{flxu}^{(m)}$ and 1397 $\mathrm{flxv}^{(m)}$. Then the adjustment is achieved by seeking a set of coefficients, 1398 $\mathrm{zcoef}_{i,j}^{(m)}$ such that: 1399 1400 \begin{equation} 1401 \label{eq:DYN_wd_continuity_coef} 1402 \begin{split} 1403 \mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\ 1404 \mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j} 1405 \end{split} 1406 \end{equation} 1407 1408 where the coefficients are $1.0$ generally but can vary between $0.0$ and $1.0$ around 1409 cells that would otherwise dry. 1410 1411 The iteration is initialised by setting 1412 1413 \begin{equation} 1414 \label{eq:DYN_wd_zzflx_initial} 1415 \mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} . 1416 \end{equation} 1417 1418 The fluxes out of cell $(i,j)$ are updated at the $m+1$th iteration if the depth of the 1419 cell on timestep $t_e$, namely $h_{i,j}(t_e)$, is less than the total flux out of the cell 1420 times the timestep divided by the cell area. Using (\autoref{eq:DYN_wd_continuity_2}) this 1421 condition is 1422 1423 \begin{equation} 1424 \label{eq:DYN_wd_continuity_if} 1425 h_{i,j}(t_e) - \mathrm{rn\_wdmin1} < \frac{\Delta t}{e_1 e_2} ( \mathrm{zzflxp}^{(m)}_{i,j} + \mathrm{zzflxn}^{(m)}_{i,j} ) . 1426 \end{equation} 1427 1428 Rearranging (\autoref{eq:DYN_wd_continuity_if}) we can obtain an expression for the maximum 1429 outward flux that can be allowed and still maintain the minimum wet depth: 1430 1431 \begin{equation} 1432 \label{eq:DYN_wd_max_flux} 1433 \begin{split} 1434 \mathrm{zzflxp}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2}) \frac{e_1 e_2}{\Delta t} \phantom{]} \\ 1435 \phantom{[} & - \mathrm{zzflxn}^{(m)}_{i,j} \Big] 1436 \end{split} 1437 \end{equation} 1438 1439 Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\itshape [Q: Why is 1440 this necessary/desirable?]}. Substituting from (\autoref{eq:DYN_wd_continuity_coef}) gives an 1441 expression for the coefficient needed to multiply the outward flux at this cell in order 1442 to avoid drying. 1443 1444 \begin{equation} 1445 \label{eq:DYN_wd_continuity_nxtcoef} 1446 \begin{split} 1447 \mathrm{zcoef}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2}) \frac{e_1 e_2}{\Delta t} \phantom{]} \\ 1448 \phantom{[} & - \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} } 1449 \end{split} 1450 \end{equation} 1451 1452 Only the outward flux components are altered but, of course, outward fluxes from one cell 1453 are inward fluxes to adjacent cells and the balance in these cells may need subsequent 1454 adjustment; hence the iterative nature of this scheme. Note, for example, that the flux 1455 across the ``eastern'' face of the $(i,j)$th cell is only updated at the $m+1$th iteration 1456 if that flux at the $m$th iteration is out of the $(i,j)$th cell. If that is the case then 1457 the flux across that face is into the $(i+1,j)$ cell and that flux will not be updated by 1458 the calculation for the $(i+1,j)$th cell. In this sense the updates to the fluxes across 1459 the faces of the cells do not ``compete'' (they do not over-write each other) and one 1460 would expect the scheme to converge relatively quickly. The scheme is flux based so 1461 conserves mass. It also conserves constant tracers for the same reason that the 1462 directional limiter does. 1463 1464 1465 % Surface pressure gradients 1466 %% ================================================================================================= 1467 \subsubsection[Modification of surface pressure gradients (\textit{dynhpg.F90})]{Modification of surface pressure gradients (\mdl{dynhpg})} 1468 \label{subsec:DYN_wd_il_spg} 1469 1470 At ``dry'' points the water depth is usually close to $\mathrm{rn\_wdmin1}$. If the 1471 topography is sloping at these points the sea-surface will have a similar slope and there 1472 will hence be very large horizontal pressure gradients at these points. The WAD modifies 1473 the magnitude but not the sign of the surface pressure gradients (zhpi and zhpj) at such 1474 points by mulitplying them by positive factors (zcpx and zcpy respectively) that lie 1475 between $0$ and $1$. 1476 1477 We describe how the scheme works for the ``eastward'' pressure gradient, zhpi, calculated 1478 at the $(i,j)$th $u$-point. The scheme uses the ht\_wd depths and surface heights at the 1479 neighbouring $(i+1,j)$ and $(i,j)$ tracer points. zcpx is calculated using two logicals 1480 variables, $\mathrm{ll\_tmp1}$ and $\mathrm{ll\_tmp2}$ which are evaluated for each grid 1481 column. The three possible combinations are illustrated in \autoref{fig:DYN_WAD_dynhpg}. 1482 1483 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1484 \begin{figure}[!ht] 1485 \centering 1486 \includegraphics[width=0.66\textwidth]{Fig_WAD_dynhpg} 1487 \caption[Combinations controlling the limiting of the horizontal pressure gradient in 1488 wetting and drying regimes]{ 1489 Three possible combinations of the logical variables controlling the 1490 limiting of the horizontal pressure gradient in wetting and drying regimes} 1491 \label{fig:DYN_WAD_dynhpg} 1492 \end{figure} 1493 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1494 1495 The first logical, $\mathrm{ll\_tmp1}$, is set to true if and only if the water depth at 1496 both neighbouring points is greater than $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ and 1497 the minimum height of the sea surface at the two points is greater than the maximum height 1498 of the topography at the two points: 1499 1500 \begin{equation} 1501 \label{eq:DYN_ll_tmp1} 1502 \begin{split} 1503 \mathrm{ll\_tmp1} = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\ 1504 & \quad \mathrm{MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj))\ .and.} \\ 1505 & \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\ 1506 & \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\ 1507 & \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 } 1508 \end{split} 1509 \end{equation} 1510 1511 The second logical, $\mathrm{ll\_tmp2}$, is set to true if and only if the maximum height 1512 of the sea surface at the two points is greater than the maximum height of the topography 1513 at the two points plus $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ 1514 1515 \begin{equation} 1516 \label{eq:DYN_ll_tmp2} 1517 \begin{split} 1518 \mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\ 1519 & \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\ 1520 & \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) . 1521 \end{split} 1522 \end{equation} 1523 1524 If $\mathrm{ll\_tmp1}$ is true then the surface pressure gradient, zhpi at the $(i,j)$ 1525 point is unmodified. If both logicals are false zhpi is set to zero. 1526 1527 If $\mathrm{ll\_tmp1}$ is true and $\mathrm{ll\_tmp2}$ is false then the surface pressure 1528 gradient is multiplied through by zcpx which is the absolute value of the difference in 1529 the water depths at the two points divided by the difference in the surface heights at the 1530 two points. Thus the sign of the sea surface height gradient is retained but the magnitude 1531 of the pressure force is determined by the difference in water depths rather than the 1532 difference in surface height between the two points. Note that dividing by the difference 1533 between the sea surface heights can be problematic if the heights approach parity. An 1534 additional condition is applied to $\mathrm{ ll\_tmp2 }$ to ensure it is .false. in such 1535 conditions. 1536 1537 %% ================================================================================================= 1538 \subsubsection[Additional considerations (\textit{usrdef\_zgr.F90})]{Additional considerations (\mdl{usrdef\_zgr})} 1539 \label{subsec:DYN_WAD_additional} 1540 1541 In the very shallow water where wetting and drying occurs the parametrisation of 1542 bottom drag is clearly very important. In order to promote stability 1543 it is sometimes useful to calculate the bottom drag using an implicit time-stepping approach. 1544 1545 Suitable specifcation of the surface heat flux in wetting and drying domains in forced and 1546 coupled simulations needs further consideration. In order to prevent freezing or boiling 1547 in uncoupled integrations the net surface heat fluxes need to be appropriately limited. 1548 1549 % The WAD test cases 1550 %% ================================================================================================= 1551 \subsection[The WAD test cases (\textit{usrdef\_zgr.F90})]{The WAD test cases (\mdl{usrdef\_zgr})} 1552 \label{subsec:DYN_WAD_test_cases} 1553 1554 See the WAD tests MY\_DOC documention for details of the WAD test cases. 1555 1556 1557 1558 %% ================================================================================================= 1559 \section[Time evolution term (\textit{dynnxt.F90})]{Time evolution term (\protect\mdl{dynnxt})} 1560 \label{sec:DYN_nxt} 1561 1562 1563 Options are defined through the \nam{dom}{dom} namelist variables. 1564 The general framework for dynamics time stepping is a leap-frog scheme, 1565 \ie\ a three level centred time scheme associated with an Asselin time filter (cf. \autoref{chap:TD}). 1566 The scheme is applied to the velocity, except when 1567 using the flux form of momentum advection (cf. \autoref{sec:DYN_adv_cor_flux}) 1568 in the variable volume case (\texttt{vvl?} defined), 1569 where it has to be applied to the thickness weighted velocity (see \autoref{sec:SCOORD_momentum}) 1570 1571 $\bullet$ vector invariant form or linear free surface 1572 (\np[=.true.]{ln_dynhpg_vec}{ln\_dynhpg\_vec} ; \texttt{vvl?} not defined): 1573 \[ 1574 % \label{eq:DYN_nxt_vec} 1575 \left\{ 1576 \begin{aligned} 1577 &u^{t+\rdt} = u_f^{t-\rdt} + 2\rdt \ \text{RHS}_u^t \\ 1578 &u_f^t \;\quad = u^t+\gamma \,\left[ {u_f^{t-\rdt} -2u^t+u^{t+\rdt}} \right] 1579 \end{aligned} 1580 \right. 1581 \] 1582 1583 $\bullet$ flux form and nonlinear free surface 1584 (\np[=.false.]{ln_dynhpg_vec}{ln\_dynhpg\_vec} ; \texttt{vvl?} defined): 1585 \[ 1586 % \label{eq:DYN_nxt_flux} 1587 \left\{ 1588 \begin{aligned} 1589 &\left(e_{3u}\,u\right)^{t+\rdt} = \left(e_{3u}\,u\right)_f^{t-\rdt} + 2\rdt \; e_{3u} \;\text{RHS}_u^t \\ 1590 &\left(e_{3u}\,u\right)_f^t \;\quad = \left(e_{3u}\,u\right)^t 1591 +\gamma \,\left[ {\left(e_{3u}\,u\right)_f^{t-\rdt} -2\left(e_{3u}\,u\right)^t+\left(e_{3u}\,u\right)^{t+\rdt}} \right] 1592 \end{aligned} 1593 \right. 1594 \] 1595 where RHS is the right hand side of the momentum equation, 1596 the subscript $f$ denotes filtered values and $\gamma$ is the Asselin coefficient. 1597 $\gamma$ is initialized as \np{nn_atfp}{nn\_atfp} (namelist parameter). 1598 Its default value is \np[=10.e-3]{nn_atfp}{nn\_atfp}. 1599 In both cases, the modified Asselin filter is not applied since perfect conservation is not an issue for 1600 the momentum equations. 1601 1602 Note that with the filtered free surface, 1603 the update of the \textit{after} velocities is done in the \mdl{dynsp\_flt} module, 1604 and only array swapping and Asselin filtering is done in \mdl{dynnxt}. 1605 1606 \onlyinsubfile{\input{../../global/epilogue}} 1607 1608 \end{document}
Note: See TracChangeset
for help on using the changeset viewer.