Index: branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 4376)
+++ branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 4380)
@@ -372,24 +372,63 @@
             h_oce     = gdepw_0(jpk)
          ELSE                                         ! bump centered in the basin
+            !IF(lwp) WRITE(numout,*)
+            !IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump'
+            !ii_bump = jpidta / 2                           ! i-index of the bump center
+            !ij_bump = jpjdta / 2                           ! j-index of the bump center
+            !r_bump  = 50000._wp                            ! bump radius (meters)       
+            !h_bump  =  2700._wp                            ! bump height (meters)
+            !h_oce   = gdepw_0(jpk)                         ! background ocean depth (meters)
+            !IF(lwp) WRITE(numout,*) '            bump characteristics: '
+            !IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump
+            !IF(lwp) WRITE(numout,*) '               bump height         = ', h_bump , ' meters'
+            !IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' index'
+            !IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters'
+            !!                                        
+            !DO jj = 1, jpjdta                              ! zdta :
+            !   DO ji = 1, jpidta
+            !      zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
+            !      zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
+            !      zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) )
+            !   END DO
+            !END DO
+
             IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump'
-            ii_bump = jpidta / 2                           ! i-index of the bump center
-            ij_bump = jpjdta / 2                           ! j-index of the bump center
-            r_bump  = 50000._wp                            ! bump radius (meters)       
-            h_bump  =  2700._wp                            ! bump height (meters)
+            IF(lwp) WRITE(numout,*) '         bathymetry field: Thacker''s parabolic basin'
+            ii_bump = jpidta / 2 + 1                       ! i-index of the basin center
+            ij_bump = jpjdta / 2 + 1                       ! j-index of the basin center
+            r_bump  = 430620._wp                           ! basin radius (meters)       
+            h_bump  =     50._wp                           ! basin depth  (meters)
             h_oce   = gdepw_0(jpk)                         ! background ocean depth (meters)
-            IF(lwp) WRITE(numout,*) '            bump characteristics: '
-            IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump
-            IF(lwp) WRITE(numout,*) '               bump height         = ', h_bump , ' meters'
-            IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' index'
+            IF(lwp) WRITE(numout,*) '           basin characteristics: '
+            IF(lwp) WRITE(numout,*) '              basin center (i,j)   = ', ii_bump, ii_bump
+            IF(lwp) WRITE(numout,*) '              basin depth          = ', h_bump , ' meters'
+            IF(lwp) WRITE(numout,*) '              basin radius         = ', r_bump , ' index'
             IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters'
-            !                                        
+                                                    
             DO jj = 1, jpjdta                              ! zdta :
                DO ji = 1, jpidta
-                  zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
-                  zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
-                  zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) )
+                  zi = FLOAT( ji - ii_bump ) * ppe1_m 
+                  zj = FLOAT( jj - ij_bump ) * ppe2_m 
+                  zdta(ji,jj) = h_bump * ( 1._wp - ( zi*zi + zj*zj ) / (r_bump * r_bump) )
                END DO
             END DO
+            !IF(lwp) WRITE(numout,*)
+            !IF(lwp) WRITE(numout,*) '         bathymetry field: Thacker''s parabolic channel'
+            !ii_bump = jpidta / 2                           ! i-index of the bump center
+            !r_bump  = 81000._wp                            ! bump radius (meters)       
+            !h_bump  =    20._wp                            ! bump height (meters)
+            !h_oce   = gdepw_0(jpk)                         ! background ocean depth (meters)
+            !IF(lwp) WRITE(numout,*) '            channel characteristics: '
+            !IF(lwp) WRITE(numout,*) '               channel center (i,j)   = ', ii_bump, ii_bump
+            !IF(lwp) WRITE(numout,*) '               channel depth          = ', h_bump , ' meters'
+            !IF(lwp) WRITE(numout,*) '               channel radius         = ', r_bump , ' index'
+            !IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters'
+            !                                        
+            !DO jj = 1, jpjdta                              ! zdta :
+            !   DO ji = 1, jpidta
+            !      zi = FLOAT( ji - ii_bump ) * ppe1_m 
+            !      zdta(ji,jj) = h_bump * ( 1._wp - (zi/r_bump) ** 2) - 10._wp
+            !   END DO
+            !END DO
             !                                              ! idta :
             IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk
@@ -1418,5 +1457,5 @@
                   IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) THEN
                     mbathy(ji,jj) = MAX( 2, jk )
-                  ELSE IF(scobot(ji,jj) <= rn_landele) THEN
+                  ELSE IF(scobot(ji,jj) <= - rn_landele) THEN
                     mbathy(ji,jj) = 0
                   ELSE
Index: branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 4376)
+++ branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 4380)
@@ -127,5 +127,5 @@
       REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum
 
-      REAL(wp), POINTER, DIMENSION(:,:) :: zwadflt
+      REAL(wp), POINTER, DIMENSION(:,:) :: zwadfltu, zwadfltv
       !!----------------------------------------------------------------------
       !
@@ -136,5 +136,5 @@
       CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum )
 
-      IF(ln_wad) CALL wrk_alloc( jpi, jpj, zwadflt)
+      IF(ln_wad) CALL wrk_alloc( jpi, jpj, zwadfltu, zwadfltv)
       !
       IF( kt == nit000 ) THEN             !* initialisation
@@ -185,5 +185,8 @@
       ENDIF
 
-      IF(ln_wad) zwadflt(:,:) = 1._wp
+      IF(ln_wad) THEN
+        zwadfltu(:,:) = 1._wp
+        zwadfltv(:,:) = 1._wp
+      END IF
 
       ! -----------------------------------------------------------------------------
@@ -412,12 +415,25 @@
          !                                                !* after ssh_e
          !                                                !  -----------
-         DO jj = 2, jpjm1                                 ! Horizontal divergence of barotropic transports
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     &
-                  &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     &
-                  &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     &
-                  &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) / ( e1t(ji,jj) * e2t(ji,jj) )
-            END DO
-         END DO
+         IF(ln_wad) THEN
+            DO jj = 2, jpjm1                                 ! Horizontal divergence of barotropic transports
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj) * zwadfltu(ji,   jj  )   &
+                     &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj) * zwadfltu(ji-1, jj  )   &
+                     &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  ) * zwadfltv(ji,   jj  )   &
+                     &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1) * zwadfltv(ji,   jj-1) ) &
+                     &         / ( e1t(ji,jj) * e2t(ji,jj) )
+               END DO
+            END DO
+         ELSE
+            DO jj = 2, jpjm1                                 ! Horizontal divergence of barotropic transports
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     &
+                     &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     &
+                     &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     &
+                     &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) &
+                     &         / ( e1t(ji,jj) * e2t(ji,jj) )
+               END DO
+            END DO
+         END IF
          !
 #if defined key_obc
@@ -433,25 +449,27 @@
 #endif
          !
+         DO jj = 2, jpjm1                                      ! leap-frog on ssh_e
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * &
+                             &   ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 
+            END DO
+         END DO
+
+         !! generate W/D filter
          IF(ln_wad) THEN
-            DO jj = 2, jpjm1                                      ! leap-frog on ssh_e
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * &
-                                &   ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 
-                  IF(ssha_e(ji,jj) <= rn_wadmin) THEN
-                    zwadflt(ji,   jj  ) = 0._wp
-                    zwadflt(ji-1, jj  ) = 0._wp
-                    zwadflt(ji,   jj-1) = 0._wp
-                    zwadflt(ji-1, jj-1) = 0._wp
+                  IF(ssha_e(ji,jj) + bathy(ji,jj) <= rn_wadmin) THEN
+                    zwadfltu(ji,   jj  ) = 0._wp
+                    zwadfltu(ji-1, jj  ) = 0._wp
+                    zwadfltv(ji,   jj  ) = 0._wp
+                    zwadfltv(ji,   jj-1) = 0._wp
                   END IF
                END DO
             END DO
-         ELSE
-            DO jj = 2, jpjm1                                      ! leap-frog on ssh_e
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * &
-                                &   ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 
-               END DO
-            END DO
+            CALL lbc_lnk(zwadfltu, 'U', 1._wp) 
+            CALL lbc_lnk(zwadfltv, 'V', 1._wp) 
          END IF
+
 
          !                                                !* after barotropic velocities (vorticity scheme dependent)
@@ -622,4 +640,13 @@
 #endif
 
+         IF(ln_wad) THEN
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                 ua_e(ji,jj) = ua_e(ji,jj) * zwadfltu(ji,jj)
+                 va_e(ji,jj) = va_e(ji,jj) * zwadfltv(ji,jj)
+               END DO
+            END DO
+         END IF
+
          !
          CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries 
@@ -693,6 +720,6 @@
       IF(ln_wad) THEN
          DO jk=1,jpkm1
-            ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b * zwadflt(:,:)
-            va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b * zwadflt(:,:)
+            ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b * zwadfltu(:,:)
+            va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b * zwadfltv(:,:)
          END DO
       ELSE
@@ -714,5 +741,5 @@
       CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum )
 
-      IF(ln_wad) CALL wrk_dealloc( jpi, jpj, zwadflt)
+      IF(ln_wad) CALL wrk_dealloc( jpi, jpj, zwadfltu, zwadfltv)
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
