Index: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90
===================================================================
--- branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90	(revision 7349)
+++ branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90	(revision 7351)
@@ -38,5 +38,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -78,6 +77,7 @@
       !!            prevent the appearance of spurious numerical oscillations
       !!
-      !! ** Action : - update (pta) with the now advective tracer trends
-      !!             - save the trends 
+      !! ** Action : - update pta  with the now advective tracer trends
+      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
+      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
       !!
       !! ** Reference : Leonard (1979, 1991)
@@ -87,5 +87,5 @@
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
-      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
+      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
@@ -105,9 +105,9 @@
       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
       !
-      ! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
+      !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
       CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 
       CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 
 
-      ! II. The vertical fluxes are computed with the 2nd order centered scheme
+      !        ! vertical fluxes are computed with the 2nd order centered scheme
       CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt )
       !
@@ -125,5 +125,5 @@
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
-      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of tracer time-step
+      REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step
       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun        ! i-velocity components
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
@@ -131,5 +131,5 @@
       !!
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp) ::   ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
+      REAL(wp) ::   ztra, zbtr, zdir, zdx, zmsk   ! local scalars
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx, zfu, zfc, zfd
       !----------------------------------------------------------------------
@@ -166,10 +166,9 @@
          !
          DO jk = 1, jpkm1  
-            zdt =  p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.   
                   zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0 
-                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk)
-                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
+                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u_n(ji,jj,jk)
+                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
                   zfc(ji,jj,jk)  = zdir * ptb(ji  ,jj,jk,jn) + ( 1. - zdir ) * ptb(ji+1,jj,jk,jn)  ! FC in the x-direction for T
                   zfd(ji,jj,jk)  = zdir * ptb(ji+1,jj,jk,jn) + ( 1. - zdir ) * ptb(ji  ,jj,jk,jn)  ! FD in the x-direction for T
@@ -216,5 +215,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.  
-                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
                   ! horizontal advective trends
                   ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
@@ -224,5 +223,5 @@
             END DO
          END DO
-         !                                 ! trend diagnostics (contribution of upstream fluxes)
+         !                                 ! trend diagnostics
          IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
          !
@@ -242,5 +241,5 @@
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
-      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of tracer time-step
+      REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step
       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
@@ -248,5 +247,5 @@
       !!
       INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
+      REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk   ! local scalars
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zwy, zfu, zfc, zfd
       !----------------------------------------------------------------------
@@ -289,10 +288,9 @@
          !
          DO jk = 1, jpkm1  
-            zdt =  p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.   
                   zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0 
-                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)
-                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
+                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v_n(ji,jj,jk)
+                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
                   zfc(ji,jj,jk)  = zdir * ptb(ji,jj  ,jk,jn) + ( 1. - zdir ) * ptb(ji,jj+1,jk,jn)  ! FC in the x-direction for T
                   zfd(ji,jj,jk)  = zdir * ptb(ji,jj+1,jk,jn) + ( 1. - zdir ) * ptb(ji,jj  ,jk,jn)  ! FD in the x-direction for T
@@ -340,5 +338,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.  
-                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
                   ! horizontal advective trends
                   ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
@@ -348,5 +346,5 @@
             END DO
          END DO
-         !                                 ! trend diagnostics (contribution of upstream fluxes)
+         !                                 ! trend diagnostics
          IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
@@ -381,7 +379,6 @@
       CALL wrk_alloc( jpi,jpj,jpk,   zwz )
       !
-      !                          ! surface & bottom values 
-      IF( lk_vvl )   zwz(:,:, 1 ) = 0._wp             ! set to zero one for all
-                     zwz(:,:,jpk) = 0._wp             ! except at the surface in linear free surface
+      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
+      zwz(:,:,jpk) = 0._wp
       !
       !                                                          ! ===========
@@ -396,5 +393,5 @@
             END DO
          END DO
-         IF(.NOT.lk_vvl ) THEN               !* top value   (only in linear free surf. as zwz is multiplied by wmask)
+         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
             IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
                DO jj = 1, jpj
@@ -403,5 +400,5 @@
                   END DO
                END DO   
-            ELSE                                   ! no ice-shelf cavities (only ocean surface)
+            ELSE                                   ! no ocean cavities (only ocean surface)
                zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
             ENDIF
@@ -412,14 +409,14 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
-                     &                                / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
-               END DO
-            END DO
-         END DO
-         !                                 ! Save the vertical advective trends for diagnostic
+                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !                                 ! Send trends for diagnostic
          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
          !
       END DO
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zwz )
+      CALL wrk_dealloc( jpi,jpj,jpk,   zwz )
       !
    END SUBROUTINE tra_adv_cen2_k
