Index: /trunk/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 2388)
+++ /trunk/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 2389)
@@ -9,4 +9,5 @@
    !!   NEMO     0.5  ! 2002-10  (G. Madec)  Free form, F90
    !!            1.0  ! 2005-10  (A. Beckmann)  correction for s-coordinates
+   !!            3.2  ! 2010-11  (F. Dupond, G. Madec)  bug correction in slopes just below the ML
    !!----------------------------------------------------------------------
 #if   defined key_ldfslp   ||   defined key_esopa
@@ -14,6 +15,6 @@
    !!   'key_ldfslp'                      Rotation of lateral mixing tensor
    !!----------------------------------------------------------------------
-   !!   ldf_slp      : compute the slopes of neutral surface
-   !!   ldf_slp_mxl  : compute the slopes of iso-neutral surface
+   !!   ldf_slp      : compute the slopes of iso-neutral surface
+   !!   ldf_slp_mxl  : compute the slopes of iso-neutral surface just below the Mixed Layer
    !!   ldf_slp_init : initialization of the slopes computation
    !!----------------------------------------------------------------------
@@ -130,19 +131,19 @@
       ENDIF
       
-      CALL ldf_slp_mxl( prd, pn2 )                ! Slopes of isopycnal surfaces just below the mixed layer
-
+
+      !                                           !* Local vertical density gradient evaluated from N^2
+      DO jk = 2, jpkm1                            !  zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
+         zwy(:,:,jk) = zmg * ( prd(:,:,jk) + 1. ) * (    pn2  (:,:,jk) + pn2  (:,:,jk+1)     )   &
+            &                                     / MAX( tmask(:,:,jk) + tmask(:,:,jk+1), 1. )
+      END DO
+      zwy(:,:,1) = 0.e0                           !  surface value set to zero
+
+      
+      CALL ldf_slp_mxl( prd, pn2 )                !* Slopes of isopycnal surfaces just below the mixed layer
+      
       
       ! 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
-      END DO
       DO jk = 2, jpkm1                            !* Slopes at u and v points
          DO jj = 2, jpjm1
@@ -349,19 +350,13 @@
       !!              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
+      !! ** Method  :   caution: use zgru, zgrv and zwy computed in ldf_slp
+      !!
+      !! ** Action  :   uslpml, wslpiml :  i- &  j-slopes of neutral surfaces
+      !!                vslpml, wslpjml    just below the mixed layer 
+      !!                omlmask         :  mixed layer mask
+      !!----------------------------------------------------------------------
+      USE oce , zgru  => ua   ! zgru, zgrv (ua, va used as workspace)  
+      USE oce , zgrv  => va   ! set to i- & j-density gradient in ldf_slp
+      USE oce , zwy   => ta   ! set to vertical density gradient at T-point in ldfslp
       !! 
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   prd   ! in situ density
@@ -374,5 +369,5 @@
       REAL(wp) ::   zau, zbu, zai, zbi   !    -         -
       REAL(wp) ::   zav, zbv, zaj, zbj   !    -         -
-      REAL(wp), DIMENSION(jpi,jpj) ::   zwy   ! 2D workspace
+      REAL(wp) ::        zbw             !    -         -
       !!----------------------------------------------------------------------
 
@@ -414,21 +409,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
@@ -443,5 +421,5 @@
             ik = MIN( ik, jpkm1 )
             zau = 1./ e1u(ji,jj) * zgru(ji,jj,ik)
-            zbu = 0.5*( zwy(ji,jj) + zwy(ji+1,jj) )
+            zbu = 0.5*( zwy(ji,jj,ik) + zwy(ji+1,jj,ik) )
             ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
             !                         kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
@@ -452,20 +430,4 @@
       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
@@ -481,5 +443,5 @@
             ik = MIN( ik,jpkm1 )
             zav = 1./ e2v(ji,jj) * zgrv(ji,jj,ik)
-            zbv = 0.5*( zwy(ji,jj) + zwy(ji,jj+1) )
+            zbv = 0.5*( zwy(ji,jj,ik) + zwy(ji,jj+1,ik) )
             ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
             !                         kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
@@ -490,21 +452,4 @@
       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
@@ -531,9 +476,11 @@
             zaj = zcoef2 * (  zgrv(ji,jj  ,ik  ) + zgrv(ji,jj  ,ikm1)   &
                &            + zgrv(ji,jj-1,ikm1) + zgrv(ji,jj-1,ik  ) ) * tmask (ji,jj,ik)
+            ! vertical density gradient at w-point (from N^2)
+            zbw = zm05g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. )
             ! 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) )
+            zbi = MIN ( zbw,- 100.*ABS(zai), -7.e+3/fse3w(ji,jj,ik)*ABS(zai) )
+            zbj = MIN ( zbw, -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)
