Index: trunk/NEMO/OPA_SRC/TRA/traadv_muscl2.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traadv_muscl2.F90	(revision 457)
+++ trunk/NEMO/OPA_SRC/TRA/traadv_muscl2.F90	(revision 503)
@@ -4,4 +4,6 @@
    !! Ocean active tracers:  horizontal & vertical advective trend
    !!==============================================================================
+   !! History :  9.0  !  02-06  (G. Madec) from traadv_muscl
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -9,5 +11,4 @@
    !!                    and vertical advection trends using MUSCL2 scheme
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and active tracers
    USE dom_oce         ! ocean space and time domain
@@ -32,7 +33,7 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
    !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -50,24 +51,18 @@
       !!
       !! ** Action  : - update (ta,sa) with the now advective tracer trends
-      !!              - save trends in (ztdta,ztdsa) ('key_trdtra')
-      !!
-      !! References :                
-      !!      Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
-      !!	IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
-      !!
-      !! History :
-      !!        !  06-00  (A.Estublier)  for passive tracers
-      !!        !  01-08  (E.Durand G.Madec)  adapted for T & S
-      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
-      !!   9.0  !  04-08  (C. Talandier) New trends organization
+      !!              - save trends in (ztrdt,ztrds) ('key_trdtra')
+      !!
+      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
+      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
       !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER , INTENT( in ) ::   kt          ! ocean time-step index
-      REAL(wp), INTENT( in ), DIMENSION(jpi,jpj,jpk) ::   &
-         pun, pvn, pwn                        ! now ocean velocity fields
-
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk               ! dummy loop indices
+      USE oce              , ztrdt => ua   ! use ua as workspace
+      USE oce              , ztrds => va   ! use va as workspace
+      !!
+      INTEGER , INTENT(in)                         ::   kt    ! ocean time-step index
+      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun   ! ocean velocity u-component
+      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pvn   ! ocean velocity v-component
+      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pwn   ! ocean velocity w-component
+      !!
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   &
          zu, zv, zw, zeu, zev,           &  
@@ -76,9 +71,8 @@
          zzt1, zzt2, zalpha,             &
          zzs1, zzs2, z2,                 &
-         zta, zsa
-      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   &
-         zt1, zt2, ztp1, ztp2,   &
-         zs1, zs2, zsp1, zsp2,   &
-         ztdta, ztdsa
+         zta, zsa,                       &
+         z_hdivn_x, z_hdivn_y, z_hdivn
+      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zt1, zt2, ztp1, ztp2   ! 3D workspace
+      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zs1, zs2, zsp1, zsp2   !  "      "
       !!----------------------------------------------------------------------
 
@@ -89,17 +83,7 @@
       ENDIF
 
-      IF( neuler == 0 .AND. kt == nit000 ) THEN
-          z2=1.
-      ELSE
-          z2=2.
-      ENDIF
-
-      ! Save ta and sa trends
-      IF( l_trdtra )   THEN
-         ztdta(:,:,:) = ta(:,:,:) 
-         ztdsa(:,:,:) = sa(:,:,:) 
-         l_adv = 'mu2'
-      ENDIF
-
+      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2 = 1.
+      ELSE                                        ;   z2 = 2.
+      ENDIF
 
       ! I. Horizontal advective fluxes
@@ -194,5 +178,5 @@
                zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
                zs1(ji,jj,jk) = zeu * ( zalpha * zzs1 + (1.-zalpha) * zzs2 )
-
+               !
                z0v = SIGN( 0.5, pvn(ji,jj,jk) )            
                zalpha = 0.5 - z0v
@@ -271,15 +255,4 @@
       CALL lbc_lnk( zt2, 'V', -1. )   ;   CALL lbc_lnk( zs2, 'V', -1. )
 
-      ! Save MUSCL fluxes to compute i- & j- horizontal 
-      ! advection trends in the MLD
-      IF( l_trdtra )   THEN
-         ! save i- terms
-         tladi(:,:,:) = zt1(:,:,:) 
-         sladi(:,:,:) = zs1(:,:,:) 
-         ! save j- terms
-         tladj(:,:,:) = zt2(:,:,:) 
-         sladj(:,:,:) = zs2(:,:,:) 
-      ENDIF
-
       ! Compute & add the horizontal advective trend
 
@@ -305,25 +278,59 @@
 
       ! Save the horizontal advective trends for diagnostic
-
-      IF( l_trdtra )   THEN
-         ! Recompute the horizontal advection zta & zsa trends computed 
-         ! at the step 2. above in making the difference between the new 
-         ! trends and the previous one ta()/sa - ztdta()/ztdsa() and add
-         ! the term tn()/sn()*hdivn() to recover the Uh gradh(T/S) trends
-         ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) + tn(:,:,:) * hdivn(:,:,:) 
-         ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) + sn(:,:,:) * hdivn(:,:,:)
-
-         CALL trd_mod(ztdta, ztdsa, jpttdlad, 'TRA', kt)
-
-         ! Save the new ta and sa trends
-         ztdta(:,:,:) = ta(:,:,:) 
-         ztdsa(:,:,:) = sa(:,:,:) 
-
-      ENDIF
-
-      IF(ln_ctl) THEN
-         CALL prt_ctl(tab3d_1=ta, clinfo1=' muscl2 had  - Ta: ', mask1=tmask, &
-            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
-      ENDIF
+      IF( l_trdtra ) THEN
+         ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0
+         !
+         ! T/S ZONAL advection trends
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  !-- Compute zonal divergence by splitting hdivn (see divcur.F90)
+                  !   N.B. This computation is not valid along OBCs (if any)
+#if defined key_zco
+                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
+                  z_hdivn_x = (  e2u(ji  ,jj) * pun(ji  ,jj,jk)                              &
+                     &         - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr
+#else
+                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  z_hdivn_x = (  e2u(ji  ,jj) * fse3u(ji  ,jj,jk) * pun(ji  ,jj,jk)          &
+                     &         - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr
+#endif
+                  ztrdt(ji,jj,jk) = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj,jk) ) + tn(ji,jj,jk) * z_hdivn_x
+                  ztrds(ji,jj,jk) = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj,jk) ) + sn(ji,jj,jk) * z_hdivn_x
+               END DO
+            END DO
+         END DO
+         CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt)
+
+         ! T/S MERIDIONAL advection trends
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  !-- Compute merid. divergence by splitting hdivn (see divcur.F90)
+                  !   N.B. This computation is not valid along OBCs (if any)
+#if defined key_zco
+                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
+                  z_hdivn_y = (  e1v(ji,jj  ) * pvn(ji,jj  ,jk)                              &
+                     &         - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr
+#else
+                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  z_hdivn_y = (  e1v(ji,  jj) * fse3v(ji,jj  ,jk) * pvn(ji,jj  ,jk)          &
+                     &         - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr
+#endif
+                  ztrdt(ji,jj,jk) = - zbtr * ( zt2(ji,jj,jk) - zt2(ji,jj-1,jk) ) + tn(ji,jj,jk) * z_hdivn_y          
+                  ztrds(ji,jj,jk) = - zbtr * ( zs2(ji,jj,jk) - zs2(ji,jj-1,jk) ) + sn(ji,jj,jk) * z_hdivn_y
+               END DO
+            END DO
+         END DO
+         CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt)
+
+         ! Save the up-to-date ta and sa trends
+         ztrdt(:,:,:) = ta(:,:,:) 
+         ztrds(:,:,:) = sa(:,:,:) 
+         !
+      ENDIF
+
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' muscl2 had  - Ta: ', mask1=tmask,   &
+         &                       tab3d_2=sa, clinfo2=              ' Sa: ', mask2=tmask, clinfo3='tra')
 
       ! "zonal" mean advective heat and salt transport
@@ -451,22 +458,35 @@
 
       ! Save the vertical advective trends for diagnostic
-
       IF( l_trdtra )   THEN
          ! Recompute the vertical advection zta & zsa trends computed 
          ! at the step 2. above in making the difference between the new 
-         ! trends and the previous one: ta()/sa - ztdta()/ztdsa() and substract
+         ! trends and the previous one: ta()/sa - ztrdt()/ztrds() and substract
          ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
-         ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) - tn(:,:,:) * hdivn(:,:,:) 
-         ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) - sn(:,:,:) * hdivn(:,:,:)
-
-         CALL trd_mod(ztdta, ztdsa, jpttdzad, 'TRA', kt)
-      ENDIF
-
-      IF(ln_ctl) THEN
-         CALL prt_ctl(tab3d_1=ta, clinfo1=' muscl2 zad  - Ta: ', mask1=tmask, &
-            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
-
-      ENDIF
-
+
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+#if defined key_zco
+                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
+                  z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk)
+                  z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk)
+#else
+                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  z_hdivn_x = e2u(ji,jj)*fse3u(ji,jj,jk)*pun(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk)*pun(ji-1,jj,jk)
+                  z_hdivn_y = e1v(ji,jj)*fse3v(ji,jj,jk)*pvn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk)*pvn(ji,jj-1,jk)
+#endif
+                  z_hdivn   = (z_hdivn_x + z_hdivn_y) * zbtr
+                  ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn 
+                  ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn
+               END DO
+            END DO
+         END DO
+         CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt)
+         !
+      ENDIF
+
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' muscl2 zad  - Ta: ', mask1=tmask,   &
+         &                       tab3d_2=sa, clinfo2=              ' Sa: ', mask2=tmask, clinfo3='tra' )
+      !
    END SUBROUTINE tra_adv_muscl2
 
