Index: NEMO/trunk/tests/BENCH/MY_SRC/zdfiwm.F90
===================================================================
--- NEMO/trunk/tests/BENCH/MY_SRC/zdfiwm.F90	(revision 12377)
+++ NEMO/trunk/tests/BENCH/MY_SRC/zdfiwm.F90	(revision 12487)
@@ -48,8 +48,10 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hcri_iwm   ! decay scale for low-mode critical slope dissipation (m)
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id: zdfiwm.F90 8093 2017-05-30 08:13:14Z gm $
-   !! Software governed by the CeCILL licence     (./LICENSE)
+   !! $Id: zdfiwm.F90 12377 2020-02-12 14:39:06Z acc $
+   !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
 CONTAINS
@@ -121,5 +123,5 @@
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zztmp        ! scalar workspace
+      REAL(wp) ::   zztmp, ztmp1, ztmp2        ! scalar workspace
       REAL(wp), DIMENSION(jpi,jpj)     ::   zfact       ! Used for vertical structure
       REAL(wp), DIMENSION(jpi,jpj)     ::   zhdep       ! Ocean depth
@@ -147,21 +149,22 @@
       !                       !* Critical slope mixing: distribute energy over the time-varying ocean depth,
       !                                                 using an exponential decay from the seafloor.
-      DO jj = 1, jpj                ! part independent of the level
-         DO ji = 1, jpi
-            zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
-            zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  )
-            IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj)
-         END DO
-      END DO
-!!gm gde3w ==>>>  check for ssh taken into account.... seem OK gde3w_n=gdept(Kmm) - ssh(Kmm)
-      DO jk = 2, jpkm1              ! complete with the level-dependent part
-         zemx_iwm(:,:,jk) = zfact(:,:) * (  EXP( ( gde3w(:,:,jk  ) - zhdep(:,:) ) / hcri_iwm(:,:) )                      &
-            &                             - EXP( ( gde3w(:,:,jk-1) - zhdep(:,:) ) / hcri_iwm(:,:) )  ) * wmask(:,:,jk)   &
-            &                          / ( gde3w(:,:,jk) - gde3w(:,:,jk-1) )
-
-!!gm delta(gde3w) = e3t(Kmm)  !!  Please verify the grid-point position w versus t-point
+      DO_2D_11_11
+         zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
+         zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  )
+         IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj)
+      END_2D
+!!gm gde3w ==>>>  check for ssh taken into account.... seem OK gde3w_n=gdept(:,:,:,Kmm) - ssh(:,:,Kmm)
+      DO_3D_11_11( 2, jpkm1 )
+         IF ( zfact(ji,jj) == 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization
+            zemx_iwm(ji,jj,jk) = 0._wp
+         ELSE
+            zemx_iwm(ji,jj,jk) = zfact(ji,jj) * (  EXP( ( gde3w(ji,jj,jk  ) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) )     &
+                 &                               - EXP( ( gde3w(ji,jj,jk-1) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) ) )   &
+                 &                            / ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) )
+         ENDIF
+      END_3D
+!!gm delta(gde3w) = e3t(:,:,:,Kmm)  !!  Please verify the grid-point position w versus t-point
 !!gm it seems to me that only 1/hcri_iwm  is used ==>  compute it one for all
 
-      END DO
 
       !                        !* Pycnocline-intensified mixing: distribute energy over the time-varying 
@@ -177,9 +180,7 @@
          END DO
          !
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
-            END DO
-         END DO
+         DO_2D_11_11
+            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
+         END_2D
          !
          DO jk = 2, jpkm1              ! complete with the level-dependent part
@@ -194,9 +195,7 @@
          END DO
          !
-         DO jj= 1, jpj
-            DO ji = 1, jpi
-               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
-            END DO
-         END DO
+         DO_2D_11_11
+            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
+         END_2D
          !
          DO jk = 2, jpkm1              ! complete with the level-dependent part
@@ -223,19 +222,18 @@
 !!gm
       !
-      DO jk = 2, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   &
-                  &                                     * wmask(ji,jj,jk) / zfact(ji,jj)
-            END DO
-         END DO
-      END DO
+      DO_3D_11_11( 2, jpkm1 )
+         IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   &
+            &                                     * wmask(ji,jj,jk) / zfact(ji,jj)
+      END_3D
       zwkb(:,:,1) = zhdep(:,:) * wmask(:,:,1)
       !
-      zweight(:,:,:) = 0._wp
-      DO jk = 2, jpkm1
-         zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_iwm(:,:) * wmask(:,:,jk)                    &
-            &   * (  EXP( -zwkb(:,:,jk) / hbot_iwm(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_iwm(:,:) )  )
-      END DO
+      DO_3D_11_11( 2, jpkm1 )
+         IF ( rn2(ji,jj,jk) <= 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization
+            zweight(ji,jj,jk) = 0._wp
+         ELSE
+            zweight(ji,jj,jk) = rn2(ji,jj,jk) * hbot_iwm(ji,jj)    &
+               &   * (  EXP( -zwkb(ji,jj,jk) / hbot_iwm(ji,jj) ) - EXP( -zwkb(ji,jj,jk-1) / hbot_iwm(ji,jj) )  )
+         ENDIF
+      END_3D
       !
       zfact(:,:) = 0._wp
@@ -244,9 +242,7 @@
       END DO
       !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
-         END DO
-      END DO
+      DO_2D_11_11
+         IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
+      END_2D
       !
       DO jk = 2, jpkm1              ! complete with the level-dependent part
@@ -276,15 +272,11 @@
       !
       IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the
-         DO jk = 2, jpkm1              ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( zReb(ji,jj,jk) > 480.00_wp ) THEN
-                     zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
-                  ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN
-                     zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
-                  ENDIF
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 2, jpkm1 )
+            IF( zReb(ji,jj,jk) > 480.00_wp ) THEN
+               zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
+            ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN
+               zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
+            ENDIF
+         END_3D
       ENDIF
       !
@@ -296,12 +288,8 @@
          zztmp = 0._wp
 !!gm used of glosum 3D....
-         DO jk = 2, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zztmp = zztmp + e3w(ji,jj,jk,Kmm) * e1e2t(ji,jj)   &
-                     &          * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 2, jpkm1 )
+            zztmp = zztmp + e3w(ji,jj,jk,Kmm) * e1e2t(ji,jj)   &
+               &          * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
+         END_3D
          CALL mpp_sum( 'zdfiwm', zztmp )
          zztmp = rau0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing 
@@ -321,13 +309,13 @@
       !      
       IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature
-         DO jk = 2, jpkm1              ! Calculate S/T diffusivity ratio as a function of Reb
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp *                                                                  &
-                      &   TANH(    0.92_wp * (   LOG10(  MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 )  ) - 0.60_wp   )    )   &
-                      &                 ) * wmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+         ztmp1 = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10( 1.e-20_wp ) - 0.60_wp ) )
+         DO_3D_11_11( 2, jpkm1 )
+            ztmp2 = zReb(ji,jj,jk) * 5._wp * r1_6
+            IF ( ztmp2 > 1.e-20_wp .AND. wmask(ji,jj,jk) == 1._wp ) THEN
+               zav_ratio(ji,jj,jk) = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10(ztmp2) - 0.60_wp ) )
+            ELSE
+               zav_ratio(ji,jj,jk) = ztmp1 * wmask(ji,jj,jk)
+            ENDIF
+         END_3D
          CALL iom_put( "av_ratio", zav_ratio )
          DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing
@@ -364,5 +352,5 @@
       CALL iom_put( "emix_iwm", zemx_iwm )
       
-      IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', kdim=jpk)
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', kdim=jpk)
       !
    END SUBROUTINE zdf_iwm
@@ -395,5 +383,4 @@
       !!              de Lavergne et al. in prep., 2017
       !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   inum         ! local integer
       INTEGER  ::   ios
