Index: NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traadv_ubs.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traadv_ubs.F90	(revision 12451)
+++ NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traadv_ubs.F90	(revision 13540)
@@ -39,4 +39,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -82,5 +83,5 @@
       !!
       !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404. 
-      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741. 
+      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731�1741.
       !!----------------------------------------------------------------------
       INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
@@ -123,6 +124,6 @@
          !                                                       ! ===========
          !                                              
-         DO jk = 1, jpkm1        !==  horizontal laplacian of before tracer ==!
-            DO_2D_10_10
+         DO jk = 1, jpkm1                !==  horizontal laplacian of before tracer ==!
+            DO_2D( 1, 0, 1, 0 )                   ! First derivative (masked gradient)
                zeeu = e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
                zeev = e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
@@ -130,5 +131,5 @@
                ztv(ji,jj,jk) = zeev * ( pt(ji  ,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
             END_2D
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )                   ! Second derivative (divergence)
                zcoef = 1._wp / ( 6._wp * e3t(ji,jj,jk,Kmm) )
                zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
@@ -137,8 +138,8 @@
             !                                    
          END DO         
-         CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1. )   ;    CALL lbc_lnk( 'traadv_ubs', zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
+         CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1.0_wp )   ;    CALL lbc_lnk( 'traadv_ubs', zltv, 'T', 1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
          !    
-         DO_3D_10_10( 1, jpkm1 )
-            zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )      ! upstream transport (x2)
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )   !==  Horizontal advective fluxes  ==!     (UBS)
+            zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )        ! upstream transport (x2)
             zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) )
             zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) )
@@ -155,8 +156,9 @@
          !
          DO jk = 1, jpkm1        !==  add the horizontal advective trend  ==!
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)                        &
                   &             - (  ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk)    &
-                  &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+                  &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) &
+                  &                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
             END_2D
             !                                             
@@ -164,5 +166,5 @@
          !
          zltu(:,:,:) = pt(:,:,:,jn,Krhs) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case
-         !                                            ! and/or in trend diagnostic (l_trd=T) 
+         !                                                ! and/or in trend diagnostic (l_trd=T) 
          !                
          IF( l_trd ) THEN                  ! trend diagnostics
@@ -185,29 +187,30 @@
             IF( l_trd )   zltv(:,:,:) = pt(:,:,:,jn,Krhs)          ! store pt(:,:,:,:,Krhs) if trend diag.
             !
-            !                          !*  upstream advection with initial mass fluxes & intermediate update  ==!
-            DO_3D_11_11( 2, jpkm1 )
+            !                               !*  upstream advection with initial mass fluxes & intermediate update  ==!
+            DO_3D( 1, 1, 1, 1, 2, jpkm1 )
                zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
                zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) )
                ztw(ji,jj,jk) = 0.5_wp * (  zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb)  ) * wmask(ji,jj,jk)
             END_3D
-            IF( ln_linssh ) THEN             ! top ocean value (only in linear free surface as ztw has been w-masked)
-               IF( ln_isfcav ) THEN                ! top of the ice-shelf cavities and at the ocean surface
-                  DO_2D_11_11
+            IF( ln_linssh ) THEN                ! top ocean value (only in linear free surface as ztw has been w-masked)
+               IF( ln_isfcav ) THEN                   ! top of the ice-shelf cavities and at the ocean surface
+                  DO_2D( 1, 1, 1, 1 )
                      ztw(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
                   END_2D
-               ELSE                                ! no cavities: only at the ocean surface
+               ELSE                                   ! no cavities: only at the ocean surface
                   ztw(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
                ENDIF
             ENDIF
             !
-            DO_3D_00_00( 1, jpkm1 )
-               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !* trend and after field with monotonic scheme
+               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )    &
+                  &     * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                pt(ji,jj,jk,jn,Krhs) =   pt(ji,jj,jk,jn,Krhs) +  ztak 
                zti(ji,jj,jk)    = ( pt(ji,jj,jk,jn,Kbb) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
             END_3D
-            CALL lbc_lnk( 'traadv_ubs', zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
+            CALL lbc_lnk( 'traadv_ubs', zti, 'T', 1.0_wp )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
             !
             !                          !*  anti-diffusive flux : high order minus low order
-            DO_3D_11_11( 2, jpkm1 )
+            DO_3D( 1, 1, 1, 1, 2, jpkm1 )
                ztw(ji,jj,jk) = (   0.5_wp * pW(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
                   &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk)
@@ -220,5 +223,5 @@
          CASE(  4  )                               ! 4th order COMPACT
             CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )         ! 4th order compact interpolation of T at w-point
-            DO_3D_00_00( 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                ztw(ji,jj,jk) = pW(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
             END_3D
@@ -227,10 +230,11 @@
          END SELECT
          !
-         DO_3D_00_00( 1, jpkm1 )
-            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !  final trend with corrected fluxes
+            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )    &
+               &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
          END_3D
          !
-         IF( l_trd )  THEN       ! vertical advective trend diagnostics
-            DO_3D_00_00( 1, jpkm1 )
+         IF( l_trd )  THEN               ! vertical advective trend diagnostics
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )                 ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
                zltv(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltv(ji,jj,jk)                          &
                   &           + pt(ji,jj,jk,jn,Kmm) * (  pW(ji,jj,jk) - pW(ji,jj,jk+1)  )   &
@@ -270,5 +274,5 @@
       !!----------------------------------------------------------------------
       !
-      zbig  = 1.e+40_wp
+      zbig  = 1.e+38_wp
       zrtrn = 1.e-15_wp
       zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
@@ -282,5 +286,5 @@
       DO jk = 1, jpkm1     ! search maximum in neighbourhood
          ikm1 = MAX(jk-1,1)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
                &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
@@ -294,5 +298,5 @@
       DO jk = 1, jpkm1     ! search minimum in neighbourhood
          ikm1 = MAX(jk-1,1)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
                &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
@@ -306,5 +310,5 @@
       ! Positive and negative part of fluxes and beta terms
       ! ---------------------------------------------------
-      DO_3D_00_00( 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          ! positive & negative part of the flux
          zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
@@ -318,8 +322,8 @@
       ! monotonic flux in the k direction, i.e. pcc
       ! -------------------------------------------
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
          zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
-         zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
+         zc = 0.5 * ( 1.e0 + SIGN( 1.0_wp, pcc(ji,jj,jk) ) )
          pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
       END_3D
