Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/domvvl.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/domvvl.F90	(revision 12150)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/domvvl.F90	(revision 12340)
@@ -65,4 +65,5 @@
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -190,22 +191,18 @@
       gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb)
       gdepw(:,:,1,Kbb) = 0.0_wp
-      DO jk = 2, jpk                               ! vertical sum
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
-               !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
-               !                             ! 0.5 where jk = mikt     
+      DO_3D_11_11( 2, jpk )
+         !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
+         !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
+         !                             ! 0.5 where jk = mikt     
 !!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ??
-               zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
-               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
-               gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  &
-                  &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm)) 
-               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
-               gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb)
-               gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  &
-                  &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb)) 
-            END DO
-         END DO
-      END DO
+         zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
+         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
+         gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  &
+            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm)) 
+         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
+         gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb)
+         gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  &
+            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb)) 
+      END_3D
       !
       !                    !==  thickness of the water column  !!   (ocean portion only)
@@ -242,27 +239,25 @@
          ENDIF
          IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator
-            DO jj = 1, jpj
-               DO ji = 1, jpi
+            DO_2D_11_11
 !!gm  case |gphi| >= 6 degrees is useless   initialized just above by default
-                  IF( ABS(gphit(ji,jj)) >= 6.) THEN
-                     ! values outside the equatorial band and transition zone (ztilde)
-                     frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp )
-                     frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
-                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star
-                     ! values inside the equatorial band (ztilde as zstar)
-                     frq_rst_e3t(ji,jj) =  0.0_wp
-                     frq_rst_hdv(ji,jj) =  1.0_wp / rdt
-                  ELSE                                      ! transition band (2.5 to 6 degrees N/S)
-                     !                                      ! (linearly transition from z-tilde to z-star)
-                     frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   &
-                        &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
-                        &                                          * 180._wp / 3.5_wp ) )
-                     frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                &
-                        &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   &
-                        &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
-                        &                                          * 180._wp / 3.5_wp ) )
-                  ENDIF
-               END DO
-            END DO
+               IF( ABS(gphit(ji,jj)) >= 6.) THEN
+                  ! values outside the equatorial band and transition zone (ztilde)
+                  frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp )
+                  frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
+               ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star
+                  ! values inside the equatorial band (ztilde as zstar)
+                  frq_rst_e3t(ji,jj) =  0.0_wp
+                  frq_rst_hdv(ji,jj) =  1.0_wp / rdt
+               ELSE                                      ! transition band (2.5 to 6 degrees N/S)
+                  !                                      ! (linearly transition from z-tilde to z-star)
+                  frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   &
+                     &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
+                     &                                          * 180._wp / 3.5_wp ) )
+                  frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                &
+                     &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   &
+                     &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
+                     &                                          * 180._wp / 3.5_wp ) )
+               ENDIF
+            END_2D
             IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN
                IF( nn_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2
@@ -413,31 +408,21 @@
          zwu(:,:) = 0._wp
          zwv(:,:) = 0._wp
-         DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
-                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
-                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
-                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
-                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
-                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         DO jj = 1, jpj          ! b - correction for last oceanic u-v points
-            DO ji = 1, jpi
-               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
-               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
-            END DO
-         END DO
-         DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
-                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
-                     &                                            ) * r1_e1e2t(ji,jj)
-               END DO
-            END DO
-         END DO
+         DO_3D_10_10( 1, jpkm1 )
+            un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
+               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
+            vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
+               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
+            zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
+            zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
+         END_3D
+         DO_2D_11_11
+            un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
+            vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
+         END_2D
+         DO_3D_00_00( 1, jpkm1 )
+            tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
+               &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
+               &                                            ) * r1_e1e2t(ji,jj)
+         END_3D
          !                       ! d - thickness diffusion transport: boundary conditions
          !                             (stored for tracer advction and continuity equation)
@@ -670,17 +655,13 @@
       gdepw(:,:,1,Kmm) = 0.0_wp
       gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
-      DO jk = 2, jpk
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
-                                                                 ! 1 for jk = mikt
-               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
-               gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  &
-                   &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) ) 
-               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
-            END DO
-         END DO
-      END DO
+      DO_3D_11_11( 2, jpk )
+        !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
+                                                           ! 1 for jk = mikt
+         zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
+         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
+         gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  &
+             &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) ) 
+         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
+      END_3D
 
       ! Local depth and Inverse of the local depth of the water
@@ -729,40 +710,28 @@
          !
       CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
-         DO jk = 1, jpk
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
-                     &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
-                     &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
-               END DO
-            END DO
-         END DO
+         DO_3D_10_10( 1, jpk )
+            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
+               &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
+               &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
+         END_3D
          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp )
          pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
          !
       CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
-         DO jk = 1, jpk
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
-                     &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
-                     &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
-               END DO
-            END DO
-         END DO
+         DO_3D_10_10( 1, jpk )
+            pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
+               &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
+               &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
+         END_3D
          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp )
          pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
          !
       CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
-         DO jk = 1, jpk
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
-                     &                       *    r1_e1e2f(ji,jj)                                                  &
-                     &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
-                     &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
-               END DO
-            END DO
-         END DO
+         DO_3D_10_10( 1, jpk )
+            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
+               &                       *    r1_e1e2f(ji,jj)                                                  &
+               &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
+               &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
+         END_3D
          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp )
          pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
@@ -926,12 +895,10 @@
                   ssh(:,:,Kbb) = -ssh_ref
 
-                  DO jj = 1, jpj
-                     DO ji = 1, jpi
-                        IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
-                           ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
-                           ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
-                        ENDIF
-                     ENDDO
-                  ENDDO
+                  DO_2D_11_11
+                     IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
+                        ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
+                        ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
+                     ENDIF
+                  END_2D
                ENDIF !If test case else
 
