Index: /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 790)
+++ /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 791)
@@ -4,5 +4,6 @@
    !! Ocean active tracers:  advection trend 
    !!==============================================================================
-   !! History :  9.0  !  05-11  (G. Madec)  Original code
+   !! History :  1.0  !  2005-11  (G. Madec)  Original code
+   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport
    !!----------------------------------------------------------------------
 
@@ -23,5 +24,4 @@
    USE ldftra_oce      ! lateral diffusion coefficient on tracers
    USE in_out_manager  ! I/O manager
-!  USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -59,59 +59,66 @@
       !! ** Method  : - Update (ua,va) with the advection term following nadv
       !!----------------------------------------------------------------------
-#if ( defined key_trabbl_adv || defined key_traldf_eiv )
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zun, zvn, zwn   ! effective velocity
-#else
-      USE oce, ONLY :                       zun => un       ! the effective velocity is the
-      USE oce, ONLY :                       zvn => vn       ! Eulerian velocity
-      USE oce, ONLY :                       zwn => wn       ! 
-#endif
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
       !!
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER ::   jk   ! dummy loop index
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zun, zvn, zwn   ! effective transports
       !!----------------------------------------------------------------------
 
       IF( kt == nit000 )   CALL tra_adv_ctl          ! initialisation & control of options
 
+
+      !                                              ! effective transport
+      DO jk = 1, jpkm1 
 #if defined key_trabbl_adv
-      zun(:,:,:) = un(:,:,:) - u_bbl(:,:,:)          ! add the bbl velocity
-      zvn(:,:,:) = vn(:,:,:) - v_bbl(:,:,:)
-      zwn(:,:,:) = wn(:,:,:) + w_bbl(:,:,:)
+         !                                                ! eulerian + bbl transport
+         zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * ( un(:,:,jk) - u_bbl(:,:,jk) )
+         zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * ( vn(:,:,jk) - v_bbl(:,:,jk) )
+         zwn(:,:,jk) = e1t(:,:) * e2t(:,:)      * ( wn(:,:,jk) + w_bbl(:,:,jk) )
+#else
+         !                                                ! eulerian transport only
+         zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) *   un(:,:,jk) 
+         zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) *   vn(:,:,jk)
+         zwn(:,:,jk) = e1t(:,:) * e2t(:,:)      *   wn(:,:,jk)
 #endif
-      IF( lk_traldf_eiv ) THEN                       ! commpute and add the eiv velocity
-         IF( .NOT. lk_trabbl_adv ) THEN 
-            zun(:,:,:) = un(:,:,:)
-            zvn(:,:,:) = vn(:,:,:)
-            zwn(:,:,:) = wn(:,:,:)
-         ENDIF
-         CALL tra_adv_eiv( kt, zun, zvn, zwn ) 
-      ENDIF
+      END DO
+      zwn(:,:,jpk) = 0.e0                                 ! no transport trough the bottom
+
+      !                                                   ! add the eiv transport (if necessary)
+      IF( lk_traldf_eiv )   CALL tra_adv_eiv( kt, zun, zvn, zwn ) 
+
 
       SELECT CASE ( nadv )                           ! compute advection trend and add it to general trend
-      CASE ( 0 )   ;   CALL tra_adv_cen2    ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! 2nd order centered
-                       CALL tra_adv_cen2    ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! 2nd order centered
-!     CASE ( 1 )   ;   CALL tra_adv_cen2_jki( kt, zun, zvn, zwn )    ! 2nd order centered scheme
+      !
+      CASE ( 1 )   ;   CALL tra_adv_cen2    ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! 2nd order centered
+                       CALL tra_adv_cen2    ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! on T & S
+      !
       CASE ( 2 )   ;   CALL tra_adv_tvd     ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! TVD      scheme
-                       CALL tra_adv_tvd     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! TVD      scheme
+                       CALL tra_adv_tvd     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! on T & S
+      !
       CASE ( 3 )   ;   CALL tra_adv_muscl   ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb    , ta )    ! MUSCL    scheme
-                       CALL tra_adv_muscl   ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb    , sa )    ! MUSCL    scheme
+                       CALL tra_adv_muscl   ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb    , sa )    ! on T & S
+      !
       CASE ( 4 )   ;   CALL tra_adv_muscl2  ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! MUSCL2   scheme
-                       CALL tra_adv_muscl2  ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! MUSCL2   scheme
+                       CALL tra_adv_muscl2  ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! on T & S
+      !
       CASE ( 5 )   ;   CALL tra_adv_ubs     ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! UBS      scheme
-                       CALL tra_adv_ubs     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! UBS      scheme
+                       CALL tra_adv_ubs     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! on T & S
+      !
       CASE ( 6 )   ;   CALL tra_adv_qck     ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! QUICKEST scheme
-                       CALL tra_adv_qck     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! QUICKEST scheme
+                       CALL tra_adv_qck     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! on T & S
       !
-      CASE (-1 )                                                     ! esopa: test all possibility with control print
+      CASE (-1 )                                     ! NEMO debug : test all the schemes at once
                        CALL tra_adv_cen2    ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! 2nd order centered
-                       CALL tra_adv_cen2    ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! 2nd order centered
+                       CALL tra_adv_cen2    ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) 
                        CALL tra_adv_tvd     ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! TVD      scheme
-                       CALL tra_adv_tvd     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! TVD      scheme
+                       CALL tra_adv_tvd     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )  
                        CALL tra_adv_muscl   ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb    , ta )    ! MUSCL    scheme
-                       CALL tra_adv_muscl   ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb    , sa )    ! MUSCL    scheme
+                       CALL tra_adv_muscl   ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb    , sa ) 
                        CALL tra_adv_muscl2  ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! MUSCL2   scheme
-                       CALL tra_adv_muscl2  ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! MUSCL2   scheme
+                       CALL tra_adv_muscl2  ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )
                        CALL tra_adv_ubs     ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! UBS      scheme
-                       CALL tra_adv_ubs     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! UBS      scheme
+                       CALL tra_adv_ubs     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )
                        CALL tra_adv_qck     ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! QUICKEST scheme
-                       CALL tra_adv_qck     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! QUICKEST scheme
+                       CALL tra_adv_qck     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )
       END SELECT
       !
@@ -127,5 +134,5 @@
       !!----------------------------------------------------------------------
       INTEGER ::   ioptio
-
+      !!
       NAMELIST/nam_traadv/ ln_traadv_cen2 , ln_traadv_tvd,    &
          &                 ln_traadv_muscl, ln_traadv_muscl2, &
@@ -163,9 +170,6 @@
          &                CALL ctl_stop( 'cross-land advection only with 2nd order advection scheme' )
 
-      !                              ! Set nadv
-      IF( ln_traadv_cen2   )   nadv =  0
-#if defined key_mpp_omp
+      !                               ! Set nadv
       IF( ln_traadv_cen2   )   nadv =  1
-#endif
       IF( ln_traadv_tvd    )   nadv =  2
       IF( ln_traadv_muscl  )   nadv =  3
@@ -175,8 +179,7 @@
       IF( lk_esopa         )   nadv = -1
 
-      IF(lwp) THEN                   ! Print the choice
+      IF(lwp) THEN                    ! Print the choice
          WRITE(numout,*)
-         IF( nadv ==  0 )   WRITE(numout,*) '         2nd order scheme is used'
-         IF( nadv ==  1 )   WRITE(numout,*) '         2nd order scheme is usedi, k-j-i case'
+         IF( nadv ==  1 )   WRITE(numout,*) '         2nd order scheme is used'
          IF( nadv ==  2 )   WRITE(numout,*) '         TVD       scheme is used'
          IF( nadv ==  3 )   WRITE(numout,*) '         MUSCL     scheme is used'
@@ -185,5 +188,5 @@
          IF( nadv ==  6 )   WRITE(numout,*) '         PPM       scheme is used'
          IF( nadv ==  7 )   WRITE(numout,*) '         QUICKEST  scheme is used'
-         IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection scheme'
+         IF( nadv == -1 )   WRITE(numout,*) '         NEMO debug: test all advection scheme at once'
       ENDIF
       !
Index: /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_cen2.F90
===================================================================
--- /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_cen2.F90	(revision 790)
+++ /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_cen2.F90	(revision 791)
@@ -4,10 +4,10 @@
    !! Ocean active tracers:  horizontal & vertical advective trend
    !!==============================================================================
-   !! History :  8.2  !  01-08  (G. Madec, E. Durand)  trahad+trazad = traadv 
-   !!            8.5  !  02-06  (G. Madec)  F90: Free form and module
-   !!   NEMO     1.0  !  05-11  (V. Garnier) Surface pressure gradient organization
-   !!             -   !  05-11  (V. Garnier) Surface pressure gradient organization
-   !!             -   !  06-04  (R. Benshila, G. Madec) Step reorganization
-   !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC
+   !! History :  8.2  !  2001-08  (G. Madec, E. Durand)  trahad+trazad = traadv 
+   !!            8.5  !  2002-06  (G. Madec)  F90: Free form and module
+   !!   NEMO     1.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
+   !!             -   !  2005-11  (V. Garnier) Surface pressure gradient organization
+   !!             -   !  2006-04  (R. Benshila, G. Madec) Step reorganization
+   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport
    !!----------------------------------------------------------------------
 
@@ -34,6 +34,4 @@
    PUBLIC   tra_adv_cen2   ! routine called by step.F90
 
-   REAL(wp), DIMENSION(jpi,jpj) ::   btr2   ! inverse of T-point surface [1/(e1t*e2t)]
-
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -41,5 +39,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 
-   !! $Id:$ 
+   !! $Id$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -62,50 +60,24 @@
       !!         Part 0 : compute the upstream / centered flag
       !!                  (3D array, zind, defined at T-point (0<zind<1))
-      !!         Part I : horizontal advection
-      !!       * centered flux:
-      !!               zcenu = e2u*e3u  un  mi(ptn)
-      !!               zcenv = e1v*e3v  vn  mj(ptn)
-      !!       * upstream flux:
-      !!               zupsu = e2u*e3u  un  (ptb(i) or ptb(i-1) ) [un>0 or <0]
-      !!               zupsv = e1v*e3v  vn  (tb(j) or tb(j-1) ) [vn>0 or <0]
-      !!       * mixed upstream / centered horizontal advection scheme
+      !!         Part I : advective fluxes
+      !!       * centered flux:   zcenu = pun  mi(ptn)         (idem for v and w)
+      !!       * upstream flux:   zupsu = pun  (ptb(i) or ptb(i-1) ) [un>0 or <0]
+      !!       * mixed upstream / centered flux:
       !!               zcofi = max(zind(i+1), zind(i))
-      !!               zcofj = max(zind(j+1), zind(j))
       !!               zwx = zcofi * zupsu + (1-zcofi) * zcenu
-      !!               zwy = zcofj * zupsv + (1-zcofj) * zcenv
+      !!         Part III : advective trend
       !!       * horizontal advective trend (divergence of the fluxes)
-      !!               zta = 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] }
-      !!       * Add this trend now to the general trend of tracer (pta):
-      !!               pta = pta + zta 
-      !!       * trend diagnostic (lk_trdtra=T): the trend is
-      !!      saved for diagnostics. The trends saved is expressed as
-      !!      Uh.gradh(T), i.e.
-      !!                     save trend = zta + ptn divn
-      !!         In addition, the advective trend in the two horizontal direc-
-      !!      tion is also re-computed as Uh gradh(T). Indeed hadt+ptn divn is
-      !!      equal to (in s-coordinates, and similarly in z-coord.):
-      !!         zta+ptn*divn=1/(e1t*e2t*e3t) { mi-1( e2u*e3u  un  di[ptn] )
-      !!                                       +mj-1( e1v*e3v  vn  mj[ptn] )  }
-      !!         NB:in z-coordinate - full step (ln_zco=T) e3u=e3v=e3t, so
-      !!      they vanish from the expression of the flux and divergence.
-      !!
-      !!         Part II : vertical advection
-      !!      the advective trend is computed as follows :
-      !!            zta = 1/e3t dk+1[ zwx ]
-      !!      where the vertical advective flux, zwx, is given by :
-      !!            zwx = zcofk * zupst + (1-zcofk) * zcent
-      !!      with
-      !!        zupsv = upstream flux = wn * (tb(k) or tb(k-1) ) [wn>0 or <0]
-      !!        zcenu = centered flux = wn * mk(ptn)
-      !!         The surface boundary condition is :
-      !!      rigid-lid (lk_dynspg_frd = T) : zero advective flux
-      !!      free-surf (lk_dynspg_fsc = T) : wn(:,:,1) * ptn(:,:,1)
-      !!         Add this trend now to the general trend of tracer (ta,sa):
-      !!            (ta,sa) = (ta,sa) + ( zta , zsa )
-      !!         Trend diagnostic ('key_trdtra' defined): the trend is
-      !!      saved for diagnostics. The trends saved is expressed as :
-      !!             save trend =  w.gradz(T) = zta - ptn divn.
-      !!
-      !! ** Action :  - update (ta,sa) with the now advective tracer trends
+      !!         added to the general trend of tracer (pta):
+      !!               pta = pta - 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] + dk[zwz] }
+      !!       * trend diagnostic (lk_trdtra=T): the advective fluxes are
+      !!      where the 3 part of advective trend are recomputed as U.grad[T]
+      !!      and then used in on-line diagnostics.
+      !!
+      !!         NB: The surface boundary condition is no flux for both
+      !!      rigid-lid (lk_dynspg_frd = T) and non-linear free-surface
+      !!      but a flux (pwn(1) * ptn(1)) cross the surface in linear
+      !!      free-surface case (lk_dynspg_fsc = T)
+      !!
+      !! ** Action :  - update pta with the 2nd order centered advective trends
       !!              - trend diagnostics (lk_trdtra=T)
       !!----------------------------------------------------------------------
@@ -118,9 +90,9 @@
       !!
       INTEGER  ::   ji, jj, jk                                    ! dummy loop indices
-      REAL(wp) ::   zbtr, zta, zhw, ze3tr                         ! temporary scalars
-      REAL(wp) ::   zcofi, zfui, zcenut, zupsut, zfp_ui, zfm_ui   !    "         "
-      REAL(wp) ::   zcofj, zfvj, zcenvt, zupsvt, zfp_vj, zfm_vj   !    "         "
-      REAL(wp) ::   zcofk,       zcent , zupst , zfp_w , zfm_w    !    "         "
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zwy, zind   ! 3D workspace 
+      REAL(wp) ::   zbtr                                          ! temporary scalars
+      REAL(wp) ::   zcofi, zcenut, zupsut, zfp_ui, zfm_ui         !    "         "
+      REAL(wp) ::   zcofj, zcenvt, zupsvt, zfp_vj, zfm_vj         !    "         "
+      REAL(wp) ::   zcofk, zcent , zupst , zfp_w , zfm_w          !    "         "
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zwy, zwz, zind   ! 3D workspace 
       !!----------------------------------------------------------------------
 
@@ -129,6 +101,4 @@
          IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
-         ! 
-         btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
       ENDIF
 
@@ -150,132 +120,81 @@
       END DO
 
-      ! I. Horizontal advection
-      ! -----------------------
-
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         !
-         DO jj = 1, jpjm1                              !  Horizontal advective fluxes
+      ! I. Advective fluxes
+      ! -------------------
+      DO jk = 1, jpkm1                                 ! Horizontal advective fluxes
+         DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
                ! upstream indicator
                zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
                zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
-               ! volume fluxes * 1/2
-#if defined key_zco
-               zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk)
-               zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk)
-#else
-               zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
-               zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
-#endif
                ! upstream scheme
-               zfp_ui = zfui + ABS( zfui )
-               zfp_vj = zfvj + ABS( zfvj )
-               zfm_ui = zfui - ABS( zfui )
-               zfm_vj = zfvj - ABS( zfvj )
+               zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
+               zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
+               zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
+               zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
                zupsut = zfp_ui * ptb(ji,jj,jk) + zfm_ui * ptb(ji+1,jj  ,jk)
                zupsvt = zfp_vj * ptb(ji,jj,jk) + zfm_vj * ptb(ji  ,jj+1,jk)
                ! centered scheme
-               zcenut = zfui * ( ptn(ji,jj,jk) + ptn(ji+1,jj  ,jk) )
-               zcenvt = zfvj * ( ptn(ji,jj,jk) + ptn(ji  ,jj+1,jk) )
+               zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji+1,jj  ,jk) )
+               zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji  ,jj+1,jk) )
                ! mixed centered / upstream scheme
-               zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut
-               zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt
-            END DO
-         END DO
-
-         DO jj = 2, jpjm1                              !  horizontal tracer flux divergence added to the general trend
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-#if defined key_zco
-               zbtr = btr2(ji,jj)
-#else
-               zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
-#endif
-               ! horizontal advective trends
-               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
-               pta(ji,jj,jk) = pta(ji,jj,jk) + zta
-            END DO
-         END DO
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-
-      !  Save the horizontal advective trends for diagnostic
-      ! -----------------------------------------------------
-      IF( l_trdtra ) THEN
-         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptn )
-         CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptn )
-      ENDIF
-
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' cen2 - had: ', mask1=tmask, clinfo3=cdtype )
-
-
-      ! "Poleward" heat and salt transport 
-      ! ----------------------------------
-      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.
-                    zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-         ENDIF
-         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zwy(:,:,:) )
-         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zwy(:,:,:) )
-      ENDIF
-
-
-      ! II. Vertical advection
-      ! ----------------------
-
-      IF( lk_dynspg_rl .OR. lk_vvl ) THEN         ! rigid lid or non-linear free surface
-         zwx(:,:, 1 ) = 0.e0                           ! Surface value : zero flux
-         zwx(:,:,jpk) = 0.e0                           ! Bottom  value : flux set to zero
-      ELSE                                        ! linear free surface
-         zwx(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1)        ! Surface :  : advection through z=0
-         zwx(:,:,jpk) = 0.e0                           ! Bottom  : flux set to zero
-      ENDIF
-
-      DO jk = 2, jpk                                   ! Vertical advective fluxes (at w-point)
+               zwx(ji,jj,jk) = 0.5 * ( zcofi * zupsut + (1.-zcofi) * zcenut )
+               zwy(ji,jj,jk) = 0.5 * ( zcofj * zupsvt + (1.-zcofj) * zcenvt )
+            END DO
+         END DO
+      END DO
+      !                                                ! Vertical advective fluxes
+      zwz(:,:,jpk) = 0.e0                                   ! Bottom  value : flux set to zero
+      !                                                     ! Surface value : 
+      IF( lk_dynspg_rl .OR. lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                      ! rigid lid or non-linear fs
+      ELSE                                  ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1)   ! linear free surface 
+      ENDIF
+      !
+      DO jk = 2, jpkm1                                      ! interior values
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
                ! upstream indicator
                zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )
-               ! velocity * 1/2
-               zhw = 0.5 * pwn(ji,jj,jk)
                ! upstream scheme
-               zfp_w = zhw + ABS( zhw )
-               zfm_w = zhw - ABS( zhw )
+               zfp_w = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
+               zfm_w = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
                zupst = zfp_w * ptb(ji,jj,jk) + zfm_w * ptb(ji,jj,jk-1)
                ! centered scheme
-               zcent = zhw * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) )
+               zcent = pwn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) )
                ! mixed centered / upstream scheme
-               zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent
-            END DO
-         END DO
-      END DO
-
-      DO jk = 1, jpkm1                                 ! Tracer flux divergence at t-point added to the general trend
-         DO jj = 2, jpjm1
+               zwz(ji,jj,jk) = 0.5 * ( zcofk * zupst + (1.-zcofk) * zcent )
+            END DO
+         END DO
+      END DO
+
+      ! II. Divergence of advective fluxes
+      ! ----------------------------------
+      DO jk = 1, jpkm1                   
+         DO jj = 2, jpjm1      
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ze3tr = 1. / fse3t(ji,jj,jk)
-               ! vertical advective trends
-               zta = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
-               ! add it to the general tracer trends
-               pta(ji,jj,jk) =  pta(ji,jj,jk) + zta
-            END DO
-         END DO
-      END DO
-
-      ! 3. Save the vertical advective trends for diagnostic
-      ! ----------------------------------------------------
-      IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwx, pwn, ptn )
-
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' cen2 - zad : ', mask1=tmask, clinfo3=cdtype )
+               zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+               ! advective trends added to the general tracer trends
+               pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                  &                                    + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                  &                                    + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  )
+            END DO
+         END DO
+      END DO
+
+
+      ! III. Diagnostics and control print
+      ! ----------------------------------
+      IF( l_trdtra ) THEN                            ! trend diagnostics
+         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptn )
+         CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptn )
+         CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwz, pwn, ptn )
+      ENDIF
+      !                                              ! "Poleward" heat and salt transport 
+      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
+         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zwy(:,:,:) )
+         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zwy(:,:,:) )
+      ENDIF
+      !                                              ! control print
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' cen2 - adv: ', mask1=tmask, clinfo3=cdtype )
       !
    END SUBROUTINE tra_adv_cen2
Index: /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_eiv.F90
===================================================================
--- /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_eiv.F90	(revision 790)
+++ /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_eiv.F90	(revision 791)
@@ -5,4 +5,5 @@
    !!======================================================================
    !! History :  9.0  !  05-11  (G. Madec)  Original code, from traldf and zdf _iso
+   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport
    !!----------------------------------------------------------------------
 #if defined key_traldf_eiv   ||   defined key_esopa
@@ -32,6 +33,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)
    !!----------------------------------------------------------------------
@@ -43,21 +44,22 @@
       !!                  ***  ROUTINE tra_adv_eiv  ***
       !! 
-      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive 
-      !!      trend and add it to the general trend of tracer equation.
-      !!
-      !! ** Method  :   The eddy induced advection is computed from the slope
+      !! ** Purpose :   Compute the eddy induced transport and add it to the
+      !!              effective transport
+      !!
+      !! ** Method  :   The eddy induced transport is computed from the slope
       !!      of iso-neutral surfaces computed in routine ldf_slp as follows:
-      !!         zu_eiv =  1/(e2u e3u)   dk[ aeiu e2u mi(wslpi) ]
-      !!         zv_eiv =  1/(e1v e3v)   dk[ aeiv e1v mj(wslpj)
-      !!         zw_eiv = -1/(e1t e2t) { di[ aeiu e2u mi(wslpi) ]
-      !!                               + dj[ aeiv e1v mj(wslpj) ] }
+      !!         zu_eiv =   dk[ aeiu e2u mi(wslpi) ]
+      !!         zv_eiv =   dk[ aeiv e1v mj(wslpj) ]
+      !!         zw_eiv = - { di[ aeiu e2u mi(wslpi) ] + dj[ aeiv e1v mj(wslpj) ] }
       !!      add the eiv component to the model velocity:
       !!         p.n = p.n + z._eiv
+      !! CAUTION : the horizontal transports not updated along jpi column and jpj row
+      !!           the vertical   transports not updated along 1 & jpi columns and 1 & jpj rows
       !!
       !! ** Action  : - add to p.n the eiv component
       !!----------------------------------------------------------------------
       INTEGER , INTENT(in   )                         ::   kt    ! ocean time-step index
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun   ! in : 3 ocean velocity components 
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pvn   ! out: 3 ocean velocity components
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun   ! in : 3 ocean transport components 
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pvn   ! out: 3 ocean transport components
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pwn   !      increased by the eiv
       !!
@@ -88,12 +90,12 @@
                zvwk1= ( wslpj(ji,jj,jk+1) + wslpj(ji,jj+1,jk+1) ) * fsaeiv(ji,jj,jk+1) * vmask(ji,jj,jk+1)
 
-               zu_eiv = 0.5 * umask(ji,jj,jk) * ( zuwk - zuwk1 ) / fse3u(ji,jj,jk)
-               zv_eiv = 0.5 * vmask(ji,jj,jk) * ( zvwk - zvwk1 ) / fse3v(ji,jj,jk)
+               zu_eiv = 0.5 * umask(ji,jj,jk) * ( zuwk - zuwk1 )
+               zv_eiv = 0.5 * vmask(ji,jj,jk) * ( zvwk - zvwk1 )
    
-               pun(ji,jj,jk) = pun(ji,jj,jk) + zu_eiv
-               pvn(ji,jj,jk) = pvn(ji,jj,jk) + zv_eiv
-# if defined key_diaeiv
-               u_eiv(ji,jj,jk) = zu_eiv
-               v_eiv(ji,jj,jk) = zv_eiv
+               pun(ji,jj,jk) = pun(ji,jj,jk) + e2u(ji,jj) * zu_eiv
+               pvn(ji,jj,jk) = pvn(ji,jj,jk) + e1v(ji,jj) * zv_eiv
+# if defined key_diaeiv
+               u_eiv(ji,jj,jk) = zu_eiv / fse3u(ji,jj,jk)
+               v_eiv(ji,jj,jk) = zv_eiv / fse3v(ji,jj,jk)
 # endif
             END DO
@@ -115,11 +117,10 @@
                   zvwj1 = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) ) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
 
-                  zw_eiv = - 0.5 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj )  &
-                     &                                                / ( e1t(ji,jj)*e2t(ji,jj) )
+                  zw_eiv = - 0.5 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj )
 # endif
                   pwn(ji,jj,jk) = pwn(ji,jj,jk) + zw_eiv
 
 # if defined key_diaeiv
-                  w_eiv(ji,jj,jk) = zw_eiv
+                  w_eiv(ji,jj,jk) = zw_eiv / ( e1t(ji,jj)*e2t(ji,jj) )
 # endif
                END DO
@@ -130,4 +131,84 @@
       !                                             ! =================
    END SUBROUTINE tra_adv_eiv
+
+
+!!gm   test tra_adv_eiv better (faster) coded?  to be verified
+
+   SUBROUTINE tra_adv_eiv2( kt, pun, pvn, pwn )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_adv_eiv  ***
+      !! 
+      !! ** Purpose :   Compute the eddy induced transport and add it to the
+      !!              effective transport
+      !!
+      !! ** Method  :   The eddy induced transport is computed from the slope
+      !!              of iso-neutral surfaces (see ldfslp.F90) as follows:
+      !!                   zu_eiv =   dk[ aeiu e2u mi(wslpi) ]
+      !!                   zv_eiv =   dk[ aeiv e1v mj(wslpj) ]
+      !!                   zw_eiv = - { di[ aeiu e2u mi(wslpi) ] + dj[ aeiv e1v mj(wslpj) ] }
+      !!              add the eiv component to the model velocity:
+      !!                   p.n = p.n + z._eiv
+      !! CAUTION : the horizontal transports not updated along jpi column and jpj row
+      !!           the vertical   transports not updated along 1 & jpi columns and 1 & jpj rows
+      !!
+      !! ** Action  : - add to p.n the eiv transport component
+      !!----------------------------------------------------------------------
+      INTEGER , INTENT(in   )                         ::   kt    ! ocean time-step index
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun   ! in : 3 ocean transport components 
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pvn   ! out: 3 ocean transport components
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pwn   !      increased by the eiv
+      !!
+      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
+      REAL(wp) ::   zuwslpi, zvwslpj           ! temporary scalar
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsfu, zsfv   ! eiv stream-function in u and v directions
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_adv_eiv : eddy induced advection :'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   add to velocity fields the eiv component'
+      ENDIF
+
+      ! eiv stream-function in u- and v-directions
+      ! NB: UW-point mask at level k is umask(:,:,k) idem form VW-point mask
+      zsfu(:,:, 1 ) = 0.e0   ;   zsfv(:,:, 1 ) = 0.e0      ! surface value set to zero
+      zsfu(:,:,jpk) = 0.e0   ;   zsfv(:,:,jpk) = 0.e0      ! bottom  value set to zero
+      DO jk = 2, jpkm1
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zuwslpi = 0.5 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
+               zvwslpj = 0.5 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) ) 
+               zsfu(ji,jj,jk) = zuwslpi * e2u(ji,jj) * fsaeiu(ji,jj,jk) * umask(ji,jj,jk)
+               zsfv(ji,jj,jk) = zvwslpj * e1v(ji,jj) * fsaeiv(ji,jj,jk) * vmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+# if defined key_diaeiv
+      ! save eiv stream function in the output
+!!gm  to be done, u_sfeiv and v_sfeiv not defined    ==> new IOM....
+!!gm  and zsfu, zfv notdefined for jpi column and jpj row
+!     u_sfeiv(:,:,:) = zsfu(:,:,:)
+!     v_sfeiv(:,:,:) = zsfu(:,:,:)
+# endif
+
+      ! increase the transport with the eiv transport
+      DO jk = 1, jpkm1
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               pun(ji,jj,jk) = pun(ji,jj,jk) + ( zsfu(ji,jj,jk) - zsfu(ji,jj,jk+1) )
+               pvn(ji,jj,jk) = pvn(ji,jj,jk) + ( zsfv(ji,jj,jk) - zsfv(ji,jj,jk+1) )
+            END DO
+         END DO
+      END DO
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               pwn(ji,jj,jk) = pwn(ji,jj,jk) - (  zsfu(ji,jj,jk) - zsfu(ji-1,jj  ,jk)                        &
+                  &                             + zsfv(ji,jj,jk) - zsfv(ji  ,jj-1,jk)  ) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+   END SUBROUTINE tra_adv_eiv2
 
 #else
Index: /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl.F90
===================================================================
--- /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl.F90	(revision 790)
+++ /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl.F90	(revision 791)
@@ -4,8 +4,8 @@
    !! Ocean active tracers:  horizontal & vertical advective trend
    !!======================================================================
-   !! History :       !  06-00  (A.Estublier)  for passive tracers
-   !!                 !  01-08  (E.Durand, G.Madec)  adapted for T & S
-   !!   NEMO     1.0  !  02-06  (G. Madec)  F90: Free form and module
-   !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC
+   !! History :  OPA  !  2006-00  (A.Estublier)  for passive tracers
+   !!            8.2  !  2001-08  (E.Durand, G.Madec)  adapted for T & S
+   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
+   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport
    !!----------------------------------------------------------------------
 
@@ -28,5 +28,5 @@
    PRIVATE
 
-   PUBLIC   tra_adv_muscl  ! routine called by step.F90
+   PUBLIC   tra_adv_muscl   ! routine called by step.F90
 
    !! * Substitutions
@@ -35,5 +35,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 
-   !! $Id:$ 
+   !! $Id$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -66,9 +66,8 @@
       !!
       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) ::   zu, z0u, zzwx
+      REAL(wp) ::   zv, z0v, zzwy 
+      REAL(wp) ::   zw, z0w 
+      REAL(wp) ::   zbtr, z2, zdt, zalpha 
       REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zwx, zwy, zslpx, zslpy   ! 3D workspace
       !!----------------------------------------------------------------------
@@ -80,13 +79,13 @@
       ENDIF
 
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
-      ELSE                                        ;    z2 = 2.
+      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.      ! euler     time-stepping
+      ELSE                                        ;    z2 = 2.      ! leap-frog time-stepping
       ENDIF
 
       ! I. Horizontal advective fluxes
       ! ------------------------------
-      ! first guess of the slopes
-      ! interior values
-      DO jk = 1, jpkm1
+      !                                             !-- first guess of the slopes
+      zwx(:,:,jpk) = 0.e0   ;   zwy(:,:,jpk) = 0.e0        ! bottom values
+      DO jk = 1, jpkm1                                     ! interior values
          DO jj = 1, jpjm1      
             DO ji = 1, fs_jpim1   ! vector opt.
@@ -96,116 +95,91 @@
          END DO
       END DO
-      ! bottom values
-      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
-      ! interior values
-      DO jk = 1, jpkm1
-         DO jj = 2, jpj
-            DO ji = fs_2, jpi   ! vector opt.
+!!gm optimisation (lbc to be suppressed)
+      CALL lbc_lnk( zwx, 'U', -1. )                        ! lateral boundary conditions on zwx, zwy   (changed sign)
+      CALL lbc_lnk( zwy, 'V', -1. )
+!!gm
+
+      !                                             !-- Slopes of tracer
+      zslpx(:,:,jpk) = 0.e0   ;   zslpy(:,:,jpk) = 0.e0    ! bottom values
+      DO jk = 1, jpkm1                                     ! interior values
+         DO jj = 2, jpjm1
+            DO ji = fs_2, jpim1   ! vector opt.
                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) ) )
+                  &            * ( 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
-      zslpx(:,:,jpk) = 0.e0    ;    zslpy(:,:,jpk) = 0.e0
-
-      ! Slopes limitation
-      DO jk = 1, jpkm1
-         DO jj = 2, jpj
-            DO ji = fs_2, jpi   ! vector opt.
-               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) ) )
+                  &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
+            END DO
+         END DO
+      END DO
+!!gm  :merge the 2 loop (above & below)
+      DO jk = 1, jpkm1                                     ! Slopes limitation
+         DO jj = 2, jpjm1
+            DO ji = fs_2, jpim1   ! vector opt.
+               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
       END DO        
-
-      ! Advection terms
-      ! interior values
-      DO jk = 1, jpkm1
-         zstep  = z2 * rdttra(jk)
-         DO jj = 2, jpjm1      
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ! volume fluxes
-#if defined key_zco
-               zeu = e2u(ji,jj)                   * pun(ji,jj,jk)
-               zev = e1v(ji,jj)                   * pvn(ji,jj,jk)
-#else
-               zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
-               zev = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
-#endif
-               ! MUSCL fluxes
+!!gm optimisation:  call lbclnk just 1 time here on zslpx & zslpy
+!!    CALL lbc_lnk( zslpx, 'U', -1. )                      ! lateral boundary conditions on zwx, zwy   (changed sign)
+!!    CALL lbc_lnk( zslpy, 'V', -1. )
+!!
+!! and loop below from 1 to jpjm1 and to jpim1
+!!gm
+
+      !                                             !-- MUSCL horizontal advective fluxes
+      DO jk = 1, jpkm1                                     ! interior values
+         zdt  = z2 * rdttra(jk)
+         DO jj = 2, jpjm1      
+            DO ji = fs_2, fs_jpim1   ! vector opt.
                z0u = SIGN( 0.5, pun(ji,jj,jk) )            
                zalpha = 0.5 - z0u
-               zu  = z0u - 0.5 * pun(ji,jj,jk) * zstep / e1u(ji,jj)
+               zu  = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
                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 )
+               zwx(ji,jj,jk) = pun(ji,jj,jk) * ( 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)
+               zv  = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
                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)
+               zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
+            END DO
+         END DO
+      END DO
+!!gm optimisation (lbc to be suppressed)
+      !                                                    ! lateral boundary conditions on zwx, zwy   (changed sign)
       CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )
-
+!!gm
+
+      !                                             !-- MUSCL advective trend
       ! Tracer flux divergence at t-point added to the general trend
       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) )
-#else
                zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-#endif
-               ! horizontal advective trends
-               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
-               pta(ji,jj,jk) = pta(ji,jj,jk) + zta
+               ! horizontal advective trends added to the general tracer trends
+               pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                  &                                   + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) ) 
             END DO
         END DO
       END DO        
 
-      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
+      IF( l_trdtra ) THEN                           !-- trend diagnostics
          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
-
+      !                                             !-- "Poleward" heat or salt transports
       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.
-                    zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-         ENDIF
          IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zwy(:,:,:) )
          IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zwy(:,:,:) )
       ENDIF
+      !                                             !-- control print
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - had: ', mask1=tmask, clinfo3=cdtype )
 
 
@@ -213,14 +187,13 @@
       ! -----------------------------
       
-      ! First guess of the slope
-      ! interior values
-      DO jk = 2, jpkm1
+      !                                             !-- first guess of the slopes
+      zwx (:,:, 1 ) = 0.e0    ;    zwx (:,:,jpk) = 0.e0    ! surface & bottom boundary conditions
+      DO jk = 2, jpkm1                                     ! interior values
          zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1) - ptb(:,:,jk) )
       END DO
-      ! surface & bottom boundary conditions
-      zwx (:,:, 1 ) = 0.e0    ;    zwx (:,:,jpk) = 0.e0
-
-      ! Slopes
-      DO jk = 2, jpkm1
+
+      !                                             !-- Slopes of tracer
+      zslpx(:,:,1) = 0.e0                                  ! surface values
+      DO jk = 2, jpkm1                                     ! interior value
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -230,60 +203,49 @@
          END DO
       END DO
-
-      ! Slopes limitation
-      DO jk = 2, jpkm1        ! interior values
+!!gm : merge the above and below loops
+      !                                             !-- Slopes limitation
+      DO jk = 2, jpkm1                                     ! interior values
          DO jj = 1, jpj
             DO ji = 1, jpi
-               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
-      DO jk = 1, jpkm1        ! interior values
-         zstep  = z2 * rdttra(jk)
-         DO jj = 2, jpjm1      
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zew = pwn(ji,jj,jk+1)
+               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
+
+      !                                             !-- vertical advective flux
+      !                                                    ! surface values  (bottom already set to zero)
+      IF( lk_dynspg_rl .OR. lk_vvl) THEN    ;   zwx(:,:, 1 ) = 0.e0                      ! rigid lid or variable volume
+      ELSE                                  ;   zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1)   ! linear free surface
+      ENDIF
+      DO jk = 1, jpkm1                                     ! interior values
+         zdt  = z2 * rdttra(jk)
+         DO jj = 2, jpjm1      
+            DO ji = fs_2, fs_jpim1   ! vector opt.
                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)
+               zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt / (e1t(ji,jj) * e2t(ji,jj) * 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
-         zwx(:,:, 1 ) = 0.e0                                 ! surface 
-         zwx(:,:,jpk) = 0.e0                                 ! bottom 
-      ELSE                                              ! free surface
-         zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1)              ! Surface
-         zwx(:,:,jpk) = 0.e0                                 ! bottom 
-
-      ENDIF
-
-      ! Compute & add the vertical advective trend
+               zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
+            END DO
+         END DO
+      END DO
+
+      !                                             !-- vertical advective trend
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1      
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zbtr = 1. / fse3t(ji,jj,jk)
-               ! horizontal advective trends
-               zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
-               ! add it to the general tracer trends
-               pta(ji,jj,jk) =  pta(ji,jj,jk) + zta
-            END DO
-         END DO
-      END DO
-
-      ! Save the vertical advective trends for diagnostic
-      ! -------------------------------------------------
+               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+               ! horizontal advective trends added to the general tracer trends
+               pta(ji,jj,jk) =  pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
+            END DO
+         END DO
+      END DO
+
+      !                                             !-- trend diagnostic
       IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwx, pwn, ptb )
-
+      !                                             !-- control print
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - zad: ', mask1=tmask, clinfo3=cdtype )
       !
Index: /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90
===================================================================
--- /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90	(revision 790)
+++ /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90	(revision 791)
@@ -4,6 +4,6 @@
    !! Ocean active tracers:  horizontal & vertical advective trend
    !!==============================================================================
-   !! History :  1.0  !  02-06  (G. Madec) from traadv_muscl
-   !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC
+   !! History :  1.0  !  2002-06  (G. Madec) from traadv_muscl
+   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport
    !!----------------------------------------------------------------------
 
@@ -26,5 +26,4 @@
    PRIVATE
 
-   !! * Accessibility
    PUBLIC tra_adv_muscl2        ! routine called by step.F90
 
@@ -34,5 +33,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 
-   !! $Id:$ 
+   !! $Id$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -66,8 +65,8 @@
       !!
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zu, zv, zw, zeu, zev
-      REAL(wp) ::   zew, zbtr, zstep, z2 
+      REAL(wp) ::   zu, zv, zw
+      REAL(wp) ::   zbtr, zstep, z2 
       REAL(wp) ::   z0u, z0v, z0w
-      REAL(wp) ::   zzwx, zzwy, zalpha, zta
+      REAL(wp) ::   zzwx, zzwy, zalpha
       REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zwx, zwy, zslpx, zslpy   ! 3D workspace
       !!----------------------------------------------------------------------
@@ -98,4 +97,5 @@
       zwx(:,:,jpk) = 0.e0    ;    zwy(:,:,jpk) = 0.e0        ! bottom values
 
+!!gm  this lbc_lnk can be omitted
       ! lateral boundary conditions on zwx, zwy   (changed sign)
       CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )
@@ -117,12 +117,10 @@
          DO jj = 2, jpj
             DO ji = fs_2, jpi   ! vector opt.
-               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) ) )
+               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
@@ -134,26 +132,17 @@
          DO jj = 2, jpjm1      
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ! volume fluxes
-#if defined key_zco
-               zeu = e2u(ji,jj)                   * pun(ji,jj,jk)
-               zev = e1v(ji,jj)                   * pvn(ji,jj,jk)
-#else
-               zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
-               zev = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
-#endif
-               ! MUSCL fluxes
-               z0u = SIGN( 0.5, pun(ji,jj,jk) )            
+               z0u  = SIGN( 0.5, pun(ji,jj,jk) )            
                zalpha = 0.5 - z0u
-               zu  = z0u - 0.5 * pun(ji,jj,jk) * zstep / e1u(ji,jj)
-               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 )
+               zu   = z0u - 0.5 * pun(ji,jj,jk) * zstep / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
+               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) = pun(ji,jj,jk) * ( 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)
-               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 )
+               zv   = z0v - 0.5 * pvn(ji,jj,jk) * zstep / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
+               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) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
             END DO
          END DO
@@ -166,43 +155,20 @@
         DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-#if defined key_zco
                IF( umask(ji,jj,jk) == 0. ) THEN
                   IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
-                     zwx(ji+1,jj,jk) = e2u(ji+1,jj) * pun(ji+1,jj,jk) * ( ptn(ji+1,jj,jk) + ptn(ji+2,jj,jk) ) * 0.5
+                     zwx(ji+1,jj,jk) = pun(ji+1,jj,jk) * ( ptn(ji+1,jj,jk) + ptn(ji+2,jj,jk) ) * 0.5
                   ENDIF
                   IF( pun(ji-1,jj,jk) < 0. ) THEN
-                     zwx(ji-1,jj,jk) = e2u(ji-1,jj) * pun(ji-1,jj,jk) * ( ptn(ji-1,jj,jk) + ptn(ji  ,jj,jk) ) * 0.5
+                     zwx(ji-1,jj,jk) = pun(ji-1,jj,jk) * ( ptn(ji-1,jj,jk) + ptn(ji  ,jj,jk) ) * 0.5
                   ENDIF
                ENDIF
                IF( vmask(ji,jj,jk) == 0. ) THEN
                   IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
-                     zwy(ji,jj+1,jk) = e1v(ji,jj+1) * pvn(ji,jj+1,jk) * ( ptn(ji,jj+1,jk) + ptn(ji,jj+2,jk) ) * 0.5
+                     zwy(ji,jj+1,jk) = pvn(ji,jj+1,jk) * ( ptn(ji,jj+1,jk) + ptn(ji,jj+2,jk) ) * 0.5
                   ENDIF
                   IF( pvn(ji,jj-1,jk) < 0. ) THEN
-                     zwy(ji,jj-1,jk) = e1v(ji,jj-1) * pvn(ji,jj-1,jk) * ( ptn(ji,jj-1,jk) + ptn(ji  ,jj,jk) ) * 0.5
+                     zwy(ji,jj-1,jk) = pvn(ji,jj-1,jk) * ( ptn(ji,jj-1,jk) + ptn(ji  ,jj,jk) ) * 0.5
                   ENDIF
                ENDIF
-#else
-               IF( umask(ji,jj,jk) == 0. ) THEN
-                  IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
-                     zwx(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk)   &
-                        &            * pun(ji+1,jj,jk) * ( ptn(ji+1,jj,jk) + ptn(ji+2,jj,jk) ) * 0.5
-                  ENDIF
-                  IF( pun(ji-1,jj,jk) < 0. ) THEN
-                     zwx(ji-1,jj,jk) = e2u(ji-1,jj)* fse3u(ji-1,jj,jk)   &
-                        &            * pun(ji-1,jj,jk) * ( ptn(ji-1,jj,jk) + ptn(ji  ,jj,jk) ) * 0.5
-                  ENDIF
-               ENDIF
-               IF( vmask(ji,jj,jk) == 0. ) THEN
-                  IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
-                     zwy(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk)   &
-                        &            * pvn(ji,jj+1,jk) * ( ptn(ji,jj+1,jk) + ptn(ji,jj+2,jk) ) * 0.5
-                  ENDIF
-                  IF( pvn(ji,jj-1,jk) < 0. ) THEN
-                     zwy(ji,jj-1,jk) = e1v(ji,jj-1)* fse3v(ji,jj-1,jk)   &
-                        &            * pvn(ji,jj-1,jk) * ( ptn(ji,jj-1,jk) + ptn(ji  ,jj,jk) ) * 0.5
-                  ENDIF
-               ENDIF
-#endif
             END DO
          END DO
@@ -217,14 +183,8 @@
          DO jj = 2, jpjm1      
             DO ji = fs_2, fs_jpim1   ! vector opt.
-#if defined key_zco
-               zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
-#else
                zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
-#endif
-               ! horizontal advective trends
-               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
-               pta(ji,jj,jk) = pta(ji,jj,jk) + zta
+               ! horizontal advective trends added to the general tracer trends
+               pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                  &                                   + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) ) 
             END DO
         END DO
@@ -270,5 +230,5 @@
             DO ji = 1, jpi
                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) ) )
+                  &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) )
             END DO
          END DO
@@ -279,8 +239,7 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               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  ) ) )
+               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
@@ -293,11 +252,10 @@
          DO jj = 2, jpjm1      
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zew = pwn(ji,jj,jk+1)
                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)
-               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 )
+               zw  = z0w - 0.5 * pwn(ji,jj,jk+1)*zstep / ( e1t(ji,jj) *e2t(ji,jj) * 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) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha)*zzwy )
             END DO
          END DO
@@ -315,12 +273,9 @@
       END DO
 
-      ! surface values
-      IF( lk_dynspg_rl .OR. lk_vvl ) THEN              ! rigid lid or variable volume: flux set to zero
-         zwx(:,:, 1 ) = 0.e0                                 ! surface
-         zwx(:,:,jpk) = 0.e0                                 ! bottom 
-      ELSE                                             ! free surface
-         zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1)               ! surface
-         zwx(:,:,jpk) = 0.e0                                 ! bottom 
-      ENDIF
+      zwx(:,:,jpk) = 0.e0                                       ! bottom value
+      !                                                         ! surface values
+      IF( lk_dynspg_rl .OR. lk_vvl) THEN    ;   zwx(:,:, 1 ) = 0.e0                      ! rigid lid or variable volume
+      ELSE                                  ;   zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1)   ! linear free surface
+      ENDIF 
 
 
@@ -330,8 +285,6 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zbtr = 1. / fse3t(ji,jj,jk)
-               ! horizontal advective trends
-               zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
-               ! add it to the general tracer trends
-               pta(ji,jj,jk) =  pta(ji,jj,jk) + zta
+               ! horizontal advective trends added to the general tracer trends
+               pta(ji,jj,jk) =  pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
             END DO
          END DO
Index: /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_qck.F90
===================================================================
--- /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_qck.F90	(revision 790)
+++ /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_qck.F90	(revision 791)
@@ -5,5 +5,5 @@
    !!======================================================================
    !! History :  1.0  !  06-09  (G. Reffray)  Original code
-   !!            2.4  !  08-01  (G. Madec)  Merge TRA-TRC
+   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport
    !!----------------------------------------------------------------------
 
@@ -30,8 +30,5 @@
    PUBLIC tra_adv_qck    ! routine called by traadv.F90
 
-   REAL(wp), DIMENSION(jpi,jpj)     ::   zbtr2
    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   sl
-   REAL(wp) ::   cst1, cst2, dt, coef1                  ! temporary scalars
-   INTEGER  ::   ji, jj, jk                             ! dummy loop indices
 
    !! * Substitutions
@@ -40,5 +37,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 
-   !! $Id:$ 
+   !! $Id$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -94,5 +91,7 @@
       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend 
       !!
-      REAL(wp) :: z2                                   ! temporary scalar 
+      INTEGER  ::   ji, jj, jk        ! dummy loop indices
+      REAL(wp) ::   zdt, z2, zcoef1   ! temporary scalars
+      REAL(wp) ::   zbtr              ! temporary scalars
       !!----------------------------------------------------------------------
 
@@ -101,12 +100,8 @@
          IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3st order quickest advection scheme'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
-
-         zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
-         cst1 = 1./12.
-         cst2 = 2./3.
       ENDIF
 
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
-      ELSE                                        ;    z2 = 2.
+      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.      ! euler     time-stepping
+      ELSE                                        ;    z2 = 2.      ! leap-frog time-stepping
       ENDIF
 
@@ -116,27 +111,27 @@
       !---------------------------------------------------------------
 
+!!gm : optimisation: sl compute twice (for t & s, and even more for trc)
+!!gm   note that sl is in permanent memory be used as workspace in the vertical part !
       sl(:,:,:) = 100.
 
-                                                      ! ===============
-      DO jk = 1, jpkm1                                ! Horizontal slab
-        !                                             ! ===============
-        dt  = z2 * rdttra(jk)
-        DO jj = 2, jpjm1
-          DO ji = 2, fs_jpim1   ! vector opt.
-            coef1 = 1.e-12
-            IF (pun(ji-1,jj  ,jk  ).LT.0.) coef1 = coef1 + ABS(pun(ji-1,jj  ,jk  ))*dt/e1t(ji,jj)
-            IF (pun(ji  ,jj  ,jk  ).GT.0.) coef1 = coef1 + ABS(pun(ji  ,jj  ,jk  ))*dt/e1t(ji,jj)
-            IF (pvn(ji  ,jj-1,jk  ).LT.0.) coef1 = coef1 + ABS(pvn(ji  ,jj-1,jk  ))*dt/e2t(ji,jj)
-            IF (pvn(ji  ,jj  ,jk  ).GT.0.) coef1 = coef1 + ABS(pvn(ji  ,jj  ,jk  ))*dt/e2t(ji,jj)
-            IF (pwn(ji  ,jj  ,jk+1).LT.0.) coef1 = coef1 + ABS(pwn(ji  ,jj  ,jk+1))*dt/(fse3t(ji,jj,jk))
-            IF (pwn(ji  ,jj  ,jk  ).GT.0.) coef1 = coef1 + ABS(pwn(ji  ,jj  ,jk  ))*dt/(fse3t(ji,jj,jk))
-            sl(ji,jj,jk) = 1./coef1
-            sl(ji,jj,jk) = MIN(100.,sl(ji,jj,jk))
-            sl(ji,jj,jk) = MAX(1.  ,sl(ji,jj,jk))
-          ENDDO
-        ENDDO
-      ENDDO
-      !--- Lateral boundary conditions on the limiter slope
-      CALL lbc_lnk(   sl(:,:,:), 'T', 1. )
+      DO jk = 1, jpkm1
+         zdt  = z2 * rdttra(jk)
+         DO jj = 2, jpjm1
+            DO ji = 2, fs_jpim1   ! vector opt.
+               zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) *fse3t(ji,jj,jk) )
+               zcoef1 = 1.e-12
+               IF( pun(ji-1,jj  ,jk  ) < 0.e0 )   zcoef1 = zcoef1 + ABS( pun(ji-1,jj  ,jk  ) ) * zdt * zbtr
+               IF( pun(ji  ,jj  ,jk  ) > 0.e0 )   zcoef1 = zcoef1 + ABS( pun(ji  ,jj  ,jk  ) ) * zdt * zbtr
+               IF( pvn(ji  ,jj-1,jk  ) < 0.e0 )   zcoef1 = zcoef1 + ABS( pvn(ji  ,jj-1,jk  ) ) * zdt * zbtr
+               IF( pvn(ji  ,jj  ,jk  ) > 0.e0 )   zcoef1 = zcoef1 + ABS( pvn(ji  ,jj  ,jk  ) ) * zdt * zbtr
+               IF( pwn(ji  ,jj  ,jk+1) < 0.e0 )   zcoef1 = zcoef1 + ABS( pwn(ji  ,jj  ,jk+1) ) * zdt * zbtr
+               IF( pwn(ji  ,jj  ,jk  ) > 0.e0 )   zcoef1 = zcoef1 + ABS( pwn(ji  ,jj  ,jk  ) ) * zdt * zbtr
+               sl(ji,jj,jk) = 1. / zcoef1
+               sl(ji,jj,jk) = MIN( 100., sl(ji,jj,jk) )
+               sl(ji,jj,jk) = MAX( 1.  , sl(ji,jj,jk) )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( sl(:,:,:), 'T', 1. )      ! Lateral boundary conditions on the limiter slope
 
 
@@ -146,4 +141,5 @@
       CALL tra_adv_qck_hor( kt, cdtype, ktra, pun, pvn, ptb, pta, z2 )
 
+      !                                       ! control print
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' qck - had: ', mask1=tmask, clinfo3=cdtype )
 
@@ -152,6 +148,7 @@
       !-------------------------------------------------------------------------
 
-      CALL tra_adv_qck_ver( pwn, ptn , pta, z2 )
-
+      CALL tra_adv_qck_ver( pwn, ptn , pta )
+
+      !                                       ! control print
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' qck - zad: ', mask1=tmask, clinfo3=cdtype )
       !
@@ -159,6 +156,9 @@
 
 
-   SUBROUTINE tra_adv_qck_hor ( kt, cdtype, ktra, pun, pvn, ptb, pta, z2 )
-      !!----------------------------------------------------------------------
+   SUBROUTINE tra_adv_qck_hor ( kt, cdtype, ktra, pun, pvn, ptb, pta, pz2 )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_adv_qck_hor  ***
+      !!
+      !! ** Purpose :   
       !!
       !!----------------------------------------------------------------------
@@ -166,13 +166,14 @@
       CHARACTER(len=3), INTENT(in   )                         ::   cdtype      ! =TRA or TRC (tracer indicator)
       INTEGER         , INTENT(in   )                         ::   ktra        ! tracer index
-      REAL(wp)        , INTENT(in   )                         ::   z2          ! ???
+      REAL(wp)        , INTENT(in   )                         ::   pz2         ! =1 or 2 (euler or leap-frog)
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! horizontal effective velocity
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptb         ! before tracer field
       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta         ! tracer trend
-
-      REAL(wp) ::   za, zbtr, e1, e2, c, dir, fu, fc, fd   ! temporary scalars
-      REAL(wp) ::   coef2, coef3, fho, mask, dx
-      REAL(wp), DIMENSION(jpi,jpj) ::  zee
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmask, zlap, dwst, lim 
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zdt                                     ! temporary scalars
+      REAL(wp) ::   zbtr, zc, zdir, zfu, zfc, zfd           ! temporary scalars
+      REAL(wp) ::   zcoef1, zcoef2, zcoef3, zfho, zmsk, zdx
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmask, zlap, zdwst, zlim 
       !!----------------------------------------------------------------------
 
@@ -182,6 +183,6 @@
       zmask(:,:,:)= 0.e0                                                
       zlap (:,:,:)= 0.e0                                                  
-      dwst (:,:,:)= 0.e0                                                  
-      lim  (:,:,:)= 0.e0     
+      zdwst(:,:,:)= 0.e0                                                  
+      zlim (:,:,:)= 0.e0     
                                             
       !----------------------------------------------------------------------
@@ -189,20 +190,5 @@
       !----------------------------------------------------------------------
 
-                                                       ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         ! Initialization of metric arrays (for z- or s-coordinates)
-         ! ---------------------------------------------------------
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-#if defined key_zco
-               ! z-coordinates, no vertical scale factors
-               zee(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
-#else
-               ! vertical scale factor are used
-               zee(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
-#endif
-            END DO
-         END DO
+      DO jk = 1, jpkm1
 
          ! Laplacian of tracers (at before time step)
@@ -211,15 +197,12 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zmask(ji,jj,jk) = zee(ji,jj) * ( ptb(ji+1,jj  ,jk) - ptb(ji,jj,jk) )
-            END DO
-         END DO
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-#if defined key_zco
-               zee(ji,jj) = e1t(ji,jj) /  e2t(ji,jj)
-#else
-               zee(ji,jj) = e1t(ji,jj) / (e2t(ji,jj) * fse3t(ji,jj,jk))
-#endif
-               zlap(ji,jj,jk) = zee(ji,jj) * (  zmask(ji,jj,jk) - zmask(ji-1,jj,jk)  )
+               zmask(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * umask(ji,jj,jk)                    &
+                  &                         * ( ptb(ji+1,jj  ,jk) - ptb(ji,jj,jk) ) / e1u(ji,jj)
+            END DO
+         END DO
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zlap(ji,jj,jk) = e1t(ji,jj) * (  zmask(ji,jj,jk) - zmask(ji-1,jj,jk)  )   &
+                  &                        / ( e2t(ji,jj) * fse3t(ji,jj,jk) )
             END DO
          END DO
@@ -231,66 +214,57 @@
             DO ji = 2, fs_jpim1   ! vector opt.
                ! Upstream in the x-direction for the tracer
-               zmask(ji,jj,jk) =ptb(ji-1,jj,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji-1,jj,jk))
+               zmask(ji,jj,jk) = ptb(ji-1,jj,jk) + sl(ji,jj,jk) *( ptb(ji,jj,jk) - ptb(ji-1,jj,jk) )
                ! Downstream in the x-direction for the tracer
-               dwst (ji,jj,jk) =ptb(ji+1,jj,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji+1,jj,jk))
-            ENDDO
-         ENDDO
-      END DO
-      !--- Lateral boundary conditions on the laplacian (unchanged sgn)
-      CALL lbc_lnk(   zlap(:,:,:), 'T', 1. ) 
-      !--- Lateral boundary conditions for the lim function
-      CALL lbc_lnk(  zmask(:,:,:), 'T', 1. )   ;   CALL lbc_lnk(  dwst(:,:,:), 'T', 1. )
-                                                       ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         ! --- lim at the U-points in function of the direction of the flow
-         ! ----------------------------------------------------------------
+               zdwst(ji,jj,jk) = ptb(ji+1,jj,jk) + sl(ji,jj,jk) * ( ptb(ji,jj,jk) - ptb(ji+1,jj,jk) )
+            END DO
+         END DO
+      END DO
+
+      !--- Lateral boundary conditions on the laplacian and lim functions (unchanged sgn)
+      CALL lbc_lnk( zlap (:,:,:), 'T', 1. ) 
+      CALL lbc_lnk( zmask(:,:,:), 'T', 1. )   ;   CALL lbc_lnk(  zdwst(:,:,:), 'T', 1. )
+            
+      ! --- lim at the U-points in function of the direction of the flow
+      ! ----------------------------------------------------------------
+      DO jk = 1, jpkm1
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               dir = 0.5 + sign(0.5,pun(ji,jj,jk))     ! if pun>0 : diru = 1 otherwise diru = 0
-               lim(ji,jj,jk)=dir*zmask(ji,jj,jk)+(1-dir)*dwst(ji+1,jj,jk)
+               zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )     ! if pun>0 : diru = 1 otherwise diru = 0
+               zlim(ji,jj,jk) = zdir * zmask(ji,jj,jk) + (1-zdir) * zdwst(ji+1,jj,jk)
                ! Mask at the T-points in the x-direction (mask=0 or mask=1)
-               zmask(ji,jj,jk)=tmask(ji-1,jj,jk)+tmask(ji,jj,jk)+tmask(ji+1,jj,jk)-2
-            END DO
-         END DO
-      END DO
-      !--- Lateral boundary conditions for the mask
-      CALL lbc_lnk(  zmask(:,:,:), 'T', 1. )  
+               zmask(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk(  zmask(:,:,:), 'T', 1. )       !--- Lateral boundary conditions for the mask
+
 
       ! Horizontal advective fluxes
       ! ---------------------------
-                                                       ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-                                                       ! ===============
-         dt  = z2 * rdttra(jk)
+      DO jk = 1, jpkm1
+         zdt  = pz2 * rdttra(jk)
          !--- tracer flux at u and v-points
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt. 
-#if defined key_zco
-               e2 = e2u(ji,jj)
-#else
-               e2 = e2u(ji,jj) * fse3u(ji,jj,jk)
-#endif
-               dir = 0.5 + sign(0.5,pun(ji,jj,jk))                ! if pun>0 : diru = 1 otherwise diru = 0
-
-               dx = dir * e1t(ji,jj) + (1-dir)* e1t(ji+1,jj)
-               c  = ABS(pun(ji,jj,jk))*dt/dx                      ! (0<cx<1 : Courant number on x-direction)
-
-               fu  = lim(ji,jj,jk)                                ! FU + sl(FC-FU) in the x-direction for T
-               fc  = dir*ptb(ji  ,jj,jk)+(1-dir)*ptb(ji+1,jj,jk)  ! FC in the x-direction for T
-               fd  = dir*ptb(ji+1,jj,jk)+(1-dir)*ptb(ji  ,jj,jk)  ! FD in the x-direction for T
+               zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )                     ! if pun>0 : diru = 1 otherwise diru = 0
+               zdx  = ( zdir * e1t(ji,jj) + (1-zdir)* e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk)
+               zc   = ABS( pun(ji,jj,jk) ) * zdt / zdx                     ! (0<cx<1 : Courant number on x-direction)
+
+               zfu  = zlim(ji,jj,jk)                                       ! FU + sl(FC-FU) in the x-direction for T
+               zfc  = zdir * ptb(ji  ,jj,jk) + (1-zdir) * ptb(ji+1,jj,jk)  ! FC in the x-direction for T
+               zfd  = zdir * ptb(ji+1,jj,jk) + (1-zdir) * ptb(ji  ,jj,jk)  ! FD in the x-direction for T
 
                !--- QUICKEST scheme
                ! Temperature on the x-direction
-               coef1 = 0.5*(fc+fd)
-               coef2 = 0.5*c*(fd-fc)
-               coef3 = ((1.-(c*c))/6.)*(dir*zlap(ji,jj,jk) + (1-dir)*zlap(ji+1,jj,jk) )
-               fho = coef1-coef2-coef3
-               fho = bound(fu,fd,fc,fho)
+               zcoef1 = 0.5      * ( zfc + zfd )
+               zcoef2 = 0.5 * zc * ( zfd - zfc )
+               zcoef3 = ( (1.-(zc*zc)) / 6. ) * ( zdir * zlap(ji,jj,jk) + (1-zdir) * zlap(ji+1,jj,jk) )
+               zfho = zcoef1 - zcoef2 - zcoef3
+               zfho = bound( zfu, zfd, zfc, zfho )
                !--- If the second ustream point is a land point
                !--- the flux is computed by the 1st order UPWIND scheme
-               mask=dir*zmask(ji,jj,jk)+(1-dir)*zmask(ji+1,jj,jk)
-               fho = mask*fho + (1-mask)*fc
-               dwst(ji,jj,jk)=e2*pun(ji,jj,jk)*fho
+               zmsk = zdir * zmask(ji,jj,jk) + (1-zdir) * zmask(ji+1,jj,jk)
+               zfho = zmsk * zfho            + (1-zmsk) * zfc
+               zdwst(ji,jj,jk) = pun(ji,jj,jk) * zfho
             END DO
          END DO
@@ -300,20 +274,13 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                !--- horizontal advective trends
-#if defined key_zco
-               zbtr = zbtr2(ji,jj)
-#else
-               zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)              
-#endif
-               za = - zbtr * ( dwst(ji,jj,jk) - dwst(ji-1,jj  ,jk) )
+               zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
                !--- add it to the general tracer trends
-               pta(ji,jj,jk) = pta(ji,jj,jk) + za
-            END DO
-         END DO
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-
-      ! Save the horizontal advective trends for diagnostic
-      IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, dwst, pun, ptb )
+               pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zdwst(ji,jj,jk) - zdwst(ji-1,jj  ,jk) )
+            END DO
+         END DO
+      END DO
+
+      !                                            !-- trend diagnostic
+      IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zdwst, pun, ptb )
 
 
@@ -321,20 +288,6 @@
       ! I. Part 2 : y-direction
       !----------------------------------------------------------------------
-                                                       ! ==============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         ! Initialization of metric arrays (for z- or s-coordinates)
-         ! ---------------------------------------------------------
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-#if defined key_zco
-               ! z-coordinates, no vertical scale factors
-               zee(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
-#else
-               ! s-coordinates, vertical scale factor are used
-               zee(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
-#endif
-            END DO
-         END DO
+
+      DO jk = 1, jpkm1
 
          ! Laplacian of tracers (at before time step)
@@ -343,5 +296,6 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zmask(ji,jj,jk) = zee(ji,jj) * ( ptb(ji  ,jj+1,jk) - ptb(ji,jj,jk) )
+               zmask(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk)   &
+                  &            * ( ptb(ji  ,jj+1,jk) - ptb(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk)
             END DO
          END DO
@@ -349,10 +303,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-#if defined key_zco
-               zee(ji,jj) = e2t(ji,jj) /  e1t(ji,jj)
-#else
-               zee(ji,jj) = e2t(ji,jj) / (e1t(ji,jj) * fse3t(ji,jj,jk))
-#endif
-               zlap(ji,jj,jk) = zee(ji,jj) * (  zmask(ji,jj,jk) - zmask(ji,jj-1,jk)  )
+               zlap(ji,jj,jk) = e2t(ji,jj) * (  zmask(ji,jj,jk) - zmask(ji,jj-1,jk)  )   &
+                  &                        / ( e1t(ji,jj) * fse3t(ji,jj,jk) )
             END DO
          END DO
@@ -362,66 +312,57 @@
             DO ji = 2, fs_jpim1   ! vector opt.
                ! Upstream in the y-direction for the tracer
-               zmask(ji,jj,jk)=ptb(ji,jj-1,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji,jj-1,jk))
+               zmask(ji,jj,jk) = ptb(ji,jj-1,jk) + sl(ji,jj,jk) *( ptb(ji,jj,jk) - ptb(ji,jj-1,jk) )
                ! Downstream in the y-direction for the tracer
-               dwst (ji,jj,jk)=ptb(ji,jj+1,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji,jj+1,jk))
+               zdwst(ji,jj,jk) = ptb(ji,jj+1,jk) + sl(ji,jj,jk) *( ptb(ji,jj,jk) - ptb(ji,jj+1,jk) )
             ENDDO
          ENDDO
       END DO
-      !--- Lateral boundary conditions on the laplacian (unchanged sgn)
-      CALL lbc_lnk(   zlap(:,:,:), 'T', 1. )
-      !--- Lateral boundary conditions for the lim function
-      CALL lbc_lnk(  zmask(:,:,:), 'T', 1. )   ;   CALL lbc_lnk(  dwst(:,:,:), 'T', 1. )
-
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
+      !--- Lateral boundary conditions on the laplacian and lim function (unchanged sgn)
+      CALL lbc_lnk( zlap (:,:,:), 'T', 1. )
+      CALL lbc_lnk( zmask(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zdwst(:,:,:), 'T', 1. )
+
+      DO jk = 1, jpkm1
+      
          ! --- lim at the V-points in function of the direction of the flow
          ! ----------------------------------------------------------------
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               dir = 0.5 + sign(0.5,pvn(ji,jj,jk))      ! if pvn>0 : dirv = 1 otherwise dirv = 0
-               lim(ji,jj,jk)=dir*zmask(ji,jj,jk)+(1-dir)*dwst(ji,jj+1,jk)
+               zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )      ! if pvn>0 : dirv = 1 otherwise dirv = 0
+               zlim(ji,jj,jk) = zdir * zmask(ji,jj,jk) + (1-zdir) * zdwst(ji,jj+1,jk)
                ! Mask at the T-points in the y-direction (mask=0 or mask=1)
-               zmask(ji,jj,jk)=tmask(ji,jj-1,jk)+tmask(ji,jj,jk)+tmask(ji,jj+1,jk)-2
-            END DO
-         END DO
-      END DO
-      !--- Lateral boundary conditions for the mask
-      CALL lbc_lnk(  zmask(:,:,:), 'T', 1. ) 
+               zmask(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk(  zmask(:,:,:), 'T', 1. )      !--- Lateral boundary conditions for the mask
 
       ! Horizontal advective fluxes
-      ! -------------------------------
-                                                       ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-                                                       ! ===============
-         dt  = z2 * rdttra(jk)
+      ! ---------------------------
+                                  
+      DO jk = 1, jpkm1           
+         zdt  = pz2 * rdttra(jk)
          !--- tracer flux at u and v-points
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-#if defined key_zco
-               e1 = e1v(ji,jj)
-#else
-               e1 = e1v(ji,jj) * fse3v(ji,jj,jk)
-#endif
-               dir = 0.5 + sign(0.5,pvn(ji,jj,jk))                ! if pvn>0 : dirv = 1 otherwise dirv = 0
-
-               dx = dir * e2t(ji,jj) + (1-dir)* e2t(ji,jj+1)
-               c  = ABS(pvn(ji,jj,jk))*dt/dx                      ! (0<cy<1 : Courant number on y-direction)
-
-               fu  = lim(ji,jj,jk)                                ! FU + sl(FC-FU) in the y-direction for T
-               fc  = dir*ptb(ji,jj  ,jk)+(1-dir)*ptb(ji,jj+1,jk)  ! FC in the y-direction for T
-               fd  = dir*ptb(ji,jj+1,jk)+(1-dir)*ptb(ji,jj  ,jk)  ! FD in the y-direction for T
+               zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )                     ! if pvn>0 : dirv = 1 otherwise dirv = 0
+               zdx  = ( zdir * e2t(ji,jj) + (1-zdir)* e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)
+               zc   = ABS( pvn(ji,jj,jk) ) * zdt / zdx                     ! (0<cy<1 : Courant number on y-direction)
+
+               zfu  = zlim(ji,jj,jk)                                       ! FU + sl(FC-FU) in the y-direction for T
+               zfc  = zdir * ptb(ji,jj  ,jk) + (1-zdir) * ptb(ji,jj+1,jk)  ! FC in the y-direction for T
+               zfd  = zdir * ptb(ji,jj+1,jk) + (1-zdir) * ptb(ji,jj  ,jk)  ! FD in the y-direction for T
 
                !--- QUICKEST scheme
                ! Temperature on the y-direction
-               coef1 = 0.5*(fc+fd)
-               coef2 = 0.5*c*(fd-fc)
-               coef3 = ((1.-(c*c))/6.)*(dir*zlap(ji,jj,jk) + (1-dir)*zlap(ji,jj+1,jk) )
-               fho = coef1-coef2-coef3
-               fho = bound(fu,fd,fc,fho)
+               zcoef1 = 0.5      * ( zfc + zfd )
+               zcoef2 = 0.5 * zc * ( zfd - zfc )
+               zcoef3 = ( (1.-(zc*zc)) / 6. ) * ( zdir * zlap(ji,jj,jk) + (1-zdir) * zlap(ji,jj+1,jk) )
+               zfho = zcoef1 - zcoef2 - zcoef3
+               zfho = bound( zfu, zfd, zfc, zfho )
                !--- If the second ustream point is a land point
                !--- the flux is computed by the 1st order UPWIND scheme
-               mask=dir*zmask(ji,jj,jk)+(1-dir)*zmask(ji,jj+1,jk)
-               fho = mask*fho + (1-mask)*fc
-               dwst(ji,jj,jk)=e1*pvn(ji,jj,jk)*fho
+               zmsk = zdir * zmask(ji,jj,jk) + (1-zdir) * zmask(ji,jj+1,jk)
+               zfho = zmsk * zfho            + (1-zmsk) * zfc
+               zdwst(ji,jj,jk) = pvn(ji,jj,jk) * zfho
             END DO
          END DO
@@ -430,35 +371,17 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               !--- horizontal advective trends
-#if defined key_zco
-               zbtr = zbtr2(ji,jj)
-#else
-               zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
-#endif
-               za = - zbtr * (  dwst(ji,jj,jk) - dwst(ji  ,jj-1,jk) )
-               !--- add it to the general tracer trends
-               pta(ji,jj,jk) = pta(ji,jj,jk) + za
-            END DO
-         END DO
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-
-      ! Save the horizontal advective trends for diagnostic
-      IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, dwst, pvn, ptb )
-
-      ! "zonal" mean advective heat and salt transport 
+               zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+               ! horizontal advective trends added to the general tracer trends
+               pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * (  zdwst(ji,jj,jk) - zdwst(ji,jj-1,jk) )
+            END DO
+         END DO
+      END DO 
+
+      !                                           !-- trend diagnostic
+      IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zdwst, pvn, ptb )
+      !                                           ! "Poleward" heat or salt transport 
       IF( ln_diaptr .AND. cdtype == 'TRA' .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
-#if defined key_zco
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  dwst(ji,jj,jk) = dwst(ji,jj,jk) * fse3v(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-# endif
-         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( dwst(:,:,:) )
-         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( dwst(:,:,:) )
+         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zdwst(:,:,:) )
+         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zdwst(:,:,:) )
       ENDIF
       !
@@ -466,14 +389,18 @@
 
 
-   SUBROUTINE tra_adv_qck_ver( pwn, ptn , pta, z2  )      
-      !!----------------------------------------------------------------------
-      !!
-      !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in   )                         ::   z2
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwn
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptn
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta
-      !!
-      REAL(wp) ::   za, ze3tr, dt, dir, fc, fd             ! temporary scalars
+   SUBROUTINE tra_adv_qck_ver( pwn, ptn , pta )      
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_adv_qck_ver  ***
+      !!
+      !! ** Purpose :   
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwn   ! vertical transport
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptn   ! now tracer
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta   ! tracer trend
+      !!
+      INTEGER  ::   ji, jj , jk                  ! dummy loop indices
+      REAL(wp) ::   zbtr, zdir, zfc, zfd   ! temporary scalars
+      !!----------------------------------------------------------------------
 
       ! Vertical advection
@@ -483,23 +410,23 @@
       ! ----------------------------
 
-      sl(:,:,jpk) = 0.e0      !Bottom value : flux set to zero
-
-      ! Surface value
-      IF( lk_dynspg_rl .OR. lk_vvl ) THEN         ! rigid lid : flux set to zero
-         sl(:,:, 1 ) = 0.e0
-      ELSE                                        ! free surface-constant volume
-         sl(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1)
+      sl(:,:,jpk) = 0.e0                          ! bottom value 
+
+      !                                           ! Surface value
+      IF( lk_dynspg_rl .OR. lk_vvl ) THEN   ;   sl(:,:, 1 ) = 0.e0                      ! rigid lid or non-linear fs
+      ELSE                                  ;   sl(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1)   ! linear free surface
       ENDIF
+
+!!gm  bug: code au true 2nd order scheme
+!!gm  : sl used as work array: not good idea for optimisation (sl compute once for all tracers...)
 
       ! Second order centered tracer flux at w-point
       DO jk = 2, jpkm1
-         dt  = z2 *  rdttra(jk)
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               dir = 0.5 + sign(0.5,pwn(ji,jj,jk))                                         ! if pwn>0 : dirw = 1 otherwise dirw = 0
-               fc = dir*ptn(ji,jj,jk  )*fse3t(ji,jj,jk-1)+(1-dir)*ptn(ji,jj,jk-1)*fse3t(ji,jj,jk  )   ! FC in the z-direction for T
-               fd = dir*ptn(ji,jj,jk-1)*fse3t(ji,jj,jk  )+(1-dir)*ptn(ji,jj,jk  )*fse3t(ji,jj,jk-1)   ! FD in the z-direction for T
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zdir = 0.5 + SIGN( 0.5, pwn(ji,jj,jk) )                        ! if pwn>0 : dirw = 1 otherwise dirw = 0
+               zfc = zdir * ptn(ji,jj,jk  ) * fse3t(ji,jj,jk-1) + (1-zdir) * ptn(ji,jj,jk-1) * fse3t(ji,jj,jk  )   ! FC 
+               zfd = zdir * ptn(ji,jj,jk-1) * fse3t(ji,jj,jk  ) + (1-zdir) * ptn(ji,jj,jk  ) * fse3t(ji,jj,jk-1)   ! FD
                !--- Second order centered scheme
-               sl(ji,jj,jk)=pwn(ji,jj,jk)*(fc+fd)/(fse3t(ji,jj,jk-1)+fse3t(ji,jj,jk))
+               sl(ji,jj,jk) = pwn(ji,jj,jk) * ( zfc + zfd ) / ( fse3t(ji,jj,jk-1) + fse3t(ji,jj,jk) )
             END DO
          END DO
@@ -511,9 +438,7 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ze3tr = 1. / fse3t(ji,jj,jk)
-               ! vertical advective trends
-               za = - ze3tr * ( sl(ji,jj,jk) - sl(ji,jj,jk+1) )
-               ! add it to the general tracer trends
-               pta(ji,jj,jk) =  pta(ji,jj,jk) + za
+               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+               ! vertical advective trends added to the general tracer trends
+               pta(ji,jj,jk) =  pta(ji,jj,jk) - zbtr * ( sl(ji,jj,jk) - sl(ji,jj,jk+1) )
             END DO
          END DO
@@ -523,9 +448,19 @@
 
 
-   REAL FUNCTION bound(fu,fd,fc,fho)
-      REAL(wp) ::  fu, fd, fc, fho, fref1, fref2
-      fref1 = fu
-      fref2 = MAX(  MIN( fc , fd ), MIN( MAX( fc , fd ), fref1 )  )
-      bound = MAX(  MIN( fho, fc ), MIN( MAX( fho, fc ), fref2 )  )
+   FUNCTION bound( pfu, pfd, pfc, pfho )   RESULT( pbound )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION bound  ***
+      !!
+      !! ** Purpose :   ???
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in) ::  pfu, pfd, pfc, pfho
+      REAL(wp)             ::  pbound
+      !!
+      REAL(wp) ::  zfref1, zfref2
+      !!----------------------------------------------------------------------
+      zfref1 = pfu
+      zfref2 = MAX(  MIN( pfc , pfd ), MIN( MAX( pfc , pfd ), zfref1 )  )
+      pbound = MAX(  MIN( pfho, pfc ), MIN( MAX( pfho, pfc ), zfref2 )  )
+      !
    END FUNCTION
 
Index: /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
===================================================================
--- /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 790)
+++ /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 791)
@@ -4,14 +4,14 @@
    !! Ocean active tracers:  horizontal & vertical advective trend
    !!==============================================================================
-   !! History :  7.0  !  95-12  (L. Mortier)  Original code
-   !!            8.0  !  00-01  (H. Loukos)  adapted to ORCA 
-   !!             -   !  00-10  (MA Foujols E.Kestenare)  include file not routine
-   !!             -   !  00-12  (E. Kestenare M. Levy)  fix bug in trtrd indexes
-   !!             -   !  01-07  (E. Durand G. Madec)  adaptation to ORCA config
-   !!            8.5  !  02-06  (G. Madec)  F90: Free form and module
-   !!   NEMO     1.0  !  04-01  (A. de Miranda, G. Madec, J.M. Molines ): advective bbl
-   !!             -   !  08-04  (S. Cravatte) add the i-, j- & k- trends computation
-   !!             -   !  05-11  (V. Garnier) Surface pressure gradient organization
-   !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC
+   !! History :  7.0  !  1995-12  (L. Mortier)  Original code
+   !!            8.0  !  2000-01  (H. Loukos)  adapted to ORCA 
+   !!             -   !  2000-10  (MA Foujols E.Kestenare)  include file not routine
+   !!             -   !  2000-12  (E. Kestenare M. Levy)  fix bug in trtrd indexes
+   !!             -   !  2001-07  (E. Durand G. Madec)  adaptation to ORCA config
+   !!            8.5  !  2002-06  (G. Madec)  F90: Free form and module
+   !!   NEMO     1.0  !  2004-01  (A. de Miranda, G. Madec, J.M. Molines ): advective bbl
+   !!             -   !  2008-04  (S. Cravatte) add the i-, j- & k- trends computation
+   !!             -   !  2005-11  (V. Garnier) Surface pressure gradient organization
+   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport
    !!----------------------------------------------------------------------
 
@@ -42,5 +42,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 
-   !! $Id:$ 
+   !! $Id$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -70,11 +70,9 @@
       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend 
       !!
-      INTEGER  ::   ji, jj, jk              ! dummy loop indices
-      REAL(wp) ::   ztai, ztaj, ztak
-      REAL(wp) ::   z2dtt, zbtr, zeu, zev   ! temporary scalar
-      REAL(wp) ::   zew, z2                          ! temporary scalar
-      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk           !    "         "
-      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk            !    "         "
-      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zti, ztu, ztv, ztw   ! temporary workspace
+      INTEGER  ::   ji, jj, jk               ! dummy loop indices
+      REAL(wp) ::   z2dtt, zbtr, z2, zzti    ! temporary scalar
+      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !    "         "
+      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !    "         "
+      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zti, ztu, ztv, ztw   ! 3D workspace
       !!----------------------------------------------------------------------
 
@@ -87,6 +85,6 @@
       ENDIF
 
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
-      ELSE                                        ;    z2 = 2.
+      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.      ! euler   time-stepping
+      ELSE                                        ;    z2 = 2.      ! leap-frog time-stepping
       ENDIF
 
@@ -103,54 +101,51 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
-               zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
-               ! upstream scheme
-               zfp_ui = zeu + ABS( zeu )
-               zfm_ui = zeu - ABS( zeu )
-               zfp_vj = zev + ABS( zev )
-               zfm_vj = zev - ABS( zev )
-               ztu(ji,jj,jk) = zfp_ui * ptb(ji,jj,jk) + zfm_ui * ptb(ji+1,jj  ,jk)
-               ztv(ji,jj,jk) = zfp_vj * ptb(ji,jj,jk) + zfm_vj * ptb(ji  ,jj+1,jk)
-            END DO
-         END DO
-      END DO
-
+               zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
+               zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
+               zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
+               zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
+               ztu(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk) + zfm_ui * ptb(ji+1,jj  ,jk) )
+               ztv(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk) + zfm_vj * ptb(ji  ,jj+1,jk) )
+            END DO
+         END DO
+      END DO
+      !
       ! upstream tracer flux in the k direction
-      ! Surface value
-      IF( lk_dynspg_rl .OR. lk_vvl ) THEN               ! rigid lid or variable volume: flux set to zero
-         ztw(:,:,1) = 0.e0
-      ELSE                                              ! free surface
-         ztw(:,:,1) = e1t(:,:) * e2t(:,:) * pwn(:,:,1) * ptb(:,:,1)
-      ENDIF
-
-      ! Interior value
-      DO jk = 2, jpkm1
+      !                                                                           ! Surface value
+      IF( lk_dynspg_rl .OR. lk_vvl ) THEN   ;   ztw(:,:,1) = 0.e0                      ! rigid lid or non-linear fs
+      ELSE                                  ;   ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1)   ! free surface
+      ENDIF
+      !
+      DO jk = 2, jpkm1                                                            ! Interior value
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
-               zfp_wk = zew + ABS( zew )
-               zfm_wk = zew - ABS( zew )
-               ztw(ji,jj,jk) = zfp_wk * ptb(ji,jj,jk) + zfm_wk * ptb(ji,jj,jk-1)
-            END DO
-         END DO
-      END DO
-
-      ! total advective trend
-      DO jk = 1, jpkm1
+               zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
+               zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
+               ztw(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk) + zfm_wk * ptb(ji,jj,jk-1) )
+            END DO
+         END DO
+      END DO
+      !
+      ! upstream advective trend
+      DO jk = 1, jpkm1
+         z2dtt = z2 * rdttra(jk)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-               ! i- j- horizontal & k- vertical advective trends
-               ztai = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  ) ) * zbtr
-               ztaj = - ( ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  ) ) * zbtr
-               ztak = - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) * zbtr
                ! total intermediate advective trends
-               zti(ji,jj,jk) = ztai + ztaj + ztak
-            END DO
-         END DO
-      END DO
-
-      !  Save the horizontal advective trends for diagnostic
-      ! -----------------------------------------------------
+               zzti = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )            &
+                  &     + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )            &
+                  &     + ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) * zbtr
+               ! update and guess with monotonic sheme
+               pta(ji,jj,jk) =   pta(ji,jj,jk) +         zzti
+               zti(ji,jj,jk) = ( ptb(ji,jj,jk) + z2dtt * zzti ) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti   (unchanged sign)
+
+
+      !                                 ! trend diagnostics (contribution of upstream fluxes)
       IF( l_trdtra ) THEN
          CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, ztu, pun, ptn )
@@ -158,90 +153,66 @@
          CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, ztw, pwn, ptn )
       ENDIF
-
-      ! update and guess with monotonic sheme
-      DO jk = 1, jpkm1
-         z2dtt = z2 * rdttra(jk)
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               pta(ji,jj,jk) =  pta(ji,jj,jk) + zti(ji,jj,jk)
-               zti (ji,jj,jk) = ( ptb(ji,jj,jk) + z2dtt * zti(ji,jj,jk) ) * tmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-
-      ! Lateral boundary conditions on zti   (unchanged sign)
-      CALL lbc_lnk( zti, 'T', 1. )
+      !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
+      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
+         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( ztv(:,:,:) ) 
+         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( ztv(:,:,:) )
+      ENDIF
 
 
       ! 3. antidiffusive flux : high order minus low order
       ! --------------------------------------------------
-      ! antidiffusive flux on i and j
+      !                                 ! anti-diffusive flux on i and j
       DO jk = 1, jpkm1
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
-               zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
-               ztu(ji,jj,jk) = zeu * ( ptn(ji,jj,jk) + ptn(ji+1,jj,jk) ) - ztu(ji,jj,jk)
-               ztv(ji,jj,jk) = zev * ( ptn(ji,jj,jk) + ptn(ji,jj+1,jk) ) - ztv(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      
-      ! antidiffusive flux on k
-      ! Surface value
-      ztw(:,:,1) = 0.e0
-
-      ! Interior value
-      DO jk = 2, jpkm1
+               ztu(ji,jj,jk) = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji+1,jj,jk) ) - ztu(ji,jj,jk)
+               ztv(ji,jj,jk) = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji,jj+1,jk) ) - ztv(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !                                 ! antidiffusive flux on k
+      ztw(:,:,1) = 0.e0                      ! Surface value
+      !
+      DO jk = 2, jpkm1                       ! Interior value
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
-               ztw(ji,jj,jk) = zew * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) - ztw(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-
-      ! Lateral bondary conditions
-      CALL lbc_lnk( ztu, 'U', -1. )
+               ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) - ztw(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      CALL lbc_lnk( ztu, 'U', -1. )     ! Lateral bondary conditions on upstream fluxes
       CALL lbc_lnk( ztv, 'V', -1. )
       CALL lbc_lnk( ztw, 'W',  1. )
 
-      ! 4. monotonicity algorithm
-      ! -------------------------
+      !                                 ! monotonicity algorithm
       CALL nonosc( ptb, ztu, ztv, ztw, zti, z2 )
 
 
-      ! 5. final trend with corrected fluxes
-      ! ------------------------------------
+      ! 4. final trend with anti-diffusive fluxes
+      ! -----------------------------------------
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.  
                zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-               ! i- j- horizontal & k- vertical advective trends
-               ztai = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )) * zbtr
-               ztaj = - ( ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )) * zbtr
-               ztak = - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1)) * zbtr
-
-               ! add them to the general tracer trends
-               pta(ji,jj,jk) = pta(ji,jj,jk) + ztai + ztaj + ztak
-            END DO
-         END DO
-      END DO
-
-!!gm  the transport computation is wrong, the upstream part is missing !
-      ! "zonal" mean advective heat and salt transport
-      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
-         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( ztv(:,:,:) )
-         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( ztv(:,:,:) )
-      ENDIF
-
-      !  Save the horizontal advective trends for diagnostic
-      ! -----------------------------------------------------
-      IF( l_trdtra ) THEN
+               ! anti-diffusive trends added to the general tracer trends
+               pta(ji,jj,jk) = pta(ji,jj,jk) - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )             &
+                  &                            + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )             &
+                  &                            + ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1)  ) * zbtr
+            END DO
+         END DO
+      END DO
+
+      IF( l_trdtra ) THEN               ! trend diagnostic (contribution of anti-diffusive fluxes)
          CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, ztu, pun, ptn, cnbpas='bis' )   ! <<< Add to iad trend
          CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, ztv, pvn, ptn, cnbpas='bis' )   ! <<< Add to jad trend
          CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, ztw, pwn, ptn, cnbpas='bis' )   ! <<< Add to zad trend
       ENDIF
-
+      !                                 ! "Poleward" heat and salt transports (contribution of anti-diffusive fluxes)
+      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
+         IF( ktra == jp_tem)   pht_adv(:) = pht_adv(:) + ptr_vj( ztv(:,:,:) )
+         IF( ktra == jp_sal)   pst_adv(:) = pst_adv(:) + ptr_vj( ztv(:,:,:) )
+      ENDIF
+      !                                 ! control print
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' tvd - adv: ', mask1=tmask, clinfo3=cdtype )
       !
@@ -276,4 +247,5 @@
       zbetup(:,:,:) = 0.e0   ;   zbetdo(:,:,:) = 0.e0  
 
+!!gm   optimisation :  add the optimal version I wrote 1 year ago
       ! Search local extrema
       ! --------------------
Index: /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
===================================================================
--- /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_ubs.F90	(revision 790)
+++ /branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_ubs.F90	(revision 791)
@@ -4,6 +4,6 @@
    !! Ocean tracers:  horizontal & vertical advective trend
    !!==============================================================================
-   !! History :  1.0  !  06-08  (L. Debreu, R. Benshila)  Original code
-   !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC
+   !! History :  1.0  !  2006-08  (L. Debreu, R. Benshila)  Original code
+   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport
    !!----------------------------------------------------------------------
 
@@ -27,6 +27,4 @@
    PUBLIC   tra_adv_ubs   ! routine called by traadv module
 
-   REAL(wp), DIMENSION(jpi,jpj) ::   e1e2tr   ! = 1/(e1t * e2t)
-
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -34,5 +32,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA & TRP 2.4 , LOCEAN-IPSL (2008) 
-   !! $Id:$ 
+   !! $Id$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    !!----------------------------------------------------------------------
@@ -71,5 +69,5 @@
       !!
       !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404. 
-      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731–1741. 
+      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731-1741. 
       !!----------------------------------------------------------------------
       INTEGER         , INTENT(in   )                         ::   kt              ! ocean time-step index
@@ -80,14 +78,13 @@
       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend 
       !!
-      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
-      REAL(wp) ::   zta, zbtr, zcoef                  ! temporary scalars
-      REAL(wp) ::   zfui, zfp_ui, zfm_ui, zcenut      !    "         "
-      REAL(wp) ::   zfvj, zfp_vj, zfm_vj, zcenvt      !    "         "
-      REAL(wp) ::   z_hdivn                           !    "         "
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zeeu, zeev     ! temporary 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zwy              ! temporary 3D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztu , ztv , zltu , zltv, ztrdt   !    "              "
-      !!----------------------------------------------------------------------
-
+      INTEGER  ::   ji, jj, jk                     ! dummy loop indices
+      REAL(wp) ::   zbtr, zcoef, z_hdivn           ! temporary scalars
+      REAL(wp) ::   zeeu, zfp_ui, zfm_ui, zcenut   !    "         "
+      REAL(wp) ::   zeev, zfp_vj, zfm_vj, zcenvt   !    "         "
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, ztu, zltu   ! temporary 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, ztv, zltv   !    "              "
+      !!----------------------------------------------------------------------
+
+!!gm  this can be optimized: we don't need zero everywhere
       zltu(:,:,:) = 0.e0
       zltv(:,:,:) = 0.e0
@@ -97,128 +94,75 @@
          IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
-         !
-         e1e2tr(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
-      ENDIF
-
-      ! store pta trends
-      ztrdt(:,:,:) = pta(:,:,:)
-
-      zcoef = 1./6.
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-
-         !  Initialization of metric arrays (for z- or s-coordinates)
+      ENDIF
+
+      !  Laplacian
+      DO jk = 1, jpkm1
+         DO jj = 1, jpjm1                               ! First derivative (gradient)
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zeeu = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
+               zeev = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
+               ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj  ,jk) - ptb(ji,jj,jk) )
+               ztv(ji,jj,jk) = zeev * ( ptb(ji  ,jj+1,jk) - ptb(ji,jj,jk) )
+            END DO
+         END DO
+         DO jj = 2, jpjm1                               ! Second derivative (divergence)
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zcoef = 1. / ( 6. * fse3t(ji,jj,jk) )
+               zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
+               zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
+            END DO
+         END DO
+      END DO
+
+      ! Lateral boundary conditions on the laplacian (zlt,zls)   (unchanged sgn)
+      CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )
+
+
+      !  Horizontal advective fluxes
+      DO jk = 1, jpkm1
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-#if defined key_zco
-               ! z-coordinates, no vertical scale factors
-               zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
-               zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
-#else
-               ! s-coordinates, vertical scale factor are used
-               zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
-               zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
-#endif
-            END DO
-         END DO
-
-         !  Laplacian
-         ! First derivative (gradient)
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               ztu(ji,jj,jk) = zeeu(ji,jj) * ( ptb(ji+1,jj  ,jk) - ptb(ji,jj,jk) )
-               ztv(ji,jj,jk) = zeev(ji,jj) * ( ptb(ji  ,jj+1,jk) - ptb(ji,jj,jk) )
-            END DO
-         END DO
-         ! Second derivative (divergence)
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-#if ! defined key_zco
-               zcoef = 1. / ( 6. * fse3t(ji,jj,jk) )
-#endif         
-               zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
-               zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
-            END DO
-         END DO
-         !                                             ! =================
-      END DO                                           !    End of slab
-      !                                                ! =================
-
-      ! Lateral boundary conditions on the laplacian (zlt,zls)   (unchanged sgn)
-      CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )
-
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         !  Horizontal advective fluxes
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               ! volume fluxes * 1/2
-#if defined key_zco
-               zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk)
-               zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk)
-#else
-               zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
-               zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
-#endif
-               ! upstream scheme
-               zfp_ui = zfui + ABS( zfui )
-               zfp_vj = zfvj + ABS( zfvj )
-               zfm_ui = zfui - ABS( zfui )
-               zfm_vj = zfvj - ABS( zfvj )
+               ! upstream transport
+               zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )      ! = 2.pun  if pun>0,  =0      if pun<0
+               zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )      ! = 0      if pun>0,  = 2.pun if pun<0
+               zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )      ! idem on pvn
+               zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )      !
                ! centered scheme
-               zcenut = zfui * ( ptn(ji,jj,jk) + ptn(ji+1,jj  ,jk) )
-               zcenvt = zfvj * ( ptn(ji,jj,jk) + ptn(ji  ,jj+1,jk) )
-               ! mixed centered / upstream scheme
-               zwx(ji,jj,jk) = zcenut - zfp_ui * zltu(ji,jj,jk) -zfm_ui * zltu(ji+1,jj,jk)
-               zwy(ji,jj,jk) = zcenvt - zfp_vj * zltv(ji,jj,jk) -zfm_vj * zltv(ji,jj+1,jk)
-            END DO
-         END DO
-
-         !  Tracer flux divergence at t-point added to the general trend
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ! horizontal advective trends
-#if defined key_zco
-               zbtr = e1e2tr(ji,jj)
-#else
-               zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk)
-#endif
-               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
-               pta(ji,jj,jk) = pta(ji,jj,jk) + zta
-            END DO
-         END DO
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-
-      ! Horizontal trend used in tra_adv_ztvd subroutine
-      zltu(:,:,:) = pta(:,:,:) - ztrdt(:,:,:) 
-
-      !  Save the horizontal advective trends for diagnostic
-      ! -----------------------------------------------------
-      IF( l_trdtra ) THEN
+               zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji+1,jj  ,jk) )
+               zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji  ,jj+1,jk) )
+               ! UBS scheme
+               zwx(ji,jj,jk) = 0.5 * (  zcenut - zfp_ui * zltu(ji,jj,jk) -zfm_ui * zltu(ji+1,jj,jk)  )
+               zwy(ji,jj,jk) = 0.5 * (  zcenvt - zfp_vj * zltv(ji,jj,jk) -zfm_vj * zltv(ji,jj+1,jk)  )
+            END DO
+         END DO
+      END DO
+
+      zltu(:,:,:) = pta(:,:,:)      ! store pta trends
+
+      ! Horizontal advective trends
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+               ! horizontal advective trends added to the general tracer trends
+               pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
+                  &                                    + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
+            END DO
+         END DO
+      END DO
+
+      zltu(:,:,:) = pta(:,:,:) - zltu(:,:,:)    ! store the Horizontal advective trend (used in tra_adv_ztvd subroutine)
+
+ 
+      IF( l_trdtra ) THEN                       ! trend diagnostics (from advective fluxes)
          CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptn )
          CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptn )
       ENDIF
-
-      ! "Poleward" heat or salt transport 
+      !                                         ! "Poleward" heat or salt transports
       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.
-                    zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-         ENDIF
-         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zwy(:,:,:) )
-         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zwy(:,:,:) )
-      ENDIF
-
+         IF( ktra == jp_tem )   pht_adv(:) = ptr_vj( zwy(:,:,:) )
+         IF( ktra == jp_sal )   pst_adv(:) = ptr_vj( zwy(:,:,:) )
+      ENDIF
+      !                                         ! control print
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' ubs - had: ', mask1=tmask, clinfo3=cdtype )
 
@@ -226,22 +170,22 @@
       ! II. Vertical advection
       ! ----------------------
-      IF( l_trdtra )   ztrdt(:,:,:) = pta(:,:,:)          ! Save ta and sa trends
+      IF( l_trdtra )   zltv(:,:,:) = pta(:,:,:)          ! store pta if trend diag.
     
       ! TVD scheme the vertical direction  
       CALL tra_adv_ztvd( kt, pwn, zltu, ptb, ptn, pta )
 
-      IF( l_trdtra )   THEN         !  vertical advective trend diagnostics
-         DO jk = 1, jpkm1
+      IF( l_trdtra )   THEN                     ! vertical advective trend diagnostics
+         DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z_hdivn = (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  ) / fse3t(ji,jj,jk)
-                  ztrdt(ji,jj,jk) = pta(ji,jj,jk) - ztrdt(ji,jj,jk+1)  +  ptn(ji,jj,jk) * z_hdivn 
+                  zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  z_hdivn = (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  ) * zbtr
+                  zltv(ji,jj,jk) = pta(ji,jj,jk) - zltv(ji,jj,jk+1)  +  ptn(ji,jj,jk) * z_hdivn 
                END DO
             END DO
          END DO
-         CALL trd_tra( kt, ktra, jpt_trd_zad, cdtype, ptrd3d=ztrdt )
-         !
-      ENDIF
-     
+         CALL trd_tra( kt, ktra, jpt_trd_zad, cdtype, ptrd3d=zltv )
+      ENDIF
+      !                                         ! control print
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' ubs - zad: ', mask1=tmask, clinfo3=cdtype )
       !
@@ -249,27 +193,26 @@
 
 
-   SUBROUTINE tra_adv_ztvd( kt, pwn, zttrd, ptb, ptn, pta )
+   SUBROUTINE tra_adv_ztvd( kt, pwn, phtrd, ptb, ptn, pta )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_ztvd  ***
       !! 
-      !! **  Purpose :   Compute the now trend due to total advection of 
-      !!       tracers and add it to the general trend of tracer equations
+      !! **  Purpose :   Compute the vertical advective trend of a tracer
+      !!               and add it to its general trend 
       !!
       !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with
       !!       corrected flux (monotonic correction)
-      !!       note: - this advection scheme needs a leap-frog time scheme
-      !!
-      !! ** Action : - update (ta,sa) with the now advective tracer trends
-      !!             - save the trends in (ztrdt,ztrds) ('key_trdtra')
-      !!----------------------------------------------------------------------
-      INTEGER , INTENT(in   )                         ::   kt      ! ocean time-step
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwn     ! verical effective velocity
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   zttrd   ! lateral advective trends on T & S 
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   ptb, ptn        ! before and now tracer fields
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend 
+      !!
+      !! ** Action : - update pta with the vertical advective trend
+      !!             - trend diagnostic (lk_trdtra=T)
+      !!----------------------------------------------------------------------
+      INTEGER , INTENT(in   )                         ::   kt         ! ocean time-step
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwn        ! verical effective velocity
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   phtrd      ! lateral advective trends on T & S 
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   ptb, ptn   ! before and now tracer fields
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta        ! tracer trend 
       !!
       INTEGER  ::   ji, jj, jk              ! dummy loop indices
-      REAL(wp) ::   z2dtt, zbtr, zew, z2    ! temporary scalar  
-      REAL(wp) ::   ztak, zfp_wk, zfm_wk            !    "         "
+      REAL(wp) ::   z2dtt, zbtr, z2         ! temporary scalar  
+      REAL(wp) ::   ztak, zfp_wk, zfm_wk    !    "         "
       REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zti, ztw   ! temporary 3D workspace
       !!----------------------------------------------------------------------
@@ -285,26 +228,19 @@
       ENDIF
 
-      !  Bottom value : flux set to zero
-      ! --------------
-      ztw(:,:,jpk) = 0.e0   ;   zti  (:,:,:) = 0.e0
-
 
       !  upstream advection with initial mass fluxes & intermediate update
-      ! -------------------------------------------------------------------
-      ! Surface value
-      IF( lk_dynspg_rl .OR. lk_vvl ) THEN                           ! rigid lid : flux set to zero
-         ztw(:,:,1) = 0.e0
-      ELSE                                                          ! free surface
-         ztw(:,:,1) = e1t(:,:) * e2t(:,:) * pwn(:,:,1) * ptb(:,:,1)
-      ENDIF
-
-      ! Interior value
-      DO jk = 2, jpkm1
+      ztw(:,:,jpk) = 0.e0                                   ! Bottom  value : flux set to zero
+      zti(:,:,jpk) = 0.e0
+      !                                                     ! Surface value : 
+      IF( lk_dynspg_rl .OR. lk_vvl ) THEN   ;   ztw(:,:, 1 ) = 0.e0                      ! rigid lid or non-linear fs
+      ELSE                                  ;   ztw(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1)   ! linear free surface 
+      ENDIF
+      !
+      DO jk = 2, jpkm1                                      ! interior values
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
-               zfp_wk = zew + ABS( zew )
-               zfm_wk = zew - ABS( zew )
-               ztw(ji,jj,jk) = zfp_wk * ptb(ji,jj,jk) + zfm_wk * ptb(ji,jj,jk-1)
+               zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
+               zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
+               ztw(ji,jj,jk) = 0.5 * (  zfp_wk * ptb(ji,jj,jk) + zfm_wk * ptb(ji,jj,jk-1)  )
             END DO
          END DO
@@ -318,42 +254,33 @@
                zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
                ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
-               pta(ji,jj,jk) =  pta(ji,jj,jk) + ztak
-               zti (ji,jj,jk) = ( ptb(ji,jj,jk) + z2dtt * ( ztak + zttrd(ji,jj,jk) ) ) * tmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-
-      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
-      CALL lbc_lnk( zti, 'T', 1. )
+               pta(ji,jj,jk) =   pta(ji,jj,jk) +           ztak
+               zti(ji,jj,jk) = ( ptb(ji,jj,jk) + z2dtt * ( ztak + phtrd(ji,jj,jk) ) ) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
 
 
       !  antidiffusive flux : high order minus low order
-      ! -------------------------------------------------      
       ztw(:,:,1) = 0.e0       ! Surface value
-
       DO jk = 2, jpkm1        ! Interior value
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
-               ztw(ji,jj,jk) = zew * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) - ztw(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-
-      !  monotonicity algorithm
-      ! ------------------------
-      CALL nonosc_z( ptb, ztw, zti, z2 )
+               ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) - ztw(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      CALL nonosc_z( ptb, ztw, zti, z2 )      !  monotonicity algorithm
 
 
       !  final trend with corrected fluxes
-      ! -----------------------------------
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.  
                zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-               ! k- vertical advective trends
-               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
-               ! add them to the general tracer trends
-               pta(ji,jj,jk) = pta(ji,jj,jk) + ztak
+               ! k- vertical advective trends added to the general tracer trends
+               pta(ji,jj,jk) = pta(ji,jj,jk) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
             END DO
          END DO
