Index: trunk/NEMO/OPA_SRC/TRA/traadv_muscl.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traadv_muscl.F90	(revision 96)
+++ trunk/NEMO/OPA_SRC/TRA/traadv_muscl.F90	(revision 97)
@@ -79,5 +79,4 @@
       REAL(wp) ::   zu, zv, zw, zeu, zev, zew, zbtr, zstep, zta, zsa
       REAL(wp) ::   z0u, z0v, z0w
-      REAL(wp) ::   z1u, z1v, z1w
       REAL(wp) ::   zzt1, zzt2, zalpha
       REAL(wp) ::   zzs1, zzs2
@@ -138,30 +137,13 @@
       DO jk = 1, jpkm1
          DO jj = 2, jpj
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               z0u = zt1(ji,jj,jk) * zt1(ji-1,jj,jk)
-               IF( z0u > 0. ) THEN
-                  ztp1(ji,jj,jk) = 0.5 * ( zt1(ji,jj,jk)+zt1(ji-1,jj,jk) )
-               ELSE
-                  ztp1(ji,jj,jk) = 0.e0
-               ENDIF
-               z1u = zs1(ji,jj,jk) * zs1(ji-1,jj,jk)
-               IF( z1u > 0. ) THEN
-                  zsp1(ji,jj,jk) = 0.5 * ( zs1(ji,jj,jk)+zs1(ji-1,jj,jk) )
-               ELSE
-                  zsp1(ji,jj,jk) = 0.e0
-               ENDIF
-
-               z0v = zt2(ji,jj,jk) * zt2(ji,jj-1,jk)
-               IF( z0v > 0. ) THEN
-                  ztp2(ji,jj,jk) = 0.5 * ( zt2(ji,jj,jk)+zt2(ji,jj-1,jk) )
-               ELSE
-                  ztp2(ji,jj,jk) = 0.e0
-               ENDIF
-               z1v = zs2(ji,jj,jk) * zs2(ji,jj-1,jk)
-               IF( z1v > 0. ) THEN
-                  zsp2(ji,jj,jk) = 0.5 * ( zs2(ji,jj,jk)+zs2(ji,jj-1,jk) )
-               ELSE
-                  zsp2(ji,jj,jk) = 0.e0
-               ENDIF
+            DO ji = fs_2, jpi   ! vector opt.
+               ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji-1,jj  ,jk) )   &
+                  &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji-1,jj  ,jk) ) )
+               zsp1(ji,jj,jk) =                    ( zs1(ji,jj,jk) + zs1(ji-1,jj  ,jk) )   &
+                  &           * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji-1,jj  ,jk) ) )
+               ztp2(ji,jj,jk) =                    ( zt2(ji,jj,jk) + zt2(ji  ,jj-1,jk) )   &
+                  &           * ( 0.25 + SIGN( 0.25, zt2(ji,jj,jk) * zt2(ji  ,jj-1,jk) ) )
+               zsp2(ji,jj,jk) =                    ( zs2(ji,jj,jk) + zs2(ji  ,jj-1,jk) )   &
+                  &           * ( 0.25 + SIGN( 0.25, zs2(ji,jj,jk) * zs2(ji  ,jj-1,jk) ) )
             END DO
          END DO
@@ -171,8 +153,8 @@
       zsp1(:,:,jpk) = 0.e0    ;    zsp2(:,:,jpk) = 0.e0
 
-! Slopes limitation
+      ! Slopes limitation
       DO jk = 1, jpkm1
          DO jj = 2, jpj
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO ji = fs_2, jpi   ! vector opt.
                ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
                   &           * MIN(    ABS( ztp1(ji  ,jj,jk) ),   &
@@ -285,4 +267,11 @@
       END DO        
 
+      IF(l_ctl) THEN
+         zta = SUM( ta(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
+         zsa = SUM( sa(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
+         WRITE(numout,*) ' had  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' muscl'
+         t_ctl = zta   ;   s_ctl = zsa
+      ENDIF
+
 #if defined key_diaptr 
       ! "zonal" mean advective heat and salt transport 
@@ -324,16 +313,8 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               z0w = zt1(ji,jj,jk) * zt1(ji,jj,jk+1) 
-               IF( z0w > 0. ) THEN
-                  ztp1(ji,jj,jk) = 0.5 * ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )
-               ELSE
-                  ztp1(ji,jj,jk) = 0.e0
-               ENDIF
-               z1w = zs1(ji,jj,jk) * zs1(ji,jj,jk+1) 
-               IF( z1w > 0. ) THEN
-                  zsp1(ji,jj,jk) = 0.5 * ( zs1(ji,jj,jk) + zs1(ji,jj,jk+1) )
-               ELSE
-                  zsp1(ji,jj,jk) = 0.e0
-               ENDIF
+               ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )   &
+                  &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) )
+               zsp1(ji,jj,jk) =                    ( zs1(ji,jj,jk) + zs1(ji,jj,jk+1) )   &
+                  &           * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji,jj,jk+1) ) )
             END DO
          END DO
@@ -357,6 +338,6 @@
       END DO
       ! surface values
-      ztp1(:,:,1) = 0. 
-      zsp1(:,:,1) = 0.
+      ztp1(:,:,1) = 0.e0
+      zsp1(:,:,1) = 0.e0
 
       ! vertical advective flux
Index: trunk/NEMO/OPA_SRC/TRA/traadv_muscl2.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traadv_muscl2.F90	(revision 96)
+++ trunk/NEMO/OPA_SRC/TRA/traadv_muscl2.F90	(revision 97)
@@ -80,5 +80,4 @@
       REAL(wp) ::   zu, zv, zw, zeu, zev, zew, zbtr, zstep, zta, zsa
       REAL(wp) ::   z0u, z0v, z0w
-      REAL(wp) ::   z1u, z1v, z1w
       REAL(wp) ::   zzt1, zzt2, zalpha
       REAL(wp) ::   zzs1, zzs2
@@ -139,30 +138,13 @@
       DO jk = 1, jpkm1
          DO jj = 2, jpj
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               z0u = zt1(ji,jj,jk) * zt1(ji-1,jj,jk)
-               IF( z0u > 0. ) THEN
-                  ztp1(ji,jj,jk) = 0.5 * ( zt1(ji,jj,jk)+zt1(ji-1,jj,jk) )
-               ELSE
-                  ztp1(ji,jj,jk) = 0.e0
-               ENDIF
-               z1u = zs1(ji,jj,jk) * zs1(ji-1,jj,jk)
-               IF( z1u > 0. ) THEN
-                  zsp1(ji,jj,jk) = 0.5 * ( zs1(ji,jj,jk)+zs1(ji-1,jj,jk) )
-               ELSE
-                  zsp1(ji,jj,jk) = 0.e0
-               ENDIF
-
-               z0v = zt2(ji,jj,jk) * zt2(ji,jj-1,jk)
-               IF( z0v > 0. ) THEN
-                  ztp2(ji,jj,jk) = 0.5 * ( zt2(ji,jj,jk)+zt2(ji,jj-1,jk) )
-               ELSE
-                  ztp2(ji,jj,jk) = 0.e0
-               ENDIF
-               z1v = zs2(ji,jj,jk) * zs2(ji,jj-1,jk)
-               IF( z1v > 0. ) THEN
-                  zsp2(ji,jj,jk) = 0.5 * ( zs2(ji,jj,jk)+zs2(ji,jj-1,jk) )
-               ELSE
-                  zsp2(ji,jj,jk) = 0.e0
-               ENDIF
+            DO ji = fs_2, jpi   ! vector opt.
+               ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji-1,jj  ,jk) )   &
+                  &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji-1,jj  ,jk) ) )
+               zsp1(ji,jj,jk) =                    ( zs1(ji,jj,jk) + zs1(ji-1,jj  ,jk) )   &
+                  &           * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji-1,jj  ,jk) ) )
+               ztp2(ji,jj,jk) =                    ( zt2(ji,jj,jk) + zt2(ji  ,jj-1,jk) )   &
+                  &           * ( 0.25 + SIGN( 0.25, zt2(ji,jj,jk) * zt2(ji  ,jj-1,jk) ) )
+               zsp2(ji,jj,jk) =                    ( zs2(ji,jj,jk) + zs2(ji  ,jj-1,jk) )   &
+                  &           * ( 0.25 + SIGN( 0.25, zs2(ji,jj,jk) * zs2(ji  ,jj-1,jk) ) )
             END DO
          END DO
@@ -172,8 +154,8 @@
       zsp1(:,:,jpk) = 0.e0    ;    zsp2(:,:,jpk) = 0.e0
 
-! Slopes limitation
+      ! Slopes limitation
       DO jk = 1, jpkm1
          DO jj = 2, jpj
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO ji = fs_2, jpi   ! vector opt.
                ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
                   &           * MIN(    ABS( ztp1(ji  ,jj,jk) ),   &
@@ -346,4 +328,11 @@
       END DO        
 
+      IF(l_ctl) THEN
+         zta = SUM( ta(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
+         zsa = SUM( sa(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
+         WRITE(numout,*) ' had  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' muscl'
+         t_ctl = zta   ;   s_ctl = zsa
+      ENDIF
+
 #if defined key_diaptr
       ! "zonal" mean advective heat and salt transport
@@ -384,16 +373,8 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               z0w = zt1(ji,jj,jk) * zt1(ji,jj,jk+1) 
-               IF( z0w > 0. ) THEN
-                  ztp1(ji,jj,jk) = 0.5 * ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )
-               ELSE
-                  ztp1(ji,jj,jk) = 0.e0
-               ENDIF
-               z1w = zs1(ji,jj,jk) * zs1(ji,jj,jk+1) 
-               IF( z1w > 0. ) THEN
-                  zsp1(ji,jj,jk) = 0.5 * ( zs1(ji,jj,jk) + zs1(ji,jj,jk+1) )
-               ELSE
-                  zsp1(ji,jj,jk) = 0.e0
-               ENDIF
+               ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )   &
+                  &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) )
+               zsp1(ji,jj,jk) =                    ( zs1(ji,jj,jk) + zs1(ji,jj,jk+1) )   &
+                  &           * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji,jj,jk+1) ) )
             END DO
          END DO
@@ -417,6 +398,6 @@
       END DO
       ! surface values
-      ztp1(:,:,1) = 0. 
-      zsp1(:,:,1) = 0.
+      ztp1(:,:,1) = 0.e0
+      zsp1(:,:,1) = 0.e0
 
       ! vertical advective flux
