Index: branches/2013/dev_r4050_NOC_WaD/NEMOGCM/CONFIG/cfg.txt
===================================================================
--- branches/2013/dev_r4050_NOC_WaD/NEMOGCM/CONFIG/cfg.txt	(revision 4356)
+++ branches/2013/dev_r4050_NOC_WaD/NEMOGCM/CONFIG/cfg.txt	(revision 4375)
@@ -9,2 +9,3 @@
 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC
 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC
+PBASIN OPA_SRC
Index: branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 4356)
+++ branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 4375)
@@ -60,4 +60,5 @@
    REAL(wp), PUBLIC ::   rn_wadmine = 0.01_wp   !: tolerrance of minimum water depth on dried cells
    REAL(wp), PUBLIC ::   rn_landele = 20.0_wp   !: land elevation below which wetting/drying will be considered
+   INTEGER , PUBLIC ::   nn_waditr  =   10      !: Maximum number of iteration for W/D limiter
 
    !                                                  !!! associated variables
Index: branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 4356)
+++ branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 4375)
@@ -280,5 +280,5 @@
       !!----------------------------------------------------------------------
       USE ioipsl
-      NAMELIST/namwad/ ln_wad , rn_wadmin, rn_wadmine, rn_landele
+      NAMELIST/namwad/ ln_wad , rn_wadmin, rn_wadmine, rn_landele, nn_waditr
       !!----------------------------------------------------------------------
 
@@ -295,4 +295,5 @@
          WRITE(numout,*) '      tollerance of minimual depth    rn_wadmine    = ', rn_wadmine
          WRITE(numout,*) ' threshold elevation for land masking rn_landele    = ', rn_landele
+         WRITE(numout,*) ' Maximum number of iteration for W/D limiter        = ', nn_waditr
       ENDIF
 
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 4356)
+++ branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 4375)
@@ -1411,12 +1411,29 @@
 !!
       ! HYBRID : 
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            DO jk = 1, jpkm1
-               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
-               IF( scobot(ji,jj) == 0._wp            )   mbathy(ji,jj) = 0
-            END DO
-         END DO
-      END DO
+      IF(ln_wad) THEN
+      !! Wetting/Drying case
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               DO jk = 1, jpkm1
+                  IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) THEN
+                    mbathy(ji,jj) = MAX( 2, jk )
+                  ELSE IF(scobot(ji,jj) <= rn_landele) THEN
+                    mbathy(ji,jj) = 0
+                  ELSE
+                    mbathy(ji,jj) = 2
+                  END IF
+               END DO
+            END DO
+         END DO
+      ELSE
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               DO jk = 1, jpkm1
+                  IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
+                  IF( scobot(ji,jj) == 0._wp            )   mbathy(ji,jj) = 0
+               END DO
+            END DO
+         END DO
+      END IF
       IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
          &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
Index: branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 4356)
+++ branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 4375)
@@ -76,5 +76,9 @@
       !!             - Save the trend (l_trddyn=T)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   kt                               ! ocean time-step index
+      INTEGER             ::   ji, jj, jk                       ! dummy loop indices
+      REAL                ::   zcpx, zcpy                       ! gravity filter for W/D
+      REAL                ::   zij, zim1j, zijm1, zim1jm1   ! local variables
+      REAL                ::   dij, dim1j, dijm1, dim1jm1   ! local variables
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
       !!----------------------------------------------------------------------
@@ -82,9 +86,11 @@
       IF( nn_timing == 1 )  CALL timing_start('dyn_hpg')
       !
-      IF( l_trddyn ) THEN                    ! Temporary saving of ua and va trends (l_trddyn)
+
+      IF( l_trddyn .OR. ln_wad ) THEN               ! Temporary saving of ua and va trends (l_trddyn)
          CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
          ztrdu(:,:,:) = ua(:,:,:)
          ztrdv(:,:,:) = va(:,:,:)
       ENDIF
+
       !
       SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation
@@ -95,4 +101,28 @@
       CASE (  4 )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme)
       END SELECT
+
+      !! Gravity filter for W/D
+      IF(ln_wad) THEN
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zij     =  sshn(ji   , jj   )
+                  zim1j   =  sshn(ji- 1, jj   )
+                  zijm1   =  sshn(ji   , jj- 1)
+                  zim1jm1 =  sshn(ji-1 , jj- 1)
+                  dij     = bathy(ji   , jj   )
+                  dim1j   = bathy(ji- 1, jj   )
+                  dijm1   = bathy(ji   , jj- 1)
+                  dim1jm1 = bathy(ji-1 , jj- 1)
+
+                  zcpx = 0.5_wp + sign(0.5_wp, min(zij,zim1j) - max(-dij,-dim1j))
+                  zcpy = 0.5_wp + sign(0.5_wp, min(zij,zijm1) - max(-dij,-dijm1))
+                  
+                  ua(ji,jj,jk) = ztrdu(ji,jj,jk) + zcpx * ( ua(ji,jj,jk) - ztrdu(ji,jj,jk))
+                  va(ji,jj,jk) = ztrdv(ji,jj,jk) + zcpy * ( va(ji,jj,jk) - ztrdv(ji,jj,jk))
+               END DO
+            END DO
+         END DO
+      END IF
       !
       IF( l_trddyn ) THEN      ! save the hydrostatic pressure gradient trends for momentum trend diagnostics
@@ -100,6 +130,10 @@
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
          CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt )
+      ENDIF
+
+      IF( l_trddyn .OR. ln_wad ) THEN 
          CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
       ENDIF
+
       !
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask,   &
Index: branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90
===================================================================
--- branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90	(revision 4356)
+++ branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90	(revision 4375)
@@ -49,10 +49,11 @@
       INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
       INTEGER  ::   zflag, z2dt         ! local scalar
-      REAL(wp) ::   zflxp, zflxn        ! local scalars
       REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
       REAL(wp) ::   ztmp                ! local scalars
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn   ! specific 2D workspace
       REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv   ! specific 2D workspace
       REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1  ! specific 2D workspace
-      REAL(wp), POINTER,  DIMENSION(:,:) ::   uwdlmt, vwdlmt  !: W/D flux limiters
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   wdlmt           !: W/D flux limiters
+
       !!----------------------------------------------------------------------
       !
@@ -62,5 +63,5 @@
       IF(ln_wad) THEN
 
-        CALL wrk_alloc( jpi, jpj, zflxu, zflxv, zflxu1, zflxv1, uwdlmt, vwdlmt)
+        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1, wdlmt )
         !
        
@@ -71,10 +72,9 @@
         IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
        
-       
+        zflxp(:,:)  = 0._wp
+        zflxn(:,:)  = 0._wp
         zflxu(:,:)  = 0._wp
         zflxv(:,:)  = 0._wp
-        uwdlmt(:,:) = 1._wp
-        vwdlmt(:,:) = 1._wp
-       
+        wdlmt(:,:)  = 1._wp
        
         ! Horizontal Flux in u and v direction
@@ -91,55 +91,67 @@
         zflxv(:,:) = zflxv(:,:) * e1v(:,:)
        
-        DO jk1 = 1, 10
-       
-        zflag = 0     ! flag marking that if we need further iteration of W/D limiters?
-       
-        zflxu1(:,:) = zflxu(:,:) * uwdlmt(:,:)
-        zflxv1(:,:) = zflxv(:,:) * vwdlmt(:,:)
-       
-       
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1   ! vector opt.
-              zflxp = 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) 
-              zflxn = 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) 
-       
+
               IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 
-       
-              !IF( .NOT. AGRIF_Root() ) THEN    !Wetting and Drying but not for AGRIF
-                 IF (((nbondi ==  1).OR.(nbondi == 2)).AND.(ji == nlci-1)) CYCLE ! east
-                 IF (((nbondi == -1).OR.(nbondi == 2)).AND.(ji == 2     )) CYCLE ! west 
-                 IF (((nbondj ==  1).OR.(nbondj == 2)).AND.(jj == nlcj-1)) CYCLE ! north
-                 IF (((nbondj == -1).OR.(nbondj == 2)).AND.(jj == 2     )) CYCLE ! south
-              !END IF
-       
-              
-              ztmp = e1t(ji,jj) * e2t(ji,jj)
-       
-              zdep1 = (zflxp + zflxn) * z2dt / ztmp
+
+              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) + sshb(ji,jj) - rn_wadmin 
-       
               IF(zdep2 < 0._wp) THEN
-                !WRITE(numout,*) 'depth less than minimum depth, cell(ji,jj):', ji,jj
                 zdep2 = 0._wp
                 sshb(ji,jj) = rn_wadmin - bathy(ji,jj)
               END IF
+           ENDDO
+        END DO
+
+      
+        !! start limiter iteration 
+        DO jk1 = 1, nn_waditr
        
-              IF(zdep1 > zdep2) THEN
-                zflag = 1
-                zcoef = ( ( zdep2 + rn_wadmine ) * ztmp + zflxn * z2dt ) / ( zflxp * z2dt )
-                IF(zflxu1(ji,  jj) >= 0._wp) uwdlmt(ji,jj  ) = zcoef
-                IF(zflxu1(ji-1,jj) <  0._wp) uwdlmt(ji-1,jj) = zcoef
-                IF(zflxv1(ji,  jj) >= 0._wp) vwdlmt(ji,jj  ) = zcoef
-                IF(zflxv1(ji,jj-1) <  0._wp) vwdlmt(ji,jj-1) = zcoef
-              END IF
-           END DO ! ji loop
-        END DO  ! jj loop
-       
-        CALL lbc_lnk( uwdlmt, 'U', 1. )   ;   CALL lbc_lnk( vwdlmt , 'V', 1. ) 
-       
-        IF(zflag == 0) EXIT
-       
+           zflag = 0     ! flag indicating if any further iteration is needed?
+          
+           zflxu1(:,:) = zflxu(:,:) * wdlmt(:,:)
+           zflxv1(:,:) = zflxv(:,:) * wdlmt(:,:)
+          
+           DO jj = 2, jpjm1
+              DO ji = fs_2, fs_jpim1   ! vector opt.
+        
+                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 
+        
+                 ztmp = e1t(ji,jj) * e2t(ji,jj)
+
+                 zflxn(ji,jj) = 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 = (zflxp(ji,jj) * wdlmt(ji,jj) + zflxn(ji,jj)) * z2dt / ztmp
+                 zdep2 = bathy(ji,jj) + sshb(ji,jj) - rn_wadmin 
+          
+                 IF(zdep1 > zdep2) THEN
+                   zflag = 1
+                   zcoef = ( ( zdep2 - rn_wadmine ) * ztmp - zflxn(ji,jj) * z2dt ) / ( zflxp(ji,jj) * z2dt )
+                   zcoef = max(zcoef, 0._wp)
+                   IF(zflxu1(ji,  jj  ) >= 0._wp) wdlmt(ji,  jj  ) = zcoef
+                   IF(zflxu1(ji-1,jj  ) <  0._wp) wdlmt(ji-1,jj  ) = zcoef
+                   IF(zflxv1(ji,  jj  ) >= 0._wp) wdlmt(ji,  jj  ) = zcoef
+                   IF(zflxv1(ji,  jj-1) <  0._wp) wdlmt(ji,  jj-1) = zcoef
+                 END IF
+              END DO ! ji loop
+           END DO  ! jj loop
+          
+           CALL lbc_lnk( wdlmt, 'T', 1. )
+          
+           IF(zflag == 0) EXIT
+          
+           IF(jk1 == nn_waditr) THEN
+              IF(zflxu1(ji,  jj  ) >= 0._wp) wdlmt(ji,  jj  ) = 0._wp
+              IF(zflxu1(ji-1,jj  ) <  0._wp) wdlmt(ji-1,jj  ) = 0._wp
+              IF(zflxv1(ji,  jj  ) >= 0._wp) wdlmt(ji,  jj  ) = 0._wp
+              IF(zflxv1(ji,  jj-1) <  0._wp) wdlmt(ji,  jj-1) = 0._wp
+           END IF
+
         END DO  ! jk1 loop
        
@@ -147,8 +159,8 @@
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
-               un(ji,jj,jk)   = 0.5_wp * ( sign(1._wp, un(ji,jj,jk))   + 1._wp ) * uwdlmt(ji,jj) 
-               vn(ji,jj,jk)   = 0.5_wp * ( sign(1._wp, vn(ji,jj,jk))   + 1._wp ) * vwdlmt(ji,jj) 
-               un(ji-1,jj,jk) = 0.5_wp * ( 1._wp - sign(1._wp, un(ji-1,jj,jk) )) * uwdlmt(ji-1,jj) 
-               vn(ji,jj-1,jk) = 0.5_wp * ( 1._wp - sign(1._wp, vn(ji,jj-1,jk) )) * vwdlmt(ji,jj-1) 
+               un(ji,  jj,  jk) = ( sign(0.5_wp, un(ji,  jj,  jk))  + 0.5_wp ) * wdlmt(ji  ,jj) 
+               vn(ji,  jj,  jk) = ( sign(0.5_wp, vn(ji,  jj,  jk))  + 0.5_wp ) * wdlmt(ji  ,jj) 
+               un(ji-1,jj,  jk) = (-sign(0.5_wp, un(ji-1,jj,  jk))  + 0.5_wp ) * wdlmt(ji-1,jj) 
+               vn(ji,  jj-1,jk) = (-sign(0.5_wp, vn(ji,  jj-1,jk))  + 0.5_wp ) * wdlmt(ji,  jj-1) 
               END DO
            END DO
@@ -161,5 +173,5 @@
         !
         !
-        CALL wrk_dealloc( jpi, jpj, zflxu, zflxv, zflxu1, zflxv1, uwdlmt, vwdlmt)
+        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1, wdlmt )
       !
       END IF
Index: branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/par_oce.F90
===================================================================
--- branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/par_oce.F90	(revision 4356)
+++ branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/par_oce.F90	(revision 4375)
@@ -86,4 +86,9 @@
    !!---------------------------------------------------------------------
 #             include "par_AMM_12km.h90"
+#elif defined key_pbasin
+   !!---------------------------------------------------------------------
+   !!   'key_pbasin':                             Parabolic Basin : PBASIN
+   !!---------------------------------------------------------------------
+#             include "par_PBASIN.h90"
 #else
    !!---------------------------------------------------------------------
Index: branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4356)
+++ branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4375)
@@ -87,9 +87,4 @@
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      !  Ocean dynamics : Wetting/Drying Limiter
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL wad_lmt( kstp )         ! wetting/drying limiter
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Update data, open boundaries, surface boundary condition (including sea-ice)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -100,4 +95,9 @@
       IF( lk_obc     )   CALL obc_rad( kstp )         ! compute phase velocities at open boundaries
       IF( lk_bdy     )   CALL bdy_dta( kstp, time_offset=+1 ) ! update dynamic and tracer data at open boundaries
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !  Ocean dynamics : Wetting/Drying Limiter
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         CALL wad_lmt( kstp )         ! wetting/drying limiter
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
