Index: /NEMO/trunk/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/trunk/cfgs/SHARED/namelist_ref	(revision 10498)
+++ /NEMO/trunk/cfgs/SHARED/namelist_ref	(revision 10499)
@@ -116,5 +116,5 @@
    ln_wd_dl    = .false.   !  T/F activation of directional limiter
    ln_wd_dl_bc = .false.   !  T/F Directional limiteer Baroclinic option
-   ln_wd_dl_rmp = .false.   !  T/F Turn on directional limiter ramp
+   ln_wd_dl_rmp = .false.  !  T/F Turn on directional limiter ramp
    rn_wdmin0   =  0.30     !  depth at which WaD starts
    rn_wdmin1   =  0.2      !  Minimum wet depth on dried cells
@@ -122,4 +122,6 @@
    rn_wdld     =  2.5      !  Land elevation below which WaD is allowed
    nn_wdit     =   20      !  Max iterations for WaD limiter
+   rn_wd_sbcdep =  5.0     !  Depth at which to taper sbc fluxes
+   rn_wd_sbcfra =  0.999   !  Fraction of SBC fluxes at taper depth (Must be <1)
 /
 !-----------------------------------------------------------------------
Index: /NEMO/trunk/doc/latex/NEMO/main/NEMO_manual.bib
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/main/NEMO_manual.bib	(revision 10498)
+++ /NEMO/trunk/doc/latex/NEMO/main/NEMO_manual.bib	(revision 10499)
@@ -2946,4 +2946,17 @@
 }
 
+@ARTICLE{Oey06,
+   author =      {Lie-Yauw Oey},
+   title =       {An OGCM with movable land-sea boundaries},
+   journal =     {Ocean Modelling},
+   volume =      "13",
+   number =      "2",
+   pages =       "176 - 195",
+   year =        "2006",
+   issn =        "1463-5003",
+   doi =         "https://doi.org/10.1016/j.ocemod.2006.01.001",
+   url =         "http://www.sciencedirect.com/science/article/pii/S1463500306000084",
+}
+
 @PHDTHESIS{Olivier_PhD01,
   author =       {F. Olivier},
@@ -3755,4 +3768,17 @@
 }
 
+@ARTICLE{WarnerEtal13,
+   author =      {John C. Warner and Zafer Defne and Kevin Haas and
+                  Hernan G. Arango},
+   title =       {A wetting and drying scheme for ROMS},
+   journal =     "Computers \& Geosciences",
+   volume =      "58", 
+   pages =       "54 - 61",
+   year =        "2013",
+   issn =        "0098-3004",
+   doi =         "https://doi.org/10.1016/j.cageo.2013.05.004",
+   url =         "http://www.sciencedirect.com/science/article/pii/S0098300413001362",
+}
+
 @ARTICLE{Weatherly_JMR84,
   author =       {G. L. Weatherly},
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex	(revision 10498)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex	(revision 10499)
@@ -1320,4 +1320,296 @@
 
 % ================================================================
+% Wetting and drying 
+% ================================================================
+\section{Wetting and drying }
+\label{sec:DYN_wetdry}
+There are two main options for wetting and drying code (wd):
+(a) an iterative limiter (il) and (b) a directional limiter (dl).
+The directional limiter is based on the scheme developed by \cite{WarnerEtal13} for RO
+MS
+which was in turn based on ideas developed for POM by \cite{Oey06}. The iterative
+limiter is a new scheme.  The iterative limiter is activated by setting $\mathrm{ln\_wd\_il} = \mathrm{.true.}$
+and $\mathrm{ln\_wd\_dl} = \mathrm{.false.}$. The directional limiter is activated
+by setting $\mathrm{ln\_wd\_dl} = \mathrm{.true.}$ and $\mathrm{ln\_wd\_il} = \mathrm{.false.}$.
+
+\nlst{namwad}
+
+The following terminology is used. The depth of the topography (positive downwards)
+at each $(i,j)$ point is the quantity stored in array $\mathrm{ht\_wd}$ in the NEMO code.
+The height of the free surface (positive upwards) is denoted by $ \mathrm{ssh}$. Given the sign
+conventions used, the water depth, $h$, is the height of the free surface plus the depth of the
+topography (i.e. $\mathrm{ssh} + \mathrm{ht\_wd}$).
+
+Both wd schemes take all points in the domain below a land elevation of $\mathrm{rn\_wdld}$ to be
+covered by water. They require the topography specified with a model
+configuration to have negative depths at points where the land is higher than the
+topography's reference sea-level. The vertical grid in NEMO is normally computed relative to an
+initial state with zero sea surface height elevation.
+The user can choose to compute the vertical grid and heights in the model relative to
+a non-zero reference height for the free surface. This choice affects the calculation of the metrics and depths
+(i.e. the $\mathrm{e3t\_0, ht\_0}$ etc. arrays).
+
+Points where the water depth is less than $\mathrm{rn\_wdmin1}$ are interpreted as ``dry''.
+$\mathrm{rn\_wdmin1}$ is usually chosen to be of order $0.05$m but extreme topographies
+with very steep slopes require larger values for normal choices of time-step. Surface fluxes
+are also switched off for dry cells to prevent freezing, boiling etc. of very thin water layers.
+The fluxes are tappered down using a $\mathrm{tanh}$ weighting function
+to no flux as the dry limit $\mathrm{rn\_wdmin1}$ is approached. Even wet cells can be very shallow.
+The depth at which to start tapering is controlled by the user by setting $\mathrm{rn\_wd\_sbcdep}$.
+The fraction $(<1)$ of sufrace fluxes to use at this depth is set by $\mathrm{rn\_wd\_sbcfra}$.
+
+Both versions of the code have been tested in six test cases provided in the WAD\_TEST\_CASES configuration
+and in ``realistic'' configurations covering parts of the north-west European shelf.
+All these configurations have used pure sigma coordinates. It is expected that
+the wetting and drying code will work in domains with more general s-coordinates provided
+the coordinates are pure sigma in the region where wetting and drying actually occurs. 
+
+The next sub-section descrbies the directional limiter and the following sub-section the iterative limiter.
+The final sub-section covers some additional considerations that are relevant to both schemes.
+
+
+%-----------------------------------------------------------------------------------------
+%   Iterative limiters
+%-----------------------------------------------------------------------------------------
+\subsection   [Directional limiter (\textit{wet\_dry})]
+         {Directional limiter (\mdl{wet\_dry})}
+\label{subsec:DYN_wd_directional_limiter}
+The principal idea of the directional limiter is that
+water should not be allowed to flow out of a dry tracer cell (i.e. one whose water depth is less than rn\_wdmin1).
+
+All the changes associated with this option are made to the barotropic solver for the non-linear
+free surface code within dynspg\_ts.
+On each barotropic sub-step the scheme determines the direction of the flow across each face of all the tracer cells
+and sets the flux across the face to zero when the flux is from a dry tracer cell. This prevents cells
+whose depth is rn\_wdmin1 or less from drying out further. The scheme does not force $h$ (the water depth) at tracer cells
+to be at least the minimum depth and hence is able to conserve mass / volume.
+
+The flux across each $u$-face of a tracer cell is multiplied by a factor zuwdmask (an array which depends on ji and jj).
+If the user sets ln\_wd\_dl\_ramp = .False. then zuwdmask is 1 when the
+flux is from a cell with water depth greater than rn\_wdmin1 and 0 otherwise. If the user sets
+ln\_wd\_dl\_ramp = .True. the flux across the face is ramped down as the water depth decreases
+from 2 * rn\_wdmin1 to rn\_wdmin1. The use of this ramp reduced grid-scale noise in idealised test cases.
+
+At the point where the flux across a $u$-face is multiplied by zuwdmask , we have chosen
+also to multiply the corresponding velocity on the ``now'' step at that face by zuwdmask. We could have
+chosen not to do that and to allow fairly large velocities to occur in these ``dry'' cells.
+The rationale for setting the velocity to zero is that it is the momentum equations that are being solved
+and the total momentum of the upstream cell (treating it as a finite volume) should be considered
+to be its depth times its velocity. This depth is considered to be zero at ``dry'' $u$-points consistent with its
+treatment in the calculation of the flux of mass across the cell face.
+
+
+\cite{WarnerEtal13} state that in their scheme the velocity masks at the cell faces for the baroclinic
+timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic sub-steps is respectively less than
+or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer
+fields (tracers independent of $x$, $y$ and $z$). Our scheme conserves constant tracers because
+the velocities used at the tracer cell faces on the baroclinic timesteps are carefully calculated by dynspg\_ts
+to equal their mean value during the barotropic steps. If the user sets ln\_wd\_dl\_bc = .True., the
+baroclinic velocities are also multiplied by a suitably weighted average of zuwdmask.  
+
+%-----------------------------------------------------------------------------------------
+%   Iterative limiters
+%-----------------------------------------------------------------------------------------
+
+\subsection   [Iterative limiter (\textit{wet\_dry})]
+         {Iterative limiter (\mdl{wet\_dry})}
+\label{subsec:DYN_wd_iterative_limiter}
+
+\subsubsection [Iterative flux limiter (\textit{wet\_dry})]
+         {Iterative flux limiter (\mdl{wet\_dry})}
+\label{subsubsec:DYN_wd_il_spg_limiter}
+
+The iterative limiter modifies the fluxes across the faces of cells that are either already ``dry''
+or may become dry within the next time-step using an iterative method.
+
+The flux limiter for the barotropic flow (devised by Hedong Liu) can be understood as follows:
+
+The continuity equation for the total water depth in a column
+\begin{equation} \label{dyn_wd_continuity}
+ \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 .
+\end{equation}
+can be written in discrete form  as
+
+\begin{align} \label{dyn_wd_continuity_2}
+\frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) ) 
+&= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j}  + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\
+&= \mathrm{zzflx}_{i,j} .
+\end{align}
+
+In the above $h$ is the depth of the water in the column at point $(i,j)$,
+$\mathrm{flxu}_{i+1,j}$ is the flux out of the ``eastern'' face of the cell and
+$\mathrm{flxv}_{i,j+1}$ the flux out of the ``northern'' face of the cell; $t_{n+1}$ is
+the new timestep, $t_e$ is the old timestep (either $t_b$ or $t_n$) and $ \Delta t =
+t_{n+1} - t_e$; $e_1 e_2$ is the area of the tracer cells centred at $(i,j)$ and
+$\mathrm{zzflx}$ is the sum of the fluxes through all the faces.
+
+The flux limiter splits the flux $\mathrm{zzflx}$ into fluxes that are out of the cell
+(zzflxp) and fluxes that are into the cell (zzflxn).  Clearly
+
+\begin{equation} \label{dyn_wd_zzflx_p_n_1}
+\mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} .  
+\end{equation}
+
+The flux limiter iteratively adjusts the fluxes $\mathrm{flxu}$ and $\mathrm{flxv}$ until
+none of the cells will ``dry out''. To be precise the fluxes are limited until none of the
+cells has water depth less than $\mathrm{rn\_wdmin1}$ on step $n+1$.
+
+Let the fluxes on the $m$th iteration step be denoted by $\mathrm{flxu}^{(m)}$ and
+$\mathrm{flxv}^{(m)}$.  Then the adjustment is achieved by seeking a set of coefficients,
+$\mathrm{zcoef}_{i,j}^{(m)}$ such that:
+
+\begin{equation} \label{dyn_wd_continuity_coef}
+\begin{split}
+\mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\
+\mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j}
+\end{split}
+\end{equation}
+
+where the coefficients are $1.0$ generally but can vary between $0.0$ and $1.0$ around
+cells that would otherwise dry.
+
+The iteration is initialised by setting
+
+\begin{equation} \label{dyn_wd_zzflx_initial}
+\mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad  \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} . 
+\end{equation}
+
+The fluxes out of cell $(i,j)$ are updated at the $m+1$th iteration if the depth of the
+cell on timestep $t_e$, namely $h_{i,j}(t_e)$, is less than the total flux out of the cell
+times the timestep divided by the cell area. Using (\ref{dyn_wd_continuity_2}) this
+condition is
+
+\begin{equation} \label{dyn_wd_continuity_if}
+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} ) .
+\end{equation}
+
+Rearranging (\ref{dyn_wd_continuity_if}) we can obtain an expression for the maximum
+outward flux that can be allowed and still maintain the minimum wet depth:
+
+\begin{equation} \label{dyn_wd_max_flux}
+\begin{split}
+\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{]} \\
+\phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big]
+\end{split}
+\end{equation}
+
+Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\it [Q: Why is
+this necessary/desirable?]}. Substituting from (\ref{dyn_wd_continuity_coef}) gives an
+expression for the coefficient needed to multiply the outward flux at this cell in order
+to avoid drying.
+
+\begin{equation} \label{dyn_wd_continuity_nxtcoef}
+\begin{split}
+\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{]} \\
+\phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} } 
+\end{split}
+\end{equation}
+
+Only the outward flux components are altered but, of course, outward fluxes from one cell
+are inward fluxes to adjacent cells and the balance in these cells may need subsequent
+adjustment; hence the iterative nature of this scheme.  Note, for example, that the flux
+across the ``eastern'' face of the $(i,j)$th cell is only updated at the $m+1$th iteration
+if that flux at the $m$th iteration is out of the $(i,j)$th cell. If that is the case then
+the flux across that face is into the $(i+1,j)$ cell and that flux will not be updated by
+the calculation for the $(i+1,j)$th cell. In this sense the updates to the fluxes across
+the faces of the cells do not ``compete'' (they do not over-write each other) and one
+would expect the scheme to converge relatively quickly. The scheme is flux based so
+conserves mass. It also conserves constant tracers for the same reason that the
+directional limiter does.
+
+
+%----------------------------------------------------------------------------------------
+%      Surface pressure gradients
+%----------------------------------------------------------------------------------------
+\subsubsection   [Modification of surface pressure gradients (\textit{dynhpg})]
+         {Modification of surface pressure gradients (\mdl{dynhpg})}
+\label{subsubsec:DYN_wd_il_spg}
+
+At ``dry'' points the water depth is usually close to $\mathrm{rn\_wdmin1}$. If the
+topography is sloping at these points the sea-surface will have a similar slope and there
+will hence be very large horizontal pressure gradients at these points. The WAD modifies
+the magnitude but not the sign of the surface pressure gradients (zhpi and zhpj) at such
+points by mulitplying them by positive factors (zcpx and zcpy respectively) that lie
+between $0$ and $1$.
+
+We describe how the scheme works for the ``eastward'' pressure gradient, zhpi, calculated
+at the $(i,j)$th $u$-point. The scheme uses the ht\_wd depths and surface heights at the
+neighbouring $(i+1,j)$ and $(i,j)$ tracer points.  zcpx is calculated using two logicals
+variables, $\mathrm{ll\_tmp1}$ and $\mathrm{ll\_tmp2}$ which are evaluated for each grid
+column.  The three possible combinations are illustrated in figure \ref{Fig_WAD_dynhpg}.
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!ht] \begin{center}
+\includegraphics[width=0.8\textwidth]{Fig_WAD_dynhpg}
+\caption{ \label{Fig_WAD_dynhpg}
+Illustrations of the three possible combinations of the logical variables controlling the
+limiting of the horizontal pressure gradient in wetting and drying regimes}
+\end{center}\end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+The first logical, $\mathrm{ll\_tmp1}$, is set to true if and only if the water depth at
+both neighbouring points is greater than $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ and
+the minimum height of the sea surface at the two points is greater than the maximum height
+of the topography at the two points:
+
+\begin{equation} \label{dyn_ll_tmp1}
+\begin{split}
+\mathrm{ll\_tmp1}  = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\
+                     & \quad \mathrm{MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj))\  .and.} \\
+& \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\
+& \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\
+& \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 }
+\end{split}
+\end{equation}
+
+The second logical, $\mathrm{ll\_tmp2}$, is set to true if and only if the maximum height
+of the sea surface at the two points is greater than the maximum height of the topography
+at the two points plus $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$
+
+\begin{equation} \label{dyn_ll_tmp2}
+\begin{split}
+\mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\
+& \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\
+& \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) .
+\end{split}
+\end{equation}
+
+If $\mathrm{ll\_tmp1}$ is true then the surface pressure gradient, zhpi at the $(i,j)$
+point is unmodified. If both logicals are false zhpi is set to zero.
+
+If $\mathrm{ll\_tmp1}$ is true and $\mathrm{ll\_tmp2}$ is false then the surface pressure
+gradient is multiplied through by zcpx which is the absolute value of the difference in
+the water depths at the two points divided by the difference in the surface heights at the
+two points. Thus the sign of the sea surface height gradient is retained but the magnitude
+of the pressure force is determined by the difference in water depths rather than the
+difference in surface height between the two points. Note that dividing by the difference
+between the sea surface heights can be problematic if the heights approach parity. An
+additional condition is applied to $\mathrm{ ll\_tmp2 }$ to ensure it is .false. in such
+conditions.
+
+\subsubsection   [Additional considerations (\textit{usrdef\_zgr})]
+         {Additional considerations (\mdl{usrdef\_zgr})}
+\label{subsubsec:WAD_additional}
+
+In the very shallow water where wetting and drying occurs the parametrisation of
+bottom drag is clearly very important. In order to promote stability
+it is sometimes useful to calculate the bottom drag using an implicit time-stepping approach.
+
+Suitable specifcation of the surface heat flux in wetting and drying domains in forced and
+coupled simulations needs further consideration. In order to prevent freezing or boiling
+in uncoupled integrations the net surface heat fluxes need to be appropriately limited.
+
+%----------------------------------------------------------------------------------------
+%      The WAD test cases
+%----------------------------------------------------------------------------------------
+\subsection   [The WAD test cases (\textit{usrdef\_zgr})]
+         {The WAD test cases (\mdl{usrdef\_zgr})}
+\label{WAD_test_cases}
+
+See the WAD tests MY\_DOC documention for details of the WAD test cases.
+
+
+
+% ================================================================
 % Time evolution term 
 % ================================================================
Index: /NEMO/trunk/doc/namelists/namwad
===================================================================
--- /NEMO/trunk/doc/namelists/namwad	(revision 10498)
+++ /NEMO/trunk/doc/namelists/namwad	(revision 10499)
@@ -1,13 +1,15 @@
 !-----------------------------------------------------------------------
-&namwad        !   Wetting and Drying (WaD)                             (default: OFF)
+&namwad  !   Wetting and drying
 !-----------------------------------------------------------------------
-   ln_wd_il    = .false    !  T/F activation of iterative   limiter
-   ln_wd_dl    = .false.   !  T/F activation of directional limiter
-   ln_wd_dl_bc = .false.   !  T/F Directional limiteer Baroclinic option
-   ln_wd_dl_rmp = .false.   !  T/F Turn on directional limiter ramp
-   rn_wdmin0   =  0.30     !  depth at which WaD starts
-   rn_wdmin1   =  0.2      !  Minimum wet depth on dried cells
-   rn_wdmin2   =  0.0001   !  Tolerance of min wet depth on dried cells
-   rn_wdld     =  2.5      !  Land elevation below which WaD is allowed
-   nn_wdit     =   20      !  Max iterations for WaD limiter
+   ln_wd_il          = .false   ! T/F activation of iterative limiter for  wetting and drying scheme
+   ln_wd_dl          = .true.   ! T/F activation of directional llimiter for wetting drying scheme
+   ln_wd_dl_bc       = .true.   ! T/F Directional limiteer Baroclinic option
+   ln_wd_dl_rmp      = .true.   ! T/F Turn on directional limiter ramp
+   rn_wdmin0         =  0.30    ! dpoth at which wetting/drying starts
+   rn_wdmin1         =  0.2     ! Minimum wet depth on dried cells
+   rn_wdmin2         =  0.0001  ! Tolerance of min wet depth on dried cells
+   rn_wdld           =  2.5     ! Land elevation below which wetting/drying is allowed
+   nn_wdit           =   20     ! Max iterations for W/D limiter
+   rn_wd_sbcdep      =  5.0     ! Depth at which to taper sbc fluxes
+   rn_wd_sbcfra      =  0.999   ! Fraction of SBC fluxes at taper depth (Must be <1)
 /
Index: /NEMO/trunk/src/OCE/DIA/dia25h.F90
===================================================================
--- /NEMO/trunk/src/OCE/DIA/dia25h.F90	(revision 10498)
+++ /NEMO/trunk/src/OCE/DIA/dia25h.F90	(revision 10499)
@@ -13,4 +13,5 @@
    USE in_out_manager  ! I/O units
    USE iom             ! I/0 library
+   USE wet_dry
 
    IMPLICIT NONE
@@ -211,5 +212,9 @@
          CALL iom_put( "salin25h", zw3d  )   ! salinity
          zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
-         CALL iom_put( "ssh25h", zw2d )   ! sea surface 
+         IF( ll_wd ) THEN
+            CALL iom_put( "ssh25h", zw2d+ssh_ref )   ! sea surface 
+         ELSE
+            CALL iom_put( "ssh25h", zw2d )   ! sea surface
+         ENDIF
          ! Write velocities (instantaneous)
          zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:))
Index: /NEMO/trunk/src/OCE/DIA/diatmb.F90
===================================================================
--- /NEMO/trunk/src/OCE/DIA/diatmb.F90	(revision 10498)
+++ /NEMO/trunk/src/OCE/DIA/diatmb.F90	(revision 10499)
@@ -12,4 +12,5 @@
    USE in_out_manager  ! I/O units
    USE iom             ! I/0 library
+   USE wet_dry
 
    IMPLICIT NONE
@@ -108,5 +109,10 @@
       CALL dia_calctmb( tsn(:,:,:,jp_tem), zwtmb )
       !ssh already output but here we output it masked
-      CALL iom_put( "sshnmasked", sshn(:,:)*tmask(:,:,1) + zmdi*(1.0 - tmask(:,:,1)) )
+      IF( ll_wd ) THEN
+         CALL iom_put( "sshnmasked", (sshn(:,:)+ssh_ref)*tmask(:,:,1) + zmdi*(1.0 - tmask(:,:,1)) )
+      ELSE
+         CALL iom_put( "sshnmasked", sshn(:,:)*tmask(:,:,1) + zmdi*(1.0 - tmask(:,:,1)) )
+      ENDIF
+
       CALL iom_put( "top_temp"  , zwtmb(:,:,1) )    ! tmb Temperature
       CALL iom_put( "mid_temp"  , zwtmb(:,:,2) )    ! tmb Temperature
Index: /NEMO/trunk/src/OCE/DOM/istate.F90
===================================================================
--- /NEMO/trunk/src/OCE/DOM/istate.F90	(revision 10498)
+++ /NEMO/trunk/src/OCE/DOM/istate.F90	(revision 10499)
@@ -100,4 +100,17 @@
             !
             sshb(:,:)   = 0._wp               ! set the ocean at rest
+            IF( ll_wd ) THEN
+               sshb(:,:) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD 
+               !
+               ! Apply minimum wetdepth criterion
+               !
+               DO jj = 1,jpj
+                  DO ji = 1,jpi
+                     IF( ht_0(ji,jj) + sshb(ji,jj)  < rn_wdmin1 ) THEN
+                        sshb(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
+                     ENDIF
+                  END DO
+               END DO 
+            ENDIF 
             ub  (:,:,:) = 0._wp
             vb  (:,:,:) = 0._wp  
Index: /NEMO/trunk/src/OCE/DYN/wet_dry.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/wet_dry.F90	(revision 10498)
+++ /NEMO/trunk/src/OCE/DYN/wet_dry.F90	(revision 10499)
@@ -45,4 +45,6 @@
    REAL(wp), PUBLIC  ::   r_rn_wdmin1 !: 1/minimum water depth on dried cells 
    REAL(wp), PUBLIC  ::   rn_wdmin2   !: tolerance of minimum water depth on dried cells
+   REAL(wp), PUBLIC  ::   rn_wd_sbcdep   !: Depth at which to taper sbc fluxes
+   REAL(wp), PUBLIC  ::   rn_wd_sbcfra   !: Fraction of SBC at taper depth
    REAL(wp), PUBLIC  ::   rn_wdld     !: land elevation below which wetting/drying will be considered
    INTEGER , PUBLIC  ::   nn_wdit     !: maximum number of iteration for W/D limiter
@@ -74,5 +76,5 @@
       !!
       NAMELIST/namwad/ ln_wd_il, ln_wd_dl   , rn_wdmin0, rn_wdmin1, rn_wdmin2, rn_wdld,   &
-         &             nn_wdit , ln_wd_dl_bc, ln_wd_dl_rmp
+         &             nn_wdit , ln_wd_dl_bc, ln_wd_dl_rmp, rn_wd_sbcdep,rn_wd_sbcfra
       !!----------------------------------------------------------------------
       !
@@ -85,4 +87,5 @@
       IF(lwm) WRITE ( numond, namwad )
       !
+      IF( rn_wd_sbcfra>=1 )   CALL ctl_stop( 'STOP', 'rn_wd_sbcfra >=1 : must be < 1' )
       IF(lwp) THEN                  ! control print
          WRITE(numout,*)
@@ -99,4 +102,6 @@
          WRITE(numout,*) '      T => baroclinic u,v=0 at dry pts: ln_wd_dl_bc = ', ln_wd_dl_bc     
          WRITE(numout,*) '      use a ramp for rwd limiter:  ln_wd_dl_rwd_rmp = ', ln_wd_dl_rmp
+         WRITE(numout,*) '      cut off depth sbc for wd   rn_wd_sbcdep       = ', rn_wd_sbcdep
+         WRITE(numout,*) '      fraction to start sbc wgt rn_wd_sbcfra        = ', rn_wd_sbcfra
       ENDIF
       IF( .NOT. ln_read_cfg ) THEN
Index: /NEMO/trunk/src/OCE/SBC/sbcmod.F90
===================================================================
--- /NEMO/trunk/src/OCE/SBC/sbcmod.F90	(revision 10498)
+++ /NEMO/trunk/src/OCE/SBC/sbcmod.F90	(revision 10499)
@@ -58,4 +58,5 @@
    USE lib_mpp        ! MPP library
    USE timing         ! Timing
+   USE wet_dry
    USE diurnal_bulk, ONLY:   ln_diurnal_only   ! diurnal SST diagnostic
 
@@ -377,4 +378,8 @@
       !
       LOGICAL ::   ll_sas, ll_opa   ! local logical
+      !
+      REAL(wp) ::     zthscl        ! wd  tanh scale
+      REAL(wp), DIMENSION(jpi,jpj) ::  zwdht, zwght  ! wd dep over wd limit, wgt  
+
       !!---------------------------------------------------------------------
       !
@@ -461,4 +466,28 @@
 !!$!clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why.
 !!$      CALL lbc_lnk( 'sbcmod', emp, 'T', 1. )
+      IF ( ll_wd ) THEN     ! If near WAD point limit the flux for now
+         zthscl = atanh(rn_wd_sbcfra)                     ! taper frac default is .999 
+         zwdht(:,:) = sshn(:,:) + ht_0(:,:) - rn_wdmin1   ! do this calc of water
+                                                     ! depth above wd limit once
+         WHERE( zwdht(:,:) <= 0.0 )
+            taum(:,:) = 0.0
+            utau(:,:) = 0.0
+            vtau(:,:) = 0.0
+            qns (:,:) = 0.0
+            qsr (:,:) = 0.0
+            emp (:,:) = min(emp(:,:),0.0) !can allow puddles to grow but not shrink
+            sfx (:,:) = 0.0
+         END WHERE
+         zwght(:,:) = tanh(zthscl*zwdht(:,:))
+         WHERE( zwdht(:,:) > 0.0  .and. zwdht(:,:) < rn_wd_sbcdep ) !  5 m hard limit here is arbitrary
+            qsr  (:,:) =  qsr(:,:)  * zwght(:,:)
+            qns  (:,:) =  qns(:,:)  * zwght(:,:)
+            taum (:,:) =  taum(:,:) * zwght(:,:)
+            utau (:,:) =  utau(:,:) * zwght(:,:)
+            vtau (:,:) =  vtau(:,:) * zwght(:,:)
+            sfx  (:,:) =  sfx(:,:)  * zwght(:,:)
+            emp  (:,:) =  emp(:,:)  * zwght(:,:)
+         END WHERE
+      ENDIF
       !
       IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
Index: /NEMO/trunk/src/OCE/TRA/trasbc.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/trasbc.F90	(revision 10498)
+++ /NEMO/trunk/src/OCE/TRA/trasbc.F90	(revision 10499)
@@ -27,5 +27,4 @@
    USE trd_oce        ! trends: ocean variables
    USE trdtra         ! trends manager: tracers 
-   USE wet_dry,  ONLY : ll_wd, rn_wdmin1, r_rn_wdmin1   ! Wetting and drying
 #if defined key_asminc   
    USE asminc         ! Assimilation increment
@@ -125,17 +124,5 @@
       DO jj = 2, jpj
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            IF ( ll_wd ) THEN     ! If near WAD point limit the flux for now
-               IF ( sshn(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN
-                  sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)   ! non solar heat flux
-               ELSE IF ( sshn(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN
-                  sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj) &
-                       &                * tanh ( 5._wp * ( ( sshn(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 ) * r_rn_wdmin1 ) )
-               ELSE
-                  sbc_tsc(ji,jj,jp_tem) = 0._wp
-               ENDIF
-            ELSE 
-               sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)   ! non solar heat flux
-            ENDIF
-
+            sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)   ! non solar heat flux
             sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)   ! salt flux due to freezing/melting
          END DO
Index: /NEMO/trunk/tests/WAD/MY_DOCS/Namelists/nam_wad
===================================================================
--- /NEMO/trunk/tests/WAD/MY_DOCS/Namelists/nam_wad	(revision 10498)
+++ /NEMO/trunk/tests/WAD/MY_DOCS/Namelists/nam_wad	(revision 10499)
@@ -11,3 +11,5 @@
    rn_wdld           =  2.5     ! Land elevation below which wetting/drying is allowed
    nn_wdit           =   20     ! Max iterations for W/D limiter
+   rn_wd_sbcdep      =  5.0     ! Depth at which to taper sbc fluxes
+   rn_wd_sbcfra      =  0.999   ! Fraction of SBC fluxes at taper depth (Must be <1)
 /
Index: /NEMO/trunk/tests/WAD/MY_DOCS/WAD_doc.tex
===================================================================
--- /NEMO/trunk/tests/WAD/MY_DOCS/WAD_doc.tex	(revision 10498)
+++ /NEMO/trunk/tests/WAD/MY_DOCS/WAD_doc.tex	(revision 10499)
@@ -25,278 +25,9 @@
 % Wetting and drying 
 % ================================================================
-\section{Wetting and drying }
-\label{DYN_wetdry}
-
-There are two main options for wetting and drying code (wd): 
-(a) an iterative limiter (il) and (b) a directional limiter (dl). 
-The directional limiter is based on the scheme developed by \cite{WarnerEtal13} for ROMS
-which was in turn based on ideas developed for POM by \cite{Oey06}. The iterative
-limiter is a new scheme.  The iterative limiter is activated by setting $\mathrm{ln\_wd\_il} = \mathrm{.true.}$
-and $\mathrm{ln\_wd\_dl} = \mathrm{.false.}$. The directional limiter is activated
-by setting $\mathrm{ln\_wd\_dl} = \mathrm{.true.}$ and $\mathrm{ln\_wd\_il} = \mathrm{.false.}$.
-
-\namdisplay{nam_wad}
-
-The following terminology is used. The depth of the topography (positive downwards)
-at each $(i,j)$ point is the quantity stored in array $\mathrm{ht\_wd}$ in the NEMO code.
-The height of the free surface (positive upwards) is denoted by $ \mathrm{ssh}$. Given the sign
-conventions used, the water depth, $h$, is the height of the free surface plus the depth of the
-topography (i.e. $\mathrm{ssh} + \mathrm{ht\_wd}$).
-
-Both wd schemes take all points in the domain below a land elevation of $\mathrm{rn\_wdld}$ to be
-covered by water. They require the topography specified with a model
-configuration to have negative depths at points where the land is higher than the
-topography's reference sea-level. The vertical grid in NEMO is normally computed relative to an
-initial state with zero sea surface height elevation. 
-The user can choose to compute the vertical grid and heights in the model relative to 
-a non-zero reference height for the free surface. This choice affects the calculation of the metrics and depths
-(i.e. the $\mathrm{e3t\_0, ht\_0}$ etc. arrays). 
-
-Points where the water depth is less than $\mathrm{rn\_wdmin1}$ are interpreted as ``dry''. 
-$\mathrm{rn\_wdmin1}$ is usually chosen to be of order $0.05$m but extreme topographies 
-with very steep slopes require larger values for normal choices of time-step. 
-
-Both versions of the code have been tested in six test cases provided in the WAD\_TEST\_CASES configuration
-and in ``realistic'' configurations covering parts of the north-west European shelf. 
-All these configurations have used pure sigma coordinates. It is expected that
-the wetting and drying code will work in domains with more general s-coordinates provided
-the coordinates are pure sigma in the region where wetting and drying actually occurs.  
-
-The next sub-section descrbies the directional limiter and the following sub-section the iterative limiter. 
-The final sub-section covers some additional considerations that are relevant to both schemes. 
-
-%-----------------------------------------------------------------------------------------
-%   Iterative limiters
-%-----------------------------------------------------------------------------------------
-\subsection   [Directional limiter (\textit{wet\_dry})]
-			{Directional limiter (\mdl{wet\_dry})}
-\label{DYN_wd_directional_limiter}
-
-The principal idea of the directional limiter is that 
-water should not be allowed to flow out of a dry tracer cell (i.e. one whose water depth is less than rn\_wdmin1).
-
-All the changes associated with this option are made to the barotropic solver for the non-linear 
-free surface code within dynspg\_ts. 
-On each barotropic sub-step the scheme determines the direction of the flow across each face of all the tracer cells
-and sets the flux across the face to zero when the flux is from a dry tracer cell. This prevents cells
-whose depth is rn\_wdmin1 or less from drying out further. The scheme does not force $h$ (the water depth) at tracer cells
-to be at least the minimum depth and hence is able to conserve mass / volume. 
-
-The flux across each $u$-face of a tracer cell is multiplied by a factor zuwdmask (an array which depends on ji and jj). 
-If the user sets ln\_wd\_dl\_ramp = .False. then zuwdmask is 1 when the
-flux is from a cell with water depth greater than rn\_wdmin1 and 0 otherwise. If the user sets
-ln\_wd\_dl\_ramp = .True. the flux across the face is ramped down as the water depth decreases
-from 2 * rn\_wdmin1 to rn\_wdmin1. The use of this ramp reduced grid-scale noise in idealised test cases. 
-
-At the point where the flux across a $u$-face is multiplied by zuwdmask , we have chosen
-also to multiply the corresponding velocity on the ``now'' step at that face by zuwdmask. We could have
-chosen not to do that and to allow fairly large velocities to occur in these ``dry'' cells. 
-The rationale for setting the velocity to zero is that it is the momentum equations that are being solved
-and the total momentum of the upstream cell (treating it as a finite volume) should be considered
-to be its depth times its velocity. This depth is considered to be zero at ``dry'' $u$-points consistent with its 
-treatment in the calculation of the flux of mass across the cell face.         
-
-\cite{WarnerEtal13} state that in their scheme the velocity masks at the cell faces for the baroclinic 
-timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic sub-steps is respectively less than 
-or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer
-fields (tracers independent of $x$, $y$ and $z$). Our scheme conserves constant tracers because
-the velocities used at the tracer cell faces on the baroclinic timesteps are carefully calculated by dynspg\_ts
-to equal their mean value during the barotropic steps. If the user sets ln\_wd\_dl\_bc = .True., the
-baroclinic velocities are also multiplied by a suitably weighted average of zuwdmask.      
-     
-%-----------------------------------------------------------------------------------------
-%   Iterative limiters
-%-----------------------------------------------------------------------------------------
-\subsection   [Iterative limiter (\textit{wet\_dry})]
-			{Iterative limiter (\mdl{wet\_dry})}
-\label{DYN_wd_iterative_limiter}
-
-\subsubsection [Iterative flux limiter (\textit{wet\_dry})]
-			{Iterative flux limiter (\mdl{wet\_dry})}
-\label{DYN_wd_il_spg_limiter}
-
-The iterative limiter modifies the fluxes across the faces of cells that are either already ``dry''
-or may become dry within the next time-step using an iterative method. 
-
-The flux limiter for the barotropic flow (devised by Hedong Liu) can be understood as follows: 
-
-The continuity equation for the total water depth in a column 
-\begin{equation} \label{dyn_wd_continuity}
- \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 .
-\end{equation} 
-can be written in discrete form  as  
-
-\begin{align} \label{dyn_wd_continuity_2}
-\frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) ) 
-&= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j}  + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\
-&= \mathrm{zzflx}_{i,j} .
-\end{align} 
-
-In the above $h$ is the depth of the water in the column at point $(i,j)$,
-$\mathrm{flxu}_{i+1,j}$ is the flux out of the ``eastern'' face of the cell and
-$\mathrm{flxv}_{i,j+1}$ the flux out of the ``northern'' face of the cell; $t_{n+1}$ is
-the new timestep, $t_e$ is the old timestep (either $t_b$ or $t_n$) and $ \Delta t =
-t_{n+1} - t_e$; $e_1 e_2$ is the area of the tracer cells centred at $(i,j)$ and
-$\mathrm{zzflx}$ is the sum of the fluxes through all the faces.
-
-The flux limiter splits the flux $\mathrm{zzflx}$ into fluxes that are out of the cell
-(zzflxp) and fluxes that are into the cell (zzflxn).  Clearly
-
-\begin{equation} \label{dyn_wd_zzflx_p_n_1}
-\mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} .  
-\end{equation} 
-
-The flux limiter iteratively adjusts the fluxes $\mathrm{flxu}$ and $\mathrm{flxv}$ until
-none of the cells will ``dry out''. To be precise the fluxes are limited until none of the
-cells has water depth less than $\mathrm{rn\_wdmin1}$ on step $n+1$.
-
-Let the fluxes on the $m$th iteration step be denoted by $\mathrm{flxu}^{(m)}$ and
-$\mathrm{flxv}^{(m)}$.  Then the adjustment is achieved by seeking a set of coefficients,
-$\mathrm{zcoef}_{i,j}^{(m)}$ such that:
-
-\begin{equation} \label{dyn_wd_continuity_coef}
-\begin{split}
-\mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\
-\mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j}
-\end{split}
-\end{equation} 
- 
-where the coefficients are $1.0$ generally but can vary between $0.0$ and $1.0$ around
-cells that would otherwise dry.
-
-The iteration is initialised by setting
-
-\begin{equation} \label{dyn_wd_zzflx_initial}
-\mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad  \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} . 
-\end{equation} 
-
-The fluxes out of cell $(i,j)$ are updated at the $m+1$th iteration if the depth of the
-cell on timestep $t_e$, namely $h_{i,j}(t_e)$, is less than the total flux out of the cell
-times the timestep divided by the cell area. Using (\ref{dyn_wd_continuity_2}) this
-condition is
-
-\begin{equation} \label{dyn_wd_continuity_if}
-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} ) .
-\end{equation} 
-
-Rearranging (\ref{dyn_wd_continuity_if}) we can obtain an expression for the maximum
-outward flux that can be allowed and still maintain the minimum wet depth:
-
-\begin{equation} \label{dyn_wd_max_flux}
-\begin{split}
-\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{]} \\
-\phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big]
-\end{split}
-\end{equation}
-
-Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\it [Q: Why is
-this necessary/desirable?]}. Substituting from (\ref{dyn_wd_continuity_coef}) gives an
-expression for the coefficient needed to multiply the outward flux at this cell in order
-to avoid drying. 
-
-\begin{equation} \label{dyn_wd_continuity_nxtcoef}
-\begin{split}
-\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{]} \\
-\phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} } 
-\end{split}
-\end{equation} 
-
-Only the outward flux components are altered but, of course, outward fluxes from one cell
-are inward fluxes to adjacent cells and the balance in these cells may need subsequent
-adjustment; hence the iterative nature of this scheme.  Note, for example, that the flux
-across the ``eastern'' face of the $(i,j)$th cell is only updated at the $m+1$th iteration
-if that flux at the $m$th iteration is out of the $(i,j)$th cell. If that is the case then
-the flux across that face is into the $(i+1,j)$ cell and that flux will not be updated by
-the calculation for the $(i+1,j)$th cell. In this sense the updates to the fluxes across
-the faces of the cells do not ``compete'' (they do not over-write each other) and one
-would expect the scheme to converge relatively quickly. The scheme is flux based so
-conserves mass. It also conserves constant tracers for the same reason that the 
-directional limiter does.  
-
-
-%----------------------------------------------------------------------------------------
-%      Surface pressure gradients
-%----------------------------------------------------------------------------------------
-\subsubsection   [Modification of surface pressure gradients (\textit{dynhpg})]
-			{Modification of surface pressure gradients (\mdl{dynhpg})}
-\label{DYN_wd_il_spg}
-
-At ``dry'' points the water depth is usually close to $\mathrm{rn\_wdmin1}$. If the
-topography is sloping at these points the sea-surface will have a similar slope and there
-will hence be very large horizontal pressure gradients at these points. The WAD modifies
-the magnitude but not the sign of the surface pressure gradients (zhpi and zhpj) at such
-points by mulitplying them by positive factors (zcpx and zcpy respectively) that lie
-between $0$ and $1$.
-
-We describe how the scheme works for the ``eastward'' pressure gradient, zhpi, calculated
-at the $(i,j)$th $u$-point. The scheme uses the ht\_wd depths and surface heights at the
-neighbouring $(i+1,j)$ and $(i,j)$ tracer points.  zcpx is calculated using two logicals
-variables, $\mathrm{ll\_tmp1}$ and $\mathrm{ll\_tmp2}$ which are evaluated for each grid
-column.  The three possible combinations are illustrated in figure \ref{Fig_WAD_dynhpg}.
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!ht] \begin{center}
-\includegraphics[width=0.8\textwidth]{Fig_WAD_dynhpg}
-\caption{ \label{Fig_WAD_dynhpg}
-Illustrations of the three possible combinations of the logical variables controlling the 
-limiting of the horizontal pressure gradient in wetting and drying regimes}
-\end{center}\end{figure}
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-
-The first logical, $\mathrm{ll\_tmp1}$, is set to true if and only if the water depth at
-both neighbouring points is greater than $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ and
-the minimum height of the sea surface at the two points is greater than the maximum height
-of the topography at the two points:
-
-\begin{equation} \label{dyn_ll_tmp1}
-\begin{split}
-\mathrm{ll\_tmp1}  = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\
-                     & \quad \mathrm{MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj))\  .and.} \\
-& \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\
-& \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\
-& \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 }
-\end{split}
-\end{equation} 
-
-The second logical, $\mathrm{ll\_tmp2}$, is set to true if and only if the maximum height
-of the sea surface at the two points is greater than the maximum height of the topography
-at the two points plus $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$
-
-\begin{equation} \label{dyn_ll_tmp2}
-\begin{split}
-\mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\
-& \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\
-& \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) .
-\end{split}
-\end{equation} 
-
-If $\mathrm{ll\_tmp1}$ is true then the surface pressure gradient, zhpi at the $(i,j)$
-point is unmodified. If both logicals are false zhpi is set to zero.
-
-If $\mathrm{ll\_tmp1}$ is true and $\mathrm{ll\_tmp2}$ is false then the surface pressure
-gradient is multiplied through by zcpx which is the absolute value of the difference in
-the water depths at the two points divided by the difference in the surface heights at the
-two points. Thus the sign of the sea surface height gradient is retained but the magnitude
-of the pressure force is determined by the difference in water depths rather than the
-difference in surface height between the two points. Note that dividing by the difference
-between the sea surface heights can be problematic if the heights approach parity. An
-additional condition is applied to $\mathrm{ ll\_tmp2 }$ to ensure it is .false. in such
-conditions.
-
-\subsection   [Additional considerations (\textit{usrdef\_zgr})]
-			{Additional considerations (\mdl{usrdef\_zgr})}
-\label{WAD_additional}
-
-In the very shallow water where wetting and drying occurs the parametrisation of 
-bottom drag is clearly very important. In order to promote stability  
-it is sometimes useful to calculate the bottom drag using an implicit time-stepping approach.  
-
-Suitable specifcation of the surface heat flux in wetting and drying domains in forced and 
-coupled simulations needs further consideration. In order to prevent freezing or boiling
-in uncoupled integrations the net surface heat fluxes need to be appropriately limited.  
  
 %----------------------------------------------------------------------------------------
 %      The WAD test cases
 %----------------------------------------------------------------------------------------
-\subsection   [The WAD test cases (\textit{usrdef\_zgr})]
+\section   [The WAD test cases (\textit{usrdef\_zgr})]
 			{The WAD test cases (\mdl{usrdef\_zgr})}
 \label{WAD_test_cases}
@@ -327,5 +58,5 @@
 
 \clearpage
-\subsubsection [WAD test case 1 : A simple linear slope]
+\subsection [WAD test case 1 : A simple linear slope]
                     {WAD test case 1 : A simple linear slope}
 \label{WAD_test_case1}
@@ -351,5 +82,5 @@
 
 \clearpage
-\subsubsection [WAD test case 2 : A parabolic channel ]
+\subsection [WAD test case 2 : A parabolic channel ]
                     {WAD test case 2 : A parabolic channel}
 \label{WAD_test_case2}
@@ -374,5 +105,5 @@
 
 \clearpage
-\subsubsection [WAD test case 3 : A parabolic channel (extreme slope) ]
+\subsection [WAD test case 3 : A parabolic channel (extreme slope) ]
                     {WAD test case 3 : A parabolic channel (extreme slope)}
 \label{WAD_test_case3}
@@ -393,5 +124,5 @@
 
 \clearpage
-\subsubsection [WAD test case 4 : A parabolic bowl ]
+\subsection [WAD test case 4 : A parabolic bowl ]
                     {WAD test case 4 : A parabolic bowl}
 \label{WAD_test_case4}
@@ -415,5 +146,5 @@
 
 \clearpage
-\subsubsection [WAD test case 5 : A double slope with shelf channel ]
+\subsection [WAD test case 5 : A double slope with shelf channel ]
                     {WAD test case 5 : A double slope with shelf channel}
 \label{WAD_test_case5}
@@ -434,5 +165,5 @@
 
 \clearpage
-\subsubsection [WAD test case 6 : A parabolic channel with central bar ]
+\subsection [WAD test case 6 : A parabolic channel with central bar ]
                     {WAD test case 6 : A parabolic channel with central bar}
 \label{WAD_test_case6}
@@ -463,5 +194,5 @@
 
 \clearpage
-\subsubsection [WAD test case 7 : A double slope with shelf, open-ended channel ]
+\subsection [WAD test case 7 : A double slope with shelf, open-ended channel ]
                     {WAD test case 7 : A double slope with shelf, open-ended channel}
 \label{WAD_test_case7}
@@ -511,6 +242,7 @@
 % ================================================================
 
-\bibliographystyle{wileyqj}
-\bibliography{references}
+%\bibliographystyle{wileyqj}
+%\bibliographystyle{../../../doc/latex/NEMO/main/ametsoc.bst}
+%\bibliography{references}
 
 \end{document}
Index: /NEMO/trunk/tests/WAD/MY_DOCS/references.bib
===================================================================
--- /NEMO/trunk/tests/WAD/MY_DOCS/references.bib	(revision 10498)
+++ /NEMO/trunk/tests/WAD/MY_DOCS/references.bib	(revision 10499)
@@ -8,16 +8,29 @@
 
 @article{Oey06,
-	author = {L-Y Oey},
-	title = {An OGCM with movable land--sea boundaries},
-	journal = {Ocean Mod.},
-	year = {2006},
-	volume = {13},
-	pages = {176--195}}
+   title = "An OGCM with movable land-sea boundaries",
+   journal = "Ocean Modelling",
+   volume = "13",
+   number = "2",
+   pages = "176 - 195",
+   year = "2006",
+   issn = "1463-5003",
+   doi = "https://doi.org/10.1016/j.ocemod.2006.01.001",
+   url = "http://www.sciencedirect.com/science/article/pii/S1463500306000084",
+   author = "Lie-Yauw Oey",
+   keywords = "Wetting and drying, Inundations, Ocean general circulation model (OGCM), Princeton Ocean Model (POM), Tides, Tsunamis, Estuarine outflows",
+   abstract = "An ocean general circulation model (OGCM) with wetting and drying (WAD) capabilities removes the vertical-wall coastal assumption and allows simultaneous modeling of open-ocean currents and water run-up (and run-down) across movable land-sea boundaries. This paper implements and tests such a WAD scheme for the Princeton Ocean Model (POM) in its most general three-dimensional setting with stratification, bathymetry and forcing. The scheme can be easily exported to other OGCM's."
+}
 
 @article{WarnerEtal13,
-	author = {J C Warner and Z Defne and K Haas and H G Arango},
-	title = {{A wetting and drying scheme for ROMS}},
-	journal = {Computers and Geosc.},
-	year = {2013},
-	volume = {58},
-	pages = {54--61}}
+   title = "A wetting and drying scheme for ROMS",
+   journal = "Computers \& Geosciences",
+   volume = "58",
+   pages = "54 - 61",
+   year = "2013",
+   issn = "0098-3004",
+   doi = "https://doi.org/10.1016/j.cageo.2013.05.004",
+   url = "http://www.sciencedirect.com/science/article/pii/S0098300413001362",
+   author = "John C. Warner and Zafer Defne and Kevin Haas and Hernan G. Arango",
+   keywords = "Wetting and drying, ROMS, Cell-face blocking",
+   abstract = "The processes of wetting and drying have many important physical and biological impacts on shallow water systems. Inundation and dewatering effects on coastal mud flats and beaches occur on various time scales ranging from storm surge, periodic rise and fall of the tide, to infragravity wave motions. To correctly simulate these physical processes with a numerical model requires the capability of the computational cells to become inundated and dewatered. In this paper, we describe a method for wetting and drying based on an approach consistent with a cell-face blocking algorithm. The method allows water to always flow into any cell, but prevents outflow from a cell when the total depth in that cell is less than a user defined critical value. We describe the method, the implementation into the three-dimensional Regional Oceanographic Modeling System (ROMS), and exhibit the new capability under three scenarios: an analytical expression for shallow water flows, a dam break test case, and a realistic application to part of a wetland area along the Georgia Coast, USA."
+}
