Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ASM/asminc.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ASM/asminc.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ASM/asminc.F90	(revision 12616)
@@ -95,4 +95,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -417,5 +418,6 @@
                      &            - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * u_bkginc(ji-1,jj,jk)    &
                      &            + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * v_bkginc(ji,jj  ,jk)    &
-                     &            - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk)  ) / e3t(ji,jj,jk,Kmm)
+                     &            - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk)  ) &
+                     &            / e3t(ji,jj,jk,Kmm)
                END_2D
                CALL lbc_lnk( 'asminc', zhdiv, 'T', 1. )   ! lateral boundary cond. (no sign change)
@@ -758,5 +760,7 @@
             !
             ssh(:,:,Kbb) = ssh(:,:,Kmm)                        ! Update before fields
+#if ! defined key_LF
             e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
+#endif
 !!gm why not e3u(:,:,:,Kbb), e3v(:,:,:,Kbb), gdept(:,:,:,Kbb) ????
             !
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdydta.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdydta.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdydta.F90	(revision 12616)
@@ -70,4 +70,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -250,5 +251,6 @@
                ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                DO ik = 1, jpkm1
-                  dta_alias%u2d(ib) = dta_alias%u2d(ib) + e3u(ii,ij,ik,Kmm) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik)
+                  dta_alias%u2d(ib) =   &
+                     & dta_alias%u2d(ib) + e3u(ii,ij,ik,Kmm) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik)
                END DO
                dta_alias%u2d(ib) =  dta_alias%u2d(ib) * r1_hu(ii,ij,Kmm)
@@ -263,5 +265,6 @@
                ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                DO ik = 1, jpkm1
-                  dta_alias%v2d(ib) = dta_alias%v2d(ib) + e3v(ii,ij,ik,Kmm) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik)
+                  dta_alias%v2d(ib) =   &
+                     & dta_alias%v2d(ib) + e3v(ii,ij,ik,Kmm) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik)
                END DO
                dta_alias%v2d(ib) =  dta_alias%v2d(ib) * r1_hv(ii,ij,Kmm)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdydyn.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdydyn.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdydyn.F90	(revision 12616)
@@ -35,4 +35,6 @@
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
+#  include "domzgr_substitute.h90"
+
 CONTAINS
 
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/CRS/crsfld.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/CRS/crsfld.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/CRS/crsfld.F90	(revision 12616)
@@ -33,4 +33,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -68,8 +69,10 @@
 
       ! Depth work arrrays
-      ze3t(:,:,:) = e3t(:,:,:,Kmm)
-      ze3u(:,:,:) = e3u(:,:,:,Kmm)
-      ze3v(:,:,:) = e3v(:,:,:,Kmm)
-      ze3w(:,:,:) = e3w(:,:,:,Kmm)
+      DO jk = 1 , jpk 
+         ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
+         ze3u(:,:,jk) = e3u(:,:,jk,Kmm)
+         ze3v(:,:,jk) = e3v(:,:,jk,Kmm)
+         ze3w(:,:,jk) = e3w(:,:,jk,Kmm)
+      END DO
 
       IF( kt == nit000  ) THEN
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/CRS/crsini.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/CRS/crsini.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/CRS/crsini.F90	(revision 12616)
@@ -33,4 +33,6 @@
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
+#  include "domzgr_substitute.h90"
+
 CONTAINS
    
@@ -174,8 +176,17 @@
      
      !
+#if ! defined key_LF
      ze3t(:,:,:) = e3t(:,:,:,Kmm)
      ze3u(:,:,:) = e3u(:,:,:,Kmm)
      ze3v(:,:,:) = e3v(:,:,:,Kmm)
      ze3w(:,:,:) = e3w(:,:,:,Kmm)
+#else
+     DO jk = 1, jpk
+        ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
+        ze3u(:,:,jk) = e3u(:,:,jk,Kmm)
+        ze3v(:,:,jk) = e3v(:,:,jk,Kmm)
+        ze3w(:,:,jk) = e3w(:,:,jk,Kmm)
+     END DO 
+#endif 
 
      !    3.d.2   Surfaces 
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domzgr_substitute.h90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domzgr_substitute.h90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domzgr_substitute.h90	(revision 12616)
@@ -11,5 +11,5 @@
 !! Software governed by the CeCILL license (see ./LICENSE)
 !!----------------------------------------------------------------------
-!#if defined key_lf
+#if defined key_LF
 #   define  e3t(i,j,k,t)   (e3t_0(i,j,k)*(1.+r3t(i,j,t)*tmask(i,j,k)))
 #   define  e3u(i,j,k,t)   (e3u_0(i,j,k)*(1.+r3u(i,j,t)*umask(i,j,k)))
@@ -19,4 +19,4 @@
 #   define  e3uw(i,j,k,t)  (e3uw_0(i,j,k)*(1.+r3u(i,j,t)))
 #   define  e3vw(i,j,k,t)  (e3vw_0(i,j,k)*(1.+r3v(i,j,t)))
-!#endif
+#endif
 !!----------------------------------------------------------------------
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynadv_cen2.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynadv_cen2.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynadv_cen2.F90	(revision 12616)
@@ -28,4 +28,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -79,7 +80,9 @@
          DO_2D_00_00
             puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
-               &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
+               &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj)   &
+               &                           / e3u(ji,jj,jk,Kmm)
             pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    &
-               &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
+               &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj)   &
+               &                           / e3v(ji,jj,jk,Kmm)
          END_2D
       END DO
@@ -115,6 +118,8 @@
       END DO
       DO_3D_00_00( 1, jpkm1 )
-         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
-         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
+         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj)   &
+            &                                      / e3u(ji,jj,jk,Kmm)
+         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj)   &
+            &                                      / e3v(ji,jj,jk,Kmm)
       END_3D
       !
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynadv_ubs.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynadv_ubs.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynadv_ubs.F90	(revision 12616)
@@ -34,4 +34,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -169,7 +170,9 @@
          DO_2D_00_00
             puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
-               &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
+               &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj)   &
+               &                           / e3u(ji,jj,jk,Kmm)
             pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    &
-               &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
+               &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj)   &
+               &                           / e3v(ji,jj,jk,Kmm)
          END_2D
       END DO
@@ -206,6 +209,8 @@
       END DO
       DO_3D_00_00( 1, jpkm1 )
-         puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
-         pvv(ji,jj,jk,Krhs) =  pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
+         puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj)   &
+            &                                       / e3u(ji,jj,jk,Kmm)
+         pvv(ji,jj,jk,Krhs) =  pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj)   &
+            &                                       / e3v(ji,jj,jk,Kmm)
       END_3D
       !
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynhpg.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynhpg.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynhpg.F90	(revision 12616)
@@ -43,5 +43,5 @@
    USE in_out_manager  ! I/O manager
    USE prtctl          ! Print control
-   USE lbclnk          ! lateral boundary condition 
+   USE lbclnk          ! lateral boundary condition
    USE lib_mpp         ! MPP library
    USE eosbn2          ! compute density
@@ -76,4 +76,6 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -204,5 +206,5 @@
       !
       IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
-      ! 
+      !
       IF(lwp) THEN
          WRITE(numout,*)
@@ -217,5 +219,5 @@
          WRITE(numout,*)
       ENDIF
-      !                          
+      !
    END SUBROUTINE dyn_hpg_init
 
@@ -427,5 +429,5 @@
             zcpx(ji,jj) = 0._wp
           END IF
-   
+
           ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                &
                &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            &
@@ -452,8 +454,12 @@
       DO_2D_00_00
          ! hydrostatic pressure gradient along s-surfaces
-         zhpi(ji,jj,1) = zcoef0 * (  e3w(ji+1,jj  ,1,Kmm) * ( znad + rhd(ji+1,jj  ,1) )    &
-            &                      - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e1u(ji,jj)
-         zhpj(ji,jj,1) = zcoef0 * (  e3w(ji  ,jj+1,1,Kmm) * ( znad + rhd(ji  ,jj+1,1) )    &
-            &                      - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e2v(ji,jj)
+         zhpi(ji,jj,1) =   &
+            &  zcoef0 * (  e3w(ji+1,jj  ,1,Kmm) * ( znad + rhd(ji+1,jj  ,1) )    &
+            &            - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) &
+            &           * r1_e1u(ji,jj)
+         zhpj(ji,jj,1) =   &
+            &  zcoef0 * (  e3w(ji  ,jj+1,1,Kmm) * ( znad + rhd(ji  ,jj+1,1) )    &
+            &            - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) &
+            &           * r1_e2v(ji,jj)
          ! s-coordinate pressure gradient correction
          zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
@@ -464,5 +470,5 @@
          IF( ln_wd_il ) THEN
             zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
-            zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
+            zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)
             zuap = zuap * zcpx(ji,jj)
             zvap = zvap * zcpy(ji,jj)
@@ -478,9 +484,9 @@
          ! hydrostatic pressure gradient along s-surfaces
          zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   &
-            &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   &
-            &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
+            &    * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )  &
+            &       - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
          zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)   &
-            &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
-            &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  )
+            &    * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
+            &       - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  )
          ! s-coordinate pressure gradient correction
          zuap = -zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
@@ -491,5 +497,5 @@
          IF( ln_wd_il ) THEN
             zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
-            zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
+            zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)
             zuap = zuap * zcpx(ji,jj)
             zvap = zvap * zcpy(ji,jj)
@@ -522,5 +528,5 @@
       !!         pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj
       !!      iceload is added
-      !!      
+      !!
       !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
       !!----------------------------------------------------------------------
@@ -540,5 +546,5 @@
       znad=1._wp                 ! To use density and not density anomaly
       !
-      !                          ! iniitialised to 0. zhpi zhpi 
+      !                          ! iniitialised to 0. zhpi zhpi
       zhpi(:,:,:) = 0._wp   ;   zhpj(:,:,:) = 0._wp
 
@@ -554,6 +560,6 @@
       CALL eos( zts_top, risfdep, zrhdtop_oce )
 
-!==================================================================================     
-!===== Compute surface value ===================================================== 
+!==================================================================================
+!===== Compute surface value =====================================================
 !==================================================================================
       DO_2D_00_00
@@ -567,10 +573,10 @@
             &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         &
             &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
-            &                                  + ( risfload(ji+1,jj) - risfload(ji,jj))                            ) 
+            &                                  + ( risfload(ji+1,jj) - risfload(ji,jj))                            )
          zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w(ji,jj+1,iktp1j,Kmm)                                    &
             &                                    * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) )   &
-            &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         & 
+            &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         &
             &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
-            &                                  + ( risfload(ji,jj+1) - risfload(ji,jj))                            ) 
+            &                                  + ( risfload(ji,jj+1) - risfload(ji,jj))                            )
          ! s-coordinate pressure gradient correction (=0 if z coordinate)
          zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
@@ -582,6 +588,6 @@
          pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1)
       END_2D
-!==================================================================================     
-!===== Compute interior value ===================================================== 
+!==================================================================================
+!===== Compute interior value =====================================================
 !==================================================================================
       ! interior value (2=<jk=<jpkm1)
@@ -589,9 +595,13 @@
          ! hydrostatic pressure gradient along s-surfaces
          zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
-            &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   &
-            &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   )
+            &           * (  e3w(ji+1,jj,jk,Kmm)                   &
+            &                  * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   &
+            &              - e3w(ji  ,jj,jk,Kmm)                   &
+            &                  * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   )
          zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
-            &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   &
-            &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   )
+            &           * (  e3w(ji,jj+1,jk,Kmm)                   &
+            &                  * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   &
+            &              - e3w(ji,jj  ,jk,Kmm)                   &
+            &                  * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   )
          ! s-coordinate pressure gradient correction
          zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
@@ -650,5 +660,5 @@
             zcpx(ji,jj) = 0._wp
           END IF
-   
+
           ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                &
                &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            &
@@ -771,5 +781,5 @@
       !-------------------------------------------------------------
 
-!!bug gm   :  e3w-gde3w = 0.5*e3w  ....  and gde3w(2)-gde3w(1)=e3w(2) ....   to be verified
+!!bug gm   :  e3w-gde3w = 0.5*e3w  ....  and gde3w(2)-gde3w(1)=e3w(:,:,2,:) ....   to be verified
 !          true if gde3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be
 
@@ -825,5 +835,5 @@
          IF( ln_wd_il ) THEN
            zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
-           zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
+           zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)
          ENDIF
          ! add to the general momentum trend
@@ -845,5 +855,5 @@
          IF( ln_wd_il ) THEN
            zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
-           zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
+           zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)
          ENDIF
          ! add to the general momentum trend
@@ -916,10 +926,10 @@
             zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
                         &    / (ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm)) )
-           
+
              zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
           ELSE
             zcpx(ji,jj) = 0._wp
           END IF
-   
+
           ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                &
                &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            &
@@ -1002,12 +1012,12 @@
 !!gm BUG ?    if it is ssh at u- & v-point then it should be:
 !          zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * &
-!                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
+!                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp
 !          zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji,jj+1) * ssh(ji,jj+1,Kmm)) * &
-!                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
+!                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp
 !!gm not this:
        zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * &
-                      & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
+                      & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp
        zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * &
-                      & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
+                      & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp
       END_2D
 
@@ -1015,5 +1025,5 @@
 
       DO_2D_00_00
-       zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad) 
+       zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad)
        zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad)
       END_2D
@@ -1098,5 +1108,5 @@
             zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj)
          ENDIF
-         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 
+         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)
       ENDIF
 
@@ -1154,6 +1164,6 @@
          ENDIF
          IF( ln_wd_il ) THEN
-            zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 
-            zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 
+            zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj)
+            zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj)
          ENDIF
 
@@ -1189,5 +1199,5 @@
       !!----------------------------------------------------------------------
       !
-!!gm  WHAT !!!!!   THIS IS VERY DANGEROUS !!!!!  
+!!gm  WHAT !!!!!   THIS IS VERY DANGEROUS !!!!!
       jpi   = size(fsp,1)
       jpj   = size(fsp,2)
@@ -1359,3 +1369,2 @@
    !!======================================================================
 END MODULE dynhpg
-
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynspg_ts.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynspg_ts.F90	(revision 12616)
@@ -87,4 +87,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -161,4 +162,5 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points
       REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV         ! fluxes
+      REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v
       !
       REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True. 
@@ -226,6 +228,11 @@
       !                                   !=  zu_frc =  1/H e3*d/dt(Ua)  =!  (Vertical mean of Ua, the 3D trends)
       !                                   !  ---------------------------  !
-      zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm)
-      zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm)
+      DO jk = 1 , jpk
+         ze3u(:,:,jk) = e3u(:,:,jk,Kmm)
+         ze3v(:,:,jk) = e3v(:,:,jk,Kmm)
+      END DO
+      !
+      zu_frc(:,:) = SUM( ze3u(:,:,:) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm)
+      zv_frc(:,:) = SUM( ze3v(:,:,:) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm)
       !
       !
@@ -1087,5 +1094,5 @@
       !
       SELECT CASE( nvor_scheme )
-      CASE( np_EEN )                != EEN scheme using e3f (energy & enstrophy scheme)
+      CASE( np_EEN )                != EEN scheme using e3f energy & enstrophy scheme
          SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point
          CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
@@ -1114,5 +1121,5 @@
          END_2D
          !
-      CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme)
+      CASE( np_EET )                  != EEN scheme using e3t energy conserving scheme
          ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
          DO_2D_01_01
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynzdf.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynzdf.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynzdf.F90	(revision 12616)
@@ -38,4 +38,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -55,8 +56,8 @@
       !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
       !!         u(after) =         u(before) + 2*dt *       u(rhs)                vector form or linear free surf.
-      !!         u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u(after)   otherwise
+      !!         u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u_after   otherwise
       !!               - update the after velocity with the implicit vertical mixing.
       !!      This requires to solver the following system: 
-      !!         u(after) = u(after) + 1/e3u(after) dk+1[ mi(avm) / e3uw(after) dk[ua] ]
+      !!         u(after) = u(after) + 1/e3u_after  dk+1[ mi(avm) / e3uw_after dk[ua] ]
       !!      with the following surface/top/bottom boundary condition:
       !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
@@ -117,8 +118,10 @@
       ELSE                                      ! applied on thickness weighted velocity
          DO jk = 1, jpkm1
-            puu(:,:,jk,Kaa) = (         e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb)  &
-               &          + r2dt * e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs)  ) / e3u(:,:,jk,Kaa) * umask(:,:,jk)
-            pvv(:,:,jk,Kaa) = (         e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb)  &
-               &          + r2dt * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs)  ) / e3v(:,:,jk,Kaa) * vmask(:,:,jk)
+            puu(:,:,jk,Kaa) = (    e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb)       &
+               &          + r2dt * e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs)  )   &
+               &                 / e3u(:,:,jk,Kaa) * umask(:,:,jk)
+            pvv(:,:,jk,Kaa) = (    e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb)       &
+               &          + r2dt * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs)  )   &
+               &                 / e3v(:,:,jk,Kaa) * vmask(:,:,jk)
          END DO
       ENDIF
@@ -136,6 +139,8 @@
             iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
             ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
-            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)
-            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)
+            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    &
+               &             + r_vvl   * e3u(ji,jj,iku,Kaa)
+            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
+               &             + r_vvl   * e3v(ji,jj,ikv,Kaa)
             puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
             pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
@@ -145,6 +150,8 @@
                iku = miku(ji,jj)         ! top ocean level at u- and v-points 
                ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
-               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)
-               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)
+               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    &
+                  &             + r_vvl   * e3u(ji,jj,iku,Kaa)
+               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
+                  &             + r_vvl   * e3v(ji,jj,ikv,Kaa)
                puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
                pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
@@ -161,5 +168,6 @@
          CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
             DO_3D_00_00( 1, jpkm1 )
-               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
+               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &
+                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
                zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
                   &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
@@ -174,7 +182,10 @@
          CASE DEFAULT               ! iso-level lateral mixing
             DO_3D_00_00( 1, jpkm1 )
-               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
-               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
-               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
+               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &    ! after scale factor at U-point
+                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)
+               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )   &
+                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
+               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )   &
+                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
                zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
                zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
@@ -186,6 +197,8 @@
          DO_2D_00_00
             zwi(ji,jj,1) = 0._wp
-            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa)
-            zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
+            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm)    &
+               &             + r_vvl   * e3u(ji,jj,1,Kaa)
+            zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) )   &
+               &         / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
             zWus = ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) / ze3ua
             zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )
@@ -196,5 +209,6 @@
          CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
             DO_3D_00_00( 1, jpkm1 )
-               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
+               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &
+                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
                zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
                   &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
@@ -207,7 +221,10 @@
          CASE DEFAULT               ! iso-level lateral mixing
             DO_3D_00_00( 1, jpkm1 )
-               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
-               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
-               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
+               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &
+                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
+               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )    &
+                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
+               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )    &
+                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
                zwi(ji,jj,jk) = zzwi
                zws(ji,jj,jk) = zzws
@@ -231,5 +248,6 @@
          DO_2D_00_00
             iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
-            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
+            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    &
+               &             + r_vvl   * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
             zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
          END_2D
@@ -238,5 +256,6 @@
                !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
                iku = miku(ji,jj)       ! ocean top level at u- and v-points 
-               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
+               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    &
+                  &             + r_vvl   * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
                zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
             END_2D
@@ -264,5 +283,6 @@
       !
       DO_2D_00_00
-         ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 
+         ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm)    &
+            &             + r_vvl   * e3u(ji,jj,1,Kaa) 
          puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
             &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
@@ -287,5 +307,6 @@
          CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
             DO_3D_00_00( 1, jpkm1 )
-               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
+               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
+                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
                zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
                   &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
@@ -300,7 +321,10 @@
          CASE DEFAULT               ! iso-level lateral mixing
             DO_3D_00_00( 1, jpkm1 )
-               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
-               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
-               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
+               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
+                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
+               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    &
+                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
+               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    &
+                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
                zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
                zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
@@ -312,6 +336,8 @@
          DO_2D_00_00
             zwi(ji,jj,1) = 0._wp
-            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa)
-            zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
+            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm)    &
+               &             + r_vvl   * e3v(ji,jj,1,Kaa)
+            zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) )    &
+               &         / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
             zWvs = ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) / ze3va
             zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
@@ -322,5 +348,6 @@
          CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
             DO_3D_00_00( 1, jpkm1 )
-               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
+               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
+                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
                zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
                   &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
@@ -333,7 +360,10 @@
          CASE DEFAULT               ! iso-level lateral mixing
             DO_3D_00_00( 1, jpkm1 )
-               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
-               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
-               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
+               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
+                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
+               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    &
+                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
+               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    &
+                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
                zwi(ji,jj,jk) = zzwi
                zws(ji,jj,jk) = zzws
@@ -356,5 +386,6 @@
          DO_2D_00_00
             ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
-            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
+            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
+               &             + r_vvl   * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
             zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
          END_2D
@@ -362,5 +393,6 @@
             DO_2D_00_00
                ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
-               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
+               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
+                  &             + r_vvl   * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
                zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va
             END_2D
@@ -388,5 +420,6 @@
       !
       DO_2D_00_00
-         ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 
+         ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm)    &
+            &             + r_vvl   * e3v(ji,jj,1,Kaa) 
          pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
             &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfcpl.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfcpl.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfcpl.F90	(revision 12616)
@@ -15,5 +15,9 @@
    USE isfutils, ONLY : debug
    USE lib_mpp , ONLY: mpp_sum, mpp_max ! mpp routine
+#if ! defined key_LF
    USE domvvl  , ONLY: dom_vvl_zgr      ! vertical scale factor interpolation
+#else
+   USE domqe   , ONLY: dom_qe_zgr      ! vertical scale factor interpolation
+#endif
    USE domngb  , ONLY: dom_ngb          ! find the closest grid point from a given lon/lat position
    !
@@ -43,4 +47,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -112,6 +117,7 @@
       vv   (:,:,:,Kbb)   = vv   (:,:,:,Kmm)
       ssh (:,:,Kbb)     = ssh (:,:,Kmm)
+#if ! defined key_LF
       e3t(:,:,:,Kbb)   = e3t(:,:,:,Kmm)
- 
+#endif 
       ! prepare writing restart
       IF( lwxios ) THEN
@@ -135,11 +141,20 @@
       INTEGER, INTENT(in) :: Kmm    ! ocean time level index
       !!----------------------------------------------------------------------
+      INTEGER :: jk                               ! loop index
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, ze3u, ze3v    ! e3t , e3u, e3v
+      !!----------------------------------------------------------------------
+      !
+      DO jk = 1, jpkm1
+         ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
+         ze3u(:,:,jk) = e3u(:,:,jk,Kmm)
+         ze3v(:,:,jk) = e3v(:,:,jk,Kmm)
+      END DO 
       !
       IF( lwxios ) CALL iom_swap( cwxios_context )
       CALL iom_rstput( kt, nitrst, numrow, 'tmask'  , tmask , ldxios = lwxios )
       CALL iom_rstput( kt, nitrst, numrow, 'ssmask' , ssmask, ldxios = lwxios )
-      CALL iom_rstput( kt, nitrst, numrow, 'e3t_n'  , e3t(:,:,:,Kmm) , ldxios = lwxios )
-      CALL iom_rstput( kt, nitrst, numrow, 'e3u_n'  , e3u(:,:,:,Kmm) , ldxios = lwxios )
-      CALL iom_rstput( kt, nitrst, numrow, 'e3v_n'  , e3v(:,:,:,Kmm) , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'e3t_n'  , ze3t , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'e3u_n'  , ze3u , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'e3v_n'  , ze3v , ldxios = lwxios )
       CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw(:,:,:,Kmm) , ldxios = lwxios )
       IF( lwxios ) CALL iom_swap( cxios_context )
@@ -209,4 +224,5 @@
       IF(lwp) write(numout,*) 'isfcpl_ssh : recompute scale factor from ssh (new wet cell,Kmm)'
       IF(lwp) write(numout,*) '~~~~~~~~~~~'
+#if ! defined key_LF
       DO jk = 1, jpk
          e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
@@ -216,4 +232,7 @@
       e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
       CALL dom_vvl_zgr(Kbb, Kmm, Kaa)
+#else
+      CALL dom_qe_zgr(Kbb, Kmm, Kaa)
+#endif
       !
    END SUBROUTINE isfcpl_ssh
@@ -412,6 +431,8 @@
          ! compute volume flux divergence after coupling
          DO_2D_00_00
-            zqvoln(ji,jj,jk) = (   e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) - e2u(ji-1,jj  ) * e3u(ji-1,jj  ,jk,Kmm) * uu(ji-1,jj  ,jk,Kmm)    &
-               &                 + e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) - e1v(ji  ,jj-1) * e3v(ji  ,jj-1,jk,Kmm) * vv(ji  ,jj-1,jk,Kmm)  ) &
+            zqvoln(ji,jj,jk) = (   e2u(ji  ,jj  ) * e3u(ji  ,jj  ,jk,Kmm) * uu(ji  ,jj  ,jk,Kmm)    &
+               &                 - e2u(ji-1,jj  ) * e3u(ji-1,jj  ,jk,Kmm) * uu(ji-1,jj  ,jk,Kmm)    &
+               &                 + e1v(ji  ,jj  ) * e3v(ji  ,jj  ,jk,Kmm) * vv(ji  ,jj  ,jk,Kmm)    &
+               &                 - e1v(ji  ,jj-1) * e3v(ji  ,jj-1,jk,Kmm) * vv(ji  ,jj-1,jk,Kmm)  ) &
                &               * tmask(ji,jj,jk)
          END_2D
@@ -555,5 +576,7 @@
             DO ji = nldi,nlei
                jip1=MIN(ji+1,jpi) ; jim1=MAX(ji-1,1) ; jjp1=MIN(jj+1,jpj) ; jjm1=MAX(jj-1,1) ;
-               IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) nisfl(narea) = nisfl(narea) + MAX(SUM(tmask(jim1:jip1,jjm1:jjp1,jk)),1._wp)
+               IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) THEN 
+                  nisfl(narea) = nisfl(narea) + MAX(SUM(tmask(jim1:jip1,jjm1:jjp1,jk)),1._wp)
+               ENDIF
             ENDDO
          ENDDO
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfdiags.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfdiags.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfdiags.F90	(revision 12616)
@@ -26,4 +26,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfdynatf.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfdynatf.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfdynatf.F90	(revision 12616)
@@ -14,5 +14,5 @@
 
    USE phycst , ONLY: r1_rau0         ! physical constant
-   USE dom_oce, ONLY: tmask, ssmask, ht, e3t, r1_e1e2t   ! time and space domain
+   USE dom_oce                        ! time and space domain
 
    USE in_out_manager
@@ -25,4 +25,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
 
 CONTAINS
@@ -81,5 +82,6 @@
       ! add the increment
       DO jk = 1, jpkm1
-         pe3t_f(:,:,jk) = pe3t_f(:,:,jk) - tmask(:,:,jk) * zfwfinc(:,:) * e3t(:,:,jk,Kmm)
+         pe3t_f(:,:,jk) = pe3t_f(:,:,jk) - tmask(:,:,jk) * zfwfinc(:,:)   &
+            &                              * e3t(:,:,jk,Kmm)
       END DO
       !
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfhdiv.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfhdiv.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfhdiv.F90	(revision 12616)
@@ -26,4 +26,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
 
 CONTAINS
@@ -134,5 +135,6 @@
       !
       DO jk=1,jpk 
-         phdiv(:,:,jk) =  phdiv(:,:,jk) + pqvol(:,:,jk) * r1_e1e2t(:,:) / e3t(:,:,jk,Kmm)
+         phdiv(:,:,jk) =  phdiv(:,:,jk) + pqvol(:,:,jk) * r1_e1e2t(:,:)   &
+            &                             / e3t(:,:,jk,Kmm)
       END DO
       !
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfstp.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfstp.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfstp.F90	(revision 12616)
@@ -21,5 +21,5 @@
    USE isfcpl , ONLY: isfcpl_rst_write, isfcpl_init ! isf variables
 
-   USE dom_oce, ONLY: ht, e3t, ln_isfcav, ln_linssh     ! ocean space and time domain
+   USE dom_oce                                          ! ocean space and time domain
    USE domvvl,  ONLY: ln_vvl_zstar                      ! zstar logical
    USE zdfdrg,  ONLY: r_Cdmin_top, r_ke0_top            ! vertical physics: top/bottom drag coef.
@@ -41,4 +41,6 @@
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
+#  include "domzgr_substitute.h90"
+
 CONTAINS
  
@@ -60,4 +62,7 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time step
       INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
+      !!----------------------------------------------------------------------
+      INTEGER :: jk                               ! loop index
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t    ! e3t 
       !!---------------------------------------------------------------------
       !
@@ -78,5 +83,8 @@
          ! 1.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
          rhisf_tbl_cav(:,:) = rn_htbl * mskisf_cav(:,:)
-         CALL isf_tbl_lvl(ht, e3t(:,:,:,Kmm), misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav)
+         DO jk = 1, jpkm1
+            ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
+         END DO 
+         CALL isf_tbl_lvl(ht, ze3t, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav)
          !
          ! 1.3: compute ice shelf melt
@@ -100,5 +108,8 @@
          ! by simplicity, we assume the top level where param applied do not change with time (done in init part)
          rhisf_tbl_par(:,:) = rhisf0_tbl_par(:,:)
-         CALL isf_tbl_lvl(ht, e3t(:,:,:,Kmm), misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par)
+         DO jk = 1, jpkm1
+            ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
+         END DO
+         CALL isf_tbl_lvl(ht, ze3t, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par)
          !
          ! 2.3: compute ice shelf melt
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isftbl.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isftbl.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isftbl.F90	(revision 12616)
@@ -25,4 +25,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
 
 CONTAINS
@@ -56,4 +57,6 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zhtbl   ! thickness of the tbl
       REAL(wp), DIMENSION(jpi,jpj) :: zfrac   ! thickness of the tbl
+      INTEGER :: jk                            ! loop index
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t,ze3u,ze3v ! e3 
       !!--------------------------------------------------------------------
       ! 
@@ -64,9 +67,12 @@
          zhtbl = phtbl
          !
+         DO jk = 1, jpkm1
+            ze3u(:,:,jk) = e3u(:,:,jk,Kmm)
+         END DO 
          ! compute tbl lvl and thickness
-         CALL isf_tbl_lvl( hu(:,:,Kmm), e3u(:,:,:,Kmm), ktop, ikbot, zhtbl, zfrac )
+         CALL isf_tbl_lvl( hu(:,:,Kmm), ze3u, ktop, ikbot, zhtbl, zfrac )
          !
          ! compute tbl property at U point
-         CALL isf_tbl_avg( miku, ikbot, zhtbl, zfrac, e3u(:,:,:,Kmm), pvarin, zvarout )
+         CALL isf_tbl_avg( miku, ikbot, zhtbl, zfrac, ze3u, pvarin, zvarout )
          !
          ! compute tbl property at T point
@@ -82,9 +88,12 @@
          zhtbl = phtbl
          !
+         DO jk = 1, jpkm1
+            ze3v(:,:,jk) = e3v(:,:,jk,Kmm)
+         END DO 
          ! compute tbl lvl and thickness
-         CALL isf_tbl_lvl( hv(:,:,Kmm), e3v(:,:,:,Kmm), ktop, ikbot, zhtbl, zfrac )
+         CALL isf_tbl_lvl( hv(:,:,Kmm), ze3v, ktop, ikbot, zhtbl, zfrac )
          !
          ! compute tbl property at V point
-         CALL isf_tbl_avg( mikv, ikbot, zhtbl, zfrac, e3v(:,:,:,Kmm), pvarin, zvarout )
+         CALL isf_tbl_avg( mikv, ikbot, zhtbl, zfrac, ze3v, pvarin, zvarout )
          !
          ! pvarout is an averaging of wet point
@@ -98,5 +107,8 @@
          !
          ! compute tbl property at T point
-         CALL isf_tbl_avg( ktop, kbot, phtbl, pfrac, e3t(:,:,:,Kmm), pvarin, pvarout )
+         DO jk = 1, jpkm1
+            ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
+         END DO 
+         CALL isf_tbl_avg( ktop, kbot, phtbl, pfrac, ze3t, pvarin, pvarout )
          !
       END SELECT
@@ -212,5 +224,5 @@
       ! phtbl need to be bounded by water column thickness before
       ! test: if htbl = water column thickness, should return mbathy
-      ! test: if htbl = 0 should return ktop (phtbl cap to e3t(ji,jj,1))
+      ! test: if htbl = 0 should return ktop (phtbl cap to pe3t(ji,jj,1))
       !
       ! get ktbl
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcice_cice.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcice_cice.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcice_cice.F90	(revision 12616)
@@ -233,39 +233,42 @@
 !!gm This should be put elsewhere....   (same remark for limsbc)
 !!gm especially here it is assumed zstar coordinate, but it can be ztilde....
-            IF( .NOT.ln_linssh )   CALL dom_vvl_zgr( Kbb, Kmm, Kaa )   ! interpolation scale factor, depth and water column
-            ! IF( .NOT.ln_linssh ) THEN
-            !    !
-            !    DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
-            !       e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-            !       e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-            !    ENDDO
-            !    e3t(:,:,:,Krhs) = e3t(:,:,:,Kbb)
-            !    ! Reconstruction of all vertical scale factors at now and before time-steps
-            !    ! =============================================================================
-            !    ! Horizontal scale factor interpolations
-            !    ! --------------------------------------
-            !    CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )
-            !    CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )
-            !    CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
-            !    CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
-            !    CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )
-            !    ! Vertical scale factor interpolations
-            !    ! ------------------------------------
-            !    CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )
-            !    CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
-            !    CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
-            !    CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
-            !    CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
-            !    ! t- and w- points depth
-            !    ! ----------------------
-            !    gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
-            !    gdepw(:,:,1,Kmm) = 0.0_wp
-            !    gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
-            !    DO jk = 2, jpk
-            !       gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk,Kmm)
-            !       gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm)
-            !       gde3w(:,:,jk)     = gdept(:,:,jk  ,Kmm) - sshn   (:,:)
-            !    END DO
-            ! ENDIF
+#if defined key_LF
+            IF( .NOT.ln_linssh )   CALL dom_qe_zgr( Kbb, Kmm, Kaa )   ! interpolation scale factor, depth and water column
+#else
+            IF( .NOT.ln_linssh ) THEN
+               !
+               DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
+                  e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
+                  e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
+               ENDDO
+               e3t(:,:,:,Krhs) = e3t(:,:,:,Kbb)
+               ! Reconstruction of all vertical scale factors at now and before time-steps
+               ! =============================================================================
+               ! Horizontal scale factor interpolations
+               ! --------------------------------------
+               CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )
+               CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )
+               CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
+               CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
+               CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )
+               ! Vertical scale factor interpolations
+               ! ------------------------------------
+               CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )
+               CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
+               CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
+               CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
+               CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
+               ! t- and w- points depth
+               ! ----------------------
+               gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
+               gdepw(:,:,1,Kmm) = 0.0_wp
+               gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
+               DO jk = 2, jpk
+                  gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk,Kmm)
+                  gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm)
+                  gde3w(:,:,jk)     = gdept(:,:,jk  ,Kmm) - sshn   (:,:)
+               END DO
+            ENDIF
+#endif
          ENDIF
       ENDIF
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdglo.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdglo.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdglo.F90	(revision 12616)
@@ -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 / rau0           ! 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_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdken.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdken.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdken.F90	(revision 12616)
@@ -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_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdmxl.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdmxl.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdmxl.F90	(revision 12616)
@@ -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_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdpen.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdpen.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdpen.F90	(revision 12616)
@@ -40,4 +40,6 @@
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
+#  include "domzgr_substitute.h90"
+
 CONTAINS
 
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdtra.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdtra.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdtra.F90	(revision 12616)
@@ -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_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfdrg.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfdrg.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfdrg.F90	(revision 12616)
@@ -74,4 +74,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfgls.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfgls.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfgls.F90	(revision 12616)
@@ -263,7 +263,9 @@
          zcof = rfact_tke * tmask(ji,jj,jk)
          !                                        ! lower diagonal, in fact not used for jk = 2 (see surface conditions)
-         zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) )
+         zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) )   &
+            &                 / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) )
          !                                        ! upper diagonal, in fact not used for jk = ibotm1 (see bottom conditions)
-         zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk,Kmm) )
+         zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) )   &
+            &                 / ( e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk,Kmm) )
          !                                        ! diagonal
          zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk)  + rdt * zdiss * wmask(ji,jj,jk) 
@@ -473,7 +475,9 @@
          zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
          !                                               ! lower diagonal
-         zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) )
+         zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) )   &
+            &            / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) )
          !                                               ! upper diagonal
-         zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk,Kmm) )
+         zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) )   &
+            &            / ( e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk,Kmm) )
          !                                               ! diagonal
          zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk)
@@ -1100,3 +1104,2 @@
    !!======================================================================
 END MODULE zdfgls
-
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfosm.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfosm.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfosm.F90	(revision 12616)
@@ -105,4 +105,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -503,5 +504,6 @@
                        & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) + zvel_max
 
-                  zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w(ji,jj,jk,Kmm) )
+                  zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ),   &
+                     &                     e3w(ji,jj,jk,Kmm) )
                   zhbl_s = MIN(zhbl_s, ht(ji,jj))
 
@@ -594,5 +596,7 @@
                      zwb_ent(ji,jj) = 0._wp
                   ENDIF
-                  inhml = MAX( INT( zari * zhbl(ji,jj) / e3t(ji,jj,ibld(ji,jj),Kmm) ) , 1 )
+                  inhml =   &
+                     & MAX( INT( zari * zhbl(ji,jj) / e3t(ji,jj,ibld(ji,jj),Kmm) ),  &
+                     &      1 )
                   imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1)
                   zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
@@ -610,5 +614,7 @@
                      zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
                        & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 )
-                     inhml = MAX( INT( zari * zhbl(ji,jj) / e3t(ji,jj,ibld(ji,jj),Kmm) ) , 1 )
+                     inhml =   &
+                        & MAX( INT( zari * zhbl(ji,jj) / e3t(ji,jj,ibld(ji,jj),Kmm) ),   &
+                        &      1 )
                      imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1)
                      zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdftke.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdftke.F90	(revision 12615)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdftke.F90	(revision 12616)
@@ -90,4 +90,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -467,6 +468,8 @@
             &            gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) )
             ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
-            zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
-            zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
+            zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk)   &
+               &            + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
+            zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk)   &
+               &            + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
          END_3D
          !
