New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 10499 for NEMO/trunk/doc/latex – NEMO

Changeset 10499 for NEMO/trunk/doc/latex


Ignore:
Timestamp:
2019-01-10T16:12:24+01:00 (5 years ago)
Author:
deazer
Message:

Fix ticket #2154

Location:
NEMO/trunk/doc/latex/NEMO
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/doc/latex/NEMO/main/NEMO_manual.bib

    r10495 r10499  
    29462946} 
    29472947 
     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 
    29482961@PHDTHESIS{Olivier_PhD01, 
    29492962  author =       {F. Olivier}, 
     
    37553768} 
    37563769 
     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 
    37573783@ARTICLE{Weatherly_JMR84, 
    37583784  author =       {G. L. Weatherly}, 
  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex

    r10468 r10499  
    13201320 
    13211321% ================================================================ 
     1322% Wetting and drying  
     1323% ================================================================ 
     1324\section{Wetting and drying } 
     1325\label{sec:DYN_wetdry} 
     1326There are two main options for wetting and drying code (wd): 
     1327(a) an iterative limiter (il) and (b) a directional limiter (dl). 
     1328The directional limiter is based on the scheme developed by \cite{WarnerEtal13} for RO 
     1329MS 
     1330which was in turn based on ideas developed for POM by \cite{Oey06}. The iterative 
     1331limiter is a new scheme.  The iterative limiter is activated by setting $\mathrm{ln\_wd\_il} = \mathrm{.true.}$ 
     1332and $\mathrm{ln\_wd\_dl} = \mathrm{.false.}$. The directional limiter is activated 
     1333by setting $\mathrm{ln\_wd\_dl} = \mathrm{.true.}$ and $\mathrm{ln\_wd\_il} = \mathrm{.false.}$. 
     1334 
     1335\nlst{namwad} 
     1336 
     1337The following terminology is used. The depth of the topography (positive downwards) 
     1338at each $(i,j)$ point is the quantity stored in array $\mathrm{ht\_wd}$ in the NEMO code. 
     1339The height of the free surface (positive upwards) is denoted by $ \mathrm{ssh}$. Given the sign 
     1340conventions used, the water depth, $h$, is the height of the free surface plus the depth of the 
     1341topography (i.e. $\mathrm{ssh} + \mathrm{ht\_wd}$). 
     1342 
     1343Both wd schemes take all points in the domain below a land elevation of $\mathrm{rn\_wdld}$ to be 
     1344covered by water. They require the topography specified with a model 
     1345configuration to have negative depths at points where the land is higher than the 
     1346topography's reference sea-level. The vertical grid in NEMO is normally computed relative to an 
     1347initial state with zero sea surface height elevation. 
     1348The user can choose to compute the vertical grid and heights in the model relative to 
     1349a 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 
     1352Points 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 
     1354with very steep slopes require larger values for normal choices of time-step. Surface fluxes 
     1355are also switched off for dry cells to prevent freezing, boiling etc. of very thin water layers. 
     1356The fluxes are tappered down using a $\mathrm{tanh}$ weighting function 
     1357to no flux as the dry limit $\mathrm{rn\_wdmin1}$ is approached. Even wet cells can be very shallow. 
     1358The depth at which to start tapering is controlled by the user by setting $\mathrm{rn\_wd\_sbcdep}$. 
     1359The fraction $(<1)$ of sufrace fluxes to use at this depth is set by $\mathrm{rn\_wd\_sbcfra}$. 
     1360 
     1361Both versions of the code have been tested in six test cases provided in the WAD\_TEST\_CASES configuration 
     1362and in ``realistic'' configurations covering parts of the north-west European shelf. 
     1363All these configurations have used pure sigma coordinates. It is expected that 
     1364the wetting and drying code will work in domains with more general s-coordinates provided 
     1365the coordinates are pure sigma in the region where wetting and drying actually occurs.  
     1366 
     1367The next sub-section descrbies the directional limiter and the following sub-section the iterative limiter. 
     1368The 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} 
     1377The principal idea of the directional limiter is that 
     1378water should not be allowed to flow out of a dry tracer cell (i.e. one whose water depth is less than rn\_wdmin1). 
     1379 
     1380All the changes associated with this option are made to the barotropic solver for the non-linear 
     1381free surface code within dynspg\_ts. 
     1382On each barotropic sub-step the scheme determines the direction of the flow across each face of all the tracer cells 
     1383and sets the flux across the face to zero when the flux is from a dry tracer cell. This prevents cells 
     1384whose depth is rn\_wdmin1 or less from drying out further. The scheme does not force $h$ (the water depth) at tracer cells 
     1385to be at least the minimum depth and hence is able to conserve mass / volume. 
     1386 
     1387The flux across each $u$-face of a tracer cell is multiplied by a factor zuwdmask (an array which depends on ji and jj). 
     1388If the user sets ln\_wd\_dl\_ramp = .False. then zuwdmask is 1 when the 
     1389flux is from a cell with water depth greater than rn\_wdmin1 and 0 otherwise. If the user sets 
     1390ln\_wd\_dl\_ramp = .True. the flux across the face is ramped down as the water depth decreases 
     1391from 2 * rn\_wdmin1 to rn\_wdmin1. The use of this ramp reduced grid-scale noise in idealised test cases. 
     1392 
     1393At the point where the flux across a $u$-face is multiplied by zuwdmask , we have chosen 
     1394also to multiply the corresponding velocity on the ``now'' step at that face by zuwdmask. We could have 
     1395chosen not to do that and to allow fairly large velocities to occur in these ``dry'' cells. 
     1396The rationale for setting the velocity to zero is that it is the momentum equations that are being solved 
     1397and the total momentum of the upstream cell (treating it as a finite volume) should be considered 
     1398to be its depth times its velocity. This depth is considered to be zero at ``dry'' $u$-points consistent with its 
     1399treatment 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 
     1403timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic sub-steps is respectively less than 
     1404or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer 
     1405fields (tracers independent of $x$, $y$ and $z$). Our scheme conserves constant tracers because 
     1406the velocities used at the tracer cell faces on the baroclinic timesteps are carefully calculated by dynspg\_ts 
     1407to equal their mean value during the barotropic steps. If the user sets ln\_wd\_dl\_bc = .True., the 
     1408baroclinic 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 
     1422The iterative limiter modifies the fluxes across the faces of cells that are either already ``dry'' 
     1423or may become dry within the next time-step using an iterative method. 
     1424 
     1425The flux limiter for the barotropic flow (devised by Hedong Liu) can be understood as follows: 
     1426 
     1427The 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} 
     1431can 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 
     1439In 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 
     1442the new timestep, $t_e$ is the old timestep (either $t_b$ or $t_n$) and $ \Delta t = 
     1443t_{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 
     1446The 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 
     1453The flux limiter iteratively adjusts the fluxes $\mathrm{flxu}$ and $\mathrm{flxv}$ until 
     1454none of the cells will ``dry out''. To be precise the fluxes are limited until none of the 
     1455cells has water depth less than $\mathrm{rn\_wdmin1}$ on step $n+1$. 
     1456 
     1457Let 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 
     1468where the coefficients are $1.0$ generally but can vary between $0.0$ and $1.0$ around 
     1469cells that would otherwise dry. 
     1470 
     1471The 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 
     1477The fluxes out of cell $(i,j)$ are updated at the $m+1$th iteration if the depth of the 
     1478cell on timestep $t_e$, namely $h_{i,j}(t_e)$, is less than the total flux out of the cell 
     1479times the timestep divided by the cell area. Using (\ref{dyn_wd_continuity_2}) this 
     1480condition is 
     1481 
     1482\begin{equation} \label{dyn_wd_continuity_if} 
     1483h_{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 
     1486Rearranging (\ref{dyn_wd_continuity_if}) we can obtain an expression for the maximum 
     1487outward 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 
     1496Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\it [Q: Why is 
     1497this necessary/desirable?]}. Substituting from (\ref{dyn_wd_continuity_coef}) gives an 
     1498expression for the coefficient needed to multiply the outward flux at this cell in order 
     1499to 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 
     1508Only the outward flux components are altered but, of course, outward fluxes from one cell 
     1509are inward fluxes to adjacent cells and the balance in these cells may need subsequent 
     1510adjustment; hence the iterative nature of this scheme.  Note, for example, that the flux 
     1511across the ``eastern'' face of the $(i,j)$th cell is only updated at the $m+1$th iteration 
     1512if that flux at the $m$th iteration is out of the $(i,j)$th cell. If that is the case then 
     1513the flux across that face is into the $(i+1,j)$ cell and that flux will not be updated by 
     1514the calculation for the $(i+1,j)$th cell. In this sense the updates to the fluxes across 
     1515the faces of the cells do not ``compete'' (they do not over-write each other) and one 
     1516would expect the scheme to converge relatively quickly. The scheme is flux based so 
     1517conserves mass. It also conserves constant tracers for the same reason that the 
     1518directional 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 
     1528At ``dry'' points the water depth is usually close to $\mathrm{rn\_wdmin1}$. If the 
     1529topography is sloping at these points the sea-surface will have a similar slope and there 
     1530will hence be very large horizontal pressure gradients at these points. The WAD modifies 
     1531the magnitude but not the sign of the surface pressure gradients (zhpi and zhpj) at such 
     1532points by mulitplying them by positive factors (zcpx and zcpy respectively) that lie 
     1533between $0$ and $1$. 
     1534 
     1535We describe how the scheme works for the ``eastward'' pressure gradient, zhpi, calculated 
     1536at the $(i,j)$th $u$-point. The scheme uses the ht\_wd depths and surface heights at the 
     1537neighbouring $(i+1,j)$ and $(i,j)$ tracer points.  zcpx is calculated using two logicals 
     1538variables, $\mathrm{ll\_tmp1}$ and $\mathrm{ll\_tmp2}$ which are evaluated for each grid 
     1539column.  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} 
     1545Illustrations of the three possible combinations of the logical variables controlling the 
     1546limiting of the horizontal pressure gradient in wetting and drying regimes} 
     1547\end{center}\end{figure} 
     1548%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     1549 
     1550The first logical, $\mathrm{ll\_tmp1}$, is set to true if and only if the water depth at 
     1551both neighbouring points is greater than $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ and 
     1552the minimum height of the sea surface at the two points is greater than the maximum height 
     1553of 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 
     1565The second logical, $\mathrm{ll\_tmp2}$, is set to true if and only if the maximum height 
     1566of the sea surface at the two points is greater than the maximum height of the topography 
     1567at 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 
     1577If $\mathrm{ll\_tmp1}$ is true then the surface pressure gradient, zhpi at the $(i,j)$ 
     1578point is unmodified. If both logicals are false zhpi is set to zero. 
     1579 
     1580If $\mathrm{ll\_tmp1}$ is true and $\mathrm{ll\_tmp2}$ is false then the surface pressure 
     1581gradient is multiplied through by zcpx which is the absolute value of the difference in 
     1582the water depths at the two points divided by the difference in the surface heights at the 
     1583two points. Thus the sign of the sea surface height gradient is retained but the magnitude 
     1584of the pressure force is determined by the difference in water depths rather than the 
     1585difference in surface height between the two points. Note that dividing by the difference 
     1586between the sea surface heights can be problematic if the heights approach parity. An 
     1587additional condition is applied to $\mathrm{ ll\_tmp2 }$ to ensure it is .false. in such 
     1588conditions. 
     1589 
     1590\subsubsection   [Additional considerations (\textit{usrdef\_zgr})] 
     1591         {Additional considerations (\mdl{usrdef\_zgr})} 
     1592\label{subsubsec:WAD_additional} 
     1593 
     1594In the very shallow water where wetting and drying occurs the parametrisation of 
     1595bottom drag is clearly very important. In order to promote stability 
     1596it is sometimes useful to calculate the bottom drag using an implicit time-stepping approach. 
     1597 
     1598Suitable specifcation of the surface heat flux in wetting and drying domains in forced and 
     1599coupled simulations needs further consideration. In order to prevent freezing or boiling 
     1600in 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 
     1609See the WAD tests MY\_DOC documention for details of the WAD test cases. 
     1610 
     1611 
     1612 
     1613% ================================================================ 
    13221614% Time evolution term  
    13231615% ================================================================ 
Note: See TracChangeset for help on using the changeset viewer.