Index: /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 5065)
+++ /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 5066)
@@ -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
    !!----------------------------------------------------------------------
 
@@ -262,4 +263,16 @@
    INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: nicoa, njcoa       !: ???
 #endif
+
+   !!----------------------------------------------------------------------
+   !! critical depths,limiters,and masks for  Wetting and Drying
+   !! ---------------------------------------------------------------------
+
+   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
 
    !!----------------------------------------------------------------------
@@ -325,5 +338,5 @@
    INTEGER FUNCTION dom_oce_alloc()
       !!----------------------------------------------------------------------
-      INTEGER, DIMENSION(11) :: ierr
+      INTEGER, DIMENSION(12) :: ierr
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -390,4 +403,7 @@
 #endif
       !
+      IF(ln_wd) &
+      ALLOCATE( wdmask(jpi,jpj), STAT=ierr(12) )
+      !
       dom_oce_alloc = MAXVAL(ierr)
       !
Index: /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 5065)
+++ /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 5066)
@@ -86,5 +86,6 @@
                              CALL dom_zgr      ! Vertical mesh and bathymetry
                              CALL dom_msk      ! Masks
-      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency
+      IF( ln_sco.AND.(.NOT.ln_wd) ) &
+                         &   CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency
       !
       ht_0(:,:) = 0.0_wp                       ! Reference ocean depth at T-points
Index: /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 5065)
+++ /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 5066)
@@ -145,4 +145,14 @@
       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
+            END IF
+          ENDDO
+        ENDDO
+      END IF
+
       ! Reconstruction of all vertical scale factors at now and before time steps
       ! =============================================================================
@@ -793,5 +803,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
       !!----------------------------------------------------------------------
@@ -867,4 +877,19 @@
             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
+                    fse3t_b(ji,jj,1:2) = 0.5_wp * rn_wdmin1
+                    fse3t_n(ji,jj,1:2) = 0.5_wp * rn_wdmin1
+                    sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj)
+                    sshn(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/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5065)
+++ /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5066)
@@ -17,4 +17,5 @@
    !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function
    !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case  
+   !!            3.6  ! 2014-09  (H. Liu) Modifications for Wetting/Drying
    !!----------------------------------------------------------------------
 
@@ -1121,4 +1122,5 @@
       REAL(wp) ::   zrmax, ztaper   ! temporary scalars
       REAL(wp) ::   zrfact
+      REAL(wp) ::   zsmth
       !
       REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
@@ -1176,9 +1178,11 @@
       bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
 
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-           IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
-         END DO
-      END DO
+      IF (.NOT. ln_wd) THEN
+        DO jj = 1, jpj
+           DO ji = 1, jpi
+             IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
+           END DO
+        END DO
+      ENDIF
       !                                        ! =============================
       !                                        ! Define the envelop bathymetry   (hbatt)
@@ -1187,19 +1191,21 @@
       zenv(:,:) = bathy(:,:)
       !
-      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-           IF( bathy(ji,jj) == 0._wp ) THEN
-             iip1 = MIN( ji+1, jpi )
-             ijp1 = MIN( jj+1, jpj )
-             iim1 = MAX( ji-1, 1 )
-             ijm1 = MAX( jj-1, 1 )
-             IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              &
-        &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
-               zenv(ji,jj) = rn_sbot_min
+      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
+           DO ji = 1, jpi
+             IF( bathy(ji,jj) == 0._wp ) THEN
+               iip1 = MIN( ji+1, jpi )
+               ijp1 = MIN( jj+1, jpj )
+               iim1 = MAX( ji-1, 1 )
+               ijm1 = MAX( jj-1, 1 )
+               IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              &
+          &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
+                 zenv(ji,jj) = rn_sbot_min
+               ENDIF
              ENDIF
-           ENDIF
-         END DO
-      END DO
+           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' )
@@ -1241,12 +1247,13 @@
          zri(:,:) = 0._wp
          zrj(:,:) = 0._wp
+         zsmth = 0._wp
          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
@@ -1271,9 +1278,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 )    ! fill 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
+      END IF
       !
       ! Envelope bathymetry saved in hbatt
@@ -1305,5 +1320,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
@@ -1318,4 +1338,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
@@ -1325,4 +1360,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)
@@ -1350,8 +1386,11 @@
 
 !!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
@@ -1395,14 +1434,16 @@
       CALL lbc_lnk( e3vw_0, 'V', 1._wp )
 
-      fsdepw(:,:,:) = gdepw_0 (:,:,:)
-      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
+      !fsdepw(:,:,:) = gdepw_0 (:,:,:)
+      !fsde3w(:,:,:) = gdep3w_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
@@ -1449,6 +1490,18 @@
                IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
 #endif
-               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
+#if key_surge
+                mbathy(ji,jj) = 1
+#else
+                mbathy(ji,jj) = 2
+#endif
+              ENDIF
+            ELSE
+              IF( scobot(ji,jj) == 0._wp            )   mbathy(ji,jj) = 0
+            END IF
          END DO
       END DO
@@ -1508,10 +1561,14 @@
          DO jj = 1, jpj
 
-            IF( hbatt(ji,jj) > 0._wp) THEN
+            IF( bathy(ji,jj) > 0._wp) THEN
                DO jk = 1, mbathy(ji,jj)
                  ! check coordinate is monotonically increasing
-                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
-                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
-                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
+                 !RF test...
+                 !IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
+                 !   WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
+                 !   WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
+                 IF (fse3w(ji,jj,jk) < 0._wp .OR. fse3t(ji,jj,jk) < 0._wp ) THEN
+                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   < 0  at point (i,j,k)= ', ji, jj, jk
+                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   < 0  at point (i,j,k)= ', ji, jj, jk
                     WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
                     WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
@@ -1572,4 +1629,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
@@ -1625,16 +1684,42 @@
       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
+
                !
                e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
@@ -1676,4 +1761,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
@@ -1754,16 +1841,41 @@
         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
 
              e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
Index: /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 5065)
+++ /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 5066)
@@ -16,4 +16,5 @@
    !!            3.4  !  2011-11  (H. Liu) hpg_prj: Original code for s-coordinates
    !!                 !           (A. Coward) suppression of hel, wdj and rot options
+   !!            3.6? !  2014-09  (H. Liu) add Wetting/Drying pressure filter
    !!----------------------------------------------------------------------
 
@@ -369,9 +370,12 @@
       !!
       INTEGER  ::   ji, jj, jk                 ! dummy loop indices
-      REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars
+      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp   ! temporary scalars
+      LOGICAL  ::   ll_tmp1, ll_tmp2           ! 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
@@ -386,4 +390,39 @@
       IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume
       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), 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))
+             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
 
@@ -402,4 +441,12 @@
             zvap = -zcoef0 * ( 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
@@ -424,4 +471,12 @@
                zvap = -zcoef0 * ( 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
@@ -445,4 +500,12 @@
             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
@@ -467,4 +530,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
@@ -476,4 +547,5 @@
       !
       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
+      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
       !
    END SUBROUTINE hpg_sco
@@ -493,8 +565,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
       !!----------------------------------------------------------------------
       !
@@ -502,4 +576,5 @@
       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 )
       !
 
@@ -508,4 +583,39 @@
          IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, density Jacobian with cubic polynomial scheme'
+      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), 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))
+             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
 
@@ -673,4 +783,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)
@@ -692,4 +806,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)
@@ -702,4 +820,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
@@ -727,4 +846,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
@@ -732,8 +852,10 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdept, zrhh
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zcpx, zcpy    !W/D pressure filter
       !!----------------------------------------------------------------------
       !
       CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
       CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh )
+      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
       !
       IF( kt == nit000 ) THEN
@@ -748,4 +870,39 @@
       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))
+             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))
+             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
@@ -907,6 +1064,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
 
@@ -964,6 +1124,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
 
@@ -975,4 +1139,5 @@
       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
       CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh )
+      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
       !
    END SUBROUTINE hpg_prj
Index: /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 5065)
+++ /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 5066)
@@ -11,4 +11,5 @@
    !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping
    !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility
+   !!             3.?  ! 2014-09  (H. Liu) Add Wetting/Drying components
    !!---------------------------------------------------------------------
 #if defined key_dynspg_ts   ||   defined key_esopa
@@ -40,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
@@ -137,14 +139,17 @@
       INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
       !
-      LOGICAL  ::   ll_fw_start			! if true, forward integration 
-      LOGICAL  ::   ll_init			    ! if true, special startup of 2d equations
-      INTEGER  ::   ji, jj, jk, jn   		! dummy loop indices
-      INTEGER  ::   ikbu, ikbv, noffset    	! local integers
-      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf  	! local scalars
-      REAL(wp) ::   zx1, zy1, zx2, zy2         !   -      -
-      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2 	  !   -      -
-      REAL(wp) ::   zu_spg, zv_spg             !   -      -
-      REAL(wp) ::   zhura, zhvra         	     !   -      -
-      REAL(wp) ::   za0, za1, za2, za3		       !   -      -
+      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
+      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars
+      REAL(wp) ::   zx1, zy1, zx2, zy2          !   -      -
+      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2     !   -      -
+      REAL(wp) ::   zu_spg, zv_spg              !   -      -
+      REAL(wp) ::   zhura, zhvra                !   -      -
+      REAL(wp) ::   za0, za1, za2, za3          !   -      -
+
+      REAL(wp) ::   ztmp                       ! temporary vaialbe
       !
       REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e
@@ -154,4 +159,5 @@
       REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a
       REAL(wp), POINTER, DIMENSION(:,:) :: zhf
+      REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef.
       !!----------------------------------------------------------------------
       !
@@ -165,4 +171,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 )
       !
       !                                         !* Local constant initialization
@@ -173,5 +180,5 @@
       zraur = 1._wp / rau0
       !
-      IF( kt == nit000 .AND. neuler == 0 ) THEN  	! reciprocal of baroclinic time step 
+      IF( kt == nit000 .AND. neuler == 0 ) THEN  ! reciprocal of baroclinic time step 
         z2dt_bf = rdt
       ELSE
@@ -180,11 +187,11 @@
       z1_2dt_b = 1.0_wp / z2dt_bf 
       !
-      ll_init = ln_bt_av                       		! if no time averaging, then no specific restart 
+      ll_init = ln_bt_av                       ! if no time averaging, then no specific restart 
       ll_fw_start = .FALSE.
       !
-                                     	                ! time offset in steps for bdy data update
+                                                     ! time offset in steps for bdy data update
       IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ;  noffset = 0 ; ENDIF
       !
-      IF( kt == nit000 ) THEN             	!* initialisation
+      IF( kt == nit000 ) THEN             !* initialisation
          !
          IF(lwp) WRITE(numout,*)
@@ -365,10 +372,54 @@
       !                                   ! ----------------------------------------------------
       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
+        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))
+                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))
+                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(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)
+              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
 
@@ -464,4 +515,11 @@
       !                                             ! ==================== !  
       ! Initialize barotropic variables:      
+
+      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 (:,:)            
@@ -537,4 +595,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(:,:)
@@ -575,4 +637,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
@@ -589,4 +652,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 )
 
@@ -652,4 +716,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
          !
@@ -717,14 +788,58 @@
          zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
          !
-         ! 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                    ! 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))
+                 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
+                 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)))
+                 ELSE
+                    zcpx(ji,jj) = 0._wp
+                 END IF
+
+                 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))
+                 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
+                 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)))
+                 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(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)
+                 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:
@@ -750,9 +865,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) ) &
@@ -770,6 +893,12 @@
          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) )
@@ -903,4 +1032,5 @@
       CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   )
       CALL wrk_dealloc( jpi, jpj, zhf )
+      IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
Index: /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 5065)
+++ /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 5066)
@@ -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.?  !  2014-09  (H. Liu) add wetting and drying  
    !!----------------------------------------------------------------------
 
@@ -42,4 +43,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE wadlmt          ! Wetting/Drying flux limting
 
    IMPLICIT NONE
@@ -94,9 +96,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
+      !
       !                                           !------------------------------!
       !                                           !   After Sea Surface Height   !
@@ -110,5 +116,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(:,:) )  ) * tmask(:,:,1)
 
Index: /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90
===================================================================
--- /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90	(revision 5066)
+++ /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90	(revision 5066)
@@ -0,0 +1,338 @@
+
+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 that ocean cells can have W/D processes
+
+       
+        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
+       
+        un(:,:,:) = un(:,:,:) * zwdlmtu(ji, jj) 
+        vn(:,:,:) = vn(:,:,:) * zwdlmtv(ji, jj) 
+
+        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(ji, jj) 
+        zflxv(:,:) = zflxv(:,:) * zwdlmtv(ji, jj) 
+
+        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/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 5065)
+++ /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 5066)
@@ -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.?! 2014-09  (H. Liu) Add Wetting and Drying 
    !!----------------------------------------------------------------------
 
@@ -233,4 +234,5 @@
       NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, &
          &             jpizoom, jpjzoom, jperio
+      NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit
       !!----------------------------------------------------------------------
       !
@@ -257,4 +259,12 @@
       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', lwp )
+
+      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', lwp )
 
 ! Force values for AGRIF zoom (cf. agrif_user.F90)
@@ -311,4 +321,5 @@
          WRITE( numond, namctl )
          WRITE( numond, namcfg )
+         WRITE( numond, namwad )
       ENDIF
 
@@ -527,4 +538,17 @@
          WRITE(numout,*) '      lateral cond. type (between 0 and 6) jperio = ', jperio   
       ENDIF
+
+      IF(lwp) THEN                  ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'namwad  : wetting and drying initialisation through namelist read'
+         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) '   Namelist namwad'
+         WRITE(numout,*) '   key to turn on/off wetting/drying              ln_wd      = ', ln_wd
+         WRITE(numout,*) '   minimum water depth on dried cells             rn_wdmin1  = ', rn_wdmin1
+         WRITE(numout,*) '   tolerrance of min water depth on dried cells   rn_wdmin2  = ', rn_wdmin2
+         WRITE(numout,*) '   land elevation below which wad considered      rn_wdld    = ', rn_wdld
+         WRITE(numout,*) '   maximum number of iteration for W/D limiter    nn_wdit    = ', nn_wdit
+      ENDIF 
+
       !                             ! Parameter control
       !
Index: /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5065)
+++ /branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5066)
@@ -118,4 +118,5 @@
    USE trcstp           ! passive tracer time-stepping      (trc_stp routine)
 #endif
+   USE wadlmt           ! wetting/drying limiter           (wad_lmt routine)
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
