Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynldf_iso.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynldf_iso.F90	(revision 12590)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynldf_iso.F90	(revision 12606)
@@ -22,5 +22,5 @@
    USE ldftra          ! lateral physics: eddy diffusivity
    USE zdf_oce         ! ocean vertical physics
-   USE ldfslp          ! iso-neutral slopes 
+   USE ldfslp          ! iso-neutral slopes
    !
    USE in_out_manager  ! I/O manager
@@ -36,10 +36,11 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akzu, akzv   !: vertical component of rotated lateral viscosity
-   
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u   ! 2D workspace (dyn_ldf_iso) 
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u   ! 2D workspace (dyn_ldf_iso)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v   !  -      -
 
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -53,5 +54,5 @@
       !!                  ***  ROUTINE dyn_ldf_iso_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( akzu(jpi,jpj,jpk) , zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) ,     & 
+      ALLOCATE( akzu(jpi,jpj,jpk) , zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) ,     &
          &      akzv(jpi,jpj,jpk) , zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc )
          !
@@ -63,5 +64,5 @@
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE dyn_ldf_iso  ***
-      !!                       
+      !!
       !! ** Purpose :   Compute the before trend of the rotated laplacian
       !!      operator of lateral momentum diffusion except the diagonal
@@ -137,7 +138,7 @@
          !
        ENDIF
-         
+
       zaht_0 = 0.5_wp * rn_Ud * rn_Ld                  ! aht_0 from namtra_ldf = zaht_max
-      
+
       !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -161,12 +162,14 @@
 
          !                               -----f-----
-         ! Horizontal fluxes on U             |  
+         ! Horizontal fluxes on U             |
          ! --------------------===        t   u   t
-         !                                    |  
+         !                                    |
          ! i-flux at t-point             -----f-----
 
          IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
             DO_2D_00_01
-               zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( e3u(ji,jj,jk,Kmm), e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj)
+               zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj)   &
+                  &    * MIN( e3u(ji  ,jj,jk,Kmm),                &
+                  &           e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj)
 
                zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)     &
@@ -181,5 +184,6 @@
          ELSE                   ! other coordinate system (zco or sco) : e3t
             DO_2D_00_01
-               zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj)
+               zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b )   &
+                  &     * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj)
 
                zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  ) + umask(ji,jj,jk+1)     &
@@ -196,5 +200,6 @@
          ! j-flux at f-point
          DO_2D_10_10
-            zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj)
+            zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b )   &
+               &     * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj)
 
             zmskf = 1._wp / MAX(   umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)     &
@@ -215,5 +220,6 @@
 
          DO_2D_00_10
-            zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj)
+            zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b )   &
+               &     * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj)
 
             zmskf = 1._wp / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)     &
@@ -230,5 +236,7 @@
          IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
             DO_2D_01_10
-               zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( e3v(ji,jj,jk,Kmm), e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj)
+               zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj)   &
+                  &     * MIN( e3v(ji,jj  ,jk,Kmm),                 &
+                  &            e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj)
 
                zmskt = 1._wp / MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)     &
@@ -243,5 +251,6 @@
          ELSE                   ! other coordinate system (zco or sco) : e3t
             DO_2D_01_10
-               zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj)
+               zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b )   &
+                  &     * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj)
 
                zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
@@ -261,7 +270,9 @@
          DO_2D_00_00
             puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (  ziut(ji+1,jj) - ziut(ji,jj  )    &
-               &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
+               &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj)   &
+               &                           / e3u(ji,jj,jk,Kmm)
             pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (  zivf(ji,jj  ) - zivf(ji-1,jj)    &
-               &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
+               &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj)   &
+               &                           / e3v(ji,jj,jk,Kmm)
          END_2D
          !                                             ! ===============
@@ -278,5 +289,5 @@
          !                                             ! ===============
 
- 
+
          ! I. vertical trends associated with the lateral mixing
          ! =====================================================
@@ -375,6 +386,8 @@
          DO jk = 1, jpkm1
             DO ji = 2, jpim1
-               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
-               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
+               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj)   &
+                  &               / e3u(ji,jj,jk,Kmm)
+               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj)   &
+                  &               / e3v(ji,jj,jk,Kmm)
             END DO
          END DO
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynldf_lap_blp.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynldf_lap_blp.F90	(revision 12590)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynldf_lap_blp.F90	(revision 12606)
@@ -14,5 +14,5 @@
    USE dom_oce        ! ocean space and time domain
    USE ldfdyn         ! lateral diffusion: eddy viscosity coef.
-   USE ldfslp         ! iso-neutral slopes 
+   USE ldfslp         ! iso-neutral slopes
    USE zdf_oce        ! ocean vertical physics
    !
@@ -28,7 +28,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)
    !!----------------------------------------------------------------------
@@ -38,10 +39,10 @@
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE dyn_ldf_lap  ***
-      !!                       
-      !! ** Purpose :   Compute the before horizontal momentum diffusive 
+      !!
+      !! ** Purpose :   Compute the before horizontal momentum diffusive
       !!      trend and add it to the general trend of momentum equation.
       !!
-      !! ** Method  :   The Laplacian operator apply on horizontal velocity is 
-      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) 
+      !! ** Method  :   The Laplacian operator apply on horizontal velocity is
+      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )
       !!
       !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv.
@@ -76,22 +77,28 @@
 !!gm open question here : e3f  at before or now ?    probably now...
 !!gm note that ahmf has already been multiplied by fmask
-            zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &
-               &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
-               &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
+            zcur(ji-1,jj-1) =  &
+               &      ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)      &
+               &  * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
+               &     - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
             !                                      ! ahm * div        (computed from 2 to jpi/jpj)
 !!gm note that ahmt has already been multiplied by tmask
-            zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)                                         &
-               &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  &
-               &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  )
+            zdiv(ji,jj)     =   &
+               &   ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)      &
+               &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk)        &
+               &        - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  &
+               &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk)        &
+               &        - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  )
          END_2D
          !
          DO_2D_00_00
-            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * (                                                 &
-               &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
-               &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     )
+            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * (                             &
+               &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj)   &
+               &              / e3u(ji,jj,jk,Kmm)   &
+               &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)        )
                !
-            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * (                                                 &
-               &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   &
-               &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     )
+            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * (                              &
+               &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj)   &
+               &              / e3v(ji,jj,jk,Kmm)   &
+               &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)       )
          END_2D
          !                                             ! ===============
@@ -105,6 +112,6 @@
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE dyn_ldf_blp  ***
-      !!                    
-      !! ** Purpose :   Compute the before lateral momentum viscous trend 
+      !!
+      !! ** Purpose :   Compute the before lateral momentum viscous trend
       !!              and add it to the general trend of momentum equation.
       !!
