Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traadv_mus.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traadv_mus.F90	(revision 12193)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traadv_mus.F90	(revision 12340)
@@ -47,4 +47,5 @@
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -131,12 +132,8 @@
          zwx(:,:,jpk) = 0._wp                   ! bottom values
          zwy(:,:,jpk) = 0._wp  
-         DO jk = 1, jpkm1                       ! interior values
-            DO jj = 1, jpjm1      
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwx(ji,jj,jk) = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
-                  zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
-               END DO
-           END DO
-         END DO
+         DO_3D_10_10( 1, jpkm1 )
+            zwx(ji,jj,jk) = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
+            zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
+         END_3D
          ! lateral boundary conditions   (changed sign)
          CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. )
@@ -144,59 +141,43 @@
          zslpx(:,:,jpk) = 0._wp                 ! bottom values
          zslpy(:,:,jpk) = 0._wp
-         DO jk = 1, jpkm1                       ! interior values
-            DO jj = 2, jpj
-               DO ji = fs_2, jpi   ! vector opt.
-                  zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
-                     &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
-                  zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   &
-                     &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
-               END DO
-            END DO
-         END DO
-         !
-         DO jk = 1, jpkm1                 !-- Slopes limitation
-            DO jj = 2, jpj
-               DO ji = fs_2, jpi   ! vector opt.
-                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
-                     &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   &
-                     &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) )
-                  zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
-                     &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   &
-                     &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) )
-               END DO
-           END DO
-         END DO
-         !
-         DO jk = 1, jpkm1                 !-- MUSCL horizontal advective fluxes
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ! MUSCL fluxes
-                  z0u = SIGN( 0.5, pU(ji,jj,jk) )
-                  zalpha = 0.5 - z0u
-                  zu  = z0u - 0.5 * pU(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
-                  zzwx = pt(ji+1,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk)
-                  zzwy = pt(ji  ,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk)
-                  zwx(ji,jj,jk) = pU(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
-                  !
-                  z0v = SIGN( 0.5, pV(ji,jj,jk) )
-                  zalpha = 0.5 - z0v
-                  zv  = z0v - 0.5 * pV(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
-                  zzwx = pt(ji,jj+1,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk)
-                  zzwy = pt(ji,jj  ,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk)
-                  zwy(ji,jj,jk) = pV(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
-               END DO
-            END DO
-         END DO
+         DO_3D_01_01( 1, jpkm1 )
+            zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
+               &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
+            zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   &
+               &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
+         END_3D
+         !
+         DO_3D_01_01( 1, jpkm1 )
+            zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
+               &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   &
+               &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) )
+            zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
+               &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   &
+               &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) )
+         END_3D
+         !
+         DO_3D_00_00( 1, jpkm1 )
+            ! MUSCL fluxes
+            z0u = SIGN( 0.5, pU(ji,jj,jk) )
+            zalpha = 0.5 - z0u
+            zu  = z0u - 0.5 * pU(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
+            zzwx = pt(ji+1,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk)
+            zzwy = pt(ji  ,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk)
+            zwx(ji,jj,jk) = pU(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
+            !
+            z0v = SIGN( 0.5, pV(ji,jj,jk) )
+            zalpha = 0.5 - z0v
+            zv  = z0v - 0.5 * pV(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
+            zzwx = pt(ji,jj+1,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk)
+            zzwy = pt(ji,jj  ,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk)
+            zwy(ji,jj,jk) = pV(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
+         END_3D
          CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. )   ! lateral boundary conditions   (changed sign)
          !
-         DO jk = 1, jpkm1                 !-- Tracer advective trend
-            DO jj = 2, jpjm1      
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  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  ) )     &
-                  &                                   * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
-               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  ) )     &
+            &                                   * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+         END_3D
          !                                ! trend diagnostics
          IF( l_trd )  THEN
@@ -219,40 +200,26 @@
          !                                !-- Slopes of tracer
          zslpx(:,:,1) = 0._wp                   ! surface values
-         DO jk = 2, jpkm1                       ! interior value
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zslpx(ji,jj,jk) =                     ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )  &
-                     &            * (  0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  )
-               END DO
-            END DO
-         END DO
-         DO jk = 2, jpkm1                 !-- Slopes limitation
-            DO jj = 1, jpj                      ! interior values
-               DO ji = 1, jpi
-                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
-                     &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   &
-                     &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  )
-               END DO
-            END DO
-         END DO
-         DO jk = 1, jpk-2                 !-- vertical advective flux
-            DO jj = 2, jpjm1      
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z0w = SIGN( 0.5, pW(ji,jj,jk+1) )
-                  zalpha = 0.5 + z0w
-                  zw  = z0w - 0.5 * pW(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm)
-                  zzwx = pt(ji,jj,jk+1,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1)
-                  zzwy = pt(ji,jj,jk  ,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  )
-                  zwx(ji,jj,jk+1) = pW(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)
-               END DO 
-            END DO
-         END DO
+         DO_3D_11_11( 2, jpkm1 )
+            zslpx(ji,jj,jk) =                     ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )  &
+               &            * (  0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  )
+         END_3D
+         DO_3D_11_11( 2, jpkm1 )
+            zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
+               &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   &
+               &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  )
+         END_3D
+         DO_3D_00_00( 1, jpk-2 )
+            z0w = SIGN( 0.5, pW(ji,jj,jk+1) )
+            zalpha = 0.5 + z0w
+            zw  = z0w - 0.5 * pW(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm)
+            zzwx = pt(ji,jj,jk+1,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1)
+            zzwy = pt(ji,jj,jk  ,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  )
+            zwx(ji,jj,jk+1) = pW(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)
+         END_3D
          IF( ln_linssh ) THEN                   ! top values, linear free surface only
             IF( ln_isfcav ) THEN                      ! ice-shelf cavities (top of the ocean)
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zwx(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)
-                  END DO
-               END DO   
+               DO_2D_11_11
+                  zwx(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)
+               END_2D
             ELSE                                      ! no cavities: only at the ocean surface
                zwx(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
@@ -260,11 +227,7 @@
          ENDIF
          !
-         DO jk = 1, jpkm1                 !-- vertical advective trend
-            DO jj = 2, jpjm1      
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
-               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,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+         END_3D
          !                                ! send trends for diagnostic
          IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwx, pW, pt(:,:,:,jn,Kbb) )
