Index: NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynzdf.F90
===================================================================
--- NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynzdf.F90	(revision 9839)
+++ NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynzdf.F90	(revision 9863)
@@ -11,4 +11,6 @@
    !!----------------------------------------------------------------------
    !!   dyn_zdf       : compute the after velocity through implicit calculation of vertical mixing
+   !!       zdf_trd   : diagnose the zdf velocity trends and the KE dissipation trend 
+!!gm                        ==>>> zdf_trd currently not used
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers variables
@@ -26,4 +28,5 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
+   USE iom             ! IOM library
    USE prtctl         ! Print control
    USE timing         ! Timing
@@ -67,10 +70,10 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
-      INTEGER  ::   ji, jj, jk         ! dummy loop indices
-      INTEGER  ::   iku, ikv           ! local integers
-      REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars
-      REAL(wp) ::   zzws, ze3va        !   -      -
-      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::  zwi, zwd, zws   ! 3D workspace 
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      -
+      INTEGER  ::   ji, jj, jk            ! dummy loop indices
+      INTEGER  ::   iku, ikv              ! local integers
+      REAL(wp) ::   zzwi, ze3ua, z2dt_2   ! local scalars
+      REAL(wp) ::   zzws, ze3va           !   -      -
+      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwd, zws   ! 3D workspace 
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv    !  -      -
       !!---------------------------------------------------------------------
       !
@@ -86,8 +89,7 @@
          ENDIF
       ENDIF
-      !                             !* set time step
-      IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdt (restart with Euler time stepping)
-      ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdt (leapfrog)
-      ENDIF
+      !
+      z2dt_2 = r2dt * 0.5_wp        !* =rdt except in 1st Euler time step where it is equal to rdt/2
+      !
       !
       !                             !* explicit top/bottom drag case
@@ -111,8 +113,6 @@
       ELSE                                      ! applied on thickness weighted velocity
          DO jk = 1, jpkm1
-            ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  &
-               &          + r2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk)
-            va(:,:,jk) = (         e3v_b(:,:,jk) * vb(:,:,jk)  &
-               &          + r2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk)
+            ua(:,:,jk) = ( e3u_b(:,:,jk)*ub(:,:,jk) + r2dt * e3u_n(:,:,jk)*ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk)
+            va(:,:,jk) = ( e3v_b(:,:,jk)*vb(:,:,jk) + r2dt * e3v_n(:,:,jk)*va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk)
          END DO
       ENDIF
@@ -133,6 +133,6 @@
                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
-               ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua
-               va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va
+               ua(ji,jj,iku) = ua(ji,jj,iku) + z2dt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua
+               va(ji,jj,ikv) = va(ji,jj,ikv) + z2dt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va
             END DO
          END DO
@@ -144,6 +144,6 @@
                   ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
-                  ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua
-                  va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va
+                  ua(ji,jj,iku) = ua(ji,jj,iku) + z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua
+                  va(ji,jj,ikv) = va(ji,jj,ikv) + z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va
                END DO
             END DO
@@ -153,16 +153,15 @@
       !              !==  Vertical diffusion on u  ==!
       !
-      !                    !* Matrix construction
-      zdt = r2dt * 0.5
-      SELECT CASE( nldf_dyn )
-      CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
+      SELECT CASE( nldf_dyn )    !* Matrix construction
+      !
+      CASE( np_lap_i )              ! rotated lateral mixing: add its vertical mixing (akzu)
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1 
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point
-                  zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
-                     &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
-                  zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
-                     &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                  zzwi = - r2dt * ( 0.5 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) + akzu(ji,jj,jk  ) )   &
+                     &            / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                  zzws = - r2dt * ( 0.5 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) + akzu(ji,jj,jk+1) )   &
+                     &            / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
                   zwi(ji,jj,jk) = zzwi
                   zws(ji,jj,jk) = zzws
@@ -171,11 +170,11 @@
             END DO
          END DO
-      CASE DEFAULT               ! iso-level lateral mixing
+      CASE DEFAULT                  ! iso-level lateral mixing
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1 
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point
-                  zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
-                  zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                  zzwi = - z2dt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                  zzws = - z2dt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
                   zwi(ji,jj,jk) = zzwi
                   zws(ji,jj,jk) = zzws
@@ -186,6 +185,6 @@
       END SELECT
       !
-      DO jj = 2, jpjm1     !* Surface boundary conditions
-         DO ji = fs_2, fs_jpim1   ! vector opt.
+      DO jj = 2, jpjm1           !* Surface boundary conditions
+         DO ji = fs_2, fs_jpim1     ! vector opt.
             zwi(ji,jj,1) = 0._wp
             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
@@ -204,5 +203,5 @@
                iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! 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
+               zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
             END DO
          END DO
@@ -213,5 +212,5 @@
                   iku = miku(ji,jj)       ! ocean top level at u- and v-points 
                   ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! 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
+                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
                END DO
             END DO
@@ -245,6 +244,5 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
-            ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
-               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
+            ua(ji,jj,1) = ua(ji,jj,1) + z2dt_2 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( ze3ua * rau0 ) * umask(ji,jj,1)
          END DO
       END DO
@@ -272,16 +270,15 @@
       !              !==  Vertical diffusion on v  ==!
       !
-      !                       !* Matrix construction
-      zdt = r2dt * 0.5
-      SELECT CASE( nldf_dyn )
-      CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
+      !                      
+      SELECT CASE( nldf_dyn )    !* Matrix construction
+      CASE( np_lap_i )              ! rotated lateral mixing: add its vertical mixing (akzu)
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1   
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
-                  zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
-                     &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
-                  zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
-                     &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                  zzwi = - r2dt * ( 0.5 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) + akzv(ji,jj,jk  ) )   &
+                     &            / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
+                  zzws = - r2dt * ( 0.5 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) + akzv(ji,jj,jk+1) )   &
+                     &            / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
                   zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
                   zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
@@ -290,11 +287,11 @@
             END DO
          END DO
-      CASE DEFAULT               ! iso-level lateral mixing
+      CASE DEFAULT                  ! iso-level lateral mixing
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1   
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
-                  zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
-                  zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                  zzwi = - z2dt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
+                  zzws = - z2dt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
                   zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
                   zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
@@ -305,6 +302,6 @@
       END SELECT
       !
-      DO jj = 2, jpjm1        !* Surface boundary conditions
-         DO ji = fs_2, fs_jpim1   ! vector opt.
+      DO jj = 2, jpjm1           !* Surface boundary conditions
+         DO ji = fs_2, fs_jpim1     ! vector opt.
             zwi(ji,jj,1) = 0._wp
             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
@@ -322,5 +319,5 @@
                ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! 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           
+               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - z2dt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
             END DO
          END DO
@@ -330,5 +327,5 @@
                   ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! 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) ) / ze3va
+                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va
                END DO
             END DO
@@ -362,6 +359,5 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.          
             ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
-            va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
-               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
+            va(ji,jj,1) = va(ji,jj,1) + z2dt_2 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( ze3va * rau0 ) * vmask(ji,jj,1)
          END DO
       END DO
@@ -387,7 +383,12 @@
       END DO
       !
-      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
-         ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:)
-         ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:)
+      IF( l_trddyn ) THEN                        ! save the vertical diffusive trends for further diagnostics
+         IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
+            ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_2dt - ztrdu(:,:,:)
+            ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_2dt - ztrdv(:,:,:)
+         ELSE                                      ! applied on thickness weighted velocity
+            ztrdu(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_2dt - ztrdu(:,:,:)
+            ztrdv(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_2dt - ztrdv(:,:,:)
+         ENDIF
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
          DEALLOCATE( ztrdu, ztrdv ) 
@@ -401,4 +402,106 @@
    END SUBROUTINE dyn_zdf
 
+!!gm currently not used : just for memory to be able to add dissipation trend through vertical mixing
+
+   SUBROUTINE zdf_trd( ptrdu, ptrdv ,kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_trd  ***
+      !!
+      !! ** Purpose :   compute the trend due to the vert. momentum diffusion
+      !!              together with the Leap-Frog time stepping using an 
+      !!              implicit scheme.
+      !!
+      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
+      !!         ua =         ub + 2*dt *       ua             vector form or linear free surf.
+      !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise
+      !!               - update the after velocity with the implicit vertical mixing.
+      !!      This requires to solver the following system: 
+      !!         ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ]
+      !!      with the following surface/top/bottom boundary condition:
+      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
+      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90)
+      !!
+      !! ** Action :   (ua,va)   after velocity 
+      !!---------------------------------------------------------------------
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   ptrdu, ptrdv   ! 3D work arrays use for zdf trends diag
+      INTEGER , INTENT(in   )                         ::   kt             ! ocean time-step index
+      !
+      INTEGER  ::   ji, jj, jk       ! dummy loop indices
+      REAL(wp) ::   zzz              ! local scalar
+      REAL(wp) ::   zavmu, zavmum1   !   -      -
+      REAL(wp) ::   zavmv, zavmvm1   !   -      -
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   z2d    !  -      -
+      !!---------------------------------------------------------------------
+      !
+      CALL lbc_lnk_multi( ua, 'U', -1., va, 'V', -1. )   ! apply lateral boundary condition on (ua,va)
+      !
+      !
+      !                 !==  momentum trend due to vertical diffusion  ==!
+      !
+      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
+         ptrdu(:,:,:) = (              ua(:,:,:) -              ub(:,:,:) )                * r1_2dt - ptrdu(:,:,:)
+         ptrdv(:,:,:) = (              va(:,:,:) -              vb(:,:,:) )                * r1_2dt - ptrdv(:,:,:)
+      ELSE                                      ! applied on thickness weighted velocity
+         ptrdu(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_2dt - ptrdu(:,:,:)
+         ptrdv(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_2dt - ptrdv(:,:,:)
+      ENDIF
+      CALL trd_dyn( ptrdu, ptrdv, jpdyn_zdf, kt )
+      !
+      !
+      !                 !==  KE dissipation trend due to vertical diffusion  ==!
+      !
+      IF( iom_use( 'dispkevfo' ) ) THEN   ! ocean kinetic energy dissipation per unit area
+         !                                ! due to v friction (v=vertical) 
+         !                                ! see NEMO_book appendix C, §C.8 (N.B. here averaged at t-points)
+         !                                ! Note that formally, in a Leap-Frog environment, the shear**2 should be the product of 
+         !                                ! now by before shears, i.e. the source term of TKE (local positivity is not ensured).
+         !                                ! Note also that now e3 scale factors are used as after one are not computed !
+         !
+         CALL wrk_alloc(jpi,jpj,   z2d )
+         z2d(:,:) = 0._wp
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  zavmu   = 0.5 * ( avm(ji+1,jj,jk) + avm(ji  ,jj,jk) )
+                  zavmum1 = 0.5 * ( avm(ji  ,jj,jk) + avm(ji-1,jj,jk) )
+                  zavmv   = 0.5 * ( avm(ji,jj+1,jk) + avm(ji,jj  ,jk) )
+                  zavmvm1 = 0.5 * ( avm(ji,jj  ,jk) + avm(ji,jj-1,jk) )
+               
+                  z2d(ji,jj) = z2d(ji,jj)  +  (                                                                                  &
+                     &   zavmu   * ( ua(ji  ,jj,jk-1) - ua(ji  ,jj,jk) )**2 / e3uw_n(ji  ,jj,jk) * wumask(ji  ,jj,jk)   &
+                     & + zavmum1 * ( ua(ji-1,jj,jk-1) - ua(ji-1,jj,jk) )**2 / e3uw_n(ji-1,jj,jk) * wumask(ji-1,jj,jk)   &
+                     & + zavmv   * ( va(ji,jj  ,jk-1) - va(ji,jj  ,jk) )**2 / e3vw_n(ji,jj  ,jk) * wvmask(ji,jj  ,jk)   &
+                     & + zavmvm1 * ( va(ji,jj-1,jk-1) - va(ji,jj-1,jk) )**2 / e3vw_n(ji,jj-1,jk) * wvmask(ji,jj-1,jk)   &
+                     &                        )
+!!gm --- This trends is in fact properly computed in zdf_sh2 but with a backward shift of one time-step  ===>>> use it ?
+!!                                                                                     No since in zdfshé only kz tke (or gls) is used
+!!
+!!gm --- formally, as done below, in a Leap-Frog environment, the shear**2 should be the product of
+!!gm     now by before shears, i.e. the source term of TKE (local positivity is not ensured).
+!!       CAUTION: requires to compute e3uw_a and e3vw_a !!!
+!                  z2d(ji,jj) = z2d(ji,jj)  + (                                                                                 &
+!                     &    avmu(ji  ,jj,jk) * ( un(ji  ,jj,jk-1) - un(ji  ,jj,jk) ) / e3uw_n(ji  ,jj,jk)                        &
+!                     &                     * ( ua(ji  ,jj,jk-1) - ua(ji  ,jj,jk) ) / e3uw_a(ji  ,jj,jk) * wumask(ji  ,jj,jk)   &
+!                     &  + avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) / e3uw_n(ji-1,jj,jk)                        &
+!                     &                       ( ua(ji-1,jj,jk-1) - ua(ji-1,jj,jk) ) / e3uw_a(ji-1,jj,jk) * wumask(ji-1,jj,jk)   &
+!                     &  + avmv(ji,jj  ,jk) * ( vn(ji,jj  ,jk-1) - vn(ji,jj  ,jk) ) / e3vw_n(ji,jj  ,jk)                        &
+!                     &                       ( va(ji,jj  ,jk-1) - va(ji,jj  ,jk) ) / e3vw_a(ji,jj  ,jk) * wvmask(ji,jj  ,jk)   &
+!                     &  + avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) / e3vw_n(ji,jj-1,jk)                        &
+!                     &                       ( va(ji,jj-1,jk-1) - va(ji,jj-1,jk) ) / e3vw_a(ji,jj-1,jk) * wvmask(ji,jj-1,jk)   &
+!                     &                       )
+!!gm end
+               END DO
+            END DO
+         END DO
+         zzz= - 0.5_wp* rau0           ! caution sign minus here
+         z2d(:,:) = zzz * z2d(:,:) 
+         CALL lbc_lnk( z2d,'T', 1. )
+         CALL iom_put( 'dispkevfo', z2d )
+      ENDIF
+      !
+   END SUBROUTINE zdf_trd
+   
+!!gm end not used
+
    !!==============================================================================
 END MODULE dynzdf
