Index: NEMO/trunk/tests/ISOMIP/MY_SRC/usrdef_hgr.F90
===================================================================
--- NEMO/trunk/tests/ISOMIP/MY_SRC/usrdef_hgr.F90	(revision 12489)
+++ NEMO/trunk/tests/ISOMIP/MY_SRC/usrdef_hgr.F90	(revision 12740)
@@ -27,4 +27,6 @@
    PUBLIC   usr_def_hgr   ! called by domhgr.F90
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -75,33 +77,30 @@
       !
       !                       !==  grid point position  ==!   (in degrees)
-      DO jj = 1, jpj
-         DO ji = 1, jpi             ! longitude   (west coast at lon=0°)
-            plamt(ji,jj) = rn_e1deg * (  - 0.5 + REAL( ji-1 + nimpp-1 , wp )  )  
-            plamu(ji,jj) = rn_e1deg * (          REAL( ji-1 + nimpp-1 , wp )  )
-            plamv(ji,jj) = plamt(ji,jj)
-            plamf(ji,jj) = plamu(ji,jj)
-            !                       ! latitude   (south coast at lat= 81°)
-            pphit(ji,jj) = rn_e2deg * (  - 0.5 + REAL( jj-1 + njmpp-1 , wp )  ) - 80._wp
-            pphiu(ji,jj) = pphit(ji,jj)
-            pphiv(ji,jj) = rn_e2deg * (          REAL( jj-1 + njmpp-1 , wp )  ) - 80_wp
-            pphif(ji,jj) = pphiv(ji,jj)
-         END DO
-      END DO
+      DO_2D_11_11
+         !                       ! longitude   (west coast at lon=0°)
+         plamt(ji,jj) = rn_e1deg * (  - 0.5 + REAL( ji-1 + nimpp-1 , wp )  )  
+         plamu(ji,jj) = rn_e1deg * (          REAL( ji-1 + nimpp-1 , wp )  )
+         plamv(ji,jj) = plamt(ji,jj)
+         plamf(ji,jj) = plamu(ji,jj)
+         !                       ! latitude   (south coast at lat= 81°)
+         pphit(ji,jj) = rn_e2deg * (  - 0.5 + REAL( jj-1 + njmpp-1 , wp )  ) - 80._wp
+         pphiu(ji,jj) = pphit(ji,jj)
+         pphiv(ji,jj) = rn_e2deg * (          REAL( jj-1 + njmpp-1 , wp )  ) - 80_wp
+         pphif(ji,jj) = pphiv(ji,jj)
+      END_2D
       !
       !                       !==  Horizontal scale factors  ==!   (in meters)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !                       ! e1   (zonal)
-            pe1t(ji,jj) = ra * rad * COS( rad * pphit(ji,jj) ) * rn_e1deg
-            pe1u(ji,jj) = ra * rad * COS( rad * pphiu(ji,jj) ) * rn_e1deg
-            pe1v(ji,jj) = ra * rad * COS( rad * pphiv(ji,jj) ) * rn_e1deg
-            pe1f(ji,jj) = ra * rad * COS( rad * pphif(ji,jj) ) * rn_e1deg
-            !                       ! e2   (meridional)
-            pe2t(ji,jj) = ra * rad * rn_e2deg
-            pe2u(ji,jj) = ra * rad * rn_e2deg
-            pe2v(ji,jj) = ra * rad * rn_e2deg
-            pe2f(ji,jj) = ra * rad * rn_e2deg
-         END DO
-      END DO
+      DO_2D_11_11
+         !                       ! e1   (zonal)
+         pe1t(ji,jj) = ra * rad * COS( rad * pphit(ji,jj) ) * rn_e1deg
+         pe1u(ji,jj) = ra * rad * COS( rad * pphiu(ji,jj) ) * rn_e1deg
+         pe1v(ji,jj) = ra * rad * COS( rad * pphiv(ji,jj) ) * rn_e1deg
+         pe1f(ji,jj) = ra * rad * COS( rad * pphif(ji,jj) ) * rn_e1deg
+         !                       ! e2   (meridional)
+         pe2t(ji,jj) = ra * rad * rn_e2deg
+         pe2u(ji,jj) = ra * rad * rn_e2deg
+         pe2v(ji,jj) = ra * rad * rn_e2deg
+         pe2f(ji,jj) = ra * rad * rn_e2deg
+      END_2D
       !                             ! NO reduction of grid size in some straits 
       ke1e2u_v    = 0               !    ==>> u_ & v_surfaces will be computed in dom_ghr routine
Index: NEMO/trunk/tests/ISOMIP/MY_SRC/usrdef_zgr.F90
===================================================================
--- NEMO/trunk/tests/ISOMIP/MY_SRC/usrdef_zgr.F90	(revision 12489)
+++ NEMO/trunk/tests/ISOMIP/MY_SRC/usrdef_zgr.F90	(revision 12740)
@@ -30,4 +30,6 @@
    PUBLIC   usr_def_zgr   ! called by domzgr.F90
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -132,48 +134,43 @@
             pe3vw(:,:,jk) = pe3w_1d (jk)
          END DO
-         DO jj = 1, jpj                      ! top scale factors and depth at T- and W-points
-            DO ji = 1, jpi
-               ik = k_top(ji,jj)
-               IF ( ik > 2 ) THEN
-                  ! pdeptw at the interface
-                  pdepw(ji,jj,ik  ) = MAX( zhisf(ji,jj) , pdepw(ji,jj,ik) )
-                  ! e3t in both side of the interface
-                  pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik)
-                  ! pdept in both side of the interface (from previous e3t)
-                  pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp
-                  pdept(ji,jj,ik-1) = pdepw(ji,jj,ik  ) - pe3t (ji,jj,ik  ) * 0.5_wp
-                  ! pe3w on both side of the interface
-                  pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik  )
-                  pe3w (ji,jj,ik  ) = pdept(ji,jj,ik  ) - pdept(ji,jj,ik-1)
-                  ! e3t into the ice shelf
-                  pe3t (ji,jj,ik-1) = pdepw(ji,jj,ik  ) - pdepw(ji,jj,ik-1)
-                  pe3w (ji,jj,ik-1) = pdept(ji,jj,ik-1) - pdept(ji,jj,ik-2)
-               END IF
-            END DO
-         END DO         
-         DO jj = 1, jpj                      ! bottom scale factors and depth at T- and W-points
-            DO ji = 1, jpi
-               ik = k_bot(ji,jj)
-               pdepw(ji,jj,ik+1) = MIN( zht(ji,jj) , pdepw_1d(ik+1) )
+         ! top scale factors and depth at T- and W-points
+         DO_2D_11_11
+            ik = k_top(ji,jj)
+            IF ( ik > 2 ) THEN
+               ! pdeptw at the interface
+               pdepw(ji,jj,ik  ) = MAX( zhisf(ji,jj) , pdepw(ji,jj,ik) )
+               ! e3t in both side of the interface
                pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik)
-               pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik  ) 
-               !
+               ! pdept in both side of the interface (from previous e3t)
                pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp
-               pdept(ji,jj,ik+1) = pdepw(ji,jj,ik+1) + pe3t (ji,jj,ik+1) * 0.5_wp
-               pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik)
-            END DO
-         END DO         
+               pdept(ji,jj,ik-1) = pdepw(ji,jj,ik  ) - pe3t (ji,jj,ik  ) * 0.5_wp
+               ! pe3w on both side of the interface
+               pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik  )
+               pe3w (ji,jj,ik  ) = pdept(ji,jj,ik  ) - pdept(ji,jj,ik-1)
+               ! e3t into the ice shelf
+               pe3t (ji,jj,ik-1) = pdepw(ji,jj,ik  ) - pdepw(ji,jj,ik-1)
+               pe3w (ji,jj,ik-1) = pdept(ji,jj,ik-1) - pdept(ji,jj,ik-2)
+            END IF
+         END_2D
+         ! bottom scale factors and depth at T- and W-points
+         DO_2D_11_11
+            ik = k_bot(ji,jj)
+            pdepw(ji,jj,ik+1) = MIN( zht(ji,jj) , pdepw_1d(ik+1) )
+            pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik)
+            pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik  ) 
+            !
+            pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp
+            pdept(ji,jj,ik+1) = pdepw(ji,jj,ik+1) + pe3t (ji,jj,ik+1) * 0.5_wp
+            pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik)
+         END_2D       
          !                                   ! bottom scale factors and depth at  U-, V-, UW and VW-points
          pe3u (:,:,:) = pe3t(:,:,:)
          pe3uw(:,:,:) = pe3w(:,:,:)
-         DO jk = 1, jpk                      ! Computed as the minimum of neighbooring scale factors
-            DO jj = 1, jpjm1
-               DO ji = 1, jpi
-                  pe3v (ji,jj,jk) = MIN( pe3t(ji,jj,jk), pe3t(ji,jj+1,jk) )
-                  pe3vw(ji,jj,jk) = MIN( pe3w(ji,jj,jk), pe3w(ji,jj+1,jk) )
-                  pe3f (ji,jj,jk) = pe3v(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+         DO_3D_00_00( 1, jpk )
+         !                                   ! Computed as the minimum of neighbooring scale factors
+            pe3v (ji,jj,jk) = MIN( pe3t(ji,jj,jk), pe3t(ji,jj+1,jk) )
+            pe3vw(ji,jj,jk) = MIN( pe3w(ji,jj,jk), pe3w(ji,jj+1,jk) )
+            pe3f (ji,jj,jk) = pe3v(ji,jj,jk)
+         END_3D
          CALL lbc_lnk( 'usrdef_zgr', pe3v , 'V', 1._wp )   ;   CALL lbc_lnk( 'usrdef_zgr', pe3vw, 'V', 1._wp )
          CALL lbc_lnk( 'usrdef_zgr', pe3f , 'F', 1._wp )
