Index: branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg	(revision 8868)
+++ branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg	(revision 8985)
@@ -444,6 +444,10 @@
 &namwad  !   Wetting and drying
 !-----------------------------------------------------------------------
-   ln_wd             = .true.   ! T/F activation of wetting and drying
-   rn_wdmin1         =  0.3     ! Minimum wet depth on dried cells
+   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
Index: branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_DOCS/Namelists/nam_wad
===================================================================
--- branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_DOCS/Namelists/nam_wad	(revision 8868)
+++ branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_DOCS/Namelists/nam_wad	(revision 8985)
@@ -2,8 +2,12 @@
 &namwad  !   Wetting and drying
 !-----------------------------------------------------------------------
-   ln_wd             = .true.   ! T/F activation of wetting and drying
-   rn_wdmin1         =  0.3     ! Minimum wet depth on dried cells
-   rn_wdmin2         =   0.0001 ! Tolerance of min wet depth on dried cells
+   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
+   nn_wdit           =   20     ! Max iterations for W/D limiter
 /
Index: branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_DOCS/WAD_doc.tex
===================================================================
--- branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_DOCS/WAD_doc.tex	(revision 8868)
+++ branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_DOCS/WAD_doc.tex	(revision 8985)
@@ -3,7 +3,7 @@
 \begin{document}
 
-\title{Draft description of NEMO wetting and drying scheme:     22 November 2016 }
-
-\author{ Hedong Liu, Jason Holt, Andrew Coward  and Michael J. Bell  }
+\title{Draft description of NEMO wetting and drying scheme:     29 November 2017 }
+
+\author{ Enda O'Dea, Hedong Liu, Jason Holt, Andrew Coward  and Michael J. Bell  }
 
 %------------------------------------------------------------------------
@@ -28,59 +28,96 @@
 \label{DYN_wetdry}
 
-This is preliminary documentation for the wetting and drying code (WAD).  The emphasis is
-on explaining the rationale for the code.  The approach used by the WAD is similar to that
-developed for POM by \cite{Oey06} and that developed for ROMS by \cite{WarnerEtal13} but
-the WAD uses schemes that have not been published.
+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) will be denoted by $ \mathrm{ssh}$. Both
-quantities are measured relative to a reference sea level at z$=$0m. Given the sign
-conventions used, the water depth is the height of the free surface plus the depth of the
+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}$).
 
-\namdisplay{nam_wad}
-
-WAD is activated by setting $\mathrm{ln\_wd} = \mathrm{.true.}$. Currently, this option
-works with six test cases provided in the WAD\_TEST\_CASES configuration. These are all
-pure sigma coordinate configurations which define their domain, surface forcing and
-initial conditions via a set of 'usrdef' routines in MY\_SRC. Extending this option to
-more realistic domains will require the derivation and provision of a suitable
-$\mathrm{ht\_wd}$ field in addition to the normal information provided in the domcfg.nc
-file. The six test cases are described in section \S\ref{WAD_test_cases}.
-
-The WAD takes all points in the domain below a land elevation of $\mathrm{rn\_wdld}$ to be
-covered by water. Points where the water depth is less than $\mathrm{rn\_wdmin1}$ are to
-be interpreted as ``dry''. The WAD requires the topography specified with a model
+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 computed relative to an
-initial state with zero sea surface height elevation. These reference metrics and depths
-(i.e. the $\mathrm{e3t\_0, ht\_0}$ etc. arrays) are unaltered by WAD.
-$\mathrm{rn\_wdmin1}$ is usually chosen to be of order $0.075$m but complex topographies
-with steep slopes may require larger values. The scheme also makes use of a second
-parameter, $\mathrm{rn\_wdmin2}$, which is intended to be much smaller than
-$\mathrm{rn\_wdmin1}$, of order $10^{-6}$m or smaller {\it (Q: What is the purpose of
-$\mathrm{rn\_wdmin2}$? Seems a non-zero value is required for the flux limiter iterations
-to converge)}.
-
-The WAD 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
-first sub-section below describes this scheme. It also briefly describes the simpler ROMS
-method that has not been implemented.
-
-The following sub-section describes how the surface pressure gradients are modified by the
-WAD. The next sub-section should describe how the WAD maintains consistency between the
-points that are ``wet'' on the barotropic sub-steps and those that are wet on the longer
-baroclinic time-step. This sub-section has not yet been written. The final sub-section
-should describe the test cases that have been used to assess the performance of the WAD.
+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. 
 
 %-----------------------------------------------------------------------------------------
-%   Flux limiters
+%   Iterative limiters
 %-----------------------------------------------------------------------------------------
-\subsection   [Flux limiters (\textit{wet\_dry})]
-			{Flux limiters (\mdl{wet\_dry})}
-\label{DYN_wd_flux_limit}
-
-The flux limiter for the barotropic flow devised by Hedong Liu can be understood as follows: 
+\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 
@@ -173,22 +210,15 @@
 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 also flux based so
-conserves mass.
-
-The ROMS scheme to prevent drying out of a cell is somewhat simpler. It specifies that if
-a tracer cell is dry (the water depth is less than $\mathrm{rn\_wdmin1}$) on the backward
-timestep, $t_e$, then any outward flux through its cell faces should be set to zero. This
-scheme has a clear physical rationale. This scheme is equivalent to setting
-$\mathrm{zcoef}^{(m+1)}_{i,j}$ to $0.0$ whenever a cell is at risk of drying.  One
-objection to the ROMS scheme is that it introduces a spurious step function in the flux
-out of a cell as the water depth in the cell passes through the ``critical'' value
-$\mathrm{rn\_wdmin1}$.
+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
 %----------------------------------------------------------------------------------------
-\subsection   [Modification of surface pressure gradients (\textit{dynhpg})]
+\subsubsection   [Modification of surface pressure gradients (\textit{dynhpg})]
 			{Modification of surface pressure gradients (\mdl{dynhpg})}
-\label{DYN_wd_spg}
+\label{DYN_wd_il_spg}
 
 At ``dry'' points the water depth is usually close to $\mathrm{rn\_wdmin1}$. If the
@@ -253,4 +283,16 @@
 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
Index: branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_DOCS/references.bib
===================================================================
--- branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_DOCS/references.bib	(revision 8985)
+++ branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_DOCS/references.bib	(revision 8985)
@@ -0,0 +1,23 @@
+@article{Oey05,
+	author = {L-Y Oey},
+	title = {{A wetting and drying scheme for POM}},
+	journal = {Ocean Mod.},
+	year = {2005},
+	volume = {9},
+	pages = {133--150}}
+
+@article{Oey06,
+	author = {L-Y Oey},
+	title = {An OGCM with movable land--sea boundaries},
+	journal = {Ocean Mod.},
+	year = {2006},
+	volume = {13},
+	pages = {176--195}}
+
+@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}}
Index: branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/domain.F90
===================================================================
--- branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/domain.F90	(revision 8868)
+++ branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/domain.F90	(revision 8985)
@@ -674,7 +674,6 @@
       ENDIF
       !
-      IF( ln_wd ) THEN              ! wetting and drying domain
+      IF( ll_wd ) THEN              ! wetting and drying domain
          CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 )
-         CALL iom_rstput( 0, 0, inum, 'ht_wd'  , ht_wd  , ktype = jp_r8 )
       ENDIF
       !
Index: branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/usrdef_istate.F90
===================================================================
--- branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/usrdef_istate.F90	(revision 8868)
+++ branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/usrdef_istate.F90	(revision 8985)
@@ -169,4 +169,6 @@
             !
       END SELECT
+
+
       !
       ! Apply minimum wetdepth criterion
@@ -174,6 +176,6 @@
       do jj = 1,jpj
          do ji = 1,jpi
-            IF( ht_wd(ji,jj) + pssh(ji,jj) < rn_wdmin1 ) THEN
-               pssh(ji,jj) = ptmask(ji,jj,1)*( rn_wdmin1 - ht_wd(ji,jj) )
+            IF( ht_0(ji,jj) + pssh(ji,jj) < rn_wdmin1 ) THEN
+               pssh(ji,jj) = ptmask(ji,jj,1)*( rn_wdmin1 - ht_0(ji,jj) )
             ENDIF
          end do
Index: branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/usrdef_zgr.F90
===================================================================
--- branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/usrdef_zgr.F90	(revision 8868)
+++ branches/2017/dev_METO_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/usrdef_zgr.F90	(revision 8985)
@@ -16,8 +16,8 @@
    USE oce            ! ocean variables
    USE dom_oce ,  ONLY: ln_zco, ln_zps, ln_sco   ! ocean space and time domain
-   USE dom_oce ,  ONLY: mi0, mi1, nimpp, njmpp,  &
-                      & mj0, mj1, glamt, gphit   ! ocean space and time domain
+   USE dom_oce ,  ONLY: ht_0, mi0, mi1, nimpp, njmpp,  &
+                      & mj0, mj1, glamt, gphit         ! ocean space and time domain
    USE usrdef_nam     ! User defined : namelist variables
-   USE wet_dry ,  ONLY: rn_wdmin1, rn_wdmin2, rn_wdld, ht_wd
+   USE wet_dry ,  ONLY: rn_wdmin1, rn_wdmin2, rn_wdld  ! Wetting and drying
    !
    USE in_out_manager ! I/O manager
@@ -232,6 +232,6 @@
          END SELECT
       END IF
-               !
-      !
+              
+
       ! at u-point: averaging zht
       DO ji = 1, jpim1
@@ -283,5 +283,5 @@
       IF ( ln_sco ) THEN      !==  s-coordinate  ==!   (terrain-following coordinate)
          !
-         ht_wd = zht
+         ht_0 = zht
          k_bot(:,:) = jpkm1 * k_top(:,:)  !* bottom ocean = jpk-1 (here use k_top as a land mask)
          DO jj = 1, jpj
