Index: branches/UKMO/r5518_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- branches/UKMO/r5518_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 7151)
+++ branches/UKMO/r5518_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 7152)
@@ -11,11 +11,10 @@
    !!            3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)  add Griffies operator
    !!             -   ! 2010-11  (F. Dupond, G. Madec)  bug correction in slopes just below the ML
+   !!            3.7  ! 2013-12  (F. Lemarie, G. Madec)  add limiter on triad slopes
    !!----------------------------------------------------------------------
-#if   defined key_ldfslp   ||   defined key_esopa
+
    !!----------------------------------------------------------------------
-   !!   'key_ldfslp'                      Rotation of lateral mixing tensor
-   !!----------------------------------------------------------------------
-   !!   ldf_slp_grif  : calculates the triads of isoneutral slopes (Griffies operator)
    !!   ldf_slp       : calculates the slopes of neutral surface   (Madec operator)
+   !!   ldf_slp_triad : calculates the triads of isoneutral slopes (Griffies operator)
    !!   ldf_slp_mxl   : calculates the slopes at the base of the mixed layer (Madec operator)
    !!   ldf_slp_init  : initialization of the slopes computation
@@ -23,6 +22,5 @@
    USE oce            ! ocean dynamics and tracers
    USE dom_oce        ! ocean space and time domain
-   USE ldftra_oce     ! lateral diffusion: traceur
-   USE ldfdyn_oce     ! lateral diffusion: dynamics
+   USE ldfdyn         ! lateral diffusion: eddy viscosity coef.
    USE phycst         ! physical constants
    USE zdfmxl         ! mixed layer depth
@@ -30,29 +28,42 @@
    !
    USE in_out_manager ! I/O manager
+   USE prtctl         ! Print control
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE prtctl         ! Print control
+   USE lib_mpp        ! distribued memory computing library
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
    USE wrk_nemo       ! work arrays
    USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   ldf_slp        ! routine called by step.F90
-   PUBLIC   ldf_slp_grif   ! routine called by step.F90
-   PUBLIC   ldf_slp_init   ! routine called by opa.F90
-
-   LOGICAL , PUBLIC, PARAMETER ::   lk_ldfslp = .TRUE.     !: slopes flag
-   !                                                                             !! Madec operator
-   !  Arrays allocated in ldf_slp_init() routine once we know whether we're using the Griffies or Madec operator
+   PUBLIC   ldf_slp         ! routine called by step.F90
+   PUBLIC   ldf_slp_triad   ! routine called by step.F90
+   PUBLIC   ldf_slp_init    ! routine called by nemogcm.F90
+
+   LOGICAL , PUBLIC ::   l_ldfslp = .FALSE.     !: slopes flag
+
+   LOGICAL , PUBLIC ::   ln_traldf_iso   = .TRUE.       !: iso-neutral direction
+   LOGICAL , PUBLIC ::   ln_traldf_triad = .FALSE.      !: griffies triad scheme
+
+   LOGICAL , PUBLIC ::   ln_triad_iso    = .FALSE.      !: pure horizontal mixing in ML
+   LOGICAL , PUBLIC ::   ln_botmix_triad = .FALSE.      !: mixing on bottom
+   REAL(wp), PUBLIC ::   rn_sw_triad     = 1._wp        !: =1 switching triads ; =0 all four triads used 
+   REAL(wp), PUBLIC ::   rn_slpmax       = 0.01_wp      !: slope limit
+
+   LOGICAL , PUBLIC ::   l_grad_zps = .FALSE.           !: special treatment for Horz Tgradients w partial steps (triad operator)
+   
+   !                                                     !! Classic operator (Madec)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   uslp, wslpi          !: i_slope at U- and W-points
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   vslp, wslpj          !: j-slope at V- and W-points
-   !                                                                !! Griffies operator
+   !                                                     !! triad operator (Griffies)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   wslp2                !: wslp**2 from Griffies quarter cells
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) ::   triadi_g, triadj_g   !: skew flux  slopes relative to geopotentials
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) ::   triadi  , triadj     !: isoneutral slopes relative to model-coordinate
-
-   !                                                              !! Madec operator
-   !  Arrays allocated in ldf_slp_init() routine once we know whether we're using the Griffies or Madec operator
+   !                                                     !! both operators
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   ah_wslp2             !: ah * slope^2 at w-point
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   akz                  !: stabilizing vertical diffusivity
+   
+   !                                                     !! Madec operator
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   omlmask           ! mask of the surface mixed layer at T-pt
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   uslpml, wslpiml   ! i_slope at U- and W-points just below the mixed layer
@@ -63,9 +74,7 @@
    !! * Substitutions
 #  include "domzgr_substitute.h90"
-#  include "ldftra_substitute.h90"
-#  include "ldfeiv_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2014)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -105,13 +114,11 @@
       INTEGER  ::   ii0, ii1, iku   ! temporary integer
       INTEGER  ::   ij0, ij1, ikv   ! temporary integer
-      REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16, zcofw ! local scalars
+      REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16, zcofw, z1_slpmax ! local scalars
       REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      -
       REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      -
       REAL(wp) ::   zck, zfk,      zbw             !   -      -
-      REAL(wp) ::   zdepv, zdepu         !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwz, zww
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdzr
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zgru, zgrv
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zhmlpu, zhmlpv
       !!----------------------------------------------------------------------
       !
@@ -119,349 +126,225 @@
       !
       CALL wrk_alloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv )
-      CALL wrk_alloc( jpi,jpj, zhmlpu, zhmlpv )
-
-      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
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj)
-                  zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj)
-               END DO
-            END DO
-         ENDIF
-         IF( ln_zps .AND. ln_isfcav ) THEN           ! partial steps correction at the bottom ocean level
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-               IF ( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj) 
-               IF ( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj)
-               END DO
-            END DO
-         ENDIF
-         !
-         !==   Local vertical density gradient at T-point   == !   (evaluated from N^2)
-         ! interior value
-         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
-         ! surface initialisation 
-         zdzr(:,:,1) = 0._wp 
-         IF ( ln_isfcav ) THEN
-            ! if isf need to overwrite the interior value at at the first ocean point
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  zdzr(ji,jj,1:mikt(ji,jj)) = 0._wp
-               END DO
-            END DO
-         END IF
-         !
-         !                          !==   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 )
-         !
-         IF ( ln_isfcav ) THEN
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  IF (miku(ji,jj) .GT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp)
-                  IF (miku(ji,jj) .LT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji+1,jj  ),                   5._wp)
-                  IF (miku(ji,jj) .EQ. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji+1,jj  ), 5._wp)
-                  IF (mikv(ji,jj) .GT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp)
-                  IF (mikv(ji,jj) .LT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj+1),                   5._wp)
-                  IF (mikv(ji,jj) .EQ. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji  ,jj+1), 5._wp)
-               ENDDO
-            ENDDO
-         ELSE
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zhmlpu(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp)
-                  zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp)
-               ENDDO
-            ENDDO
-         END IF
-         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) )
-                  ! thickness of water column between surface and level k at u/v point
-                  zdepu = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji+1,jj  ,jk) )                              &
-                                   - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj  ) ) - fse3u(ji,jj,miku(ji,jj)) )
-                  zdepv = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji  ,jj+1,jk) )                              &
-                                   - 2 * MAX( risfdep(ji,jj), risfdep(ji  ,jj+1) ) - fse3v(ji,jj,mikv(ji,jj)) )
-                  !
-                  zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps )                                          &
-                     &                 + zfi  * uslpml(ji,jj) * zdepu / zhmlpu(ji,jj)
-                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * wumask(ji,jj,jk)
-                  zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps )                                          &
-                     &                 + zfj  * vslpml(ji,jj) * zdepv / zhmlpv(ji,jj) 
-                  zww(ji,jj,jk) = zww(ji,jj,jk) * wvmask(ji,jj,jk)
-                  
-                 
+
+      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
+      z1_slpmax = 1._wp / rn_slpmax
+      !
+      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
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               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) * r1_e1u(ji,jj)
+               zav = zgrv(ji,jj,jk) * r1_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, - z1_slpmax * ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau )  )
+               zbv = MIN(  zbv, - z1_slpmax * 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)
+!               !                                         ! 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   &
-                     &                            *   umask(ji,jj,jk-1)
-                  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   &
-                     &                            *   vmask(ji,jj,jk-1)
-               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. ) * wmask(ji,jj,jk)
-                  !                                  !* 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
-                  zaj =    (  zgrv (ji,jj-1,jk  ) + zgrv (ji,jj,jk-1)           &
-                     &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj
-                  !                                        ! 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) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj), 10._wp )
-                  zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) &
-                     &            + zck * wslpiml(ji,jj) * zfk  ) * wmask(ji,jj,jk)
-                  zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) &
-                     &            + zck * wslpjml(ji,jj) * zfk  ) * wmask(ji,jj,jk)
+            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)
+!               !                                         ! 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 * wmask(ji,jj,jk)
-                  wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck * wmask(ji,jj,jk)
-               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 jj = 2, jpjm1 
-            DO ji = fs_2, fs_jpim1   ! vector opt. 
-               uslp(ji,jj,1) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,1) - fsdept_b(ji ,jj ,1) )  * umask(ji,jj,1) 
-               vslp(ji,jj,1) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,1) - fsdept_b(ji ,jj ,1) )  * vmask(ji,jj,1) 
-               wslpi(ji,jj,1) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,1) - fsdepw_b(ji-1,jj,1) ) * tmask(ji,jj,1) * 0.5 
-               wslpj(ji,jj,1) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,1) - fsdepw_b(ji,jj-1,1) ) * tmask(ji,jj,1) * 0.5 
-            END DO 
-         END DO 
-
-         DO jk = 2, 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) ) &
-                    &                              * wmask(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) ) &
-                    &                              * wmask(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 
-
+            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 = wmask(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 = wmask(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 in vicinity of topography
+         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
+
+      ! 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 )
-      CALL wrk_dealloc( jpi,jpj,     zhmlpu, zhmlpv)
       !
       IF( nn_timing == 1 )  CALL timing_stop('ldf_slp')
@@ -470,10 +353,10 @@
 
 
-   SUBROUTINE ldf_slp_grif ( kt )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE ldf_slp_grif  ***
+   SUBROUTINE ldf_slp_triad ( kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE ldf_slp_triad  ***
       !!
       !! ** Purpose :   Compute the squared slopes of neutral surfaces (slope
-      !!      of iso-pycnal surfaces referenced locally) (ln_traldf_grif=T)
+      !!      of iso-pycnal surfaces referenced locally) (ln_traldf_triad=T)
       !!      at W-points using the Griffies quarter-cells.
       !!
@@ -490,16 +373,20 @@
       REAL(wp) ::   zfacti, zfactj              ! local scalars
       REAL(wp) ::   znot_thru_surface           ! local scalars
-      REAL(wp) ::   zdit, zdis, zdjt, zdjs, zdkt, zdks, zbu, zbv, zbti, zbtj
+      REAL(wp) ::   zdit, zdis, zdkt, zbu, zbti, zisw
+      REAL(wp) ::   zdjt, zdjs, zdks, zbv, zbtj, zjsw
       REAL(wp) ::   zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_g_raw, zti_g_lim
       REAL(wp) ::   zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_g_raw, ztj_g_lim
       REAL(wp) ::   zdzrho_raw
+      REAL(wp) ::   zbeta0, ze3_e1, ze3_e2
       REAL(wp), POINTER, DIMENSION(:,:)     ::   z1_mlbw
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalbet
       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zdxrho , zdyrho, zdzrho     ! Horizontal and vertical density gradients
       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zti_mlb, ztj_mlb            ! for Griffies operator only
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('ldf_slp_grif')
+      IF( nn_timing == 1 )  CALL timing_start('ldf_slp_triad')
       !
       CALL wrk_alloc( jpi,jpj, z1_mlbw )
+      CALL wrk_alloc( jpi,jpj,jpk, zalbet )
       CALL wrk_alloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho,              klstart = 0  )
       CALL wrk_alloc( jpi,jpj,  2,2, zti_mlb, ztj_mlb,        kkstart = 0, klstart = 0  )
@@ -519,8 +406,8 @@
                   zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! j-gradient of T & S at v-point
                   zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
-                  zdxrho_raw = ( - rab_b(ji+ip,jj   ,jk,jp_tem) * zdit + rab_b(ji+ip,jj   ,jk,jp_sal) * zdis ) / e1u(ji,jj)
-                  zdyrho_raw = ( - rab_b(ji   ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji   ,jj+jp,jk,jp_sal) * zdjs ) / e2v(ji,jj)
-                  zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign
-                  zdyrho(ji   ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw )
+                  zdxrho_raw = ( - rab_b(ji+ip,jj   ,jk,jp_tem) * zdit + rab_b(ji+ip,jj   ,jk,jp_sal) * zdis ) * r1_e1u(ji,jj)
+                  zdyrho_raw = ( - rab_b(ji   ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji   ,jj+jp,jk,jp_sal) * zdjs ) * r1_e2v(ji,jj)
+                  zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN(  MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw  )   ! keep the sign
+                  zdyrho(ji   ,jj+jp,jk,1-jp) = SIGN(  MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw  )
                END DO
             END DO
@@ -533,6 +420,6 @@
                   zdit = gtsu(ji,jj,jp_tem)   ;   zdjt = gtsv(ji,jj,jp_tem)      ! i- & j-gradient of Temperature
                   zdis = gtsu(ji,jj,jp_sal)   ;   zdjs = gtsv(ji,jj,jp_sal)      ! i- & j-gradient of Salinity
-                  zdxrho_raw = ( - rab_b(ji+ip,jj   ,iku,jp_tem) * zdit + rab_b(ji+ip,jj   ,iku,jp_sal) * zdis ) / e1u(ji,jj)
-                  zdyrho_raw = ( - rab_b(ji   ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji   ,jj+jp,ikv,jp_sal) * zdjs ) / e2v(ji,jj)
+                  zdxrho_raw = ( - rab_b(ji+ip,jj   ,iku,jp_tem) * zdit + rab_b(ji+ip,jj   ,iku,jp_sal) * zdis ) * r1_e1u(ji,jj)
+                  zdyrho_raw = ( - rab_b(ji   ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji   ,jj+jp,ikv,jp_sal) * zdjs ) * r1_e2v(ji,jj)
                   zdxrho(ji+ip,jj   ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign
                   zdyrho(ji   ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw )
@@ -554,6 +441,6 @@
                      zdks = 0._wp
                   ENDIF
-                  zdzrho_raw = ( - rab_b(ji,jj,jk,jp_tem) * zdkt + rab_b(ji,jj,jk,jp_sal) * zdks ) / fse3w(ji,jj,jk+kp)
-                  zdzrho(ji,jj,jk,kp) = - MIN( - repsln, zdzrho_raw )    ! force zdzrho >= repsln
+                  zdzrho_raw = ( - zalbet(ji,jj,jk) * zdkt + zbeta0*zdks ) / fse3w(ji,jj,jk+kp)
+                  zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw )    ! force zdzrho >= repsln
                  END DO
             END DO
@@ -588,20 +475,22 @@
                   !
                   jk = nmln(ji+ip,jj) + 1
-                  IF( jk .GT. mbkt(ji+ip,jj) ) THEN  !ML reaches bottom
-                    zti_mlb(ji+ip,jj   ,1-ip,kp) = 0.0_wp
-                  ELSE
-                    ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth)
-                    zti_g_raw = (  zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp)      &
-                       &      - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj)  ) * umask(ji,jj,jk)
-                    zti_mlb(ji+ip,jj   ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw )
+                  IF( jk > mbkt(ji+ip,jj) ) THEN   ! ML reaches bottom
+                     zti_mlb(ji+ip,jj   ,1-ip,kp) = 0.0_wp
+                  ELSE                             
+                     ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth)
+                     zti_g_raw = (  zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp)      &
+                        &          - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) * r1_e1u(ji,jj)  ) * umask(ji,jj,jk)
+                     ze3_e1    =  fse3w(ji+ip,jj,jk-kp) * r1_e1u(ji,jj) 
+                     zti_mlb(ji+ip,jj   ,1-ip,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1  , ABS( zti_g_raw ) ), zti_g_raw )
                   ENDIF
                   !
                   jk = nmln(ji,jj+jp) + 1
                   IF( jk .GT. mbkt(ji,jj+jp) ) THEN  !ML reaches bottom
-                    ztj_mlb(ji   ,jj+jp,1-jp,kp) = 0.0_wp
+                     ztj_mlb(ji   ,jj+jp,1-jp,kp) = 0.0_wp
                   ELSE
-                    ztj_g_raw = (  zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp)      &
-                       &      - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj)  ) * vmask(ji,jj,jk)
-                    ztj_mlb(ji   ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw )
+                     ztj_g_raw = (  zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp)      &
+                        &      - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj)  ) * vmask(ji,jj,jk)
+                     ze3_e2    =  fse3w(ji,jj+jp,jk-kp) / e2v(ji,jj)
+                     ztj_mlb(ji   ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2  , ABS( ztj_g_raw ) ), ztj_g_raw )
                   ENDIF
                END DO
@@ -632,12 +521,16 @@
                      zti_raw   = zdxrho(ji+ip,jj   ,jk,1-ip) / zdzrho(ji+ip,jj   ,jk,kp)                   ! unmasked
                      ztj_raw   = zdyrho(ji   ,jj+jp,jk,1-jp) / zdzrho(ji   ,jj+jp,jk,kp)
-
+                     !
                      ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface
-                     zti_coord = znot_thru_surface * ( fsdept(ji+1,jj  ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj)
-                     ztj_coord = znot_thru_surface * ( fsdept(ji  ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)                  ! unmasked
+                     zti_coord = znot_thru_surface * ( fsdept(ji+1,jj  ,jk) - fsdept(ji,jj,jk) ) * r1_e1u(ji,jj)
+                     ztj_coord = znot_thru_surface * ( fsdept(ji  ,jj+1,jk) - fsdept(ji,jj,jk) ) * r1_e2v(ji,jj)     ! unmasked
                      zti_g_raw = zti_raw - zti_coord      ! ref to geopot surfaces
                      ztj_g_raw = ztj_raw - ztj_coord
-                     zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw )
-                     ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw )
+                     ! additional limit required in bilaplacian case
+                     ze3_e1    = fse3w(ji+ip,jj   ,jk+kp) * r1_e1u(ji,jj)
+                     ze3_e2    = fse3w(ji   ,jj+jp,jk+kp) * r1_e2v(ji,jj)
+                     ! NB: hard coded factor 5 (can be a namelist parameter...)
+                     zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw )
+                     ztj_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2, ABS( ztj_g_raw ) ), ztj_g_raw )
                      !
                      ! Below  ML use limited zti_g as is & mask
@@ -668,24 +561,27 @@
                      !
                      IF( ln_triad_iso ) THEN
-                        zti_raw = zti_lim**2 / zti_raw
-                        ztj_raw = ztj_lim**2 / ztj_raw
+                        zti_raw = zti_lim*zti_lim / zti_raw
+                        ztj_raw = ztj_lim*ztj_lim / ztj_raw
                         zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw )
                         ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw )
-                        zti_lim =           zfacti   * zti_lim                       &
-                        &      + ( 1._wp - zfacti ) * zti_raw
-                        ztj_lim =           zfactj   * ztj_lim                       &
-                        &      + ( 1._wp - zfactj ) * ztj_raw
+                        zti_lim = zfacti * zti_lim + ( 1._wp - zfacti ) * zti_raw
+                        ztj_lim = zfactj * ztj_lim + ( 1._wp - zfactj ) * ztj_raw
                      ENDIF
-                     triadi(ji+ip,jj   ,jk,1-ip,kp) = zti_lim
-                     triadj(ji   ,jj+jp,jk,1-jp,kp) = ztj_lim
-                    !
-                     zbu = e1u(ji    ,jj) * e2u(ji   ,jj) * fse3u(ji   ,jj,jk   )
-                     zbv = e1v(ji    ,jj) * e2v(ji   ,jj) * fse3v(ji   ,jj,jk   )
-                     zbti = e1t(ji+ip,jj) * e2t(ji+ip,jj) * fse3w(ji+ip,jj,jk+kp)
-                     zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp)
-                     !
-                     !!gm this may inhibit vectorization on Vect Computers, and even on scalar computers....  ==> to be checked
-                     wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim**2      ! masked
-                     wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim**2
+                     !                                      ! switching triad scheme 
+                     zisw = (rn_sw_triad - 1._wp ) + rn_sw_triad    &
+                        &            * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - ip ) * SIGN( 1._wp , zdxrho(ji+ip,jj,jk,1-ip) )  )
+                     zjsw = (rn_sw_triad - 1._wp ) + rn_sw_triad    &
+                        &            * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - jp ) * SIGN( 1._wp , zdyrho(ji,jj+jp,jk,1-jp) )  )
+                     !
+                     triadi(ji+ip,jj   ,jk,1-ip,kp) = zti_lim * zisw
+                     triadj(ji   ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw
+                     !
+                     zbu  = e1e2u(ji   ,jj   ) * fse3u(ji   ,jj   ,jk   )
+                     zbv  = e1e2v(ji   ,jj   ) * fse3v(ji   ,jj   ,jk   )
+                     zbti = e1e2t(ji+ip,jj   ) * fse3w(ji+ip,jj   ,jk+kp)
+                     zbtj = e1e2t(ji   ,jj+jp) * fse3w(ji   ,jj+jp,jk+kp)
+                     !
+                     wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim      ! masked
+                     wslp2(ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim*ztj_g_lim
                   END DO
                END DO
@@ -699,10 +595,11 @@
       !
       CALL wrk_dealloc( jpi,jpj, z1_mlbw )
+      CALL wrk_dealloc( jpi,jpj,jpk, zalbet )
       CALL wrk_dealloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho,              klstart = 0  )
       CALL wrk_dealloc( jpi,jpj,  2,2, zti_mlb, ztj_mlb,        kkstart = 0, klstart = 0  )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('ldf_slp_grif')
-      !
-   END SUBROUTINE ldf_slp_grif
+      IF( nn_timing == 1 )  CALL timing_stop('ldf_slp_triad')
+      !
+   END SUBROUTINE ldf_slp_triad
 
 
@@ -730,5 +627,5 @@
       INTEGER  ::   ji , jj , jk                   ! dummy loop indices
       INTEGER  ::   iku, ikv, ik, ikm1             ! local integers
-      REAL(wp) ::   zeps, zm1_g, zm1_2g            ! local scalars
+      REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_slpmax ! local scalars
       REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      -
       REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      -
@@ -741,4 +638,5 @@
       zm1_g  = -1.0_wp / grav
       zm1_2g = -0.5_wp / grav
+      z1_slpmax = 1._wp / rn_slpmax
       !
       uslpml (1,:) = 0._wp      ;      uslpml (jpi,:) = 0._wp
@@ -752,8 +650,6 @@
             DO ji = 1, jpi
                ik = nmln(ji,jj) - 1
-               IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN
-                  omlmask(ji,jj,jk) = 1._wp
-               ELSE
-                  omlmask(ji,jj,jk) = 0._wp
+               IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1._wp
+               ELSE                  ;   omlmask(ji,jj,jk) = 0._wp
                ENDIF
             END DO
@@ -777,15 +673,15 @@
             !
             !                        !- vertical density gradient for u- and v-slopes (from dzr at T-point)
-            iku = MIN(  MAX( miku(ji,jj)+1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1)
-            ikv = MIN(  MAX( mikv(ji,jj)+1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   !
+            iku = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1)
+            ikv = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   !
             zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj  ,iku) )
             zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji  ,jj+1,ikv) )
             !                        !- horizontal density gradient at u- & v-points
-            zau = p_gru(ji,jj,iku) / e1u(ji,jj)
-            zav = p_grv(ji,jj,ikv) / e2v(ji,jj)
+            zau = p_gru(ji,jj,iku) * r1_e1u(ji,jj)
+            zav = p_grv(ji,jj,ikv) * r1_e2v(ji,jj)
             !                        !- 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,iku)* ABS( zau )  )
-            zbv = MIN(  zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav )  )
+            zbu = MIN(  zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau )  )
+            zbv = MIN(  zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav )  )
             !                        !- Slope at u- & v-points (uslpml, vslpml)
             uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku)
@@ -812,6 +708,6 @@
             zbj = MIN(  zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj )  )
             !                        !- i- & j-slope at w-points (wslpiml, wslpjml)
-            wslpiml(ji,jj) = zai / ( zbi - zeps ) * wmask (ji,jj,ik)
-            wslpjml(ji,jj) = zaj / ( zbj - zeps ) * wmask (ji,jj,ik)
+            wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik)
+            wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik)
          END DO
       END DO
@@ -831,6 +727,5 @@
       !! ** Purpose :   Initialization for the isopycnal slopes computation
       !!
-      !! ** Method  :   read the nammbf namelist and check the parameter
-      !!      values called by tra_dmp at the first timestep (nit000)
+      !! ** Method  :   
       !!----------------------------------------------------------------------
       INTEGER ::   ji, jj, jk   ! dummy loop indices
@@ -845,15 +740,21 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
       ENDIF
-
-      IF( ln_traldf_grif ) THEN        ! Griffies operator : triad of slopes
-         ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , wslp2(jpi,jpj,jpk) , STAT=ierr )
-         ALLOCATE( triadi  (jpi,jpj,jpk,0:1,0:1) , triadj  (jpi,jpj,jpk,0:1,0:1)                      , STAT=ierr )
-         IF( ierr > 0             )   CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' )
-         !
+      !
+      ALLOCATE( ah_wslp2(jpi,jpj,jpk) , akz(jpi,jpj,jpk) , STAT=ierr )
+      IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate ah_slp2 or akz' )
+      !
+      IF( ln_traldf_triad ) THEN        ! Griffies operator : triad of slopes
+         IF(lwp) WRITE(numout,*) '              Griffies (triad) operator initialisation'
+         ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) ,     &
+            &      triadi  (jpi,jpj,jpk,0:1,0:1) , triadj  (jpi,jpj,jpk,0:1,0:1) ,     &
+            &      wslp2   (jpi,jpj,jpk)                                         , STAT=ierr )
+         IF( ierr > 0      )   CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' )
          IF( ln_dynldf_iso )   CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' )
          !
       ELSE                             ! Madec operator : slopes at u-, v-, and w-points
-         ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) ,                &
-            &   omlmask(jpi,jpj,jpk) , uslpml(jpi,jpj)   , vslpml(jpi,jpj)    , wslpiml(jpi,jpj)   , wslpjml(jpi,jpj) , STAT=ierr )
+         IF(lwp) WRITE(numout,*) '              Madec operator initialisation'
+         ALLOCATE( omlmask(jpi,jpj,jpk) ,                                                                        &
+            &      uslp(jpi,jpj,jpk) , uslpml(jpi,jpj) , wslpi(jpi,jpj,jpk) , wslpiml(jpi,jpj) ,     &
+            &      vslp(jpi,jpj,jpk) , vslpml(jpi,jpj) , wslpj(jpi,jpj,jpk) , wslpjml(jpi,jpj) , STAT=ierr )
          IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Madec operator slope ' )
 
@@ -865,26 +766,27 @@
          wslpj(:,:,:) = 0._wp   ;   wslpjml(:,:) = 0._wp
 
-         IF(ln_sco .AND.  (ln_traldf_hor .OR. ln_dynldf_hor )) THEN
-            IF(lwp)   WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
-
-            ! geopotential diffusion in s-coordinates on tracers and/or momentum
-            ! The slopes of s-surfaces are computed once (no call to ldfslp in step)
-            ! 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
-            CALL lbc_lnk( uslp , 'U', -1. )   ;   CALL lbc_lnk( vslp , 'V', -1. )      ! Lateral boundary conditions
-            CALL lbc_lnk( wslpi, 'W', -1. )   ;   CALL lbc_lnk( wslpj, 'W', -1. )
-         ENDIF
+         !!gm I no longer understand this.....
+!!gm         IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN
+!            IF(lwp)   WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
+!
+!            ! geopotential diffusion in s-coordinates on tracers and/or momentum
+!            ! The slopes of s-surfaces are computed once (no call to ldfslp in step)
+!            ! 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) = - ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
+!                     vslp (ji,jj,jk) = - ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
+!                     wslpi(ji,jj,jk) = - ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * r1_e1t(ji,jj) * wmask(ji,jj,jk) * 0.5
+!                     wslpj(ji,jj,jk) = - ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * r1_e2t(ji,jj) * wmask(ji,jj,jk) * 0.5
+!                  END DO
+!               END DO
+!            END DO
+!            CALL lbc_lnk( uslp , 'U', -1. )   ;   CALL lbc_lnk( vslp , 'V', -1. )      ! Lateral boundary conditions
+!            CALL lbc_lnk( wslpi, 'W', -1. )   ;   CALL lbc_lnk( wslpj, 'W', -1. )
+!!gm         ENDIF
       ENDIF
       !
@@ -892,23 +794,4 @@
       !
    END SUBROUTINE ldf_slp_init
-
-#else
-   !!------------------------------------------------------------------------
-   !!   Dummy module :                 NO Rotation of lateral mixing tensor
-   !!------------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_ldfslp = .FALSE.    !: slopes flag
-CONTAINS
-   SUBROUTINE ldf_slp( kt, prd, pn2 )   ! Dummy routine
-      INTEGER, INTENT(in) :: kt
-      REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2
-      WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1)
-   END SUBROUTINE ldf_slp
-   SUBROUTINE ldf_slp_grif( kt )        ! Dummy routine
-      INTEGER, INTENT(in) :: kt
-      WRITE(*,*) 'ldf_slp_grif: You should not have seen this print! error?', kt
-   END SUBROUTINE ldf_slp_grif
-   SUBROUTINE ldf_slp_init              ! Dummy routine
-   END SUBROUTINE ldf_slp_init
-#endif
 
    !!======================================================================
