Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/dommsk.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/dommsk.F90	(revision 12250)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/dommsk.F90	(revision 12340)
@@ -45,4 +45,5 @@
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -131,13 +132,11 @@
       !
       tmask(:,:,:) = 0._wp
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            iktop = k_top(ji,jj)
-            ikbot = k_bot(ji,jj)
-            IF( iktop /= 0 ) THEN       ! water in the column
-               tmask(ji,jj,iktop:ikbot  ) = 1._wp
-            ENDIF
-         END DO  
-      END DO
+      DO_2D_11_11
+         iktop = k_top(ji,jj)
+         ikbot = k_bot(ji,jj)
+         IF( iktop /= 0 ) THEN       ! water in the column
+            tmask(ji,jj,iktop:ikbot  ) = 1._wp
+         ENDIF
+      END_2D
       !
       ! the following call is mandatory
@@ -155,11 +154,7 @@
          CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) )
          CALL iom_close( inum )
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj)
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 1, jpkm1 )
+            tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj)
+         END_3D
       ENDIF
          
@@ -243,12 +238,10 @@
          DO jk = 1, jpk
             zwf(:,:) = fmask(:,:,jk)         
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  IF( fmask(ji,jj,jk) == 0._wp ) THEN
-                     fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
-                        &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
-                  ENDIF
-               END DO
-            END DO
+            DO_2D_00_00
+               IF( fmask(ji,jj,jk) == 0._wp ) THEN
+                  fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
+                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
+               ENDIF
+            END_2D
             DO jj = 2, jpjm1
                IF( fmask(1,jj,jk) == 0._wp ) THEN
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 12250)
+++ 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
 
Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/domwri.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/domwri.F90	(revision 12250)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/domwri.F90	(revision 12340)
@@ -34,4 +34,5 @@
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -100,23 +101,17 @@
       
       CALL dom_uniq( zprw, 'T' )
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zprt(ji,jj) = ssmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
-         END DO
-      END DO                             !    ! unique point mask
+      DO_2D_11_11
+         zprt(ji,jj) = ssmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
+      END_2D
       CALL iom_rstput( 0, 0, inum, 'tmaskutil', zprt, ktype = jp_i1 )  
       CALL dom_uniq( zprw, 'U' )
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zprt(ji,jj) = ssumask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
-         END DO
-      END DO
+      DO_2D_11_11
+         zprt(ji,jj) = ssumask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
+      END_2D
       CALL iom_rstput( 0, 0, inum, 'umaskutil', zprt, ktype = jp_i1 )  
       CALL dom_uniq( zprw, 'V' )
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zprt(ji,jj) = ssvmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
-         END DO
-      END DO
+      DO_2D_11_11
+         zprt(ji,jj) = ssvmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
+      END_2D
       CALL iom_rstput( 0, 0, inum, 'vmaskutil', zprt, ktype = jp_i1 )  
 !!gm  ssfmask has been removed  ==>> find another solution to defined fmaskutil
Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/domzgr.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/domzgr.F90	(revision 12250)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/domzgr.F90	(revision 12340)
@@ -45,4 +45,5 @@
   !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -150,12 +151,10 @@
       !
       !                                ! ice shelf draft and bathymetry
-      DO jj = 1,jpj
-         DO ji = 1,jpi
-            ikt = mikt(ji,jj)
-            ikb = mbkt(ji,jj)
-            bathy  (ji,jj) = gdepw_0(ji,jj,ikb+1)
-            risfdep(ji,jj) = gdepw_0(ji,jj,ikt  )
-         END DO
-      END DO
+      DO_2D_11_11
+         ikt = mikt(ji,jj)
+         ikb = mbkt(ji,jj)
+         bathy  (ji,jj) = gdepw_0(ji,jj,ikb+1)
+         risfdep(ji,jj) = gdepw_0(ji,jj,ikt  )
+      END_2D
       !
       !                                ! deepest/shallowest W level Above/Below ~10m
@@ -315,14 +314,12 @@
       !                                    ! N.B.  top     k-index of W-level = mikt
       !                                    !       bottom  k-index of W-level = mbkt+1
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1
-            miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  )
-            mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  )
-            mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  )
-            !
-            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
-            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
-         END DO
-      END DO
+      DO_2D_10_10
+         miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  )
+         mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  )
+         mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  )
+         !
+         mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
+         mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
+      END_2D
       ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 
       zk(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'U', 1. )   ;   miku(:,:) = MAX( NINT( zk(:,:) ), 1 )
Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/dtatsd.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/dtatsd.F90	(revision 12250)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/dtatsd.F90	(revision 12340)
@@ -35,4 +35,6 @@
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsd   ! structure of input SST (file informations, fields read)
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -184,32 +186,30 @@
          ENDIF
          !
-         DO jj = 1, jpj                         ! vertical interpolation of T & S
-            DO ji = 1, jpi
-               DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
-                  zl = gdept_0(ji,jj,jk)
-                  IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data
-                     ztp(jk) =  ptsd(ji,jj,1    ,jp_tem)
-                     zsp(jk) =  ptsd(ji,jj,1    ,jp_sal)
-                  ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data
-                     ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem)
-                     zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal)
-                  ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
-                     DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
-                        IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
-                           zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
-                           ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi 
-                           zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi
-                        ENDIF
-                     END DO
-                  ENDIF
-               END DO
-               DO jk = 1, jpkm1
-                  ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord
-                  ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk)
-               END DO
-               ptsd(ji,jj,jpk,jp_tem) = 0._wp
-               ptsd(ji,jj,jpk,jp_sal) = 0._wp
+         DO_2D_11_11
+            DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
+               zl = gdept_0(ji,jj,jk)
+               IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data
+                  ztp(jk) =  ptsd(ji,jj,1    ,jp_tem)
+                  zsp(jk) =  ptsd(ji,jj,1    ,jp_sal)
+               ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data
+                  ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem)
+                  zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal)
+               ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
+                  DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
+                     IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
+                        zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
+                        ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi 
+                        zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi
+                     ENDIF
+                  END DO
+               ENDIF
             END DO
-         END DO
+            DO jk = 1, jpkm1
+               ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord
+               ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk)
+            END DO
+            ptsd(ji,jj,jpk,jp_tem) = 0._wp
+            ptsd(ji,jj,jpk,jp_sal) = 0._wp
+         END_2D
          ! 
       ELSE                                !==   z- or zps- coordinate   ==!
@@ -219,20 +219,18 @@
          !
          IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ik = mbkt(ji,jj) 
-                  IF( ik > 1 ) THEN
-                     zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
-                     ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem)
-                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal)
-                  ENDIF
-                  ik = mikt(ji,jj)
-                  IF( ik > 1 ) THEN
-                     zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) ) 
-                     ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem)
-                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal)
-                  END IF
-               END DO
-            END DO
+            DO_2D_11_11
+               ik = mbkt(ji,jj) 
+               IF( ik > 1 ) THEN
+                  zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
+                  ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem)
+                  ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal)
+               ENDIF
+               ik = mikt(ji,jj)
+               IF( ik > 1 ) THEN
+                  zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) ) 
+                  ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem)
+                  ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal)
+               END IF
+            END_2D
          ENDIF
          !
Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/istate.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/istate.F90	(revision 12250)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/istate.F90	(revision 12340)
@@ -43,4 +43,5 @@
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -105,11 +106,9 @@
                ! Apply minimum wetdepth criterion
                !
-               DO jj = 1,jpj
-                  DO ji = 1,jpi
-                     IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN
-                        ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
-                     ENDIF
-                  END DO
-               END DO 
+               DO_2D_11_11
+                  IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN
+                     ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
+                  ENDIF
+               END_2D
             ENDIF 
             uu  (:,:,:,Kbb) = 0._wp
@@ -161,15 +160,11 @@
       !
 !!gm  the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk)
-               vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
-               !
-               uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk)
-               vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
+      DO_3D_11_11( 1, jpkm1 )
+         uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk)
+         vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
+         !
+         uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk)
+         vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk)
+      END_3D
       !
       uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
