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 5275)
+++ branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 5727)
@@ -145,13 +145,15 @@
       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
+!      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
@@ -703,5 +705,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)   * 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) ) )
@@ -721,7 +724,8 @@
             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)                                   &
-                     &                       * (   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) ) )
+!                  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)        &
+                        &                       * (   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) ) )
                END DO
             END DO
@@ -739,5 +743,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)      &
+                     &                       * 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) ) )
@@ -882,9 +888,12 @@
                 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
+                  !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
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 5275)
+++ branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5727)
@@ -300,4 +300,5 @@
       zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d )                   ! ref. depth with tolerance (10% of minimum layer thickness)
       nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
+      nlb10 = MAX(nlb10, 2)                                          ! prevent nla10 = 0 
       nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
 !!gm end bug
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 5275)
+++ branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 5727)
@@ -395,31 +395,39 @@
         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)) +&
+             IF ( tmask(ji+1,jj,1) == 0._wp) THEN
+                zcpx(ji,jj) = 1.0_wp
+             ELSE
+                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
+             ENDIF
+        
+             IF ( tmask(ji,jj+1,1) == 0._wp) THEN
+                zcpy(ji,jj) = 1.0_wp
+             ELSE
+                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
-               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
+                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
+             ENDIF
            END DO
         END DO
@@ -450,6 +458,6 @@
 
             ! add to the general momentum trend
-            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
-            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
+            ua(ji,jj,1) = ua(ji,jj,1) + ( zhpi(ji,jj,1) + zuap ) * umask(ji,jj,1)
+            va(ji,jj,1) = va(ji,jj,1) + ( zhpj(ji,jj,1) + zvap ) * vmask(ji,jj,1)
          END DO
       END DO
@@ -480,6 +488,6 @@
 
                ! add to the general momentum trend
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
-               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
+               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap ) * umask(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk)
             END DO
          END DO
@@ -546,4 +554,5 @@
 #endif
       !
+
       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
       IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
@@ -852,4 +861,5 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdept, zrhh
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   sshu_n, sshv_n
       REAL(wp), POINTER, DIMENSION(:,:)   ::   zcpx, zcpy    !W/D pressure filter
       !!----------------------------------------------------------------------
@@ -857,4 +867,5 @@
       CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
       CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh )
+      CALL wrk_alloc( jpi,jpj, sshu_n, sshv_n )
       IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
       !
@@ -967,8 +978,28 @@
 
       ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1)
+
+      ! The following modification "sshu_n -> sshn" is a big mistake, this
+      ! should never happen here, remember to correct this in NEMO v3.6
+      ! trunk.  H.L.
+
+      !prepare sshu_n and sshv_n
+      DO jj = 1, jpjm1
+        DO ji = 1, jpim1
+          sshu_n(ji,jj) = (e12u(ji,jj) * sshn(ji,jj) + e12u(ji+1, jj) * sshn(ji+1,jj)) * &
+                        & r1_e12u(ji,jj) * umask(ji,jj,1) * 0.5_wp
+          sshv_n(ji,jj) = (e12v(ji,jj) * sshn(ji,jj) + e12v(ji+1, jj) * sshn(ji,jj+1)) * &
+                        & r1_e12v(ji,jj) * vmask(ji,jj,1) * 0.5_wp
+        END DO
+      END DO
+
+      CALL lbc_lnk (sshu_n, 'U', 1.)
+      CALL lbc_lnk (sshv_n, 'V', 1.)
+
       DO jj = 2, jpjm1
         DO ji = 2, jpim1
-          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshn(ji,jj) * znad)    ! probable bug: changed from sshu_n for ztilde compilation
-          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshn(ji,jj) * znad)    ! probable bug: changed from sshv_n for ztilde compilation
+!          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshn(ji,jj) * znad)    ! probable bug: changed from sshu_n for ztilde compilation
+!          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshn(ji,jj) * znad)    ! probable bug: changed from sshv_n for ztilde compilation
+          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshu_n(ji,jj) * znad)
+          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshv_n(ji,jj) * znad)
         END DO
       END DO
@@ -1139,4 +1170,5 @@
       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
       CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh )
+      CALL wrk_dealloc( jpi,jpj, sshu_n, sshv_n )
       IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
       !
Index: branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
===================================================================
--- branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 5275)
+++ branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 5727)
@@ -99,6 +99,6 @@
          DO jj = 2, jpjm1       ! add the gradient of kinetic energy to the general momentum trends
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
-               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
+               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk)
             END DO 
          END DO
Index: branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90
===================================================================
--- branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90	(revision 5275)
+++ branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90	(revision 5727)
@@ -144,7 +144,7 @@
          
          ! Multiply by the eddy viscosity coef. (at u- and v-points)
-         zlu(:,:,jk) = zlu(:,:,jk) * ( fsahmu(:,:,jk) * (1-nkahm_smag) + nkahm_smag)
-
-         zlv(:,:,jk) = zlv(:,:,jk) * ( fsahmv(:,:,jk) * (1-nkahm_smag) + nkahm_smag)
+         zlu(:,:,jk) = zlu(:,:,jk) * ( fsahmu(:,:,jk) * (1-nkahm_smag) + nkahm_smag) * umask(:,:,jk)
+ 
+         zlv(:,:,jk) = zlv(:,:,jk) * ( fsahmv(:,:,jk) * (1-nkahm_smag) + nkahm_smag) * vmask(:,:,jk)
          
          ! Contravariant "laplacian"
@@ -174,5 +174,5 @@
                zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
                zut(ji,jj,jk) = (  zlu(ji,jj,jk) - zlu(ji-1,jj  ,jk)   &
-                  &             + zlv(ji,jj,jk) - zlv(ji  ,jj-1,jk) ) / zbt
+                  &             + zlv(ji,jj,jk) - zlv(ji  ,jj-1,jk) ) / zbt * tmask(ji,jj,jk)
             END DO
          END DO
@@ -201,6 +201,6 @@
                   &  + ( zut(ji,jj+1,jk) - zut(ji  ,jj,jk) ) / e2v(ji,jj)
                ! add it to the general momentum trends
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zua * ( fsahmu(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag ))
-               va(ji,jj,jk) = va(ji,jj,jk) + zva * ( fsahmv(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag ))
+               ua(ji,jj,jk) = ua(ji,jj,jk) + zua * ( fsahmu(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag )) * umask(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) + zva * ( fsahmv(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag )) * vmask(ji,jj,jk)
             END DO
          END DO
Index: branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 5275)
+++ branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 5727)
@@ -159,6 +159,6 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
-                  va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
+                  ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) * umask(ji,jj,jk)
+                  va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) * vmask(ji,jj,jk)
                END DO
             END DO
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 5275)
+++ branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 5727)
@@ -371,33 +371,46 @@
       !                                   !* Right-Hand-Side of the barotropic momentum equation
       !                                   ! ----------------------------------------------------
+
       IF( lk_vvl ) THEN                         ! Variable volume : remove surface pressure gradient
         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)) + &
+                IF ( tmask(ji+1,jj,1) == 0._wp ) THEN
+                   zcpx = 1._wp
+                ELSE
+                   ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji+1,jj)*tmask(ji+1,jj,1)) &
+                           & .and. MAX(sshn(ji,jj) + bathy(ji,jj)*tmask(ji,jj,1), sshn(ji+1,jj) + bathy(ji+1,jj)*tmask(ji+1,jj,1)) &
+                           &  > rn_wdmin1 + rn_wdmin2
+                   ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji+1,jj)*tmask(ji+1,jj,1)) + &
                                                           & 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)))
+                   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)*tmask(ji+1,jj,1) - sshn(ji,jj) - bathy(ji,jj)*tmask(ji,jj,1)) /&
+                                  &     (sshn(ji+1,jj) - sshn(ji,jj)))
+                   ELSE
+                      zcpx(ji,jj) = 0._wp
+                   END IF
+                ENDIF
+
+                IF ( tmask(ji,jj+1,1) == 0._wp ) THEN
+                   zcpy = 1._wp
                 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
+                   ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji,jj+1)*tmask(ji,jj+1,1)) &
+                           & .and. MAX(sshn(ji,jj) + bathy(ji,jj)*tmask(ji,jj,1), sshn(ji,jj+1) + bathy(ji,jj+1)*tmask(ji,jj+1,1)) &
+                           &  > rn_wdmin1 + rn_wdmin2
+                   ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji,jj+1)*tmask(ji,jj+1,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)*tmask(ji,jj+1,1) - sshn(ji,jj) - bathy(ji,jj)*tmask(ji,jj,1)) /&
+                                 &     (sshn(ji,jj+1) - sshn(ji,jj)))
+                   ELSE
+                     zcpy(ji,jj) = 0._wp
+                   END IF
+                ENDIF
              END DO
            END DO
@@ -453,6 +466,11 @@
       !
       ! 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
@@ -518,6 +536,6 @@
       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(:,:))
+         sshbb_e(:,:jj) = MAX( sshbb_e(:,:), (rn_wdmin1 - bathy(:,:)) ) *tmask(:,:,1)
+         sshb_e(:,:jj)  = MAX( sshb_e(:,:) , (rn_wdmin1 - bathy(:,:)) ) *tmask(:,:,1)
       ENDIF
 
@@ -637,4 +655,5 @@
          ENDIF
 #endif
+
          IF(ln_wd) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
          !
@@ -652,5 +671,5 @@
          END DO
          ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
-         IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:)) 
+         IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), (rn_wdmin1 - bathy(:,:)) ) * tmask(:,:,1)
          CALL lbc_lnk( ssha_e, 'T',  1._wp )
 
@@ -785,35 +804,52 @@
          !
          ! Add bottom stresses:
-         zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)
-         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
-         !
+         IF(ln_wd) THEN
+           zu_trd(:,:) = zu_trd(:,:) + MAX(bfrua(:,:) * hur_e(:,:), -1._wp / rdtbt) * zun_e(:,:)
+           zv_trd(:,:) = zv_trd(:,:) + MAX(bfrva(:,:) * hvr_e(:,:), -1._wp / rdtbt) * zvn_e(:,:)
+         ELSE
+           zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)
+           zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
+         END IF
+         
          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)) + &
+                 IF ( tmask(ji+1,jj,1) == 0._wp ) THEN
+                    zcpx = 1._wp
+                 ELSE
+                    ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji+1,jj)*tmask(ji+1,jj,1)) &
+                           & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj)*tmask(ji,jj,1), zsshp2_e(ji+1,jj) + bathy(ji+1,jj)*tmask(ji+1,jj,1))  &
+                           &  > rn_wdmin1 + rn_wdmin2
+                    ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji+1,jj)*tmask(ji+1,jj,1)) + &
                                                            & 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)))
+                    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)*tmask(ji+1,jj,1) - zsshp2_e(ji,jj) - bathy(ji,jj)*tmask(ji,jj,1)) /&
+                                 & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj)))
+                    ELSE
+                       zcpx(ji,jj) = 0._wp
+                    END IF
+                 ENDIF
+
+                 IF ( tmask(ji,jj+1,1) == 0._wp ) THEN
+                    zcpy = 1._wp
                  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)) + &
+                    ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji,jj+1)*tmask(ji,jj+1,1))&
+                           & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj)*tmask(ji,jj,1), zsshp2_e(ji,jj+1) + bathy(ji,jj+1)*tmask(ji,jj+1,1)) &
+                           &  > rn_wdmin1 + rn_wdmin2
+                    ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji,jj+1)*tmask(ji,jj+1,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
+                    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)*tmask(ji,jj+1,1) - zsshp2_e(ji,jj) - bathy(ji,jj)*tmask(ji,jj,1)) /&
+                                  & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj)))
+                    ELSE
+                      zcpy(ji,jj) = 0._wp
+                    END IF
+                 ENDIF
               END DO
             END DO
@@ -894,6 +930,6 @@
             !                                          !  ----------------------------------------------        
             IF(ln_wd) THEN
-              hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1)
-              hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1)
+              hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1 * umask(:,:,1) )
+              hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1 * vmask(:,:,1) )
             ELSE
               hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
@@ -1032,5 +1068,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(ln_wd) CALL wrk_dealloc( 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/dynvor.F90
===================================================================
--- branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 5275)
+++ branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 5727)
@@ -288,6 +288,6 @@
                zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
                zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
-               pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
-               pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
+               pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) * umask(ji,jj,jk)
+               pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) * vmask(ji,jj,jk)
             END DO  
          END DO  
@@ -409,6 +409,6 @@
                zcva =-zfact2 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
                ! mixed vorticity trend added to the momentum trends
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua
-               va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva
+               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zcua + zua ) * umask(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) + ( zcva + zva ) * vmask(ji,jj,jk)
             END DO  
          END DO  
@@ -529,6 +529,6 @@
                zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
                   &                         + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) )
-               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
-               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
+               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) * umask(ji,jj,jk)
+               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) * vmask(ji,jj,jk)
             END DO  
          END DO  
@@ -686,6 +686,6 @@
                zva = - zfac12 / e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
                   &                           + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
-               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
-               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
+               pua(ji,jj,jk) = pua(ji,jj,jk) + zua * umask(ji,jj,jk)
+               pva(ji,jj,jk) = pva(ji,jj,jk) + zva * vmask(ji,jj,jk)
             END DO  
          END DO  
Index: branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
===================================================================
--- branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 5275)
+++ branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 5727)
@@ -109,6 +109,6 @@
                zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
                !                         ! add the trends to the general momentum trends
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
-               va(ji,jj,jk) = va(ji,jj,jk) + zva
+               ua(ji,jj,jk) = ua(ji,jj,jk) + zua * umask(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) + zva * vmask(ji,jj,jk)
             END DO  
          END DO  
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 5275)
+++ branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90	(revision 5727)
@@ -131,5 +131,5 @@
               DO ji = 2, jpim1 
         
-                 wdmask(ji,jj) = 0
+!                 wdmask(ji,jj) = 0
                  IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 
                  IF(bathy(ji,jj) > zdepwd) CYCLE
@@ -148,5 +148,5 @@
                  IF(zdep1 > zdep2) THEN
                    zflag = 1
-                   wdmask(ji, jj) = 1
+!                   wdmask(ji, jj) = 1
                    zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
                    zcoef = max(zcoef, 0._wp)
@@ -170,6 +170,8 @@
         END DO  ! jk1 loop
        
-        un(:,:,:) = un(:,:,:) * zwdlmtu(ji, jj) 
-        vn(:,:,:) = vn(:,:,:) * zwdlmtv(ji, jj) 
+        DO jk = 1, jpkm1
+          un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :) 
+          vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :) 
+        END DO
 
         CALL lbc_lnk( un, 'U', -1. )
@@ -213,4 +215,5 @@
       REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
       REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   sum_e3u, sum_e3v         ! local 2D workspace
 
       !!----------------------------------------------------------------------
@@ -223,4 +226,5 @@
         CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
         CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
+        CALL wrk_alloc( jpi, jpj, sum_e3u, sum_e3v)
         !
        
@@ -235,6 +239,6 @@
         zflxp(:,:)   = 0._wp
         zflxn(:,:)   = 0._wp
-        zflxu(:,:)   = 0._wp
-        zflxv(:,:)   = 0._wp
+!RF bug fix!        zflxu(:,:)   = 0._wp
+!RF bug fix!        zflxv(:,:)   = 0._wp
 
         zwdlmtu(:,:)  = 1._wp
@@ -243,6 +247,6 @@
         ! Horizontal Flux in u and v direction
        
-        zflxu(:,:) = zflxu(:,:) * e2u(:,:)
-        zflxv(:,:) = zflxv(:,:) * e1v(:,:)
+!RF bug fix        zflxu(:,:) = zflxu(:,:) * e2u(:,:)
+!RF bug fix        zflxv(:,:) = zflxv(:,:) * e1v(:,:)
        
         DO jj = 2, jpjm1
@@ -276,9 +280,8 @@
               DO ji = 2, jpim1 
         
-                 wdmask(ji,jj) = 0
-                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 
-                 IF(bathy(ji,jj) > zdepwd) CYCLE
+!                 wdmask(ji,jj) = 0
+                 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 
         
-                 !ztmp = e1t(ji,jj) * e2t(ji,jj)     !there must be an array ready for this
                  ztmp = e12t(ji,jj)
 
@@ -316,6 +319,6 @@
         END DO  ! jk1 loop
        
-        zflxu(:,:) = zflxu(:,:) * zwdlmtu(ji, jj) 
-        zflxv(:,:) = zflxv(:,:) * zwdlmtv(ji, jj) 
+        zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 
+        zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 
 
         CALL lbc_lnk( zflxu, 'U', -1. )
@@ -330,4 +333,5 @@
         CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
         CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
+        CALL wrk_dealloc( jpi, jpj, sum_e3u, sum_e3v)
       !
       END IF
