Index: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/traadv_cen.F90
===================================================================
--- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/traadv_cen.F90	(revision 11139)
+++ NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/traadv_cen.F90	(revision 12928)
@@ -36,5 +36,5 @@
 
    !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -44,6 +44,6 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pun, pvn, pwn,     &
-      &                                        ptn, pta, kjpt, kn_cen_h, kn_cen_v ) 
+   SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pU, pV, pW,     &
+      &                    Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v ) 
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_cen  ***
@@ -59,17 +59,17 @@
       !!                = 4  ==>> 4th order COMPACT  scheme     -      -
       !!
-      !! ** Action : - update pta  with the now advective tracer trends
+      !! ** Action : - update pt(:,:,:,:,Krhs)  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)
+      !!             - poleward advective heat and salt transport (l_diaptr=T)
       !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
-      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
-      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
-      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
-      INTEGER                              , INTENT(in   ) ::   kn_cen_h        ! =2/4 (2nd or 4th order scheme)
-      INTEGER                              , INTENT(in   ) ::   kn_cen_v        ! =2/4 (2nd or 4th order scheme)
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn             ! now tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend 
+      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::   Kmm, Krhs       ! ocean time level indices
+      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index
+      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
+      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers
+      INTEGER                                  , INTENT(in   ) ::   kn_cen_h        ! =2/4 (2nd or 4th order scheme)
+      INTEGER                                  , INTENT(in   ) ::   kn_cen_v        ! =2/4 (2nd or 4th order scheme)
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -89,8 +89,8 @@
       l_hst = .FALSE.
       l_ptr = .FALSE.
-      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )        l_trd = .TRUE.
-      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE. 
+      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )       l_trd = .TRUE.
+      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )    l_ptr = .TRUE. 
       IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
-         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
+         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
       !
       !                    
@@ -103,40 +103,28 @@
          !
          CASE(  2  )                         !* 2nd order centered
-            DO jk = 1, jpkm1
-               DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1   ! vector opt.
-                     zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
-                     zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
-                  END DO
-               END DO
-            END DO
+            DO_3D_10_10( 1, jpkm1 )
+               zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm) )
+               zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) )
+            END_3D
             !
          CASE(  4  )                         !* 4th order centered
             ztu(:,:,jpk) = 0._wp                   ! Bottom value : flux set to zero
             ztv(:,:,jpk) = 0._wp
-            DO jk = 1, jpkm1                       ! masked gradient
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
-                     ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
+            DO_3D_00_00( 1, jpkm1 )
+               ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
+               ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
+            END_3D
             CALL lbc_lnk_multi( 'traadv_cen', ztu, 'U', -1. , ztv, 'V', -1. )   ! Lateral boundary cond.
             !
-            DO jk = 1, jpkm1                       ! Horizontal advective fluxes
-               DO jj = 2, jpjm1
-                  DO ji = 1, fs_jpim1   ! vector opt.
-                     zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! C2 interpolation of T at u- & v-points (x2)
-                     zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn)
-                     !                                                  ! C4 interpolation of T at u- & v-points (x2)
-                     zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj,jk) - ztu(ji+1,jj,jk) )
-                     zC4t_v =  zC2t_v + r1_6 * ( ztv(ji,jj-1,jk) - ztv(ji,jj+1,jk) )
-                     !                                                  ! C4 fluxes
-                     zwx(ji,jj,jk) =  0.5_wp * pun(ji,jj,jk) * zC4t_u
-                     zwy(ji,jj,jk) =  0.5_wp * pvn(ji,jj,jk) * zC4t_v
-                  END DO
-               END DO
-            END DO         
+            DO_3D_00_10( 1, jpkm1 )
+               zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! C2 interpolation of T at u- & v-points (x2)
+               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
+               !                                                  ! C4 interpolation of T at u- & v-points (x2)
+               zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj,jk) - ztu(ji+1,jj,jk) )
+               zC4t_v =  zC2t_v + r1_6 * ( ztv(ji,jj-1,jk) - ztv(ji,jj+1,jk) )
+               !                                                  ! C4 fluxes
+               zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * zC4t_u
+               zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v
+            END_3D
             !
          CASE DEFAULT
@@ -147,21 +135,13 @@
          !
          CASE(  2  )                         !* 2nd order centered
-            DO jk = 2, jpk
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-                     zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
+            DO_3D_00_00( 2, jpk )
+               zwz(ji,jj,jk) = 0.5 * pW(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) * wmask(ji,jj,jk)
+            END_3D
             !
          CASE(  4  )                         !* 4th order compact
-            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )      ! ztw = interpolated value of T at w-point
-            DO jk = 2, jpkm1
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1
-                     zwz(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
+            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )      ! ztw = interpolated value of T at w-point
+            DO_3D_00_00( 2, jpkm1 )
+               zwz(ji,jj,jk) = pW(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
+            END_3D
             !
          END SELECT
@@ -169,29 +149,23 @@
          IF( ln_linssh ) THEN                !* top value   (linear free surf. only as zwz is multiplied by wmask)
             IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn) 
-                  END DO
-               END DO   
+               DO_2D_11_11
+                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) 
+               END_2D
             ELSE                                   ! no ice-shelf cavities (only ocean surface)
-               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
+               zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kmm)
             ENDIF
          ENDIF
          !               
-         DO jk = 1, jpkm1              !--  Divergence of advective fluxes  --!
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn)    &
-                     &             - (  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)  ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+         DO_3D_00_00( 1, jpkm1 )
+            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)    &
+               &             - (  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)  ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+         END_3D
          !                             ! trend diagnostics
          IF( l_trd ) THEN
-            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
-            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
-            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
+            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
+            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
+            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
          END IF
          !                                 ! "Poleward" heat and salt transports 
