Index: /NEMO/trunk/src/ABL/ablmod.F90
===================================================================
--- /NEMO/trunk/src/ABL/ablmod.F90	(revision 14214)
+++ /NEMO/trunk/src/ABL/ablmod.F90	(revision 14215)
@@ -1263,15 +1263,15 @@
       smth_b = 1._wp / 4._wp
       !
-      DO_2D( 1, 1, 1, 0 )
+      DO_2D( 1, 0, 1, 1 )
          zumsk = msk(ji,jj) * msk(ji+1,jj)
          zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji  ,jj ) ) * zumsk
       END_2D
 
-      DO_2D( 1, 0, 1, 1 )
+      DO_2D( 1, 1, 1, 0 )
          zvmsk = msk(ji,jj) * msk(ji,jj+1)
          zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji  ,jj ) ) * zvmsk
       END_2D
 
-      DO_2D( 1, 0, 0, 0 )
+      DO_2D( 0, 0, 1, 0 )
          zFY ( ji, jj  ) =   zdY ( ji, jj   )                        &
             & +  smth_a*  ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 ))   &
@@ -1279,5 +1279,5 @@
       END_2D
 
-      DO_2D( 0, 0, 1, 0 )
+      DO_2D( 1, 0, 0, 0 )
          zFX( ji, jj  ) =    zdX( ji, jj   )                         &
            &    + smth_a*(  (zdY( ji+1, jj ) - zdY( ji+1, jj-1))     &
Index: /NEMO/trunk/src/ICE/icedyn_adv_pra.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_adv_pra.F90	(revision 14214)
+++ /NEMO/trunk/src/ICE/icedyn_adv_pra.F90	(revision 14215)
@@ -582,5 +582,5 @@
          !
          ! Limitation of moments.
-         DO_2D( 1, 1, ji0, ji0 )
+         DO_2D( ji0, ji0, 1, 1 )
             !
             zpsm  = psm (ji,jj,jl) ! optimization
@@ -642,5 +642,5 @@
          END_2D
          !
-         DO_2D( 1, 0, ji0, ji0 )
+         DO_2D( ji0, ji0, 1, 0 )
             !                                !  Flux from j+1 to j when v LT 0.
             zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)
@@ -662,5 +662,5 @@
          END_2D
 
-         DO_2D( 0, 0, ji0, ji0 )
+         DO_2D( ji0, ji0, 0, 0 )
             !                                !  Readjust moments remaining in the box.
             zbt  =         zbet(ji,jj-1)
Index: /NEMO/trunk/src/ICE/icedyn_adv_umx.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_adv_umx.F90	(revision 14214)
+++ /NEMO/trunk/src/ICE/icedyn_adv_umx.F90	(revision 14215)
@@ -491,5 +491,5 @@
       IF( pamsk == 0._wp ) THEN
          DO jl = 1, jpl
-            DO_2D( 0, 0, 1, 0 )
+            DO_2D( 1, 0, 0, 0 )
                IF( ABS( pu(ji,jj) ) > epsi10 ) THEN
                   zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc    (ji,jj,jl) / pu(ji,jj)
@@ -501,5 +501,5 @@
                !
             END_2D
-            DO_2D( 1, 0, 0, 0 )
+            DO_2D( 0, 0, 1, 0 )
                IF( ABS( pv(ji,jj) ) > epsi10 ) THEN
                   zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj)
@@ -536,9 +536,9 @@
       IF( PRESENT( pua_ho ) ) THEN
          DO jl = 1, jpl
-            DO_2D( 0, 0, 1, 0 )
+            DO_2D( 1, 0, 0, 0 )
                pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl)
                pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl)
             END_2D
-            DO_2D( 1, 0, 0, 0 )
+            DO_2D( 0, 0, 1, 0 )
                pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl)
                pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl)
@@ -594,5 +594,5 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO_2D( 1, 1, 1, 0 )
+               DO_2D( 1, 0, 1, 1 )
                   pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl)
                END_2D
@@ -600,5 +600,5 @@
             !
             DO jl = 1, jpl              !-- first guess of tracer from u-flux
-               DO_2D( 1, 1, 0, 0 )
+               DO_2D( 0, 0, 1, 1 )
                   ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) )              &
                      &   + ( pu     (ji,jj   ) - pu     (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk)
@@ -609,5 +609,5 @@
             !
             DO jl = 1, jpl              !-- flux in y-direction
-               DO_2D( 1, 0, 0, 0 )
+               DO_2D( 0, 0, 1, 0 )
                   pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl)
                END_2D
@@ -617,5 +617,5 @@
             !
             DO jl = 1, jpl              !-- flux in y-direction
-               DO_2D( 1, 0, 1, 1 )
+               DO_2D( 1, 1, 1, 0 )
                   pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl)
                END_2D
@@ -623,5 +623,5 @@
             !
             DO jl = 1, jpl              !-- first guess of tracer from v-flux
-               DO_2D( 0, 0, 1, 1 )
+               DO_2D( 1, 1, 0, 0 )
                   ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) )  &
                      &   + ( pv     (ji,jj   ) - pv     (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk)
@@ -632,5 +632,5 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO_2D( 0, 0, 1, 0 )
+               DO_2D( 1, 0, 0, 0 )
                   pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl)
                END_2D
@@ -681,8 +681,8 @@
          !
          DO jl = 1, jpl
+            DO_2D( 1, 0, 1, 1 )
+               pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj  ,jl) )
+            END_2D
             DO_2D( 1, 1, 1, 0 )
-               pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj  ,jl) )
-            END_2D
-            DO_2D( 1, 0, 1, 1 )
                pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji  ,jj+1,jl) )
             END_2D
@@ -701,5 +701,5 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO_2D( 1, 1, 1, 0 )
+               DO_2D( 1, 0, 1, 1 )
                   pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
                END_2D
@@ -708,5 +708,5 @@
 
             DO jl = 1, jpl              !-- first guess of tracer from u-flux
-               DO_2D( 1, 1, 0, 0 )
+               DO_2D( 0, 0, 1, 1 )
                   ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) )              &
                      &   + ( pu    (ji,jj   ) - pu    (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk)
@@ -717,5 +717,5 @@
 
             DO jl = 1, jpl              !-- flux in y-direction
-               DO_2D( 1, 0, 0, 0 )
+               DO_2D( 0, 0, 1, 0 )
                   pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) )
                END_2D
@@ -726,5 +726,5 @@
             !
             DO jl = 1, jpl              !-- flux in y-direction
-               DO_2D( 1, 0, 1, 1 )
+               DO_2D( 1, 1, 1, 0 )
                   pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
                END_2D
@@ -733,5 +733,5 @@
             !
             DO jl = 1, jpl              !-- first guess of tracer from v-flux
-               DO_2D( 0, 0, 1, 1 )
+               DO_2D( 1, 1, 0, 0 )
                   ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) )  &
                      &   + ( pv    (ji,jj   ) - pv    (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk)
@@ -742,5 +742,5 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO_2D( 0, 0, 1, 0 )
+               DO_2D( 1, 0, 0, 0 )
                   pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) )
                END_2D
@@ -899,5 +899,5 @@
          !
          DO jl = 1, jpl
-            DO_2D( 0, 0, 1, 0 )
+            DO_2D( 1, 0, 0, 0 )
                pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
                   &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
@@ -908,5 +908,5 @@
          !
          DO jl = 1, jpl
-            DO_2D( 0, 0, 1, 0 )
+            DO_2D( 1, 0, 0, 0 )
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
@@ -918,5 +918,5 @@
          !
          DO jl = 1, jpl
-            DO_2D( 0, 0, 1, 0 )
+            DO_2D( 1, 0, 0, 0 )
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
@@ -932,5 +932,5 @@
          !
          DO jl = 1, jpl
-            DO_2D( 0, 0, 1, 0 )
+            DO_2D( 1, 0, 0, 0 )
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
@@ -946,5 +946,5 @@
          !
          DO jl = 1, jpl
-            DO_2D( 0, 0, 1, 0 )
+            DO_2D( 1, 0, 0, 0 )
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
@@ -967,5 +967,5 @@
       IF( ll_neg ) THEN
          DO jl = 1, jpl
-            DO_2D( 0, 0, 1, 0 )
+            DO_2D( 1, 0, 0, 0 )
                IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
@@ -977,5 +977,5 @@
       !                                                     !-- High order flux in i-direction  --!
       DO jl = 1, jpl
-         DO_2D( 0, 0, 1, 0 )
+         DO_2D( 1, 0, 0, 0 )
             pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl)
          END_2D
@@ -1010,5 +1010,5 @@
       !                                                     !--  Laplacian in j-direction  --!
       DO jl = 1, jpl
-         DO_2D( 1, 0, 0, 0 )         ! First derivative (gradient)
+         DO_2D( 0, 0, 1, 0 )         ! First derivative (gradient)
             ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
          END_2D
@@ -1021,5 +1021,5 @@
       !                                                     !--  BiLaplacian in j-direction  --!
       DO jl = 1, jpl
-         DO_2D( 1, 0, 0, 0 )         ! First derivative
+         DO_2D( 0, 0, 1, 0 )         ! First derivative
             ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
          END_2D
@@ -1035,5 +1035,5 @@
       CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
          DO jl = 1, jpl
-            DO_2D( 1, 0, 0, 0 )
+            DO_2D( 0, 0, 1, 0 )
                pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   &
                   &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
@@ -1043,5 +1043,5 @@
       CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
          DO jl = 1, jpl
-            DO_2D( 1, 0, 0, 0 )
+            DO_2D( 0, 0, 1, 0 )
                zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   &
@@ -1052,5 +1052,5 @@
       CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
          DO jl = 1, jpl
-            DO_2D( 1, 0, 0, 0 )
+            DO_2D( 0, 0, 1, 0 )
                zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
@@ -1065,5 +1065,5 @@
       CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
          DO jl = 1, jpl
-            DO_2D( 1, 0, 0, 0 )
+            DO_2D( 0, 0, 1, 0 )
                zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
@@ -1078,5 +1078,5 @@
       CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
          DO jl = 1, jpl
-            DO_2D( 1, 0, 0, 0 )
+            DO_2D( 0, 0, 1, 0 )
                zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
@@ -1099,5 +1099,5 @@
       IF( ll_neg ) THEN
          DO jl = 1, jpl
-            DO_2D( 1, 0, 0, 0 )
+            DO_2D( 0, 0, 1, 0 )
                IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
                   pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
@@ -1109,5 +1109,5 @@
       !                                                     !-- High order flux in j-direction  --!
       DO jl = 1, jpl
-         DO_2D( 1, 0, 0, 0 )
+         DO_2D( 0, 0, 1, 0 )
             pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl)
          END_2D
@@ -1145,8 +1145,8 @@
       ! --------------------------------------------------
       DO jl = 1, jpl
+         DO_2D( 1, 0, 0, 0 )
+            pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)
+         END_2D
          DO_2D( 0, 0, 1, 0 )
-            pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)
-         END_2D
-         DO_2D( 1, 0, 0, 0 )
             pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)
          END_2D
@@ -1254,5 +1254,5 @@
       ! ---------------------------------
       DO jl = 1, jpl
-         DO_2D( 0, 0, 1, 0 )
+         DO_2D( 1, 0, 0, 0 )
             zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) )
             zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) )
@@ -1265,5 +1265,5 @@
          END_2D
 
-         DO_2D( 1, 0, 0, 0 )
+         DO_2D( 0, 0, 1, 0 )
             zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) )
             zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) )
Index: /NEMO/trunk/src/OCE/DIA/diaptr.F90
===================================================================
--- /NEMO/trunk/src/OCE/DIA/diaptr.F90	(revision 14214)
+++ /NEMO/trunk/src/OCE/DIA/diaptr.F90	(revision 14215)
@@ -382,5 +382,5 @@
             zts(:,:,:,:) = 0._wp
 
-            DO_3D( 1, 0, 1, 1, 1, jpkm1 )
+            DO_3D( 1, 1, 1, 0, 1, jpkm1 )
                zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
                zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc
@@ -437,5 +437,5 @@
             zts(:,:,:,:) = 0._wp
 
-            DO_3D( 1, 0, 1, 1, 1, jpkm1 )
+            DO_3D( 1, 1, 1, 0, 1, jpkm1 )
                zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
                zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc  !Tracers averaged onto V grid
Index: /NEMO/trunk/src/OCE/DOM/dommsk.F90
===================================================================
--- /NEMO/trunk/src/OCE/DOM/dommsk.F90	(revision 14214)
+++ /NEMO/trunk/src/OCE/DOM/dommsk.F90	(revision 14215)
@@ -182,5 +182,5 @@
       ssfmask(:,:) = MAXVAL( fmask(:,:,:), DIM=3 )
       IF( lk_SWE ) THEN      ! Shallow Water Eq. case : redefine ssfmask
-         DO_2D( 0,0 , 0,0 )
+         DO_2D( 0,0, 0,0  )
             ssfmask(ji,jj) = MAX(  ssmask(ji,jj+1), ssmask(ji+1,jj+1),  & 
                &                   ssmask(ji,jj  ), ssmask(ji+1,jj  )   )
Index: /NEMO/trunk/src/OCE/DYN/dynldf_iso.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/dynldf_iso.F90	(revision 14214)
+++ /NEMO/trunk/src/OCE/DYN/dynldf_iso.F90	(revision 14215)
@@ -168,5 +168,5 @@
 
          IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
-            DO_2D( 0, 0, 0, 1 )
+            DO_2D( 0, 1, 0, 0 )
                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( 0, 0, 0, 1 )
+            DO_2D( 0, 1, 0, 0 )
                zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b )   &
                   &     * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj)
@@ -219,5 +219,5 @@
          ! i-flux at f-point              |   t   |
 
-         DO_2D( 0, 0, 1, 0 )
+         DO_2D( 1, 0, 0, 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( 0, 1, 1, 0 )
+            DO_2D( 1, 0, 0, 1 )
                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( 0, 1, 1, 0 )
+            DO_2D( 1, 0, 0, 1 )
                zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b )   &
                   &     * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj)
Index: /NEMO/trunk/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/dynspg_ts.F90	(revision 14214)
+++ /NEMO/trunk/src/OCE/DYN/dynspg_ts.F90	(revision 14215)
@@ -461,18 +461,18 @@
 #if defined key_qcoTest_FluxForm
             !                                ! 'key_qcoTest_FluxForm' : simple ssh average
-            DO_2D( 1, 1, 1, 0 )   ! not jpi-column
+            DO_2D( 1, 0, 1, 1 )   ! not jpi-column
                zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * (  zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj  )  ) * ssumask(ji,jj)
             END_2D
-            DO_2D( 1, 0, 1, 1 )
+            DO_2D( 1, 1, 1, 0 )
                zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * (  zsshp2_e(ji,jj) + zsshp2_e(ji  ,jj+1)  ) * ssvmask(ji,jj)
             END_2D
 #else
             !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
-            DO_2D( 1, 1, 1, 0 )   ! not jpi-column
+            DO_2D( 1, 0, 1, 1 )   ! not jpi-column
                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( 1, 0, 1, 1 )   ! not jpj-row
+            DO_2D( 1, 1, 1, 0 )   ! not jpj-row
                zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj)                        &
                     &                              * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
@@ -538,8 +538,8 @@
 #if defined key_qcoTest_FluxForm
             !                                ! 'key_qcoTest_FluxForm' : simple ssh average
+            DO_2D( 1, 0, 1, 1 )
+               zsshu_a(ji,jj) = r1_2 * (  ssha_e(ji,jj) + ssha_e(ji+1,jj  )  ) * ssumask(ji,jj)
+            END_2D
             DO_2D( 1, 1, 1, 0 )
-               zsshu_a(ji,jj) = r1_2 * (  ssha_e(ji,jj) + ssha_e(ji+1,jj  )  ) * ssumask(ji,jj)
-            END_2D
-            DO_2D( 1, 0, 1, 1 )
                zsshv_a(ji,jj) = r1_2 * (  ssha_e(ji,jj) + ssha_e(ji  ,jj+1)  ) * ssvmask(ji,jj)
             END_2D
@@ -1253,5 +1253,5 @@
       !!----------------------------------------------------------------------
       !
-      DO_2D( 1, 1, 1, 0 )   ! not jpi-column
+      DO_2D( 1, 0, 1, 1 )   ! not jpi-column
          IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj) 
          ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj)  
@@ -1261,5 +1261,5 @@
       END_2D
       !
-      DO_2D( 1, 0, 1, 1 )   ! not jpj-row
+      DO_2D( 1, 1, 1, 0 )   ! not jpj-row
          IF ( phV(ji,jj) > 0._wp ) THEN   ;   pVmsk(ji,jj) = pTmsk(ji,jj  )
          ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1)  
Index: /NEMO/trunk/src/OCE/ISF/isftbl.F90
===================================================================
--- /NEMO/trunk/src/OCE/ISF/isftbl.F90	(revision 14214)
+++ /NEMO/trunk/src/OCE/ISF/isftbl.F90	(revision 14215)
@@ -78,5 +78,5 @@
          ! compute tbl property at T point
          pvarout(1,:) = 0._wp
-         DO_2D( 1, 1, 0, 1 )
+         DO_2D( 0, 1, 1, 1 )
             pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji-1,jj))
          END_2D
@@ -99,5 +99,5 @@
          ! pvarout is an averaging of wet point
          pvarout(:,1) = 0._wp
-         DO_2D( 0, 1, 1, 1 )
+         DO_2D( 1, 1, 0, 1 )
             pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji,jj-1))
          END_2D
Index: /NEMO/trunk/src/OCE/SBC/geo2ocean.F90
===================================================================
--- /NEMO/trunk/src/OCE/SBC/geo2ocean.F90	(revision 14214)
+++ /NEMO/trunk/src/OCE/SBC/geo2ocean.F90	(revision 14215)
@@ -160,5 +160,5 @@
       ! (computation done on the north stereographic polar plane)
       !
-      DO_2D( 0, 0, 0, 1 )
+      DO_2D( 0, 1, 0, 0 )
          !                  
          zlam = plamt(ji,jj)     ! north pole direction & modulous (at t-point)
@@ -249,5 +249,5 @@
       ! =============== !
 
-      DO_2D( 0, 0, 0, 1 )
+      DO_2D( 0, 1, 0, 0 )
          IF( MOD( ABS( plamv(ji,jj) - plamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
             gsint(ji,jj) = 0.
Index: /NEMO/trunk/src/OCE/SBC/sbcice_cice.F90
===================================================================
--- /NEMO/trunk/src/OCE/SBC/sbcice_cice.F90	(revision 14214)
+++ /NEMO/trunk/src/OCE/SBC/sbcice_cice.F90	(revision 14215)
@@ -312,5 +312,5 @@
 ! x comp of wind stress (CI_1)
 ! U point to F point
-         DO_2D( 1, 0, 1, 1 )
+         DO_2D( 1, 1, 1, 0 )
             ztmp(ji,jj) = 0.5 * (  fr_iu(ji,jj) * utau(ji,jj)      &
                                  + fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
@@ -320,5 +320,5 @@
 ! y comp of wind stress (CI_2)
 ! V point to F point
-         DO_2D( 1, 1, 1, 0 )
+         DO_2D( 1, 0, 1, 1 )
             ztmp(ji,jj) = 0.5 * (  fr_iv(ji,jj) * vtau(ji,jj)      &
                                  + fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
@@ -437,5 +437,5 @@
 ! x comp and y comp of surface ocean current
 ! U point to F point
-      DO_2D( 1, 0, 1, 1 )
+      DO_2D( 1, 1, 1, 0 )
          ztmp(ji,jj)=0.5*(ssu_m(ji,jj)+ssu_m(ji,jj+1))*fmask(ji,jj,1)
       END_2D
@@ -443,5 +443,5 @@
 
 ! V point to F point
-      DO_2D( 1, 1, 1, 0 )
+      DO_2D( 1, 0, 1, 1 )
          ztmp(ji,jj)=0.5*(ssv_m(ji,jj)+ssv_m(ji+1,jj))*fmask(ji,jj,1)
       END_2D
@@ -513,5 +513,5 @@
 ! F point to V point
 
-      DO_2D( 1, 0, 0, 0 )
+      DO_2D( 0, 0, 1, 0 )
          ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1)
       END_2D
Index: /NEMO/trunk/src/OCE/TRA/traadv_qck.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/traadv_qck.F90	(revision 14214)
+++ /NEMO/trunk/src/OCE/TRA/traadv_qck.F90	(revision 14215)
@@ -145,5 +145,5 @@
          !
 !!gm why not using a SHIFT instruction...
-         DO_3D( 0, 0, nn_hls-1, nn_hls-1, 1, jpkm1 )     !--- Computation of the ustream and downstream value of the tracer and the mask
+         DO_3D( nn_hls-1, nn_hls-1, 0, 0, 1, jpkm1 )     !--- Computation of the ustream and downstream value of the tracer and the mask
             zfc(ji,jj,jk) = pt(ji-1,jj,jk,jn,Kbb)        ! Upstream   in the x-direction for the tracer
             zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb)        ! Downstream in the x-direction for the tracer
@@ -154,10 +154,10 @@
          ! Horizontal advective fluxes
          ! ---------------------------
-         DO_3D( 0, 0, nn_hls-1, 0, 1, jpkm1 )
+         DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 )
             zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
             zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
          END_3D
          !
-         DO_3D( 0, 0, nn_hls-1, 0, 1, jpkm1 )
+         DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 )
             zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
             zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
@@ -173,5 +173,5 @@
          !
          ! Mask at the T-points in the x-direction (mask=0 or mask=1)
-         DO_3D( 0, 0, nn_hls-1, nn_hls-1, 1, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, 0, 0, 1, jpkm1 )
             zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
          END_3D
@@ -180,5 +180,5 @@
          !
          ! Tracer flux on the x-direction
-         DO_3D( 0, 0, 1, 0, 1, jpkm1 )
+         DO_3D( 1, 0, 0, 0, 1, jpkm1 )
             zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
             !--- If the second ustream point is a land point
@@ -232,5 +232,5 @@
             !
             !--- Computation of the ustream and downstream value of the tracer and the mask
-            DO_2D( nn_hls-1, nn_hls-1, 0, 0 )
+            DO_2D( 0, 0, nn_hls-1, nn_hls-1 )
                ! Upstream in the x-direction for the tracer
                zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb)
@@ -245,10 +245,10 @@
          ! ---------------------------
          !
-         DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, nn_hls-1, 0, 1, jpkm1 )
             zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
             zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
          END_3D
          !
-         DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, nn_hls-1, 0, 1, jpkm1 )
             zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
             zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
@@ -265,5 +265,5 @@
          !
          ! Mask at the T-points in the x-direction (mask=0 or mask=1)
-         DO_3D( nn_hls-1, nn_hls-1, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, nn_hls-1, nn_hls-1, 1, jpkm1 )
             zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
          END_3D
@@ -271,5 +271,5 @@
          !
          ! Tracer flux on the x-direction
-         DO_3D( 1, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 1, 0, 1, jpkm1 )
             zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
             !--- If the second ustream point is a land point
Index: /NEMO/trunk/src/OCE/TRA/trabbl.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/trabbl.F90	(revision 14214)
+++ /NEMO/trunk/src/OCE/TRA/trabbl.F90	(revision 14215)
@@ -248,5 +248,5 @@
       DO jn = 1, kjpt                                            ! tracer loop
          !                                                       ! ===========
-         DO_2D( isj, 0, isi, 0 )            ! CAUTION start from i=1 to update i=2 when cyclic east-west
+         DO_2D( isi, 0, isj, 0 )            ! CAUTION start from i=1 to update i=2 when cyclic east-west
             IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
                ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
Index: /NEMO/trunk/src/OCE/TRA/traldf_lap_blp.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/traldf_lap_blp.F90	(revision 14214)
+++ /NEMO/trunk/src/OCE/TRA/traldf_lap_blp.F90	(revision 14215)
@@ -158,5 +158,5 @@
          ENDIF
          !
-         DO_3D( isj, iej, isi, iei, 1, jpkm1 )            !== Second derivative (divergence) added to the general tracer trends  ==!
+         DO_3D( isi, iei, isj, iej, 1, jpkm1 )            !== Second derivative (divergence) added to the general tracer trends  ==!
             pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     &
                &                                      +    ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )   &
Index: /NEMO/trunk/src/OCE/TRA/traldf_triad.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/traldf_triad.F90	(revision 14214)
+++ /NEMO/trunk/src/OCE/TRA/traldf_triad.F90	(revision 14215)
@@ -387,5 +387,5 @@
          !                                !==  add the vertical 33 flux  ==!
          IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
-            DO_3D( 1, 0, 0, 0, 2, jpkm1 )
+            DO_3D( 0, 0, 1, 0, 2, jpkm1 )
                ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)   &
                   &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
@@ -395,10 +395,10 @@
             SELECT CASE( kpass )
             CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
-               DO_3D( 1, 0, 0, 0, 2, jpkm1 )
+               DO_3D( 0, 0, 1, 0, 2, jpkm1 )
                   ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)             &
                      &                            * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
                END_3D
             CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp.
-               DO_3D( 1, 0, 0, 0, 2, jpkm1 )
+               DO_3D( 0, 0, 1, 0, 2, jpkm1 )
                   ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)                      &
                      &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
Index: /NEMO/trunk/src/OCE/TRA/tranpc.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/tranpc.F90	(revision 14214)
+++ /NEMO/trunk/src/OCE/TRA/tranpc.F90	(revision 14215)
@@ -112,5 +112,5 @@
          IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; ENDIF
          !
-         DO_2D( isj, iej, isi, iei )                        ! interior column only
+         DO_2D( isi, iei, isj, iej )                        ! interior column only
             !
             IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points
Index: /NEMO/trunk/src/OCE/TRA/traqsr.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/traqsr.F90	(revision 14214)
+++ /NEMO/trunk/src/OCE/TRA/traqsr.F90	(revision 14215)
@@ -151,5 +151,5 @@
          ELSE                                           ! No restart or Euler forward at 1st time step
             z1_2 = 1._wp
-            DO_3D( isj, iej, isi, iei, 1, jpk )
+            DO_3D( isi, iei, isj, iej, 1, jpk )
                qsr_hc_b(ji,jj,jk) = 0._wp
             END_3D
@@ -157,5 +157,5 @@
       ELSE                             !==  Swap of qsr heat content  ==!
          z1_2 = 0.5_wp
-         DO_3D( isj, iej, isi, iei, 1, jpk )
+         DO_3D( isi, iei, isj, iej, 1, jpk )
             qsr_hc_b(ji,jj,jk) = qsr_hc(ji,jj,jk)
          END_3D
@@ -168,5 +168,5 @@
       CASE( np_BIO )                   !==  bio-model fluxes  ==!
          !
-         DO_3D( isj, iej, isi, iei, 1, nksr )
+         DO_3D( isi, iei, isj, iej, 1, nksr )
             qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) )
          END_3D
@@ -190,5 +190,5 @@
             ! most expensive calculations)
             !
-            DO_2D( isj, iej, isi, iei )
+            DO_2D( isi, iei, isj, iej )
                        ! zlogc = log(zchl)
                zlogc = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) )
@@ -209,5 +209,5 @@
 
 !
-            DO_3D( isj, iej, isi, iei, 1, nksr + 1 )
+            DO_3D( isi, iei, isj, iej, 1, nksr + 1 )
                ! zchl    = ALOG( ze0(ji,jj) )
                zlogc = ze0(ji,jj)
@@ -239,5 +239,5 @@
          !
          zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
-         DO_2D( isj, iej, isi, iei )
+         DO_2D( isi, iei, isj, iej )
             ze0(ji,jj) = rn_abs * qsr(ji,jj)
             ze1(ji,jj) = zcoef  * qsr(ji,jj)
@@ -250,5 +250,5 @@
          !
          !                                    !* interior equi-partition in R-G-B depending on vertical profile of Chl
-         DO_3D( isj, iej, isi, iei, 2, nksr + 1 )
+         DO_3D( isi, iei, isj, iej, 2, nksr + 1 )
             ze3t = e3t(ji,jj,jk-1,Kmm)
             irgb = NINT( ztmp3d(ji,jj,jk) )
@@ -264,5 +264,5 @@
          END_3D
          !
-         DO_3D( isj, iej, isi, iei, 1, nksr )          !* now qsr induced heat content
+         DO_3D( isi, iei, isj, iej, 1, nksr )          !* now qsr induced heat content
             qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) )
          END_3D
@@ -274,5 +274,5 @@
          zz0 =        rn_abs   * r1_rho0_rcp      ! surface equi-partition in 2-bands
          zz1 = ( 1. - rn_abs ) * r1_rho0_rcp
-         DO_3D( isj, iej, isi, iei, 1, nksr )          !* now qsr induced heat content
+         DO_3D( isi, iei, isj, iej, 1, nksr )          !* now qsr induced heat content
             zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r )
             zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r )
@@ -292,5 +292,5 @@
       !
       ! sea-ice: store the 1st ocean level attenuation coefficient
-      DO_2D( isj, iej, isi, iei )
+      DO_2D( isi, iei, isj, iej )
          IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) )
          ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
Index: /NEMO/trunk/src/OCE/TRA/trasbc.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/trasbc.F90	(revision 14214)
+++ /NEMO/trunk/src/OCE/TRA/trasbc.F90	(revision 14215)
@@ -105,5 +105,5 @@
 !!gm  This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist)
       IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
-         DO_2D( isj, iej, isi, iei )
+         DO_2D( isi, iei, isj, iej )
             qns(ji,jj) = qns(ji,jj) + qsr(ji,jj)      ! total heat flux in qns
             qsr(ji,jj) = 0._wp                        ! qsr set to zero
@@ -126,5 +126,5 @@
          ELSE                                             ! No restart or restart not found: Euler forward time stepping
             zfact = 1._wp
-            DO_2D( isj, iej, isi, iei )
+            DO_2D( isi, iei, isj, iej )
                sbc_tsc(ji,jj,:) = 0._wp
                sbc_tsc_b(ji,jj,:) = 0._wp
@@ -133,15 +133,15 @@
       ELSE                                !* other time-steps: swap of forcing fields
          zfact = 0.5_wp
-         DO_2D( isj, iej, isi, iei )
+         DO_2D( isi, iei, isj, iej )
             sbc_tsc_b(ji,jj,:) = sbc_tsc(ji,jj,:)
          END_2D
       ENDIF
       !                             !==  Now sbc tracer content fields  ==!
-      DO_2D( isj, iej, isi, iei )
+      DO_2D( isi, iei, isj, iej )
          sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj)   ! non solar heat flux
          sbc_tsc(ji,jj,jp_sal) = r1_rho0     * sfx(ji,jj)   ! salt flux due to freezing/melting
       END_2D
       IF( ln_linssh ) THEN                !* linear free surface
-         DO_2D( isj, iej, isi, iei )                    !==>> add concentration/dilution effect due to constant volume cell
+         DO_2D( isi, iei, isj, iej )                    !==>> add concentration/dilution effect due to constant volume cell
             sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm)
             sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm)
Index: /NEMO/trunk/src/OCE/ZDF/zdfosm.F90
===================================================================
--- /NEMO/trunk/src/OCE/ZDF/zdfosm.F90	(revision 14214)
+++ /NEMO/trunk/src/OCE/ZDF/zdfosm.F90	(revision 14215)
@@ -2347,5 +2347,5 @@
       ! calculate horizontal gradients at u & v points
 
-      DO_2D( 0, 0, 1, 0 )
+      DO_2D( 1, 0, 0, 0 )
          zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj)
          zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj)
@@ -2355,5 +2355,5 @@
       END_2D
 
-      DO_2D( 1, 0, 0, 0 )
+      DO_2D( 0, 0, 1, 0 )
          zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
          zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
@@ -2366,8 +2366,8 @@
       CALL eos_rab(ztsm_midv, zmld_midv, zabv, Kmm)
 
-      DO_2D( 0, 0, 1, 0 )
+      DO_2D( 1, 0, 0, 0 )
          dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal))
       END_2D
-      DO_2D( 1, 0, 0, 0 )
+      DO_2D( 0, 0, 1, 0 )
          dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal))
       END_2D
Index: /NEMO/trunk/src/OCE/do_loop_substitute.h90
===================================================================
--- /NEMO/trunk/src/OCE/do_loop_substitute.h90	(revision 14214)
+++ /NEMO/trunk/src/OCE/do_loop_substitute.h90	(revision 14215)
@@ -2,5 +2,5 @@
 ! These comments are not intended to be retained during preprocessing; i.e. do not define "show_comments"
 !!----------------------------------------------------------------------
-!! NEMO/OCE 4.0 , NEMO Consortium (2018)
+!! NEMO/OCE 4.x , NEMO Consortium (2020)
 !! Software governed by the CeCILL license (see ./LICENSE)
 !!----------------------------------------------------------------------
@@ -25,9 +25,9 @@
 ! Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)
 !
-! The macro naming convention takes the form: DO_2D( B, T, L, R) where:
+! The macro naming convention takes the form: DO_2D( L, R, B, T) where:
+!   L is the Left   offset from the PE's inner domain;
+!   R is the Right  offset from the PE's inner domain
 !   B is the Bottom offset from the PE's inner domain;
 !   T is the Top    offset from the PE's inner domain;
-!   L is the Left   offset from the PE's inner domain;
-!   R is the Right  offset from the PE's inner domain
 !
 ! So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:
@@ -42,5 +42,5 @@
 ! with:
 !
-!   DO_2D( 0, 1, 1, 0 )
+!   DO_2D( 1, 0, 0, 1 )
 !      .
 !      .
@@ -59,5 +59,5 @@
 #endif
 
-#define DO_2D(B, T, L, R) DO jj = ntsj-(B), ntej+(T)   ;   DO ji = ntsi-(L), ntei+(R)
+#define DO_2D(L, R, B, T) DO jj = ntsj-(B), ntej+(T)   ;   DO ji = ntsi-(L), ntei+(R)
 #define A1Di(H) ntsi-H:ntei+H
 #define A1Dj(H) ntsj-H:ntej+H
@@ -70,7 +70,7 @@
 #define KJPT  :
 
-#define DO_3D(B, T, L, R, ks, ke) DO jk = ks, ke   ;   DO_2D(B, T, L, R)
+#define DO_3D(L, R, B, T, ks, ke) DO jk = ks, ke   ;   DO_2D(L, R, B, T)
 
-#define DO_3DS(B, T, L, R, ks, ke, ki) DO jk = ks, ke, ki   ;   DO_2D(B, T, L, R)
+#define DO_3DS(L, R, B, T, ks, ke, ki) DO jk = ks, ke, ki   ;   DO_2D(L, R, B, T)
 
 #define END_2D   END DO   ;   END DO
Index: /NEMO/trunk/src/TOP/TRP/trcsbc.F90
===================================================================
--- /NEMO/trunk/src/TOP/TRP/trcsbc.F90	(revision 14214)
+++ /NEMO/trunk/src/TOP/TRP/trcsbc.F90	(revision 14215)
@@ -122,5 +122,5 @@
          !
          DO jn = 1, jptra
-            DO_2D( 0, 1, 0, 0 )
+            DO_2D( 0, 0, 0, 1 )
                sbc_trc(ji,jj,jn) = zsfx(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm)
             END_2D
@@ -130,5 +130,5 @@
          !
          DO jn = 1, jptra
-            DO_2D( 0, 1, 0, 0 )
+            DO_2D( 0, 0, 0, 1 )
                sbc_trc(ji,jj,jn) = ( zsfx(ji,jj) + fmmflx(ji,jj) ) * r1_rho0 * ptr(ji,jj,1,jn,Kmm)
             END_2D
@@ -138,5 +138,5 @@
          !
          DO jn = 1, jptra
-            DO_2D( 0, 1, 0, 0 )
+            DO_2D( 0, 0, 0, 1 )
                zse3t = 1. / e3t(ji,jj,1,Kmm)
                ! tracer flux at the ice/ocean interface (tracer/m2/s)
@@ -162,5 +162,5 @@
          IF( l_trdtrc )   ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs)  ! save trends
          !
-         DO_2D( 0, 1, 0, 0 )
+         DO_2D( 0, 0, 0, 1 )
             zse3t = zfact / e3t(ji,jj,1,Kmm)
             ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + ( sbc_trc_b(ji,jj,jn) + sbc_trc(ji,jj,jn) ) * zse3t
Index: /NEMO/trunk/src/TOP/trcbc.F90
===================================================================
--- /NEMO/trunk/src/TOP/trcbc.F90	(revision 14214)
+++ /NEMO/trunk/src/TOP/trcbc.F90	(revision 14215)
@@ -415,5 +415,5 @@
          ! Remove river dilution for tracers with absent river load
          IF( ln_rnf_ctl .AND. .NOT.ln_trc_cbc(jn) ) THEN
-            DO_2D( 0, 1, 0, 0 )
+            DO_2D( 0, 0, 0, 1 )
                DO jk = 1, nk_rnf(ji,jj)
                   zrnf = (rnf(ji,jj) + rnf_b(ji,jj)) * 0.5_wp * r1_rho0 / h_rnf(ji,jj)
@@ -429,5 +429,5 @@
             jl = n_trc_indsbc(jn)
             sf_trcsbc(jl)%fnow(:,:,1) = MAX( rtrn, sf_trcsbc(jl)%fnow(:,:,1) ) ! avoid nedgative value due to interpolation
-            DO_2D( 0, 1, 0, 0 )
+            DO_2D( 0, 0, 0, 1 )
                zfact = 1. / ( e3t(ji,jj,1,Kmm) * rn_sbc_time )
                ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + rf_trsfac(jl) * sf_trcsbc(jl)%fnow(ji,jj,1) * zfact
@@ -439,5 +439,5 @@
             IF( l_offline )   rn_rfact = 1._wp
             jl = n_trc_indcbc(jn)
-            DO_2D( 0, 1, 0, 0 )
+            DO_2D( 0, 0, 0, 1 )
                DO jk = 1, nk_rnf(ji,jj)
                   zfact = rn_rfact / ( e1e2t(ji,jj) * h_rnf(ji,jj) * rn_cbc_time ) * tmask(ji,jj,1)
