Index: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 4744)
+++ branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 5260)
@@ -28,6 +28,7 @@
    USE zdfmxl         ! mixed layer depth
    USE eosbn2         ! equation of states
+   !
+   USE in_out_manager ! I/O manager
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE in_out_manager ! I/O manager
    USE prtctl         ! Print control
    USE wrk_nemo       ! work arrays
@@ -108,7 +109,9 @@
       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
       !!----------------------------------------------------------------------
       !
@@ -116,4 +119,5 @@
       !
       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 
@@ -136,11 +140,6 @@
          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)
-# else
             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)
@@ -148,6 +147,15 @@
             END DO
          ENDIF
-         !
-         zdzr(:,:,1) = 0._wp        !==   Local vertical density gradient at T-point   == !   (evaluated from N^2)
+         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
@@ -159,4 +167,14 @@
                &                 * ( 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   ==!
@@ -167,4 +185,23 @@
          ! ===========================      | 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
@@ -180,14 +217,20 @@
                   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)
+                  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)
+                  
+                 
 !!gm  modif to suppress omlmask.... (as in Griffies case)
 !                  !                                         ! jk must be >= ML level for zf=1. otherwise  zf=0.
@@ -238,7 +281,9 @@
                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) + 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) + vmask(ji  ,jj,jk+1) ) * 0.5_wp   &
+                     &                            *   vmask(ji,jj,jk-1)
                END DO
             END DO
@@ -253,5 +298,5 @@
                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. )
+                  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
@@ -261,7 +306,7 @@
                      &      + 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)
+                     &      + 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 * tmask (ji,jj,jk)
+                     &      + 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)
@@ -270,7 +315,9 @@
                   !                                        ! 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)
+                  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)
 
 !!gm  modif to suppress omlmask....  (as in Griffies operator)
@@ -325,6 +372,6 @@
                   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
+                  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
@@ -377,11 +424,22 @@
          ! 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,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) ) * 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 
+                  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 
@@ -405,4 +463,5 @@
       
       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')
@@ -435,7 +494,5 @@
       REAL(wp) ::   zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_g_raw, ztj_g_lim
       REAL(wp) ::   zdzrho_raw
-      REAL(wp) ::   zbeta0
       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
@@ -445,5 +502,4 @@
       !
       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  )
@@ -452,6 +508,4 @@
       !  Some preliminary calculation  !
       !--------------------------------!
-      !
-      CALL eos_alpbet( tsb, zalbet, zbeta0 )  !==  before local thermal/haline expension ratio at T-points  ==!
       !
       DO jl = 0, 1                            !==  unmasked before density i- j-, k-gradients  ==!
@@ -465,7 +519,7 @@
                   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 = ( - zalbet(ji+ip,jj   ,jk) * zdit + zbeta0*zdis ) / e1u(ji,jj)
-                  zdyrho_raw = ( - zalbet(ji   ,jj+jp,jk) * zdjt + zbeta0*zdjs ) / e2v(ji,jj)
-                  zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN( MAX(   repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign
+                  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 )
                END DO
@@ -473,17 +527,12 @@
          END DO
          !
-         IF( ln_zps.and.l_grad_zps ) THEN     ! partial steps: correction of i- & j-grad on bottom
-# if defined key_vectopt_loop
-            DO jj = 1, 1
-               DO ji = 1, jpij-jpi            ! vector opt. (forced unrolling)
-# else
+         IF( ln_zps .AND. l_grad_zps ) THEN     ! partial steps: correction of i- & j-grad on bottom
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
-# endif
                   iku  = mbku(ji,jj)          ;   ikv  = mbkv(ji,jj)             ! last ocean level (u- & v-points)
                   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 = ( - zalbet(ji+ip,jj   ,iku) * zdit + zbeta0*zdis ) / e1u(ji,jj)
-                  zdyrho_raw = ( - zalbet(ji   ,jj+jp,ikv) * zdjt + zbeta0*zdjs ) / e2v(ji,jj)
+                  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(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 )
@@ -505,6 +554,6 @@
                      zdks = 0._wp
                   ENDIF
-                  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
+                  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
                  END DO
             END DO
@@ -650,5 +699,4 @@
       !
       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  )
@@ -701,14 +749,11 @@
       !                                            !==   surface mixed layer mask   !
       DO jk = 1, jpk                               ! =1 inside the mixed layer, =0 otherwise
-# if defined key_vectopt_loop
-         DO jj = 1, 1
-            DO ji = 1, jpij                        ! vector opt. (forced unrolling)
-# else
          DO jj = 1, jpj
             DO ji = 1, jpi
-# endif
                ik = nmln(ji,jj) - 1
-               IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1._wp
-               ELSE                  ;   omlmask(ji,jj,jk) = 0._wp
+               IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN
+                  omlmask(ji,jj,jk) = 1._wp
+               ELSE
+                  omlmask(ji,jj,jk) = 0._wp
                ENDIF
             END DO
@@ -727,16 +772,11 @@
       !-----------------------------------------------------------------------
       !
-# if defined key_vectopt_loop
-      DO jj = 1, 1
-         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
-# else
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
-# endif
             !                        !==   Slope at u- & v-points just below the Mixed Layer   ==!
             !
             !                        !- vertical density gradient for u- and v-slopes (from dzr at T-point)
-            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  )   !
+            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  )   !
             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) )
@@ -772,6 +812,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 ) * tmask (ji,jj,ik)
-            wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik)
+            wslpiml(ji,jj) = zai / ( zbi - zeps ) * wmask (ji,jj,ik)
+            wslpjml(ji,jj) = zaj / ( zbj - zeps ) * wmask (ji,jj,ik)
          END DO
       END DO
@@ -825,5 +865,5 @@
          wslpj(:,:,:) = 0._wp   ;   wslpjml(:,:) = 0._wp
 
-         IF( ln_traldf_hor .OR. ln_dynldf_hor ) THEN
+         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'
 
