Index: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icealb.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/ICE/icealb.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/ICE/icealb.F90	(revision 13470)
@@ -122,5 +122,5 @@
       !
       DO jl = 1, jpl
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             !
             !---------------------------------------------!
Index: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icecor.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/ICE/icecor.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/ICE/icecor.F90	(revision 13470)
@@ -92,5 +92,5 @@
          zzc = rhoi * r1_rdtice
          DO jl = 1, jpl
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )
                zsal = sv_i(ji,jj,jl)
                sv_i(ji,jj,jl) = MIN(  MAX( rn_simin*v_i(ji,jj,jl) , sv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl)  )
@@ -105,5 +105,5 @@
       !                             !-----------------------------------------------------
       IF( kn == 2 ) THEN            !  Ice drift case: Corrections to avoid wrong values !
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             IF ( at_i(ji,jj) == 0._wp ) THEN    ! what to do if there is no ice
                IF ( at_i(ji+1,jj) == 0._wp )   u_ice(ji  ,jj) = 0._wp   ! right side
Index: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icectl.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/ICE/icectl.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/ICE/icectl.F90	(revision 13470)
@@ -366,5 +366,5 @@
       cl_alname(ialert_id) = ' Very high salinity ' ! name of the alert
       DO jl = 1, jpl
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF( v_i(ji,jj,jl) > epsi10  ) THEN
                IF( sv_i(ji,jj,jl) / v_i(ji,jj,jl) > rn_simax ) THEN
@@ -381,5 +381,5 @@
       cl_alname(ialert_id) = ' Very low salinity ' ! name of the alert
       DO jl = 1, jpl
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF( v_i(ji,jj,jl) > epsi10  ) THEN
                IF( sv_i(ji,jj,jl) / v_i(ji,jj,jl) < rn_simin ) THEN
@@ -396,5 +396,5 @@
       cl_alname(ialert_id) = ' Very cold ice ' ! name of the alert
       DO jl = 1, jpl
-         DO_3D_11_11( 1, nlay_i )
+         DO_3D( 1, 1, 1, 1, 1, nlay_i )
             ztmelts    =  -rTmlt * sz_i(ji,jj,jk,jl) + rt0
             IF( t_i(ji,jj,jk,jl) < -50.+rt0  .AND.  v_i(ji,jj,jl) > epsi10 ) THEN
@@ -410,5 +410,5 @@
       cl_alname(ialert_id) = ' Very warm ice ' ! name of the alert
       DO jl = 1, jpl
-         DO_3D_11_11( 1, nlay_i )
+         DO_3D( 1, 1, 1, 1, 1, nlay_i )
             ztmelts    =  -rTmlt * sz_i(ji,jj,jk,jl) + rt0
             IF( t_i(ji,jj,jk,jl) > ztmelts  .AND.  v_i(ji,jj,jl) > epsi10 ) THEN
@@ -424,5 +424,5 @@
       cl_alname(ialert_id) = ' Very thick ice ' ! name of the alert
       jl = jpl 
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          IF( h_i(ji,jj,jl) > 50._wp ) THEN
             WRITE(numout,*) ' ALERTE :   Very thick ice ',h_i(ji,jj,jl)
@@ -436,5 +436,5 @@
       cl_alname(ialert_id) = ' Very thin ice ' ! name of the alert
       jl = 1 
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          IF( h_i(ji,jj,jl) < rn_himin ) THEN
             WRITE(numout,*) ' ALERTE :   Very thin ice ',h_i(ji,jj,jl)
@@ -447,5 +447,5 @@
       ialert_id = ialert_id + 1 ! reference number of this alert
       cl_alname(ialert_id) = ' Very fast ice ' ! name of the alert
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          IF( MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 2. ) THEN
             WRITE(numout,*) ' ALERTE :   Very fast ice ',MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) )
@@ -458,5 +458,5 @@
       ialert_id = ialert_id + 1 ! reference number of this alert
       cl_alname(ialert_id) = ' Ice on continents ' ! name of the alert
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          IF( tmask(ji,jj,1) == 0._wp .AND. ( at_i(ji,jj) > 0._wp .OR. vt_i(ji,jj) > 0._wp ) ) THEN 
             WRITE(numout,*) ' ALERTE :   Ice on continents ',at_i(ji,jj),vt_i(ji,jj)
@@ -469,5 +469,5 @@
       ialert_id = ialert_id + 1 ! reference number of this alert
       cl_alname(ialert_id) = ' Incompatible ice conc and vol ' ! name of the alert
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          IF(  ( vt_i(ji,jj) == 0._wp .AND. at_i(ji,jj) >  0._wp ) .OR. &
             & ( vt_i(ji,jj) >  0._wp .AND. at_i(ji,jj) == 0._wp ) ) THEN 
Index: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn.F90	(revision 13470)
@@ -127,5 +127,5 @@
          ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length
          ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s 
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. )
             zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. )
@@ -156,5 +156,5 @@
 
             ALLOCATE( zdivu_i(jpi,jpj) )
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
Index: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_adv_pra.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_adv_pra.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_adv_pra.F90	(revision 13470)
@@ -110,5 +110,5 @@
       END WHERE
       DO jl = 1, jpl
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             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), &
@@ -143,5 +143,5 @@
       END DO
       DO jl = 1, jpl
-         DO_3D_00_00( 1, nlay_i )
+         DO_3D( 0, 0, 0, 0, 1, nlay_i )
             zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj  ,jk,jl), ze_i(ji  ,jj+1,jk,jl), &
                &                                                   ze_i(ji-1,jj  ,jk,jl), ze_i(ji  ,jj-1,jk,jl), &
@@ -151,5 +151,5 @@
       END DO
       DO jl = 1, jpl
-         DO_3D_00_00( 1, nlay_s )
+         DO_3D( 0, 0, 0, 0, 1, nlay_s )
             zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj  ,jk,jl), ze_s(ji  ,jj+1,jk,jl), &
                &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), &
@@ -307,5 +307,5 @@
          ! derive open water from ice concentration
          zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                        !--- open water
                &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
@@ -363,5 +363,5 @@
          !
          ! Limitation of moments.                                           
-         DO_2D_00_11
+         DO_2D( 0, 0, 1, 1 )
             !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
             psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 )
@@ -383,5 +383,5 @@
 
          !  Calculate fluxes and moments between boxes i<-->i+1              
-         DO_2D_00_11
+         DO_2D( 0, 0, 1, 1 )
             zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
             zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl)
@@ -408,5 +408,5 @@
          END_2D
 
-         DO_2D_00_10
+         DO_2D( 0, 0, 1, 0 )
             zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 
             zalg  (ji,jj) = zalf
@@ -427,5 +427,5 @@
          END_2D
 
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zbt  =       zbet(ji-1,jj)
             zbt1 = 1.0 - zbet(ji-1,jj)
@@ -441,5 +441,5 @@
 
          !   Put the temporary moments into appropriate neighboring boxes.    
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zbt  =       zbet(ji-1,jj)
             zbt1 = 1.0 - zbet(ji-1,jj)
@@ -461,5 +461,5 @@
          END_2D
 
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zbt  =       zbet(ji,jj)
             zbt1 = 1.0 - zbet(ji,jj)
@@ -520,5 +520,5 @@
          !
          ! Limitation of moments.
-         DO_2D_11_00
+         DO_2D( 1, 1, 0, 0 )
             !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
             psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  )
@@ -540,5 +540,5 @@
  
          !  Calculate fluxes and moments between boxes j<-->j+1              
-         DO_2D_11_00
+         DO_2D( 1, 1, 0, 0 )
             zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
             zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl)
@@ -565,5 +565,5 @@
          END_2D
          !
-         DO_2D_10_00
+         DO_2D( 1, 0, 0, 0 )
             zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 
             zalg  (ji,jj) = zalf
@@ -585,5 +585,5 @@
 
          !  Readjust moments remaining in the box. 
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zbt  =         zbet(ji,jj-1)
             zbt1 = ( 1.0 - zbet(ji,jj-1) )
@@ -599,5 +599,5 @@
 
          !   Put the temporary moments into appropriate neighboring boxes.    
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zbt  =       zbet(ji,jj-1)
             zbt1 = 1.0 - zbet(ji,jj-1)
@@ -620,5 +620,5 @@
          END_2D
 
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zbt  =       zbet(ji,jj)
             zbt1 = 1.0 - zbet(ji,jj)
@@ -679,5 +679,5 @@
       !
       DO jl = 1, jpl
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
                !
@@ -726,5 +726,5 @@
       !                                           ! -- check e_i/v_i -- !
       DO jl = 1, jpl
-         DO_3D_11_11( 1, nlay_i )
+         DO_3D( 1, 1, 1, 1, 1, nlay_i )
             IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
                ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean
@@ -740,5 +740,5 @@
       !                                           ! -- check e_s/v_s -- !
       DO jl = 1, jpl
-         DO_3D_11_11( 1, nlay_s )
+         DO_3D( 1, 1, 1, 1, 1, nlay_s )
             IF ( pv_s(ji,jj,jl) > 0._wp ) THEN
                ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean
@@ -783,5 +783,5 @@
       ! -- check snow load -- !
       DO jl = 1, jpl
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
                !
Index: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_adv_umx.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_adv_umx.F90	(revision 13470)
@@ -114,5 +114,5 @@
       END WHERE
       DO jl = 1, jpl
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             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), &
@@ -147,5 +147,5 @@
       END DO
       DO jl = 1, jpl
-         DO_3D_00_00( 1, nlay_i )
+         DO_3D( 0, 0, 0, 0, 1, nlay_i )
             zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj  ,jk,jl), ze_i(ji  ,jj+1,jk,jl), &
                &                                                   ze_i(ji-1,jj  ,jk,jl), ze_i(ji  ,jj-1,jk,jl), &
@@ -155,5 +155,5 @@
       END DO
       DO jl = 1, jpl
-         DO_3D_00_00( 1, nlay_s )
+         DO_3D( 0, 0, 0, 0, 1, nlay_s )
             zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj  ,jk,jl), ze_s(ji  ,jj+1,jk,jl), &
                &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), &
@@ -191,5 +191,5 @@
       !
       ! --- define velocity for advection: u*grad(H) --- !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          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)
@@ -224,5 +224,5 @@
             IF( .NOT. ALLOCATED(jmsk_small) )   ALLOCATE( jmsk_small(jpi,jpj,jpl) ) 
             DO jl = 1, jpl
-               DO_2D_10_10
+               DO_2D( 1, 0, 1, 0 )
                   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
@@ -380,5 +380,5 @@
          !== Open water area ==!
          zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             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
@@ -490,5 +490,5 @@
       IF( pamsk == 0._wp ) THEN
          DO jl = 1, jpl
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 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)
@@ -512,5 +512,5 @@
          ! thus we calculate the upstream solution and apply a limiter again
          DO jl = 1, jpl
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) )
                !
@@ -533,5 +533,5 @@
       IF( PRESENT( pua_ho ) ) THEN
          DO jl = 1, jpl
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                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)
@@ -543,5 +543,5 @@
       ! ---------------------------------
       DO jl = 1, jpl
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) )  
             !
@@ -577,5 +577,5 @@
          !
          DO jl = 1, jpl
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                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)
@@ -588,5 +588,5 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO_2D_10_10
+               DO_2D( 1, 0, 1, 0 )
                   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
@@ -594,5 +594,5 @@
             !
             DO jl = 1, jpl              !-- first guess of tracer from u-flux
-               DO_2D_00_00
+               DO_2D( 0, 0, 0, 0 )
                   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)
@@ -604,5 +604,5 @@
             !
             DO jl = 1, jpl              !-- flux in y-direction
-               DO_2D_10_10
+               DO_2D( 1, 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
@@ -612,5 +612,5 @@
             !
             DO jl = 1, jpl              !-- flux in y-direction
-               DO_2D_10_10
+               DO_2D( 1, 0, 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
@@ -618,5 +618,5 @@
             !
             DO jl = 1, jpl              !-- first guess of tracer from v-flux
-               DO_2D_00_00
+               DO_2D( 0, 0, 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)
@@ -628,5 +628,5 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO_2D_10_10
+               DO_2D( 1, 0, 1, 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
@@ -638,5 +638,5 @@
       !
       DO jl = 1, jpl                    !-- after tracer with upstream scheme
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ztra = - (   pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj  ,jl)   &
                &       + pfv_ups(ji,jj,jl) - pfv_ups(ji  ,jj-1,jl) ) &
@@ -677,5 +677,5 @@
          !
          DO jl = 1, jpl
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                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) )
@@ -695,5 +695,5 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO_2D_10_10
+               DO_2D( 1, 0, 1, 0 )
                   pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
                END_2D
@@ -702,5 +702,5 @@
 
             DO jl = 1, jpl              !-- first guess of tracer from u-flux
-               DO_2D_00_00
+               DO_2D( 0, 0, 0, 0 )
                   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)
@@ -712,5 +712,5 @@
 
             DO jl = 1, jpl              !-- flux in y-direction
-               DO_2D_10_10
+               DO_2D( 1, 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
@@ -721,5 +721,5 @@
             !
             DO jl = 1, jpl              !-- flux in y-direction
-               DO_2D_10_10
+               DO_2D( 1, 0, 1, 0 )
                   pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
                END_2D
@@ -728,5 +728,5 @@
             !
             DO jl = 1, jpl              !-- first guess of tracer from v-flux
-               DO_2D_00_00
+               DO_2D( 0, 0, 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)
@@ -738,5 +738,5 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO_2D_10_10
+               DO_2D( 1, 0, 1, 0 )
                   pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) )
                END_2D
@@ -786,5 +786,5 @@
          !                                                        !--  advective form update in zpt  --!
          DO jl = 1, jpl
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                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) &
@@ -813,5 +813,5 @@
          !                                                        !--  advective form update in zpt  --!
          DO jl = 1, jpl
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                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) &
@@ -895,5 +895,5 @@
          !        
          DO jl = 1, jpl
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 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) ) )
@@ -904,5 +904,5 @@
          !
          DO jl = 1, jpl
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 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)   &
@@ -914,5 +914,5 @@
          !
          DO jl = 1, jpl
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
@@ -928,5 +928,5 @@
          !
          DO jl = 1, jpl
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
@@ -942,5 +942,5 @@
          !
          DO jl = 1, jpl
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
@@ -963,5 +963,5 @@
       IF( ll_neg ) THEN
          DO jl = 1, jpl
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 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)   &
@@ -973,5 +973,5 @@
       !                                                     !-- High order flux in i-direction  --!
       DO jl = 1, jpl
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl)
          END_2D
@@ -1006,8 +1006,8 @@
       !                                                     !--  Laplacian in j-direction  --!
       DO jl = 1, jpl
-         DO_2D_10_00
+         DO_2D( 1, 0, 0, 0 )
             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
+         DO_2D( 0, 0, 0, 0 )
             ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj)
          END_2D
@@ -1017,8 +1017,8 @@
       !                                                     !--  BiLaplacian in j-direction  --!
       DO jl = 1, jpl
-         DO_2D_10_00
+         DO_2D( 1, 0, 0, 0 )
             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
+         DO_2D( 0, 0, 0, 0 )
             ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj)
          END_2D
@@ -1031,5 +1031,5 @@
       CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
          DO jl = 1, jpl
-            DO_2D_10_10
+            DO_2D( 1, 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) ) )
@@ -1039,5 +1039,5 @@
       CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
          DO jl = 1, jpl
-            DO_2D_10_10
+            DO_2D( 1, 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)   &
@@ -1048,5 +1048,5 @@
       CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
          DO jl = 1, jpl
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
@@ -1061,5 +1061,5 @@
       CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
          DO jl = 1, jpl
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
@@ -1074,5 +1074,5 @@
       CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
          DO jl = 1, jpl
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
@@ -1095,5 +1095,5 @@
       IF( ll_neg ) THEN
          DO jl = 1, jpl
-            DO_2D_10_10
+            DO_2D( 1, 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) )  &
@@ -1105,5 +1105,5 @@
       !                                                     !-- High order flux in j-direction  --!
       DO jl = 1, jpl
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl)
          END_2D
@@ -1141,5 +1141,5 @@
       ! --------------------------------------------------
       DO jl = 1, jpl
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             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)
@@ -1158,5 +1158,5 @@
          
          DO jl = 1, jpl
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zti_ups(ji,jj,jl)= pt_ups(ji+1,jj  ,jl)
                ztj_ups(ji,jj,jl)= pt_ups(ji  ,jj+1,jl)
@@ -1166,5 +1166,5 @@
 
          DO jl = 1, jpl
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                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
@@ -1195,5 +1195,5 @@
       DO jl = 1, jpl
          
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF    ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
                zbup(ji,jj) = -zbig
@@ -1211,5 +1211,5 @@
          END_2D
 
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             !
             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
@@ -1248,5 +1248,5 @@
       ! ---------------------------------
       DO jl = 1, jpl
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 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) )
@@ -1259,5 +1259,5 @@
          END_2D
 
-         DO_2D_10_10
+         DO_2D( 1, 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) )
@@ -1293,5 +1293,5 @@
       !
       DO jl = 1, jpl
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1)
          END_2D
@@ -1300,5 +1300,5 @@
       
       DO jl = 1, jpl
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
             
@@ -1384,5 +1384,5 @@
       !
       DO jl = 1, jpl
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1)
          END_2D
@@ -1391,5 +1391,5 @@
 
       DO jl = 1, jpl
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
 
@@ -1488,5 +1488,5 @@
       !
       DO jl = 1, jpl
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
                !
@@ -1535,5 +1535,5 @@
       !                                           ! -- check e_i/v_i -- !
       DO jl = 1, jpl
-         DO_3D_11_11( 1, nlay_i )
+         DO_3D( 1, 1, 1, 1, 1, nlay_i )
             IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
                ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean
@@ -1549,5 +1549,5 @@
       !                                           ! -- check e_s/v_s -- !
       DO jl = 1, jpl
-         DO_3D_11_11( 1, nlay_s )
+         DO_3D( 1, 1, 1, 1, 1, nlay_s )
             IF ( pv_s(ji,jj,jl) > 0._wp ) THEN
                ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean
@@ -1592,5 +1592,5 @@
       ! -- check snow load -- !
       DO jl = 1, jpl
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
                !
Index: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_rdgrft.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_rdgrft.F90	(revision 13470)
@@ -159,5 +159,5 @@
       npti = 0   ;   nptidx(:) = 0
       ipti = 0   ;   iptidx(:) = 0
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          IF ( at_i(ji,jj) > epsi10 ) THEN
             npti           = npti + 1
@@ -775,5 +775,5 @@
       !                              !--------------------------------------------------!
       CASE( 1 )               !--- Spatial smoothing
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN 
                zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
@@ -786,5 +786,5 @@
          END_2D
          
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             strength(ji,jj) = zworka(ji,jj)
          END_2D
@@ -797,5 +797,5 @@
          ENDIF
          !
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN 
                itframe = 1 ! number of time steps for the running mean
Index: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_rhg_evp.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_rhg_evp.F90	(revision 13470)
@@ -182,5 +182,5 @@
       ! for diagnostics and convergence tests
       ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) )
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
          zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
@@ -192,5 +192,5 @@
       !------------------------------------------------------------------------------!
       ! ocean/land mask
-      DO_2D_10_10
+      DO_2D( 1, 0, 1, 0 )
          zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
       END_2D
@@ -198,5 +198,5 @@
 
       ! Lateral boundary conditions on velocity (modify zfmask)
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          IF( zfmask(ji,jj) == 0._wp ) THEN
             zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), &
@@ -266,5 +266,5 @@
       zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
 
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
 
          ! ice fraction at U-V points
@@ -317,5 +317,5 @@
       !
       IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ! ice thickness at U-V points
             zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
@@ -334,5 +334,5 @@
          !
       ELSE                               !-- no landfast
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ztaux_base(ji,jj) = 0._wp
             ztauy_base(ji,jj) = 0._wp
@@ -351,5 +351,5 @@
          ! convergence test
          IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2  ) THEN
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )
                zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step
                zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1)
@@ -358,5 +358,5 @@
 
          ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
 
             ! shear at F points
@@ -368,5 +368,5 @@
          CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. )
 
-         DO_2D_01_01
+         DO_2D( 0, 1, 0, 1 )
 
             ! shear**2 at T points (doc eq. A16)
@@ -417,10 +417,10 @@
          ! Save beta at T-points for further computations
          IF( ln_aEVP ) THEN
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )
                zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
             END_2D
          ENDIF
          
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
 
             ! alpha for aEVP
@@ -442,5 +442,5 @@
 
          ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             !                   !--- U points
             zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
@@ -470,5 +470,5 @@
          IF( MOD(jter,2) == 0 ) THEN ! even iterations
             !
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                !                 !--- tau_io/(v_oce - v_ice)
                zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
@@ -522,5 +522,5 @@
             IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
             !
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                !                 !--- tau_io/(u_oce - u_ice)
                zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
@@ -576,5 +576,5 @@
          ELSE ! odd iterations
             !
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                !                 !--- tau_io/(u_oce - u_ice)
                zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
@@ -628,5 +628,5 @@
             IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
             !
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                !                 !--- tau_io/(v_oce - v_ice)
                zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
@@ -693,5 +693,5 @@
       ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 
       !------------------------------------------------------------------------------!
-      DO_2D_10_10
+      DO_2D( 1, 0, 1, 0 )
 
          ! shear at F points
@@ -702,5 +702,5 @@
       END_2D
       
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          
          ! tension**2 at T points
@@ -766,5 +766,5 @@
          ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
          !         
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
                &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
@@ -817,5 +817,5 @@
             &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
          !
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ! 2D ice mass, snow mass, area transport arrays (X, Y)
             zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
@@ -917,5 +917,5 @@
          zresm = 0._wp
       ELSE
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
                &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj)
Index: NEMO/branches/2020/temporary_r4_trunk/src/ICE/iceistate.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/ICE/iceistate.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/ICE/iceistate.F90	(revision 13470)
@@ -287,5 +287,5 @@
          ! select ice covered grid points
          npti = 0 ; nptidx(:) = 0
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
                npti         = npti  + 1
@@ -342,5 +342,5 @@
          CALL ice_var_salprof ! for sz_i
          DO jl = 1, jpl
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )
                v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
                v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
@@ -350,5 +350,5 @@
          !
          DO jl = 1, jpl
-            DO_3D_11_11( 1, nlay_s )
+            DO_3D( 1, 1, 1, 1, 1, nlay_s )
                t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
                e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
@@ -358,5 +358,5 @@
          !
          DO jl = 1, jpl
-            DO_3D_11_11( 1, nlay_i )
+            DO_3D( 1, 1, 1, 1, 1, nlay_i )
                t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 
                ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
Index: NEMO/branches/2020/temporary_r4_trunk/src/ICE/iceitd.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/ICE/iceitd.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/ICE/iceitd.F90	(revision 13470)
@@ -97,5 +97,5 @@
       !
       npti = 0   ;   nptidx(:) = 0
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          IF ( at_i(ji,jj) > epsi10 ) THEN
             npti = npti + 1
@@ -604,5 +604,5 @@
          !                    !---------------------------------------
          npti = 0   ;   nptidx(:) = 0
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF( a_i(ji,jj,jl) > 0._wp .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN
                npti = npti + 1
@@ -643,5 +643,5 @@
          !                    !-----------------------------------------
          npti = 0 ; nptidx(:) = 0
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF( a_i(ji,jj,jl+1) > 0._wp .AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN
                npti = npti + 1
Index: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icesbc.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/ICE/icesbc.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/ICE/icesbc.F90	(revision 13470)
@@ -77,5 +77,5 @@
       IF( ln_mixcpl) THEN                                                        ! Case of a mixed Bulk/Coupled formulation
                                    CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             utau_ice(ji,jj) = utau_ice(ji,jj) * xcplmask(ji,jj,0) + zutau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
             vtau_ice(ji,jj) = vtau_ice(ji,jj) * xcplmask(ji,jj,0) + zvtau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
Index: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icethd.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/ICE/icethd.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/ICE/icethd.F90	(revision 13470)
@@ -120,5 +120,5 @@
          zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
          zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zfric(ji,jj) = rn_cio * ( 0.5_wp *  &
                &                    (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
@@ -126,5 +126,5 @@
          END_2D
       ELSE      !  if no ice dynamics => transmit directly the atmospheric stress to the ocean
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zfric(ji,jj) = r1_rau0 * SQRT( 0.5_wp *  &
                &                         (  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
@@ -137,5 +137,5 @@
       ! Partial computation of forcing for the thermodynamic sea ice model
       !--------------------------------------------------------------------!
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          rswitch  = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
          !
@@ -209,5 +209,5 @@
          ! select ice covered grid points
          npti = 0 ; nptidx(:) = 0
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF ( a_i(ji,jj,jl) > epsi10 ) THEN     
                npti         = npti  + 1
Index: NEMO/branches/2020/temporary_r4_trunk/src/ICE/iceupdate.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/ICE/iceupdate.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/ICE/iceupdate.F90	(revision 13470)
@@ -114,5 +114,5 @@
       ENDIF
       
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
 
          ! Solar heat flux reaching the ocean = zqsr (W.m-2) 
@@ -333,5 +333,5 @@
       !
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             !                                               ! 2*(U_ice-U_oce) at T-point
             zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)   
@@ -359,5 +359,5 @@
       ENDIF
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ! ice area at u and v-points 
          zat_u  = ( at_i(ji,jj) * tmask(ji,jj,1) + at_i (ji+1,jj    ) * tmask(ji+1,jj  ,1) )  &
Index: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icevar.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/ICE/icevar.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/ICE/icevar.F90	(revision 13470)
@@ -263,5 +263,5 @@
       zlay_i   = REAL( nlay_i , wp )    ! number of layers
       DO jl = 1, jpl
-         DO_3D_11_11( 1, nlay_i )
+         DO_3D( 1, 1, 1, 1, 1, nlay_i )
             IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area 
                !
@@ -368,5 +368,5 @@
          z1_dS = 1._wp / ( zsi1 - zsi0 )
          DO jl = 1, jpl
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )
                zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp )  )
                !                             ! force a constant profile when SSS too low (Baltic Sea)
@@ -377,5 +377,5 @@
          ! Computation of the profile
          DO jl = 1, jpl
-            DO_3D_11_11( 1, nlay_i )
+            DO_3D( 1, 1, 1, 1, 1, nlay_i )
                !                          ! linear profile with 0 surface value
                zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
@@ -507,5 +507,5 @@
          ! Zap ice energy and use ocean heat to melt ice
          !-----------------------------------------------------------------
-         DO_3D_11_11( 1, nlay_i )
+         DO_3D( 1, 1, 1, 1, 1, nlay_i )
             ! update exchanges with ocean
             hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
@@ -514,5 +514,5 @@
          END_3D
          !
-         DO_3D_11_11( 1, nlay_s )
+         DO_3D( 1, 1, 1, 1, 1, nlay_s )
             ! update exchanges with ocean
             hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
@@ -524,5 +524,5 @@
          ! zap ice and snow volume, add water and salt to ocean
          !-----------------------------------------------------------------
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             ! update exchanges with ocean
             sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_rdtice
@@ -597,5 +597,5 @@
          ! zap ice energy and send it to the ocean
          !----------------------------------------
-         DO_3D_11_11( 1, nlay_i )
+         DO_3D( 1, 1, 1, 1, 1, nlay_i )
             IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
                hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
@@ -604,5 +604,5 @@
          END_3D
          !
-         DO_3D_11_11( 1, nlay_s )
+         DO_3D( 1, 1, 1, 1, 1, nlay_s )
             IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
                hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
@@ -614,5 +614,5 @@
          ! zap ice and snow volume, add water and salt to ocean
          !-----------------------------------------------------
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
                wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
Index: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icewri.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/ICE/icewri.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/ICE/icewri.F90	(revision 13470)
@@ -69,5 +69,5 @@
 
       ! tresholds for outputs
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
          zmsk05(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.05_wp ) ) ! 1 if 5% ice , 0 if less
@@ -76,5 +76,5 @@
       END_2D
       DO jl = 1, jpl
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             zmsk00l(ji,jj,jl)  = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) )
             zmsksnl(ji,jj,jl)  = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi06 ) )
@@ -130,5 +130,5 @@
       !
       IF( iom_use('icevel') .OR. iom_use('fasticepres') ) THEN                                                              ! module of ice velocity
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             z2da  = u_ice(ji,jj) + u_ice(ji-1,jj)
             z2db  = v_ice(ji,jj) + v_ice(ji,jj-1)
Index: NEMO/branches/2020/temporary_r4_trunk/src/OCE/DOM/domain.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/OCE/DOM/domain.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/OCE/DOM/domain.F90	(revision 13470)
@@ -140,5 +140,5 @@
                                        ! Read in masks to define closed seas and lakes 
       !
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          ik = mikt(ji,jj)
          risfdep(ji,jj) = gdepw_0(ji,jj,ik)
Index: NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN/dynspg_ts.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN/dynspg_ts.F90	(revision 13470)
@@ -253,5 +253,5 @@
          IF( ln_wd_il ) THEN                       ! W/D : limiter applied to spgspg
             CALL wad_spg( sshn, zcpx, zcpy )          ! Calculating W/D gravity filters, zcpx and zcpy
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   &
                   &                          * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth
@@ -260,5 +260,5 @@
             END_2D
          ELSE                                      ! now suface pressure gradient
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj)
                zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj) 
@@ -268,5 +268,5 @@
       ENDIF
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
           zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
           zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
@@ -281,5 +281,5 @@
       IF( ln_apr_dyn ) THEN
          IF( ln_bt_fw ) THEN                          ! FORWARD integration: use kt+1/2 pressure (NOW+1/2)
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_frc(ji,jj) = zu_frc(ji,jj) + grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
                zv_frc(ji,jj) = zv_frc(ji,jj) + grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
@@ -287,5 +287,5 @@
          ELSE                                         ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW)
             zztmp = grav * r1_2
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)  &
                     &                                   + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
@@ -299,5 +299,5 @@
       !                                   !  ----------------------------------  !
       IF( ln_bt_fw ) THEN                        ! Add wind forcing
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zu_frc(ji,jj) =  zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj)
             zv_frc(ji,jj) =  zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj)
@@ -305,5 +305,5 @@
       ELSE
          zztmp = r1_rau0 * r1_2
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zu_frc(ji,jj) =  zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj)
             zv_frc(ji,jj) =  zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj)
@@ -443,10 +443,10 @@
             !
             !                          ! ocean u- and v-depth at mid-step   (separate DO-loops remove the need of a lbc_lnk)
-            DO_2D_11_10
+            DO_2D( 1, 1, 1, 0 )
                zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj)                        &
                     &                              * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
                     &                                 + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj)
             END_2D
-            DO_2D_10_11
+            DO_2D( 1, 0, 1, 1 )
                zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj)                        &
                     &                              * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
@@ -508,5 +508,5 @@
          !--        ssh    =  ssh   - delta_t' * [ frc + div( flux      ) ]      --!
          !-------------------------------------------------------------------------!
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zhdiv = (   zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1)   ) * r1_e1e2t(ji,jj)
             ssha_e(ji,jj) = (  sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv )  ) * ssmask(ji,jj)
@@ -533,5 +533,5 @@
          ! Sea Surface Height at u-,v-points (vvl case only)
          IF( .NOT.ln_linssh ) THEN                                
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
                   &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
@@ -553,5 +553,5 @@
          !                             ! Surface pressure gradient
          zldg = ( 1._wp - rn_scal_load ) * grav    ! local factor
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
             zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
@@ -571,5 +571,5 @@
          ! Add tidal astronomical forcing if defined
          IF ( ln_tide .AND. ln_tide_pot ) THEN
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
                zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
@@ -580,5 +580,5 @@
 !jth do implicitly instead
          IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
                zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
@@ -598,5 +598,5 @@
          !------------------------------------------------------------------------------------------------------------------------!
          IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
                          &     + rdtbt * (                   zu_spg(ji,jj)   &
@@ -613,5 +613,5 @@
             !
          ELSE                           !* Flux form
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                !                    ! hu_e, hv_e hold depth at jn,  zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
                !                    ! backward interpolated depth used in spg terms at jn+1/2
@@ -637,5 +637,5 @@
 !jth implicit bottom friction:
          IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                   ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
                   va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
@@ -703,5 +703,5 @@
       IF (ln_bt_fw) THEN
          IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )
                zun_save = un_adv(ji,jj)
                zvn_save = vn_adv(ji,jj)
@@ -734,5 +734,5 @@
       ELSE
          ! At this stage, ssha has been corrected: compute new depths at velocity points
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj) &
                &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)      &
@@ -969,5 +969,5 @@
       ! Max courant number for ext. grav. waves
       !
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
          zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
@@ -1093,5 +1093,5 @@
          SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point
          CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                    &
                     &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   ) * 0.25_wp  
@@ -1099,5 +1099,5 @@
             END_2D
          CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) =             (  ht_n  (ji  ,jj+1) + ht_n  (ji+1,jj+1)      &
                     &                      + ht_n  (ji  ,jj  ) + ht_n  (ji+1,jj  )  )   &
@@ -1110,5 +1110,5 @@
          !
          ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
-         DO_2D_01_01
+         DO_2D( 0, 1, 0, 1 )
             ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
             ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
@@ -1119,5 +1119,5 @@
       CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme)
          ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
-         DO_2D_01_01
+         DO_2D( 0, 1, 0, 1 )
             z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) )
             ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht
@@ -1152,5 +1152,5 @@
             !
             !zhf(:,:) = hbatf(:,:)
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          &
                     &            + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      &
@@ -1171,5 +1171,5 @@
          CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
          ! JC: TBC. hf should be greater than 0 
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj)
          END_2D
@@ -1194,5 +1194,5 @@
       SELECT CASE( nvor_scheme )
       CASE( np_ENT )                ! enstrophy conserving scheme (f-point)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             z1_hu = ssumask(ji,jj) / ( hu_n(ji,jj) + 1._wp - ssumask(ji,jj) )
             z1_hv = ssvmask(ji,jj) / ( hv_n(ji,jj) + 1._wp - ssvmask(ji,jj) )
@@ -1207,5 +1207,5 @@
          !         
       CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
             zy2 = ( zhV(ji,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
@@ -1218,5 +1218,5 @@
          !
       CASE( np_ENS )                ! enstrophy conserving scheme (f-point)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zy1 =   r1_8 * ( zhV(ji  ,jj-1) + zhV(ji+1,jj-1) &
               &            + zhV(ji  ,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
@@ -1228,5 +1228,5 @@
          !
       CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)         
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zhV(ji  ,jj  ) &
              &                                         + ftnw(ji+1,jj) * zhV(ji+1,jj  ) &
@@ -1262,5 +1262,5 @@
       !
       IF( ln_wd_dl_rmp ) THEN     
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF    ( pssh(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN 
                !           IF    ( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin2 ) THEN 
@@ -1273,5 +1273,5 @@
          END_2D
       ELSE  
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF ( pssh(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN   ;   ptmsk(ji,jj) = 1._wp
             ELSE                                                 ;   ptmsk(ji,jj) = 0._wp
@@ -1301,5 +1301,5 @@
       !!----------------------------------------------------------------------
       !
-      DO_2D_11_10
+      DO_2D( 1, 1, 1, 0 )
          IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj) 
          ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj)  
@@ -1309,5 +1309,5 @@
       END_2D
       !
-      DO_2D_10_11
+      DO_2D( 1, 0, 1, 1 )
          IF ( phV(ji,jj) > 0._wp ) THEN   ;   pVmsk(ji,jj) = pTmsk(ji,jj  )
          ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1)  
@@ -1331,5 +1331,5 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
       !!----------------------------------------------------------------------
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                &
               &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
@@ -1397,10 +1397,10 @@
       IF( ln_isfcav.OR.ln_drgice_imp ) THEN          ! top+bottom friction (ocean cavities)
          
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
             pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
          END_2D
       ELSE                          ! bottom friction only
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
             pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
@@ -1412,5 +1412,5 @@
       IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities
          
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ikbu = mbku(ji,jj)       
             ikbv = mbkv(ji,jj)    
@@ -1420,5 +1420,5 @@
       ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities
          
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ikbu = mbku(ji,jj)       
             ikbv = mbkv(ji,jj)    
@@ -1430,5 +1430,5 @@
       IF( ln_wd_il ) THEN      ! W/D : use the "clipped" bottom friction   !!gm   explain WHY, please !
          zztmp = -1._wp / rdtbt
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 & 
                  &                              r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp  )
@@ -1438,5 +1438,5 @@
       ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
          
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj)
             pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj)
@@ -1450,5 +1450,5 @@
          IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity
             
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                iktu = miku(ji,jj)
                iktv = mikv(ji,jj)
@@ -1458,5 +1458,5 @@
          ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity
             
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                iktu = miku(ji,jj)
                iktv = mikv(ji,jj)
@@ -1468,5 +1468,5 @@
          !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
          
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj)
             pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj)
Index: NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN/dynzdf.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN/dynzdf.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN/dynzdf.F90	(revision 13470)
@@ -131,5 +131,5 @@
             va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)
          END DO
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
             ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
@@ -140,5 +140,5 @@
          END_2D
          IF( ln_isfcav.OR.ln_drgice_imp ) THEN    ! Ocean cavities (ISF)
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                iku = miku(ji,jj)         ! top ocean level at u- and v-points 
                ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
@@ -158,5 +158,5 @@
          SELECT CASE( nldf_dyn )
          CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
                zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
@@ -171,5 +171,5 @@
             END_3D
          CASE DEFAULT               ! iso-level lateral mixing
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
                zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
@@ -182,5 +182,5 @@
             END_3D
          END SELECT
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zwi(ji,jj,1) = 0._wp
             ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)
@@ -193,5 +193,5 @@
          SELECT CASE( nldf_dyn )
          CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
                zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
@@ -204,5 +204,5 @@
             END_3D
          CASE DEFAULT               ! iso-level lateral mixing
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
                zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
@@ -213,5 +213,5 @@
             END_3D
          END SELECT
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zwi(ji,jj,1) = 0._wp
             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
@@ -227,5 +227,5 @@
       !
       IF ( ln_drgimp ) THEN      ! implicit bottom friction
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
             ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
@@ -233,5 +233,5 @@
          END_2D
          IF ( ln_isfcav.OR.ln_drgice_imp ) THEN   ! top friction (always implicit)
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
                iku = miku(ji,jj)       ! ocean top level at u- and v-points 
@@ -257,21 +257,21 @@
       !-----------------------------------------------------------------------
       !
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
       END_3D
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
          ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
             &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
       END_2D
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
       END_3D
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
       END_2D
-      DO_3D_00_00( jpk-2, 1, -1 )
+      DO_3D( 0, 0, 0, 0, jpk-2, 1, -1 )
          ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
       END_3D
@@ -284,5 +284,5 @@
          SELECT CASE( nldf_dyn )
          CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
                zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
@@ -297,5 +297,5 @@
             END_3D
          CASE DEFAULT               ! iso-level lateral mixing
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
                zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
@@ -308,5 +308,5 @@
             END_3D
          END SELECT
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zwi(ji,jj,1) = 0._wp
             ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)
@@ -319,5 +319,5 @@
          SELECT CASE( nldf_dyn )
          CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
                zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
@@ -330,5 +330,5 @@
             END_3D
          CASE DEFAULT               ! iso-level lateral mixing
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
                zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
@@ -339,5 +339,5 @@
             END_3D
          END SELECT
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zwi(ji,jj,1) = 0._wp
             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
@@ -352,5 +352,5 @@
       !
       IF( ln_drgimp ) THEN
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
             ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
@@ -358,5 +358,5 @@
          END_2D
          IF ( ln_isfcav.OR.ln_drgice_imp ) THEN
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
@@ -381,21 +381,21 @@
       !-----------------------------------------------------------------------
       !
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
       END_3D
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
          va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
             &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
       END_2D
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
       END_3D
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
       END_2D
-      DO_3D_00_00( jpk-2, 1, -1 )
+      DO_3D( 0, 0, 0, 0, jpk-2, 1, -1 )
          va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
       END_3D
Index: NEMO/branches/2020/temporary_r4_trunk/src/OCE/SBC/sbc_oce.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/OCE/SBC/sbc_oce.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/OCE/SBC/sbc_oce.F90	(revision 13470)
@@ -209,5 +209,5 @@
       !!---------------------------------------------------------------------
       zcoef = 0.5 / ( zrhoa * zcdrag ) 
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ztx = utau(ji-1,jj  ) + utau(ji,jj) 
          zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
Index: NEMO/branches/2020/temporary_r4_trunk/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/OCE/SBC/sbcblk.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/OCE/SBC/sbcblk.F90	(revision 13470)
@@ -408,10 +408,10 @@
 #if defined key_cyclone
       CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj)
          sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj)
       END_2D
 #endif
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
          zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
@@ -470,5 +470,5 @@
 !!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef.
 
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed
          taum  (ji,jj) = zztmp * wndm  (ji,jj)
@@ -485,5 +485,5 @@
       !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
       !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
-      DO_2D_10_10
+      DO_2D( 1, 0, 1, 0 )
          utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) &
             &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
@@ -625,5 +625,5 @@
       !!----------------------------------------------------------------------------------
       !
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          !
          ztmp = rt0 / ptak(ji,jj)
@@ -659,5 +659,5 @@
       !!----------------------------------------------------------------------------------
       !
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
          ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT
@@ -723,5 +723,5 @@
       ! ------------------------------------------------------------ !
       ! C-grid ice dynamics :   U & V-points (same as ocean)
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
          zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  )
@@ -749,5 +749,5 @@
       ! ------------------------------------------------------------ !
       zztmp1 = rn_vfac * 0.5_wp
-      DO_2D_01_01
+      DO_2D( 0, 1, 0, 1 )
          zztmp2 = zrhoa(ji,jj) * Cd_atm(ji,jj) * wndm_ice(ji,jj)
          utau_ice(ji,jj) = zztmp2 * ( sf(jp_wndi)%fnow(ji,jj,1) - zztmp1 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) ) )
@@ -755,5 +755,5 @@
       END_2D
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and  rheology 
          zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj  ,1) )
@@ -1007,5 +1007,5 @@
          !   
          DO jl = 1, jpl                
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )
                zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness
                IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
@@ -1022,5 +1022,5 @@
       !
       DO jl = 1, jpl
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             !                    
             zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness
@@ -1173,5 +1173,5 @@
       zqi_sat(:,:) = 0.98_wp * q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg]
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ! Virtual potential temperature [K]
          zthetav_os = zst(ji,jj)    * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean
Index: NEMO/branches/2020/temporary_r4_trunk/src/OCE/SBC/sbcblk_algo_ncar.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/OCE/SBC/sbcblk_algo_ncar.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/OCE/SBC/sbcblk_algo_ncar.F90	(revision 13470)
@@ -214,5 +214,5 @@
       !!----------------------------------------------------------------------------------
       !
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          !
          zw  = pw10(ji,jj)
@@ -279,5 +279,5 @@
       REAL(wp) :: zzeta, zx2, zx, zpsi_unst, zpsi_stab,  zstab   ! local scalars
       !!----------------------------------------------------------------------------------
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
 
          zzeta = pzeta(ji,jj)
@@ -318,5 +318,5 @@
       !!----------------------------------------------------------------------------------
       !
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          !
          zzeta = pzeta(ji,jj)
@@ -376,5 +376,5 @@
       !!-------------------------------------------------------------------
       !
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          !
          zqa = (1._wp + rctv0*pqa(ji,jj))
Index: NEMO/branches/2020/temporary_r4_trunk/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/OCE/SBC/sbccpl.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/OCE/SBC/sbccpl.F90	(revision 13470)
@@ -1193,5 +1193,5 @@
             !                              
             IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
-               DO_2D_00_00
+               DO_2D( 0, 0, 0, 0 )
                   frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1) )
                   frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji  ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) )
@@ -1217,5 +1217,5 @@
          ! => need to be done only when otx1 was changed
          IF( llnewtx ) THEN
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zzx = frcv(jpr_otx1)%z3(ji-1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1)
                zzy = frcv(jpr_oty1)%z3(ji  ,jj-1,1) + frcv(jpr_oty1)%z3(ji,jj,1)
@@ -1242,5 +1242,5 @@
          IF( llnewtau ) THEN 
             zcoef = 1. / ( zrhoa * zcdrag ) 
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )
                frcv(jpr_w10m)%z3(ji,jj,1) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
             END_2D
@@ -1580,5 +1580,5 @@
             p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
          CASE( 'T' )
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and  rheology 
                zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj  ,1) )
@@ -2463,10 +2463,10 @@
             SELECT CASE( TRIM( sn_snd_crt%cldes ) )
             CASE( 'oce only'             )      ! C-grid ==> T
-               DO_2D_00_00
+               DO_2D( 0, 0, 0, 0 )
                   zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Nii) + uu(ji-1,jj  ,1,Nii) )
                   zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Nii) + vv(ji  ,jj-1,1,Nii) ) 
                END_2D
             CASE( 'weighted oce and ice' )      ! Ocean and Ice on C-grid ==> T  
-               DO_2D_00_00
+               DO_2D( 0, 0, 0, 0 )
                   zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Nii) + uu   (ji-1,jj  ,1,Nii) ) * zfr_l(ji,jj)  
                   zoty1(ji,jj) = 0.5 * ( vv   (ji,jj,1,Nii) + vv   (ji  ,jj-1,1,Nii) ) * zfr_l(ji,jj)
@@ -2476,5 +2476,5 @@
                CALL lbc_lnk_multi( 'sbccpl', zitx1, 'T', -1., zity1, 'T', -1. )
             CASE( 'mixed oce-ice'        )      ! Ocean and Ice on C-grid ==> T
-               DO_2D_00_00
+               DO_2D( 0, 0, 0, 0 )
                   zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Nii) + uu   (ji-1,jj  ,1,Nii) ) * zfr_l(ji,jj)   &
                      &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
@@ -2540,10 +2540,10 @@
           SELECT CASE( TRIM( sn_snd_crtw%cldes ) ) 
           CASE( 'oce only'             )      ! C-grid ==> T 
-             DO_2D_00_00
+             DO_2D( 0, 0, 0, 0 )
                 zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Nii) + uu(ji-1,jj  ,1,Nii) ) 
                 zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Nii) + vv(ji , jj-1,1,Nii) )  
              END_2D
           CASE( 'weighted oce and ice' )      ! Ocean and Ice on C-grid ==> T   
-             DO_2D_00_00
+             DO_2D( 0, 0, 0, 0 )
                 zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Nii) + uu   (ji-1,jj  ,1,Nii) ) * zfr_l(ji,jj)   
                 zoty1(ji,jj) = 0.5 * ( vv   (ji,jj,1,Nii) + vv   (ji  ,jj-1,1,Nii) ) * zfr_l(ji,jj) 
@@ -2553,5 +2553,5 @@
              CALL lbc_lnk_multi( 'sbccpl', zitx1, 'T', -1.,  zity1, 'T', -1. ) 
           CASE( 'mixed oce-ice'        )      ! Ocean and Ice on C-grid ==> T  
-             DO_2D_00_00
+             DO_2D( 0, 0, 0, 0 )
                 zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Nii) + uu   (ji-1,jj  ,1,Nii) ) * zfr_l(ji,jj)   & 
                    &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj) 
Index: NEMO/branches/2020/temporary_r4_trunk/src/OCE/ZDF/zdfdrg.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/OCE/ZDF/zdfdrg.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/OCE/ZDF/zdfdrg.F90	(revision 13470)
@@ -115,5 +115,5 @@
       !
       IF( l_log_not_linssh ) THEN     !==  "log layer"  ==!   compute Cd and -Cd*|U|
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             imk = k_mk(ji,jj)          ! ocean bottom level at t-points
             zut = uu(ji,jj,imk,Nii) + uu(ji-1,jj,imk,Nii)     ! 2 x velocity at t-point
@@ -127,5 +127,5 @@
          END_2D
       ELSE                                            !==  standard Cd  ==!
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             imk = k_mk(ji,jj)    ! ocean bottom level at t-points
             zut = uu(ji,jj,imk,Nii) + uu(ji-1,jj,imk,Nii)     ! 2 x velocity at t-point
@@ -173,5 +173,5 @@
       ENDIF
 
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ikbu = mbku(ji,jj)          ! deepest wet ocean u- & v-levels
          ikbv = mbkv(ji,jj)
@@ -186,5 +186,5 @@
       !
       IF( ln_isfcav ) THEN        ! ocean cavities
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels
             ikbv = mikv(ji,jj)
@@ -434,5 +434,5 @@
             l_log_not_linssh = .FALSE.    !- don't update Cd at each time step
             !
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )
                zzz =  0.5_wp * e3t_0(ji,jj,k_mk(ji,jj))
                zcd = (  vkarmn / LOG( zzz / rn_z0 )  )**2
Index: NEMO/branches/2020/temporary_r4_trunk/src/OCE/ZDF/zdfgls.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/OCE/ZDF/zdfgls.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/OCE/ZDF/zdfgls.F90	(revision 13470)
@@ -177,5 +177,5 @@
       
       ! Compute surface, top and bottom friction at T-points
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1)
       END_2D
@@ -184,5 +184,5 @@
       !    
       IF( .NOT.ln_drg_OFF ) THEN    !== top/bottom friction   (explicit before friction)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
             zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0)
@@ -191,5 +191,5 @@
          END_2D
          IF( ln_isfcav ) THEN       !top friction
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zmsku = ( 2._wp - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
                zmskv = ( 2._wp - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0)
@@ -218,5 +218,5 @@
       zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1)  + (1._wp - tmask(:,:,1))*rn_hsro
       !
-      DO_3D_10_10( 2, jpkm1 )
+      DO_3D( 1, 0, 1, 0, 2, jpkm1 )
          eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk)
       END_3D
@@ -227,5 +227,5 @@
 
       IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             zup   = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1)
             zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) )
@@ -248,5 +248,5 @@
       ! Warning : after this step, en : right hand side of the matrix
 
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          !
          buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk)     ! stratif. destruction
@@ -342,5 +342,5 @@
          !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin
          !                      ! Balance between the production and the dissipation terms
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
 !!gm This means that bottom and ocean w-level above have a specified "en" value.   Sure ????
 !!   With thick deep ocean level thickness, this may be quite large, no ???
@@ -360,5 +360,5 @@
          !
          IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                itop   = mikt(ji,jj)       ! k   top w-point
                itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
@@ -378,5 +378,5 @@
       CASE ( 1 )             ! Neumman boundary condition
          !                      
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
             ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -392,5 +392,5 @@
          END_2D
          IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                itop   = mikt(ji,jj)       ! k   top w-point
                itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
@@ -412,11 +412,11 @@
       ! ----------------------------------------------------------
       !
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
       END_3D
-      DO_3D_00_00( 2, jpk )
+      DO_3D( 0, 0, 0, 0, 2, jpk )
          zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1)
       END_3D
-      DO_3D_00_00( jpk-1, 2, -1 )
+      DO_3D( 0, 0, 0, 0, jpk-1, 2, -1 )
          en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
       END_3D
@@ -433,20 +433,20 @@
       !
       CASE( 0 )               ! k-kl  (Mellor-Yamada)
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk)
          END_3D
          !
       CASE( 1 )               ! k-eps
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             psi(ji,jj,jk)  = eps(ji,jj,jk)
          END_3D
          !
       CASE( 2 )               ! k-w
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) )
          END_3D
          !
       CASE( 3 )               ! generic
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 
          END_3D
@@ -461,5 +461,5 @@
       ! Warning : after this step, en : right hand side of the matrix
 
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          !
          ! psi / k
@@ -557,5 +557,5 @@
          !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot
          !                      ! Balance between the production and the dissipation terms
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
             ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -576,5 +576,5 @@
       CASE ( 1 )             ! Neumman boundary condition
          !                      
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
             ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -604,11 +604,11 @@
       ! ----------------
       !
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
       END_3D
-      DO_3D_00_00( 2, jpk )
+      DO_3D( 0, 0, 0, 0, 2, jpk )
          zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1)
       END_3D
-      DO_3D_00_00( jpk-1, 2, -1 )
+      DO_3D( 0, 0, 0, 0, jpk-1, 2, -1 )
          psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
       END_3D
@@ -620,15 +620,15 @@
       !
       CASE( 0 )               ! k-kl  (Mellor-Yamada)
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin)
          END_3D
          !
       CASE( 1 )               ! k-eps
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             eps(ji,jj,jk) = psi(ji,jj,jk)
          END_3D
          !
       CASE( 2 )               ! k-w
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 
          END_3D
@@ -638,5 +638,5 @@
          zex1  =      ( 1.5_wp + rmm/rnn )
          zex2  = -1._wp / rnn
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2
          END_3D
@@ -646,5 +646,5 @@
       ! Limit dissipation rate under stable stratification
       ! --------------------------------------------------
-      DO_3D_00_00( 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          ! limitation
          eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
@@ -662,5 +662,5 @@
       !
       CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             ! zcof =  l²/q²
             zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
@@ -679,5 +679,5 @@
          !
       CASE ( 2, 3 )               ! Canuto stability functions
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             ! zcof =  l²/q²
             zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
@@ -711,5 +711,5 @@
       ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0)
       zstm(:,:,jpk) = 0.  
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj))
       END_2D
@@ -726,5 +726,5 @@
       !     later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk)
       !     for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0)
-      DO_3D_00_00( 1, jpk )
+      DO_3D( 0, 0, 0, 0, 1, jpk )
          zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk)
          zavt  = zsqen * zstt(ji,jj,jk)
Index: NEMO/branches/2020/temporary_r4_trunk/src/OCE/ZDF/zdftke.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/OCE/ZDF/zdftke.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/OCE/ZDF/zdftke.F90	(revision 13470)
@@ -231,5 +231,5 @@
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
 !! clem: this should be the right formulation but it makes the model unstable unless drags are calculated implicitly
 !!       one way around would be to increase zbbirau 
@@ -249,5 +249,5 @@
       IF( .NOT.ln_drg_OFF ) THEN    !== friction used as top/bottom boundary condition on TKE
          !
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
             zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
@@ -258,5 +258,5 @@
          END_2D
          IF( ln_isfcav ) THEN       ! top friction
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
                zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )
@@ -283,18 +283,18 @@
          zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
          imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
-         DO_3D_11_11( jpkm1, 2, -1 )
+         DO_3D( 1, 1, 1, 1, jpkm1, 2, -1 )
             zus  = zcof * taum(ji,jj)
             IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
          END_3D
          !                               ! finite LC depth
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj))
          END_2D
          zcof = 0.016 / SQRT( zrhoa * zcdrag )
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
             zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok
          END_2D
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             IF ( zus3(ji,jj) /= 0._wp ) THEN               
                ! vertical velocity due to LC   
@@ -318,5 +318,5 @@
       !
       IF( nn_pdl == 1 ) THEN      !* Prandtl number = F( Ri )
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             !                             ! local Richardson number
             zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear )
@@ -326,5 +326,5 @@
       ENDIF
       !         
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          zcof   = zfact1 * tmask(ji,jj,jk)
          !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical
@@ -346,20 +346,20 @@
       END_3D
       !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
-      DO_3D_00_00( 3, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 3, jpkm1 )
          zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
       END_3D
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
       END_2D
-      DO_3D_00_00( 3, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 3, jpkm1 )
          zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1)
       END_3D
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
       END_2D
-      DO_3D_00_00( jpk-2, 2, -1 )
+      DO_3D( 0, 0, 0, 0, jpk-2, 2, -1 )
          en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
       END_3D
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
       END_3D
@@ -373,10 +373,10 @@
       
       IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
                &                                 * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
          END_3D
       ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             jk = nmln(ji,jj)
             en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
@@ -384,5 +384,5 @@
          END_2D
       ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
             zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
@@ -459,5 +459,5 @@
          zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
 #if ! defined key_si3 && ! defined key_cice
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1)
          END_2D
@@ -467,10 +467,10 @@
          !
          CASE( 0 )                      ! No scaling under sea-ice
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)
             END_2D
             !
          CASE( 1 )                      ! scaling with constant sea-ice thickness
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zmxlm(ji,jj,1) =  ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
                   &                          fr_i(ji,jj)   * rn_mxlice           ) * tmask(ji,jj,1)
@@ -478,5 +478,5 @@
             !
          CASE( 2 )                      ! scaling with mean sea-ice thickness
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
 #if defined key_si3
                zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
@@ -490,5 +490,5 @@
             !
          CASE( 3 )                      ! scaling with max sea-ice thickness
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zmaxice = MAXVAL( h_i(ji,jj,:) )
                zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
@@ -499,5 +499,5 @@
 #endif
          !
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
          END_2D
@@ -507,5 +507,5 @@
       ENDIF
       !
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          zrn2 = MAX( rn2(ji,jj,jk), rsmall )
          zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
@@ -522,5 +522,5 @@
       ! where wmask = 0 set zmxlm == p_e3w
       CASE ( 0 )           ! bounded by the distance to surface and bottom
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   &
             &            pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) )
@@ -531,5 +531,5 @@
          !
       CASE ( 1 )           ! bounded by the vertical scale factor
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) )
             zmxlm(ji,jj,jk) = zemxl
@@ -538,8 +538,8 @@
          !
       CASE ( 2 )           ! |dk[xml]| bounded by e3t :
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
          END_3D
-         DO_3D_00_00( jpkm1, 2, -1 )
+         DO_3D( 0, 0, 0, 0, jpkm1, 2, -1 )
             zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
             zmxlm(ji,jj,jk) = zemxl
@@ -548,11 +548,11 @@
          !
       CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
          END_3D
-         DO_3D_00_00( jpkm1, 2, -1 )
+         DO_3D( 0, 0, 0, 0, jpkm1, 2, -1 )
             zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
          END_3D
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
             zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
@@ -566,5 +566,5 @@
       !                     !  Vertical eddy viscosity and diffusivity  (avm and avt)
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      DO_3D_00_00( 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          zsqen = SQRT( en(ji,jj,jk) )
          zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
@@ -576,5 +576,5 @@
       !
       IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             p_avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
          END_3D
Index: NEMO/branches/2020/temporary_r4_trunk/src/TOP/PISCES/P4Z/p4zfechem.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/TOP/PISCES/P4Z/p4zfechem.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/TOP/PISCES/P4Z/p4zfechem.F90	(revision 13470)
@@ -89,5 +89,5 @@
       ! Chemistry is supposed to be fast enough to be at equilibrium
       ! ------------------------------------------------------------
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zTL1(ji,jj,jk)  = ztotlig(ji,jj,jk)
          zkeq            = fekeq(ji,jj,jk)
@@ -104,5 +104,5 @@
 
       zdust = 0.         ! if no dust available
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          ! Scavenging rate of iron. This scavenging rate depends on the load of particles of sea water. 
          ! This parameterization assumes a simple second order kinetics (k[Particles][Fe]).
@@ -170,5 +170,5 @@
       IF( ln_ligand ) THEN
          !
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             zlam1a   = ( 0.369  * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
                 &    + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) )
Index: NEMO/branches/2020/temporary_r4_trunk/src/TOP/PISCES/P4Z/p4zsbc.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/TOP/PISCES/P4Z/p4zsbc.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/TOP/PISCES/P4Z/p4zsbc.F90	(revision 13470)
@@ -126,5 +126,5 @@
             CALL fld_read( kt, 1, sf_river )
             IF( ln_p4z ) THEN
-               DO_2D_11_11
+               DO_2D( 1, 1, 1, 1 )
                   zcoef = ryyss * e1e2t(ji,jj) * h_rnf(ji,jj) 
                   rivalk(ji,jj) =   sf_river(jr_dic)%fnow(ji,jj,1)  &
@@ -142,5 +142,5 @@
                END_2D
             ELSE    !  ln_p5z
-               DO_2D_11_11
+               DO_2D( 1, 1, 1, 1 )
                   zcoef = ryyss * e1e2t(ji,jj) * h_rnf(ji,jj) 
                   rivalk(ji,jj) =   sf_river(jr_dic)%fnow(ji,jj,1)                                    &
@@ -407,5 +407,5 @@
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) ' Level corresponding to 50m depth ',  ik50,' ', gdept_1d(ik50+1)
-         DO_3D_00_00( 1, ik50 )
+         DO_3D( 0, 0, 0, 0, 1, ik50 )
             ze3t   = e3t_0(ji,jj,jk)
             zsurfc =  e1u(ji,jj) * ( 1. - umask(ji  ,jj  ,jk) )   &
@@ -423,5 +423,5 @@
          CALL lbc_lnk( 'p4zsbc', zcmask , 'T', 1. )      ! lateral boundary conditions on cmask   (sign unchanged)
          !
-         DO_3D_11_11( 1, jpk )
+         DO_3D( 1, 1, 1, 1, 1, jpk )
             zexpide   = MIN( 8.,( gdept_n(ji,jj,jk) / 500. )**(-1.5) )
             zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
Index: NEMO/branches/2020/temporary_r4_trunk/src/TOP/PISCES/P4Z/p4zsms.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 13470)
@@ -127,5 +127,5 @@
          xnegtr(:,:,:) = 1.e0
          DO jn = jp_pcs0, jp_pcs1
-            DO_3D_11_11( 1, jpk )
+            DO_3D( 1, 1, 1, 1, 1, jpk )
                IF( ( trb(ji,jj,jk,jn) + tra(ji,jj,jk,jn) ) < 0.e0 ) THEN
                   ztra             = ABS( trb(ji,jj,jk,jn) ) / ( ABS( tra(ji,jj,jk,jn) ) + rtrn )
Index: NEMO/branches/2020/temporary_r4_trunk/tests/CANAL/MY_SRC/diawri.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/tests/CANAL/MY_SRC/diawri.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/tests/CANAL/MY_SRC/diawri.F90	(revision 13470)
@@ -150,5 +150,5 @@
       CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature
       IF ( iom_use("sbt") ) THEN
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             ikbot = mbkt(ji,jj)
             z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem)
@@ -160,5 +160,5 @@
       CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity
       IF ( iom_use("sbs") ) THEN
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             ikbot = mbkt(ji,jj)
             z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal)
@@ -170,5 +170,5 @@
          zztmp = rau0 * 0.25
          z2d(:,:) = 0._wp
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * uu(ji  ,jj,mbku(ji  ,jj),Nii)  )**2   &
                &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Nii)  )**2   &
@@ -185,5 +185,5 @@
       CALL iom_put(  "ssu", uu(:,:,1,Nii) )            ! surface i-current
       IF ( iom_use("sbu") ) THEN
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             ikbot = mbku(ji,jj)
             z2d(ji,jj) = uu(ji,jj,ikbot,Nii)
@@ -195,5 +195,5 @@
       CALL iom_put(  "ssv", vv(:,:,1,Nii) )            ! surface j-current
       IF ( iom_use("sbv") ) THEN
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             ikbot = mbkv(ji,jj)
             z2d(ji,jj) = vv(ji,jj,ikbot,Nii)
@@ -222,5 +222,5 @@
       IF ( iom_use("socegrad") .OR. iom_use("socegrad2") ) THEN
          z3d(:,:,jpk) = 0.
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             zztmp  = tsn(ji,jj,jk,jp_sal)
             zztmpx = ( tsn(ji+1,jj,jk,jp_sal) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,jk,jp_sal) ) * r1_e1u(ji-1,jj)
@@ -236,5 +236,5 @@
          
       IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zztmp  = tsn(ji,jj,1,jp_tem)
             zztmpx = ( tsn(ji+1,jj,1,jp_tem) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) * r1_e1u(ji-1,jj)
@@ -252,5 +252,5 @@
       IF( iom_use("heatc") ) THEN
          z2d(:,:)  = 0._wp 
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
          END_3D
@@ -260,5 +260,5 @@
       IF( iom_use("saltc") ) THEN
          z2d(:,:)  = 0._wp 
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
          END_3D
@@ -268,5 +268,5 @@
       IF( iom_use("salt2c") ) THEN
          z2d(:,:)  = 0._wp 
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
          END_3D
@@ -276,5 +276,5 @@
       IF ( iom_use("eken") .OR. iom_use("eken_int") ) THEN
          z3d(:,:,jpk) = 0._wp 
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             zztmpx = 0.5 * ( uu(ji-1,jj  ,jk,Nii) + uu(ji,jj,jk,Nii) )
             zztmpy = 0.5 * ( vv(ji  ,jj-1,jk,Nii) + vv(ji,jj,jk,Nii) )
@@ -285,5 +285,5 @@
 
          z2d(:,:)  = 0._wp 
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * z3d(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk)
          END_3D
@@ -296,5 +296,5 @@
          
          z3d(:,:,jpk) = 0._wp 
-         DO_3D_10_10( 1, jpkm1 )
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
             z3d(ji,jj,jk) = (  e2v(ji+1,jj  ) * vv(ji+1,jj  ,jk,Nii) - e2v(ji,jj) * vv(ji,jj,jk,Nii)    &
                &             - e1u(ji  ,jj+1) * uu(ji  ,jj+1,jk,Nii) + e1u(ji,jj) * uu(ji,jj,jk,Nii)  ) * r1_e1e2f(ji,jj)
@@ -303,10 +303,10 @@
          CALL iom_put( "relvor", z3d )                  ! relative vorticity
 
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             z3d(ji,jj,jk) = ff_f(ji,jj) + z3d(ji,jj,jk) 
          END_3D
          CALL iom_put( "absvor", z3d )                  ! absolute vorticity
 
-         DO_3D_10_10( 1, jpkm1 )
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
             ze3  = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
                &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  )
@@ -335,5 +335,5 @@
       IF( iom_use("u_heattr") ) THEN
          z2d(:,:) = 0._wp 
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
          END_3D
@@ -344,5 +344,5 @@
       IF( iom_use("u_salttr") ) THEN
          z2d(:,:) = 0.e0 
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
          END_3D
@@ -362,5 +362,5 @@
       IF( iom_use("v_heattr") ) THEN
          z2d(:,:) = 0.e0 
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
          END_3D
@@ -371,5 +371,5 @@
       IF( iom_use("v_salttr") ) THEN
          z2d(:,:) = 0._wp 
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
          END_3D
@@ -380,5 +380,5 @@
       IF( iom_use("tosmint") ) THEN
          z2d(:,:) = 0._wp
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem)
          END_3D
@@ -388,5 +388,5 @@
       IF( iom_use("somint") ) THEN
          z2d(:,:)=0._wp
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
          END_3D
Index: NEMO/branches/2020/temporary_r4_trunk/tests/CANAL/MY_SRC/usrdef_istate.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/tests/CANAL/MY_SRC/usrdef_istate.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/tests/CANAL/MY_SRC/usrdef_istate.F90	(revision 13470)
@@ -201,5 +201,5 @@
       CASE(4)    ! geostrophic zonal pulse
    
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN
                zdu = rn_uzonal
@@ -236,5 +236,5 @@
          zP0 = rau0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp)
          !
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             zx = glamt(ji,jj) * 1.e3
             zy = gphit(ji,jj) * 1.e3
@@ -267,5 +267,5 @@
          ! velocities:
          za = 2._wp * zP0 / zlambda**2
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zx = glamu(ji,jj) * 1.e3
             zy = gphiu(ji,jj) * 1.e3
@@ -282,5 +282,5 @@
          END_2D
          !
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zx = glamv(ji,jj) * 1.e3
             zy = gphiv(ji,jj) * 1.e3
