Index: NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/TOP/TRP/trcsbc.F90
===================================================================
--- NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/TOP/TRP/trcsbc.F90	(revision 11396)
+++ NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/TOP/TRP/trcsbc.F90	(revision 13463)
@@ -18,5 +18,5 @@
    USE oce_trc         ! ocean dynamics and active tracers variables
    USE trc             ! ocean  passive tracers variables
-   USE prtctl_trc      ! Print control for debbuging
+   USE prtctl          ! Print control for debbuging
    USE iom
    USE trd_oce
@@ -29,5 +29,6 @@
 
    !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -37,5 +38,5 @@
 CONTAINS
 
-   SUBROUTINE trc_sbc ( kt )
+   SUBROUTINE trc_sbc ( kt, Kmm, ptr, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trc_sbc  ***
@@ -49,14 +50,16 @@
       !!            The surface freshwater flux modify the ocean volume
       !!         and thus the concentration of a tracer as :
-      !!            tra = tra + emp * trn / e3t   for k=1
+      !!            tr(Krhs) = tr(Krhs) + emp * tr(Kmm) / e3t_   for k=1
       !!         where emp, the surface freshwater budget (evaporation minus
       !!         precipitation ) given in kg/m2/s is divided
       !!         by 1035 kg/m3 (density of ocean water) to obtain m/s.
       !!
-      !! ** Action  : - Update the 1st level of tra with the trend associated
+      !! ** Action  : - Update the 1st level of tr(:,:,:,:,Krhs) with the trend associated
       !!                with the tracer surface boundary condition 
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER,                                    INTENT(in   ) :: kt        ! ocean time-step index
+      INTEGER,                                    INTENT(in   ) :: Kmm, Krhs ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr       ! passive tracers and RHS of tracer equation
       !
       INTEGER  ::   ji, jj, jn                      ! dummy loop indices
@@ -82,8 +85,8 @@
          IF( ln_rsttr .AND. .NOT.ln_top_euler .AND.   &                     ! Restart: read in restart  file
             iom_varid( numrtr, 'sbc_'//TRIM(ctrcnm(1))//'_b', ldstop = .FALSE. ) > 0 ) THEN
-            IF(lwp) WRITE(numout,*) '          nittrc000-nn_dttrc surface tracer content forcing fields red in the restart file'
+            IF(lwp) WRITE(numout,*) '          nittrc000-1 surface tracer content forcing fields read in the restart file'
             zfact = 0.5_wp
             DO jn = 1, jptra
-               CALL iom_get( numrtr, jpdom_autoglo, 'sbc_'//TRIM(ctrcnm(jn))//'_b', sbc_trc_b(:,:,jn) )   ! before tracer content sbc
+               CALL iom_get( numrtr, jpdom_auto, 'sbc_'//TRIM(ctrcnm(jn))//'_b', sbc_trc_b(:,:,jn) )   ! before tracer content sbc
             END DO
          ELSE                                         ! No restart or restart not found: Euler forward time stepping
@@ -102,5 +105,5 @@
       ENDIF
 
-      ! Coupling online : river runoff is added to the horizontal divergence (hdivn) in the subroutine sbc_rnf_div 
+      ! Coupling online : river runoff is added to the horizontal divergence (hdiv) in the subroutine sbc_rnf_div 
       ! one only consider the concentration/dilution effect due to evaporation minus precipitation + freezing/melting of sea-ice
       ! Coupling offline : runoff are in emp which contains E-P-R
@@ -118,9 +121,7 @@
          !
          DO jn = 1, jptra
-            DO jj = 2, jpj
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  sbc_trc(ji,jj,jn) = zsfx(ji,jj) * r1_rau0 * trn(ji,jj,1,jn)
-               END DO
-            END DO
+            DO_2D( 0, 1, 0, 0 )
+               sbc_trc(ji,jj,jn) = zsfx(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm)
+            END_2D
          END DO
          !
@@ -128,9 +129,7 @@
          !
          DO jn = 1, jptra
-            DO jj = 2, jpj
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  sbc_trc(ji,jj,jn) = ( zsfx(ji,jj) + fmmflx(ji,jj) ) * r1_rau0 * trn(ji,jj,1,jn)
-               END DO
-            END DO
+            DO_2D( 0, 1, 0, 0 )
+               sbc_trc(ji,jj,jn) = ( zsfx(ji,jj) + fmmflx(ji,jj) ) * r1_rho0 * ptr(ji,jj,1,jn,Kmm)
+            END_2D
          END DO
          !
@@ -138,40 +137,36 @@
          !
          DO jn = 1, jptra
-            DO jj = 2, jpj
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zse3t = 1. / e3t_n(ji,jj,1)
-                  ! tracer flux at the ice/ocean interface (tracer/m2/s)
-                  zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice
-                  !                                         ! only used in the levitating sea ice case
-                  ! tracer flux only       : add concentration dilution term in net tracer flux, no F-M in volume flux
-                  ! tracer and mass fluxes : no concentration dilution term in net tracer flux, F-M term in volume flux
-                  ztfx  = zftra                        ! net tracer flux
-                  !
-                  zdtra = r1_rau0 * ( ztfx + ( zsfx(ji,jj) + fmmflx(ji,jj) ) * trn(ji,jj,1,jn) ) 
-                  IF ( zdtra < 0. ) THEN
-                     zdtra  = MAX(zdtra, -trn(ji,jj,1,jn) * e3t_n(ji,jj,1) / r2dttrc )   ! avoid negative concentrations to arise
-                  ENDIF
-                  sbc_trc(ji,jj,jn) =  zdtra 
-               END DO
-            END DO
+            DO_2D( 0, 1, 0, 0 )
+               zse3t = 1. / e3t(ji,jj,1,Kmm)
+               ! tracer flux at the ice/ocean interface (tracer/m2/s)
+               zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice
+               !                                         ! only used in the levitating sea ice case
+               ! tracer flux only       : add concentration dilution term in net tracer flux, no F-M in volume flux
+               ! tracer and mass fluxes : no concentration dilution term in net tracer flux, F-M term in volume flux
+               ztfx  = zftra                        ! net tracer flux
+               !
+               zdtra = r1_rho0 * ( ztfx + ( zsfx(ji,jj) + fmmflx(ji,jj) ) * ptr(ji,jj,1,jn,Kmm) ) 
+               IF ( zdtra < 0. ) THEN
+                  zdtra  = MAX(zdtra, -ptr(ji,jj,1,jn,Kmm) * e3t(ji,jj,1,Kmm) / rDt_trc )   ! avoid negative concentrations to arise
+               ENDIF
+               sbc_trc(ji,jj,jn) =  zdtra 
+            END_2D
          END DO
       END SELECT
       !
-      CALL lbc_lnk( 'trcsbc', sbc_trc(:,:,:), 'T', 1. )
+      CALL lbc_lnk( 'trcsbc', sbc_trc(:,:,:), 'T', 1.0_wp )
       !                                       Concentration dilution effect on tracers due to evaporation & precipitation 
       DO jn = 1, jptra
          !
-         IF( l_trdtrc )   ztrtrd(:,:,:) = tra(:,:,:,jn)  ! save trends
-         !
-         DO jj = 2, jpj
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zse3t = zfact / e3t_n(ji,jj,1)
-               tra(ji,jj,1,jn) = tra(ji,jj,1,jn) + ( sbc_trc_b(ji,jj,jn) + sbc_trc(ji,jj,jn) ) * zse3t
-            END DO
-         END DO
+         IF( l_trdtrc )   ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs)  ! save trends
+         !
+         DO_2D( 0, 1, 0, 0 )
+            zse3t = zfact / e3t(ji,jj,1,Kmm)
+            ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + ( sbc_trc_b(ji,jj,jn) + sbc_trc(ji,jj,jn) ) * zse3t
+         END_2D
          !
          IF( l_trdtrc ) THEN
-            ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:)
-            CALL trd_tra( kt, 'TRC', jn, jptra_nsr, ztrtrd )
+            ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs) - ztrtrd(:,:,:)
+            CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_nsr, ztrtrd )
          END IF
          !                                                       ! ===========
@@ -191,7 +186,7 @@
       ENDIF
       !
-      IF( ln_ctl )   THEN
-         WRITE(charout, FMT="('sbc ')") ;  CALL prt_ctl_trc_info(charout)
-                                           CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+      IF( sn_cfctl%l_prttrc )   THEN
+         WRITE(charout, FMT="('sbc ')") ;  CALL prt_ctl_info( charout, cdcomp = 'top' )
+                                           CALL prt_ctl( tab4d_1=ptr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm, clinfo3='trd' )
       ENDIF
       IF( l_trdtrc )  DEALLOCATE( ztrtrd )
@@ -205,7 +200,11 @@
    !!   Dummy module :                      NO passive tracer
    !!----------------------------------------------------------------------
+   USE par_oce
+   USE par_trc
 CONTAINS
-   SUBROUTINE trc_sbc (kt)              ! Empty routine
-      INTEGER, INTENT(in) :: kt
+   SUBROUTINE trc_sbc ( kt, Kmm, ptr, Krhs )      ! Empty routine
+      INTEGER,                                    INTENT(in   ) :: kt        ! ocean time-step index
+      INTEGER,                                    INTENT(in   ) :: Kmm, Krhs ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr       ! passive tracers and RHS of tracer equation
       WRITE(*,*) 'trc_sbc: You should not have seen this print! error?', kt
    END SUBROUTINE trc_sbc
