Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRA/traadv_mus.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRA/traadv_mus.F90	(revision 12529)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRA/traadv_mus.F90	(revision 13151)
@@ -29,6 +29,6 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! distribued memory computing
-   USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
 
    IMPLICIT NONE
@@ -36,9 +36,9 @@
 
    PUBLIC   tra_adv_mus   ! routine called by traadv.F90
-   
+
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   upsmsk   !: mixed upstream/centered scheme near some straits
    !                                                           !  and in closed seas (orca 2 and 1 configurations)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index
-   
+
    LOGICAL  ::   l_trd   ! flag to compute trends
    LOGICAL  ::   l_ptr   ! flag to compute poleward transport
@@ -47,7 +47,8 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$ 
+   !! $Id$
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -64,5 +65,5 @@
       !!
       !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
-      !!              ld_msc_ups=T : 
+      !!              ld_msc_ups=T :
       !!
       !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
@@ -88,5 +89,5 @@
       REAL(wp) ::   zv, z0v, zzwy, z0w           !   -      -
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zslpx   ! 3D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zslpy   ! -      - 
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zslpy   ! -      -
       !!----------------------------------------------------------------------
       !
@@ -112,13 +113,13 @@
                   &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 in some user defined area
             END DO
-         ENDIF 
-         !
-      ENDIF 
-      !      
+         ENDIF
+         !
+      ENDIF
+      !
       l_trd = .FALSE.
       l_hst = .FALSE.
       l_ptr = .FALSE.
       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
-      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE. 
+      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE.
       IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
          &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
@@ -130,5 +131,5 @@
          !                                !-- first guess of the slopes
          zwx(:,:,jpk) = 0._wp                   ! bottom values
-         zwy(:,:,jpk) = 0._wp  
+         zwy(:,:,jpk) = 0._wp
          DO_3D_10_10( 1, jpkm1 )
             zwx(ji,jj,jk) = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
@@ -176,6 +177,6 @@
          DO_3D_00_00( 1, jpkm1 )
             pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       &
-            &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     &
-            &                                   * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+            &                                           +   zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     &
+            &                           * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
          END_3D
          !                                ! trend diagnostics
@@ -184,5 +185,5 @@
             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kbb) )
          END IF
-         !                                 ! "Poleward" heat and salt transports 
+         !                                 ! "Poleward" heat and salt transports
          IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
          !                                 !  heat transport
@@ -227,5 +228,6 @@
          !
          DO_3D_00_00( 1, jpkm1 )
-            pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(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) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )   &
+               &                                      * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
          END_3D
          !                                ! send trends for diagnostic
