Index: trunk/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 1059)
+++ trunk/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 1060)
@@ -2,25 +2,25 @@
    !!======================================================================
    !!                       ***  MODULE  zdftke  ***
-   !! Ocean physics:  vertical mixing coefficient compute from the tke 
+   !! Ocean physics:  vertical mixing coefficient computed from the tke 
    !!                 turbulent closure parameterization
    !!=====================================================================
-   !! History :   6.0  !  91-03  (b. blanke)  Original code
-   !!             7.0  !  91-11  (G. Madec)   bug fix
-   !!             7.1  !  92-10  (G. Madec)   new mixing length and eav
-   !!             7.2  !  93-03  (M. Guyon)   symetrical conditions
-   !!             7.3  !  94-08  (G. Madec, M. Imbard)   npdl flag
-   !!             7.5  !  96-01  (G. Madec)   s-coordinates
-   !!             8.0  !  97-07  (G. Madec)   lbc
-   !!             8.1  !  99-01  (E. Stretta) new option for the mixing length
-   !!             8.5  !  02-06  (G. Madec) add zdf_tke_init routine
-   !!             8.5  !  02-08  (G. Madec)  ri_c and Free form, F90
-   !!             9.0  !  04-10  (C. Ethe )  1D configuration
-   !!             9.0  !  02-08  (G. Madec)  autotasking optimization
-   !!             9.0  !  06-07  (S. Masson)  distributed restart using iom
-   !!              -   !  2005-07  (C. Ethe,  G.Madec) : update TKE physics:
+   !! 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
-   !!             3.0  !  2007-11  (G. Madec)  avtb_2d, nn_avtb_2d 
+   !!              -   !  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   ||   defined key_esopa
@@ -39,14 +39,14 @@
    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
-   USE restart         ! only for lrst_oce
+   USE iom             ! I/O manager library
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   zdf_tke        ! routine called in step module
+   PUBLIC   zdf_tke    ! routine called in step module
 
    LOGICAL , PUBLIC, PARAMETER              ::   lk_zdftke = .TRUE.  !: TKE vertical mixing flag
@@ -54,40 +54,39 @@
    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
-
-   !! * Namelist (namtke)
-   LOGICAL , PUBLIC ::   ln_rstke = .FALSE.         !: =T restart with tke from a run without tke with 
-     !                                              !  a none zero initial value for en
-   INTEGER , PUBLIC ::   nitke = 50 ,            &  !: number of restart iterative loops
-      &                  nmxl  =  2 ,            &  !: = 0/1/2/3 flag for the type of mixing length used
-      &                  npdl  =  1 ,            &  !: = 0/1/2 flag on prandtl number on vert. eddy coeff.
-      &                  nave  =  1 ,            &  !: = 0/1 flag for horizontal average on avt, avmu, avmv
-      &                  navb  =  0                 !: = 0/1 flag for constant or profile background avt
-   REAL(wp), PUBLIC ::   ediff = 0.1_wp       ,  &  !: coeff. for vertical eddy coef.; avt=ediff*mxl*sqrt(e)
-      &                  ediss = 0.7_wp       ,  &  !: coef. of the Kolmogoroff dissipation 
-      &                  ebb   = 3.75_wp      ,  &  !: coef. of the surface input of tke
-      &                  efave = 1._wp        ,  &  !: coef. for the tke vert. diff. coeff.; avtke=efave*avm
-      &                  emin  = 0.7071e-6_wp ,  &  !: minimum value of tke (m2/s2)
-      &                  emin0 = 1.e-4_wp     ,  &  !: surface minimum value of tke (m2/s2)
-      &                  ri_c  = 2._wp / 9._wp      !: critic Richardson number
-   !                                      !!! ** Namelist (namtke) **
-   INTEGER  ::   nn_avtb_2d = 1            !: = 0/1 horizontal shape for avtb
-   INTEGER  ::   nn_etau    = 0            !: = 0/1/2 tke depth penetration
-   INTEGER  ::   nn_htau    = 0            !: = 0/1/2 type of tke profile of penetration
-   REAL(wp) ::   rn_lmin    = 0.4_wp       !: surface min value of mixing turbulent length scale
-   REAL(wp) ::   rn_efr     = 1.0_wp       !: fraction of TKE surface value which penetrates in the ocean
-   LOGICAL  ::   ln_lsfc    = .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
-   REAL(wp) ::   rn_lc      = 0.15_wp      !: coef to compute vertical velocity of LC
-
-   REAL(wp), DIMENSION (jpi,jpj)  :: avtb_2d ! set in tke_init, for other modif than ice
-
-# if defined key_c1d
-   !                                                                   ! 1D cfg only
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_dis, e_mix,      &  ! dissipation and mixing turbulent lengh scales
-      &                                          e_pdl, e_ric          ! prandl and local Richardson numbers
+#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
+
+   !                                       !!! ** 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.4_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
@@ -95,5 +94,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
+   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -111,23 +110,24 @@
       !! ** 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( efave eav d(en)/dz )/dz    ! diffusion of tke
-      !!                  + grav/rau0 pdl eav d(rau)/dz   ! stratif. destruc.
-      !!                  - ediss / emxl en**(2/3)        ! dissipation
+      !!         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( emin0,ebb sqrt(utau^2 + vtau^2) )
-      !!         bottom : en = emin
+      !!         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= 1/(sqrt(en)/N)  WHERE N is the brunt-vaisala frequency
+      !!         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 : 
-      !!         nmxl=0 : mxl bounded by the distance to surface and bottom.
+      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
       !!                  zmxld = zmxlm = mxl
-      !!         nmxl=1 : mxl bounded by the vertical scale factor.
+      !!         nn_mxl=1 : mxl bounded by the vertical scale factor.
       !!                  zmxld = zmxlm = mxl
-      !!         nmxl=2 : mxl bounded such that the vertical derivative of mxl 
+      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl 
       !!                  is less than 1 (|d/dz(xml)|<1). 
       !!                  zmxld = zmxlm = mxl
-      !!         nmxl=3 : lup = mxl bounded using |d/dz(xml)|<1 from the surface 
+      !!         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 
@@ -140,13 +140,13 @@
       !!      solved.
       !!         Note that - the shear production is multiplied by eboost in order
-      !!      to set the critic richardson number to ri_c (namelist parameter)
+      !!      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 npdl=1 (namelist parameter)
+      !!      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, ediff*zmxlm*en^1/2 )
-      !!              avt = max( avmb, pdl*avm )  (pdl=1 if npdl=0)
+      !!              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-
@@ -159,27 +159,27 @@
       !!                update avt, avmu, avmv (before vertical eddy coef.)
       !!
-      !! References : Gaspar et al., jgr, 95, 1990,
-      !!              Blanke and Delecluse, jpo, 1991
+      !! 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
-         &          zmxlm => ta,  &  ! use ta as workspace
-         &          zmxld => sa      ! use sa as workspace
-      USE oce     , ztkelc => va     ! use va as workspace
-      !
-      INTEGER, INTENT(in) ::   kt ! ocean time step
-      !
-      INTEGER  ::   ji, jj, jk                  ! dummy loop arguments
-      REAL(wp) ::   zmlmin, zbbrau,          &  ! temporary scalars
-         &          zfact1, zfact2, zfact3,  &  !
-         &          zrn2, zesurf,            &  !
-         &          ztx2, zty2, zav,         &  !
-         &          zcoef, zcof, zsh2,       &  !
-         &          zdku, zdkv, zpdl, zri,   &  !
-         &          zsqen, zesh2,            &  !
-         &          zemxl, zemlm, zemlp
+      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) ::   zraug, zus, zwlc, zind     ! temporary scalar
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zhtau   ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc    ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zhtau   !  -      -
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc    !  -      -
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc   ! 3D workspace
       !!--------------------------------------------------------------------
@@ -188,9 +188,8 @@
 
       !                                            ! Local constant initialization
-      zmlmin = 0.4
-      zbbrau =  .5 * ebb / rau0
-      zfact1 = -.5 * rdt * efave
-      zfact2 = 1.5 * rdt * ediss
-      zfact3 = 0.5 * rdt * ediss
+      zbbrau =  .5 * rn_ebb / rau0
+      zfact1 = -.5 * rdt * rn_efave
+      zfact2 = 1.5 * rdt * rn_ediss
+      zfact3 = 0.5 * rdt * rn_ediss
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -200,6 +199,8 @@
       ! Buoyancy length scale: l=sqrt(2*e/n**2)
       ! ---------------------
-      IF( ln_lsfc ) THEN      ! lsurf is function of wind stress : l=2.e5*sqrt(utau^2 + vtau^2)/(rau0*g)
+      IF( ln_mxl0 ) THEN         ! surface mixing length = F(stress) : l=2.e5*sqrt(utau^2 + vtau^2)/(rau0*g)
+!!gm  this should be useless
          zmxlm(:,:,1) = 0.e0
+!!gm end
          zraug = 0.5 * 2.e5 / ( rau0 * grav )
          DO jj = 2, jpjm1
@@ -208,15 +209,14 @@
                ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
                zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
-               zmxlm(ji,jj,1  ) = zraug * SQRT( ztx2 * ztx2 + zty2 * zty2 )
-               ! set the surface minimum value to rn_lmin
-               zmxlm(ji,jj,1  ) = MAX( zmxlm(ji,jj,1) , rn_lmin )
-            END DO
-         END DO
-      ELSE                    ! surface set to the minimum value
-         zmxlm(:,:,1) = rn_lmin  
+               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 minimum value
-!CDIR NOVERRCHK
-      DO jk = 2, jpkm1
+      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
@@ -224,5 +224,5 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zrn2 = MAX( rn2(ji,jj,jk), rsmall )
-               zmxlm(ji,jj,jk) = MAX( SQRT( 2. * en(ji,jj,jk) / zrn2 ), zmlmin  )
+               zmxlm(ji,jj,jk) = MAX(  rn_lmin,  SQRT( 2. * en(ji,jj,jk) / zrn2 )  )
             END DO
          END DO
@@ -232,8 +232,8 @@
       ! -------------------------------------
       zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value
-      zmxld(:,:,jpk) = rn_lmin           ! bottom  set to the minimum value
-
-      SELECT CASE ( nmxl )
-
+      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
@@ -247,5 +247,5 @@
             END DO
          END DO
-
+         !
       CASE ( 1 )           ! bounded by the vertical scale factor
          DO jk = 2, jpkm1
@@ -258,5 +258,5 @@
             END DO
          END DO
-
+         !
       CASE ( 2 )           ! |dk[xml]| bounded by e3t :
          DO jk = 2, jpkm1         ! from the surface to the bottom :
@@ -276,5 +276,5 @@
             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
@@ -305,9 +305,9 @@
             END DO
          END DO
-
+         !
       END SELECT
 
 # if defined key_c1d
-      ! save mixing and dissipation turbulent length scales
+      ! c1d configuration : save mixing and dissipation turbulent length scales
       e_dis(:,:,:) = zmxld(:,:,:)
       e_mix(:,:,:) = zmxlm(:,:,:)
@@ -346,5 +346,5 @@
          !
 # if defined key_c1d
-         hlc(:,:) = zhlc(:,:) * tmask(:,:,1)      ! save finite Langmuir Circulation  depth
+         hlc(:,:) = zhlc(:,:) * tmask(:,:,1)      ! c1d configuration: save finite Langmuir Circulation depth
 # endif
          !
@@ -382,5 +382,5 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zsqen = SQRT( en(ji,jj,jk) )
-               zav   = ediff * zmxlm(ji,jj,jk) * zsqen
+               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)
@@ -392,5 +392,5 @@
       ! 2. Surface boundary condition on tke and its eddy viscosity (zmxlm)
       ! -------------------------------------------------
-      ! en(1)   = ebb sqrt(utau^2+vtau^2) / rau0  (min value emin0)
+      ! en(1)   = rn_ebb sqrt(utau^2+vtau^2) / rau0  (min value rn_emin0)
       ! zmxlm(1) = avmb(1) and zmxlm(jpk) = 0.
 !CDIR NOVERRCHK
@@ -401,8 +401,8 @@
             zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
             zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 )
-            en (ji,jj,1) = MAX( zesurf, emin0 ) * tmask(ji,jj,1)
-            zav  =  ediff * zmxlm(ji,jj,1) * SQRT( en(ji,jj,1) )* tmask(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)
+            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
@@ -412,38 +412,36 @@
       ! -------------------------------
       ! 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 ( npdl )
-
+      ! 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.
-                  ! zesh2 = eboost * (du/dz)^2 - N^2
+                  !                                                             ! shear prod. - stratif. destruction
                   zcoef = 0.5 / fse3w(ji,jj,jk)
-                  ! shear
-                  zdku = zcoef * (   ub(ji-1, jj ,jk-1) + ub(ji,jj,jk-1)   &
-                  &                - 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  )  )
-                  ! coefficient (zesh2)
-                  zesh2 =  eboost * ( zdku*zdku + zdkv*zdkv ) - rn2(ji,jj,jk)
-
-                  ! Matrix
+                  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 ) - rn2(ji,jj,jk)        ! coefficient (zesh2)
+                  !
+                  !                                                             ! Matrix
                   zcof = zfact1 * tmask(ji,jj,jk)
-                  ! lower diagonal
+                  !                                                                  ! 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
+                     &                  / ( 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
+                     &                  / ( 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
+                  !
+                  !                                                             ! 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
+                  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          &
@@ -453,61 +451,54 @@
             END DO
          END DO
-
+         !
       CASE ( 1 )           ! Prandtl number
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ! zesh2 =  eboost * (du/dz)^2 - pdl * N^2
-                  zcoef = 0.5 / fse3w(ji,jj,jk)
-                  ! shear
-                  zdku = zcoef * (   ub(ji-1,jj  ,jk-1) + ub(ji,jj,jk-1) &
-                  &                - 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  )   )
-                  ! square of vertical shear
-                  zsh2 = zdku * zdku + zdkv * zdkv
-                  ! local Richardson number
-                  zri  = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 )
+                  !                                                             ! 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( rn2(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 )                ! local Richardson number
 # if defined key_c1d
-                  ! save masked local Richardson number in zmxlm array
-                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)
+                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri 
 # endif
-                  ! Prandtl number
-                  zpdl = 1.0
-                  IF( zri >= 0.2 ) zpdl = 0.2 / zri
+                  zpdl = 1.0                                                         ! Prandtl number
+                  IF( zri >= 0.2 )   zpdl = 0.2 / zri
                   zpdl = MAX( 0.1, zpdl )
-                  ! coefficient (esh2)
-                  zesh2 = eboost * zsh2 - zpdl * rn2(ji,jj,jk)
-
-                  ! Matrix
+                  zesh2 = eboost * zsh2 - zpdl * rn2(ji,jj,jk)                       ! coefficient (esh2)
+                  !
+                  !                                                             ! Matrix
                   zcof = zfact1 * tmask(ji,jj,jk)
-                  ! lower diagonal
+                  !                                                                 ! 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
+                  !                                                                 ! 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
+                  !                                                                 ! 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
+                  !
+                  !                                                             ! 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
+                  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
-                  ! save masked Prandlt number in zmxld array
-                  zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-
+                  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
-      !  save masked Prandlt number
-      e_pdl(:,:,2:jpkm1) = zmxld(:,:,2:jpkm1)
+      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)      
@@ -516,6 +507,5 @@
       ! 4. Matrix inversion from level 2 (tke prescribed at level 1)
       !!--------------------------------
-      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
-      DO jk = 3, jpkm1
+      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.
@@ -524,7 +514,5 @@
          END DO
       END DO
-
-      ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
-      DO jj = 2, jpjm1
+      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
@@ -538,7 +526,5 @@
          END DO
       END DO
-
-      ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
-      DO jj = 2, jpjm1
+      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)
@@ -552,78 +538,72 @@
          END DO
       END DO
-
-      ! Save the result in en and set minimum value of tke : emin
-      DO jk = 2, jpkm1
+      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), emin ) * tmask(ji,jj,jk)
+               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk)
             END DO
          END DO
       END DO
      
-      ! Modify the value of TKE by an exponential assumption
-      ! en(ji,jj,jk)=en(ji,jj,jk)+en(ji,jj,1)*EXP(-fsdepw(ji,jj,jk)/ zhtau(ji,jj) )
-
-      SELECT CASE( nn_htau )      ! Choice of H value
-      !
-      CASE( 0 )
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zhtau(ji,jj) = 10.
-            END DO
-         END DO
-         !
-      CASE( 1 )
-         DO jj = 2, jpjm1
-            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 )
-         DO jj = 2, jpjm1
-            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
-         DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  en(ji,jj,jk) = en(ji,jj,jk)   &
+      ! 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
-         END DO
-      ELSEIF( nn_etau == 2 ) THEN     !  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
+                     &                   * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
+               END DO
+            END DO
+         ENDIF
+         !
       ENDIF
 
 
-      ! Lateral boundary conditions on ( avt, en )   (sign unchanged)
-      CALL lbc_lnk( en , 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. )
+      ! 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 ( nave )
-         
-      CASE ( 0 )                ! no horizontal average
-
-         ! Vertical eddy viscosity
-
-         DO jk = 2, jpkm1                                 ! Horizontal slab
+      !
+      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.
@@ -635,86 +615,68 @@
             END DO
          END DO
-
-         ! Lateral boundary conditions (avmu,avmv) (U- and V- points, sign unchanged)
-         CALL lbc_lnk( avmu, 'U', 1. )   ;    CALL lbc_lnk( avmv, 'V', 1. )
-         
-      CASE ( 1 )                ! horizontal average
-
-         !                                                ( 1/2  1/2 )
-         ! Eddy viscosity: horizontal average: avmu = 1/4 ( 1    1   )
-         !                      ( 1/2  1 1/2 )            ( 1/2  1/2 )
-         !           avmv = 1/4 ( 1/2  1 1/2 )      
-         
-!! caution vectopt_memory change the solution (last digit of the solver stat)
+         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
-         DO jk = 2, jpkm1                                 ! Horizontal slab
-            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)   &
-                  &                 +.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)
-               END DO
-            END DO
-         END DO
+                  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
-         DO jk = 2, jpkm1                                 ! Horizontal slab
-            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)   &
-                  &              +.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) )  )
-               END DO
-            END DO
-         END DO
+                  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
-
-         ! Lateral boundary conditions (avmu,avmv) (sign unchanged)
-         CALL lbc_lnk( avmu, 'U', 1. )   ;    CALL lbc_lnk( avmv, 'V', 1. )
-
-         ! Vertical eddy diffusivity
-         ! ------------------------------
-         !                                (1 2 1)
-         ! horizontal average  avt = 1/16 (2 4 2)
-         !                                (1 2 1)
-         DO jk = 2, jpkm1                                 ! Horizontal slab
-#  if defined key_vectopt_memory
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  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)
-               END DO
-            END DO
-#  else
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  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)  )
-               END DO
-            END DO
-#  endif
-         END DO
-
+               END DO
+            END DO
+         END DO
+         !
       END SELECT
-
-      ! multiplied by the Prandtl number (npdl>1)
-      ! ----------------------------------------
-      IF( npdl == 1 ) THEN
-         DO jk = 2, jpkm1                                 ! Horizontal slab
+      !
+      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.
@@ -725,30 +687,20 @@
          END DO
       ENDIF
-
-      ! Minimum value on the eddy viscosity
-      ! ----------------------------------------
-      DO jk = 2, jpkm1                                 ! Horizontal slab
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), avmb(jk) ) * umask(ji,jj,jk)
-               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), avmb(jk) ) * vmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-
-      ! Lateral boundary conditions on avt  (sign unchanged)
-      ! ------------------------------=====
-      CALL lbc_lnk( avt, 'W', 1. )
-
-      ! write en in restart file
-      ! ------------------------
-      IF( lrst_oce )   CALL tke_rst( kt, 'WRITE' )
+      !
+      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)
+         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
 
@@ -773,5 +725,4 @@
       !
 # if defined key_vectopt_memory
-      ! caution vectopt_memory change the solution (last digit of the solver stat)
       INTEGER ::   ji, jj, jk   ! dummy loop indices
 # else
@@ -779,8 +730,8 @@
 # endif
       !!
-      NAMELIST/namtke/ ln_rstke, ediff, ediss, ebb, efave, emin, emin0,   &
-                       ri_c, nitke, nmxl, npdl, nave, navb,    &
-                       ln_lsfc, rn_lmin, nn_avtb_2d, nn_etau, nn_htau, rn_efr, &
-                       ln_lc , rn_lc 
+      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 
       !!----------------------------------------------------------------------
 
@@ -791,7 +742,8 @@
 
       ! Compute boost associated with the Richardson critic
-      !     (control values: ri_c = 0.3   ==> eboost=1.25 for npdl=1 or 2)
-      !     (                ri_c = 0.222 ==> eboost=1.                  )
-      eboost = ri_c * ( 2. + ediss / ediff ) / 2.
+      !     (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.
+
 
 
@@ -804,34 +756,36 @@
          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           ediff  = ', ediff
-         WRITE(numout,*) '             Kolmogoroff dissipation coef.  ediss  = ', ediss
-         WRITE(numout,*) '             tke surface input coef.        ebb    = ', ebb
-         WRITE(numout,*) '             tke diffusion coef.            efave  = ', efave
-         WRITE(numout,*) '             minimum value of tke           emin   = ', emin
-         WRITE(numout,*) '             surface minimum value of tke   emin0  = ', emin0
-         WRITE(numout,*) '             number of restart iter loops   nitke  = ', nitke
-         WRITE(numout,*) '             mixing length type             nmxl   = ', nmxl
-         WRITE(numout,*) '             prandl number flag             npdl   = ', npdl
-         WRITE(numout,*) '             horizontal average flag        nave   = ', nave
-         WRITE(numout,*) '             critic Richardson nb           ri_c   = ', ri_c
-         WRITE(numout,*) '                and its associated coeff.   eboost = ', eboost
-         WRITE(numout,*) '             constant background or profile navb   = ', navb
-         WRITE(numout,*) '             flag for compu.of bls as fun. of wind     ln_lsfc    = ', ln_lsfc
-         WRITE(numout,*) '             buoyancy lenght scale surface min value   rn_lmin    = ', rn_lmin
-         WRITE(numout,*) '             horizontal variation for avtb             nn_avtb_2d = ', nn_avtb_2d
-         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 emin 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,*) '             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 nmxl and npdl values
-      IF( nmxl < 0 .OR. nmxl > 3 ) CALL ctl_stop( '          bad flag: nmxl is < 0 or > 3 ' )
-      IF( npdl < 0 .OR. npdl > 1 ) CALL ctl_stop( '          bad flag: npdl is < 0 or > 1 ' )
-      IF( nn_htau < 0    .OR. nn_htau > 2   )   CALL ctl_stop( 'bad flag: nn_htau is < 0 or > 3 ' )
-      IF( rn_lc   < 0.15 .OR. rn_lc   > 0.2 )   CALL ctl_stop( 'bad value: rn_lc must be > 0.15 and < 0.2 ' )
+      ! 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 
@@ -839,7 +793,8 @@
       ! Horizontal average : initialization of weighting arrays 
       ! -------------------
-      
-      SELECT CASE ( nave )
-
+      !
+      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'
@@ -854,22 +809,17 @@
          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)  )
-
-                  evmean(ji,jj,jk) = vmask(ji,jj,jk)                     &
-                  &  / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk)  )
+                  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'
@@ -883,30 +833,21 @@
          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) )  )
+                  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
-
-      CASE DEFAULT
-         WRITE(ctmp1,*) '          bad flag value for nave = ', nave
-         CALL ctl_stop( ctmp1 )
-
+         !
       END SELECT
 
@@ -914,20 +855,17 @@
       ! Background eddy viscosity and diffusivity profil
       ! ------------------------------------------------
-      IF( navb == 0 ) THEN      ! Define avmb, avtb from namelist parameter
+      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
-!!bug  this is not valide neither in scoord
-         IF(ln_sco .AND. lwp) WRITE(numout,cform_war)
-         IF(ln_sco .AND. lwp) WRITE(numout,*) '          avtb profile not valid in sco'
-
          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_avtb_2d == 1 ) THEN      ! decrease avtb in the equatorial band
+      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.
@@ -938,11 +876,12 @@
       ENDIF
 
+!!gm  the increase is only required when using cen2 advection scheme 
+!!gm  for the equatorial upwelling.       ===>  use of TVD in ORCA  so this have to be suppressed
       !   Increase the background in the surface layers
-!!gm  the increase is only required when using cen2 advection scheme 
-!!gm  for the equatorial upwelling.
       avmb(1) = 10.  * avmb(1)      ;      avtb(1) = 10.  * avtb(1)
       avmb(2) = 10.  * avmb(2)      ;      avtb(2) = 10.  * avtb(2)
       avmb(3) =  5.  * avmb(3)      ;      avtb(3) =  5.  * avtb(3)
       avmb(4) =  2.5 * avmb(4)      ;      avtb(4) =  2.5 * avtb(4)
+!!gm end
 
 
@@ -967,8 +906,9 @@
      !!                   ***  ROUTINE ts_rst  ***
      !!                     
-     !! ** Purpose : Read or write filtered free surface arrays in restart file
+     !! ** Purpose :   Read or write TKE file (en) in restart file
      !!
-     !! ** Method  : 
-     !!
+     !! ** 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
@@ -986,13 +926,12 @@
                  &                       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 set by iterative loop'
-              IF( lwp                )   WRITE(numout,*) ' =======             ========='
-              en (:,:,:) = emin * tmask(:,:,:)
-              DO jit = 2, nitke+1
+              IF( lwp                )   WRITE(numout,*) ' ===>>>> : en is set by iterative loop'
+              en (:,:,:) = rn_emin * tmask(:,:,:)
+              DO jit = 2, nn_itke + 1
                  CALL zdf_tke( jit )
               END DO
            ENDIF
         ELSE
-           en(:,:,:) = emin * tmask(:,:,:)      ! no restart: en set to emin
+           en(:,:,:) = rn_emin * tmask(:,:,:)      ! no restart: en set to its minimum value
         ENDIF
      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
