Index: NEMO/trunk/src/OCE/USR/usrdef_sbc.F90
===================================================================
--- NEMO/trunk/src/OCE/USR/usrdef_sbc.F90	(revision 10425)
+++ NEMO/trunk/src/OCE/USR/usrdef_sbc.F90	(revision 12377)
@@ -31,5 +31,5 @@
 
    !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -39,5 +39,5 @@
 CONTAINS
 
-   SUBROUTINE usrdef_sbc_oce( kt )
+   SUBROUTINE usrdef_sbc_oce( kt, Kbb )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE usrdef_sbc  ***
@@ -55,4 +55,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
+      INTEGER, INTENT(in) ::   Kbb  ! ocean time index
       !!
       INTEGER  ::   ji, jj                 ! dummy loop indices
@@ -109,26 +110,24 @@
       ztrp= - 40.e0        ! retroaction term on heat fluxes (W/m2/K)
       zconv = 3.16e-5      ! convertion factor: 1 m/yr => 3.16e-5 mm/s
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ! domain from 15 deg to 50 deg between 27 and 28  degC at 15N, -3
-            ! and 13 degC at 50N 53.5 + or - 11 = 1/4 period :
-            ! 64.5 in summer, 42.5 in winter
-            t_star = zTstar * ( 1. + 1. / 50. * zcos_sais2 )                &
-               &                    * COS( rpi * (gphit(ji,jj) - 5.)               &
-               &                    / ( 53.5 * ( 1 + 11 / 53.5 * zcos_sais2 ) * 2.) )
-            ! 23.5 deg : tropics
-            qsr (ji,jj) =  230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais1 ) / ( 0.9 * 180 ) )
-            qns (ji,jj) = ztrp * ( tsb(ji,jj,1,jp_tem) - t_star ) - qsr(ji,jj)
-            IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj) ) THEN    ! zero at 37.8 deg, max at 24.6 deg
-               emp  (ji,jj) =   zemp_S * zconv   &
-                  &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) )  &
-                  &         * ( 1 - zemp_sais / zemp_S * zcos_sais1)
-            ELSE
-               emp (ji,jj) =  - zemp_N * zconv   &
-                  &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) )  &
-                  &         * ( 1 - zemp_sais / zemp_N * zcos_sais1 )
-            ENDIF
-         END DO
-      END DO
+      DO_2D_11_11
+         ! domain from 15 deg to 50 deg between 27 and 28  degC at 15N, -3
+         ! and 13 degC at 50N 53.5 + or - 11 = 1/4 period :
+         ! 64.5 in summer, 42.5 in winter
+         t_star = zTstar * ( 1. + 1. / 50. * zcos_sais2 )                &
+            &                    * COS( rpi * (gphit(ji,jj) - 5.)               &
+            &                    / ( 53.5 * ( 1 + 11 / 53.5 * zcos_sais2 ) * 2.) )
+         ! 23.5 deg : tropics
+         qsr (ji,jj) =  230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais1 ) / ( 0.9 * 180 ) )
+         qns (ji,jj) = ztrp * ( ts(ji,jj,1,jp_tem,Kbb) - t_star ) - qsr(ji,jj)
+         IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj) ) THEN    ! zero at 37.8 deg, max at 24.6 deg
+            emp  (ji,jj) =   zemp_S * zconv   &
+               &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) )  &
+               &         * ( 1 - zemp_sais / zemp_S * zcos_sais1)
+         ELSE
+            emp (ji,jj) =  - zemp_N * zconv   &
+               &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) )  &
+               &         * ( 1 - zemp_sais / zemp_N * zcos_sais1 )
+         ENDIF
+      END_2D
 
       zsumemp = GLOB_SUM( 'usrdef_sbc', emp  (:,:)   ) 
@@ -166,24 +165,20 @@
       ztau_sais = 0.015
       ztaun = ztau - ztau_sais * COS( (ztime - ztimemax) / (ztimemin - ztimemax) * rpi )
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-           ! domain from 15deg to 50deg and 1/2 period along 14deg
-           ! so 5/4 of half period with seasonal cycle
-           utau(ji,jj) = - ztaun * SIN( rpi * (gphiu(ji,jj) - 15.) / (29.-15.) )
-           vtau(ji,jj) =   ztaun * SIN( rpi * (gphiv(ji,jj) - 15.) / (29.-15.) )
-         END DO
-      END DO
+      DO_2D_11_11
+        ! domain from 15deg to 50deg and 1/2 period along 14deg
+        ! so 5/4 of half period with seasonal cycle
+        utau(ji,jj) = - ztaun * SIN( rpi * (gphiu(ji,jj) - 15.) / (29.-15.) )
+        vtau(ji,jj) =   ztaun * SIN( rpi * (gphiv(ji,jj) - 15.) / (29.-15.) )
+      END_2D
 
       ! module of wind stress and wind speed at T-point
       zcoef = 1. / ( zrhoa * zcdrag ) 
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vect. opt.
-            ztx = utau(ji-1,jj  ) + utau(ji,jj) 
-            zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
-            zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
-            taum(ji,jj) = zmod
-            wndm(ji,jj) = SQRT( zmod * zcoef )
-         END DO
-      END DO
+      DO_2D_00_00
+         ztx = utau(ji-1,jj  ) + utau(ji,jj) 
+         zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
+         zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
+         taum(ji,jj) = zmod
+         wndm(ji,jj) = SQRT( zmod * zcoef )
+      END_2D
       CALL lbc_lnk_multi( 'usrdef_sbc', taum(:,:), 'T', 1. , wndm(:,:), 'T', 1. )
 
