Index: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 6981)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 7646)
@@ -106,11 +106,11 @@
          WRITE(numout,*)
          WRITE(numout,*) 'dyn_adv_init : choice/control of the momentum advection scheme'
-         WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*) '       Namelist namdyn_adv : chose a advection formulation & scheme for momentum'
-         WRITE(numout,*) '          Vector/flux form (T/F)                           ln_dynadv_vec  = ', ln_dynadv_vec
-         WRITE(numout,*) '          = 0 standard scheme  ; =1 Hollingsworth scheme   nn_dynkeg      = ', nn_dynkeg
-         WRITE(numout,*) '          2nd order centred advection scheme               ln_dynadv_cen2 = ', ln_dynadv_cen2
-         WRITE(numout,*) '          3rd order UBS advection scheme                   ln_dynadv_ubs  = ', ln_dynadv_ubs
-         WRITE(numout,*) '          Sub timestepping of vertical advection           ln_dynzad_zts  = ', ln_dynzad_zts
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namdyn_adv : chose a advection formulation & scheme for momentum'
+         WRITE(numout,*) '      Vector/flux form (T/F)                           ln_dynadv_vec  = ', ln_dynadv_vec
+         WRITE(numout,*) '      = 0 standard scheme  ; =1 Hollingsworth scheme   nn_dynkeg      = ', nn_dynkeg
+         WRITE(numout,*) '      2nd order centred advection scheme               ln_dynadv_cen2 = ', ln_dynadv_cen2
+         WRITE(numout,*) '      3rd order UBS advection scheme                   ln_dynadv_ubs  = ', ln_dynadv_ubs
+         WRITE(numout,*) '      Sub timestepping of vertical advection           ln_dynzad_zts  = ', ln_dynzad_zts
       ENDIF
 
@@ -134,12 +134,12 @@
       IF(lwp) THEN                    ! Print the choice
          WRITE(numout,*)
-         IF( nadv ==  0 )   WRITE(numout,*) '         vector form : keg + zad + vor is used' 
-         IF( nadv ==  1 )   WRITE(numout,*) '         vector form : keg + zad_zts + vor is used'
+         IF( nadv ==  0 )   WRITE(numout,*) '      ===>>   vector form : keg + zad + vor is used' 
+         IF( nadv ==  1 )   WRITE(numout,*) '      ===>>   vector form : keg + zad_zts + vor is used'
          IF( nadv ==  0 .OR. nadv ==  1 ) THEN
-            IF( nn_dynkeg == nkeg_C2  )   WRITE(numout,*) 'with Centered standard keg scheme'
-            IF( nn_dynkeg == nkeg_HW  )   WRITE(numout,*) 'with Hollingsworth keg scheme'
+            IF( nn_dynkeg == nkeg_C2  )   WRITE(numout,*) '              with Centered standard keg scheme'
+            IF( nn_dynkeg == nkeg_HW  )   WRITE(numout,*) '              with Hollingsworth keg scheme'
          ENDIF
-         IF( nadv ==  2 )   WRITE(numout,*) '         flux form   : 2nd order scheme is used'
-         IF( nadv ==  3 )   WRITE(numout,*) '         flux form   : UBS       scheme is used'
+         IF( nadv ==  2 )   WRITE(numout,*) '      ===>>   flux form   : 2nd order scheme is used'
+         IF( nadv ==  3 )   WRITE(numout,*) '      ===>>   flux form   : UBS       scheme is used'
       ENDIF
       !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 6981)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 7646)
@@ -432,5 +432,5 @@
       INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices
       REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars
-      LOGICAL  ::   ll_tmp1, ll_tmp2, ll_tmp3            ! local logical variables
+      LOGICAL  ::   ll_tmp1, ll_tmp2                     ! local logical variables
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj
       REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy !W/D pressure filter
@@ -438,5 +438,5 @@
       !
       CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )
-      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
+      IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
       !
       IF( kt == nit000 ) THEN
@@ -451,47 +451,46 @@
       ENDIF
       !
-      IF(ln_wd) THEN
+      IF( ln_wd ) THEN
         DO jj = 2, jpjm1
            DO ji = 2, jpim1 
-             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
-             ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1 + rn_wdmin2
-             ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + &
-                                                       & rn_wdmin1 + rn_wdmin2
-
-             IF(ll_tmp1.AND.ll_tmp2) THEN
+             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
+                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            &
+                  &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &
+                  &                                                         > rn_wdmin1 + rn_wdmin2
+             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             &
+                  &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
+                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
+
+             IF(ll_tmp1) THEN
                zcpx(ji,jj) = 1.0_wp
-               wduflt(ji,jj) = 1.0_wp
-             ELSE IF(ll_tmp3) THEN
-               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here
-               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) / &
-                           &     (sshn(ji+1,jj) - sshn(ji,jj)))
-               wduflt(ji,jj) = 1.0_wp
+             ELSE 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) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &
+                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
              ELSE
                zcpx(ji,jj) = 0._wp
-               wduflt(ji,jj) = 0.0_wp
              END IF
       
-             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
-             ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1 + rn_wdmin2
-             ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + &
-                                                       & rn_wdmin1 + rn_wdmin2
-
-             IF(ll_tmp1.AND.ll_tmp2) THEN
+             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
+                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            &
+                  &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &
+                  &                                                         > rn_wdmin1 + rn_wdmin2
+             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             &
+                  &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
+                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
+
+             IF(ll_tmp1) THEN
                zcpy(ji,jj) = 1.0_wp
-               wdvflt(ji,jj) = 1.0_wp
-             ELSE IF(ll_tmp3) THEN
-               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here
-               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) / &
-                           &     (sshn(ji,jj+1) - sshn(ji,jj)))
-               wdvflt(ji,jj) = 1.0_wp
+             ELSE 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) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &
+                           &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) )
              ELSE
                zcpy(ji,jj) = 0._wp
-               wdvflt(ji,jj) = 0.0_wp
              END IF
            END DO
         END DO
         CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
-      ENDIF
-
+      END IF
 
       ! Surface value
@@ -510,5 +509,5 @@
 
 
-            IF(ln_wd) THEN
+            IF( ln_wd ) THEN
 
               zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
@@ -541,5 +540,5 @@
                   &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj)
 
-               IF(ln_wd) THEN
+               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) 
@@ -556,5 +555,5 @@
       !
       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
-      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
+      IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
       !
    END SUBROUTINE hpg_sco
@@ -701,38 +700,42 @@
       CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow )
       CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        )
-      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
-      !
-      !
-      IF(ln_wd) THEN
+      IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
+      !
+      !
+      IF( ln_wd ) THEN
         DO jj = 2, jpjm1
            DO ji = 2, jpim1 
-             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) &
-                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) &
-                     &  > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +&
-                                                       & rn_wdmin1 + rn_wdmin2
+             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
+                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            &
+                  &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &
+                  &                                                         > rn_wdmin1 + rn_wdmin2
+             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             &
+                  &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
+                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(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)))
+               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
+               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &
+                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
              ELSE
                zcpx(ji,jj) = 0._wp
              END IF
       
-             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) &
-                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) &
-                     &  > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +&
-                                                       & rn_wdmin1 + rn_wdmin2
+             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
+                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            &
+                  &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &
+                  &                                                         > rn_wdmin1 + rn_wdmin2
+             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             &
+                  &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
+                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(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)))
+               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
+               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &
+                           &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) )
              ELSE
                zcpy(ji,jj) = 0._wp
@@ -741,6 +744,5 @@
         END DO
         CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
-      ENDIF
-
+      END IF
 
       IF( kt == nit000 ) THEN
@@ -913,5 +915,5 @@
             zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj)
             zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj)
-            IF(ln_wd) THEN
+            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) 
@@ -936,5 +938,5 @@
                   &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    &
                   &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) * r1_e2v(ji,jj)
-               IF(ln_wd) THEN
+               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) 
@@ -950,5 +952,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 )
+      IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
       !
    END SUBROUTINE hpg_djc
@@ -987,5 +989,5 @@
       CALL wrk_alloc( jpi,jpj,jpk,   zdept, zrhh )
       CALL wrk_alloc( jpi,jpj,       zsshu_n, zsshv_n )
-      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
+      IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
       !
       IF( kt == nit000 ) THEN
@@ -1000,35 +1002,39 @@
       IF( ln_linssh )   znad = 0._wp
 
-      IF(ln_wd) THEN
+      IF( ln_wd ) THEN
         DO jj = 2, jpjm1
            DO ji = 2, jpim1 
-             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) &
-                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) &
-                     &  > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +&
-                                                       & rn_wdmin1 + rn_wdmin2
+             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
+                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            &
+                  &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &
+                  &                                                         > rn_wdmin1 + rn_wdmin2
+             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             &
+                  &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
+                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(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)))
+               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
+               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &
+                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
              ELSE
                zcpx(ji,jj) = 0._wp
              END IF
       
-             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) &
-                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) &
-                     &  > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +&
-                                                       & rn_wdmin1 + rn_wdmin2
-
-             IF(ll_tmp1.OR.ll_tmp2) THEN
+             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
+                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            &
+                  &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &
+                  &                                                         > rn_wdmin1 + rn_wdmin2
+             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             &
+                  &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
+                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(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)))
+               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
+               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &
+                           &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) )
              ELSE
                zcpy(ji,jj) = 0._wp
@@ -1037,5 +1043,5 @@
         END DO
         CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
-      ENDIF
+      END IF
 
       ! Clean 3-D work arrays
@@ -1046,5 +1052,5 @@
       DO jj = 1, jpj
         DO ji = 1, jpi
-          jk = mbathy(ji,jj)
+          jk = mbkt(ji,jj)+1
           IF(     jk <=  0   ) THEN   ;   zrhh(ji,jj,    :   ) = 0._wp
           ELSEIF( jk ==  1   ) THEN   ;   zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk)
@@ -1221,5 +1227,5 @@
                  zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)
                ENDIF
-               IF(ln_wd) THEN
+               IF( ln_wd ) THEN
                   zdpdx1 = zdpdx1 * zcpx(ji,jj)
                   zdpdx2 = zdpdx2 * zcpx(ji,jj)
@@ -1280,5 +1286,5 @@
                   zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )
                ENDIF
-               IF(ln_wd) THEN
+               IF( ln_wd ) THEN
                   zdpdy1 = zdpdy1 * zcpy(ji,jj)
                   zdpdy2 = zdpdy2 * zcpy(ji,jj)
@@ -1295,5 +1301,5 @@
       CALL wrk_dealloc( jpi,jpj,jpk,   zdept, zrhh )
       CALL wrk_dealloc( jpi,jpj,       zsshu_n, zsshv_n )
-      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
+      IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
       !
    END SUBROUTINE hpg_prj
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 6981)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 7646)
@@ -24,4 +24,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE bdy_oce         ! ocean open boundary conditions
 
    IMPLICIT NONE
@@ -78,4 +79,7 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke
       REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 
+      INTEGER  ::   jb                 ! dummy loop indices
+      INTEGER  ::   ii, ij, igrd, ib_bdy   ! local integers
+      INTEGER  ::   fu, fv
       !!----------------------------------------------------------------------
       !
@@ -98,4 +102,32 @@
       zhke(:,:,jpk) = 0._wp
       
+      IF (ln_bdy) THEN
+         ! Maria Luneva & Fred Wobus: July-2016
+         ! compensate for lack of turbulent kinetic energy on liquid bdy points
+         DO ib_bdy = 1, nb_bdy
+            IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
+               igrd = 2           ! Copying normal velocity into points outside bdy
+               DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+                  DO jk = 1, jpkm1
+                     ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
+                     ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
+                     fu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) )
+                     un(ii-fu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk)
+                  END DO
+               END DO
+               !
+               igrd = 3           ! Copying normal velocity into points outside bdy
+               DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+                  DO jk = 1, jpkm1
+                     ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
+                     ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
+                     fv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) )
+                     vn(ii,ij-fv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk)
+                  END DO
+               END DO
+            ENDIF
+         ENDDO  
+      ENDIF 
+
       SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
       !
@@ -133,4 +165,12 @@
          !
       END SELECT
+
+      IF (ln_bdy) THEN
+         ! restore velocity masks at points outside boundary
+         un(:,:,:) = un(:,:,:) * umask(:,:,:)
+         vn(:,:,:) = vn(:,:,:) * vmask(:,:,:)
+      ENDIF      
+
+
       !
       DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 6981)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 7646)
@@ -110,11 +110,11 @@
          WRITE(numout,*)
          WRITE(numout,*) 'dyn_ldf_init : Choice of the lateral diffusive operator on dynamics'
-         WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*) '       Namelist nam_dynldf : set lateral mixing parameters (type, direction, coefficients)'
-         WRITE(numout,*) '          laplacian operator          ln_dynldf_lap = ', ln_dynldf_lap
-         WRITE(numout,*) '          bilaplacian operator        ln_dynldf_blp = ', ln_dynldf_blp
-         WRITE(numout,*) '          iso-level                   ln_dynldf_lev = ', ln_dynldf_lev
-         WRITE(numout,*) '          horizontal (geopotential)   ln_dynldf_hor = ', ln_dynldf_hor
-         WRITE(numout,*) '          iso-neutral                 ln_dynldf_iso = ', ln_dynldf_iso
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist nam_dynldf : set lateral mixing parameters (type, direction, coefficients)'
+         WRITE(numout,*) '      laplacian operator          ln_dynldf_lap = ', ln_dynldf_lap
+         WRITE(numout,*) '      bilaplacian operator        ln_dynldf_blp = ', ln_dynldf_blp
+         WRITE(numout,*) '      iso-level                   ln_dynldf_lev = ', ln_dynldf_lev
+         WRITE(numout,*) '      horizontal (geopotential)   ln_dynldf_hor = ', ln_dynldf_hor
+         WRITE(numout,*) '      iso-neutral                 ln_dynldf_iso = ', ln_dynldf_iso
       ENDIF
       !                                   ! use of lateral operator or not
@@ -180,8 +180,8 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         IF( nldf == np_no_ldf )   WRITE(numout,*) '              NO lateral viscosity'
-         IF( nldf == np_lap    )   WRITE(numout,*) '              iso-level laplacian operator'
-         IF( nldf == np_lap_i  )   WRITE(numout,*) '              rotated laplacian operator with iso-level background'
-         IF( nldf == np_blp    )   WRITE(numout,*) '              iso-level bi-laplacian operator'
+         IF( nldf == np_no_ldf )   WRITE(numout,*) '      ===>>   NO lateral viscosity'
+         IF( nldf == np_lap    )   WRITE(numout,*) '      ===>>   iso-level laplacian operator'
+         IF( nldf == np_lap_i  )   WRITE(numout,*) '      ===>>   rotated laplacian operator with iso-level background'
+         IF( nldf == np_blp    )   WRITE(numout,*) '      ===>>   iso-level bi-laplacian operator'
       ENDIF
       !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 6981)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 7646)
@@ -32,5 +32,5 @@
    USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
    USE domvvl         ! variable volume
-   USE bdy_oce        ! ocean open boundary conditions
+   USE bdy_oce   , ONLY: ln_bdy
    USE bdydta         ! ocean open boundary conditions
    USE bdydyn         ! ocean open boundary conditions
@@ -77,5 +77,5 @@
       !!              * Apply lateral boundary conditions on after velocity 
       !!             at the local domain boundaries through lbc_lnk call,
-      !!             at the one-way open boundaries (lk_bdy=T),
+      !!             at the one-way open boundaries (ln_bdy=T),
       !!             at the AGRIF zoom   boundaries (lk_agrif=T)
       !!
@@ -147,12 +147,9 @@
       CALL lbc_lnk( va, 'V', -1. ) 
       !
-# if defined key_bdy
       !                                !* BDY open boundaries
-      IF( lk_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
-      IF( lk_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )
+      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
+      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )
 
 !!$   Do we need a call to bdy_vol here??
-      !
-# endif
       !
       IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 6981)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 7646)
@@ -88,5 +88,5 @@
       !
       IF(      ln_apr_dyn                                                &   ! atmos. pressure
-         .OR.  ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) )   &   ! tide potential (no time slitting)
+         .OR.  ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. ln_tide) )   &   ! tide potential (no time slitting)
          .OR.  nn_ice_embd == 2  ) THEN                                      ! embedded sea-ice
          !
@@ -111,5 +111,5 @@
          !
          !                                    !==  tide potential forcing term  ==!
-         IF( .NOT.ln_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide )  ) THEN   ! N.B. added directly at sub-time-step in ts-case
+         IF( .NOT.ln_dynspg_ts .AND. ( ln_tide_pot .AND. ln_tide )  ) THEN   ! N.B. added directly at sub-time-step in ts-case
             !
             CALL upd_tide( kt )                      ! update tide potential
@@ -216,7 +216,7 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         IF( nspg == np_EXP )   WRITE(numout,*) '     explicit free surface'
-         IF( nspg == np_TS  )   WRITE(numout,*) '     free surface with time splitting scheme'
-         IF( nspg == np_NO  )   WRITE(numout,*) '     No surface surface pressure gradient trend in momentum Eqs.'
+         IF( nspg == np_EXP )   WRITE(numout,*) '      ===>>   explicit free surface'
+         IF( nspg == np_TS  )   WRITE(numout,*) '      ===>>   free surface with time splitting scheme'
+         IF( nspg == np_NO  )   WRITE(numout,*) '      ===>>   No surface surface pressure gradient trend in momentum Eqs.'
       ENDIF
       !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 6981)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 7646)
@@ -15,4 +15,5 @@
    !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility
    !!             3.7  ! 2015-11  (J. Chanut) free surface simplification
+   !!              -   ! 2016-12  (G. Madec, E. Clementi) update for Stoke-Drift divergence
    !!---------------------------------------------------------------------
 
@@ -33,9 +34,10 @@
    USE dynvor          ! vorticity term
    USE wet_dry         ! wetting/drying flux limter
-   USE bdy_par         ! for lk_bdy
+   USE bdy_oce         ! open boundary
    USE bdytides        ! open boundary condition data
    USE bdydyn2d        ! open boundary conditions on barotropic variables
    USE sbctide         ! tides
    USE updtide         ! tide potential
+   USE sbcwave         ! surface wave
    !
    USE in_out_manager  ! I/O manager
@@ -69,5 +71,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   wgtbtp1, wgtbtp2   !: 1st & 2nd weights used in time filtering of barotropic fields
 
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz          !: ff/h at F points
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz          !: ff_f/h at F points
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftnw, ftne   !: triad of coriolis parameter
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse   !: (only used with een vorticity scheme)
@@ -156,5 +158,4 @@
       REAL(wp), POINTER, DIMENSION(:,:) :: zhf
       REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef.
-      REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1           ! Wetting/Dying velocity filter coef.
       !!----------------------------------------------------------------------
       !
@@ -168,5 +169,5 @@
       CALL wrk_alloc( jpi,jpj,   zsshu_a, zsshv_a                  )
       CALL wrk_alloc( jpi,jpj,   zhf )
-      IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 )
+      IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy )
       !
       zmdi=1.e+20                               !  missing data indicator for masking
@@ -220,5 +221,5 @@
       IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN
          IF( ln_dynvor_een ) THEN               !==  EEN scheme  ==!
-            SELECT CASE( nn_een_e3f )              !* ff/e3 at F-point
+            SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point
             CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
                DO jj = 1, jpjm1
@@ -226,5 +227,5 @@
                      zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                    &
                         &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   ) * 0.25_wp  
-                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj)
+                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
                   END DO
                END DO
@@ -236,5 +237,5 @@
                         &       / ( MAX( 1._wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    &
                         &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) )
-                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj)
+                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
                   END DO
                END DO
@@ -255,21 +256,41 @@
             zwz(:,:) = 0._wp
             zhf(:,:) = 0._wp
-            IF ( .not. ln_sco ) THEN
-
-!!gm  agree the JC comment  : this should be done in a much clear way
-
-! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
-!     Set it to zero for the time being 
-!              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
-!              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
-!              ENDIF
-!              zhf(:,:) = gdepw_0(:,:,jk+1)
-            ELSE
-               zhf(:,:) = hbatf(:,:)
-            END IF
-
-            DO jj = 1, jpjm1
-               zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
-            END DO
+            
+!!gm  assume 0 in both cases (xhich is almost surely WRONG ! ) as hvatf has been removed 
+!!gm    A priori a better value should be something like :
+!!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1) 
+!!gm                     divided by the sum of the corresponding mask 
+!!gm 
+!!            
+              IF ( .not. ln_sco ) THEN
+  
+   !!gm  agree the JC comment  : this should be done in a much clear way
+  
+   ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
+   !     Set it to zero for the time being 
+   !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
+   !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
+   !              ENDIF
+   !              zhf(:,:) = gdepw_0(:,:,jk+1)
+               ELSE
+                 !zhf(:,:) = hbatf(:,:)
+                 DO jj = 1, jpjm1
+                   DO ji = 1, jpim1
+                     zhf(ji,jj) = MAX( 0._wp,                                &
+                                & ( ht_0(ji  ,jj  )*tmask(ji  ,jj  ,1) +     &
+                                &   ht_0(ji+1,jj  )*tmask(ji+1,jj  ,1) +     &
+                                &   ht_0(ji  ,jj+1)*tmask(ji  ,jj+1,1) +     &
+                                &   ht_0(ji+1,jj+1)*tmask(ji+1,jj+1,1) ) /   &
+                                & ( tmask(ji  ,jj  ,1) + tmask(ji+1,jj  ,1) +&
+                                &   tmask(ji  ,jj+1,1) + tmask(ji+1,jj+1,1) +&
+                                &   rsmall  ) )
+                   END DO
+                 END DO
+              END IF
+  
+              DO jj = 1, jpjm1
+                 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
+              END DO
+!!gm end
 
             DO jk = 1, jpkm1
@@ -285,5 +306,5 @@
                END DO
             END DO
-            zwz(:,:) = ff(:,:) * zwz(:,:)
+            zwz(:,:) = ff_f(:,:) * zwz(:,:)
          ENDIF
       ENDIF
@@ -324,4 +345,9 @@
          END DO
       END DO
+      
+!!gm  Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....
+!!gm  Is it correct to do so ?   I think so...
+      
+      
       !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
       !                                   ! --------------------------------------------------------
@@ -374,51 +400,50 @@
       IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient
         IF( ln_wd ) THEN                        ! Calculating and applying W/D gravity filters
-          wduflt1(:,:) = 1.0_wp
-          wdvflt1(:,:) = 1.0_wp
-          DO jj = 2, jpjm1
-             DO ji = 2, jpim1
-                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))   &
-                        & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj))   &
-                        &  > rn_wdmin1 + rn_wdmin2
-                ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))   &
-                        &                                   + rn_wdmin1 + rn_wdmin2
+           DO jj = 2, jpjm1
+              DO ji = 2, jpim1 
+                ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
+                     &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            &
+                     &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &
+                     &                                                         > rn_wdmin1 + rn_wdmin2
+                ll_tmp2 = ( ABS( sshn(ji+1,jj)             -   sshn(ji  ,jj))  > 1.E-12 ).AND.( &
+                     &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
+                     &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
+   
                 IF(ll_tmp1) THEN
-                  zcpx(ji,jj)    = 1.0_wp
-                ELSEIF(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)))
+                  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) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &
+                              &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
                 ELSE
-                  zcpx(ji,jj)    = 0._wp
-                  wduflt1(ji,jj) = 0.0_wp
+                  zcpx(ji,jj) = 0._wp
                 END IF
-
-                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))   &
-                        & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1))   &
-                        &  > rn_wdmin1 + rn_wdmin2
-                ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))   &
-                        &                                   + rn_wdmin1 + rn_wdmin2
+         
+                ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
+                     &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            &
+                     &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &
+                     &                                                         > rn_wdmin1 + rn_wdmin2
+                ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1))  > 1.E-12 ).AND.( &
+                     &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
+                     &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
+   
                 IF(ll_tmp1) THEN
-                   zcpy(ji,jj)    = 1.0_wp
-                ELSEIF(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)))
+                  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) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &
+                              &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) )
                 ELSE
-                  zcpy(ji,jj)    = 0._wp
-                  wdvflt1(ji,jj) = 0.0_wp
-                ENDIF
-
-             END DO
+                  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 )
-
+ 
            DO jj = 2, jpjm1
               DO ji = 2, jpim1
-                 zu_trd(ji,jj) = ( zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   &
-                        &                        * r1_e1u(ji,jj) ) * zcpx(ji,jj) * wduflt1(ji,jj)
-                 zv_trd(ji,jj) = ( zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   &
-                        &                        * r1_e2v(ji,jj) ) * zcpy(ji,jj) * wdvflt1(ji,jj)
+                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   &
+                        &                        * r1_e1u(ji,jj) * zcpx(ji,jj)
+                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   &
+                        &                        * r1_e2v(ji,jj) * zcpy(ji,jj)
               END DO
            END DO
@@ -474,5 +499,5 @@
       !
       !                                         ! Add top stress contribution from baroclinic velocities:      
-      IF (ln_bt_fw) THEN
+      IF( ln_bt_fw ) THEN
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -538,4 +563,9 @@
                 &                        + fwfisf(:,:) + fwfisf_b(:,:)                     )
       ENDIF
+      !
+      IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary
+         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)
+      ENDIF
+      !
 #if defined key_asminc
       !                                         ! Include the IAU weighted SSH increment
@@ -567,9 +597,4 @@
       ENDIF
 
-      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                    
@@ -608,6 +633,6 @@
          ! Update only tidal forcing at open boundaries
 #if defined key_tide
-         IF( lk_bdy      .AND. lk_tide )   CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )
-         IF( ln_tide_pot .AND. lk_tide )   CALL upd_tide     ( kt, kit=jn, time_offset= noffset   )
+         IF( ln_bdy      .AND. ln_tide )   CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )
+         IF( ln_tide_pot .AND. ln_tide )   CALL upd_tide     ( kt, kit=jn, time_offset= noffset   )
 #endif
          !
@@ -646,8 +671,4 @@
             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_n(:,:)
@@ -702,11 +723,9 @@
          END DO
          ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:)
-         IF( ln_wd ) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:)) 
+         
          CALL lbc_lnk( ssha_e, 'T',  1._wp )
 
-#if defined key_bdy
          ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
-         IF( lk_bdy )   CALL bdy_ssh( ssha_e )
-#endif
+         IF( ln_bdy )   CALL bdy_ssh( ssha_e )
 #if defined key_agrif
          IF( .NOT.Agrif_Root() )   CALL agrif_ssh_ts( jn )
@@ -750,43 +769,44 @@
           &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
          IF( ln_wd ) THEN                   ! Calculating and applying W/D gravity filters
-           wduflt1(:,:) = 1._wp
-           wdvflt1(:,:) = 1._wp
            DO jj = 2, jpjm1
-              DO ji = 2, jpim1
-                 ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) &
-                        & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj) )    &
-                        &                                  > rn_wdmin1 + rn_wdmin2
-                 ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) &
-                        &                                  + rn_wdmin1 + rn_wdmin2
-                 IF(ll_tmp1) THEN
-                    zcpx(ji,jj) = 1._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
-                    wduflt1(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) ) &
-                        & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1) )    &
-                        &                                  > rn_wdmin1 + rn_wdmin2
-                 ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) &
-                        &                                  + rn_wdmin1 + rn_wdmin2
-                 IF(ll_tmp1) THEN
-                    zcpy(ji,jj) = 1._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
-                    wdvflt1(ji,jj) = 0._wp
-                 END IF
+              DO ji = 2, jpim1 
+                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                &
+                     &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji+1,jj) ) .AND.            &
+                     &    MAX( zsshp2_e(ji,jj) + ht_wd(ji,jj), zsshp2_e(ji+1,jj) + ht_wd(ji+1,jj) ) &
+                     &                                                             > rn_wdmin1 + rn_wdmin2
+                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji+1,jj))  > 1.E-12 ).AND.( &
+                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                &
+                     &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(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) +    ht_wd(ji+1,jj) - zsshp2_e(ji,jj) - ht_wd(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(   -ht_wd(ji,jj)               ,   -ht_wd(ji,jj+1) ) .AND.            &
+                     &    MAX( zsshp2_e(ji,jj) + ht_wd(ji,jj), zsshp2_e(ji,jj+1) + ht_wd(ji,jj+1) ) &
+                     &                                                             > rn_wdmin1 + rn_wdmin2
+                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji,jj+1))  > 1.E-12 ).AND.( &
+                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                &
+                     &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(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) +    ht_wd(ji,jj+1) - zsshp2_e(ji,jj) - ht_wd(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
+           END DO
+         END IF
          !
          ! Compute associated depths at U and V points:
@@ -806,9 +826,4 @@
             END DO
 
-            IF( ln_wd ) THEN
-              zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 )
-              zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 )
-            END IF
-
          ENDIF
          !
@@ -861,5 +876,5 @@
          !
          ! Add tidal astronomical forcing if defined
-         IF ( lk_tide.AND.ln_tide_pot ) THEN
+         IF ( ln_tide .AND. ln_tide_pot ) THEN
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -967,8 +982,6 @@
          CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
          !
-#if defined key_bdy  
          !                                                 ! open boundaries
-         IF( lk_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
-#endif
+         IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
 #if defined key_agrif                                                           
          IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
@@ -1086,5 +1099,5 @@
       CALL wrk_dealloc( jpi,jpj,   zsshu_a, zsshv_a                                   )
       CALL wrk_dealloc( jpi,jpj,   zhf )
-      IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 )
+      IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy )
       !
       IF ( ln_diatmb ) THEN
@@ -1248,5 +1261,5 @@
             zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
             zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
-            zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) )
+            zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
          END DO
       END DO
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 6981)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 7646)
@@ -17,4 +17,5 @@
    !!            3.7  ! 2014-04  (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity 
    !!             -   ! 2014-06  (G. Madec) suppression of velocity curl from in-core memory
+   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T)
    !!----------------------------------------------------------------------
 
@@ -32,4 +33,6 @@
    USE trd_oce        ! trends: ocean variables
    USE trddyn         ! trend manager: dynamics
+   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)
+   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force
    !
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
@@ -77,5 +80,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2016)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -108,5 +111,5 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ene( kt, nrvm, ua, va )                      ! relative vorticity or metric trend
+            CALL vor_ene( kt, nrvm, un , vn , ua, va )                    ! relative vorticity or metric trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -114,10 +117,11 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ene( kt, ncor, ua, va )                      ! planetary vorticity trend
+            CALL vor_ene( kt, ncor, un , vn , ua, va )                    ! planetary vorticity trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
-         ELSE
-            CALL vor_ene( kt, ntot, ua, va )                ! total vorticity trend
+         ELSE                                               ! total vorticity trend
+                             CALL vor_ene( kt, ntot, un , vn , ua, va )   ! total vorticity trend
+            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
          ENDIF
          !
@@ -126,5 +130,5 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ens( kt, nrvm, ua, va )                      ! relative vorticity or metric trend
+            CALL vor_ens( kt, nrvm, un , vn , ua, va )            ! relative vorticity or metric trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -132,10 +136,11 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ens( kt, ncor, ua, va )                      ! planetary vorticity trend
+            CALL vor_ens( kt, ncor, un , vn , ua, va )            ! planetary vorticity trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
-         ELSE
-            CALL vor_ens( kt, ntot, ua, va )                ! total vorticity trend
+         ELSE                                               ! total vorticity trend
+                             CALL vor_ens( kt, ntot, un , vn , ua, va )  ! total vorticity trend
+            IF( ln_stcor )   CALL vor_ens( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend
          ENDIF
          !
@@ -144,5 +149,5 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ens( kt, nrvm, ua, va )                      ! relative vorticity or metric trend (ens)
+            CALL vor_ens( kt, nrvm, un , vn , ua, va )            ! relative vorticity or metric trend (ens)
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -150,11 +155,12 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ene( kt, ncor, ua, va )                      ! planetary vorticity trend (ene)
+            CALL vor_ene( kt, ncor, un , vn , ua, va )            ! planetary vorticity trend (ene)
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
-         ELSE
-            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend (ens)
-            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend (ene)
+         ELSE                                               ! total vorticity trend
+                             CALL vor_ens( kt, nrvm, un , vn , ua, va )   ! relative vorticity or metric trend (ens)
+                             CALL vor_ene( kt, ncor, un , vn , ua, va )   ! planetary vorticity trend (ene)
+            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
         ENDIF
          !
@@ -163,5 +169,5 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_een( kt, nrvm, ua, va )                      ! relative vorticity or metric trend
+            CALL vor_een( kt, nrvm, un , vn , ua, va )            ! relative vorticity or metric trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -169,10 +175,11 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_een( kt, ncor, ua, va )                      ! planetary vorticity trend
+            CALL vor_een( kt, ncor, un , vn , ua, va )            ! planetary vorticity trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
-         ELSE
-            CALL vor_een( kt, ntot, ua, va )                ! total vorticity trend
+         ELSE                                               ! total vorticity trend
+                             CALL vor_een( kt, ntot, un , vn , ua, va )   ! total vorticity trend
+            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
          ENDIF
          !
@@ -190,5 +197,5 @@
 
 
-   SUBROUTINE vor_ene( kt, kvor, pua, pva )
+   SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE vor_ene  ***
@@ -210,9 +217,9 @@
       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
       !!----------------------------------------------------------------------
-      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
-      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
-      !                                                           ! =nrvm (relative vorticity or metric)
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
+      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index
+      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;
+      !                                                                ! =nrvm (relative vorticity or metric)
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua, pva    ! total v-trend
       !
       INTEGER  ::   ji, jj, jk           ! dummy loop indices
@@ -223,5 +230,5 @@
       IF( nn_timing == 1 )  CALL timing_start('vor_ene')
       !
-      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz ) 
+      CALL wrk_alloc( jpi,jpj,   zwx, zwy, zwz ) 
       !
       IF( kt == nit000 ) THEN
@@ -237,10 +244,10 @@
          SELECT CASE( kvor )                 !==  vorticity considered  ==!
          CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            zwz(:,:) = ff(:,:) 
+            zwz(:,:) = ff_f(:,:) 
          CASE ( np_RVO )                           !* relative vorticity
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
-                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
+                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
+                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
                END DO
             END DO
@@ -248,6 +255,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
+                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
                        &     * 0.5 * r1_e1e2f(ji,jj)
                END DO
@@ -256,6 +263,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
-                     &                      - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
+                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
+                     &                      - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &
                      &                   * r1_e1e2f(ji,jj)
                END DO
@@ -264,7 +271,7 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ff(ji,jj)                                                                        &
-                       &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                  zwz(ji,jj) = ff_f(ji,jj)                                                                        &
+                       &     + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
+                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
                        &     * 0.5 * r1_e1e2f(ji,jj)
                END DO
@@ -284,9 +291,9 @@
          IF( ln_sco ) THEN
             zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
-            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
-            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
+            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
+            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
          ELSE
-            zwx(:,:) = e2u(:,:) * un(:,:,jk)
-            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
+            zwx(:,:) = e2u(:,:) * pun(:,:,jk)
+            zwy(:,:) = e1v(:,:) * pvn(:,:,jk)
          ENDIF
          !                                   !==  compute and add the vorticity term trend  =!
@@ -311,5 +318,5 @@
 
 
-   SUBROUTINE vor_ens( kt, kvor, pua, pva )
+   SUBROUTINE vor_ens( kt, kvor, pun, pvn, pua, pva )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE vor_ens  ***
@@ -331,9 +338,9 @@
       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
       !!----------------------------------------------------------------------
-      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
-      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
-         !                                                        ! =nrvm (relative vorticity or metric)
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
+      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index
+      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;
+         !                                                             ! =nrvm (relative vorticity or metric)
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua, pva    ! total v-trend
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -344,5 +351,5 @@
       IF( nn_timing == 1 )  CALL timing_start('vor_ens')
       !
-      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz ) 
+      CALL wrk_alloc( jpi,jpj,   zwx, zwy, zwz ) 
       !
       IF( kt == nit000 ) THEN
@@ -357,10 +364,10 @@
          SELECT CASE( kvor )                 !==  vorticity considered  ==!
          CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            zwz(:,:) = ff(:,:) 
+            zwz(:,:) = ff_f(:,:) 
          CASE ( np_RVO )                           !* relative vorticity
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
-                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
+                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
+                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
                END DO
             END DO
@@ -368,6 +375,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
+                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
                        &     * 0.5 * r1_e1e2f(ji,jj)
                END DO
@@ -376,6 +383,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
-                     &                      - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
+                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
+                     &                      - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &
                      &                   * r1_e1e2f(ji,jj)
                END DO
@@ -384,7 +391,7 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ff(ji,jj)                                                                       &
-                       &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                  zwz(ji,jj) = ff_f(ji,jj)                                                                       &
+                       &     + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
+                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
                        &     * 0.5 * r1_e1e2f(ji,jj)
                END DO
@@ -404,9 +411,9 @@
          IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
             zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
-            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
-            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
+            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
+            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
          ELSE
-            zwx(:,:) = e2u(:,:) * un(:,:,jk)
-            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
+            zwx(:,:) = e2u(:,:) * pun(:,:,jk)
+            zwy(:,:) = e1v(:,:) * pvn(:,:,jk)
          ENDIF
          !                                   !==  compute and add the vorticity term trend  =!
@@ -431,5 +438,5 @@
 
 
-   SUBROUTINE vor_een( kt, kvor, pua, pva )
+   SUBROUTINE vor_een( kt, kvor, pun, pvn, pua, pva )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE vor_een  ***
@@ -448,8 +455,9 @@
       !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
       !!----------------------------------------------------------------------
-      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
-      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; =nrvm (relative or metric)
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
+      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index
+      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;
+         !                                                             ! =nrvm (relative vorticity or metric)
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua, pva    ! total v-trend
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -506,5 +514,5 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ff(ji,jj) * z1_e3f(ji,jj)
+                  zwz(ji,jj) = ff_f(ji,jj) * z1_e3f(ji,jj)
                END DO
             END DO
@@ -512,6 +520,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
-                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
+                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
+                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &
                      &       * r1_e1e2f(ji,jj) * z1_e3f(ji,jj)
                END DO
@@ -520,6 +528,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
+                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
                        &     * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj)
                END DO
@@ -528,6 +536,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (  ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
-                     &                         - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
+                  zwz(ji,jj) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
+                     &                           - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &
                      &                      * r1_e1e2f(ji,jj)    ) * z1_e3f(ji,jj)
                END DO
@@ -536,7 +544,7 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (  ff(ji,jj)                                                                        &
-                       &        + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &            - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                  zwz(ji,jj) = (  ff_f(ji,jj)                                                                        &
+                       &        + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
+                       &            - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
                        &        * 0.5 * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
                END DO
@@ -557,6 +565,6 @@
          !
          !                                   !==  horizontal fluxes  ==!
-         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
-         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
+         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
+         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
 
          !                                   !==  compute and add the vorticity term trend  =!
@@ -626,11 +634,11 @@
          WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '        Namelist namdyn_vor : choice of the vorticity term scheme'
-         WRITE(numout,*) '           energy    conserving scheme                    ln_dynvor_ene = ', ln_dynvor_ene
-         WRITE(numout,*) '           enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
-         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
-         WRITE(numout,*) '           enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
-         WRITE(numout,*) '              e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
-         WRITE(numout,*) '           masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
+         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme'
+         WRITE(numout,*) '      energy    conserving scheme                    ln_dynvor_ene = ', ln_dynvor_ene
+         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
+         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
+         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
+         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
+         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
       ENDIF
 
@@ -639,5 +647,5 @@
       ! at angles with three ocean points and one land point
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '           namlbc: change fmask value in the angles (T)   ln_vorlat = ', ln_vorlat
+      IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
       IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
          DO jk = 1, jpk
@@ -666,9 +674,9 @@
       ncor = np_COR
       IF( ln_dynadv_vec ) THEN     
-         IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity'
+         IF(lwp) WRITE(numout,*) '      ===>>   Vector form advection : vorticity = Coriolis + relative vorticity'
          nrvm = np_RVO        ! relative vorticity
          ntot = np_CRV        ! relative + planetary vorticity
       ELSE                        
-         IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term'
+         IF(lwp) WRITE(numout,*) '      ===>>   Flux form advection   : vorticity = Coriolis + metric term'
          nrvm = np_MET        ! metric term
          ntot = np_CME        ! Coriolis + metric term
@@ -677,8 +685,8 @@
       IF(lwp) THEN                   ! Print the choice
          WRITE(numout,*)
-         IF( nvor_scheme ==  np_ENE )   WRITE(numout,*) '         vorticity scheme ==>> energy conserving scheme'
-         IF( nvor_scheme ==  np_ENS )   WRITE(numout,*) '         vorticity scheme ==>> enstrophy conserving scheme'
-         IF( nvor_scheme ==  np_MIX )   WRITE(numout,*) '         vorticity scheme ==>> mixed enstrophy/energy conserving scheme'
-         IF( nvor_scheme ==  np_EEN )   WRITE(numout,*) '         vorticity scheme ==>> energy and enstrophy conserving scheme'
+         IF( nvor_scheme ==  np_ENE )   WRITE(numout,*) '      ===>>   energy conserving scheme'
+         IF( nvor_scheme ==  np_ENS )   WRITE(numout,*) '      ===>>   enstrophy conserving scheme'
+         IF( nvor_scheme ==  np_MIX )   WRITE(numout,*) '      ===>>   mixed enstrophy/energy conserving scheme'
+         IF( nvor_scheme ==  np_EEN )   WRITE(numout,*) '      ===>>   energy and enstrophy conserving scheme'
       ENDIF
       !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 6981)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 7646)
@@ -119,6 +119,6 @@
          WRITE(numout,*) 'dyn_zdf_init : vertical dynamics physics scheme'
          WRITE(numout,*) '~~~~~~~~~~~'
-         IF( nzdf ==  0 )   WRITE(numout,*) '              Explicit time-splitting scheme'
-         IF( nzdf ==  1 )   WRITE(numout,*) '              Implicit (euler backward) scheme'
+         IF( nzdf ==  0 )   WRITE(numout,*) '      ===>>   Explicit time-splitting scheme'
+         IF( nzdf ==  1 )   WRITE(numout,*) '      ===>>   Implicit (euler backward) scheme'
       ENDIF
       !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 6981)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 7646)
@@ -22,6 +22,5 @@
    USE divhor         ! horizontal divergence
    USE phycst         ! physical constants
-   USE bdy_oce        ! 
-   USE bdy_par        !
+   USE bdy_oce   , ONLY: ln_bdy, bdytmask
    USE bdydyn2d       ! bdy_ssh routine
 #if defined key_agrif
@@ -88,12 +87,17 @@
       ENDIF
       !
-      CALL div_hor( kt )                              ! Horizontal divergence
-      !
-      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
+      z2dt = 2._wp * rdt                          ! set time step size (Euler/Leapfrog)
       IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
+      zcoef = 0.5_wp * r1_rau0
 
       !                                           !------------------------------!
       !                                           !   After Sea Surface Height   !
       !                                           !------------------------------!
+      IF(ln_wd) THEN
+         CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt)
+      ENDIF
+
+      CALL div_hor( kt )                               ! Horizontal divergence
+      !
       zhdiv(:,:) = 0._wp
       DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
@@ -104,8 +108,4 @@
       ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
       ! 
-      zcoef = 0.5_wp * r1_rau0
-
-      IF(ln_wd) CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt)
-
       ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
 
@@ -116,10 +116,8 @@
          CALL agrif_ssh( kt )
 # endif
-# if defined key_bdy
-         IF( lk_bdy ) THEN
+         IF( ln_bdy ) THEN
             CALL lbc_lnk( ssha, 'T', 1. )    ! Not sure that's necessary
             CALL bdy_ssh( ssha )             ! Duplicate sea level across open boundaries
          ENDIF
-# endif
       ENDIF
 
@@ -211,11 +209,9 @@
       ENDIF
 
-#if defined key_bdy
-      IF( lk_bdy ) THEN
+      IF( ln_bdy ) THEN
          DO jk = 1, jpkm1
             wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
          END DO
       ENDIF
-#endif
       !
       IF( nn_timing == 1 )  CALL timing_stop('wzv')
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90	(revision 6981)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90	(revision 7646)
@@ -1,3 +1,2 @@
-
 MODULE wet_dry
    !!==============================================================================
@@ -7,6 +6,5 @@
    !! only effects if wetting/drying is on (ln_wd == .true.)
    !!==============================================================================
-   !! History :   
-   !!  NEMO      3.6  ! 2014-09  ((H.Liu)  Original code
+   !! History :  3.6  ! 2014-09  ((H.Liu)  Original code
    !!                 ! will add the runoff and periodic BC case later
    !!----------------------------------------------------------------------
@@ -33,6 +31,7 @@
    !! ---------------------------------------------------------------------
 
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wduflt, wdvflt !: u- and v- filter
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdmask         !: u- and v- limiter 
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   ht_wd          !: wetting and drying t-pt depths
+                                                                     !  (can include negative depths)
 
    LOGICAL,  PUBLIC  ::   ln_wd       !: Wetting/drying activation switch (T:on,F:off)
@@ -77,5 +76,5 @@
          WRITE(numout,*)
          WRITE(numout,*) 'wad_init : Wetting and drying initialization through namelist read'
-         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) '~~~~~~~~'
          WRITE(numout,*) '   Namelist namwad'
          WRITE(numout,*) '      Logical activation                 ln_wd      = ', ln_wd
@@ -84,11 +83,13 @@
          WRITE(numout,*) '      land elevation threshold         rn_wdld      = ', rn_wdld
          WRITE(numout,*) '      Max iteration for W/D limiter    nn_wdit      = ', nn_wdit
-       ENDIF
-
+      ENDIF
+      !
       IF(ln_wd) THEN
-         ALLOCATE( wduflt(jpi,jpj), wdvflt(jpi,jpj), wdmask(jpi,jpj), STAT=ierr )
+         ALLOCATE( wdmask(jpi,jpj), ht_wd(jpi,jpj),  STAT=ierr )
          IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error')
       ENDIF
+      !
    END SUBROUTINE wad_init
+
 
    SUBROUTINE wad_lmt( sshb1, sshemp, z2dt )
@@ -107,5 +108,5 @@
       !
       INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
-      INTEGER  ::   zflag               ! local scalar
+      INTEGER  ::   jflag               ! local scalar
       REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
       REAL(wp) ::   zzflxp, zzflxn      ! local scalars
@@ -116,5 +117,4 @@
       REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv            ! local 2D workspace
       REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
-
       !!----------------------------------------------------------------------
       !
@@ -131,5 +131,5 @@
         !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting'
        
-        zflag  = 0
+        jflag  = 0
         zdepwd = 50._wp   !maximum depth on which that W/D could possibly happen
 
@@ -145,6 +145,6 @@
         ! Horizontal Flux in u and v direction
         DO jk = 1, jpkm1  
-           DO jj = 1, jpjm1
-              DO ji = 1, jpim1
+           DO jj = 1, jpj
+              DO ji = 1, jpi
                  zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
                  zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
@@ -156,9 +156,10 @@
         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
+        wdmask(:,:) = 1
+        DO jj = 2, jpj
+           DO ji = 2, jpi 
+
+             IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE   ! we don't care about land cells
+             IF( ht_wd(ji,jj) > zdepwd )      CYCLE   ! and cells which are unlikely to dry
 
               zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
@@ -167,8 +168,12 @@
                            & 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)
+              zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1
+              IF(zdep2 .le. 0._wp) THEN  !add more safty, but not necessary
+                sshb1(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)
+                IF(zflxu(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = 0._wp
+                IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp
+                IF(zflxv(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = 0._wp
+                IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp 
+                wdmask(ji,jj) = 0._wp
               END IF
            ENDDO
@@ -182,11 +187,11 @@
            zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
            zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
-          
-           DO jj = 2, jpjm1
-              DO ji = 2, jpim1 
+           jflag = 0     ! flag indicating if any further iterations are needed
+          
+           DO jj = 2, jpj
+              DO ji = 2, jpi 
         
-                 wdmask(ji,jj) = 0
-                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 
-                 IF(bathy(ji,jj) > zdepwd) CYCLE
+                 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE 
+                 IF( ht_wd(ji,jj) > zdepwd )      CYCLE
         
                  ztmp = e1e2t(ji,jj)
@@ -198,16 +203,23 @@
           
                  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
+                 zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)
+          
+                 IF( zdep1 > zdep2 ) THEN
+                   wdmask(ji, jj) = 0
                    zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
-                   zcoef = max(zcoef, 0._wp)
+                   !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
+                   ! flag if the limiter has been used but stop flagging if the only
+                   ! changes have zeroed the coefficient since further iterations will
+                   ! not change anything
+                   IF( zcoef > 0._wp ) THEN
+                      jflag = 1 
+                   ELSE
+                      zcoef = 0._wp
+                   ENDIF
                    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
+                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
                  END IF
               END DO ! ji loop
@@ -217,9 +229,8 @@
            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?
+           IF(lk_mpp) CALL mpp_max(jflag)   !max over the global domain
+
+           IF(jflag == 0) EXIT
+          
         END DO  ! jk1 loop
        
@@ -231,6 +242,11 @@
         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!!!'
+      !
+        un_b(:,:) = un_b(:,:) * zwdlmtu(:, :)
+        vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :)
+        CALL lbc_lnk( un_b, 'U', -1. )
+        CALL lbc_lnk( vn_b, 'V', -1. )
+       
+        IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
        
         !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
@@ -240,9 +256,11 @@
         CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
         CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
-      !
-      END IF
-
+        !
+      ENDIF
+      !
       IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
+      !
    END SUBROUTINE wad_lmt
+
 
    SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt )
@@ -260,5 +278,5 @@
       !
       INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
-      INTEGER  ::   zflag         ! local scalar
+      INTEGER  ::   jflag               ! local scalar
       REAL(wp) ::   z2dt
       REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
@@ -269,8 +287,6 @@
       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')
 
@@ -284,5 +300,5 @@
         !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting'
        
-        zflag  = 0
+        jflag  = 0
         zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes
 
@@ -291,6 +307,4 @@
         zflxp(:,:)   = 0._wp
         zflxn(:,:)   = 0._wp
-        !zflxu(:,:)   = 0._wp
-        !zflxv(:,:)   = 0._wp
 
         zwdlmtu(:,:)  = 1._wp
@@ -299,12 +313,9 @@
         ! 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
+        DO jj = 2, jpj
+           DO ji = 2, jpi 
+
+             IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE   ! we don't care about land cells
+             IF( ht_wd(ji,jj) > zdepwd )      CYCLE   ! and cells which are unlikely to dry
 
               zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
@@ -313,8 +324,11 @@
                            & 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)
+              zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
+              IF(zdep2 .le. 0._wp) THEN  !add more safety, but not necessary
+                sshn_e(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)
+                IF(zflxu(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = 0._wp
+                IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp
+                IF(zflxv(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = 0._wp
+                IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp 
               END IF
            ENDDO
@@ -328,11 +342,11 @@
            zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
            zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
-          
-           DO jj = 2, jpjm1
-              DO ji = 2, jpim1 
+           jflag = 0     ! flag indicating if any further iterations are needed
+          
+           DO jj = 2, jpj
+              DO ji = 2, jpi 
         
-                 wdmask(ji,jj) = 0
-                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 
-                 IF(bathy(ji,jj) > zdepwd) CYCLE
+                 IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE 
+                 IF( ht_wd(ji,jj) > zdepwd )      CYCLE
         
                  ztmp = e1e2t(ji,jj)
@@ -344,17 +358,22 @@
           
                  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)
+                 zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - 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)
+                   !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
+                   ! flag if the limiter has been used but stop flagging if the only
+                   ! changes have zeroed the coefficient since further iterations will
+                   ! not change anything
+                   IF( zcoef > 0._wp ) THEN
+                      jflag = 1 
+                   ELSE
+                      zcoef = 0._wp
+                   ENDIF
                    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
+                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
                  END IF
               END DO ! ji loop
@@ -364,9 +383,8 @@
            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?
+           IF(lk_mpp) CALL mpp_max(jflag)   !max over the global domain
+
+           IF(jflag == 0) EXIT
+          
         END DO  ! jk1 loop
        
@@ -377,5 +395,5 @@
         CALL lbc_lnk( zflxv, 'V', -1. )
        
-        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
+        IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
        
         !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
@@ -385,8 +403,10 @@
         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 wet_dry
