Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90	(revision 7990)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90	(revision 7991)
@@ -9,5 +9,5 @@
    USE dom_oce         ! ocean space and time domain
    USE zdf_oce         ! ocean vertical physics    
-   USE zdfgls   , ONLY : hmxn
+   USE zdfgls   , ONLY : hmxl_n
    USE in_out_manager  ! I/O units
    USE iom             ! I/0 library
@@ -23,5 +23,6 @@
 
    ! variables for calculating 25-hourly means
-   INTEGER , SAVE ::   cnt_25h     ! Counter for 25 hour means
+   INTEGER , SAVE ::   cnt_25h           ! Counter for 25 hour means
+   REAL(wp), SAVE ::   r1_25 = 0.04_wp   ! =1/25 
    REAL(wp), SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_25h  , sn_25h
    REAL(wp), SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_25h 
@@ -94,18 +95,18 @@
       ! ------------------------- !
       cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible) 
-      tn_25h(:,:,:) = tsb(:,:,:,jp_tem)
-      sn_25h(:,:,:) = tsb(:,:,:,jp_sal)
-      sshn_25h(:,:) = sshb(:,:)
-      un_25h(:,:,:) = ub(:,:,:)
-      vn_25h(:,:,:) = vb(:,:,:)
-      wn_25h(:,:,:) = wn(:,:,:)
-      avt_25h(:,:,:) = avt(:,:,:)
-      avm_25h(:,:,:) = avm(:,:,:)
+      tn_25h  (:,:,:) = tsb (:,:,:,jp_tem)
+      sn_25h  (:,:,:) = tsb (:,:,:,jp_sal)
+      sshn_25h(:,:)   = sshb(:,:)
+      un_25h  (:,:,:) = ub  (:,:,:)
+      vn_25h  (:,:,:) = vb  (:,:,:)
+      wn_25h  (:,:,:) = wn  (:,:,:)
+      avt_25h (:,:,:) = avt (:,:,:)
+      avm_25h (:,:,:) = avm (:,:,:)
       IF( ln_zdftke ) THEN
          en_25h(:,:,:) = en(:,:,:)
       ENDIF
       IF( ln_zdfgls ) THEN
-         en_25h(:,:,:) = en(:,:,:)
-         rmxln_25h(:,:,:) = hmxn(:,:,:)
+         en_25h   (:,:,:) = en    (:,:,:)
+         rmxln_25h(:,:,:) = hmxl_n(:,:,:)
       ENDIF
 #if defined key_lim3 || defined key_lim2
@@ -156,18 +157,18 @@
          ENDIF
 
-         tn_25h(:,:,:)        = tn_25h(:,:,:) + tsn(:,:,:,jp_tem)
-         sn_25h(:,:,:)        = sn_25h(:,:,:) + tsn(:,:,:,jp_sal)
-         sshn_25h(:,:)        = sshn_25h(:,:) + sshn (:,:)
-         un_25h(:,:,:)        = un_25h(:,:,:) + un(:,:,:)
-         vn_25h(:,:,:)        = vn_25h(:,:,:) + vn(:,:,:)
-         wn_25h(:,:,:)        = wn_25h(:,:,:) + wn(:,:,:)
-         avt_25h(:,:,:)       = avt_25h(:,:,:) + avt(:,:,:)
-         avm_25h(:,:,:)       = avm_25h(:,:,:) + avm(:,:,:)
-         IF( ln_zdftke ) THEN
-            en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:)
-         ENDIF
-         IF( ln_zdfgls ) THEN
-            en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:)
-            rmxln_25h(:,:,:)      = rmxln_25h(:,:,:) + hmxn(:,:,:)
+         tn_25h  (:,:,:)     = tn_25h  (:,:,:) + tsn (:,:,:,jp_tem)
+         sn_25h  (:,:,:)     = sn_25h  (:,:,:) + tsn (:,:,:,jp_sal)
+         sshn_25h(:,:)       = sshn_25h(:,:)   + sshn(:,:)
+         un_25h  (:,:,:)     = un_25h  (:,:,:) + un  (:,:,:)
+         vn_25h  (:,:,:)     = vn_25h  (:,:,:) + vn  (:,:,:)
+         wn_25h  (:,:,:)     = wn_25h  (:,:,:) + wn  (:,:,:)
+         avt_25h (:,:,:)     = avt_25h (:,:,:) + avt (:,:,:)
+         avm_25h (:,:,:)     = avm_25h (:,:,:) + avm (:,:,:)
+         IF( ln_zdftke ) THEN
+            en_25h(:,:,:)    = en_25h  (:,:,:) + en(:,:,:)
+         ENDIF
+         IF( ln_zdfgls ) THEN
+            en_25h   (:,:,:) = en_25h   (:,:,:) + en    (:,:,:)
+            rmxln_25h(:,:,:) = rmxln_25h(:,:,:) + hmxl_n(:,:,:)
          ENDIF
          cnt_25h = cnt_25h + 1
@@ -187,18 +188,18 @@
          ENDIF
          !
-         tn_25h(:,:,:)        = tn_25h(:,:,:) / 25.0_wp
-         sn_25h(:,:,:)        = sn_25h(:,:,:) / 25.0_wp
-         sshn_25h(:,:)        = sshn_25h(:,:) / 25.0_wp
-         un_25h(:,:,:)        = un_25h(:,:,:) / 25.0_wp
-         vn_25h(:,:,:)        = vn_25h(:,:,:) / 25.0_wp
-         wn_25h(:,:,:)        = wn_25h(:,:,:) / 25.0_wp
-         avt_25h(:,:,:)       = avt_25h(:,:,:) / 25.0_wp
-         avm_25h(:,:,:)       = avm_25h(:,:,:) / 25.0_wp
-         IF( ln_zdftke ) THEN
-            en_25h(:,:,:)        = en_25h(:,:,:) / 25.0_wp
-         ENDIF
-         IF( ln_zdfgls ) THEN
-            en_25h(:,:,:)        = en_25h(:,:,:) / 25.0_wp
-            rmxln_25h(:,:,:)       = rmxln_25h(:,:,:) / 25.0_wp
+         tn_25h  (:,:,:) = tn_25h  (:,:,:) * r1_25
+         sn_25h  (:,:,:) = sn_25h  (:,:,:) * r1_25
+         sshn_25h(:,:)   = sshn_25h(:,:)   * r1_25
+         un_25h  (:,:,:) = un_25h  (:,:,:) * r1_25
+         vn_25h  (:,:,:) = vn_25h  (:,:,:) * r1_25
+         wn_25h  (:,:,:) = wn_25h  (:,:,:) * r1_25
+         avt_25h (:,:,:) = avt_25h (:,:,:) * r1_25
+         avm_25h (:,:,:) = avm_25h (:,:,:) * r1_25
+         IF( ln_zdftke ) THEN
+            en_25h(:,:,:) = en_25h(:,:,:) * r1_25
+         ENDIF
+         IF( ln_zdfgls ) THEN
+            en_25h   (:,:,:) = en_25h   (:,:,:) * r1_25
+            rmxln_25h(:,:,:) = rmxln_25h(:,:,:) * r1_25
          ENDIF
          !
@@ -236,18 +237,18 @@
          !
          ! After the write reset the values to cnt=1 and sum values equal current value 
-         tn_25h(:,:,:) = tsn(:,:,:,jp_tem)
-         sn_25h(:,:,:) = tsn(:,:,:,jp_sal)
-         sshn_25h(:,:) = sshn (:,:)
-         un_25h(:,:,:) = un(:,:,:)
-         vn_25h(:,:,:) = vn(:,:,:)
-         wn_25h(:,:,:) = wn(:,:,:)
-         avt_25h(:,:,:) = avt(:,:,:)
-         avm_25h(:,:,:) = avm(:,:,:)
+         tn_25h  (:,:,:) = tsn (:,:,:,jp_tem)
+         sn_25h  (:,:,:) = tsn (:,:,:,jp_sal)
+         sshn_25h(:,:)   = sshn(:,:)
+         un_25h  (:,:,:) = un  (:,:,:)
+         vn_25h  (:,:,:) = vn  (:,:,:)
+         wn_25h  (:,:,:) = wn  (:,:,:)
+         avt_25h (:,:,:) = avt (:,:,:)
+         avm_25h (:,:,:) = avm (:,:,:)
          IF( ln_zdftke ) THEN
             en_25h(:,:,:) = en(:,:,:)
          ENDIF
          IF( ln_zdfgls ) THEN
-            en_25h(:,:,:) = en(:,:,:)
-            rmxln_25h(:,:,:) = hmxn(:,:,:)
+            en_25h   (:,:,:) = en    (:,:,:)
+            rmxln_25h(:,:,:) = hmxl_n(:,:,:)
          ENDIF
          cnt_25h = 1
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 7990)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 7991)
@@ -19,4 +19,5 @@
    USE domvvl         ! ocean space and time domain : variable volume layer
    USE zdf_oce        ! ocean vertical physics
+   USE zdfsh2         ! vertical physics: shear production term of TKE
    USE zdfbfr         ! bottom friction (only for rn_bfrz0)
    USE sbc_oce        ! surface boundary condition: ocean
@@ -42,5 +43,5 @@
 
    !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hmxn    !: now mixing length
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hmxl_n    !: now mixing length
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zwall   !: wall function
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustars2 !: Squared surface velocity scale at T-points
@@ -114,6 +115,6 @@
       !!                ***  FUNCTION zdf_gls_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( hmxn(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     &
-         &      ustars2(jpi,jpj) , ustarb2(jpi,jpj)   , STAT= zdf_gls_alloc )
+      ALLOCATE( hmxl_n(jpi,jpj,jpk) , ustars2(jpi,jpj) ,       &
+         &      zwall (jpi,jpj,jpk) , ustarb2(jpi,jpj) ,   STAT= zdf_gls_alloc )
          !
       IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc )
@@ -141,5 +142,5 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhsro       ! Surface roughness (surface waves)
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   eb          ! tke at time before
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   mxlb        ! mixing length at time before
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   hmxl_b      ! mixing length at time before
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   shear       ! vertical shear
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   eps         ! dissipation rate
@@ -149,5 +150,5 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_b    ! element of the second matrix diagonal
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_c    ! element of the third  matrix diagonal
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3du, z3dv  ! u- and v-shears
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zstt, zstm  ! stability function on tracer and momentum
       !!--------------------------------------------------------------------
       !
@@ -155,6 +156,6 @@
       !
       CALL wrk_alloc( jpi,jpj,       zdep, zkar, zflxs, zhsro )
-      CALL wrk_alloc( jpi,jpj,jpk,   eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )
-      CALL wrk_alloc( jpi,jpj,jpk,   z3du, z3dv )
+      CALL wrk_alloc( jpi,jpj,jpk,   eb, hmxl_b, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )
+      CALL wrk_alloc( jpi,jpj,jpk,   zstt, zstm )
 
       ! Preliminary computing
@@ -198,15 +199,12 @@
       END SELECT
 
-      DO jk = 2, jpkm1              !==  Compute shear and dissipation rate  ==!
+      !                             !==  Compute shear and dissipation rate  ==!
+      CALL zdf_sh2( shear )
+
+
+      DO jk = 2, jpkm1              !==  Compute dissipation rate  ==!
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
-               z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji+1,jj,jk) )   &
-                  &                 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )   &
-                  &                 * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
-               z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji,jj+1,jk) )   &
-                  &                 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )   &
-                  &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
-                  !
-               eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / hmxn(ji,jj,jk)
+               eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk)
             END DO
          END DO
@@ -214,6 +212,6 @@
 
       ! Save tke at before time step
-      eb  (:,:,:) = en  (:,:,:)
-      mxlb(:,:,:) = hmxn(:,:,:)
+      eb    (:,:,:) = en    (:,:,:)
+      hmxl_b(:,:,:) = hmxl_n(:,:,:)
 
       IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
@@ -221,5 +219,5 @@
             DO jj = 2, jpjm1 
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zup   = hmxn(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1)
+                  zup   = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1)
                   zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) )
                   zcoef = ( zup / MAX( zdown, rsmall ) )
@@ -247,13 +245,7 @@
             DO ji = 2, jpim1
                !
-               ! shear prod. at w-point weightened by mask
-               shear(ji,jj,jk) =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
-                  &             + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
-               !
-               ! stratif. destruction
-               buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk)
-               !
-               ! shear prod. - stratif. destruction
-               diss = eps(ji,jj,jk)
+               buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk)       ! stratif. destruction
+               !
+               diss = eps(ji,jj,jk)                         ! dissipation
                !
                dir = 0.5_wp + SIGN( 0.5_wp, shear(ji,jj,jk) + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
@@ -344,5 +336,5 @@
       !
       CASE ( 0 )             ! Dirichlet 
-         !                      ! en(ibot) = u*^2 / Co2 and hmxn(ibot) = rn_lmin
+         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin
          !                      ! Balance between the production and the dissipation terms
          DO jj = 2, jpjm1
@@ -425,5 +417,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  psi(ji,jj,jk)  = eb(ji,jj,jk) * mxlb(ji,jj,jk)
+                  psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk)
                END DO
             END DO
@@ -443,5 +435,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * mxlb(ji,jj,jk) )
+                  psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) )
                END DO
             END DO
@@ -452,5 +444,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn 
+                  psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 
                END DO
             END DO
@@ -563,5 +555,5 @@
       !
       CASE ( 0 )             ! Dirichlet 
-         !                      ! en(ibot) = u*^2 / Co2 and hmxn(ibot) = vkarmn * rn_bfrz0
+         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * rn_bfrz0
          !                      ! Balance between the production and the dissipation terms
          DO jj = 2, jpjm1
@@ -686,13 +678,13 @@
       ! Limit dissipation rate under stable stratification
       ! --------------------------------------------------
-      DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxn at the same time
+      DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxl_n at the same time
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
                ! limitation
-               eps(ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
-               hmxn(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
+               eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
+               hmxl_n(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
                ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 
                zrn2 = MAX( rn2(ji,jj,jk), rsmall )
-               IF (ln_length_lim) hmxn(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxn(ji,jj,jk) )
+               IF( ln_length_lim )   hmxl_n(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) )
             END DO
          END DO
@@ -710,5 +702,5 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! zcof =  l²/q²
-                  zcof = mxlb(ji,jj,jk) * mxlb(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
+                  zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
                   ! Gh = -N²l²/q²
                   gh = - rn2(ji,jj,jk) * zcof
@@ -719,7 +711,7 @@
                   sm = ( rb1**(-1._wp/3._wp) + ( 18._wp*ra1*ra1 + 9._wp*ra1*ra2*(1._wp-rc2) )*sh*gh ) / (1._wp-9._wp*ra1*ra2*gh)
                   !
-                  ! Store stability function in z3du and z3dv
-                  z3du(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
-                  z3dv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
+                  ! Store stability function in zstt and zstm
+                  zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
+                  zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
                END DO
             END DO
@@ -731,5 +723,5 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! zcof =  l²/q²
-                  zcof = mxlb(ji,jj,jk)*mxlb(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
+                  zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
                   ! Gh = -N²l²/q²
                   gh = - rn2(ji,jj,jk) * zcof
@@ -747,7 +739,7 @@
                   sh = (rs4 - rs5*gh + rs6*gm) / rcff
                   !
-                  ! Store stability function in z3du and z3dv
-                  z3du(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
-                  z3dv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
+                  ! Store stability function in zstt and zstm
+                  zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
+                  zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
                END DO
             END DO
@@ -759,9 +751,10 @@
       ! Lines below are useless if GOTM style Dirichlet conditions are used
 
-      z3dv(:,:,1) = z3dv(:,:,2)
-
+      zstm(:,:,1) = zstm(:,:,2)
+
+!!gm should be done for ISF (top boundary cond.)
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            z3dv(ji,jj,mbkt(ji,jj)+1) = z3dv(ji,jj,mbkt(ji,jj))
+            zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj))
          END DO
       END DO
@@ -772,8 +765,8 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zsqen         = SQRT( 2._wp * en(ji,jj,jk) ) * hmxn(ji,jj,jk)
-               zav           = zsqen * z3du(ji,jj,jk)
+               zsqen         = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk)
+               zav           = zsqen * zstt(ji,jj,jk)
                avt(ji,jj,jk) = MAX( zav, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine
-               zav           = zsqen * z3dv(ji,jj,jk)
+               zav           = zsqen * zstm(ji,jj,jk)
                avm(ji,jj,jk) = MAX( zav, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom
             END DO
@@ -795,6 +788,6 @@
       !
       CALL wrk_dealloc( jpi,jpj,       zdep, zkar, zflxs, zhsro )
-      CALL wrk_dealloc( jpi,jpj,jpk,   eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )
-      CALL wrk_dealloc( jpi,jpj,jpk,   z3du, z3dv )
+      CALL wrk_dealloc( jpi,jpj,jpk,   eb, hmxl_b, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )
+      CALL wrk_dealloc( jpi,jpj,jpk,   zstt, zstm )
       !
       IF( nn_timing == 1 )   CALL timing_stop('zdf_gls')
@@ -1100,5 +1093,5 @@
 
       !                                !* read or initialize all required files  
-      CALL gls_rst( nit000, 'READ' )      ! (en, avt_k, avm_k, hmxn)
+      CALL gls_rst( nit000, 'READ' )      ! (en, avt_k, avm_k, hmxl_n)
       !
       IF( nn_timing == 1 )  CALL timing_stop('zdf_gls_init')
@@ -1129,26 +1122,26 @@
          !                                   ! ---------------
          IF( ln_rstart ) THEN                   !* Read the restart file
-            id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
-            id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
-            id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
-            id4 = iom_varid( numror, 'hmxn' , ldstop = .FALSE. )
+            id1 = iom_varid( numror, 'en'    , ldstop = .FALSE. )
+            id2 = iom_varid( numror, 'avt_k' , ldstop = .FALSE. )
+            id3 = iom_varid( numror, 'avm_k' , ldstop = .FALSE. )
+            id4 = iom_varid( numror, 'hmxl_n', ldstop = .FALSE. )
             !
             IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN        ! all required arrays exist
-               CALL iom_get( numror, jpdom_autoglo, 'en'   , en    )
-               CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k )
-               CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k )
-               CALL iom_get( numror, jpdom_autoglo, 'hmxn' , hmxn  )
+               CALL iom_get( numror, jpdom_autoglo, 'en'    , en     )
+               CALL iom_get( numror, jpdom_autoglo, 'avt_k' , avt_k  )
+               CALL iom_get( numror, jpdom_autoglo, 'avm_k' , avm_k  )
+               CALL iom_get( numror, jpdom_autoglo, 'hmxl_n', hmxl_n )
             ELSE                        
-               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and hmxn computed by iterative loop'
-               en   (:,:,:) = rn_emin
-               hmxn (:,:,:) = 0.05_wp
-               avt_k(:,:,:) = avt(:,:,:)
-               avm_k(:,:,:) = avm(:,:,:)
+               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and hmxl_n computed by iterative loop'
+               en    (:,:,:) = rn_emin
+               hmxl_n(:,:,:) = 0.05_wp
+               avt_k (:,:,:) = avt(:,:,:)
+               avm_k (:,:,:) = avm(:,:,:)
                DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_gls( jit )   ;   END DO
             ENDIF
          ELSE                                   !* Start from rest
-            IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and hmxn by background values'
-            en  (:,:,:) = rn_emin
-            hmxn(:,:,:) = 0.05_wp
+            IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and hmxl_n by background values'
+            en    (:,:,:) = rn_emin
+            hmxl_n(:,:,:) = 0.05_wp
             DO jk = 1, jpk
                avt_k(:,:,jk) = avtb(jk) * wmask(:,:,jk)
@@ -1160,8 +1153,8 @@
          !                                   ! -------------------
          IF(lwp) WRITE(numout,*) '---- gls-rst ----'
-         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    ) 
-         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k )
-         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k )
-         CALL iom_rstput( kt, nitrst, numrow, 'hmxn' , hmxn  )
+         CALL iom_rstput( kt, nitrst, numrow, 'en'    , en     ) 
+         CALL iom_rstput( kt, nitrst, numrow, 'avt_k' , avt_k  )
+         CALL iom_rstput( kt, nitrst, numrow, 'avm_k' , avm_k  )
+         CALL iom_rstput( kt, nitrst, numrow, 'hmxl_n', hmxl_n )
          !
       ENDIF
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 7990)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 7991)
@@ -21,5 +21,6 @@
    USE oce            ! ocean dynamics and tracers variables
    USE dom_oce        ! ocean space and time domain variables
-   USE zdf_oce        ! ocean vertical physics
+   USE zdf_oce        ! vertical physics: variables
+   USE zdfsh2         ! vertical physics: shear production term of TKE
    USE phycst         ! physical constants
    USE sbc_oce,  ONLY :   taum
@@ -104,64 +105,32 @@
       INTEGER  ::   ji, jj, jk                 ! dummy loop indices
       LOGICAL  ::   ll_av_weight = .TRUE.      ! local logical
-      REAL(wp) ::   zsh2, zcfRi, zav, zustar   ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   zdu2, zdv2, zh_ekm   ! 2D workspace
+      REAL(wp) ::   zcfRi, zav, zustar   ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zh_ekm   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsh2     ! 3D     -
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )   CALL timing_start('zdf_ric')
       !
-      DO jk = 2, jpkm1        !==  avm and avt = F(Richardson number)  ==!
-         !
-      IF( ll_av_weight ) THEN    !== viscosity weighted shear & stratification terms
-         !
-         DO jj = 1, jpjm1           !* viscosity weighted shear term
+      !                       !==  avm and avt = F(Richardson number)  ==!
+      !
+      !                          !* Shear production at uw- and vw-points (energy conserving form)
+      CALL zdf_sh2( zsh2 )
+      !
+      DO jk = 2, jpkm1           !* Vertical eddy viscosity and diffusivity coefficients
+         DO jj = 1, jpjm1
             DO ji = 1, jpim1
-               zdu2(ji,jj) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji+1,jj,jk) ) * wumask(ji,jj,jk)      &
-                  &              * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )                         &
-                  &              * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
-               zdv2(ji,jj) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji,jj+1,jk) ) * wumask(ji,jj,jk)      &
-                  &              * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )                         &
-                  &              * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
-            END DO
-         END DO
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1        !* Richardson number dependent coefficient
-               !                          ! shear of horizontal velocity
-               zsh2  =  ( zdu2(ji-1,jj) + zdu2(ji,jj) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
-                  &   + ( zdv2(ji,jj-1) + zdv2(ji,jj) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
-               !                          ! coefficient = F(richardson number)
-               zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , avt(ji,jj,jk) * rn2(ji,jj,jk) / ( zsh2 + 1.e-20 ) )  )
+               !                          ! coefficient = F(richardson number) (avm-weighted Ri)
+               zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( zsh2(ji,jj,jk) + 1.e-20 ) )  )
+               zav   = rn_avmri * zcfRi**nn_ric
                !
-               !                    !* Vertical eddy viscosity and diffusivity coefficients
-               zav = rn_avmri * zcfRi**nn_ric
+               !                          ! avm and avt coefficients
                avm(ji,jj,jk) = MAX(  zav         , avmb(jk)  ) * wmask(ji,jj,jk)
                avt(ji,jj,jk) = MAX(  zav * zcfRi , avtb(jk)  ) * wmask(ji,jj,jk)
             END DO
          END DO
-      ELSE                    !==  classical Richardson number definition  ==!
-         DO jj = 1, jpjm1           !* shear term
-            DO ji = 1, jpim1
-               zdu2(ji,jj) = 0.5 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )                         &
-                  &              * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
-               zdv2(ji,jj) = 0.5 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )                         &
-                  &              * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
-            END DO
-         END DO
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1        !* Richardson number dependent coefficient
-               !                          ! shear of horizontal velocity
-               zsh2  =  ( zdu2(ji-1,jj) + zdu2(ji,jj) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
-                  &   + ( zdv2(ji,jj-1) + zdv2(ji,jj) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
-               !                          ! coefficient = F(richardson number)
-               zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , rn2(ji,jj,jk) / ( zsh2 + 1.e-20 ) )  )
-               !
-               !                    !* Vertical eddy viscosity and diffusivity coefficients
-               zav = rn_avmri * zcfRi**nn_ric
-               avm(ji,jj,jk) = MAX(  zav         , avmb(jk)  ) * wmask(ji,jj,jk)
-               avt(ji,jj,jk) = MAX(  zav * zcfRi , avtb(jk)  ) * wmask(ji,jj,jk)
-            END DO
-         END DO
-      ENDIF
-         !
       END DO
+      !
+!!gm BUG <<<<====  This param can't work at low latitude 
+!!gm               it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! )
       !
       IF( ln_mldw ) THEN      !==  set a minimum value in the Ekman layer  ==!
@@ -184,5 +153,5 @@
             END DO
          END DO
-      END IF
+      ENDIF
       !
       IF( nn_timing == 1 )   CALL timing_stop('zdf_ric')
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfsh2.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfsh2.F90	(revision 7991)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfsh2.F90	(revision 7991)
@@ -0,0 +1,91 @@
+MODULE zdfsh2
+   !!======================================================================
+   !!                       ***  MODULE  zdfsh2  ***
+   !! Ocean physics:  shear production term of TKE 
+   !!=====================================================================
+   !! History :   -   !  2014-10  (A. Barthelemy, G. Madec)  original code
+   !!   NEMO     4.0  !  2017-04  (G. Madec)  remove u-,v-pts avm
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   zdf_sh2       : compute mixing the shear production term of TKE
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean: shared variables
+   USE dom_oce        ! domain: ocean
+   USE zdf_oce        ! vertical physics: variables
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE timing         ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   zdf_sh2        ! called by zdftke, zdfglf, and zdfric
+   
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
+   !! $Id: $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE zdf_sh2(  psh2  ) 
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE zdf_sh2  ***
+      !!
+      !! ** Purpose :   Compute the shear production term of a TKE equation
+      !!
+      !! ** Method  : - a stable discretization of this term is linked to the
+      !!                time-space discretization of the vertical diffusion
+      !!                of the OGCM. NEMO uses C-grid, a leap-frog environment 
+      !!                and an implicit computation of vertical mixing term,
+      !!                so the shear production at w-point is given by:
+      !!                   sh2 = mi[   mi(avm) * dk[ub]/e3ub * dk[un]/e3un   ] 
+      !!                       + mj[   mj(avm) * dk[vb]/e3vb * dk[vn]/e3vn   ] 
+      !!                NB: wet-point only horizontal averaging of shear
+      !!
+      !! ** Action  : - psh2 shear prod. term at w-point (interior ocean domain only)
+      !!
+      !! References :   Bruchard, OM 2002
+      !! ---------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   psh2   ! shear production of TKE (w-points)
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop arguments
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsh2u, zsh2v   ! 2D workspace
+      !!--------------------------------------------------------------------
+      IF( nn_timing == 1 )  CALL timing_start('zdf_sh2')
+      !
+      DO jk = 2, jpkm1
+         DO jj = 1, jpjm1        !* Shear production at uw- and vw-points (energy conserving form)
+            DO ji = 1, jpim1
+               zsh2u(ji,jj) = 0.5 * ( avm(ji+1,jj,jk) + avm(ji,jj,jk) ) * wumask(ji,jj,jk)  &
+                  &               * (  un(ji,jj,jk-1) -  un(ji,jj,jk) )   &
+                  &               * (  ub(ji,jj,jk-1) -  ub(ji,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
+               zsh2v(ji,jj) = 0.5 * ( avm(ji,jj+1,jk) + avm(ji,jj,jk) ) * wvmask(ji,jj,jk)  &
+                  &               * (  vn(ji,jj,jk-1) -  vn(ji,jj,jk) )   &
+                  &               * (  vb(ji,jj,jk-1) -  vb(ji,jj,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
+            END DO
+         END DO
+         DO jj = 2, jpjm1        !* shear production at w-point
+            DO ji = 2, jpim1
+               !
+!!gm original type of coding
+!              psh2(ji,jj,jk) = ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) / MAX( 1._wp , umask(ji-1,jj,jk)+umask(ji,jj,jk) )   &
+!                 &           + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) / MAX( 1._wp , vmask(ji,jj-1,jk)+vmask(ji,jj,jk) )    
+!
+!!gm optimized but it changes the last digits after 1 year of GYRE
+               !                       ! coast mask: =1 at the coast ; =2 otherwise
+               !                       ! no need of wmask as zsh2 are already masked
+               psh2(ji,jj,jk) = ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2._wp - umask(ji-1,jj,jk) * umask(ji,jj,jk) )   &
+                  &           + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2._wp - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )
+            END DO
+         END DO
+      END DO 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('zdf_sh2')     
+      !
+   END SUBROUTINE zdf_sh2
+
+   !!======================================================================
+END MODULE zdfsh2
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 7990)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 7991)
@@ -43,4 +43,5 @@
    USE sbc_oce        ! surface boundary condition: ocean
    USE zdf_oce        ! vertical physics: ocean variables
+   USE zdfsh2         ! vertical physics: shear production term of TKE
    USE zdfmxl         ! vertical physics: mixed layer
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
@@ -89,9 +90,4 @@
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apdlr          ! now mixing lenght of dissipation
-#if defined key_c1d
-   !                                                                        !!** 1D cfg only  **   ('key_c1d')
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_dis, e_mix   !: dissipation and mixing turbulent lengh scales
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_pdl, e_ric   !: prandl and local Richardson numbers
-#endif
 
    !! * Substitutions
@@ -108,11 +104,5 @@
       !!                ***  FUNCTION zdf_tke_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE(                                                                    &
-#if defined key_c1d
-         &      e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) ,                          &
-         &      e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) ,                          &
-#endif
-         &      htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) ,     & 
-         &      apdlr(jpi,jpj,jpk) ,                                           STAT= zdf_tke_alloc      )
+      ALLOCATE( htau(jpi,jpj) , dissl(jpi,jpj,jpk) , apdlr(jpi,jpj,jpk) ,   STAT= zdf_tke_alloc )
          !
       IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc )
@@ -211,19 +201,19 @@
       !! ---------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
-!!bfr      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! temporary scalar
-!!bfr      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          ! temporary scalar
-      REAL(wp) ::   zrhoa  = 1.22                   ! Air density kg/m3
-      REAL(wp) ::   zcdrag = 1.5e-3                 ! drag coefficient
-      REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars
-      REAL(wp) ::   zfact1, zfact2, zfact3          !    -         -
-      REAL(wp) ::   ztx2  , zty2  , zcof            !    -         -
-      REAL(wp) ::   ztau  , zdif                    !    -         -
-      REAL(wp) ::   zus   , zwlc  , zind            !    -         -
-      REAL(wp) ::   zzd_up, zzd_lw                  !    -         -
-!!bfr      REAL(wp) ::   zebot                           !    -         -
+!!bfr      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! local integers
+!!bfr      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          !   -      -
+!!bfr      REAL(wp) ::   zebot                           ! local scalars
+      REAL(wp) ::   zrhoa  = 1.22            ! Air density kg/m3
+      REAL(wp) ::   zcdrag = 1.5e-3          ! drag coefficient
+      REAL(wp) ::   zbbrau, zri              ! local scalars
+      REAL(wp) ::   zfact1, zfact2, zfact3   !   -         -
+      REAL(wp) ::   ztx2  , zty2  , zcof     !   -         -
+      REAL(wp) ::   ztau  , zdif             !   -         -
+      REAL(wp) ::   zus   , zwlc  , zind     !   -         -
+      REAL(wp) ::   zzd_up, zzd_lw           !   -         -
       INTEGER , DIMENSION(jpi,jpj) ::   imlc
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhlc
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv
-      REAL(wp)                            ::   zri  !   local Richardson number
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zpelc, zdiag, zd_up, zd_lw
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsh2
       !!--------------------------------------------------------------------
       !
@@ -231,5 +221,5 @@
       !
       CALL wrk_alloc( jpi,jpj,       zhlc ) 
-      CALL wrk_alloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv ) 
+      CALL wrk_alloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw ) 
       !
       zbbrau = rn_ebb / rau0       ! Local constant initialisation
@@ -328,36 +318,18 @@
       !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
       !
-      DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji+1,jj,jk) )   &
-                  &                 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )   &
-                  &                 * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) * wumask(ji,jj,jk) &
-                  &                 / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
-               z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji,jj+1,jk) )   &
-                  &                 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )   &
-                  &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) * wvmask(ji,jj,jk) &
-                  &                 / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
-            END DO
-         END DO
-      END DO
-      !
-      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: compute apdlr
-         ! Note that zesh2 is also computed in the next loop.
-         ! We decided to compute it twice to keep code readability and avoid an IF case in the DO loops
+      !                           !* Shear production at uw- and vw-points (energy conserving form)
+      CALL zdf_sh2( zsh2 )
+      !
+      IF( nn_pdl == 1 ) THEN      !* Prandtl number = F( Ri )
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  !                                          ! shear prod. at w-point weightened by mask
-                  zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
-                     &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
-                  !                                          ! local Richardson number
-                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zesh2 + rn_bshear )
+               DO ji = 2, jpim1
+                  !                             ! local Richardson number
+                  zri = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zsh2(ji,jj,jk) + rn_bshear )
+                  !                             ! inverse of Prandtl number
                   apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  )
-                  
-               END DO
-            END DO
-         END DO
-         !
+               END DO
+            END DO
+         END DO
       ENDIF
       !         
@@ -373,16 +345,13 @@
                   &          /    ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  )  )
                !
-               !                                   ! shear prod. at w-point weightened by mask
-               zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
-                  &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
-               !
                zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
                zd_lw(ji,jj,jk) = zzd_lw
-               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
+               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk)
                !
                !                                   ! right hand side in en
-               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    &
-                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) &
-                  &                                 * wmask(ji,jj,jk)
+               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zsh2(ji,jj,jk)                           &   ! shear
+                  &                                 - avt(ji,jj,jk) * rn2(ji,jj,jk)            &   ! stratification
+                  &                                 + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation
+                  &                                ) * wmask(ji,jj,jk)
             END DO
          END DO
@@ -470,5 +439,5 @@
       !
       CALL wrk_dealloc( jpi,jpj,       zhlc ) 
-      CALL wrk_dealloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv ) 
+      CALL wrk_dealloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('tke_tke')
@@ -631,8 +600,4 @@
       END SELECT
       !
-# if defined key_c1d
-      e_dis(:,:,:) = zmxld(:,:,:)      ! c1d configuration : save mixing and dissipation turbulent length scales
-      e_mix(:,:,:) = zmxlm(:,:,:)
-# endif
 
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -657,7 +622,4 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
-# if defined key_c1d
-                  e_pdl(ji,jj,jk) = apdlr(ji,jj,jk) * wmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number
-# endif
               END DO
             END DO
