Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv_umx.F90	(revision 12252)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv_umx.F90	(revision 12340)
@@ -52,4 +52,5 @@
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
@@ -107,20 +108,18 @@
       ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- !
       DO jl = 1, jpl
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), &
-                  &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), &
-                  &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), &
-                  &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) )
-               zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), &
-                  &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), &
-                  &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), &
-                  &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) )
-               zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), &
-                  &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), &
-                  &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), &
-                  &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) )
-            END DO
-         END DO
+         DO_2D_00_00
+            zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), &
+               &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), &
+               &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), &
+               &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) )
+            zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), &
+               &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), &
+               &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), &
+               &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) )
+            zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), &
+               &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), &
+               &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), &
+               &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) )
+         END_2D
       END DO
       CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. )
@@ -152,17 +151,15 @@
       !
       ! --- define velocity for advection: u*grad(H) --- !
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
-            IF    ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
-            ELSEIF( pu_ice(ji,jj)                   >  0._wp ) THEN   ;   zcu_box(ji,jj) = pu_ice(ji-1,jj)
-            ELSE                                                      ;   zcu_box(ji,jj) = pu_ice(ji  ,jj)
-            ENDIF
-
-            IF    ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
-            ELSEIF( pv_ice(ji,jj)                   >  0._wp ) THEN   ;   zcv_box(ji,jj) = pv_ice(ji,jj-1)
-            ELSE                                                      ;   zcv_box(ji,jj) = pv_ice(ji,jj  )
-            ENDIF
-         END DO
-      END DO
+      DO_2D_00_00
+         IF    ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
+         ELSEIF( pu_ice(ji,jj)                   >  0._wp ) THEN   ;   zcu_box(ji,jj) = pu_ice(ji-1,jj)
+         ELSE                                                      ;   zcu_box(ji,jj) = pu_ice(ji  ,jj)
+         ENDIF
+
+         IF    ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
+         ELSEIF( pv_ice(ji,jj)                   >  0._wp ) THEN   ;   zcv_box(ji,jj) = pv_ice(ji,jj-1)
+         ELSE                                                      ;   zcv_box(ji,jj) = pv_ice(ji,jj  )
+         ENDIF
+      END_2D
 
       !---------------!
@@ -187,14 +184,12 @@
             IF( .NOT. ALLOCATED(jmsk_small) )   ALLOCATE( jmsk_small(jpi,jpj,jpl) ) 
             DO jl = 1, jpl
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpim1
-                     zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) )
-                     IF( zvi_cen < epsi06) THEN   ;   imsk_small(ji,jj,jl) = 0
-                     ELSE                         ;   imsk_small(ji,jj,jl) = 1   ;   ENDIF
-                     zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) )
-                     IF( zvi_cen < epsi06) THEN   ;   jmsk_small(ji,jj,jl) = 0
-                     ELSE                         ;   jmsk_small(ji,jj,jl) = 1   ;   ENDIF
-                  END DO
-               END DO
+               DO_2D_10_10
+                  zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) )
+                  IF( zvi_cen < epsi06) THEN   ;   imsk_small(ji,jj,jl) = 0
+                  ELSE                         ;   imsk_small(ji,jj,jl) = 1   ;   ENDIF
+                  zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) )
+                  IF( zvi_cen < epsi06) THEN   ;   jmsk_small(ji,jj,jl) = 0
+                  ELSE                         ;   jmsk_small(ji,jj,jl) = 1   ;   ENDIF
+               END_2D
             END DO
          ENDIF
@@ -338,10 +333,8 @@
          !== Open water area ==!
          zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & 
-                  &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
-            END DO
-         END DO
+         DO_2D_00_00
+            pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & 
+               &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
+         END_2D
          CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1. )
          !
@@ -449,23 +442,21 @@
       IF( pamsk == 0._wp ) THEN
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  IF( ABS( pu(ji,jj) ) > epsi10 ) THEN
-                     zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc    (ji,jj,jl) / pu(ji,jj)
-                     zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj)
-                  ELSE
-                     zfu_ho (ji,jj,jl) = 0._wp
-                     zfu_ups(ji,jj,jl) = 0._wp
-                  ENDIF
-                  !
-                  IF( ABS( pv(ji,jj) ) > epsi10 ) THEN
-                     zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj)
-                     zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj)
-                  ELSE
-                     zfv_ho (ji,jj,jl) = 0._wp  
-                     zfv_ups(ji,jj,jl) = 0._wp  
-                  ENDIF
-               END DO
-            END DO
+            DO_2D_10_10
+               IF( ABS( pu(ji,jj) ) > epsi10 ) THEN
+                  zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc    (ji,jj,jl) / pu(ji,jj)
+                  zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj)
+               ELSE
+                  zfu_ho (ji,jj,jl) = 0._wp
+                  zfu_ups(ji,jj,jl) = 0._wp
+               ENDIF
+               !
+               IF( ABS( pv(ji,jj) ) > epsi10 ) THEN
+                  zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj)
+                  zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj)
+               ELSE
+                  zfv_ho (ji,jj,jl) = 0._wp  
+                  zfv_ups(ji,jj,jl) = 0._wp  
+               ENDIF
+            END_2D
          END DO
 
@@ -473,11 +464,9 @@
          ! thus we calculate the upstream solution and apply a limiter again
          DO jl = 1, jpl
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) )
-                  !
-                  zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)
-               END DO
-            END DO
+            DO_2D_00_00
+               ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) )
+               !
+               zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)
+            END_2D
          END DO
          CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T',  1. )
@@ -496,10 +485,8 @@
       IF( PRESENT( pua_ho ) ) THEN
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl)
-                  pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl)
-              END DO
-            END DO
+            DO_2D_10_10
+               pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl)
+               pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl)
+            END_2D
          END DO
       ENDIF
@@ -508,11 +495,9 @@
       ! ---------------------------------
       DO jl = 1, jpl
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1 
-               ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) )  
-               !
-               ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)               
-            END DO
-         END DO
+         DO_2D_00_00
+            ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) )  
+            !
+            ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)               
+         END_2D
       END DO
       CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1. )
@@ -544,72 +529,58 @@
          !
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
+            DO_2D_10_10
+               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)
+               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
+         END DO
+         !
+      ELSE                              !** alternate directions **!
+         !
+         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
+            !
+            DO jl = 1, jpl              !-- flux in x-direction
+               DO_2D_10_10
                   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
+            END DO
+            !
+            DO jl = 1, jpl              !-- first guess of tracer from u-flux
+               DO_2D_00_00
+                  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)
+                  !
+                  zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+               END_2D
+            END DO
+            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
+            !
+            DO jl = 1, jpl              !-- flux in y-direction
+               DO_2D_10_10
+                  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
+            END DO
+            !
+         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
+            !
+            DO jl = 1, jpl              !-- flux in y-direction
+               DO_2D_10_10
                   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 DO
-            END DO
-         END DO
-         !
-      ELSE                              !** alternate directions **!
-         !
-         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
+               END_2D
+            END DO
+            !
+            DO jl = 1, jpl              !-- first guess of tracer from v-flux
+               DO_2D_00_00
+                  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)
+                  !
+                  zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+               END_2D
+            END DO
+            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1
-                     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 DO
-               END DO
-            END DO
-            !
-            DO jl = 1, jpl              !-- first guess of tracer from u-flux
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1
-                     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)
-                     !
-                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
-                  END DO
-               END DO
-            END DO
-            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
-            !
-            DO jl = 1, jpl              !-- flux in y-direction
-               DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1
-                     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 DO
-               END DO
-            END DO
-            !
-         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
-            !
-            DO jl = 1, jpl              !-- flux in y-direction
-               DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1
-                     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 DO
-               END DO
-            END DO
-            !
-            DO jl = 1, jpl              !-- first guess of tracer from v-flux
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1
-                     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)
-                     !
-                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
-                  END DO
-               END DO
-            END DO
-            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
-            !
-            DO jl = 1, jpl              !-- flux in x-direction
-               DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1
-                     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 DO
-               END DO
+               DO_2D_10_10
+                  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
             END DO
             !
@@ -619,14 +590,12 @@
       !
       DO jl = 1, jpl                    !-- after tracer with upstream scheme
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               ztra = - (   pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj  ,jl)   &
-                  &       + pfv_ups(ji,jj,jl) - pfv_ups(ji  ,jj-1,jl) ) &
-                  &   + (   pu     (ji,jj   ) - pu     (ji-1,jj     )   &
-                  &       + pv     (ji,jj   ) - pv     (ji  ,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk)
-               !
-               pt_ups(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
-            END DO
-         END DO
+         DO_2D_00_00
+            ztra = - (   pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj  ,jl)   &
+               &       + pfv_ups(ji,jj,jl) - pfv_ups(ji  ,jj-1,jl) ) &
+               &   + (   pu     (ji,jj   ) - pu     (ji-1,jj     )   &
+               &       + pv     (ji,jj   ) - pv     (ji  ,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk)
+            !
+            pt_ups(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+         END_2D
       END DO
       CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1. )
@@ -660,10 +629,8 @@
          !
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj  ,jl) )
-                  pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji  ,jj+1,jl) )
-               END DO
-            END DO
+            DO_2D_10_10
+               pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj  ,jl) )
+               pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji  ,jj+1,jl) )
+            END_2D
          END DO
          !
@@ -680,30 +647,24 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1
-                     pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
-                  END DO
-               END DO
+               DO_2D_10_10
+                  pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
+               END_2D
             END DO
             IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
 
             DO jl = 1, jpl              !-- first guess of tracer from u-flux
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1
-                     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)
-                     !
-                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
-                  END DO
-               END DO
+               DO_2D_00_00
+                  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)
+                  !
+                  zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+               END_2D
             END DO
             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
 
             DO jl = 1, jpl              !-- flux in y-direction
-               DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1
-                     pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) )
-                  END DO
-               END DO
+               DO_2D_10_10
+                  pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) )
+               END_2D
             END DO
             IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
@@ -712,30 +673,24 @@
             !
             DO jl = 1, jpl              !-- flux in y-direction
-               DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1
-                     pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
-                  END DO
-               END DO
+               DO_2D_10_10
+                  pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
+               END_2D
             END DO
             IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
             !
             DO jl = 1, jpl              !-- first guess of tracer from v-flux
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1
-                     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)
-                     !
-                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
-                  END DO
-               END DO
+               DO_2D_00_00
+                  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)
+                  !
+                  zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+               END_2D
             END DO
             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1
-                     pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) )
-                  END DO
-               END DO
+               DO_2D_10_10
+                  pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) )
+               END_2D
             END DO
             IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
@@ -783,12 +738,10 @@
          !                                                        !--  advective form update in zpt  --!
          DO jl = 1, jpl
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  zpt(ji,jj,jl) = ( pt(ji,jj,jl) - (  pubox(ji,jj   ) * ( zt_u(ji,jj,jl) - zt_u(ji-1,jj,jl) ) * r1_e1t  (ji,jj) &
-                     &                              + pt   (ji,jj,jl) * ( pu  (ji,jj   ) - pu  (ji-1,jj   ) ) * r1_e1e2t(ji,jj) &
-                     &                                                                                        * pamsk           &
-                     &                             ) * pdt ) * tmask(ji,jj,1)
-               END DO
-            END DO
+            DO_2D_00_00
+               zpt(ji,jj,jl) = ( pt(ji,jj,jl) - (  pubox(ji,jj   ) * ( zt_u(ji,jj,jl) - zt_u(ji-1,jj,jl) ) * r1_e1t  (ji,jj) &
+                  &                              + pt   (ji,jj,jl) * ( pu  (ji,jj   ) - pu  (ji-1,jj   ) ) * r1_e1e2t(ji,jj) &
+                  &                                                                                        * pamsk           &
+                  &                             ) * pdt ) * tmask(ji,jj,1)
+            END_2D
          END DO
          CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
@@ -812,12 +765,10 @@
          !                                                        !--  advective form update in zpt  --!
          DO jl = 1, jpl
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  zpt(ji,jj,jl) = ( pt(ji,jj,jl) - (  pvbox(ji,jj   ) * ( zt_v(ji,jj,jl) - zt_v(ji,jj-1,jl) ) * r1_e2t  (ji,jj) &
-                     &                              + pt   (ji,jj,jl) * ( pv  (ji,jj   ) - pv  (ji,jj-1   ) ) * r1_e1e2t(ji,jj) &
-                     &                                                                                        * pamsk           &
-                     &                             ) * pdt ) * tmask(ji,jj,1) 
-               END DO
-            END DO
+            DO_2D_00_00
+               zpt(ji,jj,jl) = ( pt(ji,jj,jl) - (  pvbox(ji,jj   ) * ( zt_v(ji,jj,jl) - zt_v(ji,jj-1,jl) ) * r1_e2t  (ji,jj) &
+                  &                              + pt   (ji,jj,jl) * ( pv  (ji,jj   ) - pv  (ji,jj-1   ) ) * r1_e1e2t(ji,jj) &
+                  &                                                                                        * pamsk           &
+                  &                             ) * pdt ) * tmask(ji,jj,1) 
+            END_2D
          END DO
          CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
@@ -896,10 +847,8 @@
          !        
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  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) ) )
-               END DO
-            END DO
+            DO_2D_10_10
+               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) ) )
+            END_2D
          END DO
          !
@@ -907,11 +856,9 @@
          !
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  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)   &
-                     &                                                            - zcu   * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
-               END DO
-            END DO
+            DO_2D_10_10
+               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)   &
+                  &                                                            - zcu   * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
+            END_2D
          END DO
          !  
@@ -919,15 +866,13 @@
          !
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
-                  zdx2 = e1u(ji,jj) * e1u(ji,jj)
+            DO_2D_10_10
+               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
+               zdx2 = e1u(ji,jj) * e1u(ji,jj)
 !!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj)
-                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     &
-                     &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) &
-                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) *    (                      ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
-                     &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
-               END DO
-            END DO
+               pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     &
+                  &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) &
+                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) *    (                      ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
+                  &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
+            END_2D
          END DO
          !
@@ -935,15 +880,13 @@
          !
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
-                  zdx2 = e1u(ji,jj) * e1u(ji,jj)
+            DO_2D_10_10
+               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
+               zdx2 = e1u(ji,jj) * e1u(ji,jj)
 !!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj)
-                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     &
-                     &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) &
-                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) *    (                      ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
-                     &                                                   - 0.5_wp * zcu   * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
-               END DO
-            END DO
+               pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     &
+                  &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) &
+                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) *    (                      ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
+                  &                                                   - 0.5_wp * zcu   * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
+            END_2D
          END DO
          !
@@ -951,18 +894,16 @@
          !
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
-                  zdx2 = e1u(ji,jj) * e1u(ji,jj)
+            DO_2D_10_10
+               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
+               zdx2 = e1u(ji,jj) * e1u(ji,jj)
 !!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj)
-                  zdx4 = zdx2 * zdx2
-                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (        (                       pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     &
-                     &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) &
-                     &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) * (                       ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
-                     &                                                   - 0.5_wp * zcu   * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) &
-                     &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)     &
-                     &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) )
-               END DO
-            END DO
+               zdx4 = zdx2 * zdx2
+               pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (        (                       pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     &
+                  &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) &
+                  &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) * (                       ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
+                  &                                                   - 0.5_wp * zcu   * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) &
+                  &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)     &
+                  &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) )
+            END_2D
          END DO
          !
@@ -974,21 +915,17 @@
       IF( ll_neg ) THEN
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  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)   &
-                        &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
-                  ENDIF
-               END DO
-            END DO
+            DO_2D_10_10
+               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)   &
+                     &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
+               ENDIF
+            END_2D
          END DO
       ENDIF
       !                                                     !-- High order flux in i-direction  --!
       DO jl = 1, jpl
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl)
-            END DO
-         END DO
+         DO_2D_10_10
+            pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl)
+         END_2D
       END DO
       !
@@ -1021,14 +958,10 @@
       !                                                     !--  Laplacian in j-direction  --!
       DO jl = 1, jpl
-         DO jj = 1, jpjm1         ! First derivative (gradient)
-            DO ji = fs_2, fs_jpim1
-               ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
-            END DO
-         END DO
-         DO jj = 2, jpjm1         ! Second derivative (Laplacian)
-            DO ji = fs_2, fs_jpim1
-               ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj)
-            END DO
-         END DO
+         DO_2D_10_00
+            ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
+         END_2D
+         DO_2D_00_00
+            ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj)
+         END_2D
       END DO
       CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1. )
@@ -1036,14 +969,10 @@
       !                                                     !--  BiLaplacian in j-direction  --!
       DO jl = 1, jpl
-         DO jj = 1, jpjm1         ! First derivative
-            DO ji = fs_2, fs_jpim1
-               ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
-            END DO
-         END DO
-         DO jj = 2, jpjm1         ! Second derivative
-            DO ji = fs_2, fs_jpim1
-               ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj)
-            END DO
-         END DO
+         DO_2D_10_00
+            ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
+         END_2D
+         DO_2D_00_00
+            ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj)
+         END_2D
       END DO
       CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1. )
@@ -1054,69 +983,59 @@
       CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  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) ) )
-               END DO
-            END DO
+            DO_2D_10_10
+               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) ) )
+            END_2D
          END DO
          !
       CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  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)   &
-                     &                                                            - zcv *   ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
-               END DO
-            END DO
+            DO_2D_10_10
+               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)   &
+                  &                                                            - zcv *   ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
+            END_2D
          END DO
          !
       CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
-                  zdy2 = e2v(ji,jj) * e2v(ji,jj)
+            DO_2D_10_10
+               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
+               zdy2 = e2v(ji,jj) * e2v(ji,jj)
 !!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj)
-                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (      (                         pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
-                     &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
-                     &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
-                     &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
-               END DO
-            END DO
+               pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (      (                         pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
+                  &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
+                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
+                  &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
+            END_2D
          END DO
          !
       CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
-                  zdy2 = e2v(ji,jj) * e2v(ji,jj)
+            DO_2D_10_10
+               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
+               zdy2 = e2v(ji,jj) * e2v(ji,jj)
 !!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj)
-                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (      (                         pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
-                     &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
-                     &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
-                     &                                                   - 0.5_wp * zcv   * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
-               END DO
-            END DO
+               pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (      (                         pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
+                  &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
+                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
+                  &                                                   - 0.5_wp * zcv   * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
+            END_2D
          END DO
          !
       CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
-                  zdy2 = e2v(ji,jj) * e2v(ji,jj)
+            DO_2D_10_10
+               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
+               zdy2 = e2v(ji,jj) * e2v(ji,jj)
 !!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj)
-                  zdy4 = zdy2 * zdy2
-                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
-                     &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
-                     &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) * (                       ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
-                     &                                                   - 0.5_wp * zcv   * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) &
-                     &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl)     &
-                     &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) )
-               END DO
-            END DO
+               zdy4 = zdy2 * zdy2
+               pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
+                  &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
+                  &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) * (                       ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
+                  &                                                   - 0.5_wp * zcv   * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) &
+                  &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl)     &
+                  &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) )
+            END_2D
          END DO
          !
@@ -1128,21 +1047,17 @@
       IF( ll_neg ) THEN
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  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) )  &
-                        &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
-                  ENDIF
-               END DO
-            END DO
+            DO_2D_10_10
+               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) )  &
+                     &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
+               ENDIF
+            END_2D
          END DO
       ENDIF
       !                                                     !-- High order flux in j-direction  --!
       DO jl = 1, jpl
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl)
-            END DO
-         END DO
+         DO_2D_10_10
+            pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl)
+         END_2D
       END DO
       !
@@ -1178,10 +1093,8 @@
       ! --------------------------------------------------
       DO jl = 1, jpl
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)
-               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)
-            END DO
-         END DO
+         DO_2D_10_10
+            pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)
+            pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)
+         END_2D
       END DO
 
@@ -1197,34 +1110,30 @@
          
          DO jl = 1, jpl
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1 
-                  zti_ups(ji,jj,jl)= pt_ups(ji+1,jj  ,jl)
-                  ztj_ups(ji,jj,jl)= pt_ups(ji  ,jj+1,jl)
-               END DO
-            END DO
+            DO_2D_00_00
+               zti_ups(ji,jj,jl)= pt_ups(ji+1,jj  ,jl)
+               ztj_ups(ji,jj,jl)= pt_ups(ji  ,jj+1,jl)
+            END_2D
          END DO
          CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1., ztj_ups, 'T', 1. )
 
          DO jl = 1, jpl
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  IF ( pfu_ho(ji,jj,jl) * ( pt_ups(ji+1,jj  ,jl) - pt_ups(ji,jj,jl) ) <= 0._wp .AND.  &
-                     & pfv_ho(ji,jj,jl) * ( pt_ups(ji  ,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0._wp ) THEN
-                     !
-                     IF(  pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj  ,jl) - zti_ups(ji,jj,jl) ) <= 0._wp .AND.  &
-                        & pfv_ho(ji,jj,jl) * ( ztj_ups(ji  ,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0._wp ) THEN
-                        pfu_ho(ji,jj,jl)=0._wp
-                        pfv_ho(ji,jj,jl)=0._wp
-                     ENDIF
-                     !
-                     IF(  pfu_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji-1,jj  ,jl) ) <= 0._wp .AND.  &
-                        & pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji  ,jj-1,jl) ) <= 0._wp ) THEN
-                        pfu_ho(ji,jj,jl)=0._wp
-                        pfv_ho(ji,jj,jl)=0._wp
-                     ENDIF
-                     !
+            DO_2D_00_00
+               IF ( pfu_ho(ji,jj,jl) * ( pt_ups(ji+1,jj  ,jl) - pt_ups(ji,jj,jl) ) <= 0._wp .AND.  &
+                  & pfv_ho(ji,jj,jl) * ( pt_ups(ji  ,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0._wp ) THEN
+                  !
+                  IF(  pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj  ,jl) - zti_ups(ji,jj,jl) ) <= 0._wp .AND.  &
+                     & pfv_ho(ji,jj,jl) * ( ztj_ups(ji  ,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0._wp ) THEN
+                     pfu_ho(ji,jj,jl)=0._wp
+                     pfv_ho(ji,jj,jl)=0._wp
                   ENDIF
-               END DO
-            END DO
+                  !
+                  IF(  pfu_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji-1,jj  ,jl) ) <= 0._wp .AND.  &
+                     & pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji  ,jj-1,jl) ) <= 0._wp ) THEN
+                     pfu_ho(ji,jj,jl)=0._wp
+                     pfv_ho(ji,jj,jl)=0._wp
+                  ENDIF
+                  !
+               ENDIF
+            END_2D
          END DO
          CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
@@ -1238,54 +1147,50 @@
       DO jl = 1, jpl
          
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF    ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
-                  zbup(ji,jj) = -zbig
-                  zbdo(ji,jj) =  zbig
-               ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) > 0._wp ) THEN
-                  zbup(ji,jj) = pt_ups(ji,jj,jl)
-                  zbdo(ji,jj) = pt_ups(ji,jj,jl)
-               ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
-                  zbup(ji,jj) = pt(ji,jj,jl)
-                  zbdo(ji,jj) = pt(ji,jj,jl)
-               ELSE
-                  zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ups(ji,jj,jl) )
-                  zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ups(ji,jj,jl) )
-               ENDIF
-            END DO
-         END DO
-
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               !
-               zup  = MAX( zbup(ji,jj), zbup(ji-1,jj), zbup(ji+1,jj), zbup(ji,jj-1), zbup(ji,jj+1) )  ! search max/min in neighbourhood
-               zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) )
-               !
-               zpos = MAX( 0._wp, pfu_ho(ji-1,jj  ,jl) ) - MIN( 0._wp, pfu_ho(ji  ,jj  ,jl) ) &  ! positive/negative part of the flux
-                  & + MAX( 0._wp, pfv_ho(ji  ,jj-1,jl) ) - MIN( 0._wp, pfv_ho(ji  ,jj  ,jl) )
-               zneg = MAX( 0._wp, pfu_ho(ji  ,jj  ,jl) ) - MIN( 0._wp, pfu_ho(ji-1,jj  ,jl) ) &
-                  & + MAX( 0._wp, pfv_ho(ji  ,jj  ,jl) ) - MIN( 0._wp, pfv_ho(ji  ,jj-1,jl) )
-               !
-               zpos = zpos - (pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) &
-                  &          ) * ( 1. - pamsk )
-               zneg = zneg + (pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) &
-                  &          ) * ( 1. - pamsk )
-               !
-               !                                  ! up & down beta terms
-               ! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!)
-               IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt
-               ELSE                     ; zbetup(ji,jj,jl) = 0._wp ! zbig
-               ENDIF
-               !
-               IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
-               ELSE                     ; zbetdo(ji,jj,jl) = 0._wp ! zbig
-               ENDIF
-               !
-               ! if all the points are outside ice cover
-               IF( zup == -zbig )   zbetup(ji,jj,jl) = 0._wp ! zbig
-               IF( zdo ==  zbig )   zbetdo(ji,jj,jl) = 0._wp ! zbig            
-               !
-            END DO
-         END DO
+         DO_2D_11_11
+            IF    ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
+               zbup(ji,jj) = -zbig
+               zbdo(ji,jj) =  zbig
+            ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) > 0._wp ) THEN
+               zbup(ji,jj) = pt_ups(ji,jj,jl)
+               zbdo(ji,jj) = pt_ups(ji,jj,jl)
+            ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
+               zbup(ji,jj) = pt(ji,jj,jl)
+               zbdo(ji,jj) = pt(ji,jj,jl)
+            ELSE
+               zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ups(ji,jj,jl) )
+               zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ups(ji,jj,jl) )
+            ENDIF
+         END_2D
+
+         DO_2D_00_00
+            !
+            zup  = MAX( zbup(ji,jj), zbup(ji-1,jj), zbup(ji+1,jj), zbup(ji,jj-1), zbup(ji,jj+1) )  ! search max/min in neighbourhood
+            zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) )
+            !
+            zpos = MAX( 0._wp, pfu_ho(ji-1,jj  ,jl) ) - MIN( 0._wp, pfu_ho(ji  ,jj  ,jl) ) &  ! positive/negative part of the flux
+               & + MAX( 0._wp, pfv_ho(ji  ,jj-1,jl) ) - MIN( 0._wp, pfv_ho(ji  ,jj  ,jl) )
+            zneg = MAX( 0._wp, pfu_ho(ji  ,jj  ,jl) ) - MIN( 0._wp, pfu_ho(ji-1,jj  ,jl) ) &
+               & + MAX( 0._wp, pfv_ho(ji  ,jj  ,jl) ) - MIN( 0._wp, pfv_ho(ji  ,jj-1,jl) )
+            !
+            zpos = zpos - (pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) &
+               &          ) * ( 1. - pamsk )
+            zneg = zneg + (pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) &
+               &          ) * ( 1. - pamsk )
+            !
+            !                                  ! up & down beta terms
+            ! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!)
+            IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt
+            ELSE                     ; zbetup(ji,jj,jl) = 0._wp ! zbig
+            ENDIF
+            !
+            IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
+            ELSE                     ; zbetdo(ji,jj,jl) = 0._wp ! zbig
+            ENDIF
+            !
+            ! if all the points are outside ice cover
+            IF( zup == -zbig )   zbetup(ji,jj,jl) = 0._wp ! zbig
+            IF( zdo ==  zbig )   zbetdo(ji,jj,jl) = 0._wp ! zbig            
+            !
+         END_2D
       END DO
       CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
@@ -1295,29 +1200,25 @@
       ! ---------------------------------
       DO jl = 1, jpl
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               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) )
-               zcu = 0.5_wp + SIGN( 0.5_wp , pfu_ho(ji,jj,jl) )
-               !
-               zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu )
-               !
-               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_ups(ji,jj,jl)
-               !
-            END DO
-         END DO
-
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               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) )
-               zcv = 0.5_wp + SIGN( 0.5_wp , pfv_ho(ji,jj,jl) )
-               !
-               zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv )
-               !
-               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_ups(ji,jj,jl)
-               !
-            END DO
-         END DO
+         DO_2D_10_10
+            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) )
+            zcu = 0.5_wp + SIGN( 0.5_wp , pfu_ho(ji,jj,jl) )
+            !
+            zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu )
+            !
+            pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_ups(ji,jj,jl)
+            !
+         END_2D
+
+         DO_2D_10_10
+            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) )
+            zcv = 0.5_wp + SIGN( 0.5_wp , pfv_ho(ji,jj,jl) )
+            !
+            zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv )
+            !
+            pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_ups(ji,jj,jl)
+            !
+         END_2D
 
       END DO
@@ -1344,75 +1245,71 @@
       !
       DO jl = 1, jpl
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1)
-            END DO
-         END DO
+         DO_2D_00_00
+            zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1)
+         END_2D
       END DO
       CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.)   ! lateral boundary cond.
       
       DO jl = 1, jpl
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
-               
-               Rjm = zslpx(ji-1,jj,jl)
-               Rj  = zslpx(ji  ,jj,jl)
-               Rjp = zslpx(ji+1,jj,jl)
-
-               IF( np_limiter == 3 ) THEN
-
-                  IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm
-                  ELSE                        ;   Rr = Rjp
+         DO_2D_00_00
+            uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
+            
+            Rjm = zslpx(ji-1,jj,jl)
+            Rj  = zslpx(ji  ,jj,jl)
+            Rjp = zslpx(ji+1,jj,jl)
+
+            IF( np_limiter == 3 ) THEN
+
+               IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm
+               ELSE                        ;   Rr = Rjp
+               ENDIF
+
+               zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)     
+               IF( Rj > 0. ) THEN
+                  zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)),  &
+                     &        MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
+               ELSE
+                  zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)),  &
+                     &        MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
+               ENDIF
+               pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter
+
+            ELSEIF( np_limiter == 2 ) THEN
+               IF( Rj /= 0. ) THEN
+                  IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
+                  ELSE                        ;   Cr = Rjp / Rj
                   ENDIF
-
-                  zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)     
-                  IF( Rj > 0. ) THEN
-                     zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)),  &
-                        &        MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
-                  ELSE
-                     zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)),  &
-                        &        MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
-                  ENDIF
-                  pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter
-
-               ELSEIF( np_limiter == 2 ) THEN
-                  IF( Rj /= 0. ) THEN
-                     IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
-                     ELSE                        ;   Cr = Rjp / Rj
-                     ENDIF
-                  ELSE
-                     Cr = 0.
-                  ENDIF
-
-                  ! -- superbee --
-                  zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
-                  ! -- van albada 2 --
-                  !!zpsi = 2.*Cr / (Cr*Cr+1.)
-                  ! -- sweby (with beta=1) --
-                  !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
-                  ! -- van Leer --
-                  !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
-                  ! -- ospre --
-                  !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
-                  ! -- koren --
-                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
-                  ! -- charm --
-                  !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
-                  !ELSE                 ;   zpsi = 0.
-                  !ENDIF
-                  ! -- van albada 1 --
-                  !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
-                  ! -- smart --
-                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
-                  ! -- umist --
-                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
-
-                  ! high order flux corrected by the limiter
-                  pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5
-
+               ELSE
+                  Cr = 0.
                ENDIF
-            END DO
-         END DO
+
+               ! -- superbee --
+               zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
+               ! -- van albada 2 --
+               !!zpsi = 2.*Cr / (Cr*Cr+1.)
+               ! -- sweby (with beta=1) --
+               !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
+               ! -- van Leer --
+               !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
+               ! -- ospre --
+               !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
+               ! -- koren --
+               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
+               ! -- charm --
+               !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
+               !ELSE                 ;   zpsi = 0.
+               !ENDIF
+               ! -- van albada 1 --
+               !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
+               ! -- smart --
+               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
+               ! -- umist --
+               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
+
+               ! high order flux corrected by the limiter
+               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5
+
+            ENDIF
+         END_2D
       END DO
       CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.)   ! lateral boundary cond.
@@ -1439,76 +1336,72 @@
       !
       DO jl = 1, jpl
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1)
-            END DO
-         END DO
+         DO_2D_00_00
+            zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1)
+         END_2D
       END DO
       CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.)   ! lateral boundary cond.
 
       DO jl = 1, jpl
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
-
-               Rjm = zslpy(ji,jj-1,jl)
-               Rj  = zslpy(ji,jj  ,jl)
-               Rjp = zslpy(ji,jj+1,jl)
-
-               IF( np_limiter == 3 ) THEN
-
-                  IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm
-                  ELSE                        ;   Rr = Rjp
+         DO_2D_00_00
+            vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
+
+            Rjm = zslpy(ji,jj-1,jl)
+            Rj  = zslpy(ji,jj  ,jl)
+            Rjp = zslpy(ji,jj+1,jl)
+
+            IF( np_limiter == 3 ) THEN
+
+               IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm
+               ELSE                        ;   Rr = Rjp
+               ENDIF
+
+               zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)     
+               IF( Rj > 0. ) THEN
+                  zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)),  &
+                     &        MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
+               ELSE
+                  zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)),  &
+                     &        MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
+               ENDIF
+               pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter
+
+            ELSEIF( np_limiter == 2 ) THEN
+
+               IF( Rj /= 0. ) THEN
+                  IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
+                  ELSE                        ;   Cr = Rjp / Rj
                   ENDIF
-
-                  zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)     
-                  IF( Rj > 0. ) THEN
-                     zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)),  &
-                        &        MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
-                  ELSE
-                     zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)),  &
-                        &        MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
-                  ENDIF
-                  pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter
-
-               ELSEIF( np_limiter == 2 ) THEN
-
-                  IF( Rj /= 0. ) THEN
-                     IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
-                     ELSE                        ;   Cr = Rjp / Rj
-                     ENDIF
-                  ELSE
-                     Cr = 0.
-                  ENDIF
-
-                  ! -- superbee --
-                  zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
-                  ! -- van albada 2 --
-                  !!zpsi = 2.*Cr / (Cr*Cr+1.)
-                  ! -- sweby (with beta=1) --
-                  !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
-                  ! -- van Leer --
-                  !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
-                  ! -- ospre --
-                  !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
-                  ! -- koren --
-                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
-                  ! -- charm --
-                  !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
-                  !ELSE                 ;   zpsi = 0.
-                  !ENDIF
-                  ! -- van albada 1 --
-                  !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
-                  ! -- smart --
-                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
-                  ! -- umist --
-                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
-
-                  ! high order flux corrected by the limiter
-                  pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5
-
+               ELSE
+                  Cr = 0.
                ENDIF
-            END DO
-         END DO
+
+               ! -- superbee --
+               zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
+               ! -- van albada 2 --
+               !!zpsi = 2.*Cr / (Cr*Cr+1.)
+               ! -- sweby (with beta=1) --
+               !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
+               ! -- van Leer --
+               !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
+               ! -- ospre --
+               !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
+               ! -- koren --
+               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
+               ! -- charm --
+               !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
+               !ELSE                 ;   zpsi = 0.
+               !ENDIF
+               ! -- van albada 1 --
+               !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
+               ! -- smart --
+               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
+               ! -- umist --
+               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
+
+               ! high order flux corrected by the limiter
+               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5
+
+            ENDIF
+         END_2D
       END DO
       CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.)   ! lateral boundary cond.
@@ -1544,40 +1437,38 @@
       DO jl = 1, jpl
 
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
+         DO_2D_11_11
+            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
+               !
+               !                               ! -- check h_ip -- !
+               ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
+               IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
+                  zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
+                  IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
+                     pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl)
+                  ENDIF
+               ENDIF
+               !
+               !                               ! -- check h_i -- !
+               ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
+               zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl)
+               IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
+                  pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m)
+               ENDIF
+               !
+               !                               ! -- check h_s -- !
+               ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
+               zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl)
+               IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
+                  zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
                   !
-                  !                               ! -- check h_ip -- !
-                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
-                  IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
-                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
-                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
-                        pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl)
-                     ENDIF
-                  ENDIF
+                  wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt
+                  hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
                   !
-                  !                               ! -- check h_i -- !
-                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
-                  zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl)
-                  IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
-                     pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m)
-                  ENDIF
-                  !
-                  !                               ! -- check h_s -- !
-                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
-                  zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl)
-                  IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
-                     zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
-                     !
-                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt
-                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
-                     !
-                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
-                     pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl)
-                  ENDIF           
-                  !                  
-               ENDIF
-            END DO
-         END DO
+                  pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
+                  pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl)
+               ENDIF           
+               !                  
+            ENDIF
+         END_2D
       END DO 
       !
@@ -1612,23 +1503,21 @@
       ! -- check snow load -- !
       DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
-                  !
-                  zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
-                  !
-                  IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface
-                     ! put snow excess in the ocean
-                     zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
-                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
-                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
-                     ! correct snow volume and heat content
-                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
-                     pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess
-                  ENDIF
-                  !
+         DO_2D_11_11
+            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
+               !
+               zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
+               !
+               IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface
+                  ! put snow excess in the ocean
+                  zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
+                  wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
+                  hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
+                  ! correct snow volume and heat content
+                  pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
+                  pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess
                ENDIF
-            END DO
-         END DO
+               !
+            ENDIF
+         END_2D
       END DO
       !
