Index: trunk/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 1514)
+++ trunk/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 1515)
@@ -4,4 +4,10 @@
    !! Ocean physics: slopes of neutral surfaces
    !!======================================================================
+   !! History :  OPA  ! 1994-12  (G. Madec, M. Imbard)  Original code
+   !!            8.0  ! 1997-06  (G. Madec)  optimization, lbc
+   !!            8.1  ! 1999-10  (A. Jouzeau)  NEW profile in the mixed layer
+   !!   NEMO     0.5  ! 2002-10  (G. Madec)  Free form, F90
+   !!            1.0  ! 2005-10  (A. Beckmann)  correction for s-coordinates
+   !!----------------------------------------------------------------------
 #if   defined key_ldfslp   ||   defined key_esopa
    !!----------------------------------------------------------------------
@@ -12,5 +18,4 @@
    !!   ldf_slp_init : initialization of the slopes computation
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -26,21 +31,13 @@
    PRIVATE
 
-   !! * Accessibility
-   PUBLIC ldf_slp        ! routine called by step.F90
-
-   !! * Share module variables
-   LOGICAL , PUBLIC, PARAMETER ::   lk_ldfslp = .TRUE.     !: slopes flag
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !:
-      uslp, wslpi,         &  !: i_slope at U- and W-points
-      vslp, wslpj             !: j-slope at V- and W-points
+   PUBLIC   ldf_slp    ! routine called by step.F90
+
+   LOGICAL , PUBLIC, PARAMETER              ::   lk_ldfslp = .TRUE.   !: slopes flag
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   uslp, wslpi          !: i_slope at U- and W-points
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   vslp, wslpj          !: j-slope at V- and W-points
    
-   !! * Module variables
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-      omlmask                 ! mask of the surface mixed layer at T-pt
-   REAL(wp), DIMENSION(jpi,jpj) ::   &
-      uslpml, wslpiml,     &  ! i_slope at U- and W-points just below
-   !                          ! the surface mixed layer
-      vslpml, wslpjml         ! j_slope at V- and W-points just below
-   !                          ! the surface mixed layer
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   omlmask           ! mask of the surface mixed layer at T-pt
+   REAL(wp), DIMENSION(jpi,jpj)     ::   uslpml, wslpiml   ! i_slope at U- and W-points just below the mixed layer
+   REAL(wp), DIMENSION(jpi,jpj)     ::   vslpml, wslpjml   ! j_slope at V- and W-points just below the mixed layer
 
    !! * Substitutions
@@ -48,5 +45,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -59,7 +56,7 @@
       !!                 ***  ROUTINE ldf_slp  ***
       !! 
-      !! ** Purpose :   Compute the slopes of neutral surface (slope of iso-
-      !!      pycnal surfaces referenced locally) ('key_traldfiso').
-      !!
+      !! ** Purpose :   Compute the slopes of neutral surface (slope of isopycnal
+      !!              surfaces referenced locally) ('key_traldfiso').
+      !! 
       !! ** Method  :   The slope in the i-direction is computed at U- and 
       !!      W-points (uslp, wslpi) and the slope in the j-direction is 
@@ -79,56 +76,34 @@
       !! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and  j-slopes 
       !!               of now neutral surfaces at u-, w- and v- w-points, resp.
-      !!
-      !! History :
-      !!   7.0  !  94-12  (G. Madec, M. Imbard)  Original code
-      !!   8.0  !  97-06  (G. Madec)  optimization, lbc
-      !!   8.1  !  99-10  (A. Jouzeau)  NEW profile
-      !!   8.5  !  99-10  (G. Madec)  Free form, F90
-      !!   9.0  !  05-10  (A. Beckmann)  correction for s-coordinates
-      !!----------------------------------------------------------------------
-      !! * Modules used
-      USE oce            , zgru  => ua,  &  ! use ua as workspace
-                           zgrv  => va,  &  ! use va as workspace
-                           zwy   => ta,  &  ! use ta as workspace
-                           zwz   => sa      ! use sa as workspace
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   &
-         prd,                            &  ! in situ density
-         pn2                                ! Brunt-Vaisala frequency (locally ref.)
-
-      !! * Local declarations
-      INTEGER  ::   ji, jj, jk              ! dummy loop indices
-      INTEGER  ::   ii0, ii1, ij0, ij1,  &  ! temporary integer
-         &          iku, ikv                !    "          "
-      REAL(wp) ::   &
-         zeps, zmg, zm05g,               &  ! temporary scalars
-         zcoef1, zcoef2, zcoef3,         &  !
-         zau, zbu, zav, zbv,             &
-         zai, zbi, zaj, zbj,             &
-         zcofu, zcofv, zcofw,            &
-         z1u, z1v, z1wu, z1wv,           &
-         zalpha
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zww
-      !!----------------------------------------------------------------------
-
+      !!----------------------------------------------------------------------
+      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 , zwz   => sa   ! use sa as workspace
+      !!
+      INTEGER , INTENT(in)                         ::   kt    ! ocean time-step index
+      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   prd   ! in situ density
+      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pn2   ! Brunt-Vaisala frequency (locally ref.)
+      !!
+      INTEGER  ::   ji , jj , jk    ! dummy loop indices
+      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
+      !!----------------------------------------------------------------------
       
-      ! 0. Initialization (first time-step only)
-      !    --------------
-      
-      IF( kt == nit000 ) CALL ldf_slp_init
-      
-
-      ! 0. Local constant initialization
-      ! --------------------------------
-      
-      zeps  =  1.e-20
+      IF( kt == nit000 )   CALL ldf_slp_init      ! initialization (first time-step only)
+
+      zeps  =  1.e-20                             ! Local constant initialization
       zmg   = -1.0 / grav
       zm05g = -0.5 / grav
-
+      !
       zww(:,:,:) = 0.e0
       zwz(:,:,:) = 0.e0
-
-      ! horizontal density gradient computation
+      !                                           ! horizontal density gradient computation
       DO jk = 1, jpk
          DO jj = 1, jpjm1
@@ -139,56 +114,36 @@
          END DO
       END DO
-
-      IF( ln_zps ) THEN      ! partial steps correction at the bottom ocean level (zps_hde routine)
+      IF( ln_zps ) THEN                           ! partial steps correction at the bottom ocean level
 # if defined key_vectopt_loop  
-         jj = 1
-         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+         DO jj = 1, 1
+            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
 # else
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
 # endif
-               ! last ocean level
-               iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1
+               iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1               ! last ocean level
                ikv = MIN ( mbathy(ji,jj), mbathy(ji,jj+1) ) - 1
                zgru(ji,jj,iku) = gru(ji,jj) 
                zgrv(ji,jj,ikv) = grv(ji,jj)               
-# if ! defined key_vectopt_loop 
-            END DO
-# endif
+            END DO
          END DO
       ENDIF
-
-      ! Slopes of isopycnal surfaces just below the mixed layer
-      ! -------------------------------------------------------
       
-      CALL ldf_slp_mxl( prd, pn2 )
+      CALL ldf_slp_mxl( prd, pn2 )                ! Slopes of isopycnal surfaces just below the mixed layer
+
       
-      !-------------------synchro---------------------------------------------
-
-      !                                                ! ===============
-      DO jk = 2, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-
-         ! I.  slopes at u and v point
-         ! ===========================
-         
-         
-         ! I.1. Slopes of isopycnal surfaces
-         ! ---------------------------------
-         ! uslp = d/di( prd ) / d/dz( prd )
-         ! vslp = d/dj( prd ) / d/dz( prd )
-         
-         ! Local vertical density gradient evaluated from N^2
-         ! zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
-
+      ! I.  slopes at u and v point      | uslp = d/di( prd ) / d/dz( prd )
+      ! ===========================      | 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
-         
-         ! Slope at u and v points
+               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.
@@ -204,51 +159,23 @@
                ! uslp and vslp output in zwz and zww, resp.
                zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
-               zwz (ji,jj,jk) = ( zau / ( zbu - zeps ) * ( 1. - zalpha)   &
-                  &           + 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)
+               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) = ( zav / ( zbv - zeps ) * ( 1. - zalpha)   &
-                  &           + 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)
-            END DO
-         END DO
-         !                                             ! ===============
-      END DO                                           !   end of slab
-      !                                                ! ===============
-
-
-         ! lateral boundary conditions on zww and zwz
-         CALL lbc_lnk( zwz, 'U', -1. )
-         CALL lbc_lnk( zww, 'V', -1. )
-
-      !                                                ! ===============
-      DO jk = 2, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         
-         ! Shapiro filter applied in the horizontal direction
-         zcofu = 1. / 16.
-         zcofv = 1. / 16.
-         DO jj = 2, jpjm1, jpj-3   ! row jj=2 and =jpjm1 only
+               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)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( zwz, 'U', -1. )   ;   CALL lbc_lnk( zww, 'V', -1. )      ! lateral boundary conditions
+      !
+      zcofu = 1. / 16.                            !* horizontal Shapiro filter
+      zcofv = 1. / 16.
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1, jpj-3                        ! rows jj=2 and =jpjm1 only
             DO ji = 2, jpim1  
-               !uslop
-               uslp(ji,jj,jk) = zcofu * (       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)                       )
-               ! vslop
-               vslp(ji,jj,jk) = zcofv * (       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
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ! uslop
                uslp(ji,jj,jk) = zcofu * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
                   &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
@@ -256,5 +183,4 @@
                   &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
                   &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
-               ! vslop
                vslp(ji,jj,jk) = zcofv * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
                   &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
@@ -264,6 +190,19 @@
             END DO
          END DO
-
-         ! decrease along coastal boundaries
+         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)      &
+                  &                       +      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)      &
+                  &                       +      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.
@@ -276,29 +215,22 @@
             END DO
          END DO
-
-
-         ! II. Computation of slopes at w point
-         ! ====================================
-         
-         
-         ! II.1 Slopes of isopycnal surfaces
-         ! ---------------------------------
-         ! wslpi = mij( d/di( prd ) / d/dz( prd )
-         ! wslpj = mij( d/dj( prd ) / d/dz( prd )
-         
-         
-         ! Local vertical density gradient evaluated from N^2
-         !     zwy = d/dz(prd)= - mk ( prd ) / grav * pn2 -- at w point
+      END DO
+
+
+      ! II.  slopes at w point           | wslpi = mij( d/di( prd ) / d/dz( prd )
+      ! ===========================      | 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
-         
-         ! Slope at w point
+               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 jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ! horizontal density i-gradient at w-points
+               !                                        ! 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) )
@@ -306,5 +238,5 @@
                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
+               !                                        ! 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  ) )
@@ -312,10 +244,9 @@
                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)
-               ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
-               !                   static instability:
-               !                   kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
+               !                                        ! 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.
+               !                                        ! 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. )
@@ -326,22 +257,12 @@
             END DO
          END DO
-         !                                             ! ===============
-      END DO                                           !   end of slab
-      !                                                ! ===============
-
-
-      ! lateral boundary conditions on zwz and zww
-      CALL lbc_lnk( zwz, 'T', -1. )
-      CALL lbc_lnk( zww, 'T', -1. )
-
-      !                                                ! ===============
-      DO jk = 2, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-
-         ! Shapiro filter applied in the horizontal direction
-
-         DO jj = 2, jpjm1, jpj-3   ! row jj=2 and =jpjm1
+      END DO
+      CALL lbc_lnk( zwz, 'T', -1. )   ;    CALL lbc_lnk( zww, 'T', -1. )      ! lateral boundary conditions on zwz and zww
+      !
+      !                                           !* horizontal Shapiro filter
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1, jpj-3                        ! rows jj=2 and =jpjm1
             DO ji = 2, jpim1
-               zcofw = tmask(ji,jj,jk)/16.
+               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)     &
@@ -356,9 +277,8 @@
                   &                + 4.*  zww(ji  ,jj  ,jk)                        ) * zcofw
             END DO
-         END DO
-         
-         DO jj = 3, jpj-2
+         END DO  
+         DO jj = 3, jpj-2                               ! other rows
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zcofw = tmask(ji,jj,jk)/16.
+               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)     &
@@ -374,6 +294,5 @@
             END DO
          END DO
-         
-         ! decrease the slope along the boundaries
+         !                                        !* decrease along coastal boundaries
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -384,39 +303,34 @@
             END DO
          END DO
+      END DO
          
-         
-         ! III. Specific grid points
-         ! -------------------------
-
-         IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN
-            !                                        ! =======================
-            ! Horizontal diffusion in                !  ORCA_R4 configuration 
-            ! specific area                          ! =======================
-            !
-            !                                             ! Gibraltar Strait
-            ij0 =  50   ;   ij1 =  53
-            ii0 =  69   ;   ii1 =  71   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk ) = 0.e0
-            ij0 =  51   ;   ij1 =  53
-            ii0 =  68   ;   ii1 =  71   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk ) = 0.e0
-            ii0 =  69   ;   ii1 =  71   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk ) = 0.e0
-            ii0 =  69   ;   ii1 =  71   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk ) = 0.e0
-
-            !                                             ! Mediterrannean Sea
-            ij0 =  49   ;   ij1 =  56
-            ii0 =  71   ;   ii1 =  90   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk ) = 0.e0
-            ij0 =  50   ;   ij1 =  56
-            ii0 =  70   ;   ii1 =  90   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk ) = 0.e0
-            ii0 =  71   ;   ii1 =  90   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk ) = 0.e0
-            ii0 =  71   ;   ii1 =  90   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk ) = 0.e0
-         ENDIF
-         !                                             ! ===============
-      END DO                                           !   end of slab
-      !                                                ! ===============       
+
+      ! 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.e0
+         ij0 =  51   ;   ij1 =  53
+         ii0 =  68   ;   ii1 =  71   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
+         ii0 =  69   ;   ii1 =  71   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
+         ii0 =  69   ;   ii1 =  71   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
+         !
+         !                                                    ! Mediterrannean Sea
+         ij0 =  49   ;   ij1 =  56
+         ii0 =  71   ;   ii1 =  90   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
+         ij0 =  50   ;   ij1 =  56
+         ii0 =  70   ;   ii1 =  90   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
+         ii0 =  71   ;   ii1 =  90   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
+         ii0 =  71   ;   ii1 =  90   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
+      ENDIF
 
       
-      ! III Lateral boundary conditions on all slopes (uslp , vslp, 
-      ! -------------------------------                wslpi, wslpj )
+      ! 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
@@ -424,5 +338,5 @@
          CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
       ENDIF
-
+      !
    END SUBROUTINE ldf_slp
 
@@ -431,7 +345,7 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_slp_mxl  ***
-      !! ** Purpose :
-      !!     Compute the slopes of iso-neutral surface (slope of isopycnal
-      !!   surfaces referenced locally) just above the mixed layer.
+      !!
+      !! ** Purpose :   Compute the slopes of iso-neutral surface just below 
+      !!              the mixed layer.
       !!
       !! ** Method :
@@ -445,42 +359,26 @@
       !!   of 10cm/s)
       !!
-      !! ** Action :
-      !!      Compute uslp, wslpi, and vslp, wslpj, the i- and  j-slopes
+      !! ** Action :   Compute uslp, wslpi, and vslp, wslpj, the i- and  j-slopes
       !!   of now neutral surfaces at u-, w- and v- w-points, resp.
-      !!
-      !! History :
-      !!   8.1  !  99-10  (A. Jouzeau)  Original code
-      !!   8.5  !  99-10  (G. Madec)  Free form, F90
-      !!----------------------------------------------------------------------
-      !! * Modules used
-      USE oce           , zgru  => ua,  &  ! ua, va used as workspace and set to hor. 
-                          zgrv  => va      ! density gradient in ldf_slp
-
-      !! * Arguments
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   &
-         prd,                           &  ! in situ density
-         pn2                              ! Brunt-Vaisala frequency (locally ref.)
-
-      !! * Local declarations
-      INTEGER  ::   ji, jj, jk             ! dummy loop indices
-      INTEGER  ::   ik, ikm1               ! temporary integers
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         zwy                               ! temporary workspace
-      REAL(wp) ::   &
-         zeps, zmg, zm05g,              &  ! temporary scalars
-         zcoef1, zcoef2,                &  !    "         "
-         zau, zbu, zav, zbv,            &  !    "         "
-         zai, zbi, zaj, zbj                !    "         "
-      !!----------------------------------------------------------------------
-
-
-      ! 0. Local constant initialization
-      ! --------------------------------
-
-      zeps  =  1.e-20
+      !!----------------------------------------------------------------------
+      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
-
-
+      !
       uslpml (1,:) = 0.e0      ;      uslpml (jpi,:) = 0.e0
       vslpml (1,:) = 0.e0      ;      vslpml (jpi,:) = 0.e0
@@ -488,25 +386,18 @@
       wslpjml(1,:) = 0.e0      ;      wslpjml(jpi,:) = 0.e0
 
-      ! surface mixed layer mask
-
-      ! mask for mixed layer
-      DO jk = 1, jpk
+      !                                  ! surface mixed layer mask
+      DO jk = 1, jpk                          ! =1 inside the mixed layer, =0 otherwise
 # if defined key_vectopt_loop
-         jj = 1
-         DO ji = 1, jpij   ! vector opt. (forced unrolling)
+         DO jj = 1, 1
+            DO ji = 1, jpij   ! vector opt. (forced unrolling)
 # else
          DO jj = 1, jpj
             DO ji = 1, jpi
 # endif
-               ! mixed layer interior (mask = 1) and exterior (mask = 0)
                ik = nmln(ji,jj) - 1
-               IF( jk <= ik ) THEN
-                  omlmask(ji,jj,jk) = 1.e0
-               ELSE
-                  omlmask(ji,jj,jk) = 0.e0
+               IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1.e0
+               ELSE                  ;   omlmask(ji,jj,jk) = 0.e0
                ENDIF
-# if ! defined key_vectopt_loop
-            END DO
-# endif
+            END DO
          END DO
       END DO
@@ -515,37 +406,33 @@
       ! Slopes of isopycnal surfaces just before bottom of mixed layer
       ! --------------------------------------------------------------
-      ! uslpml = d/di( prd ) / d/dz( prd )
-      ! vslpml = d/dj( prd ) / d/dz( prd )
-
-      ! Local vertical density gradient evaluated from N^2
-      ! zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
-
+      ! The slope are computed as in the 3D case.
+      ! A key point here is the definition of the mixed layer at u- and v-points.
+      ! It is assumed to be the maximum of the two neighbouring T-point mixed layer depth.
+      ! Otherwise, a n2 value inside the mixed layer can be involved in the computation
+      ! of the slope, resulting in a too steep diagnosed slope and thus a spurious eddy
+      ! induce velocity field near the base of the mixed layer.
       !-----------------------------------------------------------------------
-      zwy(:,jpj) = 0.e0
+      !
+      zwy(:,jpj) = 0.e0            !* vertical density gradient for u-slope (from N^2)
       zwy(jpi,:) = 0.e0
 # if defined key_vectopt_loop
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+      DO jj = 1, 1
+         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
 # else
       DO jj = 1, jpjm1
          DO ji = 1, jpim1
 # endif
-            ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) )
-            ! if ik = jpk take jpkm1 values
-            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. )
-# if ! defined key_vectopt_loop
-         END DO
-# endif
-      END DO
-      CALL lbc_lnk( zwy, 'U', 1. )      ! lateral boundary conditions (NO sign change)
-
-
-      ! Slope at u points
+            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
-      jj = 1
-      DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
+      DO jj = 1, 1
+         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
 # else
       DO jj = 2, jpjm1
@@ -554,5 +441,5 @@
             ! horizontal and vertical density gradient at u-points
             ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) )
-            ik = MIN( ik,jpkm1 )
+            ik = MIN( ik, jpkm1 )
             zau = 1./ e1u(ji,jj) * zgru(ji,jj,ik)
             zbu = 0.5*( zwy(ji,jj) + zwy(ji+1,jj) )
@@ -562,19 +449,12 @@
             ! uslpml
             uslpml (ji,jj) = zau / ( zbu - zeps ) * umask (ji,jj,ik)
-# if ! defined key_vectopt_loop
-         END DO
-# endif
-      END DO
-
-      ! lateral boundary conditions on uslpml
-      CALL lbc_lnk( uslpml, 'U', -1. )
-
-      ! Local vertical density gradient evaluated from N^2
-      !     zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
-      zwy ( :, jpj) = 0.e0
-      zwy ( jpi, :) = 0.e0
+         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
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+      DO jj = 1, 1
+         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
 # else
       DO jj = 1, jpjm1
@@ -582,18 +462,15 @@
 # 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. )
-# if ! defined key_vectopt_loop
-         END DO
-# endif
-      END DO
-      CALL lbc_lnk( zwy, 'V', 1. )      ! lateral boundary conditions No change in sign
-
-      ! Slope at v points
+            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
-      jj = 1
-      DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
+      DO jj = 1, 1
+         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
 # else
       DO jj = 2, jpjm1
@@ -610,46 +487,36 @@
             ! vslpml
             vslpml (ji,jj) = zav / ( zbv - zeps ) * vmask (ji,jj,ik)
-# if ! defined key_vectopt_loop
-         END DO
-# endif
-      END DO
-
-      ! lateral boundary conditions on vslpml
-      CALL lbc_lnk( vslpml, 'V', -1. )
-
-      ! wslpiml = mij( d/di( prd ) / d/dz( prd )
-      ! wslpjml = mij( d/dj( prd ) / d/dz( prd )
-
-
-      ! Local vertical density gradient evaluated from N^2
-      ! zwy = d/dz(prd)= - mk ( prd ) / grav * pn2 -- at w point
+         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
-      jj = 1
-      DO ji = 1, jpij   ! vector opt. (forced unrolling)
+      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 )
+            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. )
-# if ! defined key_vectopt_loop
-         END DO
-# endif
-      END DO
-
-      ! Slope at w point
+         END DO
+      END DO
+
+      !                            !* Slopes at w points
 # if defined key_vectopt_loop
-      jj = 1
-      DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
+      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)
+            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  )   &
@@ -672,13 +539,8 @@
             wslpiml (ji,jj) = zai / ( zbi - zeps) * tmask (ji,jj,ik)
             wslpjml (ji,jj) = zaj / ( zbj - zeps) * tmask (ji,jj,ik)
-# if ! defined key_vectopt_loop
-         END DO
-# endif
-      END DO
-
-      ! lateral boundary conditions on wslpiml and wslpjml
-      CALL lbc_lnk( wslpiml, 'W', -1. )
-      CALL lbc_lnk( wslpjml, 'W', -1. )
-
+         END DO
+      END DO
+      CALL lbc_lnk( wslpiml, 'W', -1. )   ;   CALL lbc_lnk( wslpjml, 'W', -1. )      ! lateral boundary conditions
+      !
    END SUBROUTINE ldf_slp_mxl
 
@@ -692,16 +554,9 @@
       !! ** Method  :   read the nammbf namelist and check the parameter 
       !!      values called by tra_dmp at the first timestep (nit000)
-      !!
-      !! History :
-      !!  8.5  ! 02-06 (G. Madec) original code
-      !!----------------------------------------------------------------------
-      !! * local declarations
+      !!----------------------------------------------------------------------
       INTEGER ::   ji, jj, jk   ! dummy loop indices
       !!----------------------------------------------------------------------
       
-      
-      ! Parameter control and print
-      ! ---------------------------
-      IF(lwp) THEN
+      IF(lwp) THEN    
          WRITE(numout,*)
          WRITE(numout,*) 'ldf_slp : direction of lateral mixing'
@@ -744,10 +599,9 @@
             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. )
       ENDIF
-
+      !
    END SUBROUTINE ldf_slp_init
 
@@ -760,5 +614,5 @@
    SUBROUTINE ldf_slp( kt, prd, pn2 )        ! Dummy routine
       INTEGER, INTENT(in) :: kt 
-      REAL,DIMENSION(:,:,:), INTENT(in) :: prd, pn2
+      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
