Index: NEMO/trunk/src/OCE/DYN/divhor.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/divhor.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/divhor.F90	(revision 13295)
@@ -77,5 +77,5 @@
       ENDIF
       !
-      DO_3D_00_00( 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm)      &
             &               - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm)      &
Index: NEMO/trunk/src/OCE/DYN/dynadv_cen2.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynadv_cen2.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/dynadv_cen2.F90	(revision 13295)
@@ -72,5 +72,5 @@
          zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm)
          zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj  ,jk,Kmm) )
             zfv_f(ji  ,jj  ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) )
@@ -78,5 +78,5 @@
             zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji  ,jj+1,jk,Kmm) )
          END_2D
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
                &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj)   &
@@ -98,10 +98,10 @@
       !                             !==  Vertical advection  ==!
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          zfu_uw(ji,jj,jpk) = 0._wp   ;   zfv_vw(ji,jj,jpk) = 0._wp
          zfu_uw(ji,jj, 1 ) = 0._wp   ;   zfv_vw(ji,jj, 1 ) = 0._wp
       END_2D
       IF( ln_linssh ) THEN                ! linear free surface: advection through the surface
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
             zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
@@ -109,13 +109,13 @@
       ENDIF
       DO jk = 2, jpkm1                    ! interior advective fluxes
-         DO_2D_01_01
+         DO_2D( 0, 1, 0, 1 )
             zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
          END_2D
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj  ,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) )
             zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji  ,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) )
          END_2D
       END DO
-      DO_3D_00_00( 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj)   &
             &                                      / e3u(ji,jj,jk,Kmm)
Index: NEMO/trunk/src/OCE/DYN/dynadv_ubs.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynadv_ubs.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/dynadv_ubs.F90	(revision 13295)
@@ -108,5 +108,5 @@
          zfv(:,:,jk) = e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
          !            
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zlu_uu(ji,jj,jk,1) = ( puu (ji+1,jj  ,jk,Kbb) - 2.*puu (ji,jj,jk,Kbb) + puu (ji-1,jj  ,jk,Kbb) ) * umask(ji,jj,jk)
             zlv_vv(ji,jj,jk,1) = ( pvv (ji  ,jj+1,jk,Kbb) - 2.*pvv (ji,jj,jk,Kbb) + pvv (ji  ,jj-1,jk,Kbb) ) * vmask(ji,jj,jk)
@@ -136,5 +136,5 @@
          zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
          !
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj  ,jk,Kmm) )
             zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji  ,jj+1,jk,Kmm) )
@@ -168,5 +168,5 @@
                &                * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) - gamma1 * zl_v )
          END_2D
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
                &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj)   &
@@ -187,5 +187,5 @@
       !                                      !  Vertical advection  !
       !                                      ! ==================== !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          zfu_uw(ji,jj,jpk) = 0._wp
          zfv_vw(ji,jj,jpk) = 0._wp
@@ -194,5 +194,5 @@
       END_2D
       IF( ln_linssh ) THEN                         ! constant volume : advection through the surface
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
             zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
@@ -200,13 +200,13 @@
       ENDIF
       DO jk = 2, jpkm1                          ! interior fluxes
-         DO_2D_01_01
+         DO_2D( 0, 1, 0, 1 )
             zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
          END_2D
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) )
             zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) )
          END_2D
       END DO
-      DO_3D_00_00( 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj)   &
             &                                       / e3u(ji,jj,jk,Kmm)
Index: NEMO/trunk/src/OCE/DYN/dynatf.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynatf.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/dynatf.F90	(revision 13295)
@@ -197,5 +197,5 @@
          IF( ln_linssh ) THEN             ! Fixed volume !
             !                             ! =============!
-            DO_3D_11_11( 1, jpkm1 )
+            DO_3D( 1, 1, 1, 1, 1, jpkm1 )
                puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
                pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + rn_atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
@@ -233,5 +233,5 @@
                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3u(:,:,:,Kmm), 'U' )
                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3v(:,:,:,Kmm), 'V' )
-               DO_3D_11_11( 1, jpkm1 )
+               DO_3D( 1, 1, 1, 1, 1, jpkm1 )
                   puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
                   pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + rn_atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
@@ -244,5 +244,5 @@
                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3u_f, 'U' )
                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3v_f, 'V' )
-               DO_3D_11_11( 1, jpkm1 )
+               DO_3D( 1, 1, 1, 1, 1, jpkm1 )
                   zue3a = pe3u(ji,jj,jk,Kaa) * puu(ji,jj,jk,Kaa)
                   zve3a = pe3v(ji,jj,jk,Kaa) * pvv(ji,jj,jk,Kaa)
Index: NEMO/trunk/src/OCE/DYN/dynatf_qco.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynatf_qco.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/dynatf_qco.F90	(revision 13295)
@@ -140,5 +140,5 @@
          IF( ln_linssh ) THEN             ! Fixed volume !
             !                             ! =============!
-            DO_3D_11_11( 1, jpkm1 )
+            DO_3D( 1, 1, 1, 1, 1, jpkm1 )
                puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
                pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + rn_atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
@@ -150,5 +150,5 @@
             IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
                ! Before filtered scale factor at (u/v)-points
-               DO_3D_11_11( 1, jpkm1 )
+               DO_3D( 1, 1, 1, 1, 1, jpkm1 )
                   puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
                   pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + rn_atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
@@ -157,5 +157,5 @@
             ELSE                          ! Asselin filter applied on thickness weighted velocity
                !
-               DO_3D_11_11( 1, jpkm1 )
+               DO_3D( 1, 1, 1, 1, 1, jpkm1 )
                   zue3a = ( 1._wp + r3u(ji,jj,Kaa) * umask(ji,jj,jk) ) * puu(ji,jj,jk,Kaa)
                   zve3a = ( 1._wp + r3v(ji,jj,Kaa) * vmask(ji,jj,jk) ) * pvv(ji,jj,jk,Kaa)
Index: NEMO/trunk/src/OCE/DYN/dynhpg.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynhpg.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/dynhpg.F90	(revision 13295)
@@ -257,5 +257,5 @@
 
       ! Surface value
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm)
          ! hydrostatic pressure gradient
@@ -269,5 +269,5 @@
       !
       ! interior value (2=<jk=<jpkm1)
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm)
          ! hydrostatic pressure gradient
@@ -319,5 +319,5 @@
 
       !  Surface value (also valid in partial step case)
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm)
          ! hydrostatic pressure gradient
@@ -330,5 +330,5 @@
 
       ! interior value (2=<jk=<jpkm1)
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm)
          ! hydrostatic pressure gradient
@@ -346,5 +346,5 @@
 
       ! partial steps correction at the last level  (use zgru & zgrv computed in zpshde.F90)
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          iku = mbku(ji,jj)
          ikv = mbkv(ji,jj)
@@ -411,5 +411,5 @@
       !
       IF( ln_wd_il ) THEN
-        DO_2D_00_00
+        DO_2D( 0, 0, 0, 0 )
           ll_tmp1 = MIN(  ssh(ji,jj,Kmm)               ,  ssh(ji+1,jj,Kmm) ) >                &
                &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
@@ -452,5 +452,5 @@
 
       ! Surface value
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ! hydrostatic pressure gradient along s-surfaces
          zhpi(ji,jj,1) =   &
@@ -481,5 +481,5 @@
 
       ! interior value (2=<jk=<jpkm1)
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          ! hydrostatic pressure gradient along s-surfaces
          zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   &
@@ -563,5 +563,5 @@
 !===== Compute surface value ===================================================== 
 !==================================================================================
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ikt    = mikt(ji,jj)
          iktp1i = mikt(ji+1,jj)
@@ -592,5 +592,5 @@
 !==================================================================================
       ! interior value (2=<jk=<jpkm1)
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          ! hydrostatic pressure gradient along s-surfaces
          zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
@@ -643,5 +643,5 @@
       IF( ln_wd_il ) THEN
          ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) )
-        DO_2D_00_00
+        DO_2D( 0, 0, 0, 0 )
           ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji+1,jj,Kmm) ) >                &
                &    MAX( -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) .AND.            &
@@ -699,5 +699,5 @@
 !!bug gm   Not a true bug, but... dzz=e3w  for dzx, dzy verify what it is really
 
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          drhoz(ji,jj,jk) = rhd    (ji  ,jj  ,jk) - rhd    (ji,jj,jk-1)
          dzz  (ji,jj,jk) = gde3w(ji  ,jj  ,jk) - gde3w(ji,jj,jk-1)
@@ -716,5 +716,5 @@
 !!bug  gm  idem for drhox, drhoy et ji=jpi and jj=jpj
 
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          cffw = 2._wp * drhoz(ji  ,jj  ,jk) * drhoz(ji,jj,jk-1)
 
@@ -784,5 +784,5 @@
 !          true if gde3w(:,:,:) is really defined as the sum of the e3w scale factors as, it seems to me, it should be
 
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          rho_k(ji,jj,1) = -grav * ( e3w(ji,jj,1,Kmm) - gde3w(ji,jj,1) )               &
             &                   * (  rhd(ji,jj,1)                                     &
@@ -795,5 +795,5 @@
 !!bug gm    : optimisation: 1/10 and 1/12 the division should be done before the loop
 
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
 
          rho_k(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj,jk) + rhd    (ji,jj,jk-1) )                                   &
@@ -830,5 +830,5 @@
       !  Surface value
       ! ---------------
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj)
          zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj)
@@ -845,5 +845,5 @@
       !  interior value   (2=<jk=<jpkm1)
       ! ----------------
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          ! hydrostatic pressure gradient along s-surfaces
          zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)                                &
@@ -911,5 +911,5 @@
       IF( ln_wd_il ) THEN
          ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) )
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
           ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji+1,jj,Kmm) ) >                &
                &    MAX( -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) .AND.            &
@@ -960,5 +960,5 @@
 
       ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
        jk = mbkt(ji,jj)
        IF(     jk <=  1   ) THEN   ;   zrhh(ji,jj,    :   ) = 0._wp
@@ -973,9 +973,9 @@
 
       ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)"
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) * znad
       END_2D
 
-      DO_3D_11_11( 2, jpk )
+      DO_3D( 1, 1, 1, 1, 2, jpk )
          zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w(ji,jj,jk,Kmm)
       END_3D
@@ -990,5 +990,5 @@
 
       ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)"
-      DO_2D_01_01
+      DO_2D( 0, 1, 0, 1 )
        zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1),  &
           &                                              csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm)
@@ -999,5 +999,5 @@
 
       ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)"
-      DO_3D_01_01( 2, jpkm1 )
+      DO_3D( 0, 1, 0, 1, 2, jpkm1 )
       zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) +                                  &
          &             integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk),   &
@@ -1009,5 +1009,5 @@
 
       ! Prepare zsshu_n and zsshv_n
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
 !!gm BUG ?    if it is ssh at u- & v-point then it should be:
 !          zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * &
@@ -1024,20 +1024,20 @@
       CALL lbc_lnk_multi ('dynhpg', zsshu_n, 'U', 1.0_wp, zsshv_n, 'V', 1.0_wp )
 
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
        zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad) 
        zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad)
       END_2D
 
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
       zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm)
       zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm)
       END_3D
 
-      DO_3D_00_00( 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
       zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm)
       zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm)
       END_3D
 
-      DO_3D_00_00( 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
       zu(ji,jj,jk) = MIN(  zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  )
       zu(ji,jj,jk) = MAX(  zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  )
@@ -1047,5 +1047,5 @@
 
 
-      DO_3D_00_00( 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
       zpwes = 0._wp; zpwed = 0._wp
       zpnss = 0._wp; zpnsd = 0._wp
Index: NEMO/trunk/src/OCE/DYN/dynkeg.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynkeg.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/dynkeg.F90	(revision 13295)
@@ -101,5 +101,5 @@
       !
       CASE ( nkeg_C2 )                          !--  Standard scheme  --!
-         DO_3D_01_01( 1, jpkm1 )
+         DO_3D( 0, 1, 0, 1, 1, jpkm1 )
             zu =    puu(ji-1,jj  ,jk,Kmm) * puu(ji-1,jj  ,jk,Kmm)   &
                &  + puu(ji  ,jj  ,jk,Kmm) * puu(ji  ,jj  ,jk,Kmm)
@@ -109,5 +109,5 @@
          END_3D
       CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             zu = 8._wp * ( puu(ji-1,jj  ,jk,Kmm) * puu(ji-1,jj  ,jk,Kmm)    &
                &         + puu(ji  ,jj  ,jk,Kmm) * puu(ji  ,jj  ,jk,Kmm) )  &
@@ -125,5 +125,5 @@
       END SELECT 
       !
-      DO_3D_00_00( 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
          pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
Index: NEMO/trunk/src/OCE/DYN/dynldf_iso.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynldf_iso.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/dynldf_iso.F90	(revision 13295)
@@ -128,5 +128,5 @@
       IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
          !
-         DO_3D_00_00( 1, jpk )
+         DO_3D( 0, 0, 0, 0, 1, jpk )
             uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
             vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
@@ -168,5 +168,5 @@
 
          IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
-            DO_2D_00_01
+            DO_2D( 0, 0, 0, 1 )
                zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj)   &
                   &    * MIN( e3u(ji  ,jj,jk,Kmm),                &
@@ -183,5 +183,5 @@
             END_2D
          ELSE                   ! other coordinate system (zco or sco) : e3t
-            DO_2D_00_01
+            DO_2D( 0, 0, 0, 1 )
                zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b )   &
                   &     * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj)
@@ -199,5 +199,5 @@
 
          ! j-flux at f-point
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b )   &
                &     * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj)
@@ -219,5 +219,5 @@
          ! i-flux at f-point              |   t   |
 
-         DO_2D_00_10
+         DO_2D( 0, 0, 1, 0 )
             zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b )   &
                &     * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj)
@@ -235,5 +235,5 @@
          ! j-flux at t-point
          IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
-            DO_2D_01_10
+            DO_2D( 0, 1, 1, 0 )
                zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj)   &
                   &     * MIN( e3v(ji,jj  ,jk,Kmm),                 &
@@ -250,5 +250,5 @@
             END_2D
          ELSE                   ! other coordinate system (zco or sco) : e3t
-            DO_2D_01_10
+            DO_2D( 0, 1, 1, 0 )
                zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b )   &
                   &     * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj)
@@ -268,5 +268,5 @@
          ! Second derivative (divergence) and add to the general trend
          ! -----------------------------------------------------------
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (  ziut(ji+1,jj) - ziut(ji,jj  )    &
                &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj)   &
Index: NEMO/trunk/src/OCE/DYN/dynldf_lap_blp.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynldf_lap_blp.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/dynldf_lap_blp.F90	(revision 13295)
@@ -73,5 +73,5 @@
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
-         DO_2D_01_01
+         DO_2D( 0, 1, 0, 1 )
             !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
             zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &   ! ahmf already * by fmask
@@ -84,5 +84,5 @@
          END_2D
          !
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * (    &    ! * by umask is mandatory for dyn_ldf_blp use
                &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
Index: NEMO/trunk/src/OCE/DYN/dynspg.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynspg.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/dynspg.F90	(revision 13295)
@@ -95,5 +95,5 @@
          .OR.  ln_ice_embd ) THEN                                            ! embedded sea-ice
          !
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             spgu(ji,jj) = 0._wp
             spgv(ji,jj) = 0._wp
@@ -102,5 +102,5 @@
          IF( ln_apr_dyn .AND. .NOT.ln_dynspg_ts ) THEN   !==  Atmospheric pressure gradient (added later in time-split case) ==!
             zg_2 = grav * 0.5
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                spgu(ji,jj) = spgu(ji,jj) + zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    &
                   &                                + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
@@ -117,5 +117,5 @@
             CALL upd_tide(zt0step, Kmm)
             !
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
                spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
@@ -124,5 +124,5 @@
             IF (ln_scal_load) THEN
                zld = rn_scal_load * grav
-               DO_2D_00_00
+               DO_2D( 0, 0, 0, 0 )
                   spgu(ji,jj) = spgu(ji,jj) + zld * ( pssh(ji+1,jj,Kmm) - pssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)
                   spgv(ji,jj) = spgv(ji,jj) + zld * ( pssh(ji,jj+1,Kmm) - pssh(ji,jj,Kmm) ) * r1_e2v(ji,jj)
@@ -136,5 +136,5 @@
             zgrho0r     = - grav * r1_rho0
             zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrho0r
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
                spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
@@ -143,5 +143,5 @@
          ENDIF
          !
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj)
             pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + spgv(ji,jj)
Index: NEMO/trunk/src/OCE/DYN/dynspg_exp.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynspg_exp.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/dynspg_exp.F90	(revision 13295)
@@ -74,10 +74,10 @@
       IF( ln_linssh ) THEN          !* linear free surface : add the surface pressure gradient trend
          !
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             spgu(ji,jj) = - grav * ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)
             spgv(ji,jj) = - grav * ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) * r1_e2v(ji,jj)
          END_2D
          !
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj)
             pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + spgv(ji,jj)
Index: NEMO/trunk/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynspg_ts.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/dynspg_ts.F90	(revision 13295)
@@ -264,5 +264,5 @@
          IF( ln_wd_il ) THEN                       ! W/D : limiter applied to spgspg
             CALL wad_spg( pssh(:,:,Kmm), zcpx, zcpy )          ! Calculating W/D gravity filters, zcpx and zcpy
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj  ,Kmm) - pssh(ji  ,jj ,Kmm) )   &
                   &                          * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth
@@ -271,5 +271,5 @@
             END_2D
          ELSE                                      ! now suface pressure gradient
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  pssh(ji+1,jj  ,Kmm) - pssh(ji  ,jj  ,Kmm)  ) * r1_e1u(ji,jj)
                zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  pssh(ji  ,jj+1,Kmm) - pssh(ji  ,jj  ,Kmm)  ) * r1_e2v(ji,jj) 
@@ -279,5 +279,5 @@
       ENDIF
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
           zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
           zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
@@ -291,5 +291,5 @@
       IF( ln_apr_dyn ) THEN
          IF( ln_bt_fw ) THEN                          ! FORWARD integration: use kt+1/2 pressure (NOW+1/2)
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_frc(ji,jj) = zu_frc(ji,jj) + grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
                zv_frc(ji,jj) = zv_frc(ji,jj) + grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
@@ -297,5 +297,5 @@
          ELSE                                         ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW)
             zztmp = grav * r1_2
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)  &
                     &                                   + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
@@ -309,5 +309,5 @@
       !                                   !  ----------------------------------  !
       IF( ln_bt_fw ) THEN                        ! Add wind forcing
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zu_frc(ji,jj) =  zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm)
             zv_frc(ji,jj) =  zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm)
@@ -315,5 +315,5 @@
       ELSE
          zztmp = r1_rho0 * r1_2
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zu_frc(ji,jj) =  zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm)
             zv_frc(ji,jj) =  zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm)
@@ -475,10 +475,10 @@
             !
             !                          ! ocean u- and v-depth at mid-step   (separate DO-loops remove the need of a lbc_lnk)
-            DO_2D_11_10
+            DO_2D( 1, 1, 1, 0 )
                zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj)                        &
                     &                              * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
                     &                                 + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj)
             END_2D
-            DO_2D_10_11
+            DO_2D( 1, 0, 1, 1 )
                zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj)                        &
                     &                              * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
@@ -515,5 +515,5 @@
          !--        ssh    =  ssh   - delta_t' * [ frc + div( flux      ) ]      --!
          !-------------------------------------------------------------------------!
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zhdiv = (   zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1)   ) * r1_e1e2t(ji,jj)
             ssha_e(ji,jj) = (  sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv )  ) * ssmask(ji,jj)
@@ -541,5 +541,5 @@
          ! Sea Surface Height at u-,v-points (vvl case only)
          IF( .NOT.ln_linssh ) THEN                                
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
                   &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
@@ -561,5 +561,5 @@
          !                             ! Surface pressure gradient
          zldg = ( 1._wp - rn_scal_load ) * grav    ! local factor
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
             zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
@@ -579,5 +579,5 @@
          ! Add tidal astronomical forcing if defined
          IF ( ln_tide .AND. ln_tide_pot ) THEN
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
                zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
@@ -588,5 +588,5 @@
 !jth do implicitly instead
          IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
                zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
@@ -606,5 +606,5 @@
          !------------------------------------------------------------------------------------------------------------------------!
          IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
                          &     + rDt_e * (                   zu_spg(ji,jj)   &
@@ -621,5 +621,5 @@
             !
          ELSE                           !* Flux form
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                !                    ! hu_e, hv_e hold depth at jn,  zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
                !                    ! backward interpolated depth used in spg terms at jn+1/2
@@ -645,5 +645,5 @@
 !jth implicit bottom friction:
          IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                   ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj))
                   va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj))
@@ -712,5 +712,5 @@
       IF (ln_bt_fw) THEN
          IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )
                zun_save = un_adv(ji,jj)
                zvn_save = vn_adv(ji,jj)
@@ -743,5 +743,5 @@
       ELSE
          ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj) &
                &              * ( e1e2t(ji  ,jj) * pssh(ji  ,jj,Kaa)      &
@@ -975,5 +975,5 @@
       ! Max courant number for ext. grav. waves
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
          zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
@@ -1100,5 +1100,5 @@
          SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point
          CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    &
                     &           ht(ji  ,jj  ) + ht(ji+1,jj  )   ) * 0.25_wp  
@@ -1106,5 +1106,5 @@
             END_2D
          CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) =             (  ht  (ji  ,jj+1) + ht  (ji+1,jj+1)      &
                     &                    + ht  (ji  ,jj  ) + ht  (ji+1,jj  )  )   &
@@ -1117,5 +1117,5 @@
          !
          ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
-         DO_2D_01_01
+         DO_2D( 0, 1, 0, 1 )
             ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
             ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
@@ -1126,5 +1126,5 @@
       CASE( np_EET )                  != EEN scheme using e3t energy conserving scheme
          ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
-         DO_2D_01_01
+         DO_2D( 0, 1, 0, 1 )
             z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
             ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht
@@ -1159,5 +1159,5 @@
             !
             !zhf(:,:) = hbatf(:,:)
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          &
                     &            + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      &
@@ -1178,5 +1178,5 @@
          CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
          ! JC: TBC. hf should be greater than 0 
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj)
          END_2D
@@ -1201,5 +1201,5 @@
       SELECT CASE( nvor_scheme )
       CASE( np_ENT )                ! enstrophy conserving scheme (f-point)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) )
             z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) )
@@ -1214,5 +1214,5 @@
          !         
       CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
             zy2 = ( zhV(ji,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
@@ -1225,5 +1225,5 @@
          !
       CASE( np_ENS )                ! enstrophy conserving scheme (f-point)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zy1 =   r1_8 * ( zhV(ji  ,jj-1) + zhV(ji+1,jj-1) &
               &            + zhV(ji  ,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
@@ -1235,5 +1235,5 @@
          !
       CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)         
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zhV(ji  ,jj  ) &
              &                                         + ftnw(ji+1,jj) * zhV(ji+1,jj  ) &
@@ -1269,5 +1269,5 @@
       !
       IF( ln_wd_dl_rmp ) THEN     
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF    ( pssh(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN 
                !           IF    ( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin2 ) THEN 
@@ -1280,5 +1280,5 @@
          END_2D
       ELSE  
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF ( pssh(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN   ;   ptmsk(ji,jj) = 1._wp
             ELSE                                                 ;   ptmsk(ji,jj) = 0._wp
@@ -1308,5 +1308,5 @@
       !!----------------------------------------------------------------------
       !
-      DO_2D_11_10
+      DO_2D( 1, 1, 1, 0 )
          IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj) 
          ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj)  
@@ -1316,5 +1316,5 @@
       END_2D
       !
-      DO_2D_10_11
+      DO_2D( 1, 0, 1, 1 )
          IF ( phV(ji,jj) > 0._wp ) THEN   ;   pVmsk(ji,jj) = pTmsk(ji,jj  )
          ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1)  
@@ -1338,5 +1338,5 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
       !!----------------------------------------------------------------------
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ll_tmp1 = MIN(  pshn(ji,jj)               ,  pshn(ji+1,jj) ) >                &
               &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
@@ -1407,10 +1407,10 @@
       IF( ln_isfcav ) THEN          ! top+bottom friction (ocean cavities)
          
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
             pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
          END_2D
       ELSE                          ! bottom friction only
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
             pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
@@ -1422,5 +1422,5 @@
       IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities
          
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ikbu = mbku(ji,jj)       
             ikbv = mbkv(ji,jj)    
@@ -1430,5 +1430,5 @@
       ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities
          
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ikbu = mbku(ji,jj)       
             ikbv = mbkv(ji,jj)    
@@ -1440,5 +1440,5 @@
       IF( ln_wd_il ) THEN      ! W/D : use the "clipped" bottom friction   !!gm   explain WHY, please !
          zztmp = -1._wp / rDt_e
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 & 
                  &                              r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp  )
@@ -1448,5 +1448,5 @@
       ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
          
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj)
             pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj)
@@ -1460,5 +1460,5 @@
          IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity
             
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                iktu = miku(ji,jj)
                iktv = mikv(ji,jj)
@@ -1468,5 +1468,5 @@
          ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity
             
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                iktu = miku(ji,jj)
                iktv = mikv(ji,jj)
@@ -1478,5 +1478,5 @@
          !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
          
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj)
             pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj)
Index: NEMO/trunk/src/OCE/DYN/dynvor.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynvor.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/dynvor.F90	(revision 13295)
@@ -231,10 +231,10 @@
       CASE ( np_RVO )                           !* relative vorticity
          DO jk = 1, jpkm1                                 ! Horizontal slab
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
                   &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
             END_2D
             IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
-               DO_2D_10_10
+               DO_2D( 1, 0, 1, 0 )
                   zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
                END_2D
@@ -246,10 +246,10 @@
       CASE ( np_CRV )                           !* Coriolis + relative vorticity
          DO jk = 1, jpkm1                                 ! Horizontal slab
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj,jk) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   &
                   &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj)
             END_2D
             IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
-               DO_2D_10_10
+               DO_2D( 1, 0, 1, 0 )
                   zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
                END_2D
@@ -269,5 +269,5 @@
             zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm)
          CASE ( np_RVO )                           !* relative vorticity
-            DO_2D_01_01
+            DO_2D( 0, 1, 0, 1 )
                zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   &
                   &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) &
@@ -275,5 +275,5 @@
             END_2D
          CASE ( np_MET )                           !* metric term
-            DO_2D_01_01
+            DO_2D( 0, 1, 0, 1 )
                zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   &
                   &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) &
@@ -281,5 +281,5 @@
             END_2D
          CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO_2D_01_01
+            DO_2D( 0, 1, 0, 1 )
                zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    &
                   &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) &
@@ -287,5 +287,5 @@
             END_2D
          CASE ( np_CME )                           !* Coriolis + metric
-            DO_2D_01_01
+            DO_2D( 0, 1, 0, 1 )
                zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           &
                     &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  &
@@ -298,5 +298,5 @@
          !
          !                                   !==  compute and add the vorticity term trend  =!
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                    &
                &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   &
@@ -358,20 +358,20 @@
             zwz(:,:) = ff_f(:,:) 
          CASE ( np_RVO )                           !* relative vorticity
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
                   &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
             END_2D
          CASE ( np_MET )                           !* metric term
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
                   &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
             END_2D
          CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
                   &                        - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
             END_2D
          CASE ( np_CME )                           !* Coriolis + metric
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
                   &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
@@ -382,5 +382,5 @@
          !
          IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
             END_2D
@@ -396,5 +396,5 @@
          ENDIF
          !                                   !==  compute and add the vorticity term trend  =!
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
             zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
@@ -454,20 +454,20 @@
             zwz(:,:) = ff_f(:,:) 
          CASE ( np_RVO )                           !* relative vorticity
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
                   &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
             END_2D
          CASE ( np_MET )                           !* metric term
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
                   &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
             END_2D
          CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
                   &                        - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
             END_2D
          CASE ( np_CME )                           !* Coriolis + metric
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
                   &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
@@ -478,5 +478,5 @@
          !
          IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==!
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
             END_2D
@@ -492,5 +492,5 @@
          ENDIF
          !                                   !==  compute and add the vorticity term trend  =!
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
                &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
@@ -550,5 +550,5 @@
          SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point
          CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
                   &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
@@ -560,5 +560,5 @@
             END_2D
          CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
                   &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
@@ -575,19 +575,19 @@
          SELECT CASE( kvor )                 !==  vorticity considered  ==!
          CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj)
             END_2D
          CASE ( np_RVO )                           !* relative vorticity
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj,jk) = ( e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
                   &            - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
             END_2D
          CASE ( np_MET )                           !* metric term
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj,jk) = (   ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
                   &              - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
             END_2D
          CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
                   &                              - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  )   &
@@ -595,5 +595,5 @@
             END_2D
          CASE ( np_CME )                           !* Coriolis + metric
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
                   &                            - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
@@ -604,5 +604,5 @@
          !
          IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
             END_2D
@@ -635,5 +635,5 @@
             END DO
          END DO
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
                &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
@@ -695,9 +695,9 @@
          SELECT CASE( kvor )                 !==  vorticity considered  ==!
          CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj,jk) = ff_f(ji,jj)
             END_2D
          CASE ( np_RVO )                           !* relative vorticity
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj,jk) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
                   &             - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
@@ -705,10 +705,10 @@
             END_2D
          CASE ( np_MET )                           !* metric term
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj,jk) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
                   &          - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
             END_2D
          CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
                   &                              - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
@@ -716,5 +716,5 @@
             END_2D
          CASE ( np_CME )                           !* Coriolis + metric
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
                   &                        - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
@@ -725,5 +725,5 @@
          !
          IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
             END_2D
@@ -758,5 +758,5 @@
             END DO
          END DO
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
                &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
@@ -818,5 +818,5 @@
       IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
       IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
-         DO_3D_10_10( 1, jpk )
+         DO_3D( 1, 0, 1, 0, 1, jpk )
             IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              &
                & + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
@@ -857,5 +857,5 @@
          CASE( np_ENT )                      !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
             ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj  ) ) * 0.5_wp
                dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji  ,jj-1) ) * 0.5_wp
@@ -865,5 +865,5 @@
          CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f)
             ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
                dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
Index: NEMO/trunk/src/OCE/DYN/dynzad.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynzad.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/dynzad.F90	(revision 13295)
@@ -78,8 +78,8 @@
       
       DO jk = 2, jpkm1              ! Vertical momentum advection at level w and u- and v- vertical
-         DO_2D_01_01
+         DO_2D( 0, 1, 0, 1 )
             zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
          END_2D
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( puu(ji,jj,jk-1,Kmm) - puu(ji,jj,jk,Kmm) )
             zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( pvv(ji,jj,jk-1,Kmm) - pvv(ji,jj,jk,Kmm) )
@@ -88,5 +88,5 @@
       !
       ! Surface and bottom advective fluxes set to zero
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          zwuw(ji,jj, 1 ) = 0._wp
          zwvw(ji,jj, 1 ) = 0._wp
@@ -95,5 +95,5 @@
       END_2D
       !
-      DO_3D_00_00( 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj)   &
             &                                      / e3u(ji,jj,jk,Kmm)
Index: NEMO/trunk/src/OCE/DYN/dynzdf.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynzdf.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/dynzdf.F90	(revision 13295)
@@ -107,10 +107,10 @@
       !                    ! time stepping except vertical diffusion
       IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kbb) + rDt * puu(ji,jj,jk,Krhs) ) * umask(ji,jj,jk)
             pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kbb) + rDt * pvv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk)
          END_3D
       ELSE                                      ! applied on thickness weighted velocity
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             puu(ji,jj,jk,Kaa) = (         e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb )  &
                &                  + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs)  ) &
@@ -127,9 +127,9 @@
       !     G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa)
       IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
-         DO_3D_00_00( 1, jpkm1 )      ! remove barotropic velocities
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )      ! remove barotropic velocities
             puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - uu_b(ji,jj,Kaa) ) * umask(ji,jj,jk)
             pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk)
          END_3D
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
             ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
@@ -142,5 +142,5 @@
          END_2D
          IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                iku = miku(ji,jj)         ! top ocean level at u- and v-points 
                ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
@@ -162,5 +162,5 @@
          SELECT CASE( nldf_dyn )
          CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &
                   &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
@@ -176,5 +176,5 @@
             END_3D
          CASE DEFAULT               ! iso-level lateral mixing
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &    ! after scale factor at U-point
                   &             + r_vvl   * e3u(ji,jj,jk,Kaa)
@@ -190,5 +190,5 @@
             END_3D
          END SELECT
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zwi(ji,jj,1) = 0._wp
             ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm)    &
@@ -203,5 +203,5 @@
          SELECT CASE( nldf_dyn )
          CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &
                   &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
@@ -215,5 +215,5 @@
             END_3D
          CASE DEFAULT               ! iso-level lateral mixing
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &
                   &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
@@ -227,5 +227,5 @@
             END_3D
          END SELECT
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zwi(ji,jj,1) = 0._wp
             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
@@ -241,5 +241,5 @@
       !
       IF ( ln_drgimp ) THEN      ! implicit bottom friction
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
             ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    &
@@ -248,5 +248,5 @@
          END_2D
          IF ( ln_isfcav ) THEN   ! top friction (always implicit)
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
                iku = miku(ji,jj)       ! ocean top level at u- and v-points 
@@ -273,9 +273,9 @@
       !-----------------------------------------------------------------------
       !
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
       END_3D
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm)    &
             &             + r_vvl   * e3u(ji,jj,1,Kaa) 
@@ -283,12 +283,12 @@
             &                                      / ( ze3ua * rho0 ) * umask(ji,jj,1) 
       END_2D
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * puu(ji,jj,jk-1,Kaa)
       END_3D
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
       END_2D
-      DO_3DS_00_00( jpk-2, 1, -1 )
+      DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
          puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
       END_3D
@@ -301,5 +301,5 @@
          SELECT CASE( nldf_dyn )
          CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
                   &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
@@ -315,5 +315,5 @@
             END_3D
          CASE DEFAULT               ! iso-level lateral mixing
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
                   &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
@@ -329,5 +329,5 @@
             END_3D
          END SELECT
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zwi(ji,jj,1) = 0._wp
             ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm)    &
@@ -342,5 +342,5 @@
          SELECT CASE( nldf_dyn )
          CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
                   &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
@@ -354,5 +354,5 @@
             END_3D
          CASE DEFAULT               ! iso-level lateral mixing
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
                   &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
@@ -366,5 +366,5 @@
             END_3D
          END SELECT
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zwi(ji,jj,1) = 0._wp
             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
@@ -379,5 +379,5 @@
       !
       IF( ln_drgimp ) THEN
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
             ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
@@ -386,5 +386,5 @@
          END_2D
          IF ( ln_isfcav ) THEN
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
                ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
@@ -410,9 +410,9 @@
       !-----------------------------------------------------------------------
       !
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
       END_3D
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm)    &
             &             + r_vvl   * e3v(ji,jj,1,Kaa) 
@@ -420,12 +420,12 @@
             &                                      / ( ze3va * rho0 ) * vmask(ji,jj,1) 
       END_2D
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pvv(ji,jj,jk-1,Kaa)
       END_3D
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
       END_2D
-      DO_3DS_00_00( jpk-2, 1, -1 )
+      DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
          pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
       END_3D
Index: NEMO/trunk/src/OCE/DYN/sshwzv.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/sshwzv.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/sshwzv.F90	(revision 13295)
@@ -178,5 +178,5 @@
             ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
             ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
             END_2D
@@ -358,5 +358,5 @@
       zdt = 2._wp * rn_Dt                            ! 2*rn_Dt and not rDt (for restartability)
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
             Cu_adv(ji,jj,jk) =   zdt *                                                         &
@@ -375,5 +375,5 @@
          END_3D
       ELSE
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
             Cu_adv(ji,jj,jk) =   zdt *                                                      &
@@ -393,5 +393,5 @@
       !
       IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
-         DO_3DS_11_11( jpkm1, 2, -1 )
+         DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )
             !
             zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
Index: NEMO/trunk/src/OCE/DYN/wet_dry.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/wet_dry.F90	(revision 13289)
+++ NEMO/trunk/src/OCE/DYN/wet_dry.F90	(revision 13295)
@@ -174,5 +174,5 @@
       !
       wdmask(:,:) = 1._wp
-      DO_2D_01_01
+      DO_2D( 0, 1, 0, 1 )
          !
          IF( tmask(ji,jj,1)        < 0.5_wp )   CYCLE    ! we don't care about land cells
@@ -198,5 +198,5 @@
       wdramp(:,:) = min((ht_0(:,:) + psshb1(:,:) - rn_wdmin1)/(rn_wdmin0 - rn_wdmin1),1.0_wp)
       !jth assume don't need a lbc_lnk here
-      DO_2D_10_10
+      DO_2D( 1, 0, 1, 0 )
          wdrampu(ji,jj) = MIN( wdramp(ji,jj) , wdramp(ji+1,jj) )
          wdrampv(ji,jj) = MIN( wdramp(ji,jj) , wdramp(ji,jj+1) )
@@ -211,5 +211,5 @@
          jflag = 0     ! flag indicating if any further iterations are needed
          !
-         DO_2D_01_01
+         DO_2D( 0, 1, 0, 1 )
             IF( tmask(ji, jj, 1) < 0.5_wp )   CYCLE 
             IF( ht_0(ji,jj)      > zdepwd )   CYCLE
@@ -307,5 +307,5 @@
       zwdlmtv(:,:) = 1._wp
       !
-      DO_2D_01_01
+      DO_2D( 0, 1, 0, 1 )
          !
          IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE   ! we don't care about land cells
@@ -333,5 +333,5 @@
          jflag = 0     ! flag indicating if any further iterations are needed
          !
-         DO_2D_01_01
+         DO_2D( 0, 1, 0, 1 )
             !
             IF( tmask(ji, jj, 1 ) < 0.5_wp )   CYCLE 
