Index: trunk/NEMO/OPA_SRC/ZDF/zdf_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 1491)
+++ trunk/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 1492)
@@ -4,52 +4,44 @@
    !! Ocean physics : define vertical mixing variables
    !!=====================================================================
+   !! history :  1.0  !  2002-06  (G. Madec) Original code
+   !!            3.2  !  2009-07  (G.Madec) addition of avm
    !!----------------------------------------------------------------------
-   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
-   !! $Id$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+
    !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
-   !!   zdf_init    : initialization, namelist read, and parameters control
-   !!----------------------------------------------------------------------
-   !! * Modules used
-   USE par_oce         ! mesh and scale factors
+   USE par_oce         ! ocean parameters
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Share Module variables
-   LOGICAL, PARAMETER, PUBLIC ::    &   !:
 #if defined key_zdfcst   ||   defined key_esopa
-      lk_zdfcst        = .TRUE.         !: constant vertical mixing flag
+   LOGICAL, PARAMETER, PUBLIC ::   lk_zdfcst        = .TRUE.         !: constant vertical mixing flag
 #else
-      lk_zdfcst        = .FALSE.        !: constant vertical mixing flag
+   LOGICAL, PARAMETER, PUBLIC ::   lk_zdfcst        = .FALSE.        !: constant vertical mixing flag
 #endif
-   LOGICAL, PUBLIC ::                & !!! namzdf: vertical diffusion
-      ln_zdfexp        = .FALSE. ,   &  !: explicit vertical diffusion scheme flag
-      ln_zdfevd        = .TRUE.  ,   &  !: convection: enhanced vertical diffusion flag
-      ln_zdfnpc        = .FALSE.        !: convection: non-penetrative convection flag
 
-   INTEGER, PUBLIC ::    & !!: namzdf:  vertical diffusion
-      n_zdfexp = 3    ,  &  !: number of sub-time step (explicit time stepping)
-      n_evdm   = 1          !: =0/1 flag to apply enhanced avm or not
- 
-   REAL(wp), PUBLIC ::   & !!: namzdf   vertical diffusion
-      avm0  = 1.e-4_wp,  &  !: vertical eddy viscosity (m2/s)
-      avt0  = 1.e-5_wp,  &  !: vertical eddy diffusivity (m2/s)
-      avevd = 1._wp         !: vertical eddy coeff. for enhanced vert. diff. (m2/s)
+   !                                          !!* namelist namzdf: vertical diffusion
+   LOGICAL , PUBLIC ::   ln_zdfexp = .FALSE.   !: explicit vertical diffusion scheme flag
+   LOGICAL , PUBLIC ::   ln_zdfevd = .TRUE.    !: convection: enhanced vertical diffusion flag
+   LOGICAL , PUBLIC ::   ln_zdfnpc = .FALSE.   !: convection: non-penetrative convection flag
+   INTEGER , PUBLIC ::   n_zdfexp = 3          !: number of sub-time step (explicit time stepping)
+   INTEGER , PUBLIC ::   n_evdm   = 1          !: =0/1 flag to apply enhanced avm or not
+   REAL(wp), PUBLIC ::   avm0  = 1.e-4_wp      !: vertical eddy viscosity (m2/s)
+   REAL(wp), PUBLIC ::   avt0  = 1.e-5_wp      !: vertical eddy diffusivity (m2/s)
+   REAL(wp), PUBLIC ::   avevd = 1._wp         !: vertical eddy coeff. for enhanced vert. diff. (m2/s)
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
-      avmu,              &  !: vertical viscosity coeff. at uw-, vw-points
-      avmv,              &  !: vertical viscosity coeff. at uw-, vw-points
-      avt ,              &  !: vertical diffusivity coeff. at w-point
-      avt_evd,           &  !: convection: enhanced vertical diffusivity coeff. at w-point
-      avmu_evd              !: convection: enhanced vertical viscosity   coeff. at w-point
+   REAL(wp), PUBLIC, DIMENSION        (jpk) ::   avmb, avtb   !: background profile of avm and avt
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   avmu, avmv   !: vertical viscosity coeff. at uw- & vw-points   [m2/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   avm , avt    !: vertical viscosity & diffusivity coeff. at  w-point   [m2/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   avt_evd      !: enhanced vertical diffusivity coeff. at  w-point   [m2/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   avmu_evd     !: enhanced vertical viscosity   coeff. at uw-point   [m2/s]
 #if defined key_zdftmx
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   av_tide, av_tide_itf   !: Tidal mixing
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   av_tide      !: Tidal mixing
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   av_tide_itf  !: Tidal mixing in the Indonesian Through Flow
 #endif
  
-   REAL(wp), PUBLIC, DIMENSION(jpk) ::   &   !:
-      avmb, avtb            !: background profile of avm and avt
- 
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
+   !! $Id$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!======================================================================
 END MODULE zdf_oce
Index: trunk/NEMO/OPA_SRC/ZDF/zdftke2.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdftke2.F90	(revision 1491)
+++ trunk/NEMO/OPA_SRC/ZDF/zdftke2.F90	(revision 1492)
@@ -5,23 +5,25 @@
    !!                 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_tke2_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 
+   !! 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_zdftke2   ||   defined key_esopa
@@ -29,64 +31,75 @@
    !!   'key_zdftke2'                                  TKE vertical physics
    !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
    !!   zdf_tke2      : update momentum and tracer Kz from a tke scheme
-   !!   zdf_tke2_init : initialization, namelist read, and parameters control
+   !!   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
    !!   tke2_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_tke2    ! routine called in step module
-   PUBLIC   tke2_rst    ! routine called in step module
+   PUBLIC   zdf_tke2   ! routine called in step module
+   PUBLIC   tke2_rst   ! routine called in step module
 
    LOGICAL , PUBLIC, PARAMETER              ::   lk_zdftke2 = .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
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   avm                 !: now mixing coefficient of diffusion for TKE
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   dissl               !: now mixing lenght of dissipation
 
 #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
 #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 tke2_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 tke2_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
@@ -94,5 +107,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
+   !! 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)
@@ -106,53 +119,34 @@
       !!
       !! ** 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_tke2_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)
@@ -163,78 +157,89 @@
       !!              Mellor and Blumberg, JPO 2004
       !!              Axell, JGR, 2002
+      !!              Bruchard OM 2002
       !!----------------------------------------------------------------------
-
-      INTEGER, INTENT(in) ::   kt ! ocean time step
-
-      IF( kt == nit000 )   CALL zdf_tke2_init             ! Initialization (first time-step only)
-      !
-      CALL tke_tke                                        ! now tke (en)
-      CALL tke_mlmc                                       ! now avmu, avmv, avt
-
+      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_tke2
+
 
    SUBROUTINE tke_tke
       !!----------------------------------------------------------------------
-      !! Now Turbulent kinetic energy (output in en)
+      !!                   ***  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] )
       !! ---------------------------------------------------------------------
-      !! Surface boundary condition for TKE
-      !! Resolution of a tridiagonal linear system by a "methode de chasse"
-      !! computation from level 2 to jpkm1  (e(1) computed and e(jpk)=0 ).
-      !! avm represents the turbulent coefficient of the TKE
-      !!----------------------------------------------------------------------
-      USE oce,     zwd     =>   ua   ! use ua as workspace
-      USE oce,     zdiag1  =>   va   ! use va as workspace
-      USE oce,     zdiag2  =>   ta   ! use ta as workspace
-      USE oce,     ztkelc  =>   sa   ! use sa as workspace
-      !
-      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
-      REAL(wp) ::   zbbrau, zesurf, zesh2           ! temporary scalars
-      REAL(wp) ::   zfact1, ztx2, zdku              !    -         -
-      REAL(wp) ::   zfact2, zty2, zdkv              !    -         -
-      REAL(wp) ::   zfact3, zcoef, zcof             !    -         -
-      REAL(wp) ::   zus, zwlc, zind                 !    -         -
+      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
       !!--------------------------------------------------------------------
-      !                                            ! 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       * rn_ediss
-
-      ! Surface boundary condition on tke
-      ! -------------------------------------------------
-      ! en(1)   = rn_ebb sqrt(utau^2+vtau^2) / rau0  (min value rn_emin0)
-!CDIR NOVERRCHK
-      DO jj = 2, jpjm1
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  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)
+            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
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Langmuir circulation source term
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_lc ) THEN
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Bottom boundary condition on tke
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!!gm to be added soon
+      !
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke
+         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
          !
-         ! Computation of total energy produce by LC : cumulative sum over jk
+         !                        !* 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(:,:)
+         !                        !* finite Langmuir Circulation depth
+         imlc(:,:) = mbathy(:,:)         ! Initialization to the number of w ocean point 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)
+            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
@@ -242,30 +247,31 @@
             END DO
          END DO
-         !
-         ! finite LC depth
-         DO jj = 1, jpj
+         !                               ! 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
-         !
-         ! TKE Langmuir circulation source term added to en
-!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
+!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)
+                  !                                           ! 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
@@ -273,47 +279,44 @@
          !
       ENDIF
-
-
-      ! Shear production at uw- and vw-points (energy conserving form)
-      DO jk = 2, jpkm1
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  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.
-               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) )
-            ENDDO
-         ENDDO
-      ENDDO
-
-      ! Lateral boundary conditions (avmu,avmv) (sign unchanged)
-      CALL lbc_lnk( avmu, 'U', 1. )   ;    CALL lbc_lnk( avmv, 'V', 1. )
-
-      ! 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 ).
-      ! zwd : diagonal zdiag1 : upper diagonal zdiag2 : lower diagonal
-      ! Warning : after this step, en : right hand side of the matrix
-      DO jk = 2, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               !                                                             ! shear prod. at w-point weightened by mask
+               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) )    
-               !                                                             ! Matrix
-               zcof = zfact1 * tmask(ji,jj,jk)
-               !                                                                  ! lower diagonal
-               zdiag2(ji,jj,jk) = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &
-                  &                    / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
-               !                                                                  ! upper diagonal
-               zdiag1(ji,jj,jk) = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &
-                  &                    / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
-               !                                                                  ! diagonal
-               zwd(ji,jj,jk) = 1. - zdiag2(ji,jj,jk) - zdiag1(ji,jj,jk) + zfact2 * dissl(ji,jj,jk) * tmask(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
+               !                                   ! 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)
@@ -321,10 +324,9 @@
          END DO
       END DO
-      ! Matrix inversion from level 2 (tke prescribed at level 1)
-      !!--------------------------------
+      !                          !* 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) - zdiag2(ji,jj,jk) * zdiag1(ji,jj,jk-1) / zwd(ji,jj,jk-1)
+               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
@@ -332,5 +334,5 @@
       DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
          DO ji = fs_2, fs_jpim1    ! vector opt.
-            zdiag2(ji,jj,2) = en(ji,jj,2) - zdiag2(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
+            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
@@ -338,5 +340,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
-               zdiag2(ji,jj,jk) = en(ji,jj,jk) - zdiag2(ji,jj,jk) / zwd(ji,jj,jk-1) *zdiag2(ji,jj,jk-1)
+               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
@@ -344,5 +346,5 @@
       DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
          DO ji = fs_2, fs_jpim1    ! vector opt.
-            en(ji,jj,jpkm1) = zdiag2(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
+            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
          END DO
       END DO
@@ -350,5 +352,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
-               en(ji,jj,jk) = ( zdiag2(ji,jj,jk) - zdiag1(ji,jj,jk) * en(ji,jj,jk+1) ) / zwd(ji,jj,jk)
+               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
@@ -362,79 +364,85 @@
       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
-         !
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  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
-      ! Lateral boundary conditions (sign unchanged)
-      CALL lbc_lnk( en, 'W', 1. )
-
+      !
+      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
+      !
    END SUBROUTINE tke_tke
 
-   SUBROUTINE tke_mlmc
+
+   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) ::   zpdl, zri, zsqen                !    -         -
-      REAL(wp) ::   zemxl, zemlm, zemlp             !    -         -
+      !!
+      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)
-      ! ---------------------
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  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
@@ -466,10 +474,12 @@
          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 )
       !
@@ -545,21 +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  Vertical eddy viscosity on tke (put in zmxlm) and first estimate of avt !
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<!
-      ! Surface boundary condition on avt and avm : jk = 1
-      ! -------------------------------------------------
-      ! avt(1) = avmb(1) and avt(jpk) = 0.
-      ! avm(1) = avmb(1) and avm(jpk) = 0.
-      !
-!CDIR NOVERRCHK
-      DO jk = 1, 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
@@ -568,15 +572,13 @@
                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)
+               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
-
-      ! Lateral boundary conditions (sign unchanged)
-      CALL lbc_lnk( avm, 'W', 1. )
-
-      DO jk = 2, jpkm1                                 ! only vertical eddy viscosity
+      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.
@@ -587,38 +589,29 @@
       END DO
       CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
-
-
-      ! Prandtl number
-      ! ----------------------------------------
-      IF( nn_pdl == 1 ) THEN
+      !
+      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))
-                  zri = MAX( rn2b(ji,jj,jk), 0. ) * avm(ji,jj,jk)/ (zcoef * (zdku + zdkv + 1.e-20) ) ! local Richardson number
-# if defined key_cfg_1d
+                  !                                          ! 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
-                  zpdl = 1.0                                                         ! Prandtl number
-                  IF( zri >= 0.2 )   zpdl = 0.2 / zri
-                  zpdl = MAX( 0.1, zpdl )
-                  avt(ji,jj,jk)   = MAX( zpdl * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
-                  zmpdl(ji,jj,jk) = zpdl * tmask(ji,jj,jk)
               END DO
             END DO
          END DO
       ENDIF
-
-# 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
-
       CALL lbc_lnk( avt, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged)
 
@@ -629,29 +622,22 @@
       ENDIF
       !
-   END SUBROUTINE tke_mlmc
-
-   SUBROUTINE zdf_tke2_init
+   END SUBROUTINE tke_avn
+
+
+   SUBROUTINE tke_init
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_tke2_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,   &
@@ -661,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_tke2_init : tke turbulent closure scheme'
-         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) 'zdf_tke2 : 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
@@ -685,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
@@ -705,26 +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 
-      ! -------------------
-      !
-      IF(lwp) WRITE(numout,*) '          no horizontal average on avt, avmu, avmv'
-
-      ! 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
@@ -732,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.
@@ -744,6 +709,25 @@
       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)
@@ -753,10 +737,8 @@
       END DO
       dissl(:,:,:) = 1.e-12
-
-      ! read or initialize turbulent kinetic energy ( en )
-      ! -------------------------------------------------
+      !                               !* read or initialize all required files 
       CALL tke2_rst( nit000, 'READ' )
       !
-   END SUBROUTINE zdf_tke2_init
+   END SUBROUTINE tke_init
 
 
@@ -796,11 +778,11 @@
                  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_mlmc                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
+              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_mlmc                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
+              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
               DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke2( jit )   ;   END DO
            ENDIF
@@ -838,4 +820,8 @@
       WRITE(*,*) 'zdf_tke2: You should not have seen this print! error?', kt
    END SUBROUTINE zdf_tke2
+   SUBROUTINE tke2_rst( kt, cdrw )
+     CHARACTER(len=*) ::   cdrw
+     WRITE(*,*) 'tke2_rst: You should not have seen this print! error?', kt, cdwr
+   END SUBROUTINE tke2_rst
 #endif
 
Index: trunk/NEMO/OPA_SRC/step.F90
===================================================================
--- trunk/NEMO/OPA_SRC/step.F90	(revision 1491)
+++ trunk/NEMO/OPA_SRC/step.F90	(revision 1492)
@@ -5,8 +5,8 @@
    !!======================================================================
    !! History :  OPA  !  1991-03  (G. Madec)  Original code
-   !!                 !  1991-11  (G. Madec)
-   !!                 !  1992-06  (M. Imbard)  add a first output record
-   !!                 !  1996-04  (G. Madec)  introduction of dynspg
-   !!                 !  1996-04  (M.A. Foujols)  introduction of passive tracer
+   !!             -   !  1991-11  (G. Madec)
+   !!             -   !  1992-06  (M. Imbard)  add a first output record
+   !!             -   !  1996-04  (G. Madec)  introduction of dynspg
+   !!             -   !  1996-04  (M.A. Foujols)  introduction of passive tracer
    !!            8.0  !  1997-06  (G. Madec)  new architecture of call
    !!            8.2  !  1997-06  (G. Madec, M. Imbard, G. Roullet)  free surface
@@ -15,9 +15,10 @@
    !!   NEMO     1.0  !  2002-06  (G. Madec)  free form, suppress macro-tasking
    !!             -   !  2004-08  (C. Talandier) New trends organization
-   !!                 !  2005-01  (C. Ethe) Add the KPP closure scheme
-   !!                 !  2005-11  (G. Madec)  Reorganisation of tra and dyn calls
-   !!                 !  2006-01  (L. Debreu, C. Mazauric)  Agrif implementation
-   !!                 !  2006-07  (S. Masson)  restart using iom
+   !!             -   !  2005-01  (C. Ethe) Add the KPP closure scheme
+   !!             -   !  2005-11  (G. Madec)  Reorganisation of tra and dyn calls
+   !!             -   !  2006-01  (L. Debreu, C. Mazauric)  Agrif implementation
+   !!             -   !  2006-07  (S. Masson)  restart using iom
    !!            3.2  !  2009-02  (G. Madec, R. Benshila)  reintroduicing z*-coordinate
+   !!             -   !  2009-06  (S. Masson, G. Madec)  TKE2 restart compatible with key_cpl
    !!----------------------------------------------------------------------
 
@@ -123,5 +124,5 @@
    PRIVATE
 
-   PUBLIC stp            ! called by opa.F90
+   PUBLIC   stp        ! called by opa.F90
 
    !! * Substitutions
@@ -129,5 +130,5 @@
 #  include "zdfddm_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
@@ -179,92 +180,76 @@
       ! Update data, open boundaries, surface boundary condition (including sea-ice)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
       IF( lk_dtatem  )   CALL dta_tem( kstp )         ! update 3D temperature data
       IF( lk_dtasal  )   CALL dta_sal( kstp )         ! update 3D salinity data
-
                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
-
       IF( lk_obc     )   CALL obc_dta( kstp )         ! update dynamic and tracer data at open boundaries
       IF( lk_obc     )   CALL obc_rad( kstp )         ! compute phase velocities at open boundaries
-
       IF( lk_bdy     )   CALL bdy_dta( kstp )         ! update dynamic and tracer data at unstructured open boundary
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      !                Ocean dynamics : ssh, wn, hdiv, rot                   !
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                       CALL div_cur( kstp )                 ! Horizontal divergence & Relative vorticity
-      IF( n_cla == 1 ) CALL div_cla( kstp )                 ! Cross Land Advection (Update Hor. divergence)
-                       CALL ssh_wzv( kstp )                 ! after ssh & vertical velocity
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Ocean physics update
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
-      !-----------------------------------------------------------------------
-
-                       CALL bn2( tb, sb, rn2b )              ! before Brunt-Vaisala frequency
-                       CALL bn2( tn, sn, rn2  )              ! now    Brunt-Vaisala frequency
-
-      !-----------------------------------------------------------------------
-      !  VERTICAL PHYSICS
-      !-----------------------------------------------------------------------
-      !                                                     ! Vertical eddy viscosity and diffusivity coefficients
-      IF( lk_zdfric )   CALL zdf_ric( kstp )                       ! Richardson number dependent Kz
-
-      IF( lk_zdftke )   CALL zdf_tke ( kstp )                      ! TKE  closure scheme for Kz
-      IF( lk_zdftke2)   CALL zdf_tke2( kstp )                      ! TKE2 closure scheme for Kz
-
-      IF( lk_zdfkpp )   CALL zdf_kpp( kstp )                       ! KPP closure scheme for Kz
-
-      IF( lk_zdfcst )   THEN                                       ! Constant Kz (reset avt, avm[uv] to the background value)
+      !  Ocean dynamics : ssh, wn, hdiv, rot                                 !
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         CALL div_cur( kstp )         ! Horizontal divergence & Relative vorticity
+      IF( n_cla == 1 )   CALL div_cla( kstp )         ! Cross Land Advection (Update Hor. divergence)
+                         CALL ssh_wzv( kstp )         ! after ssh & vertical velocity
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Ocean physics update                (ua, va, ta, sa used as workspace)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         CALL bn2( tb, sb, rn2b )     ! before Brunt-Vaisala frequency
+                         CALL bn2( tn, sn, rn2  )     ! now    Brunt-Vaisala frequency
+      !
+      !  VERTICAL PHYSICS   
+      !                                               ! Vertical eddy viscosity and diffusivity coefficients
+      IF( lk_zdfric  )   CALL zdf_ric( kstp )         ! Richardson number dependent Kz
+      IF( lk_zdftke  )   CALL zdf_tke ( kstp )        ! TKE  closure scheme for Kz
+      IF( lk_zdftke2 )   CALL zdf_tke2( kstp )        ! TKE2 closure scheme for Kz
+      IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )         ! KPP closure scheme for Kz
+      IF( lk_zdfcst  )   THEN                         ! Constant Kz (reset avt, avm[uv] to the background value)
          avt (:,:,:) = avt0 * tmask(:,:,:)
          avmu(:,:,:) = avm0 * umask(:,:,:)
          avmv(:,:,:) = avm0 * vmask(:,:,:)
       ENDIF
-
-      IF( ln_rnf_mouth ) THEN                                ! increase diffusivity at rivers mouths
+      IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
          DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
       ENDIF
-
-      IF( ln_zdfevd )   CALL zdf_evd( kstp )                 ! enhanced vertical eddy diffusivity
-
-      IF( lk_zdftmx )   CALL zdf_tmx( kstp )                 ! tidal vertical mixing
+      IF( ln_zdfevd  )   CALL zdf_evd( kstp )         ! enhanced vertical eddy diffusivity
+
+      IF( lk_zdftmx  )   CALL zdf_tmx( kstp )         ! tidal vertical mixing
 
       IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   &
-         &              CALL zdf_ddm( kstp )                 ! double diffusive mixing
-
-
-                        CALL zdf_bfr( kstp )                 ! bottom friction
-
-                        CALL zdf_mxl( kstp )                 ! mixed layer depth
-
-      IF( lrst_oce .AND. lk_zdftke2 )   &                    ! write tke2 information in the restart file
-         &               CALL tke2_rst( kstp, 'WRITE' )
-
-      !-----------------------------------------------------------------------
-      !  LATERAL PHYSICS
-      !-----------------------------------------------------------------------
-      IF( lk_ldfslp )   THEN 
-                           CALL eos( tb, sb, rhd, rhop )          ! before in situ density
-         IF( ln_zps )      CALL zps_hde( kstp, tb, sb, rhd,  &    ! Partial steps: before horizontal gradient
-            &                                gtu, gsu, gru,  &    ! of t, s, rd at the last ocean level
-            &                                gtv, gsv, grv )
-                               CALL ldf_slp( kstp, rhd, rn2b )      ! before slope of the lateral mixing
+         &               CALL zdf_ddm( kstp )         ! double diffusive mixing
+
+                         CALL zdf_bfr( kstp )         ! bottom friction
+
+                         CALL zdf_mxl( kstp )         ! mixed layer depth
+
+                                                      ! write tke2 information in the restart file
+      IF( lrst_oce .AND. lk_zdftke2 )   CALL tke2_rst( kstp, 'WRITE' )
+      !
+      !  LATERAL  PHYSICS 
+      !
+      IF( lk_ldfslp ) THEN                            ! slope of lateral mixing
+                         CALL eos( tb, sb, rhd, rhop )          ! before in situ density
+         IF( ln_zps )    CALL zps_hde( kstp, tb, sb, rhd,  &    ! Partial steps: before horizontal gradient
+            &                              gtu, gsu, gru,  &    ! of t, s, rd at the last ocean level
+            &                              gtv, gsv, grv )
+                         CALL ldf_slp( kstp, rhd, rn2b )        ! before slope of the lateral mixing
       ENDIF
 #if defined key_traldf_c2d
-      IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )                 ! eddy induced velocity coefficient
+      IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )      ! eddy induced velocity coefficient
 #  endif
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! diagnostics and outputs
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL dia_wri( kstp, indic )         ! ocean model: outputs
-      IF( lk_floats  )   CALL flo_stp( kstp )                ! drifting Floats
-      IF( lk_diaspr  )   CALL dia_spr( kstp )                ! Surface pressure diagnostics
-      IF( lk_diahth  )   CALL dia_hth( kstp )                ! Thermocline depth (20 degres isotherm depth)
-      IF( lk_diagap  )   CALL dia_gap( kstp )                ! basin averaged diagnostics
-      IF( lk_diahdy  )   CALL dia_hdy( kstp )                ! dynamical heigh diagnostics
-      IF( lk_diafwb  )   CALL dia_fwb( kstp )                ! Fresh water budget diagnostics
-      IF( ln_diaptr  )   CALL dia_ptr( kstp )                ! Poleward TRansports diagnostics
+      ! diagnostics and outputs             (ua, va, ta, sa used as workspace)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         CALL dia_wri( kstp, indic )  ! ocean model: outputs
+      IF( lk_floats  )   CALL flo_stp( kstp )         ! drifting Floats
+      IF( lk_diaspr  )   CALL dia_spr( kstp )         ! Surface pressure diagnostics
+      IF( lk_diahth  )   CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth)
+      IF( lk_diagap  )   CALL dia_gap( kstp )         ! basin averaged diagnostics
+      IF( lk_diahdy  )   CALL dia_hdy( kstp )         ! dynamical heigh diagnostics
+      IF( lk_diafwb  )   CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
+      IF( ln_diaptr  )   CALL dia_ptr( kstp )         ! Poleward TRansports diagnostics
 
 #if defined key_top
@@ -272,14 +257,10 @@
       ! Passive Tracer Model
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
-      !-----------------------------------------------------------------------
-                             CALL trc_stp( kstp )            ! time-stepping
-#endif
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Active tracers
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      ! N.B. ua, va arrays are used as workspace in this section
-      !-----------------------------------------------------------------------
+                         CALL trc_stp( kstp )         ! time-stepping
+#endif
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Active tracers                              (ua, va used as workspace)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                              ta(:,:,:) = 0.e0               ! set tracer trends to zero
                              sa(:,:,:) = 0.e0
@@ -300,67 +281,64 @@
                              CALL tra_zdf    ( kstp )       ! vertical mixing and after tracer fields
 
-      IF( .NOT. ln_dynhpg_imp  ) THEN                       ! centered hpg (default case)
-                               CALL eos( tn, sn, rhd, rhop )          ! now (swap=before) in situ density for dynhpg module
-         IF( ln_zps    )       CALL zps_hde( kstp, tn, sn, rhd,   &   ! Partial steps: now horizontal gradient
-            &                                     gtu, gsu, gru,  &   ! of t, s, rd at the bottom ocean level
-            &                                     gtv, gsv, grv )
+      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg (time stepping then eos)
+         IF( ln_zdfnpc   )   CALL tra_npc    ( kstp )            ! update after fields by non-penetrative convection
+                             CALL tra_nxt    ( kstp )            ! tracer fields at next time step
+                             CALL eos( ta, sa, rhd, rhop )       ! Time-filtered in situ density for hpg computation
+         IF( ln_zps      )   CALL zps_hde( kstp, ta, sa, rhd,   &   ! Partial steps: time filtered hor. derivative
+            &                                   gtu, gsu, gru,  &   ! of t, s, rd at the bottom ocean level
+            &                                   gtv, gsv, grv )
+         
+      ELSE                                                  ! centered hpg  (eos then time stepping)
+                             CALL eos( tn, sn, rhd, rhop )       ! now in situ density for hpg computation
+         IF( ln_zps      )   CALL zps_hde( kstp, tn, sn, rhd,   &   ! Partial steps: now horizontal derivative
+            &                                   gtu, gsu, gru,  &   ! of t, s, rd at the bottom ocean level
+            &                                   gtv, gsv, grv )
+         IF( ln_zdfnpc   )   CALL tra_npc    ( kstp )       ! update after fields by non-penetrative convection
+                             CALL tra_nxt    ( kstp )       ! tracer fields at next time step
       ENDIF 
-      IF( ln_zdfnpc      )   CALL tra_npc    ( kstp )       ! update after fields by non-penetrative convection
-                             CALL tra_nxt    ( kstp )       ! tracer fields at next time step
-      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg
-                               CALL eos( ta, sa, rhd, rhop )          ! Time-filtered in situ density used in dynhpg module
-         IF( ln_zps    )       CALL zps_hde( kstp, ta, sa, rhd,   &   ! Partial steps: time filtered hor. gradient
-            &                                     gtu, gsu, gru,  &   ! of t, s, rd at the bottom ocean level
-            &                                     gtv, gsv, grv )
-      ENDIF
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Dynamics
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      ! N.B. ta, sa arrays are used as workspace in this section 
-      !-----------------------------------------------------------------------
-                               ua(:,:,:) = 0.e0               ! set dynamics trends to zero
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Dynamics                                    (ta, sa used as workspace)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                               ua(:,:,:) = 0.e0             ! set dynamics trends to zero
                                va(:,:,:) = 0.e0
 
-                               CALL dyn_adv( kstp )           ! advection (vector or flux form)
-                               CALL dyn_vor( kstp )           ! vorticity term including Coriolis
-                               CALL dyn_ldf( kstp )           ! lateral mixing
-#if defined key_agrif
-      IF(.NOT. Agrif_Root())   CALL Agrif_Sponge_dyn          ! momemtum sponge
-#endif
-                               CALL dyn_hpg( kstp )           ! horizontal gradient of Hydrostatic pressure
-                               CALL dyn_zdf( kstp )           ! vertical diffusion
+                               CALL dyn_adv( kstp )         ! advection (vector or flux form)
+                               CALL dyn_vor( kstp )         ! vorticity term including Coriolis
+                               CALL dyn_ldf( kstp )         ! lateral mixing
+#if defined key_agrif
+      IF(.NOT. Agrif_Root())   CALL Agrif_Sponge_dyn        ! momemtum sponge
+#endif
+                               CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
+                               CALL dyn_zdf( kstp )         ! vertical diffusion
       IF( lk_dynspg_rl ) THEN
-         IF( lk_obc    )       CALL obc_spg( kstp )           ! surface pressure gradient at open boundaries
+         IF( lk_obc    )       CALL obc_spg( kstp )         ! surface pressure gradient at open boundaries
       ENDIF
                                indic=0
-                               CALL dyn_spg( kstp, indic )    ! surface pressure gradient
-                               CALL dyn_nxt( kstp )           ! lateral velocity at next time step
-
-                               CALL ssh_nxt( kstp )           ! sea surface height at next time step
+                               CALL dyn_spg( kstp, indic )  ! surface pressure gradient
+                               CALL dyn_nxt( kstp )         ! lateral velocity at next time step
+
+                               CALL ssh_nxt( kstp )         ! sea surface height at next time step
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Control and restarts
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                                 CALL stp_ctl( kstp, indic )
-      IF( indic < 0          )   THEN
-                                 CALL ctl_stop( 'step: indic < 0' )
-                                 CALL dia_wri_state( 'output.abort', kstp )
-      ENDIF
-      IF( kstp == nit000     )   CALL iom_close( numror )             ! close input  ocean restart file
-      IF( lrst_oce           )   CALL rst_write    ( kstp )           ! write output ocean restart file
-      IF( lk_obc             )   CALL obc_rst_write( kstp )           ! write open boundary restart file
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Trends
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
-      !-----------------------------------------------------------------------
-
-      IF( nstop == 0 ) THEN                                 ! Diagnostics:
-         IF( lk_trddyn  )   CALL trd_dwr( kstp )                 ! trends: dynamics 
-         IF( lk_trdtra  )   CALL trd_twr( kstp )                 ! trends: active tracers
-         IF( lk_trdmld  )   CALL trd_mld( kstp )                 ! trends: Mixed-layer 
-         IF( lk_trdvor  )   CALL trd_vor( kstp )                 ! trends: vorticity budget
+                               CALL stp_ctl( kstp, indic )
+      IF( indic < 0        )   THEN
+                               CALL ctl_stop( 'step: indic < 0' )
+                               CALL dia_wri_state( 'output.abort', kstp )
+      ENDIF
+      IF( kstp == nit000   )   CALL iom_close( numror )     ! close input  ocean restart file
+      IF( lrst_oce         )   CALL rst_write    ( kstp )   ! write output ocean restart file
+      IF( lk_obc           )   CALL obc_rst_write( kstp )   ! write open boundary restart file
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Trends                              (ua, va, ta, sa used as workspace)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( nstop == 0 ) THEN                         
+         IF( lk_trddyn     )   CALL trd_dwr( kstp )         ! trends: dynamics 
+         IF( lk_trdtra     )   CALL trd_twr( kstp )         ! trends: active tracers
+         IF( lk_trdmld     )   CALL trd_mld( kstp )         ! trends: Mixed-layer 
+         IF( lk_trdvor     )   CALL trd_vor( kstp )         ! trends: vorticity budget
       ENDIF
 
@@ -368,6 +346,5 @@
       ! Coupled mode
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-      IF( lk_cpl )   CALL sbc_cpl_snd( kstp )                 ! coupled mode : field exchanges
+      IF( lk_cpl           )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
       !
       !
