Index: branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_cen2.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_cen2.F90	(revision 786)
+++ 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
