Changeset 10499 for NEMO/trunk/doc/latex
- Timestamp:
- 2019-01-10T16:12:24+01:00 (5 years ago)
- Location:
- NEMO/trunk/doc/latex/NEMO
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/doc/latex/NEMO/main/NEMO_manual.bib
r10495 r10499 2946 2946 } 2947 2947 2948 @ARTICLE{Oey06, 2949 author = {Lie-Yauw Oey}, 2950 title = {An OGCM with movable land-sea boundaries}, 2951 journal = {Ocean Modelling}, 2952 volume = "13", 2953 number = "2", 2954 pages = "176 - 195", 2955 year = "2006", 2956 issn = "1463-5003", 2957 doi = "https://doi.org/10.1016/j.ocemod.2006.01.001", 2958 url = "http://www.sciencedirect.com/science/article/pii/S1463500306000084", 2959 } 2960 2948 2961 @PHDTHESIS{Olivier_PhD01, 2949 2962 author = {F. Olivier}, … … 3755 3768 } 3756 3769 3770 @ARTICLE{WarnerEtal13, 3771 author = {John C. Warner and Zafer Defne and Kevin Haas and 3772 Hernan G. Arango}, 3773 title = {A wetting and drying scheme for ROMS}, 3774 journal = "Computers \& Geosciences", 3775 volume = "58", 3776 pages = "54 - 61", 3777 year = "2013", 3778 issn = "0098-3004", 3779 doi = "https://doi.org/10.1016/j.cageo.2013.05.004", 3780 url = "http://www.sciencedirect.com/science/article/pii/S0098300413001362", 3781 } 3782 3757 3783 @ARTICLE{Weatherly_JMR84, 3758 3784 author = {G. L. Weatherly}, -
NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex
r10468 r10499 1320 1320 1321 1321 % ================================================================ 1322 % Wetting and drying 1323 % ================================================================ 1324 \section{Wetting and drying } 1325 \label{sec:DYN_wetdry} 1326 There are two main options for wetting and drying code (wd): 1327 (a) an iterative limiter (il) and (b) a directional limiter (dl). 1328 The directional limiter is based on the scheme developed by \cite{WarnerEtal13} for RO 1329 MS 1330 which was in turn based on ideas developed for POM by \cite{Oey06}. The iterative 1331 limiter is a new scheme. The iterative limiter is activated by setting $\mathrm{ln\_wd\_il} = \mathrm{.true.}$ 1332 and $\mathrm{ln\_wd\_dl} = \mathrm{.false.}$. The directional limiter is activated 1333 by setting $\mathrm{ln\_wd\_dl} = \mathrm{.true.}$ and $\mathrm{ln\_wd\_il} = \mathrm{.false.}$. 1334 1335 \nlst{namwad} 1336 1337 The following terminology is used. The depth of the topography (positive downwards) 1338 at each $(i,j)$ point is the quantity stored in array $\mathrm{ht\_wd}$ in the NEMO code. 1339 The height of the free surface (positive upwards) is denoted by $ \mathrm{ssh}$. Given the sign 1340 conventions used, the water depth, $h$, is the height of the free surface plus the depth of the 1341 topography (i.e. $\mathrm{ssh} + \mathrm{ht\_wd}$). 1342 1343 Both wd schemes take all points in the domain below a land elevation of $\mathrm{rn\_wdld}$ to be 1344 covered by water. They require the topography specified with a model 1345 configuration to have negative depths at points where the land is higher than the 1346 topography's reference sea-level. The vertical grid in NEMO is normally computed relative to an 1347 initial state with zero sea surface height elevation. 1348 The user can choose to compute the vertical grid and heights in the model relative to 1349 a non-zero reference height for the free surface. This choice affects the calculation of the metrics and depths 1350 (i.e. the $\mathrm{e3t\_0, ht\_0}$ etc. arrays). 1351 1352 Points where the water depth is less than $\mathrm{rn\_wdmin1}$ are interpreted as ``dry''. 1353 $\mathrm{rn\_wdmin1}$ is usually chosen to be of order $0.05$m but extreme topographies 1354 with very steep slopes require larger values for normal choices of time-step. Surface fluxes 1355 are also switched off for dry cells to prevent freezing, boiling etc. of very thin water layers. 1356 The fluxes are tappered down using a $\mathrm{tanh}$ weighting function 1357 to no flux as the dry limit $\mathrm{rn\_wdmin1}$ is approached. Even wet cells can be very shallow. 1358 The depth at which to start tapering is controlled by the user by setting $\mathrm{rn\_wd\_sbcdep}$. 1359 The fraction $(<1)$ of sufrace fluxes to use at this depth is set by $\mathrm{rn\_wd\_sbcfra}$. 1360 1361 Both versions of the code have been tested in six test cases provided in the WAD\_TEST\_CASES configuration 1362 and in ``realistic'' configurations covering parts of the north-west European shelf. 1363 All these configurations have used pure sigma coordinates. It is expected that 1364 the wetting and drying code will work in domains with more general s-coordinates provided 1365 the coordinates are pure sigma in the region where wetting and drying actually occurs. 1366 1367 The next sub-section descrbies the directional limiter and the following sub-section the iterative limiter. 1368 The final sub-section covers some additional considerations that are relevant to both schemes. 1369 1370 1371 %----------------------------------------------------------------------------------------- 1372 % Iterative limiters 1373 %----------------------------------------------------------------------------------------- 1374 \subsection [Directional limiter (\textit{wet\_dry})] 1375 {Directional limiter (\mdl{wet\_dry})} 1376 \label{subsec:DYN_wd_directional_limiter} 1377 The principal idea of the directional limiter is that 1378 water should not be allowed to flow out of a dry tracer cell (i.e. one whose water depth is less than rn\_wdmin1). 1379 1380 All the changes associated with this option are made to the barotropic solver for the non-linear 1381 free surface code within dynspg\_ts. 1382 On each barotropic sub-step the scheme determines the direction of the flow across each face of all the tracer cells 1383 and sets the flux across the face to zero when the flux is from a dry tracer cell. This prevents cells 1384 whose depth is rn\_wdmin1 or less from drying out further. The scheme does not force $h$ (the water depth) at tracer cells 1385 to be at least the minimum depth and hence is able to conserve mass / volume. 1386 1387 The flux across each $u$-face of a tracer cell is multiplied by a factor zuwdmask (an array which depends on ji and jj). 1388 If the user sets ln\_wd\_dl\_ramp = .False. then zuwdmask is 1 when the 1389 flux is from a cell with water depth greater than rn\_wdmin1 and 0 otherwise. If the user sets 1390 ln\_wd\_dl\_ramp = .True. the flux across the face is ramped down as the water depth decreases 1391 from 2 * rn\_wdmin1 to rn\_wdmin1. The use of this ramp reduced grid-scale noise in idealised test cases. 1392 1393 At the point where the flux across a $u$-face is multiplied by zuwdmask , we have chosen 1394 also to multiply the corresponding velocity on the ``now'' step at that face by zuwdmask. We could have 1395 chosen not to do that and to allow fairly large velocities to occur in these ``dry'' cells. 1396 The rationale for setting the velocity to zero is that it is the momentum equations that are being solved 1397 and the total momentum of the upstream cell (treating it as a finite volume) should be considered 1398 to be its depth times its velocity. This depth is considered to be zero at ``dry'' $u$-points consistent with its 1399 treatment in the calculation of the flux of mass across the cell face. 1400 1401 1402 \cite{WarnerEtal13} state that in their scheme the velocity masks at the cell faces for the baroclinic 1403 timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic sub-steps is respectively less than 1404 or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer 1405 fields (tracers independent of $x$, $y$ and $z$). Our scheme conserves constant tracers because 1406 the velocities used at the tracer cell faces on the baroclinic timesteps are carefully calculated by dynspg\_ts 1407 to equal their mean value during the barotropic steps. If the user sets ln\_wd\_dl\_bc = .True., the 1408 baroclinic velocities are also multiplied by a suitably weighted average of zuwdmask. 1409 1410 %----------------------------------------------------------------------------------------- 1411 % Iterative limiters 1412 %----------------------------------------------------------------------------------------- 1413 1414 \subsection [Iterative limiter (\textit{wet\_dry})] 1415 {Iterative limiter (\mdl{wet\_dry})} 1416 \label{subsec:DYN_wd_iterative_limiter} 1417 1418 \subsubsection [Iterative flux limiter (\textit{wet\_dry})] 1419 {Iterative flux limiter (\mdl{wet\_dry})} 1420 \label{subsubsec:DYN_wd_il_spg_limiter} 1421 1422 The iterative limiter modifies the fluxes across the faces of cells that are either already ``dry'' 1423 or may become dry within the next time-step using an iterative method. 1424 1425 The flux limiter for the barotropic flow (devised by Hedong Liu) can be understood as follows: 1426 1427 The continuity equation for the total water depth in a column 1428 \begin{equation} \label{dyn_wd_continuity} 1429 \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 . 1430 \end{equation} 1431 can be written in discrete form as 1432 1433 \begin{align} \label{dyn_wd_continuity_2} 1434 \frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) ) 1435 &= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j} + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\ 1436 &= \mathrm{zzflx}_{i,j} . 1437 \end{align} 1438 1439 In the above $h$ is the depth of the water in the column at point $(i,j)$, 1440 $\mathrm{flxu}_{i+1,j}$ is the flux out of the ``eastern'' face of the cell and 1441 $\mathrm{flxv}_{i,j+1}$ the flux out of the ``northern'' face of the cell; $t_{n+1}$ is 1442 the new timestep, $t_e$ is the old timestep (either $t_b$ or $t_n$) and $ \Delta t = 1443 t_{n+1} - t_e$; $e_1 e_2$ is the area of the tracer cells centred at $(i,j)$ and 1444 $\mathrm{zzflx}$ is the sum of the fluxes through all the faces. 1445 1446 The flux limiter splits the flux $\mathrm{zzflx}$ into fluxes that are out of the cell 1447 (zzflxp) and fluxes that are into the cell (zzflxn). Clearly 1448 1449 \begin{equation} \label{dyn_wd_zzflx_p_n_1} 1450 \mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} . 1451 \end{equation} 1452 1453 The flux limiter iteratively adjusts the fluxes $\mathrm{flxu}$ and $\mathrm{flxv}$ until 1454 none of the cells will ``dry out''. To be precise the fluxes are limited until none of the 1455 cells has water depth less than $\mathrm{rn\_wdmin1}$ on step $n+1$. 1456 1457 Let the fluxes on the $m$th iteration step be denoted by $\mathrm{flxu}^{(m)}$ and 1458 $\mathrm{flxv}^{(m)}$. Then the adjustment is achieved by seeking a set of coefficients, 1459 $\mathrm{zcoef}_{i,j}^{(m)}$ such that: 1460 1461 \begin{equation} \label{dyn_wd_continuity_coef} 1462 \begin{split} 1463 \mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\ 1464 \mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j} 1465 \end{split} 1466 \end{equation} 1467 1468 where the coefficients are $1.0$ generally but can vary between $0.0$ and $1.0$ around 1469 cells that would otherwise dry. 1470 1471 The iteration is initialised by setting 1472 1473 \begin{equation} \label{dyn_wd_zzflx_initial} 1474 \mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} . 1475 \end{equation} 1476 1477 The fluxes out of cell $(i,j)$ are updated at the $m+1$th iteration if the depth of the 1478 cell on timestep $t_e$, namely $h_{i,j}(t_e)$, is less than the total flux out of the cell 1479 times the timestep divided by the cell area. Using (\ref{dyn_wd_continuity_2}) this 1480 condition is 1481 1482 \begin{equation} \label{dyn_wd_continuity_if} 1483 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} ) . 1484 \end{equation} 1485 1486 Rearranging (\ref{dyn_wd_continuity_if}) we can obtain an expression for the maximum 1487 outward flux that can be allowed and still maintain the minimum wet depth: 1488 1489 \begin{equation} \label{dyn_wd_max_flux} 1490 \begin{split} 1491 \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{]} \\ 1492 \phantom{[} & - \mathrm{zzflxn}^{(m)}_{i,j} \Big] 1493 \end{split} 1494 \end{equation} 1495 1496 Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\it [Q: Why is 1497 this necessary/desirable?]}. Substituting from (\ref{dyn_wd_continuity_coef}) gives an 1498 expression for the coefficient needed to multiply the outward flux at this cell in order 1499 to avoid drying. 1500 1501 \begin{equation} \label{dyn_wd_continuity_nxtcoef} 1502 \begin{split} 1503 \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{]} \\ 1504 \phantom{[} & - \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} } 1505 \end{split} 1506 \end{equation} 1507 1508 Only the outward flux components are altered but, of course, outward fluxes from one cell 1509 are inward fluxes to adjacent cells and the balance in these cells may need subsequent 1510 adjustment; hence the iterative nature of this scheme. Note, for example, that the flux 1511 across the ``eastern'' face of the $(i,j)$th cell is only updated at the $m+1$th iteration 1512 if that flux at the $m$th iteration is out of the $(i,j)$th cell. If that is the case then 1513 the flux across that face is into the $(i+1,j)$ cell and that flux will not be updated by 1514 the calculation for the $(i+1,j)$th cell. In this sense the updates to the fluxes across 1515 the faces of the cells do not ``compete'' (they do not over-write each other) and one 1516 would expect the scheme to converge relatively quickly. The scheme is flux based so 1517 conserves mass. It also conserves constant tracers for the same reason that the 1518 directional limiter does. 1519 1520 1521 %---------------------------------------------------------------------------------------- 1522 % Surface pressure gradients 1523 %---------------------------------------------------------------------------------------- 1524 \subsubsection [Modification of surface pressure gradients (\textit{dynhpg})] 1525 {Modification of surface pressure gradients (\mdl{dynhpg})} 1526 \label{subsubsec:DYN_wd_il_spg} 1527 1528 At ``dry'' points the water depth is usually close to $\mathrm{rn\_wdmin1}$. If the 1529 topography is sloping at these points the sea-surface will have a similar slope and there 1530 will hence be very large horizontal pressure gradients at these points. The WAD modifies 1531 the magnitude but not the sign of the surface pressure gradients (zhpi and zhpj) at such 1532 points by mulitplying them by positive factors (zcpx and zcpy respectively) that lie 1533 between $0$ and $1$. 1534 1535 We describe how the scheme works for the ``eastward'' pressure gradient, zhpi, calculated 1536 at the $(i,j)$th $u$-point. The scheme uses the ht\_wd depths and surface heights at the 1537 neighbouring $(i+1,j)$ and $(i,j)$ tracer points. zcpx is calculated using two logicals 1538 variables, $\mathrm{ll\_tmp1}$ and $\mathrm{ll\_tmp2}$ which are evaluated for each grid 1539 column. The three possible combinations are illustrated in figure \ref{Fig_WAD_dynhpg}. 1540 1541 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1542 \begin{figure}[!ht] \begin{center} 1543 \includegraphics[width=0.8\textwidth]{Fig_WAD_dynhpg} 1544 \caption{ \label{Fig_WAD_dynhpg} 1545 Illustrations of the three possible combinations of the logical variables controlling the 1546 limiting of the horizontal pressure gradient in wetting and drying regimes} 1547 \end{center}\end{figure} 1548 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1549 1550 The first logical, $\mathrm{ll\_tmp1}$, is set to true if and only if the water depth at 1551 both neighbouring points is greater than $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ and 1552 the minimum height of the sea surface at the two points is greater than the maximum height 1553 of the topography at the two points: 1554 1555 \begin{equation} \label{dyn_ll_tmp1} 1556 \begin{split} 1557 \mathrm{ll\_tmp1} = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\ 1558 & \quad \mathrm{MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj))\ .and.} \\ 1559 & \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\ 1560 & \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\ 1561 & \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 } 1562 \end{split} 1563 \end{equation} 1564 1565 The second logical, $\mathrm{ll\_tmp2}$, is set to true if and only if the maximum height 1566 of the sea surface at the two points is greater than the maximum height of the topography 1567 at the two points plus $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ 1568 1569 \begin{equation} \label{dyn_ll_tmp2} 1570 \begin{split} 1571 \mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\ 1572 & \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\ 1573 & \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) . 1574 \end{split} 1575 \end{equation} 1576 1577 If $\mathrm{ll\_tmp1}$ is true then the surface pressure gradient, zhpi at the $(i,j)$ 1578 point is unmodified. If both logicals are false zhpi is set to zero. 1579 1580 If $\mathrm{ll\_tmp1}$ is true and $\mathrm{ll\_tmp2}$ is false then the surface pressure 1581 gradient is multiplied through by zcpx which is the absolute value of the difference in 1582 the water depths at the two points divided by the difference in the surface heights at the 1583 two points. Thus the sign of the sea surface height gradient is retained but the magnitude 1584 of the pressure force is determined by the difference in water depths rather than the 1585 difference in surface height between the two points. Note that dividing by the difference 1586 between the sea surface heights can be problematic if the heights approach parity. An 1587 additional condition is applied to $\mathrm{ ll\_tmp2 }$ to ensure it is .false. in such 1588 conditions. 1589 1590 \subsubsection [Additional considerations (\textit{usrdef\_zgr})] 1591 {Additional considerations (\mdl{usrdef\_zgr})} 1592 \label{subsubsec:WAD_additional} 1593 1594 In the very shallow water where wetting and drying occurs the parametrisation of 1595 bottom drag is clearly very important. In order to promote stability 1596 it is sometimes useful to calculate the bottom drag using an implicit time-stepping approach. 1597 1598 Suitable specifcation of the surface heat flux in wetting and drying domains in forced and 1599 coupled simulations needs further consideration. In order to prevent freezing or boiling 1600 in uncoupled integrations the net surface heat fluxes need to be appropriately limited. 1601 1602 %---------------------------------------------------------------------------------------- 1603 % The WAD test cases 1604 %---------------------------------------------------------------------------------------- 1605 \subsection [The WAD test cases (\textit{usrdef\_zgr})] 1606 {The WAD test cases (\mdl{usrdef\_zgr})} 1607 \label{WAD_test_cases} 1608 1609 See the WAD tests MY\_DOC documention for details of the WAD test cases. 1610 1611 1612 1613 % ================================================================ 1322 1614 % Time evolution term 1323 1615 % ================================================================
Note: See TracChangeset
for help on using the changeset viewer.