Index: trunk/NEMO/OFF_SRC/LDF/ldfslp.F90
===================================================================
--- trunk/NEMO/OFF_SRC/LDF/ldfslp.F90	(revision 343)
+++ trunk/NEMO/OFF_SRC/LDF/ldfslp.F90	(revision 497)
@@ -20,4 +20,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -46,7 +47,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL  (2005)
-   !!   $Header$
-   !!   This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
+   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
    !!----------------------------------------------------------------------
 
@@ -69,9 +68,9 @@
       !!      of 10cm/s)
       !!        A horizontal shapiro filter is applied to the slopes
-      !!        'key_s_coord' defined: add to the previously computed slopes
+      !!        ln_sco=T, s-coordinate, add to the previously computed slopes
       !!      the slope of the model level surface.
       !!        macro-tasked on horizontal slab (jk-loop)  (2, jpk-1)
       !!      [slopes already set to zero at level 1, and to zero or the ocean
-      !!      bottom slope ('key_s_coord' defined) at level jpk in inildf]
+      !!      bottom slope (ln_sco=T) at level jpk in inildf]
       !!
       !! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and  j-slopes 
@@ -83,4 +82,5 @@
       !!   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
@@ -97,14 +97,13 @@
       !! * Local declarations
       INTEGER  ::   ji, jj, jk              ! dummy loop indices
-      INTEGER  ::   ii0, ii1, ij0, ij1      ! temporary integer
-#if defined key_partial_steps
-      INTEGER  ::   iku, ikv  ! temporary integers
-#endif
+      INTEGER  ::   ii0, ii1, ij0, ij1,  &  ! temporary integer
+         &          iku, ikv                !    "          "
       REAL(wp) ::   &
-         zeps, zmg, zm05g, zcoef1, zcoef2,   &  ! temporary scalars
-         zau, zbu, zav, zbv,                 &
-         zai, zbi, zaj, zbj,                &
-         zcofu, zcofv, zcofw,                &
-         z1u, z1v, z1wu, z1wv,               &
+         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
@@ -138,23 +137,22 @@
       END DO
 
-#if defined key_partial_steps
-      ! partial steps correction at the bottom ocean level (zps_hde routine)
-# if defined key_vectopt_loop   &&   ! defined key_autotasking
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+      IF( ln_zps ) THEN      ! partial steps correction at the bottom ocean level (zps_hde routine)
+# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
+         jj = 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
-            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   ||   defined key_autotasking
-         END DO
-# endif
-      END DO
-#endif
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+# endif
+               ! last ocean level
+               iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1
+               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   ||   defined key_mpp_omp
+            END DO
+# endif
+         END DO
+      ENDIF
 
       ! Slopes of isopycnal surfaces just below the mixed layer
@@ -203,30 +201,13 @@
                ! uslp and vslp output in zwz and zww, resp.
                zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
-#if defined key_s_coord 
                zwz (ji,jj,jk) = ( zau / ( zbu - zeps ) * ( 1. - zalpha)   &
-                  &        + zalpha * uslpml(ji,jj)   &
-                  &        * ( fsdepu(ji,jj,jk) - .5*fse3u(ji,jj,1) )   &
-                  &        / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5. ) )   &
-                  &        * umask(ji,jj,jk)
+                  &           + 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)   &
-                  &        * ( fsdepv(ji,jj,jk) - .5*fse3v(ji,jj,1) )   &
-                  &         / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) )   &
-                  &        * vmask(ji,jj,jk)
-#else
-               ! z-coord and partial steps slope computed in the same way
-               zwz (ji,jj,jk) = ( zau / ( zbu - zeps ) * ( 1. - zalpha)    &
-                  &        + zalpha * uslpml(ji,jj)    &
-                  &        * ( fsdept(ji,jj,jk) - .5*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)    &
-                  &        * ( fsdept(ji,jj,jk) - .5*fse3v(ji,jj,1))    &
-                  &        / MAX(hmlpt(ji,jj),hmlpt(ji,jj+1),5.) )    &
-                  &        * vmask (ji,jj,jk)
-#endif
+                  &           + 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
@@ -292,24 +273,4 @@
             END DO
          END DO
-
-
-         IF( lk_sco ) THEN
-            ! Add the slope of level surfaces
-            ! -----------------------------------
-            ! 'key_s_coord' defined but not 'key_traldfiso' the computation is done
-            ! in inildf, ldfslp never called
-            ! 'key_s_coord' and 'key_traldfiso' defined, the slope of level surfaces
-            ! is added to the slope of isopycnal surfaces.
-            ! c a u t i o n : minus sign as fsdep has positive value 
-         
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  uslp(ji,jj,jk) = uslp(ji,jj,jk) - 1. / e1u(ji,jj)   &
-                     &           * ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) )
-                  vslp(ji,jj,jk) = vslp(ji,jj,jk) - 1. / e2v(ji,jj)   &
-                     &           * ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) )
-               END DO
-            END DO
-         ENDIF
 
 
@@ -354,13 +315,10 @@
                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))
-               zwz(ji,jj,jk) = ( zai / ( zbi - zeps) * ( 1. - zalpha )   &
-                  &            + zalpha * wslpiml(ji,jj)   &
-                  &            * fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj),10. ) )   &
-                  &            * tmask (ji,jj,jk)
-               zww(ji,jj,jk) = ( zaj / ( zbj - zeps) * ( 1. - zalpha )   &
-                  &            + zalpha * wslpjml(ji,jj)   &
-                  &            * fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj),10. ) )   &
-                  &            * tmask (ji,jj,jk)
+               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
@@ -424,22 +382,4 @@
          END DO
          
-         IF( lk_sco ) THEN
-         
-            ! Slope of level surfaces
-            ! -----------------------
-            ! 'key_s_coord' defined but not 'key_traldfiso' the computation is done
-            ! in inildf, ldfslp never called
-            ! 'key_s_coord' and 'key_traldfiso' defined, the slope of level surfaces
-            ! is added to the slope of isopycnal surfaces.
-         
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  wslpi(ji,jj,jk) = wslpi(ji,jj,jk) - 1. / e1t(ji,jj)   &
-                     &                                   * ( fsdepuw(ji+1,jj,jk) - fsdepuw(ji,jj,jk) )
-                  wslpj(ji,jj,jk) = wslpj(ji,jj,jk) - 1. / e2t(ji,jj)   &
-                     &                                   * ( fsdepvw(ji,jj+1,jk) - fsdepvw(ji,jj,jk) )
-               END DO
-            END DO
-         ENDIF
          
          ! III. Specific grid points
@@ -474,8 +414,11 @@
       ! III Lateral boundary conditions on all slopes (uslp , vslp, 
       ! -------------------------------                wslpi, wslpj )
-      CALL lbc_lnk( uslp , 'U', -1. )
-      CALL lbc_lnk( vslp , 'V', -1. )
-      CALL lbc_lnk( wslpi, 'W', -1. )
-      CALL lbc_lnk( wslpj, 'W', -1. )
+      CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
+      CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
+
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp  - u : ', tab3d_2=vslp,  clinfo2=' v : ', kdim=jpk)
+         CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
+      ENDIF
 
    END SUBROUTINE ldf_slp
@@ -546,5 +489,5 @@
       ! mask for mixed layer
       DO jk = 1, jpk
-# if defined key_vectopt_loop   &&   ! defined key_autotasking
+# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
          jj = 1
          DO ji = 1, jpij   ! vector opt. (forced unrolling)
@@ -560,5 +503,5 @@
                   omlmask(ji,jj,jk) = 0.e0
                ENDIF
-# if ! defined key_vectopt_loop   ||   defined key_autotasking
+# if ! defined key_vectopt_loop   ||   defined key_mpp_omp
             END DO
 # endif
@@ -578,5 +521,5 @@
       zwy(:,jpj) = 0.e0
       zwy(jpi,:) = 0.e0
-# if defined key_vectopt_loop   &&   ! defined key_autotasking
+# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
       jj = 1
       DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
@@ -591,5 +534,5 @@
                &             * ( 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   ||   defined key_autotasking
+# if ! defined key_vectopt_loop   ||   defined key_mpp_omp
          END DO
 # endif
@@ -599,5 +542,5 @@
 
       ! Slope at u points
-# if defined key_vectopt_loop   &&   ! defined key_autotasking
+# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
       jj = 1
       DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
@@ -616,5 +559,5 @@
             ! uslpml
             uslpml (ji,jj) = zau / ( zbu - zeps ) * umask (ji,jj,ik)
-# if ! defined key_vectopt_loop   ||   defined key_autotasking
+# if ! defined key_vectopt_loop   ||   defined key_mpp_omp
          END DO
 # endif
@@ -628,5 +571,5 @@
       zwy ( :, jpj) = 0.e0
       zwy ( jpi, :) = 0.e0
-# if defined key_vectopt_loop   &&   ! defined key_autotasking
+# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
       jj = 1
       DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
@@ -640,5 +583,5 @@
                &             * ( 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   ||   defined key_autotasking
+# if ! defined key_vectopt_loop   ||   defined key_mpp_omp
          END DO
 # endif
@@ -649,5 +592,5 @@
 
       ! Slope at v points
-# if defined key_vectopt_loop   &&   ! defined key_autotasking
+# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
       jj = 1
       DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
@@ -666,5 +609,5 @@
             ! vslpml
             vslpml (ji,jj) = zav / ( zbv - zeps ) * vmask (ji,jj,ik)
-# if ! defined key_vectopt_loop   ||   defined key_autotasking
+# if ! defined key_vectopt_loop   ||   defined key_mpp_omp
          END DO
 # endif
@@ -680,5 +623,5 @@
       ! Local vertical density gradient evaluated from N^2
       ! zwy = d/dz(prd)= - mk ( prd ) / grav * pn2 -- at w point
-# if defined key_vectopt_loop   &&   ! defined key_autotasking
+# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
       jj = 1
       DO ji = 1, jpij   ! vector opt. (forced unrolling)
@@ -692,5 +635,5 @@
             zwy (ji,jj) = zm05g * pn2 (ji,jj,ik) *     &
                &             ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. )
-# if ! defined key_vectopt_loop   ||   defined key_autotasking
+# if ! defined key_vectopt_loop   ||   defined key_mpp_omp
          END DO
 # endif
@@ -698,5 +641,5 @@
 
       ! Slope at w point
-# if defined key_vectopt_loop   &&   ! defined key_autotasking
+# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
       jj = 1
       DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
@@ -728,5 +671,5 @@
             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   ||   defined key_autotasking
+# if ! defined key_vectopt_loop   ||   defined key_mpp_omp
          END DO
 # endif
@@ -780,4 +723,7 @@
 
       IF( ln_traldf_hor ) THEN
+         IF(lwp) THEN
+            WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
+         ENDIF
 
          ! geopotential diffusion in s-coordinates on tracers and/or momentum
@@ -790,12 +736,8 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  uslp (ji,jj,jk) = -1. / e1u(ji,jj) * umask(ji,jj,jk)   &
-                     &                               * ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) )
-                  vslp (ji,jj,jk) = -1. / e2v(ji,jj) * vmask(ji,jj,jk)   &
-                     &                               * ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) )
-                  wslpi(ji,jj,jk) = -1. / e1t(ji,jj) * tmask(ji,jj,jk)   &
-                     &                               * ( fsdepuw(ji+1,jj,jk) - fsdepuw(ji,jj,jk) )
-                  wslpj(ji,jj,jk) = -1. / e2t(ji,jj) * tmask(ji,jj,jk)   &
-                     &                               * ( fsdepvw(ji,jj+1,jk) - fsdepvw(ji,jj,jk) )
+                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)
+                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)
+                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
+                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
                END DO
             END DO
@@ -803,8 +745,6 @@
 
          ! 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. )
+         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
 
