Index: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/CONFIG/AMM12/cpp_AMM12.fcm
===================================================================
--- branches/2015/dev_r5803_NOC_WAD/NEMOGCM/CONFIG/AMM12/cpp_AMM12.fcm	(revision 5839)
+++ branches/2015/dev_r5803_NOC_WAD/NEMOGCM/CONFIG/AMM12/cpp_AMM12.fcm	(revision 5842)
@@ -1,1 +1,1 @@
- bld::tool::fppkeys  key_bdy key_tide key_dynspg_ts key_ldfslp  key_zdfgls  key_vvl key_diainstant key_mpp_mpi key_iomput
+ bld::tool::fppkeys  key_bdy key_tide key_dynspg_ts key_ldfslp  key_zdfgls  key_vvl key_diainstant key_mpp_mpi 
Index: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/CONFIG/cfg.txt
===================================================================
--- branches/2015/dev_r5803_NOC_WAD/NEMOGCM/CONFIG/cfg.txt	(revision 5839)
+++ branches/2015/dev_r5803_NOC_WAD/NEMOGCM/CONFIG/cfg.txt	(revision 5842)
@@ -5,5 +5,4 @@
 C1D_PAPA OPA_SRC
 GYRE_BFM OPA_SRC TOP_SRC
-AMM12 OPA_SRC
 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC
 GYRE OPA_SRC
@@ -11,2 +10,4 @@
 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC
 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC
+AMM12 OPA_SRC
+IRS18WD OPA_SRC
Index: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 5839)
+++ branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 5842)
@@ -10,4 +10,5 @@
    !!            3.5  ! 2012     (S. Mocavero, I. Epicoco) Add arrays associated
    !!                             to the optimization of BDY communications
+   !!            3.6.?! 2014     (H. Liu) Add arrays associated with Wetting and Drying
    !!----------------------------------------------------------------------
 
@@ -270,4 +271,23 @@
    INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: nicoa, njcoa       !: ???
 #endif
+
+   !!----------------------------------------------------------------------
+   !! critical depths,filters, limiters,and masks for  Wetting and Drying
+   !! ---------------------------------------------------------------------
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wduflt, wdvflt     !: u- and v- filter
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wdmask             !: u- and v- limiter 
+
+   LOGICAL,  PUBLIC, SAVE ::   ln_wd       !: key to turn on/off wetting/drying (T: on, F: off)
+   REAL(wp), PUBLIC, SAVE ::   rn_wdmin1   !: minimum water depth on dried cells
+   REAL(wp), PUBLIC, SAVE ::   rn_wdmin2   !: tolerrance of minimum water depth on dried cells
+   REAL(wp), PUBLIC, SAVE ::   rn_wdld     !: land elevation below which wetting/drying will be considered
+   INTEGER , PUBLIC, SAVE ::   nn_wdit     !: maximum number of iteration for W/D limiter
+
+   !LOGICAL,  PUBLIC, SAVE ::   ln_wd     =  .FALSE.  !: turn on wetting/drying (T: on, F: off)
+   !REAL(wp), PUBLIC, SAVE ::   rn_wdmin1 = 0.10_wp   !: minimum water depth on dried cells
+   !REAL(wp), PUBLIC, SAVE ::   rn_wdmin2 = 0.01_wp   !: tolerrance of minimum water depth on dried cells
+   !REAL(wp), PUBLIC, SAVE ::   rn_wdld   = 20.0_wp   !: land elevation below which wetting/drying will be considered
+   !INTEGER , PUBLIC, SAVE ::   nn_wdit   =   10      !: maximum number of iteration for W/D limiter
 
    !!----------------------------------------------------------------------
@@ -407,4 +427,7 @@
       ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(12) )
 #endif
+
+      IF(ln_wd) &
+      ALLOCATE( wduflt(jpi,jpj), wdvflt(jpi,jpj), wdmask(jpi,jpj), STAT=ierr(12) )
       !
       dom_oce_alloc = MAXVAL(ierr)
Index: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 5839)
+++ branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 5842)
@@ -147,4 +147,16 @@
       fse3t_a(:,:,jpk) = e3t_0(:,:,jpk)
 
+      !IF(ln_wd) THEN
+      !  DO jj = 1, jpj
+      !    DO ji = 1, jpi
+      !      IF(mbathy(ji,jj) == 2 .AND. e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN
+      !       fse3t_a(ji,jj,1:2) = 0.5_wp * rn_wdmin1 
+      !       fse3t_n(ji,jj,1:2) = 0.5_wp * rn_wdmin1 
+      !       fse3t_b(ji,jj,1:2) = 0.5_wp * rn_wdmin1 
+      !      END IF
+      !    ENDDO
+      !  ENDDO
+      !END IF
+
       ! Reconstruction of all vertical scale factors at now and before time steps
       ! =============================================================================
@@ -701,4 +713,7 @@
       !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
       !! * Local declarations
+      REAL(wp) ::  zwad                                                ! = 1.0 when ln_wd = .true.
+                                                                       ! = 0.0 when ln_wd = .false.
+                                                                       ! 
       INTEGER ::   ji, jj, jk                                          ! dummy loop indices
       LOGICAL ::   l_is_orca                                           ! local logical
@@ -706,4 +721,10 @@
       IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol')
          !
+      IF(ln_wd) THEN
+        zwad = 1.0_wp
+      ELSE
+        zwad = 0.0_wp
+      END IF
+
       l_is_orca = .FALSE.
       IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE.      ! ORCA R2 configuration - will need to correct some locations
@@ -717,5 +738,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   &
+                  !pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   &
+                  pe3_out(ji,jj,jk) = 0.5_wp * (umask(ji,jj,jk) * (1.0_wp - zwad) + zwad) * r1_e12u(ji,jj)        &
                      &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
                      &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
@@ -735,5 +757,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   &
+                  !pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   &
+                  pe3_out(ji,jj,jk) = 0.5_wp * (vmask(ji,jj,jk) * (1.0_wp - zwad) + zwad) * r1_e12v(ji,jj)        &
                      &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
                      &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
@@ -753,5 +776,7 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               &
+                  !pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               &
+                  pe3_out(ji,jj,jk) = 0.5_wp * (umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zwad) + zwad)     &
+                     &                       * r1_e12f(ji,jj)                                                     &
                      &                       * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
                      &                           + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
@@ -771,6 +796,10 @@
          ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
          DO jk = 2, jpk
-            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   &
-               &                            +            0.5_wp * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
+            !pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   &
+            !   &                            +            0.5_wp * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
+            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * (tmask(:,:,jk) * (1.0_wp - zwad) + zwad) ) &
+               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                         &
+               &                            +            0.5_wp * (tmask(:,:,jk) * (1.0_wp - zwad) + zwad)   &
+               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
          END DO
          !               ! -------------------------------------- !
@@ -781,6 +810,10 @@
          ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
          DO jk = 2, jpk
-            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   &
-               &                             +            0.5_wp * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
+            !pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   &
+            !   &                             +            0.5_wp * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
+            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * (umask(:,:,jk) * (1.0_wp - zwad) + zwad) ) &
+               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                         &
+               &                             +            0.5_wp * (umask(:,:,jk) * (1.0_wp - zwad) + zwad)   &
+               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
          END DO
          !               ! -------------------------------------- !
@@ -791,6 +824,10 @@
          ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
          DO jk = 2, jpk
-            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   &
-               &                             +            0.5_wp * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
+            !pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   &
+            !   &                             +            0.5_wp * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
+            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * (vmask(:,:,jk) * (1.0_wp - zwad) + zwad) ) &
+               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                         &
+               &                             +            0.5_wp * (vmask(:,:,jk) * (1.0_wp - zwad) + zwad)   &
+               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
          END DO
       END SELECT
@@ -817,5 +854,5 @@
       CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
       !! * Local declarations
-      INTEGER ::   jk
+      INTEGER ::   ji, jj, jk
       INTEGER ::   id1, id2, id3, id4, id5     ! local integers
       !!----------------------------------------------------------------------
@@ -905,4 +942,22 @@
             fse3t_n(:,:,:) = e3t_0(:,:,:)
             sshn(:,:) = 0.0_wp
+
+            IF(ln_wd) THEN
+              DO jj = 1, jpj
+                DO ji = 1, jpi
+                  !IF(e3t_0(ji,jj,1) < 0._wp) THEN
+                  !IF(mbathy(ji,jj) == 2 .AND. e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN
+                  IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN
+                    fse3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1 
+                    fse3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1 
+                    fse3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1 
+                    sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj)
+                    sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj)
+                    ssha(ji,jj) = rn_wdmin1 - bathy(ji,jj)
+                  ENDIF
+                ENDDO
+              ENDDO
+            END IF
+
             IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
                tilde_e3t_b(:,:,:) = 0.0_wp
Index: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5839)
+++ branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5842)
@@ -18,4 +18,5 @@
    !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case  
    !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capabilitye  
+   !!            3.?  ! 2015-11  (H. Liu) Modifications for Wetting/Drying
    !!----------------------------------------------------------------------
 
@@ -1810,4 +1811,5 @@
       REAL(wp) ::   zrmax, ztaper   ! temporary scalars
       REAL(wp) ::   zrfact
+      REAL(wp) ::   zsmth
       !
       REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
@@ -1865,4 +1867,5 @@
       bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
 
+      IF(.NOT.ln_wd) THEN                      
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -1870,4 +1873,5 @@
          END DO
       END DO
+      END IF
       !                                        ! =============================
       !                                        ! Define the envelop bathymetry   (hbatt)
@@ -1876,4 +1880,5 @@
       zenv(:,:) = bathy(:,:)
       !
+      IF(.NOT.ln_wd) THEN    
       ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
       DO jj = 1, jpj
@@ -1891,4 +1896,6 @@
          END DO
       END DO
+      END IF
+
       ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
       CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
@@ -1930,12 +1937,19 @@
          zri(:,:) = 0._wp
          zrj(:,:) = 0._wp
+         
+         !IF(ln_wd) THEN                     !extend the smoothed region to cover the W/D zones
+         !  zsmth = -rn_wdld    
+         !ELSE
+           zsmth = 0._wp                     ! The original form (only smooth ocean points)
+         !ENDIF
+
          DO jj = 1, nlcj
             DO ji = 1, nlci
                iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
                ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
-               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
+               IF( (zenv(ji,jj) > zsmth) .AND. (zenv(iip1,jj) > zsmth)) THEN
                   zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
                END IF
-               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
+               IF( (zenv(ji,jj) > zsmth) .AND. (zenv(ji,ijp1) > zsmth)) THEN
                   zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
                END IF
@@ -1960,9 +1974,17 @@
       END DO                                                !     End loop     !
       !                                                     ! ================ !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
-         END DO
-      END DO
+      IF(ln_wd) THEN
+        DO jj = 1, jpj
+           DO ji = 1, jpi
+              !zenv(ji,jj) = MAX( zenv(ji,jj), -rn_wdld )    ! filt out land bathy data 
+           END DO
+        END DO
+      ELSE
+        DO jj = 1, jpj
+           DO ji = 1, jpi
+              zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
+           END DO
+        END DO
+      ENDIF
       !
       ! Envelope bathymetry saved in hbatt
@@ -1994,5 +2016,10 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
+         IF(.NOT.ln_wd) THEN
+           WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
+         ELSE
+           WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min
+           WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld
+         ENDIF
       ENDIF
       hbatu(:,:) = rn_sbot_min
@@ -2007,4 +2034,19 @@
         END DO
       END DO
+
+      IF(ln_wd) THEN               !avoid the zero depth on T- (U-,V-,F-) points
+        DO jj = 1, jpj
+          DO ji = 1, jpi
+            IF(ABS(hbatt(ji,jj)) < rn_wdmin1) &
+              & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1
+            IF(ABS(hbatu(ji,jj)) < rn_wdmin1) &
+              & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1
+            IF(ABS(hbatv(ji,jj)) < rn_wdmin1) &
+              & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1
+            IF(ABS(hbatf(ji,jj)) < rn_wdmin1) &
+              & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1
+          END DO
+        END DO
+      END IF
       ! 
       ! Apply lateral boundary condition
@@ -2014,4 +2056,5 @@
          DO ji = 1, jpi
             IF( hbatu(ji,jj) == 0._wp ) THEN
+               !No worries about the following line when ln_wd == .true.
                IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
                IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
@@ -2039,8 +2082,10 @@
 
 !!bug:  key_helsinki a verifer
-      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
-      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
-      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
-      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
+      IF(.NOT.ln_wd) THEN
+       hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
+       hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
+       hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
+       hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
+      END IF
 
       IF( nprint == 1 .AND. lwp )   THEN
@@ -2087,11 +2132,13 @@
       fsde3w(:,:,:) = gdep3w_0(:,:,:)
       !
-      where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
-      where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
-      where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
-      where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
-      where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
-      where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
-      where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
+      IF(.NOT.ln_wd) THEN
+        where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
+        where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
+        where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
+        where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
+        where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
+        where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
+        where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
+      END IF
 
 #if defined key_agrif
@@ -2135,7 +2182,17 @@
                IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
             END DO
-            IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
-         END DO
-      END DO
+
+            IF(ln_wd) THEN
+              IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN
+                mbathy(ji,jj) = 0
+              ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN
+                mbathy(ji,jj) = 2
+              ENDIF
+            ELSE
+              IF( scobot(ji,jj) == 0._wp )   mbathy(ji,jj) = 0
+            ENDIF
+         END DO
+      END DO
+
       IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
          &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
@@ -2257,4 +2314,6 @@
       INTEGER  ::   ji, jj, jk           ! dummy loop argument
       REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
+      REAL(wp) ::   ztmpu,  ztmpv,  ztmpf
+      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
       !
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
@@ -2310,16 +2369,54 @@
       DO ji = 1, jpim1
          DO jj = 1, jpjm1
+            ! extended for Wetting/Drying case
+            ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj)
+            ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
+            ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
+            ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
+            ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
+            ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
+                   & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
             DO jk = 1, jpk
-               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
-               z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
-                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
-               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+               IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN
+                 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
+                 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
+               ELSE
+                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
+                        &              / ztmpu
+                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
+                        &              / ztmpu
+               END IF
+
+               IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN
+                 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
+                 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
+               ELSE
+                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
+                        &              / ztmpv
+                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
+                        &              / ztmpv
+               END IF
+
+               IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN
+                 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)  &
+                    &                            + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) )
+               ELSE
+                 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)  &
+                    &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) &
+                    &              / ztmpf
+               END IF
+
+               !
+               !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
+               !   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+               !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
+               !   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+               !z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
+               !   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
+               !   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
+               !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
+               !   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+               !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
+               !   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
                !
                e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
@@ -2361,4 +2458,6 @@
       REAL(wp) ::   zsmth               ! smoothing around critical depth
       REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
+      REAL(wp) ::   ztmpu, ztmpv, ztmpf
+      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
       !
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
@@ -2439,16 +2538,53 @@
         DO jj=1,jpj-1
 
-          DO jk = 1, jpk
-                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
-                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  &
-                                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
-                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+           ! extend to suit for Wetting/Drying case
+           ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj)
+           ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
+           ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
+           ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
+           ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
+           ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
+                  & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
+           DO jk = 1, jpk
+              IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN
+                z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
+                z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
+              ELSE
+                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
+                       &              / ztmpu
+                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
+                       &              / ztmpu
+              END IF
+
+              IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN
+                z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
+                z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
+              ELSE
+                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
+                       &              / ztmpv
+                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
+                       &              / ztmpv
+              END IF
+
+              IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN
+                z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)  &
+                   &                            + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) )
+              ELSE
+                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)  &
+                   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) &
+                   &              / ztmpf
+              END IF
+
+             !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
+             !                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+             !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
+             !                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+             !z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  &
+             !                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
+             !                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
+             !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
+             !                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+             !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
+             !                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
 
              e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
Index: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 5839)
+++ branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 5842)
@@ -17,4 +17,5 @@
    !!                 !           (A. Coward) suppression of hel, wdj and rot options
    !!            3.6  !  2014-11  (P. Mathiot) hpg_isf: original code for ice shelf cavity
+   !!            3.6? !  2015-11  (H. Liu) add Wetting/Drying pressure filter
    !!----------------------------------------------------------------------
 
@@ -378,10 +379,13 @@
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       !!
-      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
-      REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars
+      INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices
+      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars
+      LOGICAL  ::   ll_tmp1, ll_tmp2, ll_tmp3            ! local logical variables
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy !W/D pressure filter
       !!----------------------------------------------------------------------
       !
       CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )
+      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
       !
       IF( kt == nit000 ) THEN
@@ -397,4 +401,50 @@
       ELSE                 ;     znad = 0._wp         ! Fixed volume
       ENDIF
+
+      IF(ln_wd) THEN
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1 
+             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
+             ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1 + rn_wdmin2
+             ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + &
+                                                       & rn_wdmin1 + rn_wdmin2
+
+             IF(ll_tmp1.AND.ll_tmp2) THEN
+               zcpx(ji,jj) = 1.0_wp
+               wduflt(ji,jj) = 1.0_wp
+             ELSE IF(ll_tmp3) THEN
+               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here
+               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) / &
+                           &     (sshn(ji+1,jj) - sshn(ji,jj)))
+               wduflt(ji,jj) = 1.0_wp
+             ELSE
+               zcpx(ji,jj) = 0._wp
+               wduflt(ji,jj) = 0.0_wp
+             END IF
+      
+             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
+             ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1 + rn_wdmin2
+             ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + &
+                                                       & rn_wdmin1 + rn_wdmin2
+
+             IF(ll_tmp1.AND.ll_tmp2) THEN
+               zcpy(ji,jj) = 1.0_wp
+               wdvflt(ji,jj) = 1.0_wp
+             ELSE IF(ll_tmp3) THEN
+               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here
+               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) / &
+                           &     (sshn(ji,jj+1) - sshn(ji,jj)))
+               wdvflt(ji,jj) = 1.0_wp
+             ELSE
+               zcpy(ji,jj) = 0._wp
+               wdvflt(ji,jj) = 0.0_wp
+             END IF
+           END DO
+        END DO
+        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
+      ENDIF
+
+
+      !jii=jidbg-nimpp+1;jjj=jjdbg-njmpp+1
 
       ! Surface value
@@ -411,4 +461,14 @@
             zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   &
                &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
+
+
+            IF(ln_wd) THEN
+
+              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
+              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
+              zuap = zuap * zcpx(ji,jj)
+              zvap = zvap * zcpy(ji,jj)
+            ENDIF
+
             ! add to the general momentum trend
             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
@@ -433,4 +493,12 @@
                zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
                   &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj)
+
+               IF(ln_wd) THEN
+                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
+                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
+                 zuap = zuap * zcpx(ji,jj)
+                 zvap = zvap * zcpy(ji,jj)
+               ENDIF
+
                ! add to the general momentum trend
                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
@@ -441,4 +509,5 @@
       !
       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
+      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
       !
    END SUBROUTINE hpg_sco
@@ -719,8 +788,10 @@
       REAL(wp) ::   z1_10, cffu, cffx   !    "         "
       REAL(wp) ::   z1_12, cffv, cffy   !    "         "
+      LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  dzx, dzy, dzz, dzu, dzv, dzw
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  drhox, drhoy, drhoz, drhou, drhov, drhow
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  rho_i, rho_j, rho_k
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter
       !!----------------------------------------------------------------------
       !
@@ -728,5 +799,48 @@
       CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow )
       CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        )
-      !
+      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
+      !
+      !!----------------------------------------------------------------------
+      !
+      !
+      IF(ln_wd) THEN
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1 
+             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) &
+                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) &
+                     &  > rn_wdmin1 + rn_wdmin2
+             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +&
+                                                       & rn_wdmin1 + rn_wdmin2
+
+             IF(ll_tmp1) THEN
+               zcpx(ji,jj) = 1.0_wp
+             ELSE IF(ll_tmp2) THEN
+               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here
+               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /&
+                           &     (sshn(ji+1,jj) - sshn(ji,jj)))
+             ELSE
+               zcpx(ji,jj) = 0._wp
+             END IF
+      
+             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) &
+                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) &
+                     &  > rn_wdmin1 + rn_wdmin2
+             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +&
+                                                       & rn_wdmin1 + rn_wdmin2
+
+             IF(ll_tmp1) THEN
+               zcpy(ji,jj) = 1.0_wp
+             ELSE IF(ll_tmp2) THEN
+               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here
+               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /&
+                           &     (sshn(ji,jj+1) - sshn(ji,jj)))
+             ELSE
+               zcpy(ji,jj) = 0._wp
+             END IF
+           END DO
+        END DO
+        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
+      ENDIF
+
 
       IF( kt == nit000 ) THEN
@@ -899,4 +1013,8 @@
             zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj)
             zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj)
+            IF(ln_wd) THEN
+              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
+              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
+            ENDIF
             ! add to the general momentum trend
             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
@@ -918,4 +1036,8 @@
                   &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    &
                   &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) / e2v(ji,jj)
+               IF(ln_wd) THEN
+                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
+                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
+               ENDIF
                ! add to the general momentum trend
                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
@@ -928,4 +1050,5 @@
       CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow )
       CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        )
+      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
       !
    END SUBROUTINE hpg_djc
@@ -951,4 +1074,5 @@
       !! The local variables for the correction term
       INTEGER  :: jk1, jis, jid, jjs, jjd
+      LOGICAL  :: ll_tmp1, ll_tmp2                  ! local logical variables
       REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps
       REAL(wp) :: zrhdt1
@@ -957,4 +1081,5 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
       REAL(wp), POINTER, DIMENSION(:,:)   ::   zsshu_n, zsshv_n
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter
       !!----------------------------------------------------------------------
       !
@@ -962,4 +1087,5 @@
       CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh )
       CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n )
+      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
       !
       IF( kt == nit000 ) THEN
@@ -974,4 +1100,43 @@
       znad = 0.0_wp
       IF( lk_vvl ) znad = 1._wp
+
+      IF(ln_wd) THEN
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1 
+             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) &
+                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) &
+                     &  > rn_wdmin1 + rn_wdmin2
+             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +&
+                                                       & rn_wdmin1 + rn_wdmin2
+
+             IF(ll_tmp1) THEN
+               zcpx(ji,jj) = 1.0_wp
+             ELSE IF(ll_tmp2) THEN
+               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here
+               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /&
+                           &     (sshn(ji+1,jj) - sshn(ji,jj)))
+             ELSE
+               zcpx(ji,jj) = 0._wp
+             END IF
+      
+             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) &
+                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) &
+                     &  > rn_wdmin1 + rn_wdmin2
+             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +&
+                                                       & rn_wdmin1 + rn_wdmin2
+
+             IF(ll_tmp1.OR.ll_tmp2) THEN
+               zcpy(ji,jj) = 1.0_wp
+             ELSE IF(ll_tmp2) THEN
+               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here
+               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /&
+                           &     (sshn(ji,jj+1) - sshn(ji,jj)))
+             ELSE
+               zcpy(ji,jj) = 0._wp
+             END IF
+           END DO
+        END DO
+        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
+      ENDIF
 
       ! Clean 3-D work arrays
@@ -1052,4 +1217,7 @@
         END DO
       END DO
+
+      CALL lbc_lnk (zsshu_n, 'U', 1.)
+      CALL lbc_lnk (zsshv_n, 'V', 1.)
 
       DO jj = 2, jpjm1
@@ -1150,6 +1318,9 @@
                ENDIF
 
-               ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * &
-               &           umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk)
+               IF(ln_wd) THEN
+                  zdpdx1 = zdpdx1 * zcpx(ji,jj)
+                  zdpdx2 = zdpdx2 * zcpx(ji,jj)
+                ENDIF
+                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 
             ENDIF
 
@@ -1207,6 +1378,10 @@
                ENDIF
 
-               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*&
-               &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk)
+               IF(ln_wd) THEN
+                  zdpdy1 = zdpdy1 * zcpy(ji,jj)
+                  zdpdy2 = zdpdy2 * zcpy(ji,jj)
+                ENDIF
+
+               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk)
             ENDIF
 
@@ -1219,4 +1394,5 @@
       CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh )
       CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n )
+      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
       !
    END SUBROUTINE hpg_prj
Index: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 5839)
+++ branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 5842)
@@ -41,4 +41,5 @@
    USE timing          ! Timing    
    USE sbcapr          ! surface boundary condition: atmospheric pressure
+   USE wadlmt          ! wetting/drying flux limter
    USE dynadv, ONLY: ln_dynadv_vec
 #if defined key_agrif
@@ -142,4 +143,5 @@
       LOGICAL  ::   ll_fw_start			! if true, forward integration 
       LOGICAL  ::   ll_init			    ! if true, special startup of 2d equations
+      LOGICAL  ::   ll_tmp1, ll_tmp2            ! local logical variables used in W/D
       INTEGER  ::   ji, jj, jk, jn   		! dummy loop indices
       INTEGER  ::   ikbu, ikbv, noffset    	! local integers
@@ -157,5 +159,11 @@
       REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a
       REAL(wp), POINTER, DIMENSION(:,:) :: zhf
+      REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef.
+      REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1           ! Wetting/Dying velocity filter coef.
       !!----------------------------------------------------------------------
+
+      ! tempoary debugging integer
+      INTEGER :: jidbg, jjdbg
+      jidbg = 163; jjdbg = 179
       !
       IF( nn_timing == 1 )  CALL timing_start('dyn_spg_ts')
@@ -168,4 +176,5 @@
       CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a                                   )
       CALL wrk_alloc( jpi, jpj, zhf )
+      IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 )
       !
       !                                         !* Local constant initialization
@@ -379,12 +388,92 @@
       !                                   ! ----------------------------------------------------
       IF( lk_vvl ) THEN                         ! Variable volume : remove surface pressure gradient
-         DO jj = 2, jpjm1 
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj)
-               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj)
-            END DO
-         END DO
-      ENDIF
-
+        IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters
+          DO jj = 2, jpjm1
+             DO ji = 2, jpim1
+                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) &
+                        & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) &
+                        &  > rn_wdmin1 + rn_wdmin2
+                ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + &
+                                                          & rn_wdmin1 + rn_wdmin2
+                IF(ll_tmp1) THEN
+                  zcpx(ji,jj) = 1.0_wp
+                  wduflt1(ji,jj) = 1.0_wp
+                ELSE IF(ll_tmp2) THEN
+                   ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here
+                  zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /&
+                              &     (sshn(ji+1,jj) - sshn(ji,jj)))
+                  wduflt1(ji,jj) = 1.0_wp
+                ELSE
+                  zcpx(ji,jj) = 0._wp
+                  wduflt1(ji,jj) = 0.0_wp
+                END IF
+
+                !for w/d debugging, delete it when finished
+                !   zcpx(ji,jj) = 0._wp
+                ! end debugging part
+                   
+                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) &
+                        & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) &
+                        &  > rn_wdmin1 + rn_wdmin2
+                ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + &
+                                                          & rn_wdmin1 + rn_wdmin2
+                IF(ll_tmp1) THEN
+                   zcpy(ji,jj) = 1.0_wp
+                   wdvflt1(ji,jj) = 1.0_wp
+                ELSE IF(ll_tmp2) THEN
+                   ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here
+                  zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /&
+                              &     (sshn(ji,jj+1) - sshn(ji,jj)))
+                   wdvflt1(ji,jj) = 1.0_wp
+                ELSE
+                  zcpy(ji,jj) = 0._wp
+                  wdvflt1(ji,jj) = 0.0_wp
+                END IF
+
+                !for w/d debugging, delete it when finished
+                !   zcpy(ji,jj) = 0._wp
+                ! end debugging part
+             END DO
+           END DO
+           CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
+        ENDIF
+
+        IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters
+           DO jj = 2, jpjm1
+              DO ji = 2, jpim1
+                 zu_trd(ji,jj) = (zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) ) / &
+                               &                 e1u(ji,jj)) * zcpx(ji,jj)
+                 zv_trd(ji,jj) = (zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) ) / &
+                               &                 e2v(ji,jj)) * zcpy(ji,jj)
+                 IF(ji+nimpp-1 == jidbg .and. jj+njmpp-1 == jjdbg) THEN
+                   write(333,*) 'in spt_ts:  zutrdx, zutrdy, hpgx, hpgy, zcpx, zcpy'
+                   write(333,*) zu_trd(ji,jj), zv_trd(ji,jj), &
+                             &  -grav*(sshn(ji+1,jj)-sshn(ji,jj))/e1u(ji,jj)*zcpx(ji,jj), &
+                             &  -grav*(sshn(ji,jj+1)-sshn(ji,jj))/e2v(ji,jj)*zcpy(ji,jj), &
+                             &  zcpx(ji,jj), zcpy(ji,jj) 
+                   write(334,*) 'in spt_ts:  sshn_ij, sshn_ip1j, sshn_ijp1, bathy_ij, bathy_ip1j, bathy_ijp1'
+                   write(334,*) sshn(ji,jj), sshn(ji+1,jj), sshn(ji,jj+1), bathy(ji,jj), bathy(ji+1,jj), bathy(ji,jj+1) 
+                 END IF
+              END DO
+           END DO
+         ELSE
+           DO jj = 2, jpjm1
+              DO ji = fs_2, fs_jpim1   ! vector opt.
+                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj)
+                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj) 
+              END DO
+           END DO
+        END IF
+
+      ENDIF
+
+      IF(ln_wd) THEN
+      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend
+         DO ji = fs_2, fs_jpim1
+             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1) * wduflt1(ji,jj)
+             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) * wdvflt1(ji,jj)
+          END DO
+      END DO
+      ELSE
       DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend
          DO ji = fs_2, fs_jpim1
@@ -392,5 +481,6 @@
              zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)
           END DO
-      END DO 
+      END DO
+      END IF
       !
       !						! Add bottom stress contribution from baroclinic velocities:      
@@ -416,7 +506,12 @@
       !
       ! Note that the "unclipped" bottom friction parameter is used even with explicit drag
-      zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:)
-      zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:)
-      !       
+      IF(ln_wd) THEN
+        zu_frc(:,:) = zu_frc(:,:) + MAX(hur(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:)
+        zv_frc(:,:) = zv_frc(:,:) + MAX(hvr(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:)
+      ELSE
+        zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:)
+        zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:)
+      END IF
+      !
       IF (ln_bt_fw) THEN                        ! Add wind forcing
          zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)
@@ -487,8 +582,14 @@
          vb_e   (:,:) = 0._wp
       ENDIF
-      !
-      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                    
-         sshn_e(:,:) = sshn (:,:)            
-         zun_e (:,:) = un_b (:,:)            
+
+      IF(ln_wd) THEN      !preserve the positivity of water depth
+                          !ssh[b,n,a] should have already been processed for this
+         sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:))
+         sshb_e(:,:)  = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:))
+      ENDIF
+      !
+      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields
+         sshn_e(:,:) = sshn (:,:)
+         zun_e (:,:) = un_b (:,:)
          zvn_e (:,:) = vn_b (:,:)
          !
@@ -561,4 +662,8 @@
             zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points
             zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
+            IF(ln_wd) THEN
+              zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1)
+              zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1)
+            END IF
          ELSE
             zhup2_e (:,:) = hu(:,:)
@@ -599,4 +704,5 @@
          ENDIF
 #endif
+         IF(ln_wd) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
          !
          ! Sum over sub-time-steps to compute advective velocities
@@ -613,4 +719,5 @@
          END DO
          ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
+         IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:)) 
          CALL lbc_lnk( ssha_e, 'T',  1._wp )
 
@@ -660,4 +767,54 @@
           &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
 
+         IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters
+           DO jj = 2, jpjm1
+              DO ji = 2, jpim1
+                 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))&
+                        & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj)) &
+                        &  > rn_wdmin1 + rn_wdmin2
+                 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + &
+                                                           & rn_wdmin1 + rn_wdmin2
+                 IF(ll_tmp1) THEN
+                   zcpx(ji,jj) = 1.0_wp
+                   wduflt1(ji,jj) = 1.0_wp
+                 ELSE IF(ll_tmp2) THEN
+                    ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here
+                   zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) /&
+                               & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj)))
+                   wduflt1(ji,jj) = 1.0_wp
+                 ELSE
+                    zcpx(ji,jj) = 0._wp
+                    wduflt1(ji,jj) = 0.0_wp
+                 END IF
+
+                !for w/d debugging, delete it when finished
+                !   zcpx(ji,jj) = 0._wp
+                ! end debugging part
+
+                 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))&
+                        & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1)) &
+                        &  > rn_wdmin1 + rn_wdmin2
+                 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + &
+                                                           & rn_wdmin1 + rn_wdmin2
+                 IF(ll_tmp1) THEN
+                    zcpy(ji,jj) = 1.0_wp
+                    wdvflt1(ji,jj) = 1.0_wp
+                 ELSE IF(ll_tmp2) THEN
+                    ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here
+                   zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) /&
+                               & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj)))
+                    wdvflt1(ji,jj) = 1.0_wp
+                 ELSE
+                   zcpy(ji,jj) = 0._wp
+                    wdvflt1(ji,jj) = 0.0_wp
+                 END IF
+
+                !for w/d debugging, delete it when finished
+                !   zcpy(ji,jj) = 0._wp
+                ! end debugging part
+              END DO
+            END DO
+            CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
+         ENDIF
          !
          ! Compute associated depths at U and V points:
@@ -676,4 +833,11 @@
                END DO
             END DO
+
+            IF(ln_wd) THEN
+              zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 )
+              zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 )
+            END IF
+            !shall we call lbc_lnk for zhu(v)st_e() here?
+
          ENDIF
          !
@@ -742,13 +906,37 @@
          !
          ! Surface pressure trend:
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ! Add surface pressure gradient
-               zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
-               zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
-               zwx(ji,jj) = zu_spg
-               zwy(ji,jj) = zv_spg
-            END DO
-         END DO
+
+         IF(ln_wd) THEN
+           DO jj = 2, jpjm1
+              DO ji = 2, jpim1 
+                 ! Add surface pressure gradient
+                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
+                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
+                 ! for W/D debugging
+                 !zwx(ji,jj) = zu_spg * zcpx(ji,jj) * 0._wp
+                 !zwy(ji,jj) = zv_spg * zcpy(ji,jj) * 0._wp
+                 IF(ji+nimpp-1 == jidbg .and. jj+njmpp-1 == jjdbg) THEN
+                   write(444,*) 'in spt_ts_intg:  zu_spg, zv_spg, zcpx, zcpy'
+                   write(444,*) zu_spg, zv_spg, zcpx(ji,jj), zcpy(ji,jj)
+                   write(445,*) 'in spt_ts_intg:  sshn_ij, sshn_ip1j, sshn_ijp1, bathy_ij, bathy_ip1j, bathy_ijp1'
+                   write(445,*) zsshp2_e(ji,jj), zsshp2_e(ji+1,jj), zsshp2_e(ji,jj+1), &
+                    &           bathy(ji,jj), bathy(ji+1,jj), bathy(ji,jj+1) 
+                 END IF
+                 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 
+                 zwy(ji,jj) = zv_spg * zcpy(ji,jj)
+              END DO
+           END DO
+         ELSE
+           DO jj = 2, jpjm1
+              DO ji = fs_2, fs_jpim1   ! vector opt.
+                 ! Add surface pressure gradient
+                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
+                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
+                 zwx(ji,jj) = zu_spg
+                 zwy(ji,jj) = zv_spg
+              END DO
+           END DO
+         END IF
+
          !
          ! Set next velocities:
@@ -774,9 +962,17 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
 
-                  zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1))
-                  zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))
-
-                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   & 
-                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
+                  IF(ln_wd) THEN
+                    zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1)
+                    zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1)
+                  ELSE
+                    zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
+                    zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
+                  END IF
+
+                  zhura = umask(ji,jj,1)/(zhura + 1._wp - umask(ji,jj,1))
+                  zhvra = vmask(ji,jj,1)/(zhvra + 1._wp - vmask(ji,jj,1))
+
+                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   &
+                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   &
                             &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
                             &               +      hu(ji,jj)  * zu_frc(ji,jj) ) &
@@ -793,7 +989,13 @@
          !
          IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only)
-            !                                          !  ----------------------------------------------        
-            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
-            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
+            !                                          !  ----------------------------------------------
+            IF(ln_wd) THEN
+              hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1)
+              hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1)
+            ELSE
+              hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
+              hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
+            END IF
+
             hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
             hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
@@ -925,4 +1127,5 @@
       CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   )
       CALL wrk_dealloc( jpi, jpj, zhf )
+      IF(ln_wd) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
Index: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 5839)
+++ branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 5842)
@@ -9,4 +9,5 @@
    !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
    !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work
+   !!            3.?  !  2015-11  (H. Liu) add wetting and drying  
    !!----------------------------------------------------------------------
 
@@ -38,4 +39,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE wadlmt          ! Wetting/Drying flux limting
 
    IMPLICIT NONE
@@ -90,8 +92,13 @@
       ENDIF
       !
-      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
-      !
       z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
       IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
+      !
+      z1_rau0 = 0.5_wp * r1_rau0
+      !
+      IF(ln_wd) CALL wad_lmt(sshb, z1_rau0 * (emp_b(:,:) + emp(:,:)), z2dt)
+
+      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
+      !
 
       !                                           !------------------------------!
@@ -106,5 +113,4 @@
       ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
       ! 
-      z1_rau0 = 0.5_wp * r1_rau0
       ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
 
Index: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90
===================================================================
--- branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90	(revision 5842)
+++ branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90	(revision 5842)
@@ -0,0 +1,340 @@
+
+MODULE wadlmt
+   !!==============================================================================
+   !!                       ***  MODULE  wadlmt  ***
+   !! compute and apply flux limiters and preserve water depth positivity
+   !! only effects if wetting/drying is on (ln_wd == .true.)
+   !!==============================================================================
+   !! History :   
+   !!  NEMO      3.6  ! 2014-09  ((H.Liu)  Original code
+   !!                 ! will add the runoff and periodic BC case later
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   wad_lmt    : Compute the horizontal flux limiter and the limited velocity
+   !!                when wetting and drying happens 
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE sbc_oce, ONLY : ln_rnf   ! surface boundary condition: ocean
+   USE sbcrnf          ! river runoff 
+   USE cla             ! cross land advection             (cla_div routine)
+   USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp         ! MPP library
+   USE wrk_nemo        ! Memory Allocation
+   USE timing          ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+
+   PUBLIC   wad_lmt, wad_lmt_bt    ! routine called by step.F90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+CONTAINS
+
+   SUBROUTINE wad_lmt( sshb1, sshemp, z2dt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE wad_lmt  ***
+      !!                    
+      !! ** Purpose :   generate flux limiters for wetting/drying
+      !!
+      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif) 
+      !!
+      !! ** Action  : - calculate flux limiter and W/D flag
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   sshb1
+      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   sshemp
+      REAL(wp), INTENT(in) :: z2dt
+      !
+      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
+      INTEGER  ::   zflag               ! local scalar
+      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
+      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
+      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
+      REAL(wp) ::   ztmp                ! local scalars
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv            ! local 2D workspace
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
+
+      !!----------------------------------------------------------------------
+      !
+
+      IF( nn_timing == 1 )  CALL timing_start('wad_lmt')
+
+      IF(ln_wd) THEN
+
+        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
+        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
+        !
+       
+        !IF(lwp) WRITE(numout,*)
+        !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting'
+       
+        zflag  = 0
+        zdepwd = 50._wp   !maximum depth on which that W/D could possibly happen
+
+       
+        zflxp(:,:)   = 0._wp
+        zflxn(:,:)   = 0._wp
+        zflxu(:,:)   = 0._wp
+        zflxv(:,:)   = 0._wp
+
+        zwdlmtu(:,:)  = 1._wp
+        zwdlmtv(:,:)  = 1._wp
+       
+        ! Horizontal Flux in u and v direction
+        DO jk = 1, jpkm1  
+           DO jj = 1, jpjm1
+              DO ji = 1, jpim1
+                 zflxu(ji,jj) = zflxu(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
+                 zflxv(ji,jj) = zflxv(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
+              END DO  
+           END DO  
+        END DO
+       
+        zflxu(:,:) = zflxu(:,:) * e2u(:,:)
+        zflxv(:,:) = zflxv(:,:) * e1v(:,:)
+       
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1 
+
+             IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE   ! we don't care about land cells
+             IF(bathy(ji,jj) > zdepwd) CYCLE       ! and cells which will unlikely go dried out
+
+              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
+                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
+              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
+                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
+
+              zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1
+              IF(zdep2 < 0._wp) THEN  !add more safty, but not necessary
+                !zdep2 = 0._wp
+                sshb1(ji,jj) = rn_wdmin1 - bathy(ji,jj)
+              END IF
+           ENDDO
+        END DO
+
+      
+        !! start limiter iterations 
+        DO jk1 = 1, nn_wdit + 1
+       
+          
+           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
+           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
+          
+           DO jj = 2, jpjm1
+              DO ji = 2, jpim1 
+        
+                 wdmask(ji,jj) = 0
+                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 
+                 IF(bathy(ji,jj) > zdepwd) CYCLE
+        
+                 !ztmp = e1t(ji,jj) * e2t(ji,jj)     !there must be an array ready for this
+                 ztmp = e12t(ji,jj)
+
+                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
+                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
+                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
+                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
+          
+                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
+                 zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)  ! this one can be moved out of the loop
+          
+                 IF(zdep1 > zdep2) THEN
+                   zflag = 1
+                   wdmask(ji, jj) = 1
+                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
+                   zcoef = max(zcoef, 0._wp)
+                   IF(jk1 > nn_wdit) zcoef = 0._wp
+                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
+                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
+                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
+                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef
+                 END IF
+              END DO ! ji loop
+           END DO  ! jj loop
+
+           CALL lbc_lnk( zwdlmtu, 'U', 1. )
+           CALL lbc_lnk( zwdlmtv, 'V', 1. )
+
+           IF(lk_mpp) CALL mpp_max(zflag)   !max over the global domain
+
+           IF(zflag == 0) EXIT
+          
+           zflag = 0     ! flag indicating if any further iteration is needed?
+        END DO  ! jk1 loop
+       
+        DO jk = 1, jpkm1
+          un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :) 
+          vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :) 
+        END DO
+
+        CALL lbc_lnk( un, 'U', -1. )
+        CALL lbc_lnk( vn, 'V', -1. )
+       
+        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
+       
+        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
+        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
+        !
+        !
+        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
+        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
+      !
+      END IF
+
+      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
+   END SUBROUTINE wad_lmt
+
+   SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE wad_lmt  ***
+      !!                    
+      !! ** Purpose :   limiting flux in the barotropic stepping (dynspg_ts)
+      !!
+      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif) 
+      !!
+      !! ** Action  : - calculate flux limiter and W/D flag
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in)    ::   rdtbt    ! ocean time-step index
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zflxu,  zflxv, sshn_e, zssh_frc  
+      !
+      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
+      INTEGER  ::   zflag         ! local scalar
+      REAL(wp) ::   z2dt
+      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
+      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
+      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
+      REAL(wp) ::   ztmp                ! local scalars
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
+
+      !!----------------------------------------------------------------------
+      !
+
+      IF( nn_timing == 1 )  CALL timing_start('wad_lmt_bt')
+
+      IF(ln_wd) THEN
+
+        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
+        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
+        !
+       
+        !IF(lwp) WRITE(numout,*)
+        !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting'
+       
+        zflag  = 0
+        zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes
+
+        z2dt = rdtbt   
+       
+        zflxp(:,:)   = 0._wp
+        zflxn(:,:)   = 0._wp
+        !zflxu(:,:)   = 0._wp
+        !zflxv(:,:)   = 0._wp
+
+        zwdlmtu(:,:)  = 1._wp
+        zwdlmtv(:,:)  = 1._wp
+       
+        ! Horizontal Flux in u and v direction
+       
+        !zflxu(:,:) = zflxu(:,:) * e2u(:,:)
+        !zflxv(:,:) = zflxv(:,:) * e1v(:,:)
+       
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1 
+
+             IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE   ! we don't care about land cells
+             IF(bathy(ji,jj) > zdepwd) CYCLE       ! and cells which will unlikely go dried out
+
+              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
+                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
+              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
+                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
+
+              zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
+              IF(zdep2 < 0._wp) THEN  !add more safty, but not necessary
+                !zdep2 = 0._wp
+               sshn_e(ji,jj) = rn_wdmin1 - bathy(ji,jj)
+              END IF
+           ENDDO
+        END DO
+
+      
+        !! start limiter iterations 
+        DO jk1 = 1, nn_wdit + 1
+       
+          
+           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
+           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
+          
+           DO jj = 2, jpjm1
+              DO ji = 2, jpim1 
+        
+                 wdmask(ji,jj) = 0
+                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 
+                 IF(bathy(ji,jj) > zdepwd) CYCLE
+        
+                 !ztmp = e1t(ji,jj) * e2t(ji,jj)     !there must be an array ready for this
+                 ztmp = e12t(ji,jj)
+
+                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
+                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
+                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
+                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
+          
+                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
+                 zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1   ! this one can be moved out of the loop
+                 zdep2 = zdep2 - z2dt * zssh_frc(ji,jj)
+          
+                 IF(zdep1 > zdep2) THEN
+                   zflag = 1
+                   !wdmask(ji, jj) = 1
+                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
+                   zcoef = max(zcoef, 0._wp)
+                   IF(jk1 > nn_wdit) zcoef = 0._wp
+                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
+                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
+                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
+                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef
+                 END IF
+              END DO ! ji loop
+           END DO  ! jj loop
+
+           CALL lbc_lnk( zwdlmtu, 'U', 1. )
+           CALL lbc_lnk( zwdlmtv, 'V', 1. )
+
+           IF(lk_mpp) CALL mpp_max(zflag)   !max over the global domain
+
+           IF(zflag == 0) EXIT
+          
+           zflag = 0     ! flag indicating if any further iteration is needed?
+        END DO  ! jk1 loop
+       
+        zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 
+        zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 
+
+        CALL lbc_lnk( zflxu, 'U', -1. )
+        CALL lbc_lnk( zflxv, 'V', -1. )
+       
+        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
+       
+        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
+        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
+        !
+        !
+        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
+        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
+      !
+      END IF
+
+      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
+   END SUBROUTINE wad_lmt_bt
+   !!======================================================================
+END MODULE wadlmt
Index: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 5839)
+++ branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 5842)
@@ -30,4 +30,5 @@
    !!            3.4  ! 2011-11  (C. Harris) decomposition changes for running with CICE
    !!                 ! 2012-05  (C. Calone, J. Simeon, G. Madec, C. Ethe) Add grid coarsening 
+   !!            3.6.?! 2015-11  (H. Liu) Add Wetting and Drying 
    !!----------------------------------------------------------------------
 
@@ -232,4 +233,5 @@
       NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, &
          &             jpizoom, jpjzoom, jperio, ln_use_jattr
+      NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit
       !!----------------------------------------------------------------------
       !
@@ -257,4 +259,14 @@
       READ  ( numnam_cfg, namcfg, IOSTAT = ios, ERR = 904 )
 904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfg in configuration namelist', .TRUE. )   
+
+      REWIND( numnam_ref )              ! Namelist namwad in reference namelist : Parameters for Wetting/Drying
+      READ  ( numnam_ref, namwad, IOSTAT = ios, ERR = 905)
+905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', .TRUE. )
+
+      REWIND( numnam_cfg )              ! Namelist namwad in configuration namelist : Parameters for Wetting/Drying
+      READ  ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906)
+906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', .TRUE. )
+
+
 
 ! Force values for AGRIF zoom (cf. agrif_user.F90)
@@ -315,4 +327,5 @@
          WRITE( numond, namctl )
          WRITE( numond, namcfg )
+         WRITE( numond, namwad )
       ENDIF
 
