Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/C14/trcatm_c14.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/C14/trcatm_c14.F90	(revision 11960)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/C14/trcatm_c14.F90	(revision 12340)
@@ -21,4 +21,6 @@
    PUBLIC   trc_atm_c14_ini     ! called in trcini_c14.F90
    !
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -118,29 +120,27 @@
             IF( ierr3 /= 0 )   CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate fareaz' )
       !
-            DO jj = 1 , jpj                       ! from C14b package
-              DO ji = 1 , jpi
-                 IF( gphit(ji,jj) >= yn40 ) THEN
-                    fareaz(ji,jj,1) = 0.
-                    fareaz(ji,jj,2) = 0.
-                    fareaz(ji,jj,3) = 1.
-                 ELSE IF( gphit(ji,jj ) <= ys40) THEN
-                    fareaz(ji,jj,1) = 1.
-                    fareaz(ji,jj,2) = 0.
-                    fareaz(ji,jj,3) = 0.
-                 ELSE IF( gphit(ji,jj) >= yn20 ) THEN
-                    fareaz(ji,jj,1) = 0.
-                    fareaz(ji,jj,2) = 2. * ( 1. - gphit(ji,jj) / yn40 )
-                    fareaz(ji,jj,3) = 2. * gphit(ji,jj) / yn40 - 1.
-                 ELSE IF( gphit(ji,jj) <= ys20 ) THEN
-                    fareaz(ji,jj,1) = 2. * gphit(ji,jj) / ys40 - 1.
-                    fareaz(ji,jj,2) = 2. * ( 1. - gphit(ji,jj) / ys40 )
-                    fareaz(ji,jj,3) = 0.
-                 ELSE
-                    fareaz(ji,jj,1) = 0.
-                    fareaz(ji,jj,2) = 1.
-                    fareaz(ji,jj,3) = 0.
-                 ENDIF
-              END DO
-           END DO
+            DO_2D_11_11
+              IF( gphit(ji,jj) >= yn40 ) THEN
+                 fareaz(ji,jj,1) = 0.
+                 fareaz(ji,jj,2) = 0.
+                 fareaz(ji,jj,3) = 1.
+              ELSE IF( gphit(ji,jj ) <= ys40) THEN
+                 fareaz(ji,jj,1) = 1.
+                 fareaz(ji,jj,2) = 0.
+                 fareaz(ji,jj,3) = 0.
+              ELSE IF( gphit(ji,jj) >= yn20 ) THEN
+                 fareaz(ji,jj,1) = 0.
+                 fareaz(ji,jj,2) = 2. * ( 1. - gphit(ji,jj) / yn40 )
+                 fareaz(ji,jj,3) = 2. * gphit(ji,jj) / yn40 - 1.
+              ELSE IF( gphit(ji,jj) <= ys20 ) THEN
+                 fareaz(ji,jj,1) = 2. * gphit(ji,jj) / ys40 - 1.
+                 fareaz(ji,jj,2) = 2. * ( 1. - gphit(ji,jj) / ys40 )
+                 fareaz(ji,jj,3) = 0.
+              ELSE
+                 fareaz(ji,jj,1) = 0.
+                 fareaz(ji,jj,2) = 1.
+                 fareaz(ji,jj,3) = 0.
+              ENDIF
+            END_2D
       !
          ENDIF
Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/C14/trcsms_c14.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/C14/trcsms_c14.F90	(revision 11960)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/C14/trcsms_c14.F90	(revision 12340)
@@ -26,4 +26,6 @@
    PUBLIC   trc_sms_c14       ! called in trcsms.F90
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -78,40 +80,38 @@
       ! -------------------------------------------------------------------
 
-      DO jj = 1, jpj
-         DO ji = 1, jpi  
-            IF( tmask(ji,jj,1) >  0. ) THEN
-               !
-               zt   = MIN( 40. , ts(ji,jj,1,jp_tem,Kmm) )
-               !
-               !  Computation of solubility zsol in [mol/(L * atm)]
-               !   after Wanninkhof (2014) referencing Weiss (1974)
-               ztp  = ( zt + 273.16 ) * 0.01
-               zsk  = 0.027766 + ztp * ( -0.025888 + 0.0050578 * ztp )   ! [mol/(L * atm)]
-               zsol = EXP( -58.0931 + 90.5069 / ztp  + 22.2940 * LOG( ztp ) + zsk * ts(ji,jj,1,jp_sal,Kmm) )
-               ! convert solubilities [mol/(L * atm)] -> [mol/(m^3 * ppm)]
-               zsol = zsol * 1.e-03
+      DO_2D_11_11
+         IF( tmask(ji,jj,1) >  0. ) THEN
+            !
+            zt   = MIN( 40. , ts(ji,jj,1,jp_tem,Kmm) )
+            !
+            !  Computation of solubility zsol in [mol/(L * atm)]
+            !   after Wanninkhof (2014) referencing Weiss (1974)
+            ztp  = ( zt + 273.16 ) * 0.01
+            zsk  = 0.027766 + ztp * ( -0.025888 + 0.0050578 * ztp )   ! [mol/(L * atm)]
+            zsol = EXP( -58.0931 + 90.5069 / ztp  + 22.2940 * LOG( ztp ) + zsk * ts(ji,jj,1,jp_sal,Kmm) )
+            ! convert solubilities [mol/(L * atm)] -> [mol/(m^3 * ppm)]
+            zsol = zsol * 1.e-03
 
-               ! Computes the Schmidt number of CO2 in seawater
-               !               Wanninkhof-2014
-               zsch = 2116.8 + zt * ( -136.25 + zt * (4.7353 + zt * (-0.092307 + 0.0007555 * zt ) ) )
+            ! Computes the Schmidt number of CO2 in seawater
+            !               Wanninkhof-2014
+            zsch = 2116.8 + zt * ( -136.25 + zt * (4.7353 + zt * (-0.092307 + 0.0007555 * zt ) ) )
 
-               ! Wanninkhof Piston velocity: zpv in units [m/s]
-               zv2 = xkwind * (wndm(ji,jj) * wndm(ji,jj))              ! wind speed module at T points
-               ! chemical enhancement (Wanninkhof & Knox, 1996)
-               IF( ln_chemh ) zv2 = zv2 + 2.5 * ( 0.5246 + zt * (0.016256 + 0.00049946  * zt ) )
-               zv2 = zv2/360000._wp                                    ! conversion cm/h -> m/s
-               !
-               zpv  = ( zv2 * SQRT( 660./ zsch ) ) * ( 1. - fr_i(ji,jj) ) * tmask(ji,jj,1)
+            ! Wanninkhof Piston velocity: zpv in units [m/s]
+            zv2 = xkwind * (wndm(ji,jj) * wndm(ji,jj))              ! wind speed module at T points
+            ! chemical enhancement (Wanninkhof & Knox, 1996)
+            IF( ln_chemh ) zv2 = zv2 + 2.5 * ( 0.5246 + zt * (0.016256 + 0.00049946  * zt ) )
+            zv2 = zv2/360000._wp                                    ! conversion cm/h -> m/s
+            !
+            zpv  = ( zv2 * SQRT( 660./ zsch ) ) * ( 1. - fr_i(ji,jj) ) * tmask(ji,jj,1)
 
-               ! CO2 piston velocity (m/s)
-               exch_co2(ji,jj)= zpv
-               ! CO2 invasion rate (mol/ppm/m2/s) = 1st part of 14C/C exchange velocity
-               exch_c14(ji,jj)= zpv * zsol
-            ELSE
-               exch_co2(ji,jj) = 0._wp
-               exch_c14(ji,jj) = 0._wp
-            ENDIF
-         END DO
-      END DO
+            ! CO2 piston velocity (m/s)
+            exch_co2(ji,jj)= zpv
+            ! CO2 invasion rate (mol/ppm/m2/s) = 1st part of 14C/C exchange velocity
+            exch_c14(ji,jj)= zpv * zsol
+         ELSE
+            exch_co2(ji,jj) = 0._wp
+            exch_c14(ji,jj) = 0._wp
+         ENDIF
+      END_2D
 
       ! Exchange velocity for 14C/C ratio (m/s)
@@ -127,20 +127,14 @@
       !
       ! Add the surface flux to the trend of jp_c14
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            tr(ji,jj,1,jp_c14,Krhs) = tr(ji,jj,1,jp_c14,Krhs) + qtr_c14(ji,jj) / e3t(ji,jj,1,Kmm) 
-         END DO
-      END DO
+      DO_2D_11_11
+         tr(ji,jj,1,jp_c14,Krhs) = tr(ji,jj,1,jp_c14,Krhs) + qtr_c14(ji,jj) / e3t(ji,jj,1,Kmm) 
+      END_2D
       !
       ! Computation of decay effects on jp_c14
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !
-               tr(ji,jj,jk,jp_c14,Krhs) = tr(ji,jj,jk,jp_c14,Krhs) - rlam14 * tr(ji,jj,jk,jp_c14,Kbb) * tmask(ji,jj,jk) 
-               !
-            END DO
-         END DO
-      END DO
+      DO_3D_11_11( 1, jpk )
+         !
+         tr(ji,jj,jk,jp_c14,Krhs) = tr(ji,jj,jk,jp_c14,Krhs) - rlam14 * tr(ji,jj,jk,jp_c14,Kbb) * tmask(ji,jj,jk) 
+         !
+      END_3D
       !
       IF( lrst_trc ) THEN
Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/C14/trcwri_c14.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/C14/trcwri_c14.F90	(revision 11960)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/C14/trcwri_c14.F90	(revision 12340)
@@ -23,4 +23,6 @@
    !   Standard ratio: 1.176E-12 ; Avogadro's nbr = 6.022E+23 at/mol ; bomb C14 traditionally reported as 1.E+26 atoms
    REAL(wp), PARAMETER  :: atomc14 = 1.176 * 6.022E-15   ! conversion factor 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
 
 
@@ -58,14 +60,10 @@
          zz3d(:,:,:) = 0._wp
          !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( tmask(ji,jj,jk) > 0._wp) THEN
-                     z3d (ji,jj,jk) = tr(ji,jj,jk,jp_c14,Kmm)
-                     zz3d(ji,jj,jk) = LOG( z3d(ji,jj,jk) )
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDDO
+         DO_3D_11_11( 1, jpkm1 )
+            IF( tmask(ji,jj,jk) > 0._wp) THEN
+               z3d (ji,jj,jk) = tr(ji,jj,jk,jp_c14,Kmm)
+               zz3d(ji,jj,jk) = LOG( z3d(ji,jj,jk) )
+            ENDIF
+         END_3D
          zres(:,:) = z3d(:,:,1)
 
@@ -73,10 +71,8 @@
          z2d(:,:) =0._wp
          jk = 1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ztemp = zres(ji,jj) / c14sbc(ji,jj)
-               IF( ztemp > 0._wp .AND. tmask(ji,jj,jk) > 0._wp ) z2d(ji,jj) = LOG( ztemp )
-            ENDDO
-         ENDDO
+         DO_2D_11_11
+            ztemp = zres(ji,jj) / c14sbc(ji,jj)
+            IF( ztemp > 0._wp .AND. tmask(ji,jj,jk) > 0._wp ) z2d(ji,jj) = LOG( ztemp )
+         END_2D
          !
          z3d(:,:,:) = 1.d03 * ( z3d(:,:,:) - 1._wp )
@@ -131,4 +127,6 @@
 #endif
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
