Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 2400)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 2401)
@@ -10,4 +10,5 @@
    !!            1.0  ! 2005-10  (A. Beckmann)  correction for s-coordinates
    !!            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
    !!----------------------------------------------------------------------
 #if   defined key_ldfslp   ||   defined key_esopa
@@ -93,5 +94,5 @@
       USE oce , zgru  => ua   ! use ua as workspace
       USE oce , zgrv  => va   ! use va as workspace
-      USE oce , zwy   => ta   ! use ta as workspace
+      USE oce , zww   => ta   ! use ta as workspace
       USE oce , zwz   => sa   ! use sa as workspace
       !!
@@ -103,20 +104,20 @@
       INTEGER  ::   ii0, ii1, iku   ! temporary integer
       INTEGER  ::   ij0, ij1, ikv   ! temporary integer
-      REAL(wp) ::   zeps, zmg, zm05g, zalpha        ! temporary scalars
-      REAL(wp) ::   zcoef1, zcoef2, zcoef3          !    -         -
-      REAL(wp) ::   zcofu , zcofv , zcofw           !    -         -
-      REAL(wp) ::   zau, zbu, zai, zbi, z1u, z1wu   !    -         -
-      REAL(wp) ::   zav, zbv, zaj, zbj, z1v, z1wv   !
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zww     ! 3D workspace
+      REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16     ! 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), DIMENSION(jpi,jpj,jpk) ::   zdzr   ! 3D workspace
       !!----------------------------------------------------------------------
       
-      zeps  =  1.e-20                             ! Local constant initialization
-      zmg   = -1.0 / grav
-      zm05g = -0.5 / grav
+      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.e0
       zwz(:,:,:) = 0.e0
-      !                                           ! horizontal density gradient computation
-      DO jk = 1, jpk
+      !
+      DO jk = 1, jpk             !==   i- & j-gradient of density   ==!
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
@@ -141,6 +142,18 @@
          END DO
       ENDIF
-      
-      CALL ldf_slp_mxl( prd, pn2 )                ! Slopes of isopycnal surfaces just below the mixed layer
+      !
+      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
 
       
@@ -148,36 +161,36 @@
       ! ===========================      | vslp = d/dj( prd ) / d/dz( prd )
       !               
-      !                                           !* Local vertical density gradient evaluated from N^2
-      DO jk = 2, jpkm1                            !  zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zwy(ji,jj,jk) = zmg * ( prd(ji,jj,jk) + 1. ) * (    pn2  (ji,jj,jk) + pn2  (ji,jj,jk+1)     )   &
-                  &                                         / MAX( tmask(ji,jj,jk) + tmask(ji,jj,jk+1), 1. )
-            END DO
-         END DO
-      END DO
       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 = 1. / e1u(ji,jj) * zgru(ji,jj,jk)
-               zav = 1. / e2v(ji,jj) * zgrv(ji,jj,jk)
-               zbu = 0.5 * ( zwy(ji,jj,jk) + zwy(ji+1,jj  ,jk) )
-               zbv = 0.5 * ( zwy(ji,jj,jk) + zwy(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.*ABS( zau ), -7.e+3/fse3u(ji,jj,jk)*ABS( zau ) )
-               zbv = MIN( zbv, -100.*ABS( zav ), -7.e+3/fse3v(ji,jj,jk)*ABS( zav ) )
-               ! uslp and vslp output in zwz and zww, resp.
-               zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
-               zwz (ji,jj,jk) = ( ( 1. - zalpha) * zau / ( zbu - zeps )                                           &
-                  &                    + zalpha  * uslpml(ji,jj)                                                  &
-                  &                              * 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) )   &
-                  &                              / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5. ) ) * umask(ji,jj,jk)
-               zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )
-               zww (ji,jj,jk) = ( ( 1. - zalpha) * zav / ( zbv - zeps )                                           &
-                  &                    + zalpha  * vslpml(ji,jj)                                                  &
-                  &                              * 0.5 * ( 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)
+               !                                      ! 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
@@ -185,15 +198,14 @@
       CALL lbc_lnk( zwz, 'U', -1. )   ;   CALL lbc_lnk( zww, 'V', -1. )      ! lateral boundary conditions
       !
-      zcofu = 1. / 16.                            !* horizontal Shapiro filter
-      zcofv = 1. / 16.
+      !                                            !* horizontal Shapiro filter
       DO jk = 2, jpkm1
          DO jj = 2, jpjm1, jpj-3                        ! rows jj=2 and =jpjm1 only
             DO ji = 2, jpim1  
-               uslp(ji,jj,jk) = zcofu * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
+               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) = zcofv * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,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)      &
@@ -204,10 +216,10 @@
          DO jj = 3, jpj-2                               ! other rows
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               uslp(ji,jj,jk) = zcofu * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
+               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) = zcofv * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,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)      &
@@ -219,10 +231,8 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               z1u  = ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk) )*.5
-               z1v  = ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk) )*.5
-               z1wu = ( umask(ji,jj,jk)   + umask(ji,jj,jk+1) )*.5
-               z1wv = ( vmask(ji,jj,jk)   + vmask(ji,jj,jk+1) )*.5
-               uslp(ji,jj,jk) = uslp(ji,jj,jk) * z1u * z1wu
-               vslp(ji,jj,jk) = vslp(ji,jj,jk) * z1v * z1wv
+               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
@@ -233,42 +243,40 @@
       ! ===========================      | wslpj = mij( d/dj( prd ) / d/dz( prd )
       !               
-      !                                           !* Local vertical density gradient evaluated from N^2
-      DO jk = 2, jpkm1                            !  zwy = d/dz(prd)= - mk ( prd ) / grav * pn2 -- at w point
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zwy(ji,jj,jk) = zm05g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. )
-            END DO
-         END DO
-      END DO
-      DO jk = 2, jpkm1                            !* Slopes at w point
+      DO jk = 2, jpkm1
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               !                                        ! horizontal density i-gradient at w-points
-               zcoef1 = MAX( zeps, umask(ji-1,jj,jk  )+umask(ji,jj,jk  )    &
-                  &               +umask(ji-1,jj,jk-1)+umask(ji,jj,jk-1) )
-               zcoef1 = 1. / ( zcoef1 * e1t (ji,jj) )
-               zai = zcoef1 * (  zgru(ji  ,jj,jk  ) + zgru(ji  ,jj,jk-1)   &
-                  &            + zgru(ji-1,jj,jk-1) + zgru(ji-1,jj,jk  ) ) * tmask (ji,jj,jk)
-               !                                        ! horizontal density j-gradient at w-points
-               zcoef2 = MAX( zeps, vmask(ji,jj-1,jk  )+vmask(ji,jj,jk-1)   &
-                  &               +vmask(ji,jj-1,jk-1)+vmask(ji,jj,jk  ) )
-               zcoef2 = 1.0 / ( zcoef2 *  e2t (ji,jj) )
-               zaj = zcoef2 * (  zgrv(ji,jj  ,jk  ) + zgrv(ji,jj  ,jk-1)   &
-                  &            + zgrv(ji,jj-1,jk-1) + zgrv(ji,jj-1,jk  ) ) * tmask (ji,jj,jk)
+               !                                  !* 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.
-               !                                        ! static instability: kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
-               zbi = MIN( zwy (ji,jj,jk),- 100.*ABS(zai), -7.e+3/fse3w(ji,jj,jk)*ABS(zai) )
-               zbj = MIN( zwy (ji,jj,jk), -100.*ABS(zaj), -7.e+3/fse3w(ji,jj,jk)*ABS(zaj) )
-               !                                        ! wslpi and wslpj output in zwz and zww, resp.
-               zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )
-               zcoef3 = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. )
-               zwz(ji,jj,jk) = (     zai / ( zbi - zeps)  * ( 1. - zalpha )   &
-                  &             + zcoef3 * wslpiml(ji,jj) *        zalpha   ) * tmask (ji,jj,jk)
-               zww(ji,jj,jk) = (     zaj / ( zbj - zeps)  * ( 1. - zalpha )   &
-                  &             + zcoef3 * wslpjml(ji,jj) *        zalpha   ) * tmask (ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( zwz, 'T', -1. )   ;    CALL lbc_lnk( zww, 'T', -1. )      ! lateral boundary conditions on zwz and zww
+               !                                        ! + 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
@@ -276,10 +284,9 @@
          DO jj = 2, jpjm1, jpj-3                        ! rows jj=2 and =jpjm1
             DO ji = 2, jpim1
-               zcofw = tmask(ji,jj,jk) / 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
+                  &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk)
 
                wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
@@ -287,15 +294,14 @@
                   &                + 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
+                  &                + 4.*  zww(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk)
             END DO
          END DO  
          DO jj = 3, jpj-2                               ! other rows
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zcofw = tmask(ji,jj,jk) / 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
+                  &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk)
 
                wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
@@ -303,5 +309,5 @@
                   &                + 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
+                  &                + 4.*  zww(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk)
             END DO
          END DO
@@ -309,8 +315,8 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               z1u = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) *.5
-               z1v = ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) *.5
-               wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * z1u * z1v
-               wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * z1u * z1v
+               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
@@ -573,5 +579,5 @@
 
 
-   SUBROUTINE ldf_slp_mxl( prd, pn2 )
+   SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_slp_mxl  ***
@@ -580,35 +586,29 @@
       !!              the mixed layer.
       !!
-      !! ** Method :
-      !!      The slope in the i-direction is computed at u- and w-points
-      !!   (uslp, wslpi) and the slope in the j-direction is computed at
-      !!   v- and w-points (vslp, wslpj).
-      !!   They are bounded by 1/100 over the whole ocean, and within the
-      !!   surface layer they are bounded by the distance to the surface
-      !!   ( slope<= depth/l  where l is the length scale of horizontal
-      !!   diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity
-      !!   of 10cm/s)
-      !!
-      !! ** Action :   Compute uslp, wslpi, and vslp, wslpj, the i- and  j-slopes
-      !!   of now neutral surfaces at u-, w- and v- w-points, resp.
-      !!----------------------------------------------------------------------
-      USE oce , zgru  => ua   ! ua, va used as workspace and set to hor. 
-      USE oce , zgrv  => va   ! density gradient in ldf_slp
-      !! 
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   prd   ! in situ density
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   pn2   ! Brunt-Vaisala frequency (locally ref.)
-      !!
-      INTEGER  ::   ji, jj, jk           ! dummy loop indices
-      INTEGER  ::   ik, ikm1             ! temporary integers
-      REAL(wp) ::   zeps, zmg, zm05g     ! temporary scalars
-      REAL(wp) ::   zcoef1, zcoef2       !    -         -
-      REAL(wp) ::   zau, zbu, zai, zbi   !    -         -
-      REAL(wp) ::   zav, zbv, zaj, zbj   !    -         -
-      REAL(wp), DIMENSION(jpi,jpj) ::   zwy   ! 2D workspace
-      !!----------------------------------------------------------------------
-
-      zeps  =  1.e-20                    ! Local constant initialization
-      zmg   = -1.0 / grav
-      zm05g = -0.5 / grav
+      !! ** Method  :   The slope in the i-direction is computed at u- & w-points
+      !!              (uslpml, wslpiml) and the slope in the j-direction is computed
+      !!              at v- and w-points (vslpml, wslpjml) with the same bounds as
+      !!              in ldf_slp.
+      !!
+      !! ** Action  :   uslpml, wslpiml :  i- &  j-slopes of neutral surfaces
+      !!                vslpml, wslpjml    just below the mixed layer 
+      !!                omlmask         :  mixed layer mask
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   prd            ! in situ density
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   pn2            ! Brunt-Vaisala frequency (locally ref.)
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   p_gru, p_grv   ! i- & j-gradient of density (u- & v-pts)
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   p_dzr          ! z-gradient of density      (T-point)
+      !!
+      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) ::   zci, zfi, zau, zbu, zai, zbi   !   -      -
+      REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      -
+      REAL(wp) ::   zck, zfk,      zbw             !   -      -
+      !!----------------------------------------------------------------------
+
+      zeps   =  1.e-20_wp        !==   Local constant initialization   ==!
+      zm1_g  = -1.0_wp / grav
+      zm1_2g = -0.5_wp / grav
       !
       uslpml (1,:) = 0.e0      ;      uslpml (jpi,:) = 0.e0
@@ -616,7 +616,7 @@
       wslpiml(1,:) = 0.e0      ;      wslpiml(jpi,:) = 0.e0
       wslpjml(1,:) = 0.e0      ;      wslpjml(jpi,:) = 0.e0
-
-      !                                  ! surface mixed layer mask
-      DO jk = 1, jpk                          ! =1 inside the mixed layer, =0 otherwise
+      !
+      !                          !==   surface mixed layer mask   !
+      DO jk = 1, jpk                      ! =1 inside the mixed layer, =0 otherwise
 # if defined key_vectopt_loop
          DO jj = 1, 1
@@ -645,22 +645,4 @@
       !-----------------------------------------------------------------------
       !
-      zwy(:,jpj) = 0.e0            !* vertical density gradient for u-slope (from N^2)
-      zwy(jpi,:) = 0.e0
-# 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
-            ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) )       ! avoid spurious recirculation 
-            ik = MIN( ik, jpkm1 )                            ! if ik = jpk take jpkm1 values
-            zwy(ji,jj) = zmg * ( prd(ji,jj,ik) + 1. ) * (    pn2  (ji,jj,ik) + pn2  (ji,jj,ik+1)     )   &
-               &                                      / MAX( tmask(ji,jj,ik) + tmask(ji,jj,ik+1), 1. )
-         END DO
-      END DO
-      CALL lbc_lnk( zwy, 'U', 1. )      ! lateral boundary conditions NO sign change
-
-      !                            !* Slope at u points
 # if defined key_vectopt_loop
       DO jj = 1, 1
@@ -670,107 +652,49 @@
          DO ji = 2, jpim1
 # endif
-            ! horizontal and vertical density gradient at u-points
-            ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) )
-            ik = MIN( ik, jpkm1 )
-            zau = 1./ e1u(ji,jj) * zgru(ji,jj,ik)
-            zbu = 0.5*( zwy(ji,jj) + zwy(ji+1,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.*ABS(zau), -7.e+3/fse3u(ji,jj,ik)*ABS(zau) )
-            ! uslpml
-            uslpml (ji,jj) = zau / ( zbu - zeps ) * umask (ji,jj,ik)
-         END DO
-      END DO
-      CALL lbc_lnk( uslpml, 'U', -1. )      ! lateral boundary conditions (i-gradient => sign change)
-
-      !                            !* vertical density gradient for v-slope (from N^2)
-# 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
-            ik = MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) )
-            ik = MIN( ik, jpkm1 )
-            zwy(ji,jj) = zmg * ( prd(ji,jj,ik) + 1. ) * (    pn2  (ji,jj,ik) + pn2  (ji,jj,ik+1)     )   &
-               &                                      / MAX( tmask(ji,jj,ik) + tmask(ji,jj,ik+1), 1. )
-         END DO
-      END DO
-      CALL lbc_lnk( zwy, 'V', 1. )      ! lateral boundary conditions NO sign change
-
-      !                            !* Slope at v points
-# 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
-            ! horizontal and vertical density gradient at v-points
-            ik = MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) )
-            ik = MIN( ik,jpkm1 )
-            zav = 1./ e2v(ji,jj) * zgrv(ji,jj,ik)
-            zbv = 0.5*( zwy(ji,jj) + zwy(ji,jj+1) )
-            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
-            !                         kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
-            zbv = MIN( zbv, -100.*ABS(zav), -7.e+3/fse3v(ji,jj,ik)*ABS( zav ) )
-            ! vslpml
-            vslpml (ji,jj) = zav / ( zbv - zeps ) * vmask (ji,jj,ik)
-         END DO
-      END DO
-      CALL lbc_lnk( vslpml, 'V', -1. )      ! lateral boundary conditions (j-gradient => sign change)
-
-
-      !                            !* vertical density gradient for w-slope (from N^2)
-# 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
-            ik = MIN( ik, jpk )
-            ikm1 = MAX ( 1, ik-1)
-            zwy (ji,jj) = zm05g * pn2 (ji,jj,ik) *     &
-               &             ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. )
-         END DO
-      END DO
-
-      !                            !* Slopes at w points
-# 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
-            ik = nmln(ji,jj) + 1
-            ik = MIN( ik, jpk )
-            ikm1 = MAX ( 1, ik-1 )
-            ! horizontal density i-gradient at w-points
-            zcoef1 = MAX( zeps, umask(ji-1,jj,ik  )+umask(ji,jj,ik  )   &
-               &               +umask(ji-1,jj,ikm1)+umask(ji,jj,ikm1) )
-            zcoef1 = 1. / ( zcoef1 * e1t (ji,jj) )
-            zai = zcoef1 * (  zgru(ji  ,jj,ik  ) + zgru(ji  ,jj,ikm1)   &
-               &            + zgru(ji-1,jj,ikm1) + zgru(ji-1,jj,ik  ) ) * tmask (ji,jj,ik)
-            ! horizontal density j-gradient at w-points
-            zcoef2 = MAX( zeps, vmask(ji,jj-1,ik  )+vmask(ji,jj,ikm1)    &
-               &               +vmask(ji,jj-1,ikm1)+vmask(ji,jj,ik  ) )
-            zcoef2 = 1.0 / ( zcoef2 *  e2t (ji,jj) )
-            zaj = zcoef2 * (  zgrv(ji,jj  ,ik  ) + zgrv(ji,jj  ,ikm1)   &
-               &            + zgrv(ji,jj-1,ikm1) + zgrv(ji,jj-1,ik  ) ) * tmask (ji,jj,ik)
-            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
-            !                   static instability:
-            !                   kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
-            zbi = MIN ( zwy (ji,jj),- 100.*ABS(zai), -7.e+3/fse3w(ji,jj,ik)*ABS(zai) )
-            zbj = MIN ( zwy (ji,jj), -100.*ABS(zaj), -7.e+3/fse3w(ji,jj,ik)*ABS(zaj) )
-            ! wslpiml and wslpjml
-            wslpiml (ji,jj) = zai / ( zbi - zeps) * tmask (ji,jj,ik)
-            wslpjml (ji,jj) = zaj / ( zbj - zeps) * tmask (ji,jj,ik)
-         END DO
-      END DO
-      CALL lbc_lnk( wslpiml, 'W', -1. )   ;   CALL lbc_lnk( wslpjml, 'W', -1. )      ! lateral boundary conditions
+            !                    !==   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  )   ! 
+            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)
+            !                          !- 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,ik)* ABS( zau )  )
+            zbv = MIN(  zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ik)* ABS( zav )  )
+            !                          !- Slope at u- & v-points (uslpml, vslpml)
+            uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,ik)
+            vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ik)
+            !
+            !                    !==   i- & j-slopes at w-points just below the Mixed Layer   ==!
+            !
+            ik   = MIN( nmln(ji,jj) + 1, jpk )
+            ikm1 = MAX( 1, ik-1 )
+            !                          !- vertical density gradient for w-slope (from N^2)
+            zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. )
+            !                          !- horizontal density i- & j-gradient at w-points
+            zci = MAX(   umask(ji-1,jj,ik  ) + umask(ji,jj,ik  )           &
+               &       + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps  ) * e1t(ji,jj) 
+            zcj = MAX(   vmask(ji,jj-1,ik  ) + vmask(ji,jj,ik  )           &
+               &       + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps  ) * e2t(ji,jj)
+            zai =    (   p_gru(ji-1,jj,ik  ) + p_gru(ji,jj,ik)           &
+               &       + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1  )  ) / zci  * tmask(ji,jj,ik)
+            zaj =    (   p_grv(ji,jj-1,ik  ) + p_grv(ji,jj,ik  )           &
+               &       + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1)  ) / zcj  * tmask(ji,jj,ik)
+            !                          !- 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,ik)* ABS( zai )  )
+            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)
+         END DO
+      END DO
+!!gm this lbc_lnk should be useless....
+      CALL lbc_lnk( uslpml , 'U', -1. )   ;   CALL lbc_lnk( vslpml , 'V', -1. )   ! lateral boundary cond. (sign change)
+      CALL lbc_lnk( wslpiml, 'W', -1. )   ;   CALL lbc_lnk( wslpjml, 'W', -1. )   ! lateral boundary conditions
       !
    END SUBROUTINE ldf_slp_mxl
