Changeset 11596 for NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex
- Timestamp:
- 2019-09-25T19:06:37+02:00 (5 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex
r11584 r11596 2 2 3 3 \begin{document} 4 % ================================================================5 % Chapter ——— Ocean Dynamics (DYN)6 % ================================================================7 4 \chapter{Ocean Dynamics (DYN)} 8 5 \label{chap:DYN} … … 56 53 MISC correspond to "extracting tendency terms" or "vorticity balance"?} 57 54 58 % ================================================================59 % Sea Surface Height evolution & Diagnostics variables60 % ================================================================61 55 \section{Sea surface height and diagnostic variables ($\eta$, $\zeta$, $\chi$, $w$)} 62 56 \label{sec:DYN_divcur_wzv} … … 158 152 (see \autoref{subsec:DOM_Num_Index_vertical}). 159 153 160 161 % ================================================================162 % Coriolis and Advection terms: vector invariant form163 % ================================================================164 154 \section{Coriolis and advection: vector invariant form} 165 155 \label{sec:DYN_adv_cor_vect} … … 182 172 \autoref{chap:LBC}. 183 173 184 % -------------------------------------------------------------------------------------------------------------185 % Vorticity term186 % -------------------------------------------------------------------------------------------------------------187 174 \subsection[Vorticity term (\textit{dynvor.F90})]{Vorticity term (\protect\mdl{dynvor})} 188 175 \label{subsec:DYN_vor} … … 314 301 \end{equation} 315 302 316 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>317 303 \begin{figure}[!ht] 318 304 \centering … … 414 400 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}). 415 401 416 417 % ================================================================418 % Coriolis and Advection : flux form419 % ================================================================420 402 \section{Coriolis and advection: flux form} 421 403 \label{sec:DYN_adv_cor_flux} … … 430 412 At the lateral boundaries either free slip, 431 413 no slip or partial slip boundary conditions are applied following \autoref{chap:LBC}. 432 433 414 434 415 %-------------------------------------------------------------------------------------------------------------- … … 562 543 %%% 563 544 564 % ================================================================565 % Hydrostatic pressure gradient term566 % ================================================================567 545 \section[Hydrostatic pressure gradient (\textit{dynhpg.F90})]{Hydrostatic pressure gradient (\protect\mdl{dynhpg})} 568 546 \label{sec:DYN_hpg} … … 780 758 This option is controlled by \np{nn_dynhpg_rst}{nn\_dynhpg\_rst}, a namelist parameter. 781 759 782 % ================================================================783 % Surface Pressure Gradient784 % ================================================================785 760 \section[Surface pressure gradient (\textit{dynspg.F90})]{Surface pressure gradient (\protect\mdl{dynspg})} 786 761 \label{sec:DYN_spg} … … 809 784 so that the update of the next velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. 810 785 811 812 786 The form of the surface pressure gradient term depends on how the user wants to 813 787 handle the fast external gravity waves that are a solution of the analytical equation (\autoref{sec:MB_hor_pg}). … … 820 794 The extra term introduced in the filtered method is calculated implicitly, so that a solver is used to compute it. 821 795 As a consequence the update of the $next$ velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. 822 823 796 824 797 %-------------------------------------------------------------------------------------------------------------- … … 1092 1065 %>>>>>=============== 1093 1066 1094 1095 1067 %-------------------------------------------------------------------------------------------------------------- 1096 1068 % Filtered free surface formulation … … 1121 1093 It is computed once and for all and applies to all ocean time steps. 1122 1094 1123 % ================================================================1124 % Lateral diffusion term1125 % ================================================================1126 1095 \section[Lateral diffusion term and operators (\textit{dynldf.F90})]{Lateral diffusion term and operators (\protect\mdl{dynldf})} 1127 1096 \label{sec:DYN_ldf} … … 1161 1130 } 1162 1131 1163 % ================================================================1164 \subsection[Iso-level laplacian (\forcode{ln_dynldf_lap})]{Iso-level laplacian operator (\protect\np{ln_dynldf_lap}{ln\_dynldf\_lap})}1165 \label{subsec:DYN_ldf_lap}1166 1167 For lateral iso-level diffusion, the discrete operator is:1168 \begin{equation}1169 \label{eq:DYN_ldf_lap}1170 \left\{1171 \begin{aligned}1172 D_u^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm}1173 \;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta_j \left[1174 {A_f^{lm} \;e_{3f} \zeta } \right] \\ \\1175 D_v^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm}1176 \;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta_i \left[1177 {A_f^{lm} \;e_{3f} \zeta } \right]1178 \end{aligned}1179 \right.1180 \end{equation}1181 1182 As explained in \autoref{subsec:MB_ldf},1183 this formulation (as the gradient of a divergence and curl of the vorticity) preserves symmetry and1184 ensures a complete separation between the vorticity and divergence parts of the momentum diffusion.1185 1186 %--------------------------------------------------------------------------------------------------------------1187 % Rotated laplacian operator1188 %--------------------------------------------------------------------------------------------------------------1189 \subsection[Rotated laplacian (\forcode{ln_dynldf_iso})]{Rotated laplacian operator (\protect\np{ln_dynldf_iso}{ln\_dynldf\_iso})}1190 \label{subsec:DYN_ldf_iso}1191 1192 A rotation of the lateral momentum diffusion operator is needed in several cases:1193 for iso-neutral diffusion in the $z$-coordinate (\np[=.true.]{ln_dynldf_iso}{ln\_dynldf\_iso}) and1194 for either iso-neutral (\np[=.true.]{ln_dynldf_iso}{ln\_dynldf\_iso}) or1195 geopotential (\np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}) diffusion in the $s$-coordinate.1196 In the partial step case, coordinates are horizontal except at the deepest level and1197 no rotation is performed when \np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}.1198 The diffusion operator is defined simply as the divergence of down gradient momentum fluxes on1199 each momentum component.1200 It must be emphasized that this formulation ignores constraints on the stress tensor such as symmetry.1201 The resulting discrete representation is:1202 \begin{equation}1203 \label{eq:DYN_ldf_iso}1204 \begin{split}1205 D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\1206 & \left\{\quad {\delta_{i+1/2} \left[ {A_T^{lm} \left(1207 {\frac{e_{2t} \; e_{3t} }{e_{1t} } \,\delta_{i}[u]1208 -e_{2t} \; r_{1t} \,\overline{\overline {\delta_{k+1/2}[u]}}^{\,i,\,k}}1209 \right)} \right]} \right. \\1210 & \qquad +\ \delta_j \left[ {A_f^{lm} \left( {\frac{e_{1f}\,e_{3f} }{e_{2f}1211 }\,\delta_{j+1/2} [u] - e_{1f}\, r_{2f}1212 \,\overline{\overline {\delta_{k+1/2} [u]}} ^{\,j+1/2,\,k}}1213 \right)} \right] \\1214 &\qquad +\ \delta_k \left[ {A_{uw}^{lm} \left( {-e_{2u} \, r_{1uw} \,\overline{\overline1215 {\delta_{i+1/2} [u]}}^{\,i+1/2,\,k+1/2} }1216 \right.} \right. \\1217 & \ \qquad \qquad \qquad \quad\1218 - e_{1u} \, r_{2uw} \,\overline{\overline {\delta_{j+1/2} [u]}} ^{\,j,\,k+1/2} \\1219 & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\1220 +\frac{e_{1u}\, e_{2u} }{e_{3uw} }\,\left( {r_{1uw}^2+r_{2uw}^2}1221 \right)\,\delta_{k+1/2} [u]} \right)} \right]\;\;\;} \right\} \\ \\1222 D_v^{l\textbf{V}} &= \frac{1}{e_{1v} \, e_{2v} \, e_{3v} } \\1223 & \left\{\quad {\delta_{i+1/2} \left[ {A_f^{lm} \left(1224 {\frac{e_{2f} \; e_{3f} }{e_{1f} } \,\delta_{i+1/2}[v]1225 -e_{2f} \; r_{1f} \,\overline{\overline {\delta_{k+1/2}[v]}}^{\,i+1/2,\,k}}1226 \right)} \right]} \right. \\1227 & \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1t}\,e_{3t} }{e_{2t}1228 }\,\delta_{j} [v] - e_{1t}\, r_{2t}1229 \,\overline{\overline {\delta_{k+1/2} [v]}} ^{\,j,\,k}}1230 \right)} \right] \\1231 & \qquad +\ \delta_k \left[ {A_{vw}^{lm} \left( {-e_{2v} \, r_{1vw} \,\overline{\overline1232 {\delta_{i+1/2} [v]}}^{\,i+1/2,\,k+1/2} }\right.} \right. \\1233 & \ \qquad \qquad \qquad \quad\1234 - e_{1v} \, r_{2vw} \,\overline{\overline {\delta_{j+1/2} [v]}} ^{\,j+1/2,\,k+1/2} \\1235 & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\1236 +\frac{e_{1v}\, e_{2v} }{e_{3vw} }\,\left( {r_{1vw}^2+r_{2vw}^2}1237 \right)\,\delta_{k+1/2} [v]} \right)} \right]\;\;\;} \right\}1238 \end{split}1239 \end{equation}1240 where $r_1$ and $r_2$ are the slopes between the surface along which the diffusion operator acts and1241 the surface of computation ($z$- or $s$-surfaces).1242 The way these slopes are evaluated is given in the lateral physics chapter (\autoref{chap:LDF}).1243 1244 %--------------------------------------------------------------------------------------------------------------1245 % Iso-level bilaplacian operator1246 %--------------------------------------------------------------------------------------------------------------1247 \subsection[Iso-level bilaplacian (\forcode{ln_dynldf_bilap})]{Iso-level bilaplacian operator (\protect\np{ln_dynldf_bilap}{ln\_dynldf\_bilap})}1248 \label{subsec:DYN_ldf_bilap}1249 1250 The lateral fourth order operator formulation on momentum is obtained by applying \autoref{eq:DYN_ldf_lap} twice.1251 It requires an additional assumption on boundary conditions:1252 the first derivative term normal to the coast depends on the free or no-slip lateral boundary conditions chosen,1253 while the third derivative terms normal to the coast are set to zero (see \autoref{chap:LBC}).1254 %%%1255 \gmcomment{add a remark on the the change in the position of the coefficient}1256 %%%1257 1258 % ================================================================1259 1132 % Vertical diffusion term 1260 % ================================================================1261 \section[Vertical diffusion term (\textit{dynzdf.F90})]{Vertical diffusion term (\protect\mdl{dynzdf})}1262 \label{sec:DYN_zdf}1263 %----------------------------------------------namzdf------------------------------------------------------1264 1265 %-------------------------------------------------------------------------------------------------------------1266 1267 Options are defined through the \nam{zdf}{zdf} namelist variables.1268 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.1269 Two time stepping schemes can be used for the vertical diffusion term:1270 $(a)$ a forward time differencing scheme1271 (\np[=.true.]{ln_zdfexp}{ln\_zdfexp}) using a time splitting technique (\np{nn_zdfexp}{nn\_zdfexp} $>$ 1) or1272 $(b)$ a backward (or implicit) time differencing scheme (\np[=.false.]{ln_zdfexp}{ln\_zdfexp})1273 (see \autoref{chap:TD}).1274 Note that namelist variables \np{ln_zdfexp}{ln\_zdfexp} and \np{nn_zdfexp}{nn\_zdfexp} apply to both tracers and dynamics.1275 1276 The formulation of the vertical subgrid scale physics is the same whatever the vertical coordinate is.1277 The vertical diffusion operators given by \autoref{eq:MB_zdf} take the following semi-discrete space form:1278 \[1279 % \label{eq:DYN_zdf}1280 \left\{1281 \begin{aligned}1282 D_u^{vm} &\equiv \frac{1}{e_{3u}} \ \delta_k \left[ \frac{A_{uw}^{vm} }{e_{3uw} }1283 \ \delta_{k+1/2} [\,u\,] \right] \\1284 \\1285 D_v^{vm} &\equiv \frac{1}{e_{3v}} \ \delta_k \left[ \frac{A_{vw}^{vm} }{e_{3vw} }1286 \ \delta_{k+1/2} [\,v\,] \right]1287 \end{aligned}1288 \right.1289 \]1290 where $A_{uw}^{vm} $ and $A_{vw}^{vm} $ are the vertical eddy viscosity and diffusivity coefficients.1291 The way these coefficients are evaluated depends on the vertical physics used (see \autoref{chap:ZDF}).1292 1293 The surface boundary condition on momentum is the stress exerted by the wind.1294 At the surface, the momentum fluxes are prescribed as the boundary condition on1295 the vertical turbulent momentum fluxes,1296 \begin{equation}1297 \label{eq:DYN_zdf_sbc}1298 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1}1299 = \frac{1}{\rho_o} \binom{\tau_u}{\tau_v }1300 \end{equation}1301 where $\left( \tau_u ,\tau_v \right)$ are the two components of the wind stress vector in1302 the (\textbf{i},\textbf{j}) coordinate system.1303 The high mixing coefficients in the surface mixed layer ensure that the surface wind stress is distributed in1304 the vertical over the mixed layer depth.1305 If the vertical mixing coefficient is small (when no mixed layer scheme is used)1306 the surface stress enters only the top model level, as a body force.1307 The surface wind stress is calculated in the surface module routines (SBC, see \autoref{chap:SBC}).1308 1309 The turbulent flux of momentum at the bottom of the ocean is specified through a bottom friction parameterisation1310 (see \autoref{sec:ZDF_drg})1311 1312 % ================================================================1313 1133 % External Forcing 1314 % ================================================================1315 \section{External forcings}1316 \label{sec:DYN_forcing}1317 1318 Besides the surface and bottom stresses (see the above section)1319 which are introduced as boundary conditions on the vertical mixing,1320 three other forcings may enter the dynamical equations by affecting the surface pressure gradient.1321 1322 (1) When \np[=.true.]{ln_apr_dyn}{ln\_apr\_dyn} (see \autoref{sec:SBC_apr}),1323 the atmospheric pressure is taken into account when computing the surface pressure gradient.1324 1325 (2) When \np[=.true.]{ln_tide_pot}{ln\_tide\_pot} and \np[=.true.]{ln_tide}{ln\_tide} (see \autoref{sec:SBC_tide}),1326 the tidal potential is taken into account when computing the surface pressure gradient.1327 1328 (3) When \np[=2]{nn_ice_embd}{nn\_ice\_embd} and LIM or CICE is used1329 (\ie\ when the sea-ice is embedded in the ocean),1330 the snow-ice mass is taken into account when computing the surface pressure gradient.1331 1332 1333 \gmcomment{ missing : the lateral boundary condition !!! another external forcing1334 }1335 1336 % ================================================================1337 1134 % Wetting and drying 1338 % ================================================================1339 \section{Wetting and drying }1340 \label{sec:DYN_wetdry}1341 1342 There are two main options for wetting and drying code (wd):1343 (a) an iterative limiter (il) and (b) a directional limiter (dl).1344 The directional limiter is based on the scheme developed by \cite{warner.defne.ea_CG13} for RO1345 MS1346 which was in turn based on ideas developed for POM by \cite{oey_OM06}. The iterative1347 limiter is a new scheme. The iterative limiter is activated by setting $\mathrm{ln\_wd\_il} = \mathrm{.true.}$1348 and $\mathrm{ln\_wd\_dl} = \mathrm{.false.}$. The directional limiter is activated1349 by setting $\mathrm{ln\_wd\_dl} = \mathrm{.true.}$ and $\mathrm{ln\_wd\_il} = \mathrm{.false.}$.1350 1351 \begin{listing}1352 \nlst{namwad}1353 \caption{\forcode{&namwad}}1354 \label{lst:namwad}1355 \end{listing}1356 1357 The following terminology is used. The depth of the topography (positive downwards)1358 at each $(i,j)$ point is the quantity stored in array $\mathrm{ht\_wd}$ in the \NEMO\ code.1359 The height of the free surface (positive upwards) is denoted by $ \mathrm{ssh}$. Given the sign1360 conventions used, the water depth, $h$, is the height of the free surface plus the depth of the1361 topography (i.e. $\mathrm{ssh} + \mathrm{ht\_wd}$).1362 1363 Both wd schemes take all points in the domain below a land elevation of $\mathrm{rn\_wdld}$ to be1364 covered by water. They require the topography specified with a model1365 configuration to have negative depths at points where the land is higher than the1366 topography's reference sea-level. The vertical grid in \NEMO\ is normally computed relative to an1367 initial state with zero sea surface height elevation.1368 The user can choose to compute the vertical grid and heights in the model relative to1369 a non-zero reference height for the free surface. This choice affects the calculation of the metrics and depths1370 (i.e. the $\mathrm{e3t\_0, ht\_0}$ etc. arrays).1371 1372 Points where the water depth is less than $\mathrm{rn\_wdmin1}$ are interpreted as ``dry''.1373 $\mathrm{rn\_wdmin1}$ is usually chosen to be of order $0.05$m but extreme topographies1374 with very steep slopes require larger values for normal choices of time-step. Surface fluxes1375 are also switched off for dry cells to prevent freezing, boiling etc. of very thin water layers.1376 The fluxes are tappered down using a $\mathrm{tanh}$ weighting function1377 to no flux as the dry limit $\mathrm{rn\_wdmin1}$ is approached. Even wet cells can be very shallow.1378 The depth at which to start tapering is controlled by the user by setting $\mathrm{rn\_wd\_sbcdep}$.1379 The fraction $(<1)$ of sufrace fluxes to use at this depth is set by $\mathrm{rn\_wd\_sbcfra}$.1380 1381 Both versions of the code have been tested in six test cases provided in the WAD\_TEST\_CASES configuration1382 and in ``realistic'' configurations covering parts of the north-west European shelf.1383 All these configurations have used pure sigma coordinates. It is expected that1384 the wetting and drying code will work in domains with more general s-coordinates provided1385 the coordinates are pure sigma in the region where wetting and drying actually occurs.1386 1387 The next sub-section descrbies the directional limiter and the following sub-section the iterative limiter.1388 The final sub-section covers some additional considerations that are relevant to both schemes.1389 1390 1391 %-----------------------------------------------------------------------------------------1392 % Iterative limiters1393 %-----------------------------------------------------------------------------------------1394 \subsection[Directional limiter (\textit{wet\_dry.F90})]{Directional limiter (\mdl{wet\_dry})}1395 \label{subsec:DYN_wd_directional_limiter}1396 1397 The principal idea of the directional limiter is that1398 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}).1399 1400 All the changes associated with this option are made to the barotropic solver for the non-linear1401 free surface code within dynspg\_ts.1402 On each barotropic sub-step the scheme determines the direction of the flow across each face of all the tracer cells1403 and sets the flux across the face to zero when the flux is from a dry tracer cell. This prevents cells1404 whose depth is rn\_wdmin1 or less from drying out further. The scheme does not force $h$ (the water depth) at tracer cells1405 to be at least the minimum depth and hence is able to conserve mass / volume.1406 1407 The flux across each $u$-face of a tracer cell is multiplied by a factor zuwdmask (an array which depends on ji and jj).1408 If the user sets \np[=.false.]{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp} then zuwdmask is 1 when the1409 flux is from a cell with water depth greater than \np{rn_wdmin1}{rn\_wdmin1} and 0 otherwise. If the user sets1410 \np[=.true.]{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp} the flux across the face is ramped down as the water depth decreases1411 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.1412 1413 At the point where the flux across a $u$-face is multiplied by zuwdmask , we have chosen1414 also to multiply the corresponding velocity on the ``now'' step at that face by zuwdmask. We could have1415 chosen not to do that and to allow fairly large velocities to occur in these ``dry'' cells.1416 The rationale for setting the velocity to zero is that it is the momentum equations that are being solved1417 and the total momentum of the upstream cell (treating it as a finite volume) should be considered1418 to be its depth times its velocity. This depth is considered to be zero at ``dry'' $u$-points consistent with its1419 treatment in the calculation of the flux of mass across the cell face.1420 1421 1422 \cite{warner.defne.ea_CG13} state that in their scheme the velocity masks at the cell faces for the baroclinic1423 timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic sub-steps is respectively less than1424 or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer1425 fields (tracers independent of $x$, $y$ and $z$). Our scheme conserves constant tracers because1426 the velocities used at the tracer cell faces on the baroclinic timesteps are carefully calculated by dynspg\_ts1427 to equal their mean value during the barotropic steps. If the user sets \np[=.true.]{ln_wd_dl_bc}{ln\_wd\_dl\_bc}, the1428 baroclinic velocities are also multiplied by a suitably weighted average of zuwdmask.1429 1430 %-----------------------------------------------------------------------------------------1431 % Iterative limiters1432 %-----------------------------------------------------------------------------------------1433 1434 \subsection[Iterative limiter (\textit{wet\_dry.F90})]{Iterative limiter (\mdl{wet\_dry})}1435 \label{subsec:DYN_wd_iterative_limiter}1436 1437 \subsubsection[Iterative flux limiter (\textit{wet\_dry.F90})]{Iterative flux limiter (\mdl{wet\_dry})}1438 \label{subsec:DYN_wd_il_spg_limiter}1439 1440 The iterative limiter modifies the fluxes across the faces of cells that are either already ``dry''1441 or may become dry within the next time-step using an iterative method.1442 1443 The flux limiter for the barotropic flow (devised by Hedong Liu) can be understood as follows:1444 1445 The continuity equation for the total water depth in a column1446 \begin{equation}1447 \label{eq:DYN_wd_continuity}1448 \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 .1449 \end{equation}1450 can be written in discrete form as1451 1452 \begin{align}1453 \label{eq:DYN_wd_continuity_2}1454 \frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) )1455 &= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j} + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\1456 &= \mathrm{zzflx}_{i,j} .1457 \end{align}1458 1459 In the above $h$ is the depth of the water in the column at point $(i,j)$,1460 $\mathrm{flxu}_{i+1,j}$ is the flux out of the ``eastern'' face of the cell and1461 $\mathrm{flxv}_{i,j+1}$ the flux out of the ``northern'' face of the cell; $t_{n+1}$ is1462 the new timestep, $t_e$ is the old timestep (either $t_b$ or $t_n$) and $ \Delta t =1463 t_{n+1} - t_e$; $e_1 e_2$ is the area of the tracer cells centred at $(i,j)$ and1464 $\mathrm{zzflx}$ is the sum of the fluxes through all the faces.1465 1466 The flux limiter splits the flux $\mathrm{zzflx}$ into fluxes that are out of the cell1467 (zzflxp) and fluxes that are into the cell (zzflxn). Clearly1468 1469 \begin{equation}1470 \label{eq:DYN_wd_zzflx_p_n_1}1471 \mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} .1472 \end{equation}1473 1474 The flux limiter iteratively adjusts the fluxes $\mathrm{flxu}$ and $\mathrm{flxv}$ until1475 none of the cells will ``dry out''. To be precise the fluxes are limited until none of the1476 cells has water depth less than $\mathrm{rn\_wdmin1}$ on step $n+1$.1477 1478 Let the fluxes on the $m$th iteration step be denoted by $\mathrm{flxu}^{(m)}$ and1479 $\mathrm{flxv}^{(m)}$. Then the adjustment is achieved by seeking a set of coefficients,1480 $\mathrm{zcoef}_{i,j}^{(m)}$ such that:1481 1482 \begin{equation}1483 \label{eq:DYN_wd_continuity_coef}1484 \begin{split}1485 \mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\1486 \mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j}1487 \end{split}1488 \end{equation}1489 1490 where the coefficients are $1.0$ generally but can vary between $0.0$ and $1.0$ around1491 cells that would otherwise dry.1492 1493 The iteration is initialised by setting1494 1495 \begin{equation}1496 \label{eq:DYN_wd_zzflx_initial}1497 \mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} .1498 \end{equation}1499 1500 The fluxes out of cell $(i,j)$ are updated at the $m+1$th iteration if the depth of the1501 cell on timestep $t_e$, namely $h_{i,j}(t_e)$, is less than the total flux out of the cell1502 times the timestep divided by the cell area. Using (\autoref{eq:DYN_wd_continuity_2}) this1503 condition is1504 1505 \begin{equation}1506 \label{eq:DYN_wd_continuity_if}1507 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} ) .1508 \end{equation}1509 1510 Rearranging (\autoref{eq:DYN_wd_continuity_if}) we can obtain an expression for the maximum1511 outward flux that can be allowed and still maintain the minimum wet depth:1512 1513 \begin{equation}1514 \label{eq:DYN_wd_max_flux}1515 \begin{split}1516 \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{]} \\1517 \phantom{[} & - \mathrm{zzflxn}^{(m)}_{i,j} \Big]1518 \end{split}1519 \end{equation}1520 1521 Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\itshape [Q: Why is1522 this necessary/desirable?]}. Substituting from (\autoref{eq:DYN_wd_continuity_coef}) gives an1523 expression for the coefficient needed to multiply the outward flux at this cell in order1524 to avoid drying.1525 1526 \begin{equation}1527 \label{eq:DYN_wd_continuity_nxtcoef}1528 \begin{split}1529 \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{]} \\1530 \phantom{[} & - \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} }1531 \end{split}1532 \end{equation}1533 1534 Only the outward flux components are altered but, of course, outward fluxes from one cell1535 are inward fluxes to adjacent cells and the balance in these cells may need subsequent1536 adjustment; hence the iterative nature of this scheme. Note, for example, that the flux1537 across the ``eastern'' face of the $(i,j)$th cell is only updated at the $m+1$th iteration1538 if that flux at the $m$th iteration is out of the $(i,j)$th cell. If that is the case then1539 the flux across that face is into the $(i+1,j)$ cell and that flux will not be updated by1540 the calculation for the $(i+1,j)$th cell. In this sense the updates to the fluxes across1541 the faces of the cells do not ``compete'' (they do not over-write each other) and one1542 would expect the scheme to converge relatively quickly. The scheme is flux based so1543 conserves mass. It also conserves constant tracers for the same reason that the1544 directional limiter does.1545 1546 1547 %----------------------------------------------------------------------------------------1548 % Surface pressure gradients1549 %----------------------------------------------------------------------------------------1550 \subsubsection[Modification of surface pressure gradients (\textit{dynhpg.F90})]{Modification of surface pressure gradients (\mdl{dynhpg})}1551 \label{subsec:DYN_wd_il_spg}1552 1553 At ``dry'' points the water depth is usually close to $\mathrm{rn\_wdmin1}$. If the1554 topography is sloping at these points the sea-surface will have a similar slope and there1555 will hence be very large horizontal pressure gradients at these points. The WAD modifies1556 the magnitude but not the sign of the surface pressure gradients (zhpi and zhpj) at such1557 points by mulitplying them by positive factors (zcpx and zcpy respectively) that lie1558 between $0$ and $1$.1559 1560 We describe how the scheme works for the ``eastward'' pressure gradient, zhpi, calculated1561 at the $(i,j)$th $u$-point. The scheme uses the ht\_wd depths and surface heights at the1562 neighbouring $(i+1,j)$ and $(i,j)$ tracer points. zcpx is calculated using two logicals1563 variables, $\mathrm{ll\_tmp1}$ and $\mathrm{ll\_tmp2}$ which are evaluated for each grid1564 column. The three possible combinations are illustrated in \autoref{fig:DYN_WAD_dynhpg}.1565 1566 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>1567 \begin{figure}[!ht]1568 \centering1569 \includegraphics[width=0.66\textwidth]{Fig_WAD_dynhpg}1570 \caption[Combinations controlling the limiting of the horizontal pressure gradient in1571 wetting and drying regimes]{1572 Three possible combinations of the logical variables controlling the1573 limiting of the horizontal pressure gradient in wetting and drying regimes}1574 \label{fig:DYN_WAD_dynhpg}1575 \end{figure}1576 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>1577 1578 The first logical, $\mathrm{ll\_tmp1}$, is set to true if and only if the water depth at1579 both neighbouring points is greater than $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ and1580 the minimum height of the sea surface at the two points is greater than the maximum height1581 of the topography at the two points:1582 1583 \begin{equation}1584 \label{eq:DYN_ll_tmp1}1585 \begin{split}1586 \mathrm{ll\_tmp1} = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\1587 & \quad \mathrm{MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj))\ .and.} \\1588 & \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\1589 & \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\1590 & \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 }1591 \end{split}1592 \end{equation}1593 1594 The second logical, $\mathrm{ll\_tmp2}$, is set to true if and only if the maximum height1595 of the sea surface at the two points is greater than the maximum height of the topography1596 at the two points plus $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$1597 1598 \begin{equation}1599 \label{eq:DYN_ll_tmp2}1600 \begin{split}1601 \mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\1602 & \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\1603 & \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) .1604 \end{split}1605 \end{equation}1606 1607 If $\mathrm{ll\_tmp1}$ is true then the surface pressure gradient, zhpi at the $(i,j)$1608 point is unmodified. If both logicals are false zhpi is set to zero.1609 1610 If $\mathrm{ll\_tmp1}$ is true and $\mathrm{ll\_tmp2}$ is false then the surface pressure1611 gradient is multiplied through by zcpx which is the absolute value of the difference in1612 the water depths at the two points divided by the difference in the surface heights at the1613 two points. Thus the sign of the sea surface height gradient is retained but the magnitude1614 of the pressure force is determined by the difference in water depths rather than the1615 difference in surface height between the two points. Note that dividing by the difference1616 between the sea surface heights can be problematic if the heights approach parity. An1617 additional condition is applied to $\mathrm{ ll\_tmp2 }$ to ensure it is .false. in such1618 conditions.1619 1620 \subsubsection[Additional considerations (\textit{usrdef\_zgr.F90})]{Additional considerations (\mdl{usrdef\_zgr})}1621 \label{subsec:DYN_WAD_additional}1622 1623 In the very shallow water where wetting and drying occurs the parametrisation of1624 bottom drag is clearly very important. In order to promote stability1625 it is sometimes useful to calculate the bottom drag using an implicit time-stepping approach.1626 1627 Suitable specifcation of the surface heat flux in wetting and drying domains in forced and1628 coupled simulations needs further consideration. In order to prevent freezing or boiling1629 in uncoupled integrations the net surface heat fluxes need to be appropriately limited.1630 1631 %----------------------------------------------------------------------------------------1632 % The WAD test cases1633 %----------------------------------------------------------------------------------------1634 \subsection[The WAD test cases (\textit{usrdef\_zgr.F90})]{The WAD test cases (\mdl{usrdef\_zgr})}1635 \label{subsec:DYN_WAD_test_cases}1636 1637 See the WAD tests MY\_DOC documention for details of the WAD test cases.1638 1639 1640 1641 % ================================================================1642 1135 % Time evolution term 1643 % ================================================================1644 \section[Time evolution term (\textit{dynnxt.F90})]{Time evolution term (\protect\mdl{dynnxt})}1645 \label{sec:DYN_nxt}1646 1647 %----------------------------------------------namdom----------------------------------------------------1648 1649 %-------------------------------------------------------------------------------------------------------------1650 1651 Options are defined through the \nam{dom}{dom} namelist variables.1652 The general framework for dynamics time stepping is a leap-frog scheme,1653 \ie\ a three level centred time scheme associated with an Asselin time filter (cf. \autoref{chap:TD}).1654 The scheme is applied to the velocity, except when1655 using the flux form of momentum advection (cf. \autoref{sec:DYN_adv_cor_flux})1656 in the variable volume case (\texttt{vvl?} defined),1657 where it has to be applied to the thickness weighted velocity (see \autoref{sec:SCOORD_momentum})1658 1659 $\bullet$ vector invariant form or linear free surface1660 (\np[=.true.]{ln_dynhpg_vec}{ln\_dynhpg\_vec} ; \texttt{vvl?} not defined):1661 \[1662 % \label{eq:DYN_nxt_vec}1663 \left\{1664 \begin{aligned}1665 &u^{t+\rdt} = u_f^{t-\rdt} + 2\rdt \ \text{RHS}_u^t \\1666 &u_f^t \;\quad = u^t+\gamma \,\left[ {u_f^{t-\rdt} -2u^t+u^{t+\rdt}} \right]1667 \end{aligned}1668 \right.1669 \]1670 1671 $\bullet$ flux form and nonlinear free surface1672 (\np[=.false.]{ln_dynhpg_vec}{ln\_dynhpg\_vec} ; \texttt{vvl?} defined):1673 \[1674 % \label{eq:DYN_nxt_flux}1675 \left\{1676 \begin{aligned}1677 &\left(e_{3u}\,u\right)^{t+\rdt} = \left(e_{3u}\,u\right)_f^{t-\rdt} + 2\rdt \; e_{3u} \;\text{RHS}_u^t \\1678 &\left(e_{3u}\,u\right)_f^t \;\quad = \left(e_{3u}\,u\right)^t1679 +\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]1680 \end{aligned}1681 \right.1682 \]1683 where RHS is the right hand side of the momentum equation,1684 the subscript $f$ denotes filtered values and $\gamma$ is the Asselin coefficient.1685 $\gamma$ is initialized as \np{nn_atfp}{nn\_atfp} (namelist parameter).1686 Its default value is \np[=10.e-3]{nn_atfp}{nn\_atfp}.1687 In both cases, the modified Asselin filter is not applied since perfect conservation is not an issue for1688 the momentum equations.1689 1690 Note that with the filtered free surface,1691 the update of the \textit{after} velocities is done in the \mdl{dynsp\_flt} module,1692 and only array swapping and Asselin filtering is done in \mdl{dynnxt}.1693 1694 \onlyinsubfile{\input{../../global/epilogue}}1695 1696 \end{document}
Note: See TracChangeset
for help on using the changeset viewer.
