Index: branches/2015/dev_agrif_v3_6_STABLE/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
===================================================================
--- branches/2015/dev_agrif_v3_6_STABLE/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 6099)
+++ branches/2015/dev_agrif_v3_6_STABLE/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 6102)
@@ -316,5 +316,5 @@
    !
 # if defined key_zdftke
-   CALL Agrif_Update_tke(0)
+!   CALL Agrif_Update_tke(0)
 # endif
    !
Index: branches/2015/dev_agrif_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90
===================================================================
--- branches/2015/dev_agrif_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 6099)
+++ branches/2015/dev_agrif_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 6102)
@@ -43,8 +43,5 @@
    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avmu , avmv    !: vertical viscosity coef at uw- & vw-pts       [m2/s]
    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avm  , avt     !: vertical viscosity & diffusivity coef at w-pt [m2/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_k , avm_k  ! not enhanced Kz
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avmu_k, avmv_k ! not enhanced Kz
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy   [m2/s2]
-
+ 
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011)
@@ -62,9 +59,6 @@
          &     avtb(jpk) , bfrva(jpi,jpj) , avtb_2d(jpi,jpj) ,      &
          &     tfrua(jpi, jpj), tfrva(jpi, jpj)              ,      &
-         &     avmu  (jpi,jpj,jpk), avm   (jpi,jpj,jpk),            &
-         &     avmv  (jpi,jpj,jpk), avt   (jpi,jpj,jpk),            &
-         &     avt_k (jpi,jpj,jpk), avm_k (jpi,jpj,jpk),            & 
-         &     avmu_k(jpi,jpj,jpk), avmv_k(jpi,jpj,jpk),            & 
-         &     en    (jpi,jpj,jpk), STAT = zdf_oce_alloc )
+         &     avmu(jpi,jpj,jpk), avm(jpi,jpj,jpk)           ,      &
+         &     avmv(jpi,jpj,jpk), avt(jpi,jpj,jpk)           , STAT = zdf_oce_alloc )
          !
       IF( zdf_oce_alloc /= 0 )   CALL ctl_warn('zdf_oce_alloc: failed to allocate arrays')
Index: branches/2015/dev_agrif_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- branches/2015/dev_agrif_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 6099)
+++ branches/2015/dev_agrif_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 6102)
@@ -42,6 +42,11 @@
    LOGICAL , PUBLIC, PARAMETER ::   lk_zdfgls = .TRUE.   !: TKE vertical mixing flag
    !
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en      !: now turbulent kinetic energy
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mxln    !: now mixing length
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zwall   !: wall function
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_k   ! not enhanced Kz
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avm_k   ! not enhanced Kz
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avmu_k  ! not enhanced Kz
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avmv_k  ! not enhanced Kz
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustars2 !: Squared surface velocity scale at T-points
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustarb2 !: Squared bottom  velocity scale at T-points
@@ -115,5 +120,7 @@
       !!                ***  FUNCTION zdf_gls_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     &
+      ALLOCATE( en(jpi,jpj,jpk),  mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     &
+         &      avt_k (jpi,jpj,jpk) , avm_k (jpi,jpj,jpk),                    &
+         &      avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk),                    &
          &      ustars2(jpi,jpj), ustarb2(jpi,jpj)                      , STAT= zdf_gls_alloc )
          !
@@ -322,5 +329,6 @@
       ! 
       ! One level below
-      en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)
+      en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) &
+          &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)
       en(:,:,2) = MAX(en(:,:,2), rn_emin )
       z_elem_a(:,:,2) = 0._wp 
@@ -343,5 +351,6 @@
       z_elem_a(:,:,2) = 0._wp
       zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) ))
-      zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf)
+      zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) &
+           &                      * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf)
 
       en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2)
Index: branches/2015/dev_agrif_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- branches/2015/dev_agrif_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 6099)
+++ branches/2015/dev_agrif_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 6102)
@@ -53,8 +53,4 @@
    USE timing         ! Timing
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-#if defined key_agrif
-   USE agrif_opa_interp
-   USE agrif_opa_update
-#endif
 
    IMPLICIT NONE
@@ -89,7 +85,9 @@
    REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3)
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy   [m2/s2]
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau)
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation
-   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apdlr          ! now mixing lenght of dissipation
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_k , avm_k  ! not enhanced Kz
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avmu_k, avmv_k ! not enhanced Kz
 #if defined key_c1d
    !                                                                        !!** 1D cfg only  **   ('key_c1d')
@@ -102,5 +100,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.6 , NEMO Consortium (2011)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -117,6 +115,7 @@
          &      e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) ,                          &
 #endif
-         &      apdlr(jpi,jpj,jpk) , htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) ,     & 
-         &      STAT= zdf_tke_alloc      )
+         &      en    (jpi,jpj,jpk) , htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) ,     & 
+         &      avt_k (jpi,jpj,jpk) , avm_k (jpi,jpj,jpk),                          &
+         &      avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), STAT= zdf_tke_alloc      )
          !
       IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc )
@@ -174,9 +173,4 @@
       !!----------------------------------------------------------------------
       !
-#if defined key_agrif 
-      ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2)
-      IF( .NOT.Agrif_Root() )   CALL Agrif_Tke
-#endif
-      !
       IF( kt /= nit000 ) THEN   ! restore before value to compute tke
          avt (:,:,:) = avt_k (:,:,:) 
@@ -195,10 +189,5 @@
       avmv_k(:,:,:) = avmv(:,:,:) 
       !
-#if defined key_agrif
-      ! Update child grid f => parent grid 
-      IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only
-#endif      
-     ! 
-  END SUBROUTINE zdf_tke
+   END SUBROUTINE zdf_tke
 
 
@@ -234,6 +223,5 @@
       INTEGER , POINTER, DIMENSION(:,:  ) :: 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
       !!--------------------------------------------------------------------
       !
@@ -242,5 +230,5 @@
       CALL wrk_alloc( jpi,jpj, imlc )    ! integer
       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
@@ -344,37 +332,18 @@
       !
       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) &
-                  &                 / (  fse3uw_n(ji,jj,jk) * fse3uw_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) &
-                  &                 / (  fse3vw_n(ji,jj,jk) * fse3vw_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
-         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 )
-                  apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  )
-                  
-               END DO
-            END DO
-         END DO
-         !
-      ENDIF
-         !         
+         DO jj = 1, jpj                 ! here avmu, avmv used as workspace
+            DO ji = 1, jpi
+               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
+                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   & 
+                  &                            / (  fse3uw_n(ji,jj,jk)               &
+                  &                              *  fse3uw_b(ji,jj,jk)  )
+               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   &
+                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   &
+                  &                            / (  fse3vw_n(ji,jj,jk)               &
+                  &                              *  fse3vw_b(ji,jj,jk)  )
+            END DO
+         END DO
+      END DO
+      !
       DO jk = 2, jpkm1           !* Matrix and right hand side in en
          DO jj = 2, jpjm1
@@ -385,8 +354,8 @@
                zzd_lw = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &  ! lower diagonal
                   &          / ( fse3t(ji,jj,jk-1) * fse3w(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) )    
-               !
+                  !                                                           ! shear prod. at w-point weightened by mask
+               zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
+                  &    + ( avmv(ji,jj-1,jk) + avmv(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
@@ -486,5 +455,5 @@
       CALL wrk_dealloc( jpi,jpj, imlc )    ! integer
       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')
@@ -530,5 +499,5 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zrn2, zraug, zcoef, zav     ! local scalars
-      REAL(wp) ::   zdku, zri, zsqen     !   -      -
+      REAL(wp) ::   zdku, zpdlr, zri, zsqen     !   -      -
       REAL(wp) ::   zdkv, zemxl, zemlm, zemlp   !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld
@@ -691,8 +660,19 @@
             DO jj = 2, jpjm1
                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)
+                  zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk)
+                  !                                          ! shear
+                  zdku = avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) * ( ub(ji-1,jj,jk-1) - ub(ji-1,jj,jk) )   &
+                    &  + avmu(ji  ,jj,jk) * ( un(ji  ,jj,jk-1) - un(ji  ,jj,jk) ) * ( ub(ji  ,jj,jk-1) - ub(ji  ,jj,jk) )
+                  zdkv = avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) * ( vb(ji,jj-1,jk-1) - vb(ji,jj-1,jk) )   &
+                    &  + avmv(ji,jj  ,jk) * ( vn(ji,jj  ,jk-1) - vn(ji,jj  ,jk) ) * ( vb(ji,jj  ,jk-1) - vb(ji,jj  ,jk) )
+                  !                                          ! local Richardson number
+                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear )
+                  zpdlr = MAX(  0.1_wp,  0.2 / MAX( 0.2 , zri )  )
+!!gm and even better with the use of the "true" ri_crit=0.22222...  (this change the results!)
+!!gm              zpdlr = MAX(  0.1_wp,  ri_crit / MAX( ri_crit , zri )  )
+                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
 # if defined key_c1d
-                  e_pdl(ji,jj,jk) = apdlr(ji,jj,jk) * tmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number
-                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri
+                  e_pdl(ji,jj,jk) = zpdlr * wmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number
+                  e_ric(ji,jj,jk) = zri   * wmask(ji,jj,jk)  ! c1d config. : save Ri
 # endif
               END DO
