Index: NEMO/trunk/src/OCE/TRD/trdtra.F90
===================================================================
--- NEMO/trunk/src/OCE/TRD/trdtra.F90	(revision 10425)
+++ NEMO/trunk/src/OCE/TRD/trdtra.F90	(revision 12377)
@@ -41,5 +41,5 @@
 
    !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -60,5 +60,5 @@
 
 
-   SUBROUTINE trd_tra( kt, ctype, ktra, ktrd, ptrd, pun, ptra )
+   SUBROUTINE trd_tra( kt, Kmm, Krhs, ctype, ktra, ktrd, ptrd, pu, ptra )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_tra  ***
@@ -77,6 +77,7 @@
       INTEGER                         , INTENT(in)           ::   ktra    ! tracer index
       INTEGER                         , INTENT(in)           ::   ktrd    ! tracer trend index
+      INTEGER                         , INTENT(in)           ::   Kmm, Krhs ! time level indices
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)           ::   ptrd    ! tracer trend  or flux
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pun     ! now velocity 
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pu      ! now velocity 
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   ptra    ! now tracer variable
       !
@@ -94,11 +95,11 @@
          SELECT CASE( ktrd )
          !                            ! advection: transform the advective flux into a trend
-         CASE( jptra_xad )   ;   CALL trd_tra_adv( ptrd, pun, ptra, 'X', trdtx ) 
-         CASE( jptra_yad )   ;   CALL trd_tra_adv( ptrd, pun, ptra, 'Y', trdty ) 
-         CASE( jptra_zad )   ;   CALL trd_tra_adv( ptrd, pun, ptra, 'Z', trdt  ) 
+         CASE( jptra_xad )   ;   CALL trd_tra_adv( ptrd, pu, ptra, 'X', trdtx, Kmm ) 
+         CASE( jptra_yad )   ;   CALL trd_tra_adv( ptrd, pu, ptra, 'Y', trdty, Kmm ) 
+         CASE( jptra_zad )   ;   CALL trd_tra_adv( ptrd, pu, ptra, 'Z', trdt, Kmm )
          CASE( jptra_bbc,    &        ! qsr, bbc: on temperature only, send to trd_tra_mng
             &  jptra_qsr )   ;   trdt(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
                                  ztrds(:,:,:) = 0._wp
-                                 CALL trd_tra_mng( trdt, ztrds, ktrd, kt )
+                                 CALL trd_tra_mng( trdt, ztrds, ktrd, kt, Kmm )
  !!gm Gurvan, verify the jptra_evd trend please !
          CASE( jptra_evd )   ;   avt_evd(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
@@ -114,10 +115,10 @@
          !                            ! advection: transform the advective flux into a trend
          !                            !            and send T & S trends to trd_tra_mng
-         CASE( jptra_xad  )   ;   CALL trd_tra_adv( ptrd , pun  , ptra, 'X'  , ztrds ) 
-                                  CALL trd_tra_mng( trdtx, ztrds, ktrd, kt   )
-         CASE( jptra_yad  )   ;   CALL trd_tra_adv( ptrd , pun  , ptra, 'Y'  , ztrds ) 
-                                  CALL trd_tra_mng( trdty, ztrds, ktrd, kt   )
-         CASE( jptra_zad  )   ;   CALL trd_tra_adv( ptrd , pun  , ptra, 'Z'  , ztrds ) 
-                                  CALL trd_tra_mng( trdt , ztrds, ktrd, kt   )
+         CASE( jptra_xad  )   ;   CALL trd_tra_adv( ptrd , pu  , ptra, 'X'  , ztrds, Kmm ) 
+                                  CALL trd_tra_mng( trdtx, ztrds, ktrd, kt, Kmm   )
+         CASE( jptra_yad  )   ;   CALL trd_tra_adv( ptrd , pu  , ptra, 'Y'  , ztrds, Kmm ) 
+                                  CALL trd_tra_mng( trdty, ztrds, ktrd, kt, Kmm   )
+         CASE( jptra_zad  )   ;   CALL trd_tra_adv( ptrd , pu  , ptra, 'Z'  , ztrds, Kmm ) 
+                                  CALL trd_tra_mng( trdt , ztrds, ktrd, kt, Kmm   )
          CASE( jptra_zdfp )           ! diagnose the "PURE" Kz trend (here: just before the swap)
             !                         ! iso-neutral diffusion case otherwise jptra_zdf is "PURE"
@@ -127,28 +128,28 @@
             zwt(:,:,jpk) = 0._wp   ;   zws(:,:,jpk) = 0._wp
             DO jk = 2, jpk
-               zwt(:,:,jk) = avt(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / e3w_n(:,:,jk) * tmask(:,:,jk)
-               zws(:,:,jk) = avs(:,:,jk) * ( tsa(:,:,jk-1,jp_sal) - tsa(:,:,jk,jp_sal) ) / e3w_n(:,:,jk) * tmask(:,:,jk)
+               zwt(:,:,jk) = avt(:,:,jk) * ( ts(:,:,jk-1,jp_tem,Krhs) - ts(:,:,jk,jp_tem,Krhs) ) / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
+               zws(:,:,jk) = avs(:,:,jk) * ( ts(:,:,jk-1,jp_sal,Krhs) - ts(:,:,jk,jp_sal,Krhs) ) / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
             END DO
             !
             ztrdt(:,:,jpk) = 0._wp   ;   ztrds(:,:,jpk) = 0._wp
             DO jk = 1, jpkm1
-               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / e3t_n(:,:,jk)
-               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / e3t_n(:,:,jk) 
+               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / e3t(:,:,jk,Kmm)
+               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / e3t(:,:,jk,Kmm) 
             END DO
-            CALL trd_tra_mng( ztrdt, ztrds, jptra_zdfp, kt )  
+            CALL trd_tra_mng( ztrdt, ztrds, jptra_zdfp, kt, Kmm )  
             !
             !                         ! Also calculate EVD trend at this point. 
             zwt(:,:,:) = 0._wp   ;   zws(:,:,:) = 0._wp            ! vertical diffusive fluxes
             DO jk = 2, jpk
-               zwt(:,:,jk) = avt_evd(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / e3w_n(:,:,jk) * tmask(:,:,jk)
-               zws(:,:,jk) = avt_evd(:,:,jk) * ( tsa(:,:,jk-1,jp_sal) - tsa(:,:,jk,jp_sal) ) / e3w_n(:,:,jk) * tmask(:,:,jk)
+               zwt(:,:,jk) = avt_evd(:,:,jk) * ( ts(:,:,jk-1,jp_tem,Krhs) - ts(:,:,jk,jp_tem,Krhs) ) / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
+               zws(:,:,jk) = avt_evd(:,:,jk) * ( ts(:,:,jk-1,jp_sal,Krhs) - ts(:,:,jk,jp_sal,Krhs) ) / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
             END DO
             !
             ztrdt(:,:,jpk) = 0._wp   ;   ztrds(:,:,jpk) = 0._wp
             DO jk = 1, jpkm1
-               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / e3t_n(:,:,jk)
-               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / e3t_n(:,:,jk) 
+               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / e3t(:,:,jk,Kmm)
+               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / e3t(:,:,jk,Kmm) 
             END DO
-            CALL trd_tra_mng( ztrdt, ztrds, jptra_evd, kt )  
+            CALL trd_tra_mng( ztrdt, ztrds, jptra_evd, kt, Kmm )  
             !
             DEALLOCATE( zwt, zws, ztrdt )
@@ -156,5 +157,5 @@
          CASE DEFAULT                 ! other trends: mask and send T & S trends to trd_tra_mng
             ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
-            CALL trd_tra_mng( trdt, ztrds, ktrd, kt )  
+            CALL trd_tra_mng( trdt, ztrds, ktrd, kt, Kmm )  
          END SELECT
       ENDIF
@@ -164,12 +165,12 @@
          SELECT CASE( ktrd )
          !                            ! advection: transform the advective flux into a masked trend
-         CASE( jptra_xad )   ;   CALL trd_tra_adv( ptrd , pun , ptra, 'X', ztrds ) 
-         CASE( jptra_yad )   ;   CALL trd_tra_adv( ptrd , pun , ptra, 'Y', ztrds ) 
-         CASE( jptra_zad )   ;   CALL trd_tra_adv( ptrd , pun , ptra, 'Z', ztrds ) 
+         CASE( jptra_xad )   ;   CALL trd_tra_adv( ptrd , pu , ptra, 'X', ztrds, Kmm ) 
+         CASE( jptra_yad )   ;   CALL trd_tra_adv( ptrd , pu , ptra, 'Y', ztrds, Kmm ) 
+         CASE( jptra_zad )   ;   CALL trd_tra_adv( ptrd , pu , ptra, 'Z', ztrds, Kmm ) 
          CASE DEFAULT                 ! other trends: just masked 
                                  ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
          END SELECT
          !                            ! send trend to trd_trc
-         CALL trd_trc( ztrds, ktra, ktrd, kt ) 
+         CALL trd_trc( ztrds, ktra, ktrd, kt, Kmm ) 
          !
       ENDIF
@@ -178,5 +179,5 @@
 
 
-   SUBROUTINE trd_tra_adv( pf, pun, ptn, cdir, ptrd )
+   SUBROUTINE trd_tra_adv( pf, pu, pt, cdir, ptrd, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_tra_adv  ***
@@ -191,8 +192,9 @@
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pf      ! advective flux in one direction
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pun     ! now velocity   in one direction
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   ptn     ! now or before tracer 
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu      ! now velocity   in one direction
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pt      ! now or before tracer 
       CHARACTER(len=1)                , INTENT(in   ) ::   cdir    ! X/Y/Z direction
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   ptrd    ! advective trend in one direction
+      INTEGER,  INTENT(in)                            ::   Kmm     ! time level index
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -211,18 +213,14 @@
       ptrd(:,:,jpk) = 0._wp
       !
-      DO jk = 1, jpkm1         ! advective trend
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ptrd(ji,jj,jk) = - (     pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik)                        &
-                 &                  - ( pun(ji,jj,jk) - pun(ji-ii,jj-ij,jk-ik) ) * ptn(ji,jj,jk)  )   &
-                 &              * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
+      DO_3D_00_00( 1, jpkm1 )
+         ptrd(ji,jj,jk) = - (     pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik)                        &
+           &                  - ( pu(ji,jj,jk) - pu(ji-ii,jj-ij,jk-ik) ) * pt(ji,jj,jk)  )   &
+           &              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
+      END_3D
       !
    END SUBROUTINE trd_tra_adv
 
 
-   SUBROUTINE trd_tra_mng( ptrdx, ptrdy, ktrd, kt )
+   SUBROUTINE trd_tra_mng( ptrdx, ptrdy, ktrd, kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_tra_mng  ***
@@ -236,4 +234,5 @@
       INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
       INTEGER                   , INTENT(in   ) ::   kt      ! time step
+      INTEGER                   , INTENT(in   ) ::   Kmm     ! time level index
       !!----------------------------------------------------------------------
 
@@ -243,11 +242,11 @@
 
       !                   ! 3D output of tracers trends using IOM interface
-      IF( ln_tra_trd )   CALL trd_tra_iom ( ptrdx, ptrdy, ktrd, kt )
+      IF( ln_tra_trd )   CALL trd_tra_iom ( ptrdx, ptrdy, ktrd, kt, Kmm )
 
       !                   ! Integral Constraints Properties for tracers trends                                       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_glo_trd )   CALL trd_glo( ptrdx, ptrdy, ktrd, 'TRA', kt )
+      IF( ln_glo_trd )   CALL trd_glo( ptrdx, ptrdy, ktrd, 'TRA', kt, Kmm )
 
       !                   ! Potential ENergy trends
-      IF( ln_PE_trd  )   CALL trd_pen( ptrdx, ptrdy, ktrd, kt, r2dt )
+      IF( ln_PE_trd  )   CALL trd_pen( ptrdx, ptrdy, ktrd, kt, r2dt, Kmm )
 
       !                   ! Mixed layer trends for active tracers
@@ -290,5 +289,5 @@
 
 
-   SUBROUTINE trd_tra_iom( ptrdx, ptrdy, ktrd, kt )
+   SUBROUTINE trd_tra_iom( ptrdx, ptrdy, ktrd, kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_tra_iom  ***
@@ -300,4 +299,5 @@
       INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
       INTEGER                   , INTENT(in   ) ::   kt      ! time step
+      INTEGER                   , INTENT(in   ) ::   Kmm     ! time level index
       !!
       INTEGER ::   ji, jj, jk   ! dummy loop indices
@@ -326,6 +326,6 @@
                                   IF( ln_linssh ) THEN                   ! cst volume : adv flux through z=0 surface
                                      ALLOCATE( z2dx(jpi,jpj), z2dy(jpi,jpj) )
-                                     z2dx(:,:) = wn(:,:,1) * tsn(:,:,1,jp_tem) / e3t_n(:,:,1)
-                                     z2dy(:,:) = wn(:,:,1) * tsn(:,:,1,jp_sal) / e3t_n(:,:,1)
+                                     z2dx(:,:) = ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) / e3t(:,:,1,Kmm)
+                                     z2dy(:,:) = ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) / e3t(:,:,1,Kmm)
                                      CALL iom_put( "ttrd_sad", z2dx )
                                      CALL iom_put( "strd_sad", z2dy )
