Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traqsr.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traqsr.F90	(revision 12236)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traqsr.F90	(revision 12340)
@@ -68,4 +68,5 @@
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -197,47 +198,37 @@
          !
          zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               ze0(ji,jj,1) = rn_abs * qsr(ji,jj)
-               ze1(ji,jj,1) = zcoef  * qsr(ji,jj)
-               ze2(ji,jj,1) = zcoef  * qsr(ji,jj)
-               ze3(ji,jj,1) = zcoef  * qsr(ji,jj)
-               zea(ji,jj,1) =          qsr(ji,jj)
-            END DO
+         DO_2D_00_00
+            ze0(ji,jj,1) = rn_abs * qsr(ji,jj)
+            ze1(ji,jj,1) = zcoef  * qsr(ji,jj)
+            ze2(ji,jj,1) = zcoef  * qsr(ji,jj)
+            ze3(ji,jj,1) = zcoef  * qsr(ji,jj)
+            zea(ji,jj,1) =          qsr(ji,jj)
+         END_2D
+         !
+         DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B depending of vertical profile of Chl
+            DO_2D_00_00
+               zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) )
+               irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
+               zekb(ji,jj) = rkrgb(1,irgb)
+               zekg(ji,jj) = rkrgb(2,irgb)
+               zekr(ji,jj) = rkrgb(3,irgb)
+            END_2D
+
+            DO_2D_00_00
+               zc0 = ze0(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * xsi0r       )
+               zc1 = ze1(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekb(ji,jj) )
+               zc2 = ze2(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekg(ji,jj) )
+               zc3 = ze3(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekr(ji,jj) )
+               ze0(ji,jj,jk) = zc0
+               ze1(ji,jj,jk) = zc1
+               ze2(ji,jj,jk) = zc2
+               ze3(ji,jj,jk) = zc3
+               zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
+            END_2D
          END DO
          !
-         DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B depending of vertical profile of Chl
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) )
-                  irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
-                  zekb(ji,jj) = rkrgb(1,irgb)
-                  zekg(ji,jj) = rkrgb(2,irgb)
-                  zekr(ji,jj) = rkrgb(3,irgb)
-               END DO
-            END DO
-
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  zc0 = ze0(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * xsi0r       )
-                  zc1 = ze1(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekb(ji,jj) )
-                  zc2 = ze2(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekg(ji,jj) )
-                  zc3 = ze3(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekr(ji,jj) )
-                  ze0(ji,jj,jk) = zc0
-                  ze1(ji,jj,jk) = zc1
-                  ze2(ji,jj,jk) = zc2
-                  ze3(ji,jj,jk) = zc3
-                  zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         !
-         DO jk = 1, nksr                     !* now qsr induced heat content
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  qsr_hc(ji,jj,jk) = r1_rau0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) )
-               END DO
-            END DO
-         END DO
+         DO_3D_00_00( 1, nksr )
+            qsr_hc(ji,jj,jk) = r1_rau0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) )
+         END_3D
          !
          DEALLOCATE( zekb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d ) 
@@ -247,34 +238,24 @@
          zz0 =        rn_abs   * r1_rau0_rcp      ! surface equi-partition in 2-bands
          zz1 = ( 1. - rn_abs ) * r1_rau0_rcp
-         DO jk = 1, nksr                          ! solar heat absorbed at T-point in the top 400m 
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r )
-                  zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r )
-                  qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 
-               END DO
-            END DO
-         END DO
+         DO_3D_00_00( 1, nksr )
+            zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r )
+            zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r )
+            qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 
+         END_3D
          !
       END SELECT
       !
       !                          !-----------------------------!
-      DO jk = 1, nksr            !  update to the temp. trend  !
-         DO jj = 2, jpjm1        !-----------------------------!
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
-                  &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,Kmm)
-            END DO
-         END DO
-      END DO
+      DO_3D_00_00( 1, nksr )
+         pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
+            &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,Kmm)
+      END_3D
       !
       ! sea-ice: store the 1st ocean level attenuation coefficient
-      DO jj = 2, jpjm1 
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) )
-            ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
-            ENDIF
-         END DO
-      END DO
+      DO_2D_00_00
+         IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) )
+         ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
+         ENDIF
+      END_2D
       CALL lbc_lnk( 'traqsr', fraqsr_1lev(:,:), 'T', 1._wp )
       !
