Index: trunk/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 1531)
+++ trunk/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 1532)
@@ -1,93 +1,105 @@
-MODULE zdftke_old
+MODULE zdftke
    !!======================================================================
-   !!                       ***  MODULE  zdftke_old  ***
+   !!                       ***  MODULE  zdftke  ***
    !! Ocean physics:  vertical mixing coefficient computed from the tke 
    !!                 turbulent closure parameterization
    !!=====================================================================
-   !! History :   OPA  !  1991-03  (b. blanke)  Original code
-   !!             7.0  !  1991-11  (G. Madec)   bug fix
-   !!             7.1  !  1992-10  (G. Madec)   new mixing length and eav
-   !!             7.2  !  1993-03  (M. Guyon)   symetrical conditions
-   !!             7.3  !  1994-08  (G. Madec, M. Imbard)  nn_pdl flag
-   !!             7.5  !  1996-01  (G. Madec)   s-coordinates
-   !!             8.0  !  1997-07  (G. Madec)   lbc
-   !!             8.1  !  1999-01  (E. Stretta) new option for the mixing length
-   !!  NEMO       1.0  !  2002-06  (G. Madec) add zdf_tke_init routine
-   !!              -   !  2002-08  (G. Madec)  rn_cri and Free form, F90
-   !!              -   !  2004-10  (C. Ethe )  1D configuration
-   !!             2.0  !  2006-07  (S. Masson)  distributed restart using iom
-   !!             3.0  !  2008-05  (C. Ethe,  G.Madec) : update TKE physics:
-   !!                              - tke penetration (wind steering)
-   !!                              - suface condition for tke & mixing length
-   !!                              - Langmuir cells
-   !!              -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
-   !!              -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
+   !! History :  OPA  !  1991-03  (b. blanke)  Original code
+   !!            7.0  !  1991-11  (G. Madec)   bug fix
+   !!            7.1  !  1992-10  (G. Madec)   new mixing length and eav
+   !!            7.2  !  1993-03  (M. Guyon)   symetrical conditions
+   !!            7.3  !  1994-08  (G. Madec, M. Imbard)  nn_pdl flag
+   !!            7.5  !  1996-01  (G. Madec)   s-coordinates
+   !!            8.0  !  1997-07  (G. Madec)   lbc
+   !!            8.1  !  1999-01  (E. Stretta) new option for the mixing length
+   !!  NEMO      1.0  !  2002-06  (G. Madec) add tke_init routine
+   !!             -   !  2002-08  (G. Madec)  rn_cri and Free form, F90
+   !!             -   !  2004-10  (C. Ethe )  1D configuration
+   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
+   !!            3.0  !  2008-05  (C. Ethe,  G.Madec) : update TKE physics:
+   !!                 !           - tke penetration (wind steering)
+   !!                 !           - suface condition for tke & mixing length
+   !!                 !           - Langmuir cells
+   !!             -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
+   !!             -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
+   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term 
+   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl 
+   !!                 !                                + cleaning of the parameters + bugs correction
    !!----------------------------------------------------------------------
-#if defined key_zdftke_old   ||   defined key_esopa
+#if defined key_zdftke   ||   defined key_esopa
    !!----------------------------------------------------------------------
-   !!   'key_zdftke_old'                               TKE vertical physics
+   !!   'key_zdftke'                                   TKE vertical physics
    !!----------------------------------------------------------------------
+   !!   zdf_tke       : update momentum and tracer Kz from a tke scheme
+   !!   tke_tke       : tke time stepping: update tke at now time step (en)
+   !!   tke_avn       : compute mixing length scale and deduce avm and avt
+   !!   tke_init      : initialization, namelist read, and parameters control
+   !!   tke_rst       : read/write tke restart in ocean restart file
    !!----------------------------------------------------------------------
-   !!   zdf_tke_old  : update momentum and tracer Kz from a tke scheme
-   !!   zdf_tke_init : initialization, namelist read, and parameters control
-   !!   tke_rst      : read/write tke restart in ocean restart file
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and active tracers 
-   USE dom_oce         ! ocean space and time domain
-   USE zdf_oce         ! ocean vertical physics
-   USE sbc_oce         ! surface boundary condition: ocean
-   USE phycst          ! physical constants
-   USE zdfmxl          ! mixed layer
-   USE restart         ! only for lrst_oce
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE prtctl          ! Print control
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O manager library
+   USE oce            ! ocean dynamics and active tracers 
+   USE dom_oce        ! ocean space and time domain
+   USE domvvl         ! ocean space and time domain : variable volume layer
+   USE zdf_oce        ! ocean vertical physics
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE phycst         ! physical constants
+   USE zdfmxl         ! mixed layer
+   USE restart        ! only for lrst_oce
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl         ! Print control
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   zdf_tke_old   ! routine called in step module
-
-   LOGICAL , PUBLIC, PARAMETER              ::   lk_zdftke_old = .TRUE.  !: TKE vertical mixing flag
-   REAL(wp), PUBLIC                         ::   eboost              !: multiplicative coeff of the shear product.
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   en                  !: now turbulent kinetic energy
-# if defined key_vectopt_memory
-   !                                                                !!! key_vectopt_memory
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   etmean              !: coefficient used for horizontal smoothing
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   eumean, evmean      !: at t-, u- and v-points
-# endif
+   PUBLIC   zdf_tke    ! routine called in step module
+   PUBLIC   tke_rst    ! routine called in step module
+
+   LOGICAL , PUBLIC, PARAMETER              ::   lk_zdftke = .TRUE.  !: TKE vertical mixing flag
+
 #if defined key_c1d
-   !                                                                !!! 1D cfg only
+   !                                                                !!* 1D cfg only *
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_dis, e_mix        !: dissipation and mixing turbulent lengh scales
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_pdl, e_ric        !: prandl and local Richardson numbers
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   hlc                 !: save finite Langmuir Circulation depth
 #endif
 
+!!gm  - variables to be suppressed from the namelist:
+   LOGICAL  ::   ln_rstke = .FALSE.         ! =T restart with tke from a run without tke
+   INTEGER  ::   nn_itke  = 50              ! number of restart iterative loops
+   INTEGER  ::   nn_ave   =  1              ! horizontal average or not on avt, avmu, avmv (=0/1)
+   REAL(wp) ::   rn_cri   = 2._wp / 9._wp   ! critic Richardson number
+   REAL(wp) ::   rn_efave = 1._wp           ! coefficient for ave : ave=rn_efave*avm
+!!gm end
+
+!!gm  - variables to be added in the namelist (with a larger default value)
+   REAL(wp) ::   rn_bshear = 1.e-20   ! backgrounf shear (>0)
+!!gm
+
    !                                       !!! ** Namelist  namtke  **
-   LOGICAL  ::   ln_rstke = .FALSE.         ! =T restart with tke from a run without tke
    LOGICAL  ::   ln_mxl0  = .FALSE.         ! mixing length scale surface value as function of wind stress or not
-   LOGICAL  ::   ln_lc    = .FALSE.         ! Langmuir cells (LC) as a source term of TKE or not
-   INTEGER  ::   nn_itke  = 50              ! number of restart iterative loops
    INTEGER  ::   nn_mxl   =  2              ! type of mixing length (=0/1/2/3)
+!!gm to be change: the ref value of lmin0 and lmin
+   REAL(wp) ::   rn_lmin0 = 0.4_wp          ! surface  min value of mixing length   [m]
+   REAL(wp) ::   rn_lmin  = 0.1_wp          ! interior min value of mixing length   [m]
    INTEGER  ::   nn_pdl   =  1              ! Prandtl number or not (ratio avt/avm) (=0/1)
-   INTEGER  ::   nn_ave   =  1              ! horizontal average or not on avt, avmu, avmv (=0/1)
    INTEGER  ::   nn_avb   =  0              ! constant or profile background on avt (=0/1)
    REAL(wp) ::   rn_ediff = 0.1_wp          ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
    REAL(wp) ::   rn_ediss = 0.7_wp          ! coefficient of the Kolmogoroff dissipation 
    REAL(wp) ::   rn_ebb   = 3.75_wp         ! coefficient of the surface input of tke
-   REAL(wp) ::   rn_efave = 1._wp           ! coefficient for ave : ave=rn_efave*avm
-   REAL(wp) ::   rn_emin  = 0.7071e-6_wp    ! minimum value of tke (m2/s2)
-   REAL(wp) ::   rn_emin0 = 1.e-4_wp        ! surface minimum value of tke (m2/s2)
-   REAL(wp) ::   rn_cri   = 2._wp / 9._wp   ! critic Richardson number
+   REAL(wp) ::   rn_emin  = 0.7071e-6_wp    ! minimum value of tke           [m2/s2]
+   REAL(wp) ::   rn_emin0 = 1.e-4_wp        ! surface minimum value of tke   [m2/s2]
    INTEGER  ::   nn_havtb = 1               ! horizontal shape or not for avtb (=0/1)
    INTEGER  ::   nn_etau  = 0               ! type of depth penetration of surface tke (=0/1/2)
    INTEGER  ::   nn_htau  = 0               ! type of tke profile of penetration (=0/1/2)
-   REAL(wp) ::   rn_lmin0 = 0.4_wp          ! surface  min value of mixing length
-   REAL(wp) ::   rn_lmin  = 0.1_wp          ! interior min value of mixing length
    REAL(wp) ::   rn_efr   = 1.0_wp          ! fraction of TKE surface value which penetrates in the ocean
+   LOGICAL  ::   ln_lc    = .FALSE.         ! Langmuir cells (LC) as a source term of TKE or not
    REAL(wp) ::   rn_lc    = 0.15_wp         ! coef to compute vertical velocity of Langmuir cells
 
-   REAL(wp), DIMENSION (jpi,jpj) ::   avtb_2d   ! set in tke_init, for other modif than ice
+   REAL(wp) ::   ri_cri   ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
+
+   REAL(wp), DIMENSION(jpi,jpj)     ::   htau      ! depth of tke penetration (nn_htau)
+   REAL(wp), DIMENSION(jpi,jpj)     ::   avtb_2d   ! set in tke_init, for other modif than ice
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   en        ! now turbulent kinetic energy   [m2/s2]
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   dissl     ! now mixing lenght of dissipation
 
    !! * Substitutions
@@ -95,6 +107,6 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
-   !! $Id$
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
+   !! $Id: zdftke2.F90 1201 2008-09-24 13:24:21Z rblod $
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -102,58 +114,39 @@
 CONTAINS
 
-   SUBROUTINE zdf_tke_old( kt )
+   SUBROUTINE zdf_tke( kt )
       !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE zdf_tke_old  ***
+      !!                   ***  ROUTINE zdf_tke  ***
       !!
       !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
-      !!      coefficients using a 1.5 turbulent closure scheme.
-      !!
-      !! ** Method  :   The time evolution of the turbulent kinetic energy 
-      !!      (tke) is computed from a prognostic equation :
-      !!         d(en)/dt = eboost eav (d(u)/dz)**2         ! shear production
-      !!                  + d( rn_efave eav d(en)/dz )/dz   ! diffusion of tke
-      !!                  + grav/rau0 pdl eav d(rau)/dz     ! stratif. destruc.
-      !!                  - rn_ediss / emxl en**(2/3)       ! dissipation
+      !!              coefficients using a turbulent closure scheme (TKE).
+      !!
+      !! ** Method  :   The time evolution of the turbulent kinetic energy (tke)
+      !!              is computed from a prognostic equation :
+      !!         d(en)/dt = avm (d(u)/dz)**2             ! shear production
+      !!                  + d( avm d(en)/dz )/dz         ! diffusion of tke
+      !!                  + avt N^2                      ! stratif. destruc.
+      !!                  - rn_ediss / emxl en**(2/3)    ! Kolmogoroff dissipation
       !!      with the boundary conditions:
       !!         surface: en = max( rn_emin0, rn_ebb sqrt(utau^2 + vtau^2) )
       !!         bottom : en = rn_emin
-      !!      -1- The dissipation and mixing turbulent lengh scales are computed
-      !!         from the usual diagnostic buoyancy length scale:  
-      !!         mxl= sqrt(2*en)/N  where N is the brunt-vaisala frequency
-      !!         with mxl = rn_lmin at the bottom minimum value of 0.4
-      !!      Four cases : 
-      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
-      !!                  zmxld = zmxlm = mxl
-      !!         nn_mxl=1 : mxl bounded by the vertical scale factor.
-      !!                  zmxld = zmxlm = mxl
-      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl 
-      !!                  is less than 1 (|d/dz(xml)|<1). 
-      !!                  zmxld = zmxlm = mxl
-      !!         nn_mxl=3 : lup = mxl bounded using |d/dz(xml)|<1 from the surface 
-      !!                        to the bottom
-      !!                  ldown = mxl bounded using |d/dz(xml)|<1 from the bottom 
-      !!                        to the surface
-      !!                  zmxld = sqrt (lup*ldown) ; zmxlm = min(lup,ldown)
-      !!      -2- Compute the now Turbulent kinetic energy. The time differencing
-      !!      is implicit for vertical diffusion term, linearized for kolmo-
-      !!      goroff dissipation term, and explicit forward for both buoyancy
-      !!      and dynamic production terms. Thus a tridiagonal linear system is
-      !!      solved.
-      !!         Note that - the shear production is multiplied by eboost in order
-      !!      to set the critic richardson number to rn_cri (namelist parameter)
-      !!                   - the destruction by stratification term is multiplied
-      !!      by the Prandtl number (defined by an empirical funtion of the local 
-      !!      Richardson number) if nn_pdl=1 (namelist parameter)
-      !!      coefficient (zesh2):
-      !!      -3- Compute the now vertical eddy vicosity and diffusivity
-      !!      coefficients from en (before the time stepping) and zmxlm:
-      !!              avm = max( avtb, rn_ediff*zmxlm*en^1/2 )
-      !!              avt = max( avmb, pdl*avm )  (pdl=1 if nn_pdl=0)
+      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 
+      !!
+      !!        The now Turbulent kinetic energy is computed using the following 
+      !!      time stepping: implicit for vertical diffusion term, linearized semi
+      !!      implicit for kolmogoroff dissipation term, and explicit forward for 
+      !!      both buoyancy and shear production terms. Therefore a tridiagonal 
+      !!      linear system is solved. Note that buoyancy and shear terms are
+      !!      discretized in a energy conserving form (Bruchard 2002).
+      !!
+      !!        The dissipative and mixing length scale are computed from en and
+      !!      the stratification (see tke_avn)
+      !!
+      !!        The now vertical eddy vicosity and diffusivity coefficients are
+      !!      given by: 
+      !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
+      !!              avt = max( avmb, pdl * avm                 )  
       !!              eav = max( avmb, avm )
-      !!      avt and avm are horizontally averaged to avoid numerical insta-
-      !!      bilities.
-      !!        N.B. The computation is done from jk=2 to jpkm1 except for
-      !!      en. Surface value of avt avmu avmv are set once a time to
-      !!      their background value in routine zdf_tke_init.
+      !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
+      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1 
       !!
       !! ** Action  :   compute en (now turbulent kinetic energy)
@@ -164,40 +157,292 @@
       !!              Mellor and Blumberg, JPO 2004
       !!              Axell, JGR, 2002
+      !!              Bruchard OM 2002
       !!----------------------------------------------------------------------
-      USE oce,     zwd    =>   ua   ! use ua as workspace
-      USE oce,     zmxlm  =>   va   ! use va as workspace
-      USE oce,     zmxld  =>   ta   ! use ta as workspace
-      USE oce,     ztkelc =>   sa   ! use sa as workspace
-      !
       INTEGER, INTENT(in) ::   kt   ! ocean time step
-      !
-      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
-      REAL(wp) ::   zbbrau, zrn2, zesurf            ! temporary scalars
-      REAL(wp) ::   zfact1, ztx2, zdku              !    -         -
-      REAL(wp) ::   zfact2, zty2, zdkv              !    -         -
-      REAL(wp) ::   zfact3, zcoef, zcof, zav        !    -         -
-      REAL(wp) ::   zsh2, zpdl, zri, zsqen, zesh2   !    -         -
-      REAL(wp) ::   zemxl, zemlm, zemlp             !    -         -
-      REAL(wp) ::   zraug, zus, zwlc, zind          !    -         -
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == nit000 )   CALL tke_init     ! initialisation 
+                           !
+                           CALL tke_tke      ! now tke (en)
+                           !
+                           CALL tke_avn      ! now avt, avm, avmu, avmv
+      !
+   END SUBROUTINE zdf_tke
+
+
+   SUBROUTINE tke_tke
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE tke_tke  ***
+      !!
+      !! ** Purpose :   Compute the now Turbulente Kinetic Energy (TKE)
+      !!
+      !! ** Method  : - TKE surface boundary condition
+      !!              - source term due to Langmuir cells (ln_lc=T)
+      !!              - source term due to shear (saved in avmu, avmv arrays)
+      !!              - Now TKE : resolution of the TKE equation by inverting
+      !!                a tridiagonal linear system by a "methode de chasse"
+      !!              - increase TKE due to surface and internal wave breaking
+      !!
+      !! ** Action  : - en : now turbulent kinetic energy)
+      !!              - avmu, avmv : production of TKE by shear at u and v-points
+      !!                (= Kz dz[Ub] * dz[Un] )
+      !! ---------------------------------------------------------------------
+      USE oce,   zdiag  =>   ua   ! use ua as workspace
+      USE oce,   zd_up  =>   va   ! use va as workspace
+      USE oce,   zd_lw  =>   ta   ! use ta as workspace
+      !!
+      INTEGER  ::   ji, jj, jk                ! dummy loop arguments
+      REAL(wp) ::   zbbrau, zesurf, zesh2     ! temporary scalars
+      REAL(wp) ::   zfact1, zfact2, zfact3    !    -         -
+      REAL(wp) ::   ztx2  , zty2  , zcof      !    -         -
+      REAL(wp) ::   zus   , zwlc  , zind      !    -         -
+      REAL(wp) ::   zzd_up, zzd_lw            !    -         -
       INTEGER , DIMENSION(jpi,jpj)     ::   imlc    ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zhtau   !  -      -
       REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc    !  -      -
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc   ! 3D workspace
       !!--------------------------------------------------------------------
-
-      IF( kt == nit000  )   CALL zdf_tke_init      ! Initialization (first time-step only)
-
-      !                                            ! Local constant initialization
-      zbbrau =  .5 * rn_ebb / rau0
-      zfact1 = -.5 * rdt * rn_efave
+      !
+      zbbrau =  .5 * rn_ebb / rau0       ! Local constant initialisation
+      zfact1 = -.5 * rdt 
       zfact2 = 1.5 * rdt * rn_ediss
-      zfact3 = 0.5 * rdt * rn_ediss
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! I.  Mixing length
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-      ! Buoyancy length scale: l=sqrt(2*e/n**2)
-      ! ---------------------
+      zfact3 = 0.5       * rn_ediss
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Surface boundary condition on tke
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!CDIR NOVERRCHK
+      DO jj = 2, jpjm1         ! en(1)   = rn_ebb sqrt(utau^2+vtau^2) / rau0  (min value rn_emin0)
+!CDIR NOVERRCHK
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ztx2   = utau(ji-1,jj  ) + utau(ji,jj)
+            zty2   = vtau(ji  ,jj-1) + vtau(ji,jj)
+            zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 )
+            en(ji,jj,1) = MAX( zesurf, rn_emin0 ) * tmask(ji,jj,1)
+         END DO
+      END DO
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Bottom boundary condition on tke
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!!gm to be added soon
+      !
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke
+         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+         !
+         !                        !* total energy produce by LC : cumulative sum over jk
+         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1)
+         DO jk = 2, jpk
+            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
+         END DO
+         !                        !* finite Langmuir Circulation depth
+         imlc(:,:) = mbathy(:,:)         ! Initialization to the number of w ocean point mbathy
+         DO jk = jpkm1, 2, -1
+            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us 
+               DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
+                  zus  = 0.000128 * wndm(ji,jj) * wndm(ji,jj)
+                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
+               END DO
+            END DO
+         END DO
+         !                               ! finite LC depth
+# if defined key_vectopt_loop
+         DO jj = 1, 1
+            DO ji = 1, jpij   ! vector opt. (forced unrolling)
+# else
+         DO jj = 1, jpj 
+            DO ji = 1, jpi
+# endif
+               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
+            END DO
+         END DO
+# if defined key_c1d
+         hlc(:,:) = zhlc(:,:) * tmask(:,:,1)      ! c1d configuration: save finite Langmuir Circulation depth
+# endif
+!CDIR NOVERRCHK
+         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en
+!CDIR NOVERRCHK
+            DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+!!gm replace here wndn by a formulation with the stress module
+                  zus  = 0.016 * wndm(ji,jj)                  ! Stokes drift
+                  !                                           ! vertical velocity due to LC
+                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
+                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
+                  !                                           ! TKE Langmuir circulation source term
+                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
+      ENDIF
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Now Turbulent kinetic energy (output in en)
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
+      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
+      !                     ! 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, 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
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zcof   = zfact1 * tmask(ji,jj,jk)
+               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal
+                  &          / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk  ) )
+               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  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
+                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , 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.e0 - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(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)  ) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
+      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1    ! vector opt.
+               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+         DO ji = fs_2, fs_jpim1    ! vector opt.
+            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
+         END DO
+      END DO
+      DO jk = 3, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1    ! vector opt.
+               zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+         DO ji = fs_2, fs_jpim1    ! vector opt.
+            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
+         END DO
+      END DO
+      DO jk = jpk-2, 2, -1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1    ! vector opt.
+               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      DO jk = 2, jpkm1                             ! set the minimum value of tke
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  TKE due to surface and internal wave breaking
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( nn_etau == 1 ) THEN           !*  penetration throughout the water column
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
+                     &                                               * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               jk = nmln(ji,jj)
+               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
+                  &                                               * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
+            END DO
+         END DO
+      ENDIF
+      !
+      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
+      !
+   END SUBROUTINE tke_tke
+
+
+   SUBROUTINE tke_avn
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE tke_avn  ***
+      !!
+      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
+      !!
+      !! ** Method  :   At this stage, en, the now TKE, is known (computed in 
+      !!              the tke_tke routine). First, the now mixing lenth is 
+      !!      computed from en and the strafification (N^2), then the mixings
+      !!      coefficients are computed.
+      !!              - Mixing length : a first evaluation of the mixing lengh
+      !!      scales is:
+      !!                      mxl = sqrt(2*en) / N  
+      !!      where N is the brunt-vaisala frequency, with a minimum value set
+      !!      to rn_lmin (rn_lmin0) in the interior (surface) ocean.
+      !!        The mixing and dissipative length scale are bound as follow : 
+      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
+      !!                        zmxld = zmxlm = mxl
+      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
+      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 
+      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
+      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
+      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to 
+      !!                    the surface to obtain ldown. the resulting length 
+      !!                    scales are:
+      !!                        zmxld = sqrt( lup * ldown ) 
+      !!                        zmxlm = min ( lup , ldown )
+      !!              - Vertical eddy viscosity and diffusivity:
+      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
+      !!                      avt = max( avmb, pdlr * avm )  
+      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
+      !!
+      !! ** Action  : - avt : now vertical eddy diffusivity (w-point)
+      !!              - avmu, avmv : now vertical eddy viscosity at uw- and vw-points
+      !!----------------------------------------------------------------------
+      USE oce,     zmpdl  =>   ua   ! use ua as workspace
+      USE oce,     zmxlm  =>   va   ! use va as workspace
+      USE oce,     zmxld  =>   ta   ! use ta as workspace
+      !!
+      INTEGER  ::   ji, jj, jk            ! dummy loop arguments
+      REAL(wp) ::   zrn2, zraug           ! temporary scalars
+      REAL(wp) ::   ztx2, zdku            !    -         -
+      REAL(wp) ::   zty2, zdkv            !    -         -
+      REAL(wp) ::   zcoef, zav            !    -         -
+      REAL(wp) ::   zpdlr, zri, zsqen     !    -         -
+      REAL(wp) ::   zemxl, zemlm, zemlp   !    -         -
+      !!--------------------------------------------------------------------
+
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Mixing length
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
+      !
       IF( ln_mxl0 ) THEN         ! surface mixing length = F(stress) : l=vkarmn*2.e5*sqrt(utau^2 + vtau^2)/(rau0*g)
 !!gm  this should be useless
@@ -224,15 +469,17 @@
 !CDIR NOVERRCHK
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zrn2 = MAX( rn2b(ji,jj,jk), rsmall )
+               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
                zmxlm(ji,jj,jk) = MAX(  rn_lmin,  SQRT( 2. * en(ji,jj,jk) / zrn2 )  )
             END DO
          END DO
       END DO
-
-      ! Physical limits for the mixing length
-      ! -------------------------------------
+      !
+!!gm  CAUTION: to be added here the bottom boundary condition on zmxlm
+      !
+      !                     !* Physical limits for the mixing length
+      !
       zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value
       zmxld(:,:,jpk) = rn_lmin        ! bottom  set to the minimum value
-
+      !
       SELECT CASE ( nn_mxl )
       !
@@ -308,74 +555,15 @@
          !
       END SELECT
-
+      !
 # if defined key_c1d
-      ! c1d configuration : save mixing and dissipation turbulent length scales
-      e_dis(:,:,:) = zmxld(:,:,:)
+      e_dis(:,:,:) = zmxld(:,:,:)      ! c1d configuration : save mixing and dissipation turbulent length scales
       e_mix(:,:,:) = zmxlm(:,:,:)
 # endif
 
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! II  TKE Langmuir circulation source term
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_lc ) THEN
-         !
-         ! Computation of total energy produce by LC : cumulative sum over jk
-         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1)
-         DO jk = 2, jpk
-            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
-         END DO
-         !
-         ! Computation of finite Langmuir Circulation depth
-         ! Initialization to the number of w ocean point mbathy
-         imlc(:,:) = mbathy(:,:)
-         DO jk = jpkm1, 2, -1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ! Last w-level at which zpelc>=0.5*us*us with us=0.016*wind(starting from jpk-1)
-                  zus  = 0.000128 * wndm(ji,jj) * wndm(ji,jj)
-                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
-               END DO
-            END DO
-         END DO
-         !
-         ! finite LC depth
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
-            END DO
-         END DO
-         !
-# if defined key_c1d
-         hlc(:,:) = zhlc(:,:) * tmask(:,:,1)      ! c1d configuration: save finite Langmuir Circulation depth
-# endif
-         !
-         ! TKE Langmuir circulation source term
-!CDIR NOVERRCHK
-         DO jk = 2, jpkm1
-!CDIR NOVERRCHK
-            DO jj = 2, jpjm1
-!CDIR NOVERRCHK
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ! Stokes drift
-                  zus  = 0.016 * wndm(ji,jj)
-                  ! computation of vertical velocity due to LC
-                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
-                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
-                  ! TKE Langmuir circulation source term
-                  ztkelc(ji,jj,jk) = ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)
-               END DO
-            END DO
-         END DO
-         !
-      ENDIF
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! III  Tubulent kinetic energy time stepping
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-      ! 1. Vertical eddy viscosity on tke (put in zmxlm) and first estimate of avt
-      ! ---------------------------------------------------------------------
-!CDIR NOVERRCHK
-      DO jk = 2, jpkm1
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt)
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!CDIR NOVERRCHK
+      DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
 !CDIR NOVERRCHK
          DO jj = 2, jpjm1
@@ -384,317 +572,47 @@
                zsqen = SQRT( en(ji,jj,jk) )
                zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
+               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * tmask(ji,jj,jk)
                avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
-               zmxlm(ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)
-               zmxld(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-
-      ! 2. Surface boundary condition on tke and its eddy viscosity (zmxlm)
-      ! -------------------------------------------------
-      ! en(1)   = rn_ebb sqrt(utau^2+vtau^2) / rau0  (min value rn_emin0)
-      ! zmxlm(1) = avmb(1) and zmxlm(jpk) = 0.
-!CDIR NOVERRCHK
-      DO jj = 2, jpjm1
-!CDIR NOVERRCHK
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
-            zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
-            zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 )
-            en (ji,jj,1) = MAX( zesurf, rn_emin0 ) * tmask(ji,jj,1)
-            zav  =  rn_ediff * zmxlm(ji,jj,1) * SQRT( en(ji,jj,1) )
-            zmxlm(ji,jj,1  ) = MAX( zav, avmb(1)                  ) * tmask(ji,jj,1)
-            avt  (ji,jj,1  ) = MAX( zav, avtb(1) * avtb_2d(ji,jj) ) * tmask(ji,jj,1)
-            zmxlm(ji,jj,jpk) = 0.e0
-         END DO
-      END DO
-
-      ! 3. Now Turbulent kinetic energy (output in en)
-      ! -------------------------------
-      ! Resolution of a tridiagonal linear system by a "methode de chasse"
-      ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ). 
-
-      SELECT CASE ( nn_pdl )
-      !
-      CASE ( 0 )           ! No Prandtl number
+               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
+      !
+      DO jk = 2, jpkm1            !* vertical eddy viscosity at u- and v-points
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
+               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
+      !
+      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  !                                                             ! shear prod. - stratif. destruction
-                  zcoef = 0.5 / fse3w(ji,jj,jk)
-                  zdku = zcoef * (  ub(ji-1, jj ,jk-1) + ub(ji,jj,jk-1)   &          ! shear
-                     &            - ub(ji-1, jj ,jk  ) - ub(ji,jj,jk  )  )
-                  zdkv = zcoef * (  vb( ji ,jj-1,jk-1) + vb(ji,jj,jk-1)   &
-                     &            - vb( ji ,jj-1,jk  ) - vb(ji,jj,jk  )  )
-                  zesh2 =  eboost * ( zdku*zdku + zdkv*zdkv ) - rn2b(ji,jj,jk)        ! coefficient (zesh2)
-                  !
-                  !                                                             ! Matrix
-                  zcof = zfact1 * tmask(ji,jj,jk)
-                  !                                                                  ! lower diagonal
-                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   &
-                     &                  / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
-                  !                                                                  ! upper diagonal
-                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   &
-                     &                  / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
-                  !                                                                  ! diagonal
-                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
-                  !
-                  !                                                             ! right hand side in en 
-                  IF( .NOT. ln_lc ) THEN                                             ! No Langmuir cells
-                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en(ji,jj,jk)   &
-                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2
-                  ELSE                                                               ! Langmuir cell source term
-                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en(ji,jj,jk)   &
-                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2          &
-                        &                        +   rdt  * ztkelc(ji,jj,jk)
-                  ENDIF
-               END DO
-            END DO
-         END DO
-         !
-      CASE ( 1 )           ! Prandtl number
-         DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  !                                                             ! shear prod. - stratif. destruction
-                  zcoef = 0.5 / fse3w(ji,jj,jk)                                     
-                  zdku = zcoef * (  ub(ji-1,jj  ,jk-1) + ub(ji,jj,jk-1)   &          ! shear
-                  &               - ub(ji-1,jj  ,jk  ) - ub(ji,jj,jk  )  )
-                  zdkv = zcoef * (  vb(ji  ,jj-1,jk-1) + vb(ji,jj,jk-1)   &
-                  &               - vb(ji  ,jj-1,jk  ) - vb(ji,jj,jk  )  )
-                  zsh2 = zdku * zdku + zdkv * zdkv                                   ! square of shear
-                  zri  = MAX( rn2b(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 )                ! local Richardson number
+                  zcoef = 0.5 / ( 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. ) * avm(ji,jj,jk) / (zcoef * (zdku + zdkv + rn_bshear ) )
+                  zpdlr = MAX(  0.1,  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,  ri_crit / MAX( ri_crit , zri )  )
+                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
 # if defined key_c1d
-                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri 
+                  e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number
+                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri
 # endif
-                  zpdl = 1.0                                                         ! Prandtl number
-                  IF( zri >= 0.2 )   zpdl = 0.2 / zri
-                  zpdl = MAX( 0.1, zpdl )
-                  zesh2 = eboost * zsh2 - zpdl * rn2b(ji,jj,jk)                       ! coefficient (esh2)
-                  !
-                  !                                                             ! Matrix
-                  zcof = zfact1 * tmask(ji,jj,jk)
-                  !                                                                 ! lower diagonal
-                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   &
-                  &                     / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
-                  !                                                                 ! upper diagonal
-                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   &
-                  &                     / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
-                  !                                                                 ! diagonal
-                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
-                  !
-                  !                                                             ! right hand side in en 
-                  IF( .NOT. ln_lc ) THEN                                             ! No Langmuir cells
-                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en   (ji,jj,jk)   &
-                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2
-                  ELSE                                                               ! Langmuir cell source term
-                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en   (ji,jj,jk)   &
-                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2             &
-                        &                        +   rdt  * ztkelc(ji,jj,jk)
-                  ENDIF
-                  zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk)                      ! store masked Prandlt number in zmxld array
-               END DO
-            END DO
-         END DO
-         !
-      END SELECT
-
-# if defined key_c1d
-      e_pdl(:,:,2:jpkm1) = zmxld(:,:,2:jpkm1)      ! c1d configuration : save masked Prandlt number
-      e_pdl(:,:,      1) = e_pdl(:,:,      2)
-      e_pdl(:,:,    jpk) = e_pdl(:,:,  jpkm1)      
-# endif
-
-      ! 4. Matrix inversion from level 2 (tke prescribed at level 1)
-      !!--------------------------------
-      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
-               zwd(ji,jj,jk) = zwd(ji,jj,jk) - avmv(ji,jj,jk) * avmu(ji,jj,jk-1) / zwd(ji,jj,jk-1)
-            END DO
-         END DO
-      END DO
-      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
-         DO ji = fs_2, fs_jpim1    ! vector opt.
-            avmv(ji,jj,2) = en(ji,jj,2) - avmv(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
-         END DO
-      END DO
-      DO jk = 3, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
-               avmv(ji,jj,jk) = en(ji,jj,jk) - avmv(ji,jj,jk) / zwd(ji,jj,jk-1) *avmv(ji,jj,jk-1)
-            END DO
-         END DO
-      END DO
-      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
-         DO ji = fs_2, fs_jpim1    ! vector opt.
-            en(ji,jj,jpkm1) = avmv(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
-         END DO
-      END DO
-      DO jk = jpk-2, 2, -1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
-               en(ji,jj,jk) = ( avmv(ji,jj,jk) - avmu(ji,jj,jk) * en(ji,jj,jk+1) ) / zwd(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      DO jk = 2, jpkm1                             ! set the minimum value of tke
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-     
-      ! 5. Add extra TKE due to surface and internal wave breaking (nn_etau /= 0)
-      !!----------------------------------------------------------
-      IF( nn_etau /= 0 ) THEN        ! extra tke : en = en + rn_efr * en(1) * exp( -z/zhtau )
-         !
-         SELECT CASE( nn_htau )           ! Choice of the depth of penetration
-         CASE( 0 )                                    ! constant depth penetration (here 10 meters)
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zhtau(ji,jj) = 10.
-               END DO
-            END DO
-         CASE( 1 )                                    ! meridional profile  1
-            DO jj = 2, jpjm1                          ! ( 5m in the tropics to a maximum of 40 m at high lat.)
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zhtau(ji,jj) = MAX( 5., MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
-               END DO
-            END DO
-         CASE( 2 )                                    ! meridional profile 2
-            DO jj = 2, jpjm1                          ! ( 5m in the tropics to a maximum of 60 m at high lat.)
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zhtau(ji,jj) =  MAX( 5.,6./4.* MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
-               END DO
-            END DO
-         END SELECT
-         ! 
-         IF( nn_etau == 1 ) THEN          ! extra term throughout the water column
-            DO jk = 2, jpkm1
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-                     en(ji,jj,jk) = en(ji,jj,jk)   &
-                        &         + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) &
-                        &                  * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-         ELSEIF( nn_etau == 2 ) THEN      ! extra term only at the base of the mixed layer
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  jk = nmln(ji,jj)
-                  en(ji,jj,jk) = en(ji,jj,jk)    &
-                     &         + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) &
-                     &                   * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
-               END DO
-            END DO
-         ENDIF
-         !
+              END DO
+            END DO
+         END DO
       ENDIF
-
-
-      ! Lateral boundary conditions (sign unchanged)
-      CALL lbc_lnk( en, 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. )
-
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! IV.  Before vertical eddy vicosity and diffusivity coefficients
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !
-      SELECT CASE ( nn_ave )
-      CASE ( 0 )                                     ! no horizontal average 
-         DO jk = 2, jpkm1                                 ! only vertical eddy viscosity 
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  avmu(ji,jj,jk) = ( avt  (ji,jj,jk) + avt  (ji+1,jj  ,jk) ) * umask(ji,jj,jk)   &
-                  &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk) )
-                  avmv(ji,jj,jk) = ( avt  (ji,jj,jk) + avt  (ji  ,jj+1,jk) ) * vmask(ji,jj,jk)   &
-                  &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk) )
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
-         !
-      CASE ( 1 )                                     ! horizontal average                         ( 1/2  1/2 )
-         !                                                ! Vertical eddy viscosity    avmu = 1/4 ( 1    1   )
-         !                                                !                                       ( 1/2  1/2 )
-         !                                                !
-         !                                                !                                       ( 1/2  1   1/2 )            
-         !                                                !                            avmv = 1/4 ( 1/2  1   1/2 )      
-         DO jk = 2, jpkm1                                
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-# if defined key_vectopt_memory
-                  !                                       ! caution : vectopt_memory change the solution
-                  !                                       !           (last digit of the solver stat)
-                  avmu(ji,jj,jk) = (      avt(ji,jj  ,jk) + avt(ji+1,jj  ,jk)   &
-                     &              +.5*( avt(ji,jj-1,jk) + avt(ji+1,jj-1,jk)   &
-                     &                   +avt(ji,jj+1,jk) + avt(ji+1,jj+1,jk) ) ) * eumean(ji,jj,jk)
-
-                  avmv(ji,jj,jk) = (      avt(ji  ,jj,jk) + avt(ji  ,jj+1,jk)   &
-                     &              +.5*( avt(ji-1,jj,jk) + avt(ji-1,jj+1,jk)   &
-                     &                   +avt(ji+1,jj,jk) + avt(ji+1,jj+1,jk) ) ) * evmean(ji,jj,jk)
-# else
-                  avmu(ji,jj,jk) = (   avt  (ji,jj  ,jk) + avt  (ji+1,jj  ,jk)   &
-                     &           +.5*( avt  (ji,jj-1,jk) + avt  (ji+1,jj-1,jk)   &
-                     &                +avt  (ji,jj+1,jk) + avt  (ji+1,jj+1,jk) ) ) * umask(ji,jj,jk)   &
-                     &    / MAX( 1.,   tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk)   &
-                     &           +.5*( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk)   &
-                     &                +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) )  )
-
-                  avmv(ji,jj,jk) = (   avt  (ji  ,jj,jk) + avt  (ji  ,jj+1,jk)   &
-                     &           +.5*( avt  (ji-1,jj,jk) + avt  (ji-1,jj+1,jk)   &
-                     &                +avt  (ji+1,jj,jk) + avt  (ji+1,jj+1,jk) ) ) * vmask(ji,jj,jk)   &
-                     &   /  MAX( 1.,   tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
-                     &           +.5*( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk)   &
-                     &                +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) )  )
-# endif
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions 
-         !
-         !                                                ! Vertical eddy diffusivity             (1 2 1)
-         !                                                !                            avt = 1/16 (2 4 2)
-         !                                                !                                       (1 2 1)
-         DO jk = 2, jpkm1         
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-#  if defined key_vectopt_memory
-                  avt(ji,jj,jk) = ( avmu(ji,jj,jk) + avmu(ji-1,jj  ,jk)                          &
-                     &            + avmv(ji,jj,jk) + avmv(ji  ,jj-1,jk)  ) * etmean(ji,jj,jk)
-#  else
-                  avt(ji,jj,jk) = ( avmu (ji,jj,jk) + avmu (ji-1,jj  ,jk)                        &
-                     &            + avmv (ji,jj,jk) + avmv (ji  ,jj-1,jk)  ) * tmask(ji,jj,jk)   &
-                     &  / MAX( 1.,  umask(ji,jj,jk) + umask(ji-1,jj  ,jk)                        &
-                     &            + vmask(ji,jj,jk) + vmask(ji  ,jj-1,jk)  )
-#  endif
-               END DO
-            END DO
-         END DO
-         !
-      END SELECT
-      !
-      IF( nn_pdl == 1 ) THEN                            ! Ponderation by the Prandtl number (nn_pdl=1)
-         DO jk = 2, jpkm1 
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zpdl = zmxld(ji,jj,jk)
-                  avt(ji,jj,jk) = MAX( zpdl * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-      ENDIF
-      !
-      DO jk = 2, jpkm1                                  ! Minimum value on the eddy viscosity
-         avmu(:,:,jk) = MAX( avmu(:,:,jk), avmb(jk) ) * umask(:,:,jk)
-         avmv(:,:,jk) = MAX( avmv(:,:,jk), avmb(jk) ) * vmask(:,:,jk)
-      END DO
-
       CALL lbc_lnk( avt, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged)
-
-      IF( lrst_oce )   CALL tke_rst( kt, 'WRITE' )      ! write en in restart file
 
       IF(ln_ctl) THEN
@@ -704,30 +622,22 @@
       ENDIF
       !
-   END SUBROUTINE zdf_tke_old
-
-
-   SUBROUTINE zdf_tke_init
+   END SUBROUTINE tke_avn
+
+
+   SUBROUTINE tke_init
       !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE zdf_tke_init  ***
+      !!                  ***  ROUTINE tke_init  ***
       !!                     
       !! ** Purpose :   Initialization of the vertical eddy diffivity and 
-      !!      viscosity when using a tke turbulent closure scheme
+      !!              viscosity when using a tke turbulent closure scheme
       !!
       !! ** Method  :   Read the namtke namelist and check the parameters
-      !!      called at the first timestep (nit000)
+      !!              called at the first timestep (nit000)
       !!
       !! ** input   :   Namlist namtke
       !!
       !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
-      !!
       !!----------------------------------------------------------------------
-      USE dynzdf_exp
-      USE trazdf_exp
-      !
-# if defined key_vectopt_memory
       INTEGER ::   ji, jj, jk   ! dummy loop indices
-# else
-      INTEGER ::           jk   ! dummy loop indices
-# endif
       !!
       NAMELIST/namtke/ ln_rstke, rn_ediff, rn_ediss, rn_ebb  , rn_efave, rn_emin,   &
@@ -737,23 +647,13 @@
       !!----------------------------------------------------------------------
 
-      ! Read Namelist namtke : Turbulente Kinetic Energy
-      ! --------------------
-      REWIND ( numnam )
+      REWIND ( numnam )               !* Read Namelist namtke : Turbulente Kinetic Energy
       READ   ( numnam, namtke )
-
-      ! Compute boost associated with the Richardson critic
-      !     (control values: rn_cri = 0.3   ==> eboost=1.25 for nn_pdl=1)
-      !     (                rn_cri = 0.222 ==> eboost=1.               )
-      eboost = rn_cri * ( 2. + rn_ediss / rn_ediff ) / 2.
-
-
-
-      ! Parameter control and print
-      ! ---------------------------
-      ! Control print
-      IF(lwp) THEN
+      
+      ri_cri = 2. / ( 2. + rn_ediss / rn_ediff )      ! resulting critical Richardson number
+
+      IF(lwp) THEN                    !* Control print
          WRITE(numout,*)
-         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme (old scheme)'
-         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) 'zdf_tke : tke turbulent closure scheme - initialisation'
+         WRITE(numout,*) '~~~~~~~~'
          WRITE(numout,*) '          Namelist namtke : set tke mixing parameters'
          WRITE(numout,*) '             restart with tke from no tke              ln_rstke = ', ln_rstke
@@ -761,13 +661,8 @@
          WRITE(numout,*) '             Kolmogoroff dissipation coef.             rn_ediss = ', rn_ediss
          WRITE(numout,*) '             tke surface input coef.                   rn_ebb   = ', rn_ebb
-         WRITE(numout,*) '             tke diffusion coef.                       rn_efave = ', rn_efave
          WRITE(numout,*) '             minimum value of tke                      rn_emin  = ', rn_emin
          WRITE(numout,*) '             surface minimum value of tke              rn_emin0 = ', rn_emin0
-         WRITE(numout,*) '             number of restart iter loops              nn_itke  = ', nn_itke
          WRITE(numout,*) '             mixing length type                        nn_mxl   = ', nn_mxl
          WRITE(numout,*) '             prandl number flag                        nn_pdl   = ', nn_pdl
-         WRITE(numout,*) '             horizontal average flag                   nn_ave   = ', nn_ave
-         WRITE(numout,*) '             critic Richardson nb                      rn_cri   = ', rn_cri
-         WRITE(numout,*) '                and its associated coeff.              eboost   = ', eboost
          WRITE(numout,*) '             constant background or profile            nn_avb   = ', nn_avb
          WRITE(numout,*) '             surface mixing length = F(stress) or not  ln_mxl0  = ', ln_mxl0
@@ -781,83 +676,20 @@
          WRITE(numout,*) '             coef to compute verticla velocity of LC   rn_lc    = ', rn_lc
          WRITE(numout,*)
+         WRITE(numout,*) '             critical Richardson nb with your choice of coefs.  = ', ri_cri
       ENDIF
 
-      ! Check of some namelist values
+      !                               !* Check of some namelist values
       IF( nn_mxl  < 0    .OR. nn_mxl  > 3   )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
       IF( nn_pdl  < 0    .OR. nn_pdl  > 1   )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
-      IF( nn_ave  < 0    .OR. nn_ave  > 1   )   CALL ctl_stop( 'bad flag: nn_ave is  0 or 1    ' )
       IF( nn_htau < 0    .OR. nn_htau > 2   )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
       IF( rn_lc   < 0.15 .OR. rn_lc   > 0.2 )   CALL ctl_stop( 'bad value: rn_lc must be between 0.15 and 0.2 ' )
 
-      IF( nn_etau == 2 )   CALL zdf_mxl( nit000 )      ! Initialization of nmln 
-
-      ! Horizontal average : initialization of weighting arrays 
-      ! -------------------
-      !
-      SELECT CASE ( nn_ave )
-      ! Notice: mean arrays have not to by defined at local domain boundaries due to the vertical nature of TKE
-      !
-      CASE ( 0 )                ! no horizontal average
-         IF(lwp) WRITE(numout,*) '          no horizontal average on avt, avmu, avmv'
-         IF(lwp) WRITE(numout,*) '          only in very high horizontal resolution !'
-# if defined key_vectopt_memory
-         ! caution vectopt_memory change the solution (last digit of the solver stat)
-         ! weighting mean arrays etmean, eumean and evmean
-         !           ( 1  1 )                                          ( 1 )
-         ! avt = 1/4 ( 1  1 )     avmu = 1/2 ( 1  1 )       avmv=  1/2 ( 1 )
-         !                          
-         etmean(:,:,:) = 0.e0
-         eumean(:,:,:) = 0.e0
-         evmean(:,:,:) = 0.e0
-         !
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  etmean(ji,jj,jk) = tmask(ji,jj,jk) / MAX( 1.,  umask(ji,jj,jk) + umask(ji-1,jj  ,jk)   &
-                  &                                            + vmask(ji,jj,jk) + vmask(ji  ,jj-1,jk) )
-                  eumean(ji,jj,jk) = umask(ji,jj,jk) / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk) )
-                  evmean(ji,jj,jk) = vmask(ji,jj,jk) / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk) )
-               END DO
-            END DO
-         END DO
-# endif
-         !
-      CASE ( 1 )                ! horizontal average 
-         IF(lwp) WRITE(numout,*) '          horizontal average on avt, avmu, avmv'
-# if defined key_vectopt_memory
-         ! caution vectopt_memory change the solution (last digit of the solver stat)
-         ! weighting mean arrays etmean, eumean and evmean
-         !           ( 1  1 )              ( 1/2  1/2 )             ( 1/2  1  1/2 )
-         ! avt = 1/4 ( 1  1 )   avmu = 1/4 ( 1    1   )   avmv= 1/4 ( 1/2  1  1/2 )
-         !                                 ( 1/2  1/2 )
-         etmean(:,:,:) = 0.e0
-         eumean(:,:,:) = 0.e0
-         evmean(:,:,:) = 0.e0
-         !
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  etmean(ji,jj,jk) = tmask(ji,jj,jk) / MAX( 1.,   umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk)   &
-                  &                                       +       vmask(ji  ,jj-1,jk) + vmask(ji  ,jj  ,jk) )
-                  eumean(ji,jj,jk) = umask(ji,jj,jk) / MAX( 1.,   tmask(ji  ,jj  ,jk) + tmask(ji+1,jj  ,jk)   &
-                  &                                       +.5 * ( tmask(ji  ,jj-1,jk) + tmask(ji+1,jj-1,jk)   &
-                  &                                             + tmask(ji  ,jj+1,jk) + tmask(ji+1,jj+1,jk) )  )
-                  evmean(ji,jj,jk) = vmask(ji,jj,jk) / MAX( 1.,   tmask(ji  ,jj  ,jk) + tmask(ji  ,jj+1,jk)   &
-                  &                                       +.5 * ( tmask(ji-1,jj  ,jk) + tmask(ji-1,jj+1,jk)   &
-                  &                                             + tmask(ji+1,jj  ,jk) + tmask(ji+1,jj+1,jk) )  )
-               END DO
-            END DO
-         END DO
-# endif
-         !
-      END SELECT
-
-
-      ! Background eddy viscosity and diffusivity profil
-      ! ------------------------------------------------
-      IF( nn_avb == 0 ) THEN      ! Define avmb, avtb from namelist parameter
+      IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln 
+
+      !                               !* Background eddy viscosity and diffusivity profil
+      IF( nn_avb == 0 ) THEN                ! Define avmb, avtb from namelist parameter
          avmb(:) = avm0
          avtb(:) = avt0
-      ELSE                      ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990) 
+      ELSE                                  ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990) 
          avmb(:) = avm0
          avtb(:) = avt0 + ( 3.0e-4 - 2 * avt0 ) * 1.0e-4 * gdepw_0(:)   ! m2/s
@@ -865,8 +697,8 @@
       ENDIF
       !
-      !                         ! 2D shape of the avtb
-      avtb_2d(:,:) = 1.e0             ! uniform 
-      !
-      IF( nn_havtb == 1 ) THEN      ! decrease avtb in the equatorial band
+      !                                     ! 2D shape of the avtb
+      avtb_2d(:,:) = 1.e0                        ! uniform 
+      !
+      IF( nn_havtb == 1 ) THEN                   ! decrease avtb in the equatorial band
            !  -15S -5S : linear decrease from avt0 to avt0/10.
            !  -5S  +5N : cst value avt0/10.
@@ -877,23 +709,41 @@
       ENDIF
 
-      ! Initialization of vertical eddy coef. to the background value
-      ! -------------------------------------------------------------
+      !                               !* depth of penetration of surface tke
+      IF( nn_etau /= 0 ) THEN      
+         SELECT CASE( nn_htau )           ! Choice of the depth of penetration
+         CASE( 0 )                                    ! constant depth penetration (here 10 meters)
+            htau(:,:) = 10.e0
+         CASE( 1 )                                    ! F(latitude) : 5m to 40m at high lat.
+            DO jj = 1, jpj
+               DO ji = 1, jpi  
+                  htau(ji,jj) = MAX( 5., MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
+               END DO
+            END DO
+         CASE( 2 )                                    ! F(latitude) : 5m to 60m at high lat.
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  htau(ji,jj) =  MAX( 5.,6./4.* MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
+               END DO
+            END DO
+         END SELECT
+      ENDIF
+
+      !                               !* set vertical eddy coef. to the background value
       DO jk = 1, jpk
          avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
+         avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
          avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
          avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
       END DO
-
-
-      ! read or initialize turbulent kinetic energy ( en )
-      ! -------------------------------------------------
+      dissl(:,:,:) = 1.e-12
+      !                               !* read or initialize all required files 
       CALL tke_rst( nit000, 'READ' )
       !
-   END SUBROUTINE zdf_tke_init
+   END SUBROUTINE tke_init
 
 
    SUBROUTINE tke_rst( kt, cdrw )
      !!---------------------------------------------------------------------
-     !!                   ***  ROUTINE ts_rst  ***
+     !!                   ***  ROUTINE tke_rst  ***
      !!                     
      !! ** Purpose :   Read or write TKE file (en) in restart file
@@ -906,26 +756,55 @@
      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
      !
-     INTEGER ::   jit   ! dummy loop indices
+     INTEGER ::   jit, jk   ! dummy loop indices
+     INTEGER ::   id1, id2, id3, id4, id5, id6
      !!----------------------------------------------------------------------
      !
-     IF( TRIM(cdrw) == 'READ' ) THEN
-        IF( ln_rstart ) THEN
-           IF( iom_varid( numror, 'en', ldstop = .FALSE. ) > 0 .AND. .NOT.(ln_rstke) ) THEN 
+     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+        !                                   ! ---------------
+        IF( ln_rstart ) THEN                   !* Read the restart file
+           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
+           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
+           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
+           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
+           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
+           id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
+           !
+           IF( id1 > 0 ) THEN                       ! 'en' exists
               CALL iom_get( numror, jpdom_autoglo, 'en', en )
-           ELSE
-              IF( lwp .AND. iom_varid( numror, 'en', ldstop = .FALSE. ) > 0 )   &
-                 &                       WRITE(numout,*) ' ===>>>> : previous run without tke scheme'
-              IF( lwp .AND. ln_rstke )   WRITE(numout,*) ' ===>>>> : We do not use en from the restart file'
-              IF( lwp                )   WRITE(numout,*) ' ===>>>> : en is set by iterative loop'
+              IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
+                 CALL iom_get( numror, jpdom_autoglo, 'avt'  , avt   )
+                 CALL iom_get( numror, jpdom_autoglo, 'avm'  , avm   )
+                 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu  )
+                 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv  )
+                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
+              ELSE                                                 ! one at least array is missing
+                 CALL tke_avn                                          ! compute avt, avm, avmu, avmv and dissl (approximation)
+              ENDIF
+           ELSE                                     ! No TKE array found: initialisation
+              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
               en (:,:,:) = rn_emin * tmask(:,:,:)
-              DO jit = 2, nn_itke + 1
-                 CALL zdf_tke_old( jit )
-              END DO
+              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
+              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
            ENDIF
-        ELSE
-           en(:,:,:) = rn_emin * tmask(:,:,:)      ! no restart: en set to its minimum value
+        ELSE                                   !* Start from rest
+           en(:,:,:) = rn_emin * tmask(:,:,:)
+           DO jk = 1, jpk                           ! set the Kz to the background value
+              avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
+              avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
+              avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
+              avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
+           END DO
         ENDIF
-     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
-        CALL iom_rstput( kt, nitrst, numrow, 'en', en )
+        !
+     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+        !                                   ! -------------------
+        IF(lwp) WRITE(numout,*) '---- tke-rst ----'
+        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
+        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt    )
+        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm    )
+        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu   )
+        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv   )
+        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl  )
+        !
      ENDIF
      !
@@ -936,11 +815,15 @@
    !!   Dummy module :                                        NO TKE scheme
    !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke_old = .FALSE.   !: TKE flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
 CONTAINS
-   SUBROUTINE zdf_tke_old( kt )          ! Empty routine
-      WRITE(*,*) 'zdf_tke_old: You should not have seen this print! error?', kt
-   END SUBROUTINE zdf_tke_old
+   SUBROUTINE zdf_tke( kt )          ! Empty routine
+      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
+   END SUBROUTINE zdf_tke
+   SUBROUTINE tke_rst( kt, cdrw )
+     CHARACTER(len=*) ::   cdrw
+     WRITE(*,*) 'tke_rst: You should not have seen this print! error?', kt, cdwr
+   END SUBROUTINE tke_rst
 #endif
 
    !!======================================================================
-END MODULE zdftke_old
+END MODULE zdftke
Index: trunk/NEMO/OPA_SRC/ZDF/zdftke2.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdftke2.F90	(revision 1531)
+++ 	(revision )
@@ -1,829 +1,0 @@
-MODULE zdftke
-   !!======================================================================
-   !!                       ***  MODULE  zdftke  ***
-   !! Ocean physics:  vertical mixing coefficient computed from the tke 
-   !!                 turbulent closure parameterization
-   !!=====================================================================
-   !! History :  OPA  !  1991-03  (b. blanke)  Original code
-   !!            7.0  !  1991-11  (G. Madec)   bug fix
-   !!            7.1  !  1992-10  (G. Madec)   new mixing length and eav
-   !!            7.2  !  1993-03  (M. Guyon)   symetrical conditions
-   !!            7.3  !  1994-08  (G. Madec, M. Imbard)  nn_pdl flag
-   !!            7.5  !  1996-01  (G. Madec)   s-coordinates
-   !!            8.0  !  1997-07  (G. Madec)   lbc
-   !!            8.1  !  1999-01  (E. Stretta) new option for the mixing length
-   !!  NEMO      1.0  !  2002-06  (G. Madec) add tke_init routine
-   !!             -   !  2002-08  (G. Madec)  rn_cri and Free form, F90
-   !!             -   !  2004-10  (C. Ethe )  1D configuration
-   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
-   !!            3.0  !  2008-05  (C. Ethe,  G.Madec) : update TKE physics:
-   !!                 !           - tke penetration (wind steering)
-   !!                 !           - suface condition for tke & mixing length
-   !!                 !           - Langmuir cells
-   !!             -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
-   !!             -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
-   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term 
-   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl 
-   !!                 !                                + cleaning of the parameters + bugs correction
-   !!----------------------------------------------------------------------
-#if defined key_zdftke   ||   defined key_esopa
-   !!----------------------------------------------------------------------
-   !!   'key_zdftke'                                   TKE vertical physics
-   !!----------------------------------------------------------------------
-   !!   zdf_tke       : update momentum and tracer Kz from a tke scheme
-   !!   tke_tke       : tke time stepping: update tke at now time step (en)
-   !!   tke_avn       : compute mixing length scale and deduce avm and avt
-   !!   tke_init      : initialization, namelist read, and parameters control
-   !!   tke_rst       : read/write tke restart in ocean restart file
-   !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and active tracers 
-   USE dom_oce        ! ocean space and time domain
-   USE domvvl         ! ocean space and time domain : variable volume layer
-   USE zdf_oce        ! ocean vertical physics
-   USE sbc_oce        ! surface boundary condition: ocean
-   USE phycst         ! physical constants
-   USE zdfmxl         ! mixed layer
-   USE restart        ! only for lrst_oce
-   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE prtctl         ! Print control
-   USE in_out_manager ! I/O manager
-   USE iom            ! I/O manager library
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   zdf_tke    ! routine called in step module
-   PUBLIC   tke_rst    ! routine called in step module
-
-   LOGICAL , PUBLIC, PARAMETER              ::   lk_zdftke = .TRUE.  !: TKE vertical mixing flag
-
-#if defined key_c1d
-   !                                                                !!* 1D cfg only *
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_dis, e_mix        !: dissipation and mixing turbulent lengh scales
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_pdl, e_ric        !: prandl and local Richardson numbers
-#endif
-
-!!gm  - variables to be suppressed from the namelist:
-   LOGICAL  ::   ln_rstke = .FALSE.         ! =T restart with tke from a run without tke
-   INTEGER  ::   nn_itke  = 50              ! number of restart iterative loops
-   INTEGER  ::   nn_ave   =  1              ! horizontal average or not on avt, avmu, avmv (=0/1)
-   REAL(wp) ::   rn_cri   = 2._wp / 9._wp   ! critic Richardson number
-   REAL(wp) ::   rn_efave = 1._wp           ! coefficient for ave : ave=rn_efave*avm
-!!gm end
-
-!!gm  - variables to be added in the namelist (with a larger default value)
-   REAL(wp) ::   rn_bshear = 1.e-20   ! backgrounf shear (>0)
-!!gm
-
-   !                                       !!! ** Namelist  namtke  **
-   LOGICAL  ::   ln_mxl0  = .FALSE.         ! mixing length scale surface value as function of wind stress or not
-   INTEGER  ::   nn_mxl   =  2              ! type of mixing length (=0/1/2/3)
-!!gm to be change: the ref value of lmin0 and lmin
-   REAL(wp) ::   rn_lmin0 = 0.4_wp          ! surface  min value of mixing length   [m]
-   REAL(wp) ::   rn_lmin  = 0.1_wp          ! interior min value of mixing length   [m]
-   INTEGER  ::   nn_pdl   =  1              ! Prandtl number or not (ratio avt/avm) (=0/1)
-   INTEGER  ::   nn_avb   =  0              ! constant or profile background on avt (=0/1)
-   REAL(wp) ::   rn_ediff = 0.1_wp          ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
-   REAL(wp) ::   rn_ediss = 0.7_wp          ! coefficient of the Kolmogoroff dissipation 
-   REAL(wp) ::   rn_ebb   = 3.75_wp         ! coefficient of the surface input of tke
-   REAL(wp) ::   rn_emin  = 0.7071e-6_wp    ! minimum value of tke           [m2/s2]
-   REAL(wp) ::   rn_emin0 = 1.e-4_wp        ! surface minimum value of tke   [m2/s2]
-   INTEGER  ::   nn_havtb = 1               ! horizontal shape or not for avtb (=0/1)
-   INTEGER  ::   nn_etau  = 0               ! type of depth penetration of surface tke (=0/1/2)
-   INTEGER  ::   nn_htau  = 0               ! type of tke profile of penetration (=0/1/2)
-   REAL(wp) ::   rn_efr   = 1.0_wp          ! fraction of TKE surface value which penetrates in the ocean
-   LOGICAL  ::   ln_lc    = .FALSE.         ! Langmuir cells (LC) as a source term of TKE or not
-   REAL(wp) ::   rn_lc    = 0.15_wp         ! coef to compute vertical velocity of Langmuir cells
-
-   REAL(wp) ::   ri_cri   ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
-
-   REAL(wp), DIMENSION(jpi,jpj)     ::   htau      ! depth of tke penetration (nn_htau)
-   REAL(wp), DIMENSION(jpi,jpj)     ::   avtb_2d   ! set in tke_init, for other modif than ice
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   en        ! now turbulent kinetic energy   [m2/s2]
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   dissl     ! now mixing lenght of dissipation
-
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
-   !! $Id: zdftke2.F90 1201 2008-09-24 13:24:21Z rblod $
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE zdf_tke( kt )
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE zdf_tke  ***
-      !!
-      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
-      !!              coefficients using a turbulent closure scheme (TKE).
-      !!
-      !! ** Method  :   The time evolution of the turbulent kinetic energy (tke)
-      !!              is computed from a prognostic equation :
-      !!         d(en)/dt = avm (d(u)/dz)**2             ! shear production
-      !!                  + d( avm d(en)/dz )/dz         ! diffusion of tke
-      !!                  + avt N^2                      ! stratif. destruc.
-      !!                  - rn_ediss / emxl en**(2/3)    ! Kolmogoroff dissipation
-      !!      with the boundary conditions:
-      !!         surface: en = max( rn_emin0, rn_ebb sqrt(utau^2 + vtau^2) )
-      !!         bottom : en = rn_emin
-      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 
-      !!
-      !!        The now Turbulent kinetic energy is computed using the following 
-      !!      time stepping: implicit for vertical diffusion term, linearized semi
-      !!      implicit for kolmogoroff dissipation term, and explicit forward for 
-      !!      both buoyancy and shear production terms. Therefore a tridiagonal 
-      !!      linear system is solved. Note that buoyancy and shear terms are
-      !!      discretized in a energy conserving form (Bruchard 2002).
-      !!
-      !!        The dissipative and mixing length scale are computed from en and
-      !!      the stratification (see tke_avn)
-      !!
-      !!        The now vertical eddy vicosity and diffusivity coefficients are
-      !!      given by: 
-      !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
-      !!              avt = max( avmb, pdl * avm                 )  
-      !!              eav = max( avmb, avm )
-      !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
-      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1 
-      !!
-      !! ** Action  :   compute en (now turbulent kinetic energy)
-      !!                update avt, avmu, avmv (before vertical eddy coef.)
-      !!
-      !! References : Gaspar et al., JGR, 1990,
-      !!              Blanke and Delecluse, JPO, 1991
-      !!              Mellor and Blumberg, JPO 2004
-      !!              Axell, JGR, 2002
-      !!              Bruchard OM 2002
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time step
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == nit000 )   CALL tke_init     ! initialisation 
-                           !
-                           CALL tke_tke      ! now tke (en)
-                           !
-                           CALL tke_avn      ! now avt, avm, avmu, avmv
-      !
-   END SUBROUTINE zdf_tke
-
-
-   SUBROUTINE tke_tke
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE tke_tke  ***
-      !!
-      !! ** Purpose :   Compute the now Turbulente Kinetic Energy (TKE)
-      !!
-      !! ** Method  : - TKE surface boundary condition
-      !!              - source term due to Langmuir cells (ln_lc=T)
-      !!              - source term due to shear (saved in avmu, avmv arrays)
-      !!              - Now TKE : resolution of the TKE equation by inverting
-      !!                a tridiagonal linear system by a "methode de chasse"
-      !!              - increase TKE due to surface and internal wave breaking
-      !!
-      !! ** Action  : - en : now turbulent kinetic energy)
-      !!              - avmu, avmv : production of TKE by shear at u and v-points
-      !!                (= Kz dz[Ub] * dz[Un] )
-      !! ---------------------------------------------------------------------
-      USE oce,   zdiag  =>   ua   ! use ua as workspace
-      USE oce,   zd_up  =>   va   ! use va as workspace
-      USE oce,   zd_lw  =>   ta   ! use ta as workspace
-      !!
-      INTEGER  ::   ji, jj, jk                ! dummy loop arguments
-      REAL(wp) ::   zbbrau, zesurf, zesh2     ! temporary scalars
-      REAL(wp) ::   zfact1, zfact2, zfact3    !    -         -
-      REAL(wp) ::   ztx2  , zty2  , zcof      !    -         -
-      REAL(wp) ::   zus   , zwlc  , zind      !    -         -
-      REAL(wp) ::   zzd_up, zzd_lw            !    -         -
-      INTEGER , DIMENSION(jpi,jpj)     ::   imlc    ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc    !  -      -
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc   ! 3D workspace
-      !!--------------------------------------------------------------------
-      !
-      zbbrau =  .5 * rn_ebb / rau0       ! Local constant initialisation
-      zfact1 = -.5 * rdt 
-      zfact2 = 1.5 * rdt * rn_ediss
-      zfact3 = 0.5       * rn_ediss
-      !
-      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !                     !  Surface boundary condition on tke
-      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-!CDIR NOVERRCHK
-      DO jj = 2, jpjm1         ! en(1)   = rn_ebb sqrt(utau^2+vtau^2) / rau0  (min value rn_emin0)
-!CDIR NOVERRCHK
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            ztx2   = utau(ji-1,jj  ) + utau(ji,jj)
-            zty2   = vtau(ji  ,jj-1) + vtau(ji,jj)
-            zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 )
-            en(ji,jj,1) = MAX( zesurf, rn_emin0 ) * tmask(ji,jj,1)
-         END DO
-      END DO
-      !
-      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !                     !  Bottom boundary condition on tke
-      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-!!gm to be added soon
-      !
-      !
-      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke
-         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-         !
-         !                        !* total energy produce by LC : cumulative sum over jk
-         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1)
-         DO jk = 2, jpk
-            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
-         END DO
-         !                        !* finite Langmuir Circulation depth
-         imlc(:,:) = mbathy(:,:)         ! Initialization to the number of w ocean point mbathy
-         DO jk = jpkm1, 2, -1
-            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us 
-               DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
-                  zus  = 0.000128 * wndm(ji,jj) * wndm(ji,jj)
-                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
-               END DO
-            END DO
-         END DO
-         !                               ! finite LC depth
-# if defined key_vectopt_loop
-         DO jj = 1, 1
-            DO ji = 1, jpij   ! vector opt. (forced unrolling)
-# else
-         DO jj = 1, jpj 
-            DO ji = 1, jpi
-# endif
-               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
-            END DO
-         END DO
-# if defined key_c1d
-         hlc(:,:) = zhlc(:,:) * tmask(:,:,1)      ! c1d configuration: save finite Langmuir Circulation depth
-# endif
-!CDIR NOVERRCHK
-         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en
-!CDIR NOVERRCHK
-            DO jj = 2, jpjm1
-!CDIR NOVERRCHK
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-!!gm replace here wndn by a formulation with the stress module
-                  zus  = 0.016 * wndm(ji,jj)                  ! Stokes drift
-                  !                                           ! vertical velocity due to LC
-                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
-                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
-                  !                                           ! TKE Langmuir circulation source term
-                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         !
-      ENDIF
-      !
-      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !                     !  Now Turbulent kinetic energy (output in en)
-      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
-      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
-      !                     ! 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, 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
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zcof   = zfact1 * tmask(ji,jj,jk)
-               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal
-                  &          / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk  ) )
-               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  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
-                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , 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.e0 - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(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)  ) * tmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
-      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
-               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
-            END DO
-         END DO
-      END DO
-      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
-         DO ji = fs_2, fs_jpim1    ! vector opt.
-            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
-         END DO
-      END DO
-      DO jk = 3, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
-               zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1)
-            END DO
-         END DO
-      END DO
-      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
-         DO ji = fs_2, fs_jpim1    ! vector opt.
-            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
-         END DO
-      END DO
-      DO jk = jpk-2, 2, -1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
-               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      DO jk = 2, jpkm1                             ! set the minimum value of tke
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-
-      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !                            !  TKE due to surface and internal wave breaking
-      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( nn_etau == 1 ) THEN           !*  penetration throughout the water column
-         DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
-                     &                                               * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               jk = nmln(ji,jj)
-               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
-                  &                                               * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
-            END DO
-         END DO
-      ENDIF
-      !
-      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
-      !
-   END SUBROUTINE tke_tke
-
-
-   SUBROUTINE tke_avn
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE tke_avn  ***
-      !!
-      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
-      !!
-      !! ** Method  :   At this stage, en, the now TKE, is known (computed in 
-      !!              the tke_tke routine). First, the now mixing lenth is 
-      !!      computed from en and the strafification (N^2), then the mixings
-      !!      coefficients are computed.
-      !!              - Mixing length : a first evaluation of the mixing lengh
-      !!      scales is:
-      !!                      mxl = sqrt(2*en) / N  
-      !!      where N is the brunt-vaisala frequency, with a minimum value set
-      !!      to rn_lmin (rn_lmin0) in the interior (surface) ocean.
-      !!        The mixing and dissipative length scale are bound as follow : 
-      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
-      !!                        zmxld = zmxlm = mxl
-      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
-      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 
-      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
-      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
-      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to 
-      !!                    the surface to obtain ldown. the resulting length 
-      !!                    scales are:
-      !!                        zmxld = sqrt( lup * ldown ) 
-      !!                        zmxlm = min ( lup , ldown )
-      !!              - Vertical eddy viscosity and diffusivity:
-      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
-      !!                      avt = max( avmb, pdlr * avm )  
-      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
-      !!
-      !! ** Action  : - avt : now vertical eddy diffusivity (w-point)
-      !!              - avmu, avmv : now vertical eddy viscosity at uw- and vw-points
-      !!----------------------------------------------------------------------
-      USE oce,     zmpdl  =>   ua   ! use ua as workspace
-      USE oce,     zmxlm  =>   va   ! use va as workspace
-      USE oce,     zmxld  =>   ta   ! use ta as workspace
-      !!
-      INTEGER  ::   ji, jj, jk            ! dummy loop arguments
-      REAL(wp) ::   zrn2, zraug           ! temporary scalars
-      REAL(wp) ::   ztx2, zdku            !    -         -
-      REAL(wp) ::   zty2, zdkv            !    -         -
-      REAL(wp) ::   zcoef, zav            !    -         -
-      REAL(wp) ::   zpdlr, zri, zsqen     !    -         -
-      REAL(wp) ::   zemxl, zemlm, zemlp   !    -         -
-      !!--------------------------------------------------------------------
-
-      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !                     !  Mixing length
-      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !
-      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
-      !
-      IF( ln_mxl0 ) THEN         ! surface mixing length = F(stress) : l=vkarmn*2.e5*sqrt(utau^2 + vtau^2)/(rau0*g)
-!!gm  this should be useless
-         zmxlm(:,:,1) = 0.e0
-!!gm end
-         zraug = 0.5 * vkarmn * 2.e5 / ( rau0 * grav )
-         DO jj = 2, jpjm1
-!CDIR NOVERRCHK
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
-               zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
-               zmxlm(ji,jj,1) = MAX(  rn_lmin0,  zraug * SQRT( ztx2 * ztx2 + zty2 * zty2 )  )
-            END DO
-         END DO
-      ELSE                       ! surface set to the minimum value
-         zmxlm(:,:,1) = rn_lmin0
-      ENDIF
-      zmxlm(:,:,jpk) = rn_lmin   ! bottom set to the interior minium value
-      !
-!CDIR NOVERRCHK
-      DO jk = 2, jpkm1           ! interior value : l=sqrt(2*e/n**2)
-!CDIR NOVERRCHK
-         DO jj = 2, jpjm1
-!CDIR NOVERRCHK
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
-               zmxlm(ji,jj,jk) = MAX(  rn_lmin,  SQRT( 2. * en(ji,jj,jk) / zrn2 )  )
-            END DO
-         END DO
-      END DO
-      !
-!!gm  CAUTION: to be added here the bottom boundary condition on zmxlm
-      !
-      !                     !* Physical limits for the mixing length
-      !
-      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value
-      zmxld(:,:,jpk) = rn_lmin        ! bottom  set to the minimum value
-      !
-      SELECT CASE ( nn_mxl )
-      !
-      CASE ( 0 )           ! bounded by the distance to surface and bottom
-         DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk),   &
-                  &            fsdepw(ji,jj,mbathy(ji,jj)) - fsdepw(ji,jj,jk) )
-                  zmxlm(ji,jj,jk) = zemxl
-                  zmxld(ji,jj,jk) = zemxl
-               END DO
-            END DO
-         END DO
-         !
-      CASE ( 1 )           ! bounded by the vertical scale factor
-         DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
-                  zmxlm(ji,jj,jk) = zemxl
-                  zmxld(ji,jj,jk) = zemxl
-               END DO
-            END DO
-         END DO
-         !
-      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
-         DO jk = 2, jpkm1         ! from the surface to the bottom :
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
-               END DO
-            END DO
-         END DO
-         DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
-                  zmxlm(ji,jj,jk) = zemxl
-                  zmxld(ji,jj,jk) = zemxl
-               END DO
-            END DO
-         END DO
-         !
-      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
-         DO jk = 2, jpkm1         ! from the surface to the bottom : lup
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
-               END DO
-            END DO
-         END DO
-         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
-               END DO
-            END DO
-         END DO
-!CDIR NOVERRCHK
-         DO jk = 2, jpkm1
-!CDIR NOVERRCHK
-            DO jj = 2, jpjm1
-!CDIR NOVERRCHK
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
-                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
-                  zmxlm(ji,jj,jk) = zemlm
-                  zmxld(ji,jj,jk) = zemlp
-               END DO
-            END DO
-         END DO
-         !
-      END SELECT
-      !
-# if defined key_c1d
-      e_dis(:,:,:) = zmxld(:,:,:)      ! c1d configuration : save mixing and dissipation turbulent length scales
-      e_mix(:,:,:) = zmxlm(:,:,:)
-# endif
-
-      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt)
-      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-!CDIR NOVERRCHK
-      DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
-!CDIR NOVERRCHK
-         DO jj = 2, jpjm1
-!CDIR NOVERRCHK
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zsqen = SQRT( en(ji,jj,jk) )
-               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
-               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * tmask(ji,jj,jk)
-               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
-               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
-      !
-      DO jk = 2, jpkm1            !* vertical eddy viscosity at u- and v-points
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
-               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
-      !
-      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
-         DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zcoef = 0.5 / ( 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. ) * avm(ji,jj,jk) / (zcoef * (zdku + zdkv + rn_bshear ) )
-                  zpdlr = MAX(  0.1,  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,  ri_crit / MAX( ri_crit , zri )  )
-                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
-# if defined key_c1d
-                  e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number
-                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri
-# endif
-              END DO
-            END DO
-         END DO
-      ENDIF
-      CALL lbc_lnk( avt, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged)
-
-      IF(ln_ctl) THEN
-         CALL prt_ctl( tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
-         CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask,                   &
-            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
-      ENDIF
-      !
-   END SUBROUTINE tke_avn
-
-
-   SUBROUTINE tke_init
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE tke_init  ***
-      !!                     
-      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
-      !!              viscosity when using a tke turbulent closure scheme
-      !!
-      !! ** Method  :   Read the namtke namelist and check the parameters
-      !!              called at the first timestep (nit000)
-      !!
-      !! ** input   :   Namlist namtke
-      !!
-      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
-      !!----------------------------------------------------------------------
-      INTEGER ::   ji, jj, jk   ! dummy loop indices
-      !!
-      NAMELIST/namtke/ ln_rstke, rn_ediff, rn_ediss, rn_ebb  , rn_efave, rn_emin,   &
-         &             rn_emin0, rn_cri  , nn_itke , nn_mxl  , nn_pdl  , nn_ave ,   &
-         &             nn_avb  , ln_mxl0 , rn_lmin , rn_lmin0, nn_havtb, nn_etau,   &
-         &             nn_htau , rn_efr  , ln_lc   , rn_lc 
-      !!----------------------------------------------------------------------
-
-      REWIND ( numnam )               !* Read Namelist namtke : Turbulente Kinetic Energy
-      READ   ( numnam, namtke )
-      
-      ri_cri = 2. / ( 2. + rn_ediss / rn_ediff )      ! resulting critical Richardson number
-
-      IF(lwp) THEN                    !* Control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'zdf_tke : tke turbulent closure scheme - initialisation'
-         WRITE(numout,*) '~~~~~~~~'
-         WRITE(numout,*) '          Namelist namtke : set tke mixing parameters'
-         WRITE(numout,*) '             restart with tke from no tke              ln_rstke = ', ln_rstke
-         WRITE(numout,*) '             coef. to compute avt                      rn_ediff = ', rn_ediff
-         WRITE(numout,*) '             Kolmogoroff dissipation coef.             rn_ediss = ', rn_ediss
-         WRITE(numout,*) '             tke surface input coef.                   rn_ebb   = ', rn_ebb
-         WRITE(numout,*) '             minimum value of tke                      rn_emin  = ', rn_emin
-         WRITE(numout,*) '             surface minimum value of tke              rn_emin0 = ', rn_emin0
-         WRITE(numout,*) '             mixing length type                        nn_mxl   = ', nn_mxl
-         WRITE(numout,*) '             prandl number flag                        nn_pdl   = ', nn_pdl
-         WRITE(numout,*) '             constant background or profile            nn_avb   = ', nn_avb
-         WRITE(numout,*) '             surface mixing length = F(stress) or not  ln_mxl0  = ', ln_mxl0
-         WRITE(numout,*) '             surface  mixing length minimum value      rn_lmin0 = ', rn_lmin0
-         WRITE(numout,*) '             interior mixing length minimum value      rn_lmin0 = ', rn_lmin0
-         WRITE(numout,*) '             horizontal variation for avtb             nn_havtb = ', nn_havtb
-         WRITE(numout,*) '             test param. to add tke induced by wind    nn_etau  = ', nn_etau
-         WRITE(numout,*) '             flag for computation of exp. tke profile  nn_htau  = ', nn_htau
-         WRITE(numout,*) '             % of rn_emin0 which pene. the thermocline rn_efr   = ', rn_efr
-         WRITE(numout,*) '             flag to take into acc.  Langmuir circ.    ln_lc    = ', ln_lc
-         WRITE(numout,*) '             coef to compute verticla velocity of LC   rn_lc    = ', rn_lc
-         WRITE(numout,*)
-         WRITE(numout,*) '             critical Richardson nb with your choice of coefs.  = ', ri_cri
-      ENDIF
-
-      !                               !* Check of some namelist values
-      IF( nn_mxl  < 0    .OR. nn_mxl  > 3   )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
-      IF( nn_pdl  < 0    .OR. nn_pdl  > 1   )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
-      IF( nn_htau < 0    .OR. nn_htau > 2   )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
-      IF( rn_lc   < 0.15 .OR. rn_lc   > 0.2 )   CALL ctl_stop( 'bad value: rn_lc must be between 0.15 and 0.2 ' )
-
-      IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln 
-
-      !                               !* Background eddy viscosity and diffusivity profil
-      IF( nn_avb == 0 ) THEN                ! Define avmb, avtb from namelist parameter
-         avmb(:) = avm0
-         avtb(:) = avt0
-      ELSE                                  ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990) 
-         avmb(:) = avm0
-         avtb(:) = avt0 + ( 3.0e-4 - 2 * avt0 ) * 1.0e-4 * gdepw_0(:)   ! m2/s
-         IF(ln_sco .AND. lwp)   CALL ctl_warn( '          avtb profile not valid in sco' )
-      ENDIF
-      !
-      !                                     ! 2D shape of the avtb
-      avtb_2d(:,:) = 1.e0                        ! uniform 
-      !
-      IF( nn_havtb == 1 ) THEN                   ! decrease avtb in the equatorial band
-           !  -15S -5S : linear decrease from avt0 to avt0/10.
-           !  -5S  +5N : cst value avt0/10.
-           !   5N  15N : linear increase from avt0/10, to avt0
-           WHERE(-15. <= gphit .AND. gphit < -5 )   avtb_2d = (1.  - 0.09 * (gphit + 15.))
-           WHERE( -5. <= gphit .AND. gphit <  5 )   avtb_2d =  0.1
-           WHERE(  5. <= gphit .AND. gphit < 15 )   avtb_2d = (0.1 + 0.09 * (gphit -  5.))
-      ENDIF
-
-      !                               !* depth of penetration of surface tke
-      IF( nn_etau /= 0 ) THEN      
-         SELECT CASE( nn_htau )           ! Choice of the depth of penetration
-         CASE( 0 )                                    ! constant depth penetration (here 10 meters)
-            htau(:,:) = 10.e0
-         CASE( 1 )                                    ! F(latitude) : 5m to 40m at high lat.
-            DO jj = 1, jpj
-               DO ji = 1, jpi  
-                  htau(ji,jj) = MAX( 5., MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
-               END DO
-            END DO
-         CASE( 2 )                                    ! F(latitude) : 5m to 60m at high lat.
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  htau(ji,jj) =  MAX( 5.,6./4.* MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
-               END DO
-            END DO
-         END SELECT
-      ENDIF
-
-      !                               !* set vertical eddy coef. to the background value
-      DO jk = 1, jpk
-         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
-         avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
-         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
-         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
-      END DO
-      dissl(:,:,:) = 1.e-12
-      !                               !* read or initialize all required files 
-      CALL tke_rst( nit000, 'READ' )
-      !
-   END SUBROUTINE tke_init
-
-
-   SUBROUTINE tke_rst( kt, cdrw )
-     !!---------------------------------------------------------------------
-     !!                   ***  ROUTINE tke_rst  ***
-     !!                     
-     !! ** Purpose :   Read or write TKE file (en) in restart file
-     !!
-     !! ** Method  :   use of IOM library
-     !!                if the restart does not contain TKE, en is either 
-     !!                set to rn_emin or recomputed (nn_itke/=0)
-     !!----------------------------------------------------------------------
-     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
-     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
-     !
-     INTEGER ::   jit, jk   ! dummy loop indices
-     INTEGER ::   id1, id2, id3, id4, id5, id6
-     !!----------------------------------------------------------------------
-     !
-     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
-        !                                   ! ---------------
-        IF( ln_rstart ) THEN                   !* Read the restart file
-           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
-           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
-           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
-           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
-           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
-           id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
-           !
-           IF( id1 > 0 ) THEN                       ! 'en' exists
-              CALL iom_get( numror, jpdom_autoglo, 'en', en )
-              IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
-                 CALL iom_get( numror, jpdom_autoglo, 'avt'  , avt   )
-                 CALL iom_get( numror, jpdom_autoglo, 'avm'  , avm   )
-                 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu  )
-                 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv  )
-                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
-              ELSE                                                 ! one at least array is missing
-                 CALL tke_avn                                          ! compute avt, avm, avmu, avmv and dissl (approximation)
-              ENDIF
-           ELSE                                     ! No TKE array found: initialisation
-              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
-              en (:,:,:) = rn_emin * tmask(:,:,:)
-              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
-              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
-           ENDIF
-        ELSE                                   !* Start from rest
-           en(:,:,:) = rn_emin * tmask(:,:,:)
-           DO jk = 1, jpk                           ! set the Kz to the background value
-              avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
-              avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
-              avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
-              avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
-           END DO
-        ENDIF
-        !
-     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
-        !                                   ! -------------------
-        IF(lwp) WRITE(numout,*) '---- tke-rst ----'
-        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
-        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt    )
-        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm    )
-        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu   )
-        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv   )
-        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl  )
-        !
-     ENDIF
-     !
-   END SUBROUTINE tke_rst
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Dummy module :                                        NO TKE scheme
-   !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
-CONTAINS
-   SUBROUTINE zdf_tke( kt )          ! Empty routine
-      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
-   END SUBROUTINE zdf_tke
-   SUBROUTINE tke_rst( kt, cdrw )
-     CHARACTER(len=*) ::   cdrw
-     WRITE(*,*) 'tke_rst: You should not have seen this print! error?', kt, cdwr
-   END SUBROUTINE tke_rst
-#endif
-
-   !!======================================================================
-END MODULE zdftke
Index: trunk/NEMO/OPA_SRC/ZDF/zdftke_old.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdftke_old.F90	(revision 1532)
+++ trunk/NEMO/OPA_SRC/ZDF/zdftke_old.F90	(revision 1532)
@@ -0,0 +1,946 @@
+MODULE zdftke_old
+   !!======================================================================
+   !!                       ***  MODULE  zdftke_old  ***
+   !! Ocean physics:  vertical mixing coefficient computed from the tke 
+   !!                 turbulent closure parameterization
+   !!=====================================================================
+   !! History :   OPA  !  1991-03  (b. blanke)  Original code
+   !!             7.0  !  1991-11  (G. Madec)   bug fix
+   !!             7.1  !  1992-10  (G. Madec)   new mixing length and eav
+   !!             7.2  !  1993-03  (M. Guyon)   symetrical conditions
+   !!             7.3  !  1994-08  (G. Madec, M. Imbard)  nn_pdl flag
+   !!             7.5  !  1996-01  (G. Madec)   s-coordinates
+   !!             8.0  !  1997-07  (G. Madec)   lbc
+   !!             8.1  !  1999-01  (E. Stretta) new option for the mixing length
+   !!  NEMO       1.0  !  2002-06  (G. Madec) add zdf_tke_init routine
+   !!              -   !  2002-08  (G. Madec)  rn_cri and Free form, F90
+   !!              -   !  2004-10  (C. Ethe )  1D configuration
+   !!             2.0  !  2006-07  (S. Masson)  distributed restart using iom
+   !!             3.0  !  2008-05  (C. Ethe,  G.Madec) : update TKE physics:
+   !!                              - tke penetration (wind steering)
+   !!                              - suface condition for tke & mixing length
+   !!                              - Langmuir cells
+   !!              -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
+   !!              -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
+   !!----------------------------------------------------------------------
+#if defined key_zdftke_old   ||   defined key_esopa
+   !!----------------------------------------------------------------------
+   !!   'key_zdftke_old'                               TKE vertical physics
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   zdf_tke_old  : update momentum and tracer Kz from a tke scheme
+   !!   zdf_tke_init : initialization, namelist read, and parameters control
+   !!   tke_rst      : read/write tke restart in ocean restart file
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and active tracers 
+   USE dom_oce         ! ocean space and time domain
+   USE zdf_oce         ! ocean vertical physics
+   USE sbc_oce         ! surface boundary condition: ocean
+   USE phycst          ! physical constants
+   USE zdfmxl          ! mixed layer
+   USE restart         ! only for lrst_oce
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl          ! Print control
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O manager library
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   zdf_tke_old   ! routine called in step module
+
+   LOGICAL , PUBLIC, PARAMETER              ::   lk_zdftke_old = .TRUE.  !: TKE vertical mixing flag
+   REAL(wp), PUBLIC                         ::   eboost              !: multiplicative coeff of the shear product.
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   en                  !: now turbulent kinetic energy
+# if defined key_vectopt_memory
+   !                                                                !!! key_vectopt_memory
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   etmean              !: coefficient used for horizontal smoothing
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   eumean, evmean      !: at t-, u- and v-points
+# endif
+#if defined key_c1d
+   !                                                                !!! 1D cfg only
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_dis, e_mix        !: dissipation and mixing turbulent lengh scales
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_pdl, e_ric        !: prandl and local Richardson numbers
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   hlc                 !: save finite Langmuir Circulation depth
+#endif
+
+   !                                       !!! ** Namelist  namtke  **
+   LOGICAL  ::   ln_rstke = .FALSE.         ! =T restart with tke from a run without tke
+   LOGICAL  ::   ln_mxl0  = .FALSE.         ! mixing length scale surface value as function of wind stress or not
+   LOGICAL  ::   ln_lc    = .FALSE.         ! Langmuir cells (LC) as a source term of TKE or not
+   INTEGER  ::   nn_itke  = 50              ! number of restart iterative loops
+   INTEGER  ::   nn_mxl   =  2              ! type of mixing length (=0/1/2/3)
+   INTEGER  ::   nn_pdl   =  1              ! Prandtl number or not (ratio avt/avm) (=0/1)
+   INTEGER  ::   nn_ave   =  1              ! horizontal average or not on avt, avmu, avmv (=0/1)
+   INTEGER  ::   nn_avb   =  0              ! constant or profile background on avt (=0/1)
+   REAL(wp) ::   rn_ediff = 0.1_wp          ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
+   REAL(wp) ::   rn_ediss = 0.7_wp          ! coefficient of the Kolmogoroff dissipation 
+   REAL(wp) ::   rn_ebb   = 3.75_wp         ! coefficient of the surface input of tke
+   REAL(wp) ::   rn_efave = 1._wp           ! coefficient for ave : ave=rn_efave*avm
+   REAL(wp) ::   rn_emin  = 0.7071e-6_wp    ! minimum value of tke (m2/s2)
+   REAL(wp) ::   rn_emin0 = 1.e-4_wp        ! surface minimum value of tke (m2/s2)
+   REAL(wp) ::   rn_cri   = 2._wp / 9._wp   ! critic Richardson number
+   INTEGER  ::   nn_havtb = 1               ! horizontal shape or not for avtb (=0/1)
+   INTEGER  ::   nn_etau  = 0               ! type of depth penetration of surface tke (=0/1/2)
+   INTEGER  ::   nn_htau  = 0               ! type of tke profile of penetration (=0/1/2)
+   REAL(wp) ::   rn_lmin0 = 0.4_wp          ! surface  min value of mixing length
+   REAL(wp) ::   rn_lmin  = 0.1_wp          ! interior min value of mixing length
+   REAL(wp) ::   rn_efr   = 1.0_wp          ! fraction of TKE surface value which penetrates in the ocean
+   REAL(wp) ::   rn_lc    = 0.15_wp         ! coef to compute vertical velocity of Langmuir cells
+
+   REAL(wp), DIMENSION (jpi,jpj) ::   avtb_2d   ! set in tke_init, for other modif than ice
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE zdf_tke_old( kt )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE zdf_tke_old  ***
+      !!
+      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
+      !!      coefficients using a 1.5 turbulent closure scheme.
+      !!
+      !! ** Method  :   The time evolution of the turbulent kinetic energy 
+      !!      (tke) is computed from a prognostic equation :
+      !!         d(en)/dt = eboost eav (d(u)/dz)**2         ! shear production
+      !!                  + d( rn_efave eav d(en)/dz )/dz   ! diffusion of tke
+      !!                  + grav/rau0 pdl eav d(rau)/dz     ! stratif. destruc.
+      !!                  - rn_ediss / emxl en**(2/3)       ! dissipation
+      !!      with the boundary conditions:
+      !!         surface: en = max( rn_emin0, rn_ebb sqrt(utau^2 + vtau^2) )
+      !!         bottom : en = rn_emin
+      !!      -1- The dissipation and mixing turbulent lengh scales are computed
+      !!         from the usual diagnostic buoyancy length scale:  
+      !!         mxl= sqrt(2*en)/N  where N is the brunt-vaisala frequency
+      !!         with mxl = rn_lmin at the bottom minimum value of 0.4
+      !!      Four cases : 
+      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
+      !!                  zmxld = zmxlm = mxl
+      !!         nn_mxl=1 : mxl bounded by the vertical scale factor.
+      !!                  zmxld = zmxlm = mxl
+      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl 
+      !!                  is less than 1 (|d/dz(xml)|<1). 
+      !!                  zmxld = zmxlm = mxl
+      !!         nn_mxl=3 : lup = mxl bounded using |d/dz(xml)|<1 from the surface 
+      !!                        to the bottom
+      !!                  ldown = mxl bounded using |d/dz(xml)|<1 from the bottom 
+      !!                        to the surface
+      !!                  zmxld = sqrt (lup*ldown) ; zmxlm = min(lup,ldown)
+      !!      -2- Compute the now Turbulent kinetic energy. The time differencing
+      !!      is implicit for vertical diffusion term, linearized for kolmo-
+      !!      goroff dissipation term, and explicit forward for both buoyancy
+      !!      and dynamic production terms. Thus a tridiagonal linear system is
+      !!      solved.
+      !!         Note that - the shear production is multiplied by eboost in order
+      !!      to set the critic richardson number to rn_cri (namelist parameter)
+      !!                   - the destruction by stratification term is multiplied
+      !!      by the Prandtl number (defined by an empirical funtion of the local 
+      !!      Richardson number) if nn_pdl=1 (namelist parameter)
+      !!      coefficient (zesh2):
+      !!      -3- Compute the now vertical eddy vicosity and diffusivity
+      !!      coefficients from en (before the time stepping) and zmxlm:
+      !!              avm = max( avtb, rn_ediff*zmxlm*en^1/2 )
+      !!              avt = max( avmb, pdl*avm )  (pdl=1 if nn_pdl=0)
+      !!              eav = max( avmb, avm )
+      !!      avt and avm are horizontally averaged to avoid numerical insta-
+      !!      bilities.
+      !!        N.B. The computation is done from jk=2 to jpkm1 except for
+      !!      en. Surface value of avt avmu avmv are set once a time to
+      !!      their background value in routine zdf_tke_init.
+      !!
+      !! ** Action  :   compute en (now turbulent kinetic energy)
+      !!                update avt, avmu, avmv (before vertical eddy coef.)
+      !!
+      !! References : Gaspar et al., JGR, 1990,
+      !!              Blanke and Delecluse, JPO, 1991
+      !!              Mellor and Blumberg, JPO 2004
+      !!              Axell, JGR, 2002
+      !!----------------------------------------------------------------------
+      USE oce,     zwd    =>   ua   ! use ua as workspace
+      USE oce,     zmxlm  =>   va   ! use va as workspace
+      USE oce,     zmxld  =>   ta   ! use ta as workspace
+      USE oce,     ztkelc =>   sa   ! use sa as workspace
+      !
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      !
+      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
+      REAL(wp) ::   zbbrau, zrn2, zesurf            ! temporary scalars
+      REAL(wp) ::   zfact1, ztx2, zdku              !    -         -
+      REAL(wp) ::   zfact2, zty2, zdkv              !    -         -
+      REAL(wp) ::   zfact3, zcoef, zcof, zav        !    -         -
+      REAL(wp) ::   zsh2, zpdl, zri, zsqen, zesh2   !    -         -
+      REAL(wp) ::   zemxl, zemlm, zemlp             !    -         -
+      REAL(wp) ::   zraug, zus, zwlc, zind          !    -         -
+      INTEGER , DIMENSION(jpi,jpj)     ::   imlc    ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zhtau   !  -      -
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc    !  -      -
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc   ! 3D workspace
+      !!--------------------------------------------------------------------
+
+      IF( kt == nit000  )   CALL zdf_tke_init      ! Initialization (first time-step only)
+
+      !                                            ! Local constant initialization
+      zbbrau =  .5 * rn_ebb / rau0
+      zfact1 = -.5 * rdt * rn_efave
+      zfact2 = 1.5 * rdt * rn_ediss
+      zfact3 = 0.5 * rdt * rn_ediss
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! I.  Mixing length
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+      ! Buoyancy length scale: l=sqrt(2*e/n**2)
+      ! ---------------------
+      IF( ln_mxl0 ) THEN         ! surface mixing length = F(stress) : l=vkarmn*2.e5*sqrt(utau^2 + vtau^2)/(rau0*g)
+!!gm  this should be useless
+         zmxlm(:,:,1) = 0.e0
+!!gm end
+         zraug = 0.5 * vkarmn * 2.e5 / ( rau0 * grav )
+         DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
+               zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
+               zmxlm(ji,jj,1) = MAX(  rn_lmin0,  zraug * SQRT( ztx2 * ztx2 + zty2 * zty2 )  )
+            END DO
+         END DO
+      ELSE                       ! surface set to the minimum value
+         zmxlm(:,:,1) = rn_lmin0
+      ENDIF
+      zmxlm(:,:,jpk) = rn_lmin   ! bottom set to the interior minium value
+      !
+!CDIR NOVERRCHK
+      DO jk = 2, jpkm1           ! interior value : l=sqrt(2*e/n**2)
+!CDIR NOVERRCHK
+         DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zrn2 = MAX( rn2b(ji,jj,jk), rsmall )
+               zmxlm(ji,jj,jk) = MAX(  rn_lmin,  SQRT( 2. * en(ji,jj,jk) / zrn2 )  )
+            END DO
+         END DO
+      END DO
+
+      ! Physical limits for the mixing length
+      ! -------------------------------------
+      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value
+      zmxld(:,:,jpk) = rn_lmin        ! bottom  set to the minimum value
+
+      SELECT CASE ( nn_mxl )
+      !
+      CASE ( 0 )           ! bounded by the distance to surface and bottom
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk),   &
+                  &            fsdepw(ji,jj,mbathy(ji,jj)) - fsdepw(ji,jj,jk) )
+                  zmxlm(ji,jj,jk) = zemxl
+                  zmxld(ji,jj,jk) = zemxl
+               END DO
+            END DO
+         END DO
+         !
+      CASE ( 1 )           ! bounded by the vertical scale factor
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
+                  zmxlm(ji,jj,jk) = zemxl
+                  zmxld(ji,jj,jk) = zemxl
+               END DO
+            END DO
+         END DO
+         !
+      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
+         DO jk = 2, jpkm1         ! from the surface to the bottom :
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+         DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
+                  zmxlm(ji,jj,jk) = zemxl
+                  zmxld(ji,jj,jk) = zemxl
+               END DO
+            END DO
+         END DO
+         !
+      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
+         DO jk = 2, jpkm1         ! from the surface to the bottom : lup
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+!CDIR NOVERRCHK
+         DO jk = 2, jpkm1
+!CDIR NOVERRCHK
+            DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
+                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
+                  zmxlm(ji,jj,jk) = zemlm
+                  zmxld(ji,jj,jk) = zemlp
+               END DO
+            END DO
+         END DO
+         !
+      END SELECT
+
+# if defined key_c1d
+      ! c1d configuration : save mixing and dissipation turbulent length scales
+      e_dis(:,:,:) = zmxld(:,:,:)
+      e_mix(:,:,:) = zmxlm(:,:,:)
+# endif
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! II  TKE Langmuir circulation source term
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_lc ) THEN
+         !
+         ! Computation of total energy produce by LC : cumulative sum over jk
+         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1)
+         DO jk = 2, jpk
+            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
+         END DO
+         !
+         ! Computation of finite Langmuir Circulation depth
+         ! Initialization to the number of w ocean point mbathy
+         imlc(:,:) = mbathy(:,:)
+         DO jk = jpkm1, 2, -1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ! Last w-level at which zpelc>=0.5*us*us with us=0.016*wind(starting from jpk-1)
+                  zus  = 0.000128 * wndm(ji,jj) * wndm(ji,jj)
+                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
+               END DO
+            END DO
+         END DO
+         !
+         ! finite LC depth
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
+            END DO
+         END DO
+         !
+# if defined key_c1d
+         hlc(:,:) = zhlc(:,:) * tmask(:,:,1)      ! c1d configuration: save finite Langmuir Circulation depth
+# endif
+         !
+         ! TKE Langmuir circulation source term
+!CDIR NOVERRCHK
+         DO jk = 2, jpkm1
+!CDIR NOVERRCHK
+            DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ! Stokes drift
+                  zus  = 0.016 * wndm(ji,jj)
+                  ! computation of vertical velocity due to LC
+                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
+                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
+                  ! TKE Langmuir circulation source term
+                  ztkelc(ji,jj,jk) = ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)
+               END DO
+            END DO
+         END DO
+         !
+      ENDIF
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! III  Tubulent kinetic energy time stepping
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+      ! 1. Vertical eddy viscosity on tke (put in zmxlm) and first estimate of avt
+      ! ---------------------------------------------------------------------
+!CDIR NOVERRCHK
+      DO jk = 2, jpkm1
+!CDIR NOVERRCHK
+         DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zsqen = SQRT( en(ji,jj,jk) )
+               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
+               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
+               zmxlm(ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)
+               zmxld(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      ! 2. Surface boundary condition on tke and its eddy viscosity (zmxlm)
+      ! -------------------------------------------------
+      ! en(1)   = rn_ebb sqrt(utau^2+vtau^2) / rau0  (min value rn_emin0)
+      ! zmxlm(1) = avmb(1) and zmxlm(jpk) = 0.
+!CDIR NOVERRCHK
+      DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
+            zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
+            zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 )
+            en (ji,jj,1) = MAX( zesurf, rn_emin0 ) * tmask(ji,jj,1)
+            zav  =  rn_ediff * zmxlm(ji,jj,1) * SQRT( en(ji,jj,1) )
+            zmxlm(ji,jj,1  ) = MAX( zav, avmb(1)                  ) * tmask(ji,jj,1)
+            avt  (ji,jj,1  ) = MAX( zav, avtb(1) * avtb_2d(ji,jj) ) * tmask(ji,jj,1)
+            zmxlm(ji,jj,jpk) = 0.e0
+         END DO
+      END DO
+
+      ! 3. Now Turbulent kinetic energy (output in en)
+      ! -------------------------------
+      ! Resolution of a tridiagonal linear system by a "methode de chasse"
+      ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ). 
+
+      SELECT CASE ( nn_pdl )
+      !
+      CASE ( 0 )           ! No Prandtl number
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  !                                                             ! shear prod. - stratif. destruction
+                  zcoef = 0.5 / fse3w(ji,jj,jk)
+                  zdku = zcoef * (  ub(ji-1, jj ,jk-1) + ub(ji,jj,jk-1)   &          ! shear
+                     &            - ub(ji-1, jj ,jk  ) - ub(ji,jj,jk  )  )
+                  zdkv = zcoef * (  vb( ji ,jj-1,jk-1) + vb(ji,jj,jk-1)   &
+                     &            - vb( ji ,jj-1,jk  ) - vb(ji,jj,jk  )  )
+                  zesh2 =  eboost * ( zdku*zdku + zdkv*zdkv ) - rn2b(ji,jj,jk)        ! coefficient (zesh2)
+                  !
+                  !                                                             ! Matrix
+                  zcof = zfact1 * tmask(ji,jj,jk)
+                  !                                                                  ! lower diagonal
+                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   &
+                     &                  / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
+                  !                                                                  ! upper diagonal
+                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   &
+                     &                  / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
+                  !                                                                  ! diagonal
+                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
+                  !
+                  !                                                             ! right hand side in en 
+                  IF( .NOT. ln_lc ) THEN                                             ! No Langmuir cells
+                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en(ji,jj,jk)   &
+                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2
+                  ELSE                                                               ! Langmuir cell source term
+                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en(ji,jj,jk)   &
+                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2          &
+                        &                        +   rdt  * ztkelc(ji,jj,jk)
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         !
+      CASE ( 1 )           ! Prandtl number
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  !                                                             ! shear prod. - stratif. destruction
+                  zcoef = 0.5 / fse3w(ji,jj,jk)                                     
+                  zdku = zcoef * (  ub(ji-1,jj  ,jk-1) + ub(ji,jj,jk-1)   &          ! shear
+                  &               - ub(ji-1,jj  ,jk  ) - ub(ji,jj,jk  )  )
+                  zdkv = zcoef * (  vb(ji  ,jj-1,jk-1) + vb(ji,jj,jk-1)   &
+                  &               - vb(ji  ,jj-1,jk  ) - vb(ji,jj,jk  )  )
+                  zsh2 = zdku * zdku + zdkv * zdkv                                   ! square of shear
+                  zri  = MAX( rn2b(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 )                ! local Richardson number
+# if defined key_c1d
+                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri 
+# endif
+                  zpdl = 1.0                                                         ! Prandtl number
+                  IF( zri >= 0.2 )   zpdl = 0.2 / zri
+                  zpdl = MAX( 0.1, zpdl )
+                  zesh2 = eboost * zsh2 - zpdl * rn2b(ji,jj,jk)                       ! coefficient (esh2)
+                  !
+                  !                                                             ! Matrix
+                  zcof = zfact1 * tmask(ji,jj,jk)
+                  !                                                                 ! lower diagonal
+                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   &
+                  &                     / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
+                  !                                                                 ! upper diagonal
+                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   &
+                  &                     / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
+                  !                                                                 ! diagonal
+                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
+                  !
+                  !                                                             ! right hand side in en 
+                  IF( .NOT. ln_lc ) THEN                                             ! No Langmuir cells
+                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en   (ji,jj,jk)   &
+                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2
+                  ELSE                                                               ! Langmuir cell source term
+                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en   (ji,jj,jk)   &
+                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2             &
+                        &                        +   rdt  * ztkelc(ji,jj,jk)
+                  ENDIF
+                  zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk)                      ! store masked Prandlt number in zmxld array
+               END DO
+            END DO
+         END DO
+         !
+      END SELECT
+
+# if defined key_c1d
+      e_pdl(:,:,2:jpkm1) = zmxld(:,:,2:jpkm1)      ! c1d configuration : save masked Prandlt number
+      e_pdl(:,:,      1) = e_pdl(:,:,      2)
+      e_pdl(:,:,    jpk) = e_pdl(:,:,  jpkm1)      
+# endif
+
+      ! 4. Matrix inversion from level 2 (tke prescribed at level 1)
+      !!--------------------------------
+      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1    ! vector opt.
+               zwd(ji,jj,jk) = zwd(ji,jj,jk) - avmv(ji,jj,jk) * avmu(ji,jj,jk-1) / zwd(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+         DO ji = fs_2, fs_jpim1    ! vector opt.
+            avmv(ji,jj,2) = en(ji,jj,2) - avmv(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
+         END DO
+      END DO
+      DO jk = 3, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1    ! vector opt.
+               avmv(ji,jj,jk) = en(ji,jj,jk) - avmv(ji,jj,jk) / zwd(ji,jj,jk-1) *avmv(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+         DO ji = fs_2, fs_jpim1    ! vector opt.
+            en(ji,jj,jpkm1) = avmv(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
+         END DO
+      END DO
+      DO jk = jpk-2, 2, -1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1    ! vector opt.
+               en(ji,jj,jk) = ( avmv(ji,jj,jk) - avmu(ji,jj,jk) * en(ji,jj,jk+1) ) / zwd(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      DO jk = 2, jpkm1                             ! set the minimum value of tke
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+     
+      ! 5. Add extra TKE due to surface and internal wave breaking (nn_etau /= 0)
+      !!----------------------------------------------------------
+      IF( nn_etau /= 0 ) THEN        ! extra tke : en = en + rn_efr * en(1) * exp( -z/zhtau )
+         !
+         SELECT CASE( nn_htau )           ! Choice of the depth of penetration
+         CASE( 0 )                                    ! constant depth penetration (here 10 meters)
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zhtau(ji,jj) = 10.
+               END DO
+            END DO
+         CASE( 1 )                                    ! meridional profile  1
+            DO jj = 2, jpjm1                          ! ( 5m in the tropics to a maximum of 40 m at high lat.)
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zhtau(ji,jj) = MAX( 5., MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
+               END DO
+            END DO
+         CASE( 2 )                                    ! meridional profile 2
+            DO jj = 2, jpjm1                          ! ( 5m in the tropics to a maximum of 60 m at high lat.)
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zhtau(ji,jj) =  MAX( 5.,6./4.* MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
+               END DO
+            END DO
+         END SELECT
+         ! 
+         IF( nn_etau == 1 ) THEN          ! extra term throughout the water column
+            DO jk = 2, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     en(ji,jj,jk) = en(ji,jj,jk)   &
+                        &         + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) &
+                        &                  * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         ELSEIF( nn_etau == 2 ) THEN      ! extra term only at the base of the mixed layer
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  jk = nmln(ji,jj)
+                  en(ji,jj,jk) = en(ji,jj,jk)    &
+                     &         + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) &
+                     &                   * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
+               END DO
+            END DO
+         ENDIF
+         !
+      ENDIF
+
+
+      ! Lateral boundary conditions (sign unchanged)
+      CALL lbc_lnk( en, 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. )
+
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! IV.  Before vertical eddy vicosity and diffusivity coefficients
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      SELECT CASE ( nn_ave )
+      CASE ( 0 )                                     ! no horizontal average 
+         DO jk = 2, jpkm1                                 ! only vertical eddy viscosity 
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  avmu(ji,jj,jk) = ( avt  (ji,jj,jk) + avt  (ji+1,jj  ,jk) ) * umask(ji,jj,jk)   &
+                  &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk) )
+                  avmv(ji,jj,jk) = ( avt  (ji,jj,jk) + avt  (ji  ,jj+1,jk) ) * vmask(ji,jj,jk)   &
+                  &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk) )
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
+         !
+      CASE ( 1 )                                     ! horizontal average                         ( 1/2  1/2 )
+         !                                                ! Vertical eddy viscosity    avmu = 1/4 ( 1    1   )
+         !                                                !                                       ( 1/2  1/2 )
+         !                                                !
+         !                                                !                                       ( 1/2  1   1/2 )            
+         !                                                !                            avmv = 1/4 ( 1/2  1   1/2 )      
+         DO jk = 2, jpkm1                                
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+# if defined key_vectopt_memory
+                  !                                       ! caution : vectopt_memory change the solution
+                  !                                       !           (last digit of the solver stat)
+                  avmu(ji,jj,jk) = (      avt(ji,jj  ,jk) + avt(ji+1,jj  ,jk)   &
+                     &              +.5*( avt(ji,jj-1,jk) + avt(ji+1,jj-1,jk)   &
+                     &                   +avt(ji,jj+1,jk) + avt(ji+1,jj+1,jk) ) ) * eumean(ji,jj,jk)
+
+                  avmv(ji,jj,jk) = (      avt(ji  ,jj,jk) + avt(ji  ,jj+1,jk)   &
+                     &              +.5*( avt(ji-1,jj,jk) + avt(ji-1,jj+1,jk)   &
+                     &                   +avt(ji+1,jj,jk) + avt(ji+1,jj+1,jk) ) ) * evmean(ji,jj,jk)
+# else
+                  avmu(ji,jj,jk) = (   avt  (ji,jj  ,jk) + avt  (ji+1,jj  ,jk)   &
+                     &           +.5*( avt  (ji,jj-1,jk) + avt  (ji+1,jj-1,jk)   &
+                     &                +avt  (ji,jj+1,jk) + avt  (ji+1,jj+1,jk) ) ) * umask(ji,jj,jk)   &
+                     &    / MAX( 1.,   tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk)   &
+                     &           +.5*( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk)   &
+                     &                +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) )  )
+
+                  avmv(ji,jj,jk) = (   avt  (ji  ,jj,jk) + avt  (ji  ,jj+1,jk)   &
+                     &           +.5*( avt  (ji-1,jj,jk) + avt  (ji-1,jj+1,jk)   &
+                     &                +avt  (ji+1,jj,jk) + avt  (ji+1,jj+1,jk) ) ) * vmask(ji,jj,jk)   &
+                     &   /  MAX( 1.,   tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
+                     &           +.5*( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk)   &
+                     &                +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) )  )
+# endif
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions 
+         !
+         !                                                ! Vertical eddy diffusivity             (1 2 1)
+         !                                                !                            avt = 1/16 (2 4 2)
+         !                                                !                                       (1 2 1)
+         DO jk = 2, jpkm1         
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+#  if defined key_vectopt_memory
+                  avt(ji,jj,jk) = ( avmu(ji,jj,jk) + avmu(ji-1,jj  ,jk)                          &
+                     &            + avmv(ji,jj,jk) + avmv(ji  ,jj-1,jk)  ) * etmean(ji,jj,jk)
+#  else
+                  avt(ji,jj,jk) = ( avmu (ji,jj,jk) + avmu (ji-1,jj  ,jk)                        &
+                     &            + avmv (ji,jj,jk) + avmv (ji  ,jj-1,jk)  ) * tmask(ji,jj,jk)   &
+                     &  / MAX( 1.,  umask(ji,jj,jk) + umask(ji-1,jj  ,jk)                        &
+                     &            + vmask(ji,jj,jk) + vmask(ji  ,jj-1,jk)  )
+#  endif
+               END DO
+            END DO
+         END DO
+         !
+      END SELECT
+      !
+      IF( nn_pdl == 1 ) THEN                            ! Ponderation by the Prandtl number (nn_pdl=1)
+         DO jk = 2, jpkm1 
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zpdl = zmxld(ji,jj,jk)
+                  avt(ji,jj,jk) = MAX( zpdl * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ENDIF
+      !
+      DO jk = 2, jpkm1                                  ! Minimum value on the eddy viscosity
+         avmu(:,:,jk) = MAX( avmu(:,:,jk), avmb(jk) ) * umask(:,:,jk)
+         avmv(:,:,jk) = MAX( avmv(:,:,jk), avmb(jk) ) * vmask(:,:,jk)
+      END DO
+
+      CALL lbc_lnk( avt, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged)
+
+      IF( lrst_oce )   CALL tke_rst( kt, 'WRITE' )      ! write en in restart file
+
+      IF(ln_ctl) THEN
+         CALL prt_ctl( tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
+         CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask,                   &
+            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
+      ENDIF
+      !
+   END SUBROUTINE zdf_tke_old
+
+
+   SUBROUTINE zdf_tke_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_tke_init  ***
+      !!                     
+      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
+      !!      viscosity when using a tke turbulent closure scheme
+      !!
+      !! ** Method  :   Read the namtke namelist and check the parameters
+      !!      called at the first timestep (nit000)
+      !!
+      !! ** input   :   Namlist namtke
+      !!
+      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
+      !!
+      !!----------------------------------------------------------------------
+      USE dynzdf_exp
+      USE trazdf_exp
+      !
+# if defined key_vectopt_memory
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+# else
+      INTEGER ::           jk   ! dummy loop indices
+# endif
+      !!
+      NAMELIST/namtke/ ln_rstke, rn_ediff, rn_ediss, rn_ebb  , rn_efave, rn_emin,   &
+         &             rn_emin0, rn_cri  , nn_itke , nn_mxl  , nn_pdl  , nn_ave ,   &
+         &             nn_avb  , ln_mxl0 , rn_lmin , rn_lmin0, nn_havtb, nn_etau,   &
+         &             nn_htau , rn_efr  , ln_lc   , rn_lc 
+      !!----------------------------------------------------------------------
+
+      ! Read Namelist namtke : Turbulente Kinetic Energy
+      ! --------------------
+      REWIND ( numnam )
+      READ   ( numnam, namtke )
+
+      ! Compute boost associated with the Richardson critic
+      !     (control values: rn_cri = 0.3   ==> eboost=1.25 for nn_pdl=1)
+      !     (                rn_cri = 0.222 ==> eboost=1.               )
+      eboost = rn_cri * ( 2. + rn_ediss / rn_ediff ) / 2.
+
+
+
+      ! Parameter control and print
+      ! ---------------------------
+      ! Control print
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme (old scheme)'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '          Namelist namtke : set tke mixing parameters'
+         WRITE(numout,*) '             restart with tke from no tke              ln_rstke = ', ln_rstke
+         WRITE(numout,*) '             coef. to compute avt                      rn_ediff = ', rn_ediff
+         WRITE(numout,*) '             Kolmogoroff dissipation coef.             rn_ediss = ', rn_ediss
+         WRITE(numout,*) '             tke surface input coef.                   rn_ebb   = ', rn_ebb
+         WRITE(numout,*) '             tke diffusion coef.                       rn_efave = ', rn_efave
+         WRITE(numout,*) '             minimum value of tke                      rn_emin  = ', rn_emin
+         WRITE(numout,*) '             surface minimum value of tke              rn_emin0 = ', rn_emin0
+         WRITE(numout,*) '             number of restart iter loops              nn_itke  = ', nn_itke
+         WRITE(numout,*) '             mixing length type                        nn_mxl   = ', nn_mxl
+         WRITE(numout,*) '             prandl number flag                        nn_pdl   = ', nn_pdl
+         WRITE(numout,*) '             horizontal average flag                   nn_ave   = ', nn_ave
+         WRITE(numout,*) '             critic Richardson nb                      rn_cri   = ', rn_cri
+         WRITE(numout,*) '                and its associated coeff.              eboost   = ', eboost
+         WRITE(numout,*) '             constant background or profile            nn_avb   = ', nn_avb
+         WRITE(numout,*) '             surface mixing length = F(stress) or not  ln_mxl0  = ', ln_mxl0
+         WRITE(numout,*) '             surface  mixing length minimum value      rn_lmin0 = ', rn_lmin0
+         WRITE(numout,*) '             interior mixing length minimum value      rn_lmin0 = ', rn_lmin0
+         WRITE(numout,*) '             horizontal variation for avtb             nn_havtb = ', nn_havtb
+         WRITE(numout,*) '             test param. to add tke induced by wind    nn_etau  = ', nn_etau
+         WRITE(numout,*) '             flag for computation of exp. tke profile  nn_htau  = ', nn_htau
+         WRITE(numout,*) '             % of rn_emin0 which pene. the thermocline rn_efr   = ', rn_efr
+         WRITE(numout,*) '             flag to take into acc.  Langmuir circ.    ln_lc    = ', ln_lc
+         WRITE(numout,*) '             coef to compute verticla velocity of LC   rn_lc    = ', rn_lc
+         WRITE(numout,*)
+      ENDIF
+
+      ! Check of some namelist values
+      IF( nn_mxl  < 0    .OR. nn_mxl  > 3   )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
+      IF( nn_pdl  < 0    .OR. nn_pdl  > 1   )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
+      IF( nn_ave  < 0    .OR. nn_ave  > 1   )   CALL ctl_stop( 'bad flag: nn_ave is  0 or 1    ' )
+      IF( nn_htau < 0    .OR. nn_htau > 2   )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
+      IF( rn_lc   < 0.15 .OR. rn_lc   > 0.2 )   CALL ctl_stop( 'bad value: rn_lc must be between 0.15 and 0.2 ' )
+
+      IF( nn_etau == 2 )   CALL zdf_mxl( nit000 )      ! Initialization of nmln 
+
+      ! Horizontal average : initialization of weighting arrays 
+      ! -------------------
+      !
+      SELECT CASE ( nn_ave )
+      ! Notice: mean arrays have not to by defined at local domain boundaries due to the vertical nature of TKE
+      !
+      CASE ( 0 )                ! no horizontal average
+         IF(lwp) WRITE(numout,*) '          no horizontal average on avt, avmu, avmv'
+         IF(lwp) WRITE(numout,*) '          only in very high horizontal resolution !'
+# if defined key_vectopt_memory
+         ! caution vectopt_memory change the solution (last digit of the solver stat)
+         ! weighting mean arrays etmean, eumean and evmean
+         !           ( 1  1 )                                          ( 1 )
+         ! avt = 1/4 ( 1  1 )     avmu = 1/2 ( 1  1 )       avmv=  1/2 ( 1 )
+         !                          
+         etmean(:,:,:) = 0.e0
+         eumean(:,:,:) = 0.e0
+         evmean(:,:,:) = 0.e0
+         !
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  etmean(ji,jj,jk) = tmask(ji,jj,jk) / MAX( 1.,  umask(ji,jj,jk) + umask(ji-1,jj  ,jk)   &
+                  &                                            + vmask(ji,jj,jk) + vmask(ji  ,jj-1,jk) )
+                  eumean(ji,jj,jk) = umask(ji,jj,jk) / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk) )
+                  evmean(ji,jj,jk) = vmask(ji,jj,jk) / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk) )
+               END DO
+            END DO
+         END DO
+# endif
+         !
+      CASE ( 1 )                ! horizontal average 
+         IF(lwp) WRITE(numout,*) '          horizontal average on avt, avmu, avmv'
+# if defined key_vectopt_memory
+         ! caution vectopt_memory change the solution (last digit of the solver stat)
+         ! weighting mean arrays etmean, eumean and evmean
+         !           ( 1  1 )              ( 1/2  1/2 )             ( 1/2  1  1/2 )
+         ! avt = 1/4 ( 1  1 )   avmu = 1/4 ( 1    1   )   avmv= 1/4 ( 1/2  1  1/2 )
+         !                                 ( 1/2  1/2 )
+         etmean(:,:,:) = 0.e0
+         eumean(:,:,:) = 0.e0
+         evmean(:,:,:) = 0.e0
+         !
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  etmean(ji,jj,jk) = tmask(ji,jj,jk) / MAX( 1.,   umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk)   &
+                  &                                       +       vmask(ji  ,jj-1,jk) + vmask(ji  ,jj  ,jk) )
+                  eumean(ji,jj,jk) = umask(ji,jj,jk) / MAX( 1.,   tmask(ji  ,jj  ,jk) + tmask(ji+1,jj  ,jk)   &
+                  &                                       +.5 * ( tmask(ji  ,jj-1,jk) + tmask(ji+1,jj-1,jk)   &
+                  &                                             + tmask(ji  ,jj+1,jk) + tmask(ji+1,jj+1,jk) )  )
+                  evmean(ji,jj,jk) = vmask(ji,jj,jk) / MAX( 1.,   tmask(ji  ,jj  ,jk) + tmask(ji  ,jj+1,jk)   &
+                  &                                       +.5 * ( tmask(ji-1,jj  ,jk) + tmask(ji-1,jj+1,jk)   &
+                  &                                             + tmask(ji+1,jj  ,jk) + tmask(ji+1,jj+1,jk) )  )
+               END DO
+            END DO
+         END DO
+# endif
+         !
+      END SELECT
+
+
+      ! Background eddy viscosity and diffusivity profil
+      ! ------------------------------------------------
+      IF( nn_avb == 0 ) THEN      ! Define avmb, avtb from namelist parameter
+         avmb(:) = avm0
+         avtb(:) = avt0
+      ELSE                      ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990) 
+         avmb(:) = avm0
+         avtb(:) = avt0 + ( 3.0e-4 - 2 * avt0 ) * 1.0e-4 * gdepw_0(:)   ! m2/s
+         IF(ln_sco .AND. lwp)   CALL ctl_warn( '          avtb profile not valid in sco' )
+      ENDIF
+      !
+      !                         ! 2D shape of the avtb
+      avtb_2d(:,:) = 1.e0             ! uniform 
+      !
+      IF( nn_havtb == 1 ) THEN      ! decrease avtb in the equatorial band
+           !  -15S -5S : linear decrease from avt0 to avt0/10.
+           !  -5S  +5N : cst value avt0/10.
+           !   5N  15N : linear increase from avt0/10, to avt0
+           WHERE(-15. <= gphit .AND. gphit < -5 )   avtb_2d = (1.  - 0.09 * (gphit + 15.))
+           WHERE( -5. <= gphit .AND. gphit <  5 )   avtb_2d =  0.1
+           WHERE(  5. <= gphit .AND. gphit < 15 )   avtb_2d = (0.1 + 0.09 * (gphit -  5.))
+      ENDIF
+
+      ! Initialization of vertical eddy coef. to the background value
+      ! -------------------------------------------------------------
+      DO jk = 1, jpk
+         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
+         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
+         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
+      END DO
+
+
+      ! read or initialize turbulent kinetic energy ( en )
+      ! -------------------------------------------------
+      CALL tke_rst( nit000, 'READ' )
+      !
+   END SUBROUTINE zdf_tke_init
+
+
+   SUBROUTINE tke_rst( kt, cdrw )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE ts_rst  ***
+     !!                     
+     !! ** Purpose :   Read or write TKE file (en) in restart file
+     !!
+     !! ** Method  :   use of IOM library
+     !!                if the restart does not contain TKE, en is either 
+     !!                set to rn_emin or recomputed (nn_itke/=0)
+     !!----------------------------------------------------------------------
+     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
+     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
+     !
+     INTEGER ::   jit   ! dummy loop indices
+     !!----------------------------------------------------------------------
+     !
+     IF( TRIM(cdrw) == 'READ' ) THEN
+        IF( ln_rstart ) THEN
+           IF( iom_varid( numror, 'en', ldstop = .FALSE. ) > 0 .AND. .NOT.(ln_rstke) ) THEN 
+              CALL iom_get( numror, jpdom_autoglo, 'en', en )
+           ELSE
+              IF( lwp .AND. iom_varid( numror, 'en', ldstop = .FALSE. ) > 0 )   &
+                 &                       WRITE(numout,*) ' ===>>>> : previous run without tke scheme'
+              IF( lwp .AND. ln_rstke )   WRITE(numout,*) ' ===>>>> : We do not use en from the restart file'
+              IF( lwp                )   WRITE(numout,*) ' ===>>>> : en is set by iterative loop'
+              en (:,:,:) = rn_emin * tmask(:,:,:)
+              DO jit = 2, nn_itke + 1
+                 CALL zdf_tke_old( jit )
+              END DO
+           ENDIF
+        ELSE
+           en(:,:,:) = rn_emin * tmask(:,:,:)      ! no restart: en set to its minimum value
+        ENDIF
+     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
+        CALL iom_rstput( kt, nitrst, numrow, 'en', en )
+     ENDIF
+     !
+   END SUBROUTINE tke_rst
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module :                                        NO TKE scheme
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke_old = .FALSE.   !: TKE flag
+CONTAINS
+   SUBROUTINE zdf_tke_old( kt )          ! Empty routine
+      WRITE(*,*) 'zdf_tke_old: You should not have seen this print! error?', kt
+   END SUBROUTINE zdf_tke_old
+#endif
+
+   !!======================================================================
+END MODULE zdftke_old
