Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trddyn.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trddyn.F90	(revision 12529)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trddyn.F90	(revision 13151)
@@ -37,4 +37,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdglo.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdglo.F90	(revision 12529)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdglo.F90	(revision 13151)
@@ -52,4 +52,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -116,7 +117,7 @@
             DO_3D_10_10( 1, jpkm1 )
                zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
-                  &                                     * e1e2u  (ji,jj) * e3u(ji,jj,jk,Kmm)
+                  &                  * e1e2u  (ji  ,jj) * e3u(ji,jj,jk,Kmm)
                zvs = ptrdy(ji,jj,jk) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
-                  &                                     * e1e2v  (ji,jj) * e3u(ji,jj,jk,Kmm)
+                  &                  * e1e2v  (ji,jj  ) * e3u(ji,jj,jk,Kmm)
                umo(ktrd) = umo(ktrd) + zvt
                vmo(ktrd) = vmo(ktrd) + zvs
@@ -211,6 +212,8 @@
          zcof   = 0.5_wp / rho0           ! Density flux at u and v-points
          DO_3D_10_10( 1, jpkm1 )
-            zkx(ji,jj,jk) = zcof * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
-            zky(ji,jj,jk) = zcof * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
+            zkx(ji,jj,jk) = zcof * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)   &
+               &                              *  uu(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
+            zky(ji,jj,jk) = zcof * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)   &
+               &                              *  vv(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
          END_3D
          
@@ -219,5 +222,5 @@
                &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
                &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
-               &              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
+               &     / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
          END_3D
 
@@ -226,5 +229,6 @@
          peke = 0._wp
          DO jk = 1, jpkm1
-            peke = peke + SUM( zkepe(:,:,jk) * gdept(:,:,jk,Kmm) * e1e2t(:,:) * e3t(:,:,jk,Kmm) )
+            peke = peke + SUM( zkepe(:,:,jk) * gdept(:,:,jk,Kmm) * e1e2t(:,:)   &
+               &                               * e3t(:,:,jk,Kmm) )
          END DO
          peke = grav * peke
@@ -524,6 +528,8 @@
 
       DO_3D_00_00( 1, jpk )
-         tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
-         tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
+         tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)   &
+            &                                       * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
+         tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * e3v(ji,jj,jk,Kmm)   &
+            &                                       * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
       END_3D
       CALL mpp_sum( 'trdglo', tvolu )   ! sums over the global domain
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdken.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdken.F90	(revision 12529)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdken.F90	(revision 13151)
@@ -41,4 +41,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdmxl.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdmxl.F90	(revision 12529)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdmxl.F90	(revision 13151)
@@ -70,4 +70,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -120,5 +121,7 @@
          wkx(:,:,:) = 0._wp         !==  now ML weights for vertical averaging  ==!
          DO_3D_11_11( 1, jpktrd )
-            IF( jk - kmxln(ji,jj) < 0 )   wkx(ji,jj,jk) = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
+            IF( jk - kmxln(ji,jj) < 0 )   THEN
+               wkx(ji,jj,jk) = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
+            ENDIF
          END_3D
          hmxl(:,:) = 0._wp               ! NOW mixed-layer depth
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdpen.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdpen.F90	(revision 12529)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdpen.F90	(revision 13151)
@@ -35,4 +35,6 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   rab_pe   ! partial derivatives of PE anomaly with respect to T and S
 
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -40,4 +42,5 @@
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
+
 CONTAINS
 
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdtra.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdtra.F90	(revision 12529)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdtra.F90	(revision 13151)
@@ -42,4 +42,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -128,6 +129,8 @@
             zwt(:,:,jpk) = 0._wp   ;   zws(:,:,jpk) = 0._wp
             DO jk = 2, jpk
-               zwt(:,:,jk) = avt(:,:,jk) * ( ts(:,:,jk-1,jp_tem,Krhs) - ts(:,:,jk,jp_tem,Krhs) ) / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
-               zws(:,:,jk) = avs(:,:,jk) * ( ts(:,:,jk-1,jp_sal,Krhs) - ts(:,:,jk,jp_sal,Krhs) ) / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
+               zwt(:,:,jk) = avt(:,:,jk) * ( ts(:,:,jk-1,jp_tem,Krhs) - ts(:,:,jk,jp_tem,Krhs) )   &
+                  &        / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
+               zws(:,:,jk) = avs(:,:,jk) * ( ts(:,:,jk-1,jp_sal,Krhs) - ts(:,:,jk,jp_sal,Krhs) )   &
+                  &        / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
             END DO
             !
@@ -142,6 +145,8 @@
             zwt(:,:,:) = 0._wp   ;   zws(:,:,:) = 0._wp            ! vertical diffusive fluxes
             DO jk = 2, jpk
-               zwt(:,:,jk) = avt_evd(:,:,jk) * ( ts(:,:,jk-1,jp_tem,Krhs) - ts(:,:,jk,jp_tem,Krhs) ) / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
-               zws(:,:,jk) = avt_evd(:,:,jk) * ( ts(:,:,jk-1,jp_sal,Krhs) - ts(:,:,jk,jp_sal,Krhs) ) / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
+               zwt(:,:,jk) = avt_evd(:,:,jk) * ( ts(:,:,jk-1,jp_tem,Krhs) - ts(:,:,jk,jp_tem,Krhs) )   &
+                  &            / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
+               zws(:,:,jk) = avt_evd(:,:,jk) * ( ts(:,:,jk-1,jp_sal,Krhs) - ts(:,:,jk,jp_sal,Krhs) )   &
+                  &            / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
             END DO
             !
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdvor.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdvor.F90	(revision 12529)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRD/trdvor.F90	(revision 13151)
@@ -57,4 +57,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -192,5 +193,6 @@
          DO jj = 1, jpjm1
             vortrd(ji,jj,ktrd) = (    zvdpvor(ji+1,jj) - zvdpvor(ji,jj)       &
-                 &                - ( zudpvor(ji,jj+1) - zudpvor(ji,jj) )   ) / ( e1f(ji,jj) * e2f(ji,jj) )
+                 &                - ( zudpvor(ji,jj+1) - zudpvor(ji,jj) )   ) &
+                 &                  / ( e1f(ji,jj) * e2f(ji,jj) )
          END DO
       END DO
@@ -268,5 +270,6 @@
             DO jj = 1, jpjm1
                vortrd(ji,jj,jpvor_bev) = (    zvbet(ji+1,jj) - zvbet(ji,jj)     &
-                  &                       - ( zubet(ji,jj+1) - zubet(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) )
+                  &                       - ( zubet(ji,jj+1) - zubet(ji,jj) ) ) &
+                  &                           / ( e1f(ji,jj) * e2f(ji,jj) )
             END DO
          END DO
@@ -283,5 +286,6 @@
          DO jj=1,jpjm1
             vortrd(ji,jj,ktrd) = (    zvdpvor(ji+1,jj) - zvdpvor(ji,jj)     &
-               &                  - ( zudpvor(ji,jj+1) - zudpvor(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) )
+               &                  - ( zudpvor(ji,jj+1) - zudpvor(ji,jj) ) ) &
+               &                         / ( e1f(ji,jj) * e2f(ji,jj) )
          END DO
       END DO
@@ -345,5 +349,6 @@
          DO jj = 1, jpjm1
             vor_avr(ji,jj) = (  ( zvv(ji+1,jj) - zvv(ji,jj) )    &
-               &              - ( zuu(ji,jj+1) - zuu(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,1)
+               &              - ( zuu(ji,jj+1) - zuu(ji,jj) ) )  &
+               &             / ( e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,1)
          END DO
       END DO
