Index: NEMO/branches/2020/test_12905_xios_restart/src/OCE/TRA/traadv_fct.F90
===================================================================
--- NEMO/branches/2020/test_12905_xios_restart/src/OCE/TRA/traadv_fct.F90	(revision 12987)
+++ NEMO/branches/2020/test_12905_xios_restart/src/OCE/TRA/traadv_fct.F90	(revision 13751)
@@ -46,4 +46,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -96,4 +97,14 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
       ENDIF
+      !! -- init to 0
+      zwi(:,:,:) = 0._wp
+      zwx(:,:,:) = 0._wp
+      zwy(:,:,:) = 0._wp
+      zwz(:,:,:) = 0._wp
+      ztu(:,:,:) = 0._wp
+      ztv(:,:,:) = 0._wp
+      zltu(:,:,:) = 0._wp
+      zltv(:,:,:) = 0._wp
+      ztw(:,:,:) = 0._wp
       !
       l_trd = .FALSE.            ! set local switches
@@ -128,6 +139,7 @@
       IF( ll_zAimp ) THEN
          ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk))
-         DO_3D_00_00( 1, jpkm1 )
-            zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t(ji,jj,jk,Krhs)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )   &
+            &                               / e3t(ji,jj,jk,Krhs)
             zwinf(ji,jj,jk) =  p2dt * MIN( wi(ji,jj,jk  ) , 0._wp ) / e3t(ji,jj,jk,Krhs)
             zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs)
@@ -139,5 +151,5 @@
          !        !==  upstream advection with initial mass fluxes & intermediate update  ==!
          !                    !* upstream tracer flux in the i and j direction 
-         DO_3D_10_10( 1, jpkm1 )
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
             ! upstream scheme
             zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )
@@ -148,28 +160,32 @@
             zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji  ,jj+1,jk,jn,Kbb) )
          END_3D
-         !                    !* upstream tracer flux in the k direction *!
-         DO_3D_11_11( 2, jpkm1 )
+         !                               !* upstream tracer flux in the k direction *!
+         DO_3D( 1, 1, 1, 1, 2, jpkm1 )      ! Interior value ( multiplied by wmask)
             zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
             zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) )
             zwz(ji,jj,jk) = 0.5 * ( zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb) ) * wmask(ji,jj,jk)
          END_3D
-         IF( ln_linssh ) THEN    ! top ocean value (only in linear free surface as zwz has been w-masked)
-            IF( ln_isfcav ) THEN             ! top of the ice-shelf cavities and at the ocean surface
-               DO_2D_11_11
+         IF( ln_linssh ) THEN               ! top ocean value (only in linear free surface as zwz has been w-masked)
+            IF( ln_isfcav ) THEN                        ! top of the ice-shelf cavities and at the ocean surface
+               DO_2D( 1, 1, 1, 1 )
                   zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
                END_2D
-            ELSE                             ! no cavities: only at the ocean surface
-               zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
+            ELSE                                        ! no cavities: only at the ocean surface
+               DO_2D( 1, 1, 1, 1 )
+                  zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb)
+               END_2D
             ENDIF
          ENDIF
          !               
-         DO_3D_00_00( 1, jpkm1 )
-            !                             ! total intermediate advective trends
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !* trend and after field with monotonic scheme
+            !                               ! total intermediate advective trends
             ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
                &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
                &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
-            !                             ! update and guess with monotonic sheme
-            pt(ji,jj,jk,jn,Krhs) =                     pt(ji,jj,jk,jn,Krhs) +        ztra   / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
-            zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
+            !                               ! update and guess with monotonic sheme
+            pt(ji,jj,jk,jn,Krhs) =                   pt(ji,jj,jk,jn,Krhs) +       ztra   &
+               &                                  / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk)
+            zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) &
+               &                                  / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
          END_3D
          
@@ -178,5 +194,5 @@
             !
             ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
-            DO_3D_00_00( 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! Interior value ( multiplied by wmask)
                zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
                zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
@@ -184,5 +200,5 @@
                zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes
             END_3D
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) &
                   &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
@@ -202,5 +218,5 @@
          !
          CASE(  2  )                   !- 2nd order centered
-            DO_3D_10_10( 1, jpkm1 )
+            DO_3D( 1, 0, 1, 0, 1, jpkm1 )
                zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx(ji,jj,jk)
                zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy(ji,jj,jk)
@@ -211,19 +227,19 @@
             zltv(:,:,jpk) = 0._wp
             DO jk = 1, jpkm1                 ! Laplacian
-               DO_2D_10_10
+               DO_2D( 1, 0, 1, 0 )                 ! 1st derivative (gradient)
                   ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
                   ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
                END_2D
-               DO_2D_00_00
+               DO_2D( 0, 0, 0, 0 )                 ! 2nd derivative * 1/ 6
                   zltu(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6
                   zltv(ji,jj,jk) = (  ztv(ji,jj,jk) + ztv(ji,jj-1,jk)  ) * r1_6
                END_2D
             END DO
-            CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1. , zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
-            !
-            DO_3D_10_10( 1, jpkm1 )
+            CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1.0_wp , zltv, 'T', 1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
+            !
+            DO_3D( 1, 0, 1, 0, 1, jpkm1 )    ! Horizontal advective fluxes
                zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points
                zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
-               !                                                  ! C4 minus upstream advective fluxes 
+               !                                                        ! C4 minus upstream advective fluxes 
                zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)
                zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
@@ -233,11 +249,11 @@
             ztu(:,:,jpk) = 0._wp             ! Bottom value : flux set to zero
             ztv(:,:,jpk) = 0._wp
-            DO_3D_10_10( 1, jpkm1 )
+            DO_3D( 1, 0, 1, 0, 1, jpkm1 )    ! 1st derivative (gradient)
                ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
                ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
             END_3D
-            CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1. , ztv, 'V', -1. )   ! Lateral boundary cond. (unchanged sgn)
-            !
-            DO_3D_00_00( 1, jpkm1 )
+            CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
+            !
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )    ! Horizontal advective fluxes
                zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points (x2)
                zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
@@ -255,5 +271,5 @@
          !
          CASE(  2  )                   !- 2nd order centered
-            DO_3D_00_00( 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                zwz(ji,jj,jk) =  (  pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
                   &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk)
@@ -262,5 +278,5 @@
          CASE(  4  )                   !- 4th order COMPACT
             CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )   ! zwt = COMPACT interpolation of T at w-point
-            DO_3D_00_00( 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
             END_3D
@@ -272,6 +288,6 @@
          !         
          IF ( ll_zAimp ) THEN
-            DO_3D_00_00( 1, jpkm1 )
-               !                             ! total intermediate advective trends
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !* trend and after field with monotonic scheme
+               !                                                ! total intermediate advective trends
                ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
                   &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
@@ -282,5 +298,5 @@
             CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
             !
-            DO_3D_00_00( 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! Interior value ( multiplied by wmask)
                zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
                zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
@@ -289,5 +305,5 @@
          END IF
          !
-         CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1., zwx, 'U', -1. , zwy, 'V', -1.,  zwz, 'W',  1. )
+         CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1.0_wp, zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp,  zwz, 'W',  1.0_wp )
          !
          !        !==  monotonicity algorithm  ==!
@@ -297,5 +313,5 @@
          !        !==  final trend with corrected fluxes  ==!
          !
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
                &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
@@ -308,5 +324,5 @@
             !
             ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp
-            DO_3D_00_00( 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )      ! Interior value ( multiplied by wmask)
                zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
                zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
@@ -314,5 +330,5 @@
                zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic
             END_3D
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) &
                   &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
@@ -374,12 +390,12 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   ikm1         ! local integer
-      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars
-      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo
-      !!----------------------------------------------------------------------
-      !
-      zbig  = 1.e+40_wp
-      zrtrn = 1.e-15_wp
-      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
+      REAL(dp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars
+      REAL(dp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
+      REAL(dp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo
+      !!----------------------------------------------------------------------
+      !
+      zbig  = 1.e+40_dp
+      zrtrn = 1.e-15_dp
+      zbetup(:,:,:) = 0._dp   ;   zbetdo(:,:,:) = 0._dp
 
       ! Search local extrema
@@ -393,5 +409,5 @@
       DO jk = 1, jpkm1
          ikm1 = MAX(jk-1,1)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
 
             ! search maximum in neighbourhood
@@ -423,27 +439,27 @@
          END_2D
       END DO
-      CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1. , zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
+      CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp )   ! lateral boundary cond. (unchanged sign)
 
       ! 3. monotonic flux in the i & j direction (paa & pbb)
       ! ----------------------------------------
-      DO_3D_00_00( 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
          zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
-         zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) )
+         zcu =       ( 0.5  + SIGN( 0.5_wp , paa(ji,jj,jk) ) )
          paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
 
          zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
          zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
-         zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) )
+         zcv =       ( 0.5  + SIGN( 0.5_wp , pbb(ji,jj,jk) ) )
          pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
 
-! monotonic flux in the k direction, i.e. pcc
-! -------------------------------------------
+      ! monotonic flux in the k direction, i.e. pcc
+      ! -------------------------------------------
          za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
          zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
-         zc =       ( 0.5  + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
+         zc =       ( 0.5  + SIGN( 0.5_wp , pcc(ji,jj,jk+1) ) )
          pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
       END_3D
-      CALL lbc_lnk_multi( 'traadv_fct', paa, 'U', -1. , pbb, 'V', -1. )   ! lateral boundary condition (changed sign)
+      CALL lbc_lnk_multi( 'traadv_fct', paa, 'U', -1.0_wp , pbb, 'V', -1.0_wp )   ! lateral boundary condition (changed sign)
       !
    END SUBROUTINE nonosc
@@ -465,5 +481,5 @@
       !!----------------------------------------------------------------------
       
-      DO_3D_11_11( 3, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 3, jpkm1 )       !==  build the three diagonal matrix  ==!
          zwd (ji,jj,jk) = 4._wp
          zwi (ji,jj,jk) = 1._wp
@@ -479,6 +495,6 @@
       END_3D
       !
-      jk = 2                                          ! Switch to second order centered at top
-      DO_2D_11_11
+      jk = 2                                    ! Switch to second order centered at top
+      DO_2D( 1, 1, 1, 1 )
          zwd (ji,jj,jk) = 1._wp
          zwi (ji,jj,jk) = 0._wp
@@ -488,22 +504,22 @@
       !
       !                       !==  tridiagonal solve  ==!
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )           ! first recurrence
          zwt(ji,jj,2) = zwd(ji,jj,2)
       END_2D
-      DO_3D_11_11( 3, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 3, jpkm1 )
          zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
       END_3D
       !
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )           ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
          pt_out(ji,jj,2) = zwrm(ji,jj,2)
       END_2D
-      DO_3D_11_11( 3, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 3, jpkm1 )
          pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
       END_3D
 
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )           ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
          pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
       END_2D
-      DO_3DS_11_11( jpk-2, 2, -1 )
+      DO_3DS( 1, 1, 1, 1, jpk-2, 2, -1 )
          pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
       END_3D
@@ -530,5 +546,5 @@
       !                      !==  build the three diagonal matrix & the RHS  ==!
       !
-      DO_3D_00_00( 3, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 3, jpkm1 )    ! interior (from jk=3 to jpk-1)
          zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp                 !       diagonal
          zwi (ji,jj,jk) =         wmask(ji,jj,jk)                         ! lower diagonal
@@ -549,5 +565,5 @@
       END IF
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )              ! 2nd order centered at top & bottom
          ikt = mikt(ji,jj) + 1            ! w-point below the 1st  wet point
          ikb = MAX(mbkt(ji,jj), 2)        !     -   above the last wet point
@@ -566,22 +582,22 @@
       !                       !==  tridiagonal solver  ==!
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )           !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1
          zwt(ji,jj,2) = zwd(ji,jj,2)
       END_2D
-      DO_3D_00_00( 3, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 3, jpkm1 )
          zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
       END_3D
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
          pt_out(ji,jj,2) = zwrm(ji,jj,2)
       END_2D
-      DO_3D_00_00( 3, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 3, jpkm1 )
          pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
       END_3D
 
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
          pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
       END_2D
-      DO_3DS_00_00( jpk-2, 2, -1 )
+      DO_3DS( 0, 0, 0, 0, jpk-2, 2, -1 )
          pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
       END_3D
@@ -622,22 +638,22 @@
       kstart =  1  + klev
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )                         !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1
          zwt(ji,jj,kstart) = pD(ji,jj,kstart)
       END_2D
-      DO_3D_00_00( kstart+1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, kstart+1, jpkm1 )
          zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1)
       END_3D
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )                        !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
          pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart)
       END_2D
-      DO_3D_00_00( kstart+1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, kstart+1, jpkm1 )
          pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
       END_3D
 
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )                       !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
          pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
       END_2D
-      DO_3DS_00_00( jpk-2, kstart, -1 )
+      DO_3DS( 0, 0, 0, 0, jpk-2, kstart, -1 )
          pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
       END_3D
