Index: branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg	(revision 4487)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg	(revision 4488)
@@ -323,7 +323,10 @@
    !                       !  Type of the operator :
    ln_dynldf_bilap  =  .true.   !  bilaplacian operator
+   ln_dynldf_lap    =  .false.  !  bilaplacian operator
+   !                       !  Direction of action  :
+   ln_dynldf_level  =  .true.   !  iso-level
+   ln_dynldf_hor    =  .false.  !  horizontal (geopotential)            (require "key_ldfslp" in s-coord.)
                            !  Coefficient
    rn_ahm_0_lap     = 60.0      !  horizontal laplacian eddy viscosity   [m2/s]
-   rn_ahmb_0        =  0.0      !  background eddy viscosity for ldf_iso [m2/s]
    rn_ahm_0_blp     = -1.0e+10  !  horizontal bilaplacian eddy viscosity [m4/s]
 /
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 4487)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 4488)
@@ -191,4 +191,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdep3w_n           !: now depth of T-points (sum of e3w) (m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_n, gdepw_n   !: now depth at T-W  points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_b, gdepw_b   !: before depth at T-W  points (m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_n              !: now    vertical scale factors at  t       point  (m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3u_n  , e3v_n     !:            -      -      -    -   u --v   points (m)
@@ -196,4 +197,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3uw_n , e3vw_n    !:            -      -      -    -   uw--vw  points (m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_b              !: before     -      -      -    -   t       points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3w_b              !: before     -      -      -    -   t       points (m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3u_b  , e3v_b     !:   -        -      -      -    -   u --v   points (m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3uw_b , e3vw_b    !:   -        -      -      -    -   uw--vw  points (m)
@@ -350,5 +352,6 @@
          &      e3t_b   (jpi,jpj,jpk) , e3u_b (jpi,jpj,jpk) , e3v_b (jpi,jpj,jpk) ,                           &
          &      e3uw_b  (jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) ,                                                 &
-         &      e3t_a   (jpi,jpj,jpk) , e3u_a (jpi,jpj,jpk) , e3v_a (jpi,jpj,jpk),                            &
+         &      gdept_b (jpi,jpj,jpk) ,gdepw_b(jpi,jpj,jpk) , e3w_b (jpi,jpj,jpk) ,                           &
+         &      e3t_a   (jpi,jpj,jpk) , e3u_a (jpi,jpj,jpk) , e3v_a (jpi,jpj,jpk) ,                           &
          &      ehu_a    (jpi,jpj)    , ehv_a  (jpi,jpj),                                                     &
          &      ehur_a   (jpi,jpj)    , ehvr_a (jpi,jpj),                                                     &
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 4487)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 4488)
@@ -159,4 +159,5 @@
       CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
       CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
+      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W'  )
       CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
       CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
@@ -166,8 +167,12 @@
       fsdepw_n(:,:,1) = 0.0_wp
       fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
+      fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1)
+      fsdepw_b(:,:,1) = 0.0_wp
       DO jk = 2, jpk
          fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
          fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
          fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
+         fsdept_b(:,:,jk) = fsdept_b(:,:,jk-1) + fse3w_b(:,:,jk)
+         fsdepw_b(:,:,jk) = fsdepw_b(:,:,jk-1) + fse3t_b(:,:,jk-1)
       END DO
       ! Reference water column height at t-, u- and v- point
@@ -600,4 +605,7 @@
          tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
       ENDIF
+      fsdept_b(:,:,:) = fsdept_n(:,:,:)
+      fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
+
       fse3t_n(:,:,:) = fse3t_a(:,:,:)
       fse3u_n(:,:,:) = fse3u_a(:,:,:)
@@ -616,4 +624,5 @@
       CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
       CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
+      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W'  )
       CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
       CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 4487)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 4488)
@@ -28,4 +28,8 @@
 #   define  fse3uw_n(i,j,k)  e3uw_n(i,j,k)
 #   define  fse3vw_n(i,j,k)  e3vw_n(i,j,k)
+
+#   define  fsdept_b(i,j,k)  gdept_b(i,j,k)
+#   define  fsdepw_b(i,j,k)  gdepw_b(i,j,k)
+#   define  fse3w_b(i,j,k)   e3w_b(i,j,k)
 
 #   define  fse3t_a(i,j,k)   e3t_a(i,j,k)
@@ -78,4 +82,8 @@
 #   define  fse3vw_n(i,j,k)  e3vw_0(i,j,k)
 
+#   define  fsdept_b(i,j,k)  gdept_0(i,j,k)
+#   define  fsdepw_b(i,j,k)  gdepw_0(i,j,k)
+#   define  fse3w_b(i,j,k)   e3w_0(i,j,k)
+
 #   define  fse3t_a(i,j,k)   e3t_0(i,j,k)
 #   define  fse3u_a(i,j,k)   e3u_0(i,j,k)
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 4487)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 4488)
@@ -20,4 +20,5 @@
    USE dynldf_iso     ! lateral mixing            (dyn_ldf_iso    routine)
    USE dynldf_lap     ! lateral mixing            (dyn_ldf_lap    routine)
+   USE ldftra_oce, ONLY: ln_traldf_hor     ! ocean tracers lateral physics
    USE trdmod         ! ocean dynamics and tracer trends
    USE trdmod_oce     ! ocean variables trends
@@ -152,4 +153,6 @@
       IF( ioptio >  1 ) CALL ctl_stop( '          use only ONE direction (level/hor/iso)' )
 
+      IF( ln_dynldf_iso .AND. ln_traldf_hor ) CALL ctl_stop( 'Not sensible to use geopotential diffusion for tracers with isoneutral diffusion for dynamics' )
+
       !                                   ! Set nldf, the type of lateral diffusion, from ln_dynldf_... logicals
       ierr = 0
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90	(revision 4487)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90	(revision 4488)
@@ -19,4 +19,5 @@
    USE dom_oce         ! ocean space and time domain
    USE ldfdyn_oce      ! ocean dynamics lateral physics
+   USE ldftra_oce, ONLY: ln_traldf_iso
    USE zdf_oce         ! ocean vertical physics
    USE trdmod          ! ocean dynamics trends 
@@ -29,5 +30,4 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
-
 
    IMPLICIT NONE
@@ -104,5 +104,30 @@
          IF( dyn_ldf_bilapg_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_bilapg: failed to allocate arrays')
       ENDIF
-      !
+
+      ! s-coordinate: Iso-level diffusion on tracer, but geopotential level diffusion on momentum
+      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
+         !
+         DO jk = 1, jpk         ! set the slopes of iso-level
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk)
+                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)
+                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
+                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
+               END DO
+            END DO
+         END DO
+         ! Lateral boundary conditions on the slopes
+         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
+         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
+ 
+!!bug
+         IF( kt == nit000 ) then
+            IF(lwp) WRITE(numout,*) ' max slop: u', SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)),  &
+               &                             ' wi', sqrt(MAXVAL(wslpi))     , ' wj', sqrt(MAXVAL(wslpj))
+         endif
+!!end
+      ENDIF
+
       zwk1(:,:,:) = 0.e0   ;   zwk3(:,:,:) = 0.e0
       zwk2(:,:,:) = 0.e0   ;   zwk4(:,:,:) = 0.e0
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90	(revision 4487)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90	(revision 4488)
@@ -131,14 +131,14 @@
       ENDIF
 
-      ! s-coordinate: Iso-level diffusion on momentum but not on tracer
+      ! s-coordinate: Iso-level diffusion on tracer, but geopotential level diffusion on momentum
       IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
          !
          DO jk = 1, jpk         ! set the slopes of iso-level
             DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)
-                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)
-                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
-                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
+               DO ji = 2, jpim1
+                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk)
+                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)
+                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
+                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
                END DO
             END DO
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 4487)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 4488)
@@ -117,249 +117,291 @@
       CALL wrk_alloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv )
 
-      zeps   =  1.e-20_wp        !==   Local constant initialization   ==!
-      z1_16  =  1.0_wp / 16._wp
-      zm1_g  = -1.0_wp / grav
-      zm1_2g = -0.5_wp / grav
-      !
-      zww(:,:,:) = 0._wp
-      zwz(:,:,:) = 0._wp
-      !
-      DO jk = 1, jpk             !==   i- & j-gradient of density   ==!
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj  ,jk) - prd(ji,jj,jk) )
-               zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji  ,jj+1,jk) - prd(ji,jj,jk) )
-            END DO
-         END DO
-      END DO
-      IF( ln_zps ) THEN                           ! partial steps correction at the bottom ocean level
+      IF ( ln_traldf_iso .OR. ln_dynldf_iso ) THEN 
+     
+         zeps   =  1.e-20_wp        !==   Local constant initialization   ==!
+         z1_16  =  1.0_wp / 16._wp
+         zm1_g  = -1.0_wp / grav
+         zm1_2g = -0.5_wp / grav
+         !
+         zww(:,:,:) = 0._wp
+         zwz(:,:,:) = 0._wp
+         !
+         DO jk = 1, jpk             !==   i- & j-gradient of density   ==!
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj  ,jk) - prd(ji,jj,jk) )
+                  zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji  ,jj+1,jk) - prd(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+         IF( ln_zps ) THEN                           ! partial steps correction at the bottom ocean level
 # if defined key_vectopt_loop
-         DO jj = 1, 1
-            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+            DO jj = 1, 1
+               DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
 # else
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
 # endif
-               zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj)
-               zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj)
-            END DO
-         END DO
+                  zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj)
+                  zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj)
+               END DO
+            END DO
+         ENDIF
+         !
+         zdzr(:,:,1) = 0._wp        !==   Local vertical density gradient at T-point   == !   (evaluated from N^2)
+         DO jk = 2, jpkm1
+            !                                ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
+            !                                !   trick: tmask(ik  )  = 0   =>   all pn2   = 0   =>   zdzr = 0
+            !                                !    else  tmask(ik+1)  = 0   =>   pn2(ik+1) = 0   =>   zdzr divides by 1
+            !                                !          umask(ik+1) /= 0   =>   all pn2  /= 0   =>   zdzr divides by 2
+            !                                ! NB: 1/(tmask+1) = (1-.5*tmask)  substitute a / by a *  ==> faster
+            zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp )              &
+               &                 * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) )
+         END DO
+         !
+         !                          !==   Slopes just below the mixed layer   ==!
+         CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr )        ! output: uslpml, vslpml, wslpiml, wslpjml
+
+
+         ! I.  slopes at u and v point      | uslp = d/di( prd ) / d/dz( prd )
+         ! ===========================      | vslp = d/dj( prd ) / d/dz( prd )
+         !
+         DO jk = 2, jpkm1                            !* Slopes at u and v points
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  !                                      ! horizontal and vertical density gradient at u- and v-points
+                  zau = zgru(ji,jj,jk) / e1u(ji,jj)
+                  zav = zgrv(ji,jj,jk) / e2v(ji,jj)
+                  zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj  ,jk) )
+                  zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji  ,jj+1,jk) )
+                  !                                      ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
+                  !                                      ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
+                  zbu = MIN(  zbu, -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau )  )
+                  zbv = MIN(  zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav )  )
+                  !                                      ! uslp and vslp output in zwz and zww, resp.
+                  zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
+                  zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )
+                  zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps )                                              &
+                     &                   + zfi  * uslpml(ji,jj)                                                     &
+                     &                          * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) )   &
+                     &                          / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk)
+                  zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps )                                              &
+                     &                   + zfj  * vslpml(ji,jj)                                                     &
+                     &                          * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) )   &
+                     &                          / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk)
+!!gm  modif to suppress omlmask.... (as in Griffies case)
+!                  !                                         ! jk must be >= ML level for zf=1. otherwise  zf=0.
+!                  zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp )
+!                  zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp )
+!                  zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )
+!                  zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) )
+!                  zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk)
+!                  zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk)
+!!gm end modif
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( zwz, 'U', -1. )   ;   CALL lbc_lnk( zww, 'V', -1. )      ! lateral boundary conditions
+         !
+         !                                            !* horizontal Shapiro filter
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only
+               DO ji = 2, jpim1
+                  uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
+                     &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
+                     &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
+                     &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
+                     &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
+                  vslp(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
+                     &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
+                     &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
+                     &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
+                     &                       + 4.*  zww(ji,jj    ,jk)                       )
+               END DO
+            END DO
+            DO jj = 3, jpj-2                               ! other rows
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
+                     &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
+                     &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
+                     &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
+                     &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
+                  vslp(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
+                     &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
+                     &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)         &
+                     &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
+                     &                       + 4.*  zww(ji,jj    ,jk)                       )
+               END DO
+            END DO
+            !                                        !* decrease along coastal boundaries
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk  ) ) * 0.5_wp   &
+                     &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5_wp
+                  vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk  ) ) * 0.5_wp   &
+                     &                            * ( vmask(ji  ,jj,jk) + vmask(ji  ,jj,jk+1) ) * 0.5_wp
+               END DO
+            END DO
+         END DO
+
+
+         ! II.  slopes at w point           | wslpi = mij( d/di( prd ) / d/dz( prd )
+         ! ===========================      | wslpj = mij( d/dj( prd ) / d/dz( prd )
+         !
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  !                                  !* Local vertical density gradient evaluated from N^2
+                  zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. )
+                  !                                  !* Slopes at w point
+                  !                                        ! i- & j-gradient of density at w-points
+                  zci = MAX(  umask(ji-1,jj,jk  ) + umask(ji,jj,jk  )           &
+                     &      + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps  ) * e1t(ji,jj)
+                  zcj = MAX(  vmask(ji,jj-1,jk  ) + vmask(ji,jj,jk-1)           &
+                     &      + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk  ) , zeps  ) *  e2t(ji,jj)
+                  zai =    (  zgru (ji-1,jj,jk  ) + zgru (ji,jj,jk-1)           &
+                     &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci * tmask (ji,jj,jk)
+                  zaj =    (  zgrv (ji,jj-1,jk  ) + zgrv (ji,jj,jk-1)           &
+                     &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj * tmask (ji,jj,jk)
+                  !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
+                  !                                        ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
+                  zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zai )  )
+                  zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj )  )
+                  !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.)
+                  zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0
+                  zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp )
+                  zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk  ) * tmask(ji,jj,jk)
+                  zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk  ) * tmask(ji,jj,jk)
+
+!!gm  modif to suppress omlmask....  (as in Griffies operator)
+!                  !                                         ! jk must be >= ML level for zfk=1. otherwise  zfk=0.
+!                  zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp )
+!                  zck = fsdepw(ji,jj,jk)    / MAX( hmlp(ji,jj), 10. )
+!                  zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk)
+!                  zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk)
+!!gm end modif
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( zwz, 'T', -1. )   ;    CALL lbc_lnk( zww, 'T', -1. )      ! lateral boundary conditions
+         !
+         !                                           !* horizontal Shapiro filter
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only
+               DO ji = 2, jpim1
+                  zcofw = tmask(ji,jj,jk) * z1_16
+                  wslpi(ji,jj,jk) = (          zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
+                       &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
+                       &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
+                       &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
+                       &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * zcofw
+
+                  wslpj(ji,jj,jk) = (          zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
+                       &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
+                       &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
+                       &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
+                       &                + 4.*  zww(ji  ,jj  ,jk)                         ) * zcofw
+               END DO
+            END DO
+            DO jj = 3, jpj-2                               ! other rows
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zcofw = tmask(ji,jj,jk) * z1_16
+                  wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
+                       &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
+                       &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
+                       &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
+                       &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * zcofw
+
+                  wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
+                       &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
+                       &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
+                       &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
+                       &                + 4.*  zww(ji  ,jj  ,jk)                         ) * zcofw
+               END DO
+            END DO
+            !                                        !* decrease along coastal boundaries
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zck =   ( umask(ji,jj,jk) + umask(ji-1,jj,jk) )   &
+                     &  * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25
+                  wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck
+                  wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck
+               END DO
+            END DO
+         END DO
+
+         ! III.  Specific grid points
+         ! ===========================
+         !
+         IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN     !  ORCA_R4 configuration: horizontal diffusion in specific area
+            !                                                    ! Gibraltar Strait
+            ij0 =  50   ;   ij1 =  53
+            ii0 =  69   ;   ii1 =  71   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
+            ij0 =  51   ;   ij1 =  53
+            ii0 =  68   ;   ii1 =  71   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
+            ii0 =  69   ;   ii1 =  71   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
+            ii0 =  69   ;   ii1 =  71   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
+            !
+            !                                                    ! Mediterrannean Sea
+            ij0 =  49   ;   ij1 =  56
+            ii0 =  71   ;   ii1 =  90   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
+            ij0 =  50   ;   ij1 =  56
+            ii0 =  70   ;   ii1 =  90   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
+            ii0 =  71   ;   ii1 =  90   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
+            ii0 =  71   ;   ii1 =  90   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
+         ENDIF
+
+
+         ! IV. Lateral boundary conditions
+         ! ===============================
+         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
+         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
+
+
+         IF(ln_ctl) THEN
+            CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp  - u : ', tab3d_2=vslp,  clinfo2=' v : ', kdim=jpk)
+            CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
+         ENDIF
+         !
+
+      ELSEIF ( lk_vvl ) THEN 
+ 
+         IF(lwp) THEN 
+            WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 
+         ENDIF 
+
+         ! geopotential diffusion in s-coordinates on tracers and/or momentum 
+         ! The slopes of s-surfaces are computed at each time step due to vvl 
+         ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 
+
+         ! set the slope of diffusion to the slope of s-surfaces 
+         !      ( c a u t i o n : minus sign as fsdep has positive value ) 
+         DO jk = 1, jpk 
+            DO jj = 2, jpjm1 
+               DO ji = fs_2, fs_jpim1   ! vector opt. 
+                  uslp(ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) 
+                  vslp(ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 
+                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 
+                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 
+               END DO 
+            END DO 
+         END DO 
+
+         ! Lateral boundary conditions on the slopes 
+         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. ) 
+         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. ) 
+  
+         if( kt == nit000 ) then 
+            IF(lwp) WRITE(numout,*) ' max slop: u',SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)),  & 
+               &                             ' wi', sqrt(MAXVAL(wslpi)), ' wj', sqrt(MAXVAL(wslpj)) 
+         endif 
+  
+         IF(ln_ctl) THEN 
+            CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp  - u : ', tab3d_2=vslp,  clinfo2=' v : ', kdim=jpk) 
+            CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) 
+         ENDIF 
+
       ENDIF
-      !
-      zdzr(:,:,1) = 0._wp        !==   Local vertical density gradient at T-point   == !   (evaluated from N^2)
-      DO jk = 2, jpkm1
-         !                                ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
-         !                                !   trick: tmask(ik  )  = 0   =>   all pn2   = 0   =>   zdzr = 0
-         !                                !    else  tmask(ik+1)  = 0   =>   pn2(ik+1) = 0   =>   zdzr divides by 1
-         !                                !          umask(ik+1) /= 0   =>   all pn2  /= 0   =>   zdzr divides by 2
-         !                                ! NB: 1/(tmask+1) = (1-.5*tmask)  substitute a / by a *  ==> faster
-         zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp )              &
-            &                 * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) )
-      END DO
-      !
-      !                          !==   Slopes just below the mixed layer   ==!
-      CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr )        ! output: uslpml, vslpml, wslpiml, wslpjml
-
-
-      ! I.  slopes at u and v point      | uslp = d/di( prd ) / d/dz( prd )
-      ! ===========================      | vslp = d/dj( prd ) / d/dz( prd )
-      !
-      DO jk = 2, jpkm1                            !* Slopes at u and v points
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               !                                      ! horizontal and vertical density gradient at u- and v-points
-               zau = zgru(ji,jj,jk) / e1u(ji,jj)
-               zav = zgrv(ji,jj,jk) / e2v(ji,jj)
-               zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj  ,jk) )
-               zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji  ,jj+1,jk) )
-               !                                      ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
-               !                                      ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
-               zbu = MIN(  zbu, -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau )  )
-               zbv = MIN(  zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav )  )
-               !                                      ! uslp and vslp output in zwz and zww, resp.
-               zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
-               zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )
-               zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps )                                              &
-                  &                   + zfi  * uslpml(ji,jj)                                                     &
-                  &                          * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) )   &
-                  &                          / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk)
-               zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps )                                              &
-                  &                   + zfj  * vslpml(ji,jj)                                                     &
-                  &                          * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) )   &
-                  &                          / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk)
-!!gm  modif to suppress omlmask.... (as in Griffies case)
-!               !                                         ! jk must be >= ML level for zf=1. otherwise  zf=0.
-!               zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp )
-!               zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp )
-!               zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )
-!               zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) )
-!               zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk)
-!               zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk)
-!!gm end modif
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( zwz, 'U', -1. )   ;   CALL lbc_lnk( zww, 'V', -1. )      ! lateral boundary conditions
-      !
-      !                                            !* horizontal Shapiro filter
-      DO jk = 2, jpkm1
-         DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only
-            DO ji = 2, jpim1
-               uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
-                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
-                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
-                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
-                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
-               vslp(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
-                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
-                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
-                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
-                  &                       + 4.*  zww(ji,jj    ,jk)                       )
-            END DO
-         END DO
-         DO jj = 3, jpj-2                               ! other rows
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
-                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
-                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
-                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
-                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
-               vslp(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
-                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
-                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
-                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
-                  &                       + 4.*  zww(ji,jj    ,jk)                       )
-            END DO
-         END DO
-         !                                        !* decrease along coastal boundaries
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk  ) ) * 0.5_wp   &
-                  &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5_wp
-               vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk  ) ) * 0.5_wp   &
-                  &                            * ( vmask(ji  ,jj,jk) + vmask(ji  ,jj,jk+1) ) * 0.5_wp
-            END DO
-         END DO
-      END DO
-
-
-      ! II.  slopes at w point           | wslpi = mij( d/di( prd ) / d/dz( prd )
-      ! ===========================      | wslpj = mij( d/dj( prd ) / d/dz( prd )
-      !
-      DO jk = 2, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               !                                  !* Local vertical density gradient evaluated from N^2
-               zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. )
-               !                                  !* Slopes at w point
-               !                                        ! i- & j-gradient of density at w-points
-               zci = MAX(  umask(ji-1,jj,jk  ) + umask(ji,jj,jk  )           &
-                  &      + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps  ) * e1t(ji,jj)
-               zcj = MAX(  vmask(ji,jj-1,jk  ) + vmask(ji,jj,jk-1)           &
-                  &      + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk  ) , zeps  ) *  e2t(ji,jj)
-               zai =    (  zgru (ji-1,jj,jk  ) + zgru (ji,jj,jk-1)           &
-                  &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci * tmask (ji,jj,jk)
-               zaj =    (  zgrv (ji,jj-1,jk  ) + zgrv (ji,jj,jk-1)           &
-                  &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj * tmask (ji,jj,jk)
-               !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
-               !                                        ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
-               zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zai )  )
-               zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj )  )
-               !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.)
-               zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0
-               zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp )
-               zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk  ) * tmask(ji,jj,jk)
-               zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk  ) * tmask(ji,jj,jk)
-
-!!gm  modif to suppress omlmask....  (as in Griffies operator)
-!               !                                         ! jk must be >= ML level for zfk=1. otherwise  zfk=0.
-!               zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp )
-!               zck = fsdepw(ji,jj,jk)    / MAX( hmlp(ji,jj), 10. )
-!               zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk)
-!               zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk)
-!!gm end modif
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( zwz, 'T', -1. )   ;    CALL lbc_lnk( zww, 'T', -1. )      ! lateral boundary conditions
-      !
-      !                                           !* horizontal Shapiro filter
-      DO jk = 2, jpkm1
-         DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only
-            DO ji = 2, jpim1
-               zcofw = tmask(ji,jj,jk) * z1_16
-               wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
-                    &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
-                    &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
-                    &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
-                    &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * zcofw
-
-               wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
-                    &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
-                    &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
-                    &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
-                    &                + 4.*  zww(ji  ,jj  ,jk)                         ) * zcofw
-            END DO
-         END DO
-         DO jj = 3, jpj-2                               ! other rows
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zcofw = tmask(ji,jj,jk) * z1_16
-               wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
-                    &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
-                    &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
-                    &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
-                    &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * zcofw
-
-               wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
-                    &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
-                    &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
-                    &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
-                    &                + 4.*  zww(ji  ,jj  ,jk)                         ) * zcofw
-            END DO
-         END DO
-         !                                        !* decrease along coastal boundaries
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zck =   ( umask(ji,jj,jk) + umask(ji-1,jj,jk) )   &
-                  &  * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25
-               wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck
-               wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck
-            END DO
-         END DO
-      END DO
-
-      ! III.  Specific grid points
-      ! ===========================
-      !
-      IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN     !  ORCA_R4 configuration: horizontal diffusion in specific area
-         !                                                    ! Gibraltar Strait
-         ij0 =  50   ;   ij1 =  53
-         ii0 =  69   ;   ii1 =  71   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
-         ij0 =  51   ;   ij1 =  53
-         ii0 =  68   ;   ii1 =  71   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
-         ii0 =  69   ;   ii1 =  71   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
-         ii0 =  69   ;   ii1 =  71   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
-         !
-         !                                                    ! Mediterrannean Sea
-         ij0 =  49   ;   ij1 =  56
-         ii0 =  71   ;   ii1 =  90   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
-         ij0 =  50   ;   ij1 =  56
-         ii0 =  70   ;   ii1 =  90   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
-         ii0 =  71   ;   ii1 =  90   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
-         ii0 =  71   ;   ii1 =  90   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
-      ENDIF
-
-
-      ! IV. Lateral boundary conditions
-      ! ===============================
-      CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
-      CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
-
-
-      IF(ln_ctl) THEN
-         CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp  - u : ', tab3d_2=vslp,  clinfo2=' v : ', kdim=jpk)
-         CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
-      ENDIF
-      !
+      
       CALL wrk_dealloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv )
       !
@@ -783,6 +825,5 @@
          wslpj(:,:,:) = 0._wp   ;   wslpjml(:,:) = 0._wp
 
-         !!gm I no longer understand this.....
-         IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN
+         IF( ln_traldf_hor .OR. ln_dynldf_hor ) THEN
             IF(lwp)   WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
 
@@ -796,8 +837,8 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)
-                     vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)
-                     wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
-                     wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
+                     uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk)
+                     vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)
+                     wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
+                     wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
                   END DO
                END DO
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 4487)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 4488)
@@ -204,4 +204,5 @@
       ENDIF
 
+      IF( nldf == 3 )   CALL ctl_warn( 'geopotential bilaplacian tracer diffusion in s-coords not thoroughly tested' )
       IF( ierr == 1 )   CALL ctl_stop( ' iso-level in z-coordinate - partial step, not allowed' )
       IF( ierr == 2 )   CALL ctl_stop( ' isoneutral bilaplacian operator does not exist' )
