Index: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/DYN/dynzdf.F90
===================================================================
--- NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/DYN/dynzdf.F90	(revision 12939)
+++ NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/DYN/dynzdf.F90	(revision 13248)
@@ -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)
@@ -112,8 +113,10 @@
       ELSE                                      ! applied on thickness weighted velocity
          DO_3D_00_00( 1, jpkm1 )
-            puu(ji,jj,jk,Kaa) = (         e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb)  &
-               &                  + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs)  ) / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk)
-            pvv(ji,jj,jk,Kaa) = (         e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb)  &
-               &                  + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs)  ) / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk)
+            puu(ji,jj,jk,Kaa) = (         e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb )  &
+               &                  + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs)  ) &
+               &                        / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk)
+            pvv(ji,jj,jk,Kaa) = (         e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb )  &
+               &                  + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs)  ) &
+               &                        / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk)
          END_3D
       ENDIF
@@ -131,6 +134,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) + rDt * 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) + rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
@@ -140,6 +145,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) + rDt * 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) + rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
@@ -156,5 +163,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  )
@@ -169,7 +177,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
@@ -181,6 +192,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 )
@@ -191,5 +204,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  )
@@ -202,7 +216,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
@@ -226,5 +243,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) - rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
          END_2D
@@ -233,5 +251,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) - rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
             END_2D
@@ -259,5 +278,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) + rDt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
             &                                      / ( ze3ua * rho0 ) * umask(ji,jj,1) 
@@ -282,5 +302,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  )
@@ -295,7 +316,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
@@ -307,6 +331,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 )
@@ -317,5 +343,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  )
@@ -328,7 +355,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
@@ -351,5 +381,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) - rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
          END_2D
@@ -357,5 +388,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) - rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va
             END_2D
@@ -383,5 +415,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) + rDt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
             &                                      / ( ze3va * rho0 ) * vmask(ji,jj,1) 
