ecearth3/sources/nemo-3.6/NEMO/OPA_SRC/ZDF/zdftke.F90

MODULE zdftke
   !!======================================================================
   !!                       ***  MODULE  zdftke  ***
   !! Ocean physics:  vertical mixing coefficient computed from the tke 
   !!                 turbulent closure parameterization
   !!=====================================================================
   !! History :  OPA  !  1991-03  (b. blanke)  Original code
   !!            7.0  !  1991-11  (G. Madec)   bug fix
   !!            7.1  !  1992-10  (G. Madec)   new mixing length and eav
   !!            7.2  !  1993-03  (M. Guyon)   symetrical conditions
   !!            7.3  !  1994-08  (G. Madec, M. Imbard)  nn_pdl flag
   !!            7.5  !  1996-01  (G. Madec)   s-coordinates
   !!            8.0  !  1997-07  (G. Madec)   lbc
   !!            8.1  !  1999-01  (E. Stretta) new option for the mixing length
   !!  NEMO      1.0  !  2002-06  (G. Madec) add tke_init routine
   !!             -   !  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
   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
   !!----------------------------------------------------------------------
#if defined key_zdftke   ||   defined key_esopa
   !!----------------------------------------------------------------------
   !!   'key_zdftke'                                   TKE vertical physics
   !!----------------------------------------------------------------------
   !!   zdf_tke       : update momentum and tracer Kz from a tke scheme
   !!   tke_tke       : tke time stepping: update tke at now time step (en)
   !!   tke_avn       : compute mixing length scale and deduce avm and avt
   !!   zdf_tke_init  : initialization, namelist read, and parameters control
   !!   tke_rst       : read/write tke restart in ocean restart file
   !!----------------------------------------------------------------------
   USE oce            ! ocean: dynamics and active tracers variables
   USE phycst         ! physical constants
   USE dom_oce        ! domain: ocean
   USE domvvl         ! domain: variable volume layer
   USE sbc_oce        ! surface boundary condition: ocean
   USE zdf_oce        ! vertical physics: ocean variables
   USE zdfmxl         ! vertical physics: mixed layer
   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 lib_mpp        ! MPP library
   USE wrk_nemo       ! work arrays
   USE timing         ! Timing
   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
#if defined key_agrif
   USE agrif_opa_interp
   USE agrif_opa_update
#endif



   IMPLICIT NONE
   PRIVATE

   PUBLIC   zdf_tke        ! routine called in step module
   PUBLIC   zdf_tke_init   ! routine called in opa module
   PUBLIC   tke_rst        ! routine called in step module

   LOGICAL , PUBLIC, PARAMETER ::   lk_zdftke = .TRUE.  !: TKE vertical mixing flag

   !                      !!** Namelist  namzdf_tke  **
   LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not
   INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3)
   REAL(wp) ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m]
   INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1)
   REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
   REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation 
   REAL(wp) ::   rn_ebb    ! coefficient of the surface input of tke
   REAL(wp) ::   rn_emin   ! minimum value of tke           [m2/s2]
   REAL(wp) ::   rn_emin0  ! surface minimum value of tke   [m2/s2]
   REAL(wp) ::   rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
   INTEGER  ::   nn_etau   ! type of depth penetration of surface tke (=0/1/2/3)
   INTEGER  ::   nn_htau   ! type of tke profile of penetration (=0/1)
   REAL(wp) ::   rn_efr    ! fraction of TKE surface value which penetrates in the ocean
   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
   REAL(wp) ::   rn_lc     ! coef to compute vertical velocity of Langmuir cells

   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
   REAL(wp) ::   rmxl_min  ! minimum mixing length value (deduced from rn_ediff and rn_emin values)  [m]
   REAL(wp) ::   rhftau_add = 1.e-3_wp     ! add offset   applied to HF part of taum  (nn_etau=3)
   REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3)

   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau)
   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation
#if defined key_c1d
   !                                                                        !!** 1D cfg only  **   ('key_c1d')
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_dis, e_mix   !: dissipation and mixing turbulent lengh scales
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_pdl, e_ric   !: prandl and local Richardson numbers
#endif

   !! * Substitutions
#  include "domzgr_substitute.h90"
#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
   !! $Id: zdftke.F90 4990 2014-12-15 16:42:49Z timgraham $
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   INTEGER FUNCTION zdf_tke_alloc()
      !!----------------------------------------------------------------------
      !!                ***  FUNCTION zdf_tke_alloc  ***
      !!----------------------------------------------------------------------
      ALLOCATE(                                                                    &
#if defined key_c1d
         &      e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) ,                          &
         &      e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) ,                          &
#endif
         &      htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) , STAT= zdf_tke_alloc      )
         !
      IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc )
      IF( zdf_tke_alloc /= 0 )   CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays')
      !
   END FUNCTION zdf_tke_alloc


   SUBROUTINE zdf_tke( kt )
      !!----------------------------------------------------------------------
      !!                   ***  ROUTINE zdf_tke  ***
      !!
      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
      !!              coefficients using a turbulent closure scheme (TKE).
      !!
      !! ** Method  :   The time evolution of the turbulent kinetic energy (tke)
      !!              is computed from a prognostic equation :
      !!         d(en)/dt = avm (d(u)/dz)**2             ! shear production
      !!                  + d( avm d(en)/dz )/dz         ! diffusion of tke
      !!                  + avt N^2                      ! stratif. destruc.
      !!                  - rn_ediss / emxl en**(2/3)    ! Kolmogoroff dissipation
      !!      with the boundary conditions:
      !!         surface: en = max( rn_emin0, rn_ebb * taum )
      !!         bottom : en = rn_emin
      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 
      !!
      !!        The now Turbulent kinetic energy is computed using the following 
      !!      time stepping: implicit for vertical diffusion term, linearized semi
      !!      implicit for kolmogoroff dissipation term, and explicit forward for 
      !!      both buoyancy and shear production terms. Therefore a tridiagonal 
      !!      linear system is solved. Note that buoyancy and shear terms are
      !!      discretized in a energy conserving form (Bruchard 2002).
      !!
      !!        The dissipative and mixing length scale are computed from en and
      !!      the stratification (see tke_avn)
      !!
      !!        The now vertical eddy vicosity and diffusivity coefficients are
      !!      given by: 
      !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
      !!              avt = max( avmb, pdl * avm                 )  
      !!              eav = max( avmb, avm )
      !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1 
      !!
      !! ** Action  :   compute en (now turbulent kinetic energy)
      !!                update avt, avmu, avmv (before vertical eddy coef.)
      !!
      !! References : Gaspar et al., JGR, 1990,
      !!              Blanke and Delecluse, JPO, 1991
      !!              Mellor and Blumberg, JPO 2004
      !!              Axell, JGR, 2002
      !!              Bruchard OM 2002
      !!----------------------------------------------------------------------
      INTEGER, INTENT(in) ::   kt   ! ocean time step
      !!----------------------------------------------------------------------
      !
      IF( kt /= nit000 ) THEN   ! restore before value to compute tke
         avt (:,:,:) = avt_k (:,:,:) 
         avm (:,:,:) = avm_k (:,:,:) 
         avmu(:,:,:) = avmu_k(:,:,:) 
         avmv(:,:,:) = avmv_k(:,:,:) 
      ENDIF 
      !
      CALL tke_tke      ! now tke (en)
      !
      CALL tke_avn      ! now avt, avm, avmu, avmv
      !
      avt_k (:,:,:) = avt (:,:,:) 
      avm_k (:,:,:) = avm (:,:,:) 
      avmu_k(:,:,:) = avmu(:,:,:) 
      avmv_k(:,:,:) = avmv(:,:,:) 
      !
#if defined key_agrif
      ! Update child grid f => parent grid 
      IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only
#endif      
     ! 
   END SUBROUTINE zdf_tke


   SUBROUTINE tke_tke
      !!----------------------------------------------------------------------
      !!                   ***  ROUTINE tke_tke  ***
      !!
      !! ** Purpose :   Compute the now Turbulente Kinetic Energy (TKE)
      !!
      !! ** Method  : - TKE surface boundary condition
      !!              - source term due to Langmuir cells (Axell JGR 2002) (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] )
      !! ---------------------------------------------------------------------
      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
!!bfr      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! temporary scalar
!!bfr      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          ! temporary scalar
      REAL(wp) ::   zrhoa  = 1.22                   ! Air density kg/m3
      REAL(wp) ::   zcdrag = 1.5e-3                 ! drag coefficient
      REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars
      REAL(wp) ::   zfact1, zfact2, zfact3          !    -         -
      REAL(wp) ::   ztx2  , zty2  , zcof            !    -         -
      REAL(wp) ::   ztau  , zdif                    !    -         -
      REAL(wp) ::   zus   , zwlc  , zind            !    -         -
      REAL(wp) ::   zzd_up, zzd_lw                  !    -         -
!!bfr      REAL(wp) ::   zebot                           !    -         -
      INTEGER , POINTER, DIMENSION(:,:  ) :: imlc
      REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc
      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw
      !!--------------------------------------------------------------------
      !
      IF( nn_timing == 1 )  CALL timing_start('tke_tke')
      !
      CALL wrk_alloc( jpi,jpj, imlc )    ! integer
      CALL wrk_alloc( jpi,jpj, zhlc ) 
      CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
      !
      zbbrau = rn_ebb / rau0       ! Local constant initialisation
      zfact1 = -.5_wp * rdt 
      zfact2 = 1.5_wp * rdt * rn_ediss
      zfact3 = 0.5_wp       * rn_ediss
      !
      !
      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      !                     !  Surface boundary condition on tke
      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      IF ( ln_isfcav ) THEN
         DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin
            DO ji = fs_2, fs_jpim1   ! vector opt.
               en(ji,jj,mikt(ji,jj))=rn_emin * tmask(ji,jj,1)
            END DO
         END DO
      END IF
      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
         DO ji = fs_2, fs_jpim1   ! vector opt.
            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
         END DO
      END DO
      
!!bfr   - start commented area
      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      !                     !  Bottom boundary condition on tke
      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      !
      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      ! Tests to date have found the bottom boundary condition on tke to have very little effect.
      ! The condition is coded here for completion but commented out until there is proof that the
      ! computational cost is justified
      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      !                     en(bot)   = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
!CDIR NOVERRCHK
!!    DO jj = 2, jpjm1
!CDIR NOVERRCHK
!!       DO ji = fs_2, fs_jpim1   ! vector opt.
!!          ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + &
!!                 bfrua(ji  ,jj) * ub(ji  ,jj,mbku(ji  ,jj) )
!!          zty2 = bfrva(ji,jj  ) * vb(ji,jj  ,mbkv(ji,jj  )) + &
!!                 bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) )
!!          zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )   !  where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.
!!          en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1)
!!       END DO
!!    END DO
!!bfr   - end commented area
      !
      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke       (Axell JGR 2002)
         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
         !
         !                        !* total energy produce by LC : cumulative sum over jk
         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * fsdepw(:,:,1) * fse3w(:,:,1)
         DO jk = 2, jpk
            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
         END DO
         !                        !* finite Langmuir Circulation depth
         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
         DO jk = jpkm1, 2, -1
            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us 
               DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
                  zus  = zcof * taum(ji,jj)
                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
               END DO
            END DO
         END DO
         !                               ! finite LC depth
         DO jj = 1, jpj 
            DO ji = 1, jpi
               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
            END DO
         END DO
         zcof = 0.016 / SQRT( zrhoa * zcdrag )
!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.
                  zus  = zcof * SQRT( taum(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 * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) /   &
                     &   zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1)
               END DO
            END DO
         END DO
         !
      ENDIF
      !
      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      !                     !  Now Turbulent kinetic energy (output in en)
      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
      !
      DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
         DO jj = 1, jpj                 ! here avmu, avmv used as workspace
            DO ji = 1, jpi
               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   & 
                  &                            / (  fse3uw_n(ji,jj,jk)               &
                  &                              *  fse3uw_b(ji,jj,jk)  )
               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   &
                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   &
                  &                            / (  fse3vw_n(ji,jj,jk)               &
                  &                              *  fse3vw_b(ji,jj,jk)  )
            END DO
         END DO
      END DO
      !
      DO jk = 2, jpkm1           !* Matrix and right hand side in en
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1   ! vector opt.
               zcof   = zfact1 * tmask(ji,jj,jk)
# if defined key_zdftmx_new
               ! key_zdftmx_new: New internal wave-driven param: set a minimum value for Kz on TKE (ensure numerical stability)
               zzd_up = zcof * ( MAX( avm(ji,jj,jk+1) + avm(ji,jj,jk), 2.e-5_wp ) )   &  ! upper diagonal
                  &          / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk  ) )
               zzd_lw = zcof * ( MAX( avm(ji,jj,jk) + avm(ji,jj,jk-1), 2.e-5_wp ) )   &  ! lower diagonal
                  &          / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
# else
               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  ) )
# endif
                  !                                                           ! shear prod. at w-point weightened by mask
               zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
                  !
               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
               zd_lw(ji,jj,jk) = zzd_lw
               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
               !
               !                                   ! right hand side in en
               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    &
                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) &
                  &                                 * wmask(ji,jj,jk)
            END DO
         END DO
      END DO
      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1    ! vector opt.
               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
            END DO
         END DO
      END DO
      !
      ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
      DO jj = 2, jpjm1
         DO ji = fs_2, fs_jpim1   ! vector opt.
            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
         END DO
      END DO
      DO jk = 3, jpkm1
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1    ! vector opt.
               zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1)
            END DO
         END DO
      END DO
      !
      ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
      DO jj = 2, jpjm1
         DO ji = fs_2, fs_jpim1   ! vector opt.
            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
         END DO
      END DO
      DO jk = jpk-2, 2, -1
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1    ! vector opt.
               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
            END DO
         END DO
      END DO
      DO jk = 2, jpkm1                             ! set the minimum value of tke
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1   ! vector opt.
               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
            END DO
         END DO
      END DO

      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      !                            !  TKE due to surface and internal wave breaking
      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
         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) )   &
                     &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
               END DO
            END DO
         END DO
      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
         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) )   &
                  &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
            END DO
         END DO
      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
!CDIR NOVERRCHK
         DO jk = 2, jpkm1
!CDIR NOVERRCHK
            DO jj = 2, jpjm1
!CDIR NOVERRCHK
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
                  zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress 
                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean 
                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
                     &                        * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
               END DO
            END DO
         END DO
      ENDIF
      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
      !
      CALL wrk_dealloc( jpi,jpj, imlc )    ! integer
      CALL wrk_dealloc( jpi,jpj, zhlc ) 
      CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
      !
      IF( nn_timing == 1 )  CALL timing_stop('tke_tke')
      !
   END SUBROUTINE tke_tke


   SUBROUTINE tke_avn
      !!----------------------------------------------------------------------
      !!                   ***  ROUTINE tke_avn  ***
      !!
      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
      !!
      !! ** Method  :   At this stage, en, the now TKE, is known (computed in 
      !!              the tke_tke routine). First, the now mixing lenth is 
      !!      computed from en and the strafification (N^2), then the mixings
      !!      coefficients are computed.
      !!              - Mixing length : a first evaluation of the mixing lengh
      !!      scales is:
      !!                      mxl = sqrt(2*en) / N  
      !!      where N is the brunt-vaisala frequency, with a minimum value set
      !!      to rmxl_min (rn_mxl0) 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
      !!----------------------------------------------------------------------
      INTEGER  ::   ji, jj, jk   ! dummy loop indices
      REAL(wp) ::   zrn2, zraug, zcoef, zav     ! local scalars
      REAL(wp) ::   zdku, zpdlr, zri, zsqen     !   -      -
      REAL(wp) ::   zdkv, zemxl, zemlm, zemlp   !   -      -
      REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld
      !!--------------------------------------------------------------------
      !
      IF( nn_timing == 1 )  CALL timing_start('tke_avn')

      CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 

      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      !                     !  Mixing length
      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      !
      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
      !
      ! initialisation of interior minimum value (avoid a 2d loop with mikt)
      zmxlm(:,:,:)  = rmxl_min    
      zmxld(:,:,:)  = rmxl_min
      !
      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1
               zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
               zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) )
            END DO
         END DO
      ELSE 
         zmxlm(:,:,1) = rn_mxl0
      ENDIF
      !
!CDIR NOVERRCHK
      DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2)
!CDIR NOVERRCHK
         DO jj = 2, jpjm1
!CDIR NOVERRCHK
            DO ji = fs_2, fs_jpim1   ! vector opt.
               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) )
            END DO
         END DO
      END DO
      !
      !                     !* Physical limits for the mixing length
      !
      zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value 
      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
      !
      SELECT CASE ( nn_mxl )
      !
      ! where wmask = 0 set zmxlm == fse3w
      CASE ( 0 )           ! bounded by the distance to surface and bottom
         DO jk = 2, jpkm1
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   &
                  &            fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) )
                  ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
               END DO
            END DO
         END DO
         !
      CASE ( 1 )           ! bounded by the vertical scale factor
         DO jk = 2, jpkm1
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
                  zmxlm(ji,jj,jk) = zemxl
                  zmxld(ji,jj,jk) = zemxl
               END DO
            END DO
         END DO
         !
      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
         DO jk = 2, jpkm1         ! from the surface to the bottom :
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
               END DO
            END DO
         END DO
         DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
                  zmxlm(ji,jj,jk) = zemxl
                  zmxld(ji,jj,jk) = zemxl
               END DO
            END DO
         END DO
         !
      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
         DO jk = 2, jpkm1         ! from the surface to the bottom : lup
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
               END DO
            END DO
         END DO
         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
               END DO
            END DO
         END DO
!CDIR NOVERRCHK
         DO jk = 2, jpkm1
!CDIR NOVERRCHK
            DO jj = 2, jpjm1
!CDIR NOVERRCHK
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
                  zmxlm(ji,jj,jk) = zemlm
                  zmxld(ji,jj,jk) = zemlp
               END DO
            END DO
         END DO
         !
      END SELECT
      !
# if defined key_c1d
      e_dis(:,:,:) = zmxld(:,:,:)      ! c1d configuration : save mixing and dissipation turbulent length scales
      e_mix(:,:,:) = zmxlm(:,:,:)
# endif

      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt)
      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!CDIR NOVERRCHK
      DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
!CDIR NOVERRCHK
         DO jj = 2, jpjm1
!CDIR NOVERRCHK
            DO ji = fs_2, fs_jpim1   ! vector opt.
               zsqen = SQRT( en(ji,jj,jk) )
               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk)
               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
            END DO
         END DO
      END DO
      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
      !
      DO jk = 2, jpkm1            !* vertical eddy viscosity at wu- and wv-points
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1   ! vector opt.
               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk)
               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk)
            END DO
         END DO
      END DO
      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
      !
      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
         DO jk = 2, jpkm1
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk)
                  !                                          ! shear
                  zdku = avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) * ( ub(ji-1,jj,jk-1) - ub(ji-1,jj,jk) )   &
                    &  + avmu(ji  ,jj,jk) * ( un(ji  ,jj,jk-1) - un(ji  ,jj,jk) ) * ( ub(ji  ,jj,jk-1) - ub(ji  ,jj,jk) )
                  zdkv = avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) * ( vb(ji,jj-1,jk-1) - vb(ji,jj-1,jk) )   &
                    &  + avmv(ji,jj  ,jk) * ( vn(ji,jj  ,jk-1) - vn(ji,jj  ,jk) ) * ( vb(ji,jj  ,jk-1) - vb(ji,jj  ,jk) )
                  !                                          ! local Richardson number
                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear )
                  zpdlr = MAX(  0.1_wp,  0.2 / MAX( 0.2 , zri )  )
!!gm and even better with the use of the "true" ri_crit=0.22222...  (this change the results!)
!!gm              zpdlr = MAX(  0.1_wp,  ri_crit / MAX( ri_crit , zri )  )
                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
# if defined key_c1d
                  e_pdl(ji,jj,jk) = zpdlr * wmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number
                  e_ric(ji,jj,jk) = zri   * wmask(ji,jj,jk)  ! c1d config. : save Ri
# endif
              END DO
            END DO
         END DO
      ENDIF
      CALL lbc_lnk( avt, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged)

      IF(ln_ctl) THEN
         CALL prt_ctl( tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
         CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask,                   &
            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
      ENDIF
      !
      CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
      !
      IF( nn_timing == 1 )  CALL timing_stop('tke_avn')
      !
   END SUBROUTINE tke_avn


   SUBROUTINE zdf_tke_init
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE zdf_tke_init  ***
      !!                     
      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
      !!              viscosity when using a tke turbulent closure scheme
      !!
      !! ** Method  :   Read the namzdf_tke namelist and check the parameters
      !!              called at the first timestep (nit000)
      !!
      !! ** input   :   Namlist namzdf_tke
      !!
      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
      !!----------------------------------------------------------------------
      INTEGER ::   ji, jj, jk   ! dummy loop indices
      INTEGER ::   ios, ierr
      !!
      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,   &
         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   &
         &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   &
         &                 nn_etau , nn_htau  , rn_efr   
      !!----------------------------------------------------------------------
      !
      REWIND( numnam_ref )              ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy
      READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist', lwp )

      REWIND( numnam_cfg )              ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
      READ  ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 )
902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist', lwp )
      IF(lwm) WRITE ( numond, namzdf_tke )
      !
      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
# if defined key_zdftmx_new
      ! key_zdftmx_new: New internal wave-driven param: specified value of rn_emin & rmxl_min are used
      rn_emin  = 1.e-10_wp
      rmxl_min = 1.e-03_wp
      IF(lwp) THEN                  ! Control print
         WRITE(numout,*)
         WRITE(numout,*) 'zdf_tke_init :  New tidal mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 '
         WRITE(numout,*) '~~~~~~~~~~~~'
      ENDIF
# else
      rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
# endif
      !
      IF(lwp) THEN                    !* Control print
         WRITE(numout,*)
         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
         WRITE(numout,*) '~~~~~~~~~~~~'
         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters'
         WRITE(numout,*) '      coef. to compute avt                        rn_ediff  = ', rn_ediff
         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss  = ', rn_ediss
         WRITE(numout,*) '      tke surface input coef.                     rn_ebb    = ', rn_ebb
         WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin
         WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0
         WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear
         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
         WRITE(numout,*) '      surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
         WRITE(numout,*) '      surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
         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,*) '      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,*) '      fraction of en which pene. the thermocline  rn_efr    = ', rn_efr
         WRITE(numout,*)
         WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri
      ENDIF
      !
      !                              ! allocate tke arrays
      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
      !
      !                               !* Check of some namelist values
      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )

      IF( ln_mxl0 ) THEN
         IF(lwp) WRITE(numout,*) '   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
         rn_mxl0 = rmxl_min
      ENDIF
      
      IF( nn_etau == 2  ) THEN
          ierr = zdf_mxl_alloc()
          nmln(:,:) = nlb10           ! Initialization of nmln
      ENDIF

      !                               !* depth of penetration of surface tke
      IF( nn_etau /= 0 ) THEN      
         SELECT CASE( nn_htau )             ! Choice of the depth of penetration
         CASE( 0 )                                 ! constant depth penetration (here 10 meters)
            htau(:,:) = 10._wp
         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )            
         END SELECT
      ENDIF
      !                               !* set vertical eddy coef. to the background value
      DO jk = 1, jpk
         avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
         avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
         avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
         avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
      END DO
      dissl(:,:,:) = 1.e-12_wp
      !                              
      CALL tke_rst( nit000, 'READ' )  !* read or initialize all required files
      !
   END SUBROUTINE zdf_tke_init


   SUBROUTINE tke_rst( kt, cdrw )
     !!---------------------------------------------------------------------
     !!                   ***  ROUTINE tke_rst  ***
     !!                     
     !! ** Purpose :   Read or write TKE file (en) in restart file
     !!
     !! ** Method  :   use of IOM library
     !!                if the restart does not contain TKE, en is either 
     !!                set to rn_emin or recomputed 
     !!----------------------------------------------------------------------
     INTEGER         , INTENT(in) ::   kt     ! ocean time-step
     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
     !
     INTEGER ::   jit, jk   ! dummy loop indices
     INTEGER ::   id1, id2, id3, id4, id5, id6   ! local integers
     !!----------------------------------------------------------------------
     !
     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
        !                                   ! ---------------
        IF( ln_rstart ) THEN                   !* Read the restart file
           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
           id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
           !
           IF( id1 > 0 ) THEN                       ! 'en' exists
              CALL iom_get( numror, jpdom_autoglo, 'en', en )
              IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
                 CALL iom_get( numror, jpdom_autoglo, 'avt'  , avt   )
                 CALL iom_get( numror, jpdom_autoglo, 'avm'  , avm   )
                 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu  )
                 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv  )
                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
              ELSE                                                 ! one at least array is missing
                 CALL tke_avn                                          ! compute avt, avm, avmu, avmv and dissl (approximation)
              ENDIF
           ELSE                                     ! No TKE array found: initialisation
              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
              en (:,:,:) = rn_emin * tmask(:,:,:)
              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
              !
              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
           ENDIF
        ELSE                                   !* Start from rest
           en(:,:,:) = rn_emin * tmask(:,:,:)
        ENDIF
           !                                   ! Initialize av*_k
           avt_k (:,:,:) = avt (:,:,:)
           avm_k (:,:,:) = avm (:,:,:)
           avmu_k(:,:,:) = avmu(:,:,:)
           avmv_k(:,:,:) = avmv(:,:,:)
           !
     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
        !                                   ! -------------------
        IF(lwp) WRITE(numout,*) '---- tke-rst ----'
        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  )
        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  )
        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )
        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl  )
        !
     ENDIF
     !
   END SUBROUTINE tke_rst

#else
   !!----------------------------------------------------------------------
   !!   Dummy module :                                        NO TKE scheme
   !!----------------------------------------------------------------------
   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
CONTAINS
   SUBROUTINE zdf_tke_init           ! Dummy routine
   END SUBROUTINE zdf_tke_init
   SUBROUTINE zdf_tke( kt )          ! Dummy routine
      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
   END SUBROUTINE zdf_tke
   SUBROUTINE tke_rst( kt, cdrw )
     CHARACTER(len=*) ::   cdrw
     WRITE(*,*) 'tke_rst: You should not have seen this print! error?', kt, cdwr
   END SUBROUTINE tke_rst
#endif

   !!======================================================================
END MODULE zdftke