Index: branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl.F90	(revision 740)
+++ branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl.F90	(revision 786)
@@ -6,5 +6,6 @@
    !! History :       !  06-00  (A.Estublier)  for passive tracers
    !!                 !  01-08  (E.Durand, G.Madec)  adapted for T & S
-   !!            8.5  !  02-06  (G. Madec)  F90: Free form and module
+   !!   NEMO     1.0  !  02-06  (G. Madec)  F90: Free form and module
+   !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC
    !!----------------------------------------------------------------------
 
@@ -13,5 +14,4 @@
    !!                   and vertical advection trends using MUSCL scheme
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and active tracers
    USE dom_oce         ! ocean space and time domain
    USE trdmod          ! ocean active tracers trends 
@@ -34,6 +34,6 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
-   !! $Header$ 
+   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 
+   !! $Id:$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -41,5 +41,6 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_muscl( kt, pun, pvn, pwn )
+   SUBROUTINE tra_adv_muscl( kt, cdtype, ktra, pun, pvn, pwn,   &
+      &                                        ptb     , pta )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE tra_adv_muscl  ***
@@ -51,5 +52,5 @@
       !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
       !!
-      !! ** Action  : - update (ta,sa) with the now advective tracer trends
+      !! ** Action  : - update (pta,sa) with the now advective tracer trends
       !!              - save trends in (ztrdt,ztrds) ('key_trdtra')
       !!
@@ -57,23 +58,18 @@
       !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
       !!----------------------------------------------------------------------
-      USE oce, ONLY :   ztrdt => ua   ! use ua as workspace
-      USE oce, ONLY :   ztrds => va   ! use va as workspace
-      !!
-      INTEGER , INTENT(in)                         ::   kt    ! ocean time-step index
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun   ! ocean velocity u-component
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pvn   ! ocean velocity v-component
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pwn   ! ocean velocity w-component
-      !!
-      INTEGER ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   &
-         zu, zv, zw, zeu, zev,           &  
-         zew, zbtr, zstep,               &
-         z0u, z0v, z0w,                  &
-         zzt1, zzt2, zalpha,             &
-         zzs1, zzs2, z2,                 &
-         zta, zsa,                       &
-        z_hdivn_x, z_hdivn_y, z_hdivn
-      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zt1, zt2, ztp1, ztp2   ! 3D workspace
-      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zs1, zs2, zsp1, zsp2   !  "      "
+      INTEGER         , INTENT(in   )                         ::   kt              ! ocean time-step index
+      CHARACTER(len=3), INTENT(in   )                         ::   cdtype          ! =TRA or TRC (tracer indicator)
+      INTEGER         , INTENT(in   )                         ::   ktra            ! tracer index
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn, pwn   ! 3 ocean velocity components
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptb             ! before tracer fields
+      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend 
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zu, zv, zw, zeu, zev  
+      REAL(wp) ::   zew, zbtr, z2, zstep 
+      REAL(wp) ::   z0u, z0v, z0w 
+      REAL(wp) ::   zzwx, zzwy, zalpha 
+      REAL(wp) ::   zta
+      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zwx, zwy, zslpx, zslpy   ! 3D workspace
       !!----------------------------------------------------------------------
 
@@ -95,18 +91,14 @@
          DO jj = 1, jpjm1      
             DO ji = 1, fs_jpim1   ! vector opt.
-               zt1(ji,jj,jk) = umask(ji,jj,jk) * ( tb(ji+1,jj,jk) - tb(ji,jj,jk) )
-               zs1(ji,jj,jk) = umask(ji,jj,jk) * ( sb(ji+1,jj,jk) - sb(ji,jj,jk) )
-               zt2(ji,jj,jk) = vmask(ji,jj,jk) * ( tb(ji,jj+1,jk) - tb(ji,jj,jk) )
-               zs2(ji,jj,jk) = vmask(ji,jj,jk) * ( sb(ji,jj+1,jk) - sb(ji,jj,jk) )
+               zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk) - ptb(ji,jj,jk) )
+               zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk) - ptb(ji,jj,jk) )
             END DO
          END DO
       END DO
       ! bottom values
-      zt1(:,:,jpk) = 0.e0    ;    zt2(:,:,jpk) = 0.e0
-      zs1(:,:,jpk) = 0.e0    ;    zs2(:,:,jpk) = 0.e0
-
-      ! lateral boundary conditions on zt1, zt2 ; zs1, zs2   (changed sign)
-      CALL lbc_lnk( zt1, 'U', -1. )   ;   CALL lbc_lnk( zs1, 'U', -1. )
-      CALL lbc_lnk( zt2, 'V', -1. )   ;   CALL lbc_lnk( zs2, 'V', -1. )
+      zwx(:,:,jpk) = 0.e0    ;    zwy(:,:,jpk) = 0.e0
+
+      ! lateral boundary conditions on zwx, zwy   (changed sign)
+      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )
 
       ! Slopes
@@ -115,18 +107,13 @@
          DO jj = 2, jpj
             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) ) )
+               zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
+                  &           * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
+               zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   &
+                  &           * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
             END DO
          END DO
       END DO
       ! bottom values
-      ztp1(:,:,jpk) = 0.e0    ;    ztp2(:,:,jpk) = 0.e0
-      zsp1(:,:,jpk) = 0.e0    ;    zsp2(:,:,jpk) = 0.e0
+      zslpx(:,:,jpk) = 0.e0    ;    zslpy(:,:,jpk) = 0.e0
 
       ! Slopes limitation
@@ -134,20 +121,12 @@
          DO jj = 2, jpj
             DO ji = fs_2, jpi   ! vector opt.
-               ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
-                  &           * MIN(    ABS( ztp1(ji  ,jj,jk) ),   &
-                  &                  2.*ABS( zt1 (ji-1,jj,jk) ),   &
-                  &                  2.*ABS( zt1 (ji  ,jj,jk) ) )
-               zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) )   &
-                  &           * MIN(    ABS( zsp1(ji  ,jj,jk) ),   &
-                  &                  2.*ABS( zs1 (ji-1,jj,jk) ),   &
-                  &                  2.*ABS( zs1 (ji  ,jj,jk) ) )
-               ztp2(ji,jj,jk) = SIGN( 1., ztp2(ji,jj,jk) )   &
-                  &           * MIN(    ABS( ztp2(ji,jj  ,jk) ),   &
-                  &                  2.*ABS( zt2 (ji,jj-1,jk) ),   &
-                  &                  2.*ABS( zt2 (ji,jj  ,jk) ) )
-               zsp2(ji,jj,jk) = SIGN( 1., zsp2(ji,jj,jk) )   &
-                  &           * MIN(    ABS( zsp2(ji,jj  ,jk) ),   &
-                  &                  2.*ABS( zs2 (ji,jj-1,jk) ),   &
-                  &                  2.*ABS( zs2 (ji,jj  ,jk) ) )
+               zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) )   &
+                  &            * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
+                  &                   2.*ABS( zwx (ji-1,jj,jk) ),   &
+                  &                   2.*ABS( zwx (ji  ,jj,jk) ) )
+               zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) )   &
+                  &            * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
+                  &                   2.*ABS( zwy (ji,jj-1,jk) ),   &
+                  &                   2.*ABS( zwy (ji,jj  ,jk) ) )
             END DO
          END DO
@@ -172,27 +151,21 @@
                zalpha = 0.5 - z0u
                zu  = z0u - 0.5 * pun(ji,jj,jk) * zstep / e1u(ji,jj)
-               zzt1 = tb(ji+1,jj,jk) + zu*ztp1(ji+1,jj,jk)
-               zzt2 = tb(ji  ,jj,jk) + zu*ztp1(ji  ,jj,jk)
-               zzs1 = sb(ji+1,jj,jk) + zu*zsp1(ji+1,jj,jk)
-               zzs2 = sb(ji  ,jj,jk) + zu*zsp1(ji  ,jj,jk)
-               zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
-               zs1(ji,jj,jk) = zeu * ( zalpha * zzs1 + (1.-zalpha) * zzs2 )
+               zzwx = ptb(ji+1,jj,jk) + zu * zslpx(ji+1,jj,jk)
+               zzwy = ptb(ji  ,jj,jk) + zu * zslpx(ji  ,jj,jk)
+               zwx(ji,jj,jk) = zeu * ( zalpha * zzwx + (1.-zalpha) * zzwy )
                !
                z0v = SIGN( 0.5, pvn(ji,jj,jk) )            
                zalpha = 0.5 - z0v
                zv  = z0v - 0.5 * pvn(ji,jj,jk) * zstep / e2v(ji,jj)
-               zzt1 = tb(ji,jj+1,jk) + zv*ztp2(ji,jj+1,jk)
-               zzt2 = tb(ji,jj  ,jk) + zv*ztp2(ji,jj  ,jk)
-               zzs1 = sb(ji,jj+1,jk) + zv*zsp2(ji,jj+1,jk)
-               zzs2 = sb(ji,jj  ,jk) + zv*zsp2(ji,jj  ,jk)
-               zt2(ji,jj,jk) = zev * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
-               zs2(ji,jj,jk) = zev * ( zalpha * zzs1 + (1.-zalpha) * zzs2 )
-            END DO
-         END DO
-      END DO
-
-      ! lateral boundary conditions on zt1, zt2 ; zs1, zs2   (changed sign)
-      CALL lbc_lnk( zt1, 'U', -1. )   ;   CALL lbc_lnk( zs1, 'U', -1. ) 
-      CALL lbc_lnk( zt2, 'V', -1. )   ;   CALL lbc_lnk( zs2, 'V', -1. )
+               zzwx = ptb(ji,jj+1,jk) + zv * zslpy(ji,jj+1,jk)
+               zzwy = ptb(ji,jj  ,jk) + zv * zslpy(ji,jj  ,jk)
+               zwy(ji,jj,jk) = zev * ( zalpha * zzwx + (1.-zalpha) * zzwy )
+            END DO
+         END DO
+      END DO
+
+!!gm bug?   there is too many lbc: this have to be checked
+      ! lateral boundary conditions on zwx, zwy   (changed sign)
+      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )
 
       ! Tracer flux divergence at t-point added to the general trend
@@ -201,92 +174,39 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
 #if defined key_zco
-               zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
+               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
 #else
-               zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
+               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
 #endif
                ! horizontal advective trends
-               zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk  )   &
-                  &           + zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk  ) )
-               zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj  ,jk  )   &
-                  &           + zs2(ji,jj,jk) - zs2(ji  ,jj-1,jk  ) ) 
+               zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                  &           + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) ) 
                ! add it to the general tracer trends
-               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
-               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
+               pta(ji,jj,jk) = pta(ji,jj,jk) + zta
             END DO
         END DO
       END DO        
 
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' muscl had  - Ta: ', mask1=tmask ,  &
-         &                       tab3d_2=sa, clinfo2=             ' Sa: ', mask2=tmask, clinfo3='tra' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - had: ', mask1=tmask, clinfo3=cdtype )
 
       ! Save the horizontal advective trends for diagnostics
       IF( l_trdtra ) THEN
-         ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0
-         !
-         ! T/S ZONAL advection trends
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  !-- Compute zonal divergence by splitting hdivn (see divcur.F90)
-                  !   N.B. This computation is not valid along OBCs (if any)
-#if defined key_zco
-                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
-                  z_hdivn_x = (  e2u(ji  ,jj) * pun(ji  ,jj,jk)                              &
-                     &         - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr
-#else
-                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-                  z_hdivn_x = (  e2u(ji  ,jj) * fse3u(ji  ,jj,jk) * pun(ji  ,jj,jk)          &
-                     &         - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr
-#endif
-                  ztrdt(ji,jj,jk) = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj,jk) ) + tn(ji,jj,jk) * z_hdivn_x
-                  ztrds(ji,jj,jk) = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj,jk) ) + sn(ji,jj,jk) * z_hdivn_x
-               END DO
-            END DO
-         END DO
-         CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt)
-
-         ! T/S MERIDIONAL advection trends
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  !-- Compute merid. divergence by splitting hdivn (see divcur.F90)
-                  !   N.B. This computation is not valid along OBCs (if any)
-#if defined key_zco
-                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
-                  z_hdivn_y = (  e1v(ji,jj  ) * pvn(ji,jj  ,jk)                              &
-                     &         - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr
-#else
-                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-                  z_hdivn_y = (  e1v(ji,  jj) * fse3v(ji,jj  ,jk) * pvn(ji,jj  ,jk)          &
-                     &         - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr
-#endif
-                  ztrdt(ji,jj,jk) = - zbtr * ( zt2(ji,jj,jk) - zt2(ji,jj-1,jk) ) + tn(ji,jj,jk) * z_hdivn_y          
-                  ztrds(ji,jj,jk) = - zbtr * ( zs2(ji,jj,jk) - zs2(ji,jj-1,jk) ) + sn(ji,jj,jk) * z_hdivn_y
-               END DO
-            END DO
-         END DO
-         CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt)
-
-         ! Save the up-to-date ta and sa trends
-         ztrdt(:,:,:) = ta(:,:,:) 
-         ztrds(:,:,:) = sa(:,:,:) 
-         !
-      ENDIF
-
-      ! "zonal" mean advective heat and salt transport
-      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
+         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptb )
+         CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptb )
+      ENDIF
+
+      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
          IF( lk_zco ) THEN
             DO jk = 1, jpkm1
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                    zt2(ji,jj,jk) = zt2(ji,jj,jk) * fse3v(ji,jj,jk)
-                    zs2(ji,jj,jk) = zs2(ji,jj,jk) * fse3v(ji,jj,jk)
+                    zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)
                   END DO
                END DO
             END DO
          ENDIF
-         pht_adv(:) = ptr_vj( zt2(:,:,:) )
-         pst_adv(:) = ptr_vj( zs2(:,:,:) )
-      ENDIF
+         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zwy(:,:,:) )
+         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zwy(:,:,:) )
+      ENDIF
+
 
       ! II. Vertical advective fluxes
@@ -296,10 +216,8 @@
       ! interior values
       DO jk = 2, jpkm1
-         zt1(:,:,jk) = tmask(:,:,jk) * ( tb(:,:,jk-1) - tb(:,:,jk) )
-         zs1(:,:,jk) = tmask(:,:,jk) * ( sb(:,:,jk-1) - sb(:,:,jk) )
+         zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1) - ptb(:,:,jk) )
       END DO
       ! surface & bottom boundary conditions
-      zt1 (:,:, 1 ) = 0.e0    ;    zt1 (:,:,jpk) = 0.e0
-      zs1 (:,:, 1 ) = 0.e0    ;    zs1 (:,:,jpk) = 0.e0
+      zwx (:,:, 1 ) = 0.e0    ;    zwx (:,:,jpk) = 0.e0
 
       ! Slopes
@@ -307,8 +225,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               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) ) )
+               zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )   &
+                  &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) )
             END DO
          END DO
@@ -316,26 +232,18 @@
 
       ! Slopes limitation
-      ! interior values
-      DO jk = 2, jpkm1
+      DO jk = 2, jpkm1        ! interior values
          DO jj = 1, jpj
             DO ji = 1, jpi
-               ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
-                  &           * MIN(    ABS( ztp1(ji,jj,jk  ) ),   &
-                  &                  2.*ABS( zt1 (ji,jj,jk+1) ),   &
-                  &                  2.*ABS( zt1 (ji,jj,jk  ) ) )
-               zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) )   &
-                  &           * MIN(    ABS( zsp1(ji,jj,jk  ) ),   &
-                  &                  2.*ABS( zs1 (ji,jj,jk+1) ),   &
-                  &                  2.*ABS( zs1 (ji,jj,jk  ) ) )
-            END DO
-         END DO
-      END DO
-      ! surface values
-      ztp1(:,:,1) = 0.e0
-      zsp1(:,:,1) = 0.e0
+               zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) )   &
+                  &            * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
+                  &                   2.*ABS( zwx (ji,jj,jk+1) ),   &
+                  &                   2.*ABS( zwx (ji,jj,jk  ) ) )
+            END DO
+         END DO
+      END DO
+      zslpx(:,:,1) = 0.e0      ! surface values
 
       ! vertical advective flux
-      ! interior values
-      DO jk = 1, jpkm1
+      DO jk = 1, jpkm1        ! interior values
          zstep  = z2 * rdttra(jk)
          DO jj = 2, jpjm1      
@@ -344,30 +252,22 @@
                z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
                zalpha = 0.5 + z0w
-               zw  = z0w - 0.5 * pwn(ji,jj,jk+1)*zstep / fse3w(ji,jj,jk+1)
-               zzt1 = tb(ji,jj,jk+1) + zw*ztp1(ji,jj,jk+1)
-               zzt2 = tb(ji,jj,jk  ) + zw*ztp1(ji,jj,jk  )
-               zzs1 = sb(ji,jj,jk+1) + zw*zsp1(ji,jj,jk+1)
-               zzs2 = sb(ji,jj,jk  ) + zw*zsp1(ji,jj,jk  )
-               zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 )
-               zs1(ji,jj,jk+1) = zew * ( zalpha * zzs1 + (1.-zalpha)*zzs2 )
-            END DO
-         END DO
-      END DO
-      ! surface values
+               zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1)
+               zzwx = ptb(ji,jj,jk+1) + zw * zslpx(ji,jj,jk+1)
+               zzwy = ptb(ji,jj,jk  ) + zw * zslpx(ji,jj,jk  )
+               zwx(ji,jj,jk+1) = zew * ( zalpha * zzwx + (1.-zalpha) * zzwy )
+            END DO
+         END DO
+      END DO
+      !                       ! surface values
       IF( lk_dynspg_rl .OR. lk_vvl) THEN                ! rigid lid or variable volume: flux set to zero
-         zt1(:,:, 1 ) = 0.e0
-         zs1(:,:, 1 ) = 0.e0
+         zwx(:,:, 1 ) = 0.e0                                 ! surface 
+         zwx(:,:,jpk) = 0.e0                                 ! bottom 
       ELSE                                              ! free surface
-         zt1(:,:, 1 ) = pwn(:,:,1) * tb(:,:,1)
-         zs1(:,:, 1 ) = pwn(:,:,1) * sb(:,:,1)
-      ENDIF
-
-      ! bottom values
-      zt1(:,:,jpk) = 0.e0
-      zs1(:,:,jpk) = 0.e0
-
+         zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1)              ! Surface
+         zwx(:,:,jpk) = 0.e0                                 ! bottom 
+
+      ENDIF
 
       ! Compute & add the vertical advective trend
-
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1      
@@ -375,9 +275,7 @@
                zbtr = 1. / fse3t(ji,jj,jk)
                ! horizontal advective trends
-               zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )
-               zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji,jj,jk+1) )
+               zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
                ! add it to the general tracer trends
-               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta
-               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa
+               pta(ji,jj,jk) =  pta(ji,jj,jk) + zta
             END DO
          END DO
@@ -386,34 +284,7 @@
       ! Save the vertical advective trends for diagnostic
       ! -------------------------------------------------
-      IF( l_trdtra )   THEN
-         ! Recompute the vertical advection zta & zsa trends computed 
-         ! at the step 2. above in making the difference between the new 
-         ! trends and the previous one: ta()/sa - ztrdt()/ztrds() and substract
-         ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
-
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-#if defined key_zco
-                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
-                  z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk)
-                  z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk)
-#else
-                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-                  z_hdivn_x = e2u(ji,jj)*fse3u(ji,jj,jk)*pun(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk)*pun(ji-1,jj,jk)
-                  z_hdivn_y = e1v(ji,jj)*fse3v(ji,jj,jk)*pvn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk)*pvn(ji,jj-1,jk)
-#endif
-                  z_hdivn   = (z_hdivn_x + z_hdivn_y) * zbtr
-                  ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn 
-                  ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn
-               END DO
-            END DO
-         END DO
-         CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt)
-         !
-      ENDIF
-
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' muscl zad  - Ta: ', mask1=tmask ,   &
-         &                       tab3d_2=sa, clinfo2=             ' Sa: ', mask2=tmask, clinfo3='tra' )
+      IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwx, pwn, ptb )
+
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - zad: ', mask1=tmask, clinfo3=cdtype )
       !
    END SUBROUTINE tra_adv_muscl
