Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p4zbc.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p4zbc.F90	(revision 12258)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p4zbc.F90	(revision 12340)
@@ -48,4 +48,5 @@
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -111,13 +112,11 @@
       IF( ll_river ) THEN
           jl = n_trc_indcbc(jpno3)
-          DO jj = 1, jpj
-             DO ji = 1, jpi
-                DO jk = 1, nk_rnf(ji,jj)
-                   zcoef = rn_rfact / ( e1e2t(ji,jj) * h_rnf(ji,jj) * rn_cbc_time ) * tmask(ji,jj,1)
-                   zrivdin = rf_trcfac(jl) * sf_trccbc(jl)%fnow(ji,jj,1) * zcoef
-                   tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) - rno3 * zrivdin * rfact
-               ENDDO
-             END DO
-          END DO
+          DO_2D_11_11
+             DO jk = 1, nk_rnf(ji,jj)
+                zcoef = rn_rfact / ( e1e2t(ji,jj) * h_rnf(ji,jj) * rn_cbc_time ) * tmask(ji,jj,1)
+                zrivdin = rf_trcfac(jl) * sf_trccbc(jl)%fnow(ji,jj,1) * zcoef
+                tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) - rno3 * zrivdin * rfact
+            ENDDO
+          END_2D
       ENDIF
       
@@ -146,11 +145,9 @@
          ALLOCATE( zironice(jpi,jpj) )
          !
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zdep    = rfact / e3t(ji,jj,1,Kmm)
-               zwflux  = fmmflx(ji,jj) / 1000._wp
-               zironice(ji,jj) =  MAX( -0.99 * tr(ji,jj,1,jpfer,Kbb), -zwflux * icefeinput * zdep )
-            END DO
-         END DO
+         DO_2D_11_11
+            zdep    = rfact / e3t(ji,jj,1,Kmm)
+            zwflux  = fmmflx(ji,jj) / 1000._wp
+            zironice(ji,jj) =  MAX( -0.99 * tr(ji,jj,1,jpfer,Kbb), -zwflux * icefeinput * zdep )
+         END_2D
          !
          tr(:,:,1,jpfer,Krhs) = tr(:,:,1,jpfer,Krhs) + zironice(:,:)
@@ -300,33 +297,25 @@
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) ' Level corresponding to 50m depth ',  ik50,' ', gdept_1d(ik50+1)
-         DO jk = 1, ik50
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  ze3t   = e3t_0(ji,jj,jk)
-                  zsurfc =  e1u(ji,jj) * ( 1. - umask(ji  ,jj  ,jk) )   &
-                          + e1u(ji,jj) * ( 1. - umask(ji-1,jj  ,jk) )   &
-                          + e2v(ji,jj) * ( 1. - vmask(ji  ,jj  ,jk) )   &
-                          + e2v(ji,jj) * ( 1. - vmask(ji  ,jj-1,jk) )
-                  zsurfp = zsurfc * ze3t / e1e2t(ji,jj)
-                  ! estimation of the coastal slope : 5 km off the coast
-                  ze3t2 = ze3t * ze3t
-                  zcslp = SQRT( ( distcoast*distcoast + ze3t2 ) / ze3t2 )
-                  !
-                  zcmask(ji,jj,jk) = zcmask(ji,jj,jk) + zcslp * zsurfp
-               END DO
-            END DO
-         END DO
+         DO_3D_00_00( 1, ik50 )
+            ze3t   = e3t_0(ji,jj,jk)
+            zsurfc =  e1u(ji,jj) * ( 1. - umask(ji  ,jj  ,jk) )   &
+                    + e1u(ji,jj) * ( 1. - umask(ji-1,jj  ,jk) )   &
+                    + e2v(ji,jj) * ( 1. - vmask(ji  ,jj  ,jk) )   &
+                    + e2v(ji,jj) * ( 1. - vmask(ji  ,jj-1,jk) )
+            zsurfp = zsurfc * ze3t / e1e2t(ji,jj)
+            ! estimation of the coastal slope : 5 km off the coast
+            ze3t2 = ze3t * ze3t
+            zcslp = SQRT( ( distcoast*distcoast + ze3t2 ) / ze3t2 )
+            !
+            zcmask(ji,jj,jk) = zcmask(ji,jj,jk) + zcslp * zsurfp
+         END_3D
          !
          CALL lbc_lnk( 'p4zbc', zcmask , 'T', 1. )      ! lateral boundary conditions on cmask   (sign unchanged)
          !
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zexpide   = MIN( 8.,( gdept(ji,jj,jk,Kmm) / 500. )**(-1.5) )
-                  zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
-                  zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 )
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 1, jpk )
+            zexpide   = MIN( 8.,( gdept(ji,jj,jk,Kmm) / 500. )**(-1.5) )
+            zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
+            zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 )
+         END_3D
          ! Coastal supply of iron
          ! -------------------------
