Index: NEMO/trunk/src/OCE/SBC/sbcwave.F90
===================================================================
--- NEMO/trunk/src/OCE/SBC/sbcwave.F90	(revision 11536)
+++ NEMO/trunk/src/OCE/SBC/sbcwave.F90	(revision 12377)
@@ -72,5 +72,5 @@
 
    !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -80,5 +80,5 @@
 CONTAINS
 
-   SUBROUTINE sbc_stokes( )
+   SUBROUTINE sbc_stokes( Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE sbc_stokes  ***
@@ -92,4 +92,5 @@
       !! ** action  
       !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kmm ! ocean time level index
       INTEGER  ::   jj, ji, jk   ! dummy loop argument
       INTEGER  ::   ik           ! local integer 
@@ -111,105 +112,87 @@
       IF( ll_st_bv_li ) THEN   ! (Eq. (19) in Breivik et al. (2014) )
          zfac = 2.0_wp * rpi / 16.0_wp
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ! Stokes drift velocity estimated from Hs and Tmean
-               ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp )
-               ! Stokes surface speed
-               tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj))
-               ! Wavenumber scale
-               zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp )
-            END DO
-         END DO
-         DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points
-            DO ji = 1, jpim1
-               zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) )
-               zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) )
-               !
-               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
-               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
-            END DO
-         END DO
+         DO_2D_11_11
+            ! Stokes drift velocity estimated from Hs and Tmean
+            ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp )
+            ! Stokes surface speed
+            tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj))
+            ! Wavenumber scale
+            zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp )
+         END_2D
+         DO_2D_10_10
+            zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) )
+            zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) )
+            !
+            zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
+            zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
+         END_2D
       ELSE IF( ll_st_peakfr ) THEN    ! peak wave number calculated from the peak frequency received by the wave model
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zk_t(ji,jj) = ( 2.0_wp * rpi * wfreq(ji,jj) ) * ( 2.0_wp * rpi * wfreq(ji,jj) ) / grav
-            END DO
-         END DO
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) )
-               zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) )
-               !
-               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
-               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
-            END DO
-         END DO
+         DO_2D_11_11
+            zk_t(ji,jj) = ( 2.0_wp * rpi * wfreq(ji,jj) ) * ( 2.0_wp * rpi * wfreq(ji,jj) ) / grav
+         END_2D
+         DO_2D_10_10
+            zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) )
+            zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) )
+            !
+            zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
+            zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
+         END_2D
       ENDIF
       !
       !                       !==  horizontal Stokes Drift 3D velocity  ==!
       IF( ll_st_bv2014 ) THEN
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) )
-                  zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) )
-                  !                          
-                  zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth
-                  zkh_v = zk_v(ji,jj) * zdep_v
-                  !                                ! Depth attenuation
-                  zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u )
-                  zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v )
-                  !
-                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
-                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+         DO_3D_00_00( 1, jpkm1 )
+            zdep_u = 0.5_wp * ( gdept(ji,jj,jk,Kmm) + gdept(ji+1,jj,jk,Kmm) )
+            zdep_v = 0.5_wp * ( gdept(ji,jj,jk,Kmm) + gdept(ji,jj+1,jk,Kmm) )
+            !                          
+            zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth
+            zkh_v = zk_v(ji,jj) * zdep_v
+            !                                ! Depth attenuation
+            zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u )
+            zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v )
+            !
+            usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
+            vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
+         END_3D
       ELSE IF( ll_st_li2017 .OR. ll_st_peakfr ) THEN
          ALLOCATE( zstokes_psi_u_top(jpi,jpj), zstokes_psi_v_top(jpi,jpj) )
-         DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points
-            DO ji = 1, jpim1
-               zstokes_psi_u_top(ji,jj) = 0._wp
-               zstokes_psi_v_top(ji,jj) = 0._wp
-            END DO
-         END DO
+         DO_2D_10_10
+            zstokes_psi_u_top(ji,jj) = 0._wp
+            zstokes_psi_v_top(ji,jj) = 0._wp
+         END_2D
          zsqrtpi = SQRT(rpi)
          z_two_thirds = 2.0_wp / 3.0_wp
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  zbot_u = ( gdepw_n(ji,jj,jk+1) + gdepw_n(ji+1,jj,jk+1) )  ! 2 * bottom depth
-                  zbot_v = ( gdepw_n(ji,jj,jk+1) + gdepw_n(ji,jj+1,jk+1) )  ! 2 * bottom depth
-                  zkb_u  = zk_u(ji,jj) * zbot_u                             ! 2 * k * bottom depth
-                  zkb_v  = zk_v(ji,jj) * zbot_v                             ! 2 * k * bottom depth
-                  !
-                  zke3_u = MAX(1.e-8_wp, 2.0_wp * zk_u(ji,jj) * e3u_n(ji,jj,jk))     ! 2k * thickness
-                  zke3_v = MAX(1.e-8_wp, 2.0_wp * zk_v(ji,jj) * e3v_n(ji,jj,jk))     ! 2k * thickness
-
-                  ! Depth attenuation .... do u component first..
-                  zdepth      = zkb_u
-                  zsqrt_depth = SQRT(zdepth)
-                  zexp_depth  = EXP(-zdepth)
-                  zstokes_psi_u_bot = 1.0_wp - zexp_depth  &
-                       &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) &
-                       &              + 1.0_wp - (1.0_wp + zdepth)*zexp_depth )
-                  zda_u                    = ( zstokes_psi_u_bot - zstokes_psi_u_top(ji,jj) ) / zke3_u
-                  zstokes_psi_u_top(ji,jj) =   zstokes_psi_u_bot
-
-                  !         ... and then v component
-                  zdepth      =zkb_v
-                  zsqrt_depth = SQRT(zdepth)
-                  zexp_depth  = EXP(-zdepth)
-                  zstokes_psi_v_bot = 1.0_wp - zexp_depth  &
-                       &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) &
-                       &              + 1.0_wp - (1.0_wp + zdepth)*zexp_depth )
-                  zda_v                    = ( zstokes_psi_v_bot - zstokes_psi_v_top(ji,jj) ) / zke3_v
-                  zstokes_psi_v_top(ji,jj) =   zstokes_psi_v_bot
-                  !
-                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
-                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+         DO_3D_00_00( 1, jpkm1 )
+            zbot_u = ( gdepw(ji,jj,jk+1,Kmm) + gdepw(ji+1,jj,jk+1,Kmm) )  ! 2 * bottom depth
+            zbot_v = ( gdepw(ji,jj,jk+1,Kmm) + gdepw(ji,jj+1,jk+1,Kmm) )  ! 2 * bottom depth
+            zkb_u  = zk_u(ji,jj) * zbot_u                             ! 2 * k * bottom depth
+            zkb_v  = zk_v(ji,jj) * zbot_v                             ! 2 * k * bottom depth
+            !
+            zke3_u = MAX(1.e-8_wp, 2.0_wp * zk_u(ji,jj) * e3u(ji,jj,jk,Kmm))     ! 2k * thickness
+            zke3_v = MAX(1.e-8_wp, 2.0_wp * zk_v(ji,jj) * e3v(ji,jj,jk,Kmm))     ! 2k * thickness
+
+            ! Depth attenuation .... do u component first..
+            zdepth      = zkb_u
+            zsqrt_depth = SQRT(zdepth)
+            zexp_depth  = EXP(-zdepth)
+            zstokes_psi_u_bot = 1.0_wp - zexp_depth  &
+                 &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) &
+                 &              + 1.0_wp - (1.0_wp + zdepth)*zexp_depth )
+            zda_u                    = ( zstokes_psi_u_bot - zstokes_psi_u_top(ji,jj) ) / zke3_u
+            zstokes_psi_u_top(ji,jj) =   zstokes_psi_u_bot
+
+            !         ... and then v component
+            zdepth      =zkb_v
+            zsqrt_depth = SQRT(zdepth)
+            zexp_depth  = EXP(-zdepth)
+            zstokes_psi_v_bot = 1.0_wp - zexp_depth  &
+                 &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) &
+                 &              + 1.0_wp - (1.0_wp + zdepth)*zexp_depth )
+            zda_v                    = ( zstokes_psi_v_bot - zstokes_psi_v_top(ji,jj) ) / zke3_v
+            zstokes_psi_v_top(ji,jj) =   zstokes_psi_v_bot
+            !
+            usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
+            vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
+         END_3D
          DEALLOCATE( zstokes_psi_u_top, zstokes_psi_v_top )
       ENDIF
@@ -220,14 +203,10 @@
       !                       !==  vertical Stokes Drift 3D velocity  ==!
       !
-      DO jk = 1, jpkm1               ! Horizontal e3*divergence
-         DO jj = 2, jpj
-            DO ji = fs_2, jpi
-               ze3divh(ji,jj,jk) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * usd(ji  ,jj,jk)    &
-                  &                 - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * usd(ji-1,jj,jk)    &
-                  &                 + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vsd(ji,jj  ,jk)    &
-                  &                 - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vsd(ji,jj-1,jk)  ) * r1_e1e2t(ji,jj)
-            END DO
-         END DO
-      END DO
+      DO_3D_01_01( 1, jpkm1 )
+         ze3divh(ji,jj,jk) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * usd(ji  ,jj,jk)    &
+            &                 - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * usd(ji-1,jj,jk)    &
+            &                 + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vsd(ji,jj  ,jk)    &
+            &                 - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vsd(ji,jj-1,jk)  ) * r1_e1e2t(ji,jj)
+      END_3D
       !
 #if defined key_agrif
@@ -291,14 +270,12 @@
       !
       IF( ln_tauw ) THEN
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               ! Stress components at u- & v-points
-               utau(ji,jj) = 0.5_wp * ( tauw_x(ji,jj) + tauw_x(ji+1,jj) )
-               vtau(ji,jj) = 0.5_wp * ( tauw_y(ji,jj) + tauw_y(ji,jj+1) )
-               !
-               ! Stress module at t points
-               taum(ji,jj) = SQRT( tauw_x(ji,jj)*tauw_x(ji,jj) + tauw_y(ji,jj)*tauw_y(ji,jj) )
-            END DO
-         END DO
+         DO_2D_10_10
+            ! Stress components at u- & v-points
+            utau(ji,jj) = 0.5_wp * ( tauw_x(ji,jj) + tauw_x(ji+1,jj) )
+            vtau(ji,jj) = 0.5_wp * ( tauw_y(ji,jj) + tauw_y(ji,jj+1) )
+            !
+            ! Stress module at t points
+            taum(ji,jj) = SQRT( tauw_x(ji,jj)*tauw_x(ji,jj) + tauw_y(ji,jj)*tauw_y(ji,jj) )
+         END_2D
          CALL lbc_lnk_multi( 'sbcwave', utau(:,:), 'U', -1. , vtau(:,:), 'V', -1. , taum(:,:) , 'T', -1. )
       ENDIF
@@ -307,5 +284,5 @@
 
 
-   SUBROUTINE sbc_wave( kt )
+   SUBROUTINE sbc_wave( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE sbc_wave  ***
@@ -322,4 +299,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in   ) ::   kt   ! ocean time step
+      INTEGER, INTENT(in   ) ::   Kmm  ! ocean time index
       !!---------------------------------------------------------------------
       !
@@ -361,5 +339,5 @@
          !
          IF( ( ll_st_bv_li   .AND. jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0 ) .OR. &
-           & ( ll_st_peakfr  .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0                ) ) CALL sbc_stokes()
+           & ( ll_st_peakfr  .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0                ) ) CALL sbc_stokes( Kmm )
          !
       ENDIF
@@ -395,9 +373,7 @@
       !!---------------------------------------------------------------------
       !
-      REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model
       READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_wave in reference namelist' )
          
-      REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model
       READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist' )
