MODULE zdfkpp
   !!======================================================================
   !!                       ***  MODULE  zdfkpp  ***
   !! Ocean physics:  vertical mixing coefficient compute from the KPP 
   !!                 turbulent closure parameterization
   !!=====================================================================
   !! History :  OPA  ! 2000-03 (W.G. Large, J. Chanut) Original code
   !!            8.1  ! 2002-06 (J.M. Molines) for real case CLIPPER  
   !!            8.2  ! 2003-10 (Chanut J.) re-writting
   !!   NEMO     1.0  ! 2005-01 (C. Ethe, G. Madec) Free form, F90 + creation of tra_kpp routine
   !!            3.3  ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
   !!----------------------------------------------------------------------
#if defined key_zdfkpp   ||   defined key_esopa
   !!----------------------------------------------------------------------
   !!   'key_zdfkpp'                                             KPP scheme
   !!----------------------------------------------------------------------
   !!   zdf_kpp       : update momentum and tracer Kz from a kpp scheme
   !!   zdf_kpp_init  : initialization, namelist read, and parameters control
   !!   tra_kpp       : compute and add to the T & S trend the non-local flux
   !!   trc_kpp       : compute and add to the passive tracer trend the non-local flux (lk_top=T)
   !!----------------------------------------------------------------------
   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 eosbn2         ! equation of state
   USE zdfddm         ! double diffusion mixing (avs array)
   USE lib_mpp        ! MPP library
   USE trd_oce        ! ocean trends definition
   USE trdtra         ! tracers trends
   !
   USE in_out_manager ! I/O manager
   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
   USE prtctl         ! Print control
   USE wrk_nemo       ! work arrays
   USE timing         ! Timing
   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)

   IMPLICIT NONE
   PRIVATE

   PUBLIC   zdf_kpp       ! routine called by step.F90
   PUBLIC   zdf_kpp_init  ! routine called by opa.F90
   PUBLIC   tra_kpp       ! routine called by step.F90
   PUBLIC   trc_kpp       ! routine called by trcstp.F90

   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfkpp = .TRUE.    !: KPP vertical mixing flag

   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghats    !: non-local scalar mixing term (gamma/<ws>o)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wt0      !: surface temperature flux for non local flux
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ws0      !: surface salinity flux for non local flux
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hkpp     !: boundary layer depth

   !                          !!* Namelist namzdf_kpp *
   REAL(wp) ::   rn_difmiw     ! constant internal wave viscosity (m2/s)
   REAL(wp) ::   rn_difsiw     ! constant internal wave diffusivity (m2/s)
   REAL(wp) ::   rn_riinfty    ! local Richardson Number limit for shear instability
   REAL(wp) ::   rn_difri      ! maximum shear mixing at Rig = 0    (m2/s)
   REAL(wp) ::   rn_bvsqcon    ! Brunt-Vaisala squared (1/s**2) for maximum convection
   REAL(wp) ::   rn_difcon     ! maximum mixing in interior convection (m2/s)
   INTEGER  ::   nn_ave        ! = 0/1 flag for horizontal average on avt, avmu, avmv

#if defined key_zdfddm
   !                                        !!! ** Double diffusion Mixing
   REAL(wp) ::   difssf  = 1.e-03_wp         ! maximum salt fingering mixing 
   REAL(wp) ::   Rrho0   = 1.9_wp            ! limit for salt  fingering mixing 
   REAL(wp) ::   difsdc  = 1.5e-06_wp        ! maximum diffusive convection mixing
#endif
   LOGICAL  ::   ln_kpprimix       ! Shear instability mixing 

   !                                        !!! ** General constants  **
   REAL(wp) ::   epsln   = 1.0e-20_wp        ! a small positive number    
   REAL(wp) ::   pthird  = 1._wp/3._wp       ! 1/3
   REAL(wp) ::   pfourth = 1._wp/4._wp       ! 1/4

   !                                        !!! ** Boundary Layer Turbulence Parameters  **
   REAL(wp) ::   vonk     = 0.4_wp           ! von Karman's constant
   REAL(wp) ::   epsilon  = 0.1_wp           ! nondimensional extent of the surface layer
   REAL(wp) ::   rconc1   = 5.0_wp           ! standard flux profile function parmaeters
   REAL(wp) ::   rconc2   = 16.0_wp          !         "        "
   REAL(wp) ::   rconcm   = 8.38_wp          ! momentum flux profile fit
   REAL(wp) ::   rconam   = 1.26_wp          !         "       "
   REAL(wp) ::   rzetam   = -.20_wp          !         "       "       
   REAL(wp) ::   rconcs   = 98.96_wp         !  scalar  flux profile fit
   REAL(wp) ::   rconas   = -28.86_wp        !         "       "
   REAL(wp) ::   rzetas   = -1.0_wp          !         "       "  
   
   !                                        !!! ** Boundary Layer Depth Diagnostic  **
   REAL(wp) ::   Ricr     = 0.3_wp           ! critical bulk Richardson Number
   REAL(wp) ::   rcekman  = 0.7_wp           ! coefficient for ekman depth  
   REAL(wp) ::   rcmonob  = 1.0_wp           ! coefficient for Monin-Obukhov depth 
   REAL(wp) ::   rconcv   = 1.7_wp           ! ratio of interior buoyancy frequency to its value at entrainment depth
   REAL(wp) ::   hbf      = 1.0_wp           ! fraction of bound. layer depth to which absorbed solar 
      !                                      ! rad. and contributes to surf. buo. forcing
   REAL(wp) ::   Vtc                         ! function of rconcv,rconcs,epsilon,vonk,Ricr
   
   !                                        !!! ** Nonlocal Boundary Layer Mixing **
   REAL(wp) ::   rcstar   = 5.0_wp           ! coefficient for convective nonlocal transport
   REAL(wp) ::   rcs      = 1.0e-3_wp        ! conversion: mm/s ==> m/s   
   REAL(wp) ::   rcg                         ! non-dimensional coefficient for nonlocal transport

#if ! defined key_kppcustom
   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   del     ! array for reference mean values of vertical integration 
#endif

#if defined key_kpplktb
   !                                         !!! ** Parameters for lookup table for turbulent velocity scales ** 
   INTEGER, PARAMETER ::   nilktb   = 892     ! number of values for zehat in KPP lookup table
   INTEGER, PARAMETER ::   njlktb   = 482     ! number of values for ustar in KPP lookup table
   INTEGER, PARAMETER ::   nilktbm1 = nilktb-1   !
   INTEGER, PARAMETER ::   njlktbm1 = njlktb-1   !

   REAL(wp), DIMENSION(nilktb,njlktb) ::   wmlktb   ! lookup table for the turbulent vertical velocity scale (momentum)
   REAL(wp), DIMENSION(nilktb,njlktb) ::   wslktb   ! lookup table for the turbulent vertical velocity scale (tracers)

   REAL(wp) ::   dehatmin = -4.e-7_wp    ! minimum limit for zhat in lookup table (m3/s3) 
   REAL(wp) ::   dehatmax = 0._wp        ! maximum limit for zhat in lookup table (m3/s3)
   REAL(wp) ::   ustmin   = 0._wp        ! minimum limit for ustar in lookup table (m/s)
   REAL(wp) ::   ustmax   = 0.04_wp      ! maximum limit for ustar in lookup table (m/s)    
   REAL(wp) ::   dezehat                 ! delta zhat in lookup table
   REAL(wp) ::   deustar                 ! delta ustar in lookup table
#endif
   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   ratt   ! attenuation coef  (already defines in module traqsr, 
   !                                    ! but only if the solar radiation penetration is considered)
   
   !                                    !!! * penetrative solar radiation coefficient *
   REAL(wp) ::   rabs = 0.58_wp          ! fraction associated with xsi1
   REAL(wp) ::   xsi1 = 0.35_wp          ! first depth of extinction 
   REAL(wp) ::   xsi2 = 23.0_wp          ! second depth of extinction 
      !                           ! (default values: water type Ib) 

   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   etmean, eumean, evmean   ! coeff. used for hor. smoothing at t-, u- & v-points
        
#if defined key_c1d
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rig    !: gradient Richardson number
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rib    !: bulk Richardson number
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   buof   !: buoyancy forcing
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mols   !: moning-Obukhov length scale 
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ekdp   !: Ekman depth
#endif

   INTEGER  ::   jip = 62 , jjp = 111

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

   INTEGER FUNCTION zdf_kpp_alloc()
      !!----------------------------------------------------------------------
      !!                 ***  FUNCTION zdf_kpp_alloc  ***
      !!----------------------------------------------------------------------
      ALLOCATE( ghats(jpi,jpj,jpk), wt0(jpi,jpj), ws0(jpi,jpj), hkpp(jpi,jpj), &
#if ! defined key_kpplktb
         &      del(jpk,jpk),                                                  &
#endif
         &      ratt(jpk),                                                     &
         &      etmean(jpi,jpj,jpk), eumean(jpi,jpj,jpk), evmean(jpi,jpj,jpk), &
#if defined key_c1d
         &      rig (jpi,jpj,jpk), rib(jpi,jpj,jpk), buof(jpi,jpj,jpk),        &
         &      mols(jpi,jpj,jpk), ekdp(jpi,jpj),                              &
#endif
         &      STAT= zdf_kpp_alloc )
         !
      IF( lk_mpp             )   CALL mpp_sum ( zdf_kpp_alloc )
      IF( zdf_kpp_alloc /= 0 )   CALL ctl_warn('zdf_kpp_alloc: failed to allocate arrays')
   END FUNCTION zdf_kpp_alloc


   SUBROUTINE zdf_kpp( kt )
      !!----------------------------------------------------------------------
      !!                   ***  ROUTINE zdf_kpp  ***
      !!
      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
      !!      coefficients and non local mixing using K-profile parameterization
      !!
      !! ** Method :   The boundary layer depth hkpp is diagnosed at tracer points
      !!      from profiles of buoyancy, and shear, and the surface forcing.
      !!      Above hbl (sigma=-z/hbl <1) the mixing coefficients are computed from 
      !!
      !!                      Kx =  hkpp  Wx(sigma) G(sigma)  
      !!
      !!             and the non local term ghat = Cs / Ws(sigma) / hkpp
      !!      Below hkpp  the coefficients are the sum of mixing due to internal waves
      !!      shear instability and double diffusion.
      !!
      !!      -1- Compute the now interior vertical mixing coefficients at all depths. 
      !!      -2- Diagnose the boundary layer depth. 
      !!      -3- Compute the now boundary layer vertical mixing coefficients. 
      !!      -4- Compute the now vertical eddy vicosity and diffusivity.
      !!      -5- Smoothing
      !!
      !!        N.B. The computation is done from jk=2 to jpkm1 
      !!             Surface value of avt avmu avmv are set once a time to zero
      !!             in routine zdf_kpp_init.
      !!
      !! ** Action  :   update the non-local terms ghats
      !!                update avt, avmu, avmv (before vertical eddy coef.)
      !!
      !! References : Large W.G., Mc Williams J.C. and Doney S.C.              
      !!         Reviews of Geophysics, 32, 4, November 1994
      !!         Comments in the code refer to this paper, particularly 
      !!         the equation number. (LMD94, here after)
      !!----------------------------------------------------------------------
      USE oce     , zviscos => ua   ! temp. array for viscosities use ua as workspace
      USE oce     , zdiffut => va   ! temp. array for diffusivities use sa as workspace
      !!
      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
      !!
      INTEGER ::   ji, jj, jk              ! dummy loop indices
      INTEGER ::   ikbot, jkmax, jkm1, jkp2   !

      REAL(wp) ::   ztx, zty, zflageos, zstabl, zbuofdep,zucube     !
      REAL(wp) ::   zrhos, zalbet, zbeta, zthermal, zhalin, zatt1   !
      REAL(wp) ::   zref, zt, zs, zh, zu, zv, zrh                   ! Bulk richardson number
      REAL(wp) ::   zrib, zrinum, zdVsq, zVtsq                      !
      REAL(wp) ::   zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm   ! Velocity scales
#if defined key_kpplktb
      INTEGER ::    il, jl                                          ! Lookup table or Analytical functions 
      REAL(wp) ::   ud, zfrac, ufrac, zwam, zwbm, zwas, zwbs        !
#else
      REAL(wp) ::   zwsun, zwmun, zcons, zconm, zwcons, zwconm      !
#endif
      REAL(wp) ::   zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed   ! In situ density
#if ! defined key_kppcustom     
      INTEGER  ::   jm                          ! dummy loop indices
      REAL(wp) ::   zr1, zr2, zr3, zr4, zrhop   ! Compression terms
#endif
      REAL(wp) ::   zflag, ztemp, zrn2, zdep21, zdep32, zdep43
      REAL(wp) ::   zdku2, zdkv2, ze3sqr, zsh2, zri, zfri          ! Interior richardson mixing
      REAL(wp), POINTER, DIMENSION(:,:) ::   zmoek                 ! Moning-Obukov limitation
      REAL(wp), POINTER, DIMENSION(:)   ::   zmoa, zekman                
      REAL(wp)                          ::   zmob, zek
      REAL(wp), POINTER, DIMENSION(:,:) ::   zdepw, zdift, zvisc   ! The pipe 
      REAL(wp), POINTER, DIMENSION(:,:) ::   zdept
      REAL(wp), POINTER, DIMENSION(:,:) ::   zriblk
      REAL(wp), POINTER, DIMENSION(:)   ::   zhmax, zria, zhbl 
      REAL(wp) ::   zflagri, zflagek, zflagmo, zflagh, zflagkb   !
      REAL(wp), POINTER, DIMENSION(:)   ::   za2m, za3m, zkmpm, za2t, za3t, zkmpt   ! Shape function (G)
      REAL(wp) ::   zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t
#if defined key_zdfddm
      REAL(wp) ::   zrw, zkm1s                    ! local scalars
      REAL(wp) ::   zrrau, zdt, zds, zavdds, zavddt, zinr   ! double diffusion mixing
      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdifs
      REAL(wp), POINTER, DIMENSION(:)     ::   za2s, za3s, zkmps
      REAL(wp), POINTER, DIMENSION(:,:)   ::   zblcs
      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdiffus
#endif
      REAL(wp), POINTER, DIMENSION(:,:) ::   zBo, zBosol, zustar         ! Surface buoyancy forcing, friction velocity
      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask, zblcm, zblct
      !!--------------------------------------------------------------------
      !
      IF( nn_timing == 1 )  CALL timing_start('zdf_kpp')
      !
      CALL wrk_alloc( jpi, zmoa, zekman, zhmax, zria, zhbl )
      CALL wrk_alloc( jpi, za2m, za3m, zkmpm, za2t, za3t, zkmpt )
      CALL wrk_alloc( jpi,2, zriblk )
      CALL wrk_alloc( jpi,3, zmoek, kjstart = 0 )
      CALL wrk_alloc( jpi,3, zdept )
      CALL wrk_alloc( jpi,4, zdepw, zdift, zvisc )
      CALL wrk_alloc( jpi,jpj, zBo, zBosol, zustar )
      CALL wrk_alloc( jpi,jpk, zmask, zblcm, zblct )
#if defined key_zdfddm
      CALL wrk_alloc( jpi,4, zdifs )
      CALL wrk_alloc( jpi, zmoa, za2s, za3s, zkmps )
      CALL wrk_alloc( jpi,jpk, zblcs )
      CALL wrk_alloc( jpi,jpi,jpk, zdiffus )
#endif

      zviscos(:,:,:) = 0._wp
      zblcm  (:,:  ) = 0._wp
      zdiffut(:,:,:) = 0._wp
      zblct  (:,:  ) = 0._wp 
#if defined key_zdfddm
      zdiffus(:,:,:) = 0._wp
      zblcs  (:,:  ) = 0._wp 
#endif
      ghats  (:,:,:) = 0._wp
      zBo    (:,:  ) = 0._wp
      zBosol (:,:  ) = 0._wp
      zustar (:,:  ) = 0._wp
      !
      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      ! I. Interior diffusivity and viscosity at w points ( T interfaces)
      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  
      DO jk = 2, jpkm1
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1 
               ! Mixing due to internal waves breaking
               ! -------------------------------------
               avmu(ji,jj,jk)  = rn_difmiw 
               avt (ji,jj,jk)  = rn_difsiw             
               ! Mixing due to vertical shear instability
               ! -------------------------------------               
               IF( ln_kpprimix ) THEN          
                  ! Compute the gradient Richardson  number at interfaces (zri):
                  ! LMD94, eq. 27 (is vertical smoothing needed : Rig=N^2 / (dz(u))^2 + (dz(v))^2
                  zdku2 =   ( un(ji - 1,jj,jk - 1) - un(ji - 1,jj,jk) ) &
                     &    * ( un(ji - 1,jj,jk - 1) - un(ji - 1,jj,jk) ) &
                     &    + ( un(ji    ,jj,jk - 1) - un(ji    ,jj,jk) ) &
                     &    * ( un(ji    ,jj,jk - 1) - un(ji    ,jj,jk) )  
                  
                  zdkv2 =   ( vn(ji,jj - 1,jk - 1) - vn(ji,jj - 1,jk) ) &
                     &    * ( vn(ji,jj - 1,jk - 1) - vn(ji,jj - 1,jk) ) &
                     &    + ( vn(ji,    jj,jk - 1) - vn(ji,    jj,jk) ) &
                     &    * ( vn(ji,    jj,jk - 1) - vn(ji,    jj,jk) )  

                  ze3sqr = 1. / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) )
                  ! Square of vertical shear  at interfaces
                  zsh2   = 0.5 * ( zdku2 + zdkv2 ) * ze3sqr
                  zri    = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + epsln ) 
#if defined key_c1d
                  ! save the gradient richardson number
                  rig(ji,jj,jk) = zri * tmask(ji,jj,jk)
#endif                  
                  ! Evaluate f of Ri (zri) for shear instability store in zfri
                  ! LMD94, eq. 28a,b,c, figure 3 ; Rem: p1 is 3, hard coded
                  zfri  = MAX( zri , 0. )
                  zfri  = MIN( zfri / rn_riinfty , 1.0 )
                  zfri  = ( 1.0 - zfri * zfri )
                  zfri  = zfri * zfri  * zfri
                  ! add shear contribution to mixing coef. 
                  avmu(ji,jj,jk) =  avmu(ji,jj,jk) + rn_difri * zfri   
                  avt (ji,jj,jk) =  avt (ji,jj,jk) + rn_difri * zfri    
               ENDIF
               !
#if defined key_zdfddm 
               !
               !  Double diffusion mixing ; NOT IN ROUTINE ZDFDDM.F90
               ! -------------------------
               avs (ji,jj,jk) =  avt (ji,jj,jk)   

               ! R=zrau = (alpha / beta) (dk[t] / dk[s])
               zrw =   ( fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk) )   &
                  &  / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) ) 
               !
               zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  ) * tmask(ji,jj,jk)
               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  ) * tmask(ji,jj,jk)
               !
               zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )
               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) 
               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
               zrrau = MAX(  epsln , zdt / zds  )    ! only retains positive value of zrau
               !
               IF( zrrau > 1. .AND. zds > 0.) THEN                        ! Salt fingering case.
                  !                                                       !---------------------
                  ! Compute interior diffusivity for double diffusive mixing of salinity. 
                  ! Upper bound "zrrau" by "Rrho0"; (Rrho0=1.9, difcoefnuf=0.001).
                  ! After that set interior diffusivity for double diffusive mixing of temperature
                  zavdds = MIN( zrrau, Rrho0 )
                  zavdds = ( zavdds - 1.0 ) / ( Rrho0 - 1.0 )
                  zavdds = 1.0 - zavdds * zavdds 
                  zavdds = zavdds * zavdds * zavdds 
                  zavdds = difssf * zavdds 
                  zavddt = 0.7 * zavdds
                  !
               ELSEIF( zrrau < 1. .AND. zrrau > 0. .AND. zds < 0.) THEN   ! Diffusive convection case.
                  !                                                       !---------------------------
                  ! Compute interior diffusivity for double diffusive mixing of temperature (Marmorino and Caldwell, 1976); 
                  ! Compute interior diffusivity for double diffusive mixing of salinity 
                  zinr   = 1. / zrrau
                  zavddt = 0.909 * EXP( 4.6 * EXP( -0.54* ( zinr - 1. ) ) ) 
                  zavddt = difsdc * zavddt
                  IF( zrrau < 0.5) THEN   ;   zavdds = zavddt * 0.15 * zrrau
                  ELSE                    ;   zavdds = zavddt * (1.85 * zrrau - 0.85 ) 
                  ENDIF
               ELSE
                  zavddt = 0.
                  zavdds = 0.
               ENDIF 
               ! Add double diffusion contribution to temperature and salinity  mixing coefficients.
               avt (ji,jj,jk) =  avt (ji,jj,jk) +  zavddt 
               avs (ji,jj,jk) =  avs (ji,jj,jk) +  zavdds         
#endif                      
            END DO
         END DO
      END DO


      ! Radiative (zBosol) and non radiative (zBo) surface buoyancy
      !JMM at the time zdfkpp is called, q still holds the sum q + qsr
      !---------------------------------------------------------------------
      DO jj = 2, jpjm1
         DO ji = fs_2, fs_jpim1           
            zrhos    = rau0 * ( 1._wp + rhd(ji,jj,1) ) * tmask(ji,jj,1)
            zthermal = rab_n(ji,jj,1,jp_tem) / ( rcp * zrhos + epsln )
            zbeta    = rab_n(ji,jj,1,jp_sal)
            zhalin   = zbeta * tsn(ji,jj,1,jp_sal) * rcs
            !
            ! Radiative surface buoyancy force
            zBosol(ji,jj) = grav * zthermal * qsr(ji,jj)
            ! Non radiative surface buoyancy force
            zBo   (ji,jj) = grav * zthermal * qns(ji,jj) -  grav * zhalin * ( emp(ji,jj)-rnf(ji,jj) )  &
               &                                         -  grav * zbeta * rcs * sfx(ji,jj)
            ! Surface Temperature flux for non-local term
            wt0(ji,jj) = - ( qsr(ji,jj) + qns(ji,jj) )* r1_rau0_rcp * tmask(ji,jj,1)
            ! Surface salinity flux for non-local term
            ws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal)                          &
               &             + sfx(ji,jj)                                     ) * rcs * tmask(ji,jj,1) 
         END DO
      END DO

      zflageos = 0.5 + SIGN( 0.5, nn_eos - 1. ) 
      !  Compute surface buoyancy forcing, Monin Obukhov and Ekman depths  
      !------------------------------------------------------------------    
      DO jj = 2, jpjm1
         DO ji = fs_2, fs_jpim1
            !  Reference surface density = density at first T point level   
            zrhos         = rhop(ji,jj,1) + zflageos * rau0 * ( 1. - tmask(ji,jj,1) )  
            ! Friction velocity (zustar), at T-point : LMD94 eq. 2
            zustar(ji,jj) = SQRT( taum(ji,jj) / ( zrhos +  epsln ) )
         END DO
      END DO

!CDIR NOVERRCHK  
      !                                               ! ===============
      DO jj = 2, jpjm1                                 !  Vertical slab
         !                                             ! ===============
         
         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
         ! II Compute Boundary layer mixing coef. and diagnose the new boundary layer depth
         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
         
         ! Initialization
         jkmax       = 0
         zdept (:,:) = 0.
         zdepw (:,:) = 0.
         zriblk(:,:) = 0.
         zmoek (:,:) = 0.
         zvisc (:,:) = 0.
         zdift (:,:) = 0.
#if defined key_zdfddm
         zdifs (:,:) = 0.
#endif
         zmask (:,:) = 0.
         DO ji = fs_2, fs_jpim1
            zria(ji ) = 0.
            ! Maximum boundary layer depth
            ikbot     = mbkt(ji,jj)    ! ikbot is the last T point in the water
            zhmax(ji) = fsdept(ji,jj,ikbot) - 0.001      
            ! Compute Monin obukhov length scale at the surface and Ekman depth:
            zbuofdep   = zBo(ji,jj) + zBosol(ji,jj) * ratt(1)
            zekman(ji) = rcekman * zustar(ji,jj) / ( ABS( ff(ji,jj) ) + epsln )
            zucube     = zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) 
            zmoa(ji)   = zucube / ( vonk * ( zbuofdep + epsln ) )   
#if defined key_c1d
            ! store the surface buoyancy forcing
            zstabl        = 0.5 + SIGN( 0.5, zbuofdep )
            buof(ji,jj,1) = zbuofdep * tmask(ji,jj,1)
            ! store the moning-oboukov length scale at surface
            zmob          = zstabl * zmoa(ji) + ( 1.0 - zstabl ) * fsdept(ji,jj,1)
            mols(ji,jj,1) = MIN( zmob , zhmax(ji) ) * tmask(ji,jj,1)
            ! store Ekman depth
            zek           = zstabl * zekman(ji) + ( 1.0 - zstabl ) * fsdept(ji,jj,1)  
            ekdp(ji,jj )  = MIN( zek , zhmax(ji) ) * tmask(ji,jj,1)  
#endif 
         END DO     
         ! Compute the pipe
         ! ---------------------
         DO jk = 2, jpkm1
            DO ji = fs_2, fs_jpim1
               ! Compute bfsfc = Bo + radiative contribution down to hbf*depht
               zbuofdep = zBo(ji,jj) + zBosol(ji,jj) * ratt(jk)
               ! Flag (zstabl  = 1) if positive forcing
               zstabl   =  0.5 + SIGN(  0.5, zbuofdep)

               !   Compute bulk richardson number zrib at depht 
               !-------------------------------------------------------
               !                           [Br - B(d)] * d         zrinum
               !             Rib(z) = ----------------------- = -------------
               !                       |Vr - V(d)|^2 + Vt(d)^2   zdVsq + zVtsq
               !
               ! First compute zt,zs,zu,zv = means in the surface layer < epsilon*depht  
               ! Else surface values are taken at the first T level.
               ! For stability, resolved vertical shear is computed with "before velocities".
               zref = epsilon * fsdept(ji,jj,jk)
#if defined key_kppcustom
               ! zref = gdept(1)
               zref = fsdept(ji,jj,1)
               zt   = tsn(ji,jj,1,jp_tem)
               zs   = tsn(ji,jj,1,jp_sal)
               zrh  = rhop(ji,jj,1)
               zu   = ( ub(ji,jj,1) + ub(ji - 1,jj    ,1) ) / MAX( 1. , umask(ji,jj,1) + umask(ji - 1,jj   ,1) )
               zv   = ( vb(ji,jj,1) + vb(ji    ,jj - 1,1) ) / MAX( 1. , vmask(ji,jj,1) + vmask(ji   ,jj - 1,1) )
#else
               zt   = 0.
               zs   = 0.
               zu   = 0.
               zv   = 0.
               zrh  = 0.
               ! vertically integration over the upper epsilon*gdept(jk) ; del () array is computed once in zdf_kpp_init
               DO jm = 1, jpkm1
                  zt   = zt  + del(jk,jm) * tsn(ji,jj,jm,jp_tem)
                  zs   = zs  + del(jk,jm) * tsn(ji,jj,jm,jp_sal)
                  zu   = zu  + 0.5 * del(jk,jm) &
                     &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) &
                     &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
                  zv   = zv  + 0.5 * del(jk,jm) &
                     &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) &
                     &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
                  zrh  = zrh + del(jk,jm) * rhop(ji,jj,jm)
               END DO
#endif
               zsr = SQRT( ABS( tsn(ji,jj,jk,jp_sal) ) )
               ! depth
               zh = fsdept(ji,jj,jk)
               ! compute compression terms on density
               ze  = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6
               zbw = (  1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4
               zb  = zbw + ze * zs
               
               zd  = -2.042967e-2
               zc  =   (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896
               zaw = ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt - 4.721788
               za  = ( zd*zsr + zc ) *zs + zaw
               
               zb1 =   (-0.1909078*zt+7.390729 ) *zt-55.87545
               za1 = ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077
               zkw = ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt + 2098.925 ) *zt+190925.6
               zk0 = ( zb1*zsr + za1 )*zs + zkw
               zcomp =   1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )
               
#if defined key_kppcustom
               ! potential density of water(zrh = zt,zs at level jk):
               zrhdr = zrh / zcomp
#else
               ! potential density of water(ztref,zsref at level jk):
               ! compute volumic mass pure water at atm pressure
               IF ( nn_eos < 1 ) THEN
                  zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4)*zt   &
                     &       -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594
                  ! seawater volumic mass atm pressure
                  zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt   &
                     &   -4.0899e-3 ) *zt+0.824493
                  zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3
                  zr4= 4.8314e-4              
                  ! potential volumic mass (reference to the surface)
                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1                 
                  zrhdr = zrhop / zcomp
               ELSE
                  zrhdr = zrh / zcomp
               ENDIF
#endif
               
               ! potential density of ambiant water at level jk :
               zrhd   = ( rhd(ji,jj,jk) * rau0 + rau0 )  
               
               ! And now the Rib number numerator .
               zrinum = grav * ( zrhd - zrhdr ) / rau0
               zrinum = zrinum * ( fsdept(ji,jj,jk) - zref ) * tmask(ji,jj,jk)
           
               ! Resolved shear contribution to Rib at depth T-point (zdVsq)
               ztx    =   ( ub( ji , jj ,jk)   +  ub(ji - 1, jj ,jk) ) &
                  &     / MAX( 1. , umask( ji , jj ,jk) + umask(ji - 1, jj ,jk) )   
               zty    =   ( vb( ji , jj ,jk)   +  vb(ji  ,jj - 1,jk) ) &
                  &     / MAX( 1., vmask( ji , jj ,jk) + vmask(ji  ,jj - 1,jk) ) 
               
               zdVsq  = ( zu - ztx ) * ( zu - ztx ) + ( zv - zty ) * ( zv - zty )
               
               ! Scalar turbulent velocity scale zws for hbl=gdept
               zscale = zstabl + ( 1.0 - zstabl ) * epsilon
               zehat  = vonk * zscale * fsdept(ji,jj,jk) * zbuofdep
               zucube = zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj)              
               zeta   = zehat / ( zucube + epsln )
               
               IF( zehat > 0. ) THEN
                  ! Stable case
                  zws  = vonk * zustar(ji,jj) / ( 1.0 + rconc1 * zeta )
               ELSE
                  ! Unstable case
#if defined key_kpplktb
                  ! use lookup table
                  zd     = zehat - dehatmin
                  il     = INT( zd / dezehat )
                  il     = MIN( il, nilktbm1 )
                  il     = MAX( il, 1 )
                  
                  ud     = zustar(ji,jj) - ustmin
                  jl     = INT( ud / deustar )
                  jl     = MIN( jl, njlktbm1 )
                  jl     = MAX( jl, 1 )
                  
                  zfrac  = zd / dezehat - FLOAT( il )  
                  ufrac  = ud / deustar - FLOAT( jl )
                  zwas   = ( 1. - zfrac ) * wslktb(il,jl+1) + zfrac * wslktb(il+1,jl+1)
                  zwbs   = ( 1. - zfrac ) * wslktb(il,jl  ) + zfrac * wslktb(il+1,jl  )
                  !
                  zws    = ( 1. - ufrac ) * zwbs + ufrac * zwas
#else
                  ! use analytical functions:
                  zcons  = 0.5 + SIGN( 0.5 , ( rzetas - zeta ) )
                  zwcons = vonk * zustar(ji,jj) * ( ( ABS( rconas - rconcs * zeta ) )**pthird ) 
                  zwsun  = vonk * zustar(ji,jj) * SQRT( ABS ( 1.0 - rconc2 * zeta ) )
                  !
                  zws    = zcons * zwcons +  ( 1.0 - zcons) * zwsun
#endif
               ENDIF
               
               ! Turbulent shear contribution to Rib (zVtsq) bv frequency at levels  ( ie T-point jk)
               zrn2   = 0.5 * ( rn2(ji,jj,jk) + rn2(ji,jj,jk+1) )   
               zbvzed = SQRT( ABS( zrn2 ) ) 
               zVtsq  = fsdept(ji,jj,jk) * zws * zbvzed  * Vtc
               
               ! Finally, the bulk Richardson number at depth fsdept(i,j,k) 
               zrib  = zrinum   / ( zdVsq + zVtsq + epsln )
 
               ! Find subscripts around the boundary layer depth, build the pipe
               ! ----------------------------------------------------------------

               ! Flag (zflagri = 1) if zrib < Ricr  
               zflagri = 0.5 + SIGN( 0.5, ( Ricr - zrib ) )
               !  Flag (zflagh  = 1) if still within overall boundary layer
               zflagh  = 0.5 + SIGN( 0.5, ( fsdept(ji,jj,1) - zdept(ji,2) ) )
               
               ! Ekman layer depth
               zek     = zstabl * zekman(ji) + ( 1.0 - zstabl ) * zhmax(ji)
               zflag   = 0.5 + SIGN( 0.5, ( zek - fsdept(ji,jj,jk-1) ) )
               zek     = zflag * zek + ( 1.0 - zflag ) * zhmax(ji)
               zflagek = 0.5 + SIGN( 0.5, ( zek - fsdept(ji,jj,jk) ) )
               ! Flag (zflagmo = 1) if still within stable Monin-Obukhov and in water
               zmob    = zucube / ( vonk * ( zbuofdep + epsln ) )  
               ztemp   = zstabl * zmob + ( 1.0 - zstabl) * zhmax(ji)
               ztemp   = MIN( ztemp , zhmax(ji) ) 
               zflagmo = 0.5 + SIGN( 0.5, ( ztemp - fsdept(ji,jj,jk) ) )             

               ! No limitation by Monin Obukhov or Ekman depths:
!               zflagek = 1.0
!               zflagmo = 0.5 + SIGN( 0.5, ( zhmax(ji) - fsdept(ji,jj,jk) ) )

               ! Load  pipe via zflagkb  for later calculations
               ! Flag (zflagkb = 1) if zflagh = 1 and (zflagri = 0 or zflagek = 0 or zflagmo = 0)
               zflagkb = zflagh * ( 1.0 - ( zflagri * zflagek * zflagmo ) )

               zmask(ji,jk) = zflagh
               jkp2         = MIN( jk+2 , ikbot )
               jkm1         = MAX( jk-1 , 2 )
               jkmax        = MAX( jkmax, jk * INT( REAL( zflagh+epsln ) ) ) 

               zdept(ji,1)  = zdept(ji,1) + zflagkb * fsdept(ji,jj,jk-1) 
               zdept(ji,2)  = zdept(ji,2) + zflagkb * fsdept(ji,jj,jk  ) 
               zdept(ji,3)  = zdept(ji,3) + zflagkb * fsdept(ji,jj,jk+1) 

               zdepw(ji,1)  = zdepw(ji,1) + zflagkb * fsdepw(ji,jj,jk-1) 
               zdepw(ji,2)  = zdepw(ji,2) + zflagkb * fsdepw(ji,jj,jk  ) 
               zdepw(ji,3)  = zdepw(ji,3) + zflagkb * fsdepw(ji,jj,jk+1)
               zdepw(ji,4)  = zdepw(ji,4) + zflagkb * fsdepw(ji,jj,jkp2)  

               zriblk(ji,1) = zriblk(ji,1) + zflagkb * zria(ji)
               zriblk(ji,2) = zriblk(ji,2) + zflagkb * zrib

               zmoek (ji,0) = zmoek (ji,0) + zflagkb * zek
               zmoek (ji,1) = zmoek (ji,1) + zflagkb * zmoa(ji)
               zmoek (ji,2) = zmoek (ji,2) + zflagkb * ztemp  
               ! Save Monin Obukhov depth
               zmoa  (ji)   = zmob
           
               zvisc(ji,1) = zvisc(ji,1) + zflagkb * avmu(ji,jj,jkm1)
               zvisc(ji,2) = zvisc(ji,2) + zflagkb * avmu(ji,jj,jk  )
               zvisc(ji,3) = zvisc(ji,3) + zflagkb * avmu(ji,jj,jk+1)
               zvisc(ji,4) = zvisc(ji,4) + zflagkb * avmu(ji,jj,jkp2)
               
               zdift(ji,1) = zdift(ji,1) + zflagkb * avt (ji,jj,jkm1)
               zdift(ji,2) = zdift(ji,2) + zflagkb * avt (ji,jj,jk  )
               zdift(ji,3) = zdift(ji,3) + zflagkb * avt (ji,jj,jk+1)
               zdift(ji,4) = zdift(ji,4) + zflagkb * avt (ji,jj,jkp2)

#if defined key_zdfddm
               zdifs(ji,1) = zdifs(ji,1) + zflagkb * avs (ji,jj,jkm1)
               zdifs(ji,2) = zdifs(ji,2) + zflagkb * avs (ji,jj,jk  )
               zdifs(ji,3) = zdifs(ji,3) + zflagkb * avs (ji,jj,jk+1)
               zdifs(ji,4) = zdifs(ji,4) + zflagkb * avs (ji,jj,jkp2)
#endif               
               ! Save the Richardson number 
               zria  (ji)   = zrib  
#if defined key_c1d
               ! store buoyancy length scale
               buof(ji,jj,jk) = zbuofdep * tmask(ji,jj,jk) 
               ! store Monin Obukhov
               zmob           = zstabl * zmob + ( 1.0 - zstabl) * fsdept(ji,jj,1)
               mols(ji,jj,jk) = MIN( zmob , zhmax(ji) ) * tmask(ji,jj,jk) 
               ! Bulk Richardson number
               rib(ji,jj,jk)  = zrib * tmask(ji,jj,jk)             
#endif               
            END DO
         END DO
         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
         ! III PROCESS THE PIPE
         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
         
         DO ji = fs_2, fs_jpim1 
            
            ! Find the boundary layer depth zhbl
            ! ----------------------------------------
            
            ! Interpolate monin Obukhov and critical Ri mumber depths   
            ztemp = zdept(ji,2) - zdept(ji,1)
            zflag = ( Ricr - zriblk(ji,1) ) / ( zriblk(ji,2) - zriblk(ji,1)  + epsln )
            zhrib = zdept(ji,1) + zflag * ztemp      

            IF( zriblk(ji,2) < Ricr ) zhrib = zhmax(ji) 
         
            IF( zmoek(ji,2) < zdept(ji,2) ) THEN
               IF ( zmoek(ji,1) < 0. ) THEN
                  zmob = zdept(ji,2) - epsln
               ELSE
                  zmob = ztemp + zmoek(ji,1) - zmoek(ji,2)
                  zmob = ( zmoek(ji,1) * zdept(ji,2) - zmoek(ji,2) * zdept(ji,1) ) / zmob
                  zmob = MAX( zmob , zdept(ji,1) + epsln )                
               ENDIF
            ELSE           
               zmob = zhmax(ji) 
            ENDIF
            ztemp   = MIN( zmob , zmoek(ji,0) )
                         
            ! Finally, the boundary layer depth, zhbl 
            zhbl(ji) = MAX( fsdept(ji,jj,1) + epsln, MIN( zhrib , ztemp ) )
            
            ! Save hkpp for further diagnostics (optional)
            hkpp(ji,jj) = zhbl(ji) * tmask(ji,jj,1) 
          
            ! Correct mask if zhbl < fsdepw(ji,jj,2) for no viscosity/diffusivity enhancement at fsdepw(ji,jj,2)
            !     zflag = 1 if zhbl(ji) > fsdepw(ji,jj,2)
            IF( zhbl(ji) < fsdepw(ji,jj,2) ) zmask(ji,2) = 0.
          
            
            !  Velocity scales at depth zhbl
            ! -----------------------------------
            
            !  Compute bouyancy forcing down to zhbl
            ztemp    = -hbf * zhbl(ji)
            zatt1    = 1.0 - ( rabs * EXP( ztemp / xsi1 ) + ( 1.0 - rabs ) * EXP( ztemp / xsi2 ) )
            zbuofdep = zBo(ji,jj) + zBosol(ji,jj) * zatt1
            zstabl   = 0.5 + SIGN( 0.5 , zbuofdep ) 

            zbuofdep = zbuofdep + zstabl * epsln

            zscale = zstabl + ( 1.0 - zstabl ) * epsilon          
            zehat  = vonk * zscale * zhbl(ji) * zbuofdep
            zucube = zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj)              
            zeta   = zehat / ( zucube + epsln )
            
            IF( zehat > 0. ) THEN
               ! Stable case
               zws  = vonk * zustar(ji,jj) / ( 1.0 + rconc1 * zeta )
               zwm  = zws
            ELSE
               ! Unstable case
#if defined key_kpplktb
               ! use lookup table
               zd     = zehat - dehatmin
               il     = INT( zd / dezehat )
               il     = MIN( il, nilktbm1 )
               il     = MAX( il, 1 )
               
               ud     = zustar(ji,jj) - ustmin
               jl     = INT( ud / deustar )
               jl     = MIN( jl, njlktbm1 )
               jl     = MAX( jl, 1 )
               
               zfrac  = zd / dezehat - FLOAT( il )  
               ufrac  = ud / deustar - FLOAT( jl )
               zwas   = ( 1. - zfrac ) * wslktb(il,jl+1) + zfrac * wslktb(il+1,jl+1)
               zwbs   = ( 1. - zfrac ) * wslktb(il,jl  ) + zfrac * wslktb(il+1,jl  )
               zwam   = ( 1. - zfrac ) * wmlktb(il,jl+1) + zfrac * wmlktb(il+1,jl+1)
               zwbm   = ( 1. - zfrac ) * wmlktb(il,jl  ) + zfrac * wmlktb(il+1,jl  )
               !
               zws    = ( 1. - ufrac ) * zwbs + ufrac * zwas
               zwm    = ( 1. - ufrac ) * zwbm + ufrac * zwam
#else
               ! use analytical functions
               zconm  = 0.5 + SIGN( 0.5, ( rzetam - zeta) )
               zcons  = 0.5 + SIGN( 0.5, ( rzetas - zeta) )
               
               ! Momentum : zeta < rzetam (zconm = 1)
               ! Scalars  : zeta < rzetas (zcons = 1) 
               zwconm = zustar(ji,jj) * vonk * ( ( ABS( rconam - rconcm * zeta) )**pthird )
               zwcons = zustar(ji,jj) * vonk * ( ( ABS( rconas - rconcs * zeta) )**pthird )
               
               ! Momentum : rzetam <= zeta < 0 (zconm = 0)
               ! Scalars  : rzetas <= zeta < 0 (zcons = 0)	
               zwmun  = SQRT( ABS( 1.0 - rconc2 * zeta ) )
               zwsun  = vonk * zustar(ji,jj) * zwmun
               zwmun  = vonk * zustar(ji,jj) * SQRT(zwmun)
               !
               zwm    = zconm * zwconm + ( 1.0 - zconm ) * zwmun
               zws    = zcons * zwcons + ( 1.0 - zcons ) * zwsun
               
#endif
            ENDIF
            
            
            ! Viscosity, diffusivity values and derivatives at h
            ! --------------------------------------------------------
            
            ! check between at which interfaces is located zhbl(ji)
            ! ztemp = 1, zdepw(ji,2) < zhbl <  zdepw(ji,3)
            ! ztemp = 0, zdepw(ji,1) < zhbl <  zdepw(ji,2)
            ztemp  =  0.5 + SIGN( 0.5, ( zhbl(ji) - zdepw(ji,2) ) )  
            zdep21 =   zdepw(ji,2) - zdepw(ji,1) + epsln
            zdep32 =   zdepw(ji,3) - zdepw(ji,2) + epsln
            zdep43 =   zdepw(ji,4) - zdepw(ji,3) + epsln  
            
            ! Compute R as in LMD94, eq D5b
            zdelta =  ( zhbl(ji) - zdepw(ji,2) ) *         ztemp   / zdep32   &
               &    + ( zhbl(ji) - zdepw(ji,1) ) * ( 1.0 - ztemp ) / zdep21 
            
            ! Compute the vertical derivative of viscosities (zdzh) at z=zhbl(ji)
            zdzup  =  ( zvisc(ji,2) - zvisc(ji,3) ) *         ztemp   / zdep32 &
               &    + ( zvisc(ji,1) - zvisc(ji,2) ) * ( 1.0 - ztemp ) / zdep21
            
            zdzdn  =  ( zvisc(ji,3) - zvisc(ji,4) ) *         ztemp   / zdep43 &
               &    + ( zvisc(ji,2) - zvisc(ji,3) ) * ( 1.0 - ztemp ) / zdep32
            
            ! LMD94, eq D5b :          
            zdzh   = ( 1.0 - zdelta ) * zdzup + zdelta * zdzdn
            zdzh   = MAX( zdzh , 0. )
            
            ! Compute viscosities (zvath) at z=zhbl(ji), LMD94 eq D5a
            zvath  =          ztemp   * ( zvisc(ji,3) + zdzh * ( zdepw(ji,3) - zhbl(ji) ) ) &
               &    + ( 1.0 - ztemp ) * ( zvisc(ji,2) + zdzh * ( zdepw(ji,2) - zhbl(ji) ) )
            
            ! Compute G (zgat1) and its derivative (zdat1) at z=hbl(ji), LMD94 eq 18
            
            ! Vertical derivative of velocity scale divided by velocity scale squared at z=hbl(ji) 
            ! (non zero only in stable conditions)
            zflag  =  -zstabl * rconc1 * zbuofdep / ( zucube * zustar(ji,jj) + epsln )
            
            ! G at its derivative at z=hbl:
            zgat1  = zvath  / ( zhbl(ji) * ( zwm + epsln )  )
            zdat1  = -zdzh  / ( zwm + epsln ) -  zflag * zvath / zhbl(ji)
            
            ! G coefficients, LMD94 eq 17
            za2m(ji) = -2.0 + 3.0 * zgat1 - zdat1
            za3m(ji) =  1.0 - 2.0 * zgat1 + zdat1

            
            ! Compute the vertical derivative of temperature diffusivities (zdzh) at z=zhbl(ji)
            zdzup  =  ( zdift(ji,2) - zdift(ji,3) ) *         ztemp   / zdep32 &
               &    + ( zdift(ji,1) - zdift(ji,2) ) * ( 1.0 - ztemp ) / zdep21
            
            zdzdn  =  ( zdift(ji,3) - zdift(ji,4) ) *         ztemp   / zdep43 &
               &    + ( zdift(ji,2) - zdift(ji,3) ) * ( 1.0 - ztemp ) / zdep32
            
            ! LMD94, eq D5b :          
            zdzh   = ( 1.0 - zdelta ) * zdzup + zdelta * zdzdn
            zdzh   = MAX( zdzh , 0. )
            
            
            ! Compute diffusivities (zvath) at z=zhbl(ji), LMD94 eq D5a
            zvath  =          ztemp   * ( zdift(ji,3) + zdzh * ( zdepw(ji,3) - zhbl(ji) ) ) &
               &    + ( 1.0 - ztemp ) * ( zdift(ji,2) + zdzh * ( zdepw(ji,2) - zhbl(ji) ) )
                        
            ! G at its derivative at z=hbl:
            zgat1  = zvath  / ( zhbl(ji) * ( zws + epsln )  )
            zdat1  = -zdzh  / ( zws + epsln ) -  zflag * zvath / zhbl(ji)
            
            ! G coefficients, LMD94 eq 17
            za2t(ji) = -2.0 + 3.0 * zgat1 - zdat1
            za3t(ji) =  1.0 - 2.0 * zgat1 + zdat1

#if defined key_zdfddm
            ! Compute the vertical derivative of salinities diffusivities (zdzh) at z=zhbl(ji)
            zdzup  =  ( zdifs(ji,2) - zdifs(ji,3) ) *         ztemp   / zdep32 &
               &    + ( zdifs(ji,1) - zdifs(ji,2) ) * ( 1.0 - ztemp ) / zdep21
            
            zdzdn  =  ( zdifs(ji,3) - zdifs(ji,4) ) *         ztemp   / zdep43 &
               &    + ( zdifs(ji,2) - zdifs(ji,3) ) * ( 1.0 - ztemp ) / zdep32
            
            ! LMD94, eq D5b :          
            zdzh   = ( 1.0 - zdelta ) * zdzup + zdelta * zdzdn
            zdzh   = MAX( zdzh , 0. )           
            
            ! Compute diffusivities (zvath) at z=zhbl(ji), LMD94 eq D5a
            zvath  =          ztemp   * ( zdifs(ji,3) + zdzh * ( zdepw(ji,3) - zhbl(ji) ) ) &
               &    + ( 1.0 - ztemp ) * ( zdifs(ji,2) + zdzh * ( zdepw(ji,2) - zhbl(ji) ) )
                        
            ! G at its derivative at z=hbl:
            zgat1  = zvath  / ( zhbl(ji) * ( zws + epsln )  )
            zdat1  = -zdzh  / ( zws + epsln ) -  zflag * zvath / zhbl(ji)
            
            ! G coefficients, LMD94 eq 17
            za2s(ji) = -2.0 + 3.0 * zgat1 - zdat1
            za3s(ji) =  1.0 - 2.0 * zgat1 + zdat1
#endif

            !-------------------turn off interior matching here------
            !          za2(ji,1) = -2.0
            !          za3(ji,1) =  1.0
            !          za2(ji,2) = -2.0
            !          za3(ji,2) =  1.0
            !--------------------------------------------------------
            
            !  Compute Enhanced Mixing Coefficients (LMD94,eq D6)
            ! ---------------------------------------------------------------
            
            ! Delta 
            zdelta  = ( zhbl(ji)  - zdept(ji,1) ) / ( zdept(ji,2) - zdept(ji,1) + epsln )
            zdelta2 = zdelta * zdelta
            
            !  Mixing coefficients at first level above h (zdept(ji,1)) 
            ! and at first interface in the pipe (zdepw(ji,2))
            
            ! At first T level above h (zdept(ji,1)) (always in the boundary layer)
            zsig    = zdept(ji,1) / zhbl(ji)
            ztemp   = zstabl * zsig  + ( 1.0 - zstabl ) * MIN( zsig , epsilon )
            zehat   = vonk * ztemp * zhbl(ji) * zbuofdep
            zeta    = zehat / ( zucube + epsln)
            zwst    = vonk * zustar(ji,jj) / ( ABS( 1.0 + rconc1 * zeta ) + epsln)
            zwm     = zstabl * zwst + ( 1.0 - zstabl ) * zwm
            zws     = zstabl * zwst + ( 1.0 - zstabl ) * zws

            zkm1m  = zhbl(ji) * zwm * zsig * ( 1.0 + zsig * ( za2m(ji) + zsig * za3m(ji) ) )
            zkm1t  = zhbl(ji) * zws * zsig * ( 1.0 + zsig * ( za2t(ji) + zsig * za3t(ji) ) )
#if defined key_zdfddm
            zkm1s  = zhbl(ji) * zws * zsig * ( 1.0 + zsig * ( za2s(ji) + zsig * za3s(ji) ) )
#endif                        
            ! At first W level in the pipe (zdepw(ji,2)) (not always in the boundary layer ):
            zsig    = MIN( zdepw(ji,2) / zhbl(ji) , 1.0 )
            ztemp   = zstabl * zsig + ( 1.0 - zstabl ) * MIN( zsig , epsilon )
            zehat   = vonk * ztemp * zhbl(ji) * zbuofdep
            zeta    = zehat / ( zucube + epsln )
            zwst    = vonk * zustar(ji,jj) / ( ABS( 1.0 + rconc1 * zeta ) + epsln)
            zws     = zstabl * zws + ( 1.0 - zstabl ) * zws
            zwm     = zstabl * zws + ( 1.0 - zstabl ) * zwm

            zkmpm(ji) = zhbl(ji) * zwm * zsig * ( 1.0 + zsig * ( za2m(ji) + zsig * za3m(ji) ) )
            zkmpt(ji) = zhbl(ji) * zws * zsig * ( 1.0 + zsig * ( za2t(ji) + zsig * za3t(ji) ) )
#if defined key_zdfddm
            zkmps(ji) = zhbl(ji) * zws * zsig * ( 1.0 + zsig * ( za2s(ji) + zsig * za3s(ji) ) )
#endif  
      
            ! check if this point is in the boundary layer,else take interior viscosity/diffusivity:
            zflag       = 0.5 + SIGN( 0.5, ( zhbl(ji) - zdepw(ji,2) ) )
            zkmpm(ji) = zkmpm(ji) * zflag + ( 1.0 - zflag ) * zvisc(ji,2)
            zkmpt(ji) = zkmpt(ji) * zflag + ( 1.0 - zflag ) * zdift(ji,2)
#if defined key_zdfddm
            zkmps(ji) = zkmps(ji) * zflag + ( 1.0 - zflag ) * zdifs(ji,2)
#endif

            ! Enhanced viscosity/diffusivity at zdepw(ji,2)
            ztemp     = ( 1.0 - 2.0 * zdelta + zdelta2 ) * zkm1m + zdelta2 * zkmpm(ji)
            zkmpm(ji) = ( 1.0 - zdelta ) * zvisc(ji,2) + zdelta * ztemp
            ztemp     = ( 1.0 - 2.0 * zdelta + zdelta2 ) * zkm1t + zdelta2 * zkmpt(ji)
            zkmpt(ji) = ( 1.0 - zdelta ) * zdift(ji,2) + zdelta * ztemp
#if defined key_zdfddm
            ztemp     = ( 1.0 - 2.0 * zdelta + zdelta2 ) * zkm1s + zdelta2 * zkmps(ji)
            zkmps(ji) = ( 1.0 - zdelta ) * zdifs(ji,2) + zdelta * ztemp
#endif            

         END DO
         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
         ! IV. Compute vertical eddy viscosity and diffusivity coefficients
         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
         
         DO jk  = 2, jkmax
            
            ! Compute turbulent velocity scales on the interfaces
            ! --------------------------------------------------------
            DO  ji = fs_2, fs_jpim1
               zbuofdep = zBo(ji,jj) + zBosol(ji,jj) * zatt1
               zstabl   = 0.5 + SIGN( 0.5 , zbuofdep ) 
               zbuofdep = zbuofdep + zstabl * epsln          
               zsig    = fsdepw(ji,jj,jk) / zhbl(ji)
               ztemp   = zstabl * zsig + ( 1. - zstabl ) * MIN( zsig , epsilon )
               zehat   = vonk * ztemp * zhbl(ji) * zbuofdep
               zucube  = zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj)
               zeta    = zehat / ( zucube + epsln )

               IF( zehat > 0. ) THEN
                  ! Stable case
                  zws  = vonk * zustar(ji,jj) / ( 1.0 + rconc1 * zeta )
                  zwm  = zws
               ELSE
                  ! Unstable case
#if defined key_kpplktb
                  ! use lookup table
                  zd     = zehat - dehatmin
                  il     = INT( zd / dezehat )
                  il     = MIN( il, nilktbm1 )
                  il     = MAX( il, 1 )
                  
                  ud     = zustar(ji,jj) - ustmin
                  jl     = INT( ud / deustar )
                  jl     = MIN( jl, njlktbm1 )
                  jl     = MAX( jl, 1 )
                  
                  zfrac  = zd / dezehat - FLOAT( il )  
                  ufrac  = ud / deustar - FLOAT( jl )
                  zwas   = ( 1. - zfrac ) * wslktb(il,jl+1) + zfrac * wslktb(il+1,jl+1)
                  zwbs   = ( 1. - zfrac ) * wslktb(il,jl  ) + zfrac * wslktb(il+1,jl  )
                  zwam   = ( 1. - zfrac ) * wmlktb(il,jl+1) + zfrac * wmlktb(il+1,jl+1)
                  zwbm   = ( 1. - zfrac ) * wmlktb(il,jl  ) + zfrac * wmlktb(il+1,jl  )
                  !
                  zws    = ( 1. - ufrac ) * zwbs + ufrac * zwas
                  zwm    = ( 1. - ufrac ) * zwbm + ufrac * zwam
#else
                  ! use analytical functions
                  zconm  = 0.5 + SIGN( 0.5, ( rzetam - zeta) )
                  zcons  = 0.5 + SIGN( 0.5, ( rzetas - zeta) )
                  
                  ! Momentum : zeta < rzetam (zconm = 1)
                  ! Scalars  : zeta < rzetas (zcons = 1) 
                  zwconm = zustar(ji,jj) * vonk * ( ( ABS( rconam - rconcm * zeta) )**pthird )
                  zwcons = zustar(ji,jj) * vonk * ( ( ABS( rconas - rconcs * zeta) )**pthird )
                  
                  ! Momentum : rzetam <= zeta < 0 (zconm = 0)
                  ! Scalars  : rzetas <= zeta < 0 (zcons = 0)	
                  zwmun  = SQRT( ABS( 1.0 - rconc2 * zeta ) )
                  zwsun  = vonk * zustar(ji,jj) * zwmun
                  zwmun  = vonk * zustar(ji,jj) * SQRT(zwmun)
                  !
                  zwm    = zconm * zwconm + ( 1.0 - zconm ) * zwmun
                  zws    = zcons * zwcons + ( 1.0 - zcons ) * zwsun
                  
#endif
               ENDIF
               
               zblcm(ji,jk) = zhbl(ji) * zwm * zsig  * ( 1.0 + zsig * ( za2m(ji) + zsig * za3m(ji) ) )
               zblct(ji,jk) = zhbl(ji) * zws * zsig  * ( 1.0 + zsig * ( za2t(ji) + zsig * za3t(ji) ) )
#if defined key_zdfddm
               zblcs(ji,jk) = zhbl(ji) * zws * zsig  * ( 1.0 + zsig * ( za2s(ji) + zsig * za3s(ji) ) )
#endif              
               !  Compute Nonlocal transport term = ghats * <ws>o
               ! ----------------------------------------------------
               ghats(ji,jj,jk-1) = ( 1. - zstabl ) * rcg / ( zws * zhbl(ji) + epsln ) * tmask(ji,jj,jk)

            END DO
         END DO     
         !  Combine interior and boundary layer coefficients and nonlocal term
         ! -----------------------------------------------------------------------
         DO jk = 2, jpkm1   
            DO ji = fs_2, fs_jpim1
               zflag = zmask(ji,jk) * zmask(ji,jk+1)
               zviscos(ji,jj,jk) = ( 1.0 - zmask(ji,jk) )         * avmu (ji,jj,jk) & ! interior viscosities
                  &              +                        zflag   * zblcm(ji,jk    ) & ! boundary layer viscosities
                  &              + zmask(ji,jk) * ( 1.0 - zflag ) * zkmpm(ji       )   ! viscosity enhancement at W_level near zhbl
               
               zviscos(ji,jj,jk) = zviscos(ji,jj,jk) * tmask(ji,jj,jk)   

            
               zdiffut(ji,jj,jk) = ( 1.0 - zmask(ji,jk) )          * avt (ji,jj,jk) & ! interior diffusivities 
                  &              +                        zflag   * zblct(ji,jk   ) & ! boundary layer diffusivities
                  &              + zmask(ji,jk) * ( 1.0 - zflag ) * zkmpt(ji      )   ! diffusivity enhancement at W_level near zhbl
                       
               zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) * tmask(ji,jj,jk) 
#if defined key_zdfddm
               zdiffus(ji,jj,jk) = ( 1.0 - zmask(ji,jk) )         * avs (ji,jj,jk) & ! interior diffusivities 
                  &              +                        zflag   * zblcs(ji,jk   ) & ! boundary layer diffusivities
                  &              + zmask(ji,jk) * ( 1.0 - zflag ) * zkmps(ji      )   ! diffusivity enhancement at W_level near zhbl
               zdiffus(ji,jj,jk) = zdiffus(ji,jj,jk) * tmask(ji,jj,jk) 
#endif               
               ! Non local flux in the boundary layer only
               ghats(ji,jj,jk-1) = zmask(ji,jk) * ghats(ji,jj,jk-1)

            ENDDO
         END DO
         !                                             ! ===============
      END DO                                           !   End of slab
      !                                                ! ===============

      ! Lateral boundary conditions on zvicos and zdiffus  (sign unchanged)
      CALL lbc_lnk( zviscos(:,:,:), 'U', 1. )  ; CALL lbc_lnk( zdiffut(:,:,:), 'W', 1. )  
#if defined key_zdfddm  
      CALL lbc_lnk( zdiffus(:,:,:), 'W', 1. ) 
#endif

      SELECT CASE ( nn_ave )
         !
         CASE ( 0 )             ! no viscosity and diffusivity smoothing

            DO jk = 2, jpkm1
               DO jj = 2, jpjm1
                  DO ji = fs_2, fs_jpim1
                     avmu(ji,jj,jk) = ( zviscos(ji,jj,jk) + zviscos(ji+1,jj,jk) ) &
                        &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk)
                     
                     avmv(ji,jj,jk) = ( zviscos(ji,jj,jk) + zviscos(ji,jj+1,jk) ) &
                        &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk)
                     
                     avt (ji,jj,jk) =  zdiffut(ji,jj,jk) * tmask(ji,jj,jk)  
#if defined key_zdfddm     
                     avs (ji,jj,jk) =  zdiffus(ji,jj,jk) * tmask(ji,jj,jk)  
#endif
                  END DO
               END DO
            END DO
            
         CASE ( 1 )                ! viscosity and diffusivity smoothing
            !                      
            !           ( 1/2  1  1/2 )              ( 1/2  1/2 )             ( 1/2  1  1/2 )
            ! avt = 1/8 ( 1    2  1   )   avmu = 1/4 ( 1    1   )   avmv= 1/4 ( 1/2  1  1/2 )
            !           ( 1/2  1  1/2 )              ( 1/2  1/2 )
  
            DO jk = 2, jpkm1
               DO jj = 2, jpjm1
                  DO ji = fs_2, fs_jpim1

                     avmu(ji,jj,jk) = (      zviscos(ji  ,jj  ,jk) + zviscos(ji+1,jj  ,jk)   &
                        &              +.5*( zviscos(ji  ,jj-1,jk) + zviscos(ji+1,jj-1,jk)   &
                        &                   +zviscos(ji  ,jj+1,jk) + zviscos(ji+1,jj+1,jk) ) ) * eumean(ji,jj,jk)
                     
                     avmv(ji,jj,jk) = (      zviscos(ji  ,jj  ,jk) + zviscos(ji  ,jj+1,jk)   &
                        &              +.5*( zviscos(ji-1,jj  ,jk) + zviscos(ji-1,jj+1,jk)   &
                        &                   +zviscos(ji+1,jj  ,jk) + zviscos(ji+1,jj+1,jk) ) ) * evmean(ji,jj,jk)
 
                     avt (ji,jj,jk) = ( .5*( zdiffut(ji-1,jj+1,jk) + zdiffut(ji-1,jj-1,jk)    &
                        &                   +zdiffut(ji+1,jj+1,jk) + zdiffut(ji+1,jj-1,jk) )  &
                        &              +1.*( zdiffut(ji-1,jj  ,jk) + zdiffut(ji  ,jj+1,jk)    &
                        &                   +zdiffut(ji  ,jj-1,jk) + zdiffut(ji+1,jj  ,jk) )  &
                        &              +2.*  zdiffut(ji  ,jj  ,jk)                          ) * etmean(ji,jj,jk)
#if defined key_zdfddm   
                     avs (ji,jj,jk) = ( .5*( zdiffus(ji-1,jj+1,jk) + zdiffus(ji-1,jj-1,jk)    &
                        &                   +zdiffus(ji+1,jj+1,jk) + zdiffus(ji+1,jj-1,jk) )  &
                        &              +1.*( zdiffus(ji-1,jj  ,jk) + zdiffus(ji  ,jj+1,jk)    &
                        &                   +zdiffus(ji  ,jj-1,jk) + zdiffus(ji+1,jj  ,jk) )  &
                        &              +2.*  zdiffus(ji  ,jj  ,jk)                          ) * etmean(ji,jj,jk)  
#endif                
                  END DO
               END DO
            END DO
         
         END SELECT

         DO jk = 2, jpkm1                       ! vertical slab
            !
            !  Minimum value on the eddy diffusivity
            ! ----------------------------------------
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  avt(ji,jj,jk) = MAX( avt(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
#if defined key_zdfddm  
                  avs(ji,jj,jk) = MAX( avs(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
#endif
               END DO
            END DO

            !
            ! Minimum value on the eddy viscosity
            ! ----------------------------------------
            DO jj = 1, jpj
               DO ji = 1, jpi
                  avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), avmb(jk) ) * umask(ji,jj,jk)
                  avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), avmb(jk) ) * vmask(ji,jj,jk)
               END DO
            END DO
            !
         END DO

         ! Lateral boundary conditions on avt  (sign unchanged)
         CALL lbc_lnk( hkpp(:,:), 'T', 1. )

         ! Lateral boundary conditions on avt  (sign unchanged)
         CALL lbc_lnk( avt(:,:,:), 'W', 1. )
#if defined key_zdfddm  
         CALL lbc_lnk( avs(:,:,:), 'W', 1. )  
#endif
         ! Lateral boundary conditions (avmu,avmv) (U- and V- points, sign unchanged)
         CALL lbc_lnk( avmu(:,:,:), 'U', 1. )   ;    CALL lbc_lnk( avmv(:,:,:), 'V', 1. )  
 
         IF(ln_ctl) THEN
#if defined key_zdfddm
            CALL prt_ctl(tab3d_1=avt , clinfo1=' kpp - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
#else
            CALL prt_ctl(tab3d_1=avt , clinfo1=' kpp - t: ', ovlap=1, kdim=jpk)
#endif
            CALL prt_ctl(tab3d_1=avmu, clinfo1=' kpp - u: ', mask1=umask,  &
               &         tab3d_2=avmv, clinfo2=      ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
         ENDIF

      CALL wrk_dealloc( jpi, zmoa, zekman, zhmax, zria, zhbl )
      CALL wrk_dealloc( jpi, za2m, za3m, zkmpm, za2t, za3t, zkmpt )
      CALL wrk_dealloc( jpi,2, zriblk )
      CALL wrk_dealloc( jpi,3, zmoek, kjstart = 0 )
      CALL wrk_dealloc( jpi,3, zdept )
      CALL wrk_dealloc( jpi,4, zdepw, zdift, zvisc )
      CALL wrk_dealloc( jpi,jpj, zBo, zBosol, zustar )
      CALL wrk_dealloc( jpi,jpk, zmask, zblcm, zblct )
#if defined key_zdfddm
      CALL wrk_dealloc( jpi,4, zdifs )
      CALL wrk_dealloc( jpi, zmoa, za2s, za3s, zkmps )
      CALL wrk_dealloc( jpi,jpk, zblcs )
      CALL wrk_dealloc( jpi,jpi,jpk, zdiffus )
#endif
      !
      IF( nn_timing == 1 )  CALL timing_stop('zdf_kpp')
      !
   END SUBROUTINE zdf_kpp


   SUBROUTINE tra_kpp( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE tra_kpp  ***
      !!
      !! ** Purpose :   compute and add to the tracer trend the non-local tracer flux
      !!
      !! ** Method  :   ???
      !!----------------------------------------------------------------------
      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
      !!----------------------------------------------------------------------
      INTEGER, INTENT(in) :: kt
      INTEGER :: ji, jj, jk
      !
      IF( nn_timing == 1 )  CALL timing_start('tra_kpp')
      !
      IF( kt == nit000 ) THEN
         IF(lwp) WRITE(numout,*) 
         IF(lwp) WRITE(numout,*) 'tra_kpp : KPP non-local tracer fluxes'
         IF(lwp) WRITE(numout,*) '~~~~~~~   '
      ENDIF

      IF( l_trdtra )   THEN                    !* Save ta and sa trends
         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsa(:,:,:,jp_sal)
      ENDIF

      ! add non-local temperature and salinity flux ( in convective case only)
      DO jk = 1, jpkm1
         DO jj = 2, jpjm1 
            DO ji = fs_2, fs_jpim1
               tsa(ji,jj,jk,jp_tem) =  tsa(ji,jj,jk,jp_tem)                      &
                  &                 - (  ghats(ji,jj,jk  ) * avt  (ji,jj,jk  )   & 
                  &                    - ghats(ji,jj,jk+1) * avt  (ji,jj,jk+1) ) * wt0(ji,jj) / fse3t(ji,jj,jk)
               tsa(ji,jj,jk,jp_sal) =  tsa(ji,jj,jk,jp_sal)                      &
                  &                 - (  ghats(ji,jj,jk  ) * fsavs(ji,jj,jk  )   & 
                  &                    - ghats(ji,jj,jk+1) * fsavs(ji,jj,jk+1) ) * ws0(ji,jj) / fse3t(ji,jj,jk)
            END DO
         END DO
      END DO

      ! save the non-local tracer flux trends for diagnostic
      IF( l_trdtra )   THEN
         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
!!bug gm jpttdzdf ==> jpttkpp
         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )
         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )
         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds )
      ENDIF

      IF(ln_ctl) THEN
         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' kpp  - Ta: ', mask1=tmask,   &
         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
      ENDIF
      !
      IF( nn_timing == 1 )  CALL timing_stop('tra_kpp')
      !
   END SUBROUTINE tra_kpp

#if defined key_top
   !!----------------------------------------------------------------------
   !!   'key_top'                                                TOP models
   !!----------------------------------------------------------------------
   SUBROUTINE trc_kpp( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE trc_kpp  ***
      !!
      !! ** Purpose :   compute and add to the tracer trend the non-local
      !!                tracer flux
      !!
      !! ** Method  :   ???
      !!
      !! history :
      !!            9.0  ! 2005-11 (G. Madec)  Original code
      !!       NEMO 3.3  ! 2010-06 (C. Ethe )  Adapted to passive tracers
      !!----------------------------------------------------------------------
      USE trc
      USE prtctl_trc          ! Print control
      !
      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
      !
      INTEGER  ::   ji, jj, jk, jn      ! Dummy loop indices
      CHARACTER (len=35) :: charout
      REAL(wp) ::   ztra, zflx
      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrtrd
      !!----------------------------------------------------------------------

      IF( kt == nit000 ) THEN
         IF(lwp) WRITE(numout,*) 
         IF(lwp) WRITE(numout,*) 'trc_kpp : KPP non-local tracer fluxes'
         IF(lwp) WRITE(numout,*) '~~~~~~~   '
      ENDIF

      IF( l_trdtrc )  ALLOCATE( ztrtrd(jpi,jpj,jpk) )
      !
      DO jn = 1, jptra
         !
         IF( l_trdtrc )  ztrtrd(:,:,:)  = tra(:,:,:,jn)
         ! add non-local on passive tracer flux ( in convective case only)
         DO jk = 1, jpkm1
            DO jj = 2, jpjm1 
               DO ji = fs_2, fs_jpim1
                  ! Surface tracer flux for non-local term 
                  zflx = - ( sfx (ji,jj) * tra(ji,jj,1,jn) * rcs ) * tmask(ji,jj,1)
                  ! compute the trend
                  ztra = - ( ghats(ji,jj,jk  ) * fsavs(ji,jj,jk  )   &
                  &        - ghats(ji,jj,jk+1) * fsavs(ji,jj,jk+1) ) * zflx / fse3t(ji,jj,jk)
                  ! add the trend to the general trend
                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn)  + ztra
               END DO
            END DO
         END DO
         !
         IF( l_trdtrc ) THEN         ! save the non-local tracer flux trends for diagnostic
            ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:)
            CALL trd_tra( kt, 'TRC', jn, jptra_zdf, ztrtrd(:,:,:) )
         ENDIF
         !
      END DO
      IF( l_trdtrc )  DEALLOCATE( ztrtrd )
      IF( ln_ctl )   THEN
         WRITE(charout, FMT="(' kpp')")  ;  CALL prt_ctl_trc_info(charout)
         CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
      ENDIF
      !
   END SUBROUTINE trc_kpp

#else
   !!----------------------------------------------------------------------
   !!   NO 'key_top'           DUMMY routine                  No TOP models
   !!----------------------------------------------------------------------
   SUBROUTINE trc_kpp( kt )         ! Dummy routine
      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
      WRITE(*,*) 'tra_kpp: You should not have seen this print! error?', kt
   END SUBROUTINE trc_kpp
#endif

   SUBROUTINE zdf_kpp_init
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE zdf_kpp_init  ***
      !!                     
      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
      !!      viscosity when using a kpp turbulent closure scheme
      !!
      !! ** Method  :   Read the namkpp namelist and check the parameters
      !!      called at the first timestep (nit000)
      !!
      !! ** input   :   Namlist namkpp
      !!----------------------------------------------------------------------
      INTEGER  ::   ji, jj, jk     ! dummy loop indices
      INTEGER  ::   ios            ! local integer
#if ! defined key_kppcustom
      INTEGER  ::   jm             ! dummy loop indices     
      REAL(wp) ::   zref, zdist    ! tempory scalars
#endif
#if defined key_kpplktb
      REAL(wp) ::   zustar, zucube, zustvk, zeta, zehat   ! tempory scalars
#endif
      REAL(wp) ::   zhbf           ! tempory scalars
      LOGICAL  ::   ll_kppcustom   ! 1st ocean level taken as surface layer
      LOGICAL  ::   ll_kpplktb     ! Lookup table for turbul. velocity scales 
      !!
      NAMELIST/namzdf_kpp/ ln_kpprimix, rn_difmiw, rn_difsiw, rn_riinfty, rn_difri, rn_bvsqcon, rn_difcon, nn_ave
      !!----------------------------------------------------------------------
      !
      IF( nn_timing == 1 )  CALL timing_start('zdf_kpp_init')
      !
      REWIND( numnam_ref )              ! Namelist namzdf_kpp in reference namelist : Vertical eddy diffivity and viscosity using kpp turbulent closure scheme
      READ  ( numnam_ref, namzdf_kpp, IOSTAT = ios, ERR = 901)
901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_kpp in reference namelist', lwp )

      REWIND( numnam_cfg )              ! Namelist namzdf_kpp in configuration namelist : Vertical eddy diffivity and viscosity using kpp turbulent closure scheme
      READ  ( numnam_cfg, namzdf_kpp, IOSTAT = ios, ERR = 902 )
902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_kpp in configuration namelist', lwp )
      IF(lwm) WRITE ( numond, namzdf_kpp )

      IF(lwp) THEN                    ! Control print
         WRITE(numout,*)
         WRITE(numout,*) 'zdf_kpp_init : K-Profile Parameterisation'
         WRITE(numout,*) '~~~~~~~~~~~~'
         WRITE(numout,*) '   Namelist namzdf_kpp : set tke mixing parameters'
         WRITE(numout,*) '     Shear instability mixing                      ln_kpprimix = ', ln_kpprimix
         WRITE(numout,*) '     max. internal wave viscosity                  rn_difmiw   = ', rn_difmiw
         WRITE(numout,*) '     max. internal wave diffusivity                rn_difsiw   = ', rn_difsiw
         WRITE(numout,*) '     Richardson Number limit for shear instability rn_riinfty  = ', rn_riinfty
         WRITE(numout,*) '     max. shear mixing at Rig = 0                  rn_difri    = ', rn_difri
         WRITE(numout,*) '     Brunt-Vaisala squared for max. convection     rn_bvsqcon  = ', rn_bvsqcon
         WRITE(numout,*) '     max. mix. in interior convec.                 rn_difcon   = ', rn_difcon
         WRITE(numout,*) '     horizontal average flag                       nn_ave      = ', nn_ave
      ENDIF

      !                              ! allocate zdfkpp arrays
      IF( zdf_kpp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_kpp_init : unable to allocate arrays' )

      ll_kppcustom = .FALSE.
      ll_kpplktb   = .FALSE.

#if defined key_kppcustom
      ll_kppcustom = .TRUE.
#endif
#if defined key_kpplktb
      ll_kpplktb   = .TRUE.
#endif
      IF(lwp) THEN
         WRITE(numout,*) '     Lookup table for turbul. velocity scales ll_kpplktb   = ', ll_kpplktb
         WRITE(numout,*) '     1st ocean level taken as surface layer   ll_kppcustom = ', ll_kppcustom
      ENDIF
      
      IF( lk_zdfddm) THEN
         IF(lwp) THEN
            WRITE(numout,*)
            WRITE(numout,*) '    Double diffusion mixing on temperature and salinity '
            WRITE(numout,*) '    CAUTION : done in routine zdfkpp, not in routine zdfddm '
         ENDIF
      ENDIF
      

      !set constants not in namelist
      !-----------------------------
      Vtc  = rconcv * SQRT( 0.2 / ( rconcs * epsilon ) ) / ( vonk * vonk * Ricr )
      rcg  = rcstar * vonk * ( rconcs * vonk * epsilon )**pthird

      IF(lwp) THEN
         WRITE(numout,*)
         WRITE(numout,*) '     Constant value for unreso. turbul. velocity shear Vtc = ', Vtc
         WRITE(numout,*) '     Non-dimensional coef. for nonlocal transport      rcg = ', rcg
       ENDIF

      ! ratt is the attenuation coefficient for solar flux
      ! Should be different is s_coordinate
      DO jk = 1, jpk
         zhbf     = - fsdept(1,1,jk) * hbf
         ratt(jk) = 1.0 - ( rabs * EXP( zhbf / xsi1 ) + ( 1.0 - rabs ) * EXP( zhbf / xsi2 ) )        
      ENDDO

      ! Horizontal average : initialization of weighting arrays 
      ! -------------------
      
      SELECT CASE ( nn_ave )

      CASE ( 0 )                ! no horizontal average
         IF(lwp) WRITE(numout,*) '          no horizontal average on avt, avmu, avmv'
         IF(lwp) WRITE(numout,*) '          only in very high horizontal resolution !'
         ! weighting mean arrays etmean, eumean and evmean
         !           ( 1  1 )                                          ( 1 )
         ! avt = 1/4 ( 1  1 )     avmu = 1/2 ( 1  1 )       avmv=  1/2 ( 1 )
         !                          
         etmean(:,:,:) = 0.e0
         eumean(:,:,:) = 0.e0
         evmean(:,:,:) = 0.e0
         
         DO jk = 1, jpkm1
            DO jj = 2, jpjm1
               DO ji = 2, jpim1   ! vector opt.
                  etmean(ji,jj,jk) = tmask(ji,jj,jk)                     &
                  &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
                  &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
                  
                  eumean(ji,jj,jk) = umask(ji,jj,jk)                     &
                  &  / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk)  )

                  evmean(ji,jj,jk) = vmask(ji,jj,jk)                     &
                  &  / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk)  )
               END DO
            END DO
         END DO

      CASE ( 1 )                ! horizontal average 
         IF(lwp) WRITE(numout,*) '          horizontal average on avt, avmu, avmv'
         ! weighting mean arrays etmean, eumean and evmean
         !           ( 1/2  1  1/2 )              ( 1/2  1/2 )             ( 1/2  1  1/2 )
         ! avt = 1/8 ( 1    2  1   )   avmu = 1/4 ( 1    1   )   avmv= 1/4 ( 1/2  1  1/2 )
         !           ( 1/2  1  1/2 )              ( 1/2  1/2 )
         etmean(:,:,:) = 0.e0
         eumean(:,:,:) = 0.e0
         evmean(:,:,:) = 0.e0
         
         DO jk = 1, jpkm1
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  etmean(ji,jj,jk) = tmask(ji, jj,jk)                           &
                     & / MAX( 1., 2.* tmask(ji,jj,jk)                           &
                     &      +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk)   &
                     &             +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) &
                     &      +1. * ( tmask(ji-1,jj  ,jk) + tmask(ji  ,jj+1,jk)   &
                     &             +tmask(ji  ,jj-1,jk) + tmask(ji+1,jj  ,jk) ) )
                  
                  eumean(ji,jj,jk) = umask(ji,jj,jk)                        &
                     &  / MAX( 1.,   tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk)   &
                     &       +.5 * ( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk)   &
                     &              +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) )  )
                  
                  evmean(ji,jj,jk) = vmask(ji,jj,jk)                        &
                     &  / MAX( 1.,   tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
                     &       +.5 * ( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk)   &
                     &              +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) )  )
               END DO
            END DO
         END DO

      CASE DEFAULT
         WRITE(ctmp1,*) '          bad flag value for nn_ave = ', nn_ave
         CALL ctl_stop( ctmp1 )

      END SELECT
 
      ! Initialization of vertical eddy coef. to the background value
      ! -------------------------------------------------------------
      DO jk = 1, jpk
         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
      END DO

      ! zero the surface flux for non local term and kpp mixed layer depth
      ! ------------------------------------------------------------------
      ghats(:,:,:) = 0.
      wt0  (:,:  ) = 0.
      ws0  (:,:  ) = 0.
      hkpp (:,:  ) = 0. ! just a diagnostic (not essential)

#if ! defined key_kppcustom
      ! compute arrays (del, wz) for reference mean values 
      ! (increase speed for vectorization key_kppcustom not defined)
      del(1:jpk, 1:jpk) = 0.
      DO jk = 1, jpk
         zref = epsilon * fsdept(1,1,jk)    
         DO jm = 1 , jpk
            zdist = zref - fsdepw(1,1,jm)   
            IF( zdist > 0.  ) THEN
               del(jk,jm) = MIN( zdist, fse3t(1,1,jm) ) / zref   
            ELSE
               del(jk,jm) = 0.
            ENDIF
         ENDDO
      ENDDO
#endif

#if defined key_kpplktb
      ! build lookup table for turbulent velocity scales
      dezehat = ( dehatmax - dehatmin ) / nilktbm1
      deustar = ( ustmax   - ustmin   ) / njlktbm1
 
      DO jj = 1, njlktb
         zustar = ( jj - 1) * deustar + ustmin
         zustvk = vonk * zustar 
         zucube = zustar * zustar * zustar 
         DO ji = 1 , nilktb
            zehat = ( ji - 1 ) * dezehat + dehatmin
            zeta   = zehat / ( zucube + epsln )
            IF( zehat >= 0 ) THEN             ! Stable case
               wmlktb(ji,jj) = zustvk / ABS( 1.0 + rconc1 * zeta + epsln )                        
               wslktb(ji,jj) = wmlktb(ji,jj)
            ELSE                                ! Unstable case 
               IF( zeta > rzetam ) THEN
                  wmlktb(ji,jj) = zustvk * ABS( 1.0    - rconc2 * zeta )**pfourth
               ELSE
                  wmlktb(ji,jj) = zustvk * ABS( rconam - rconcm * zeta )**pthird
               ENDIF
               
               IF( zeta > rzetas ) THEN
                  wslktb(ji,jj) = zustvk * SQRT( ABS( 1.0 - rconc2 * zeta ) )
               ELSE
                  wslktb(ji,jj) = zustvk * ABS( rconas - rconcs * zeta )**pthird
               ENDIF
            ENDIF
         END DO
      END DO
#endif
      !
      IF( nn_timing == 1 )  CALL timing_stop('zdf_kpp_init')
      !
   END SUBROUTINE zdf_kpp_init

#else
   !!----------------------------------------------------------------------
   !!   Dummy module :                                        NO KPP scheme
   !!----------------------------------------------------------------------
   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfkpp = .FALSE.   !: KPP flag
CONTAINS
   SUBROUTINE zdf_kpp_init           ! Dummy routine
      WRITE(*,*) 'zdf_kpp_init: You should not have seen this print! error?'
   END SUBROUTINE zdf_kpp_init
   SUBROUTINE zdf_kpp( kt )          ! Dummy routine
      WRITE(*,*) 'zdf_kpp: You should not have seen this print! error?', kt
   END SUBROUTINE zdf_kpp
   SUBROUTINE tra_kpp( kt )          ! Dummy routine
      WRITE(*,*) 'tra_kpp: You should not have seen this print! error?', kt
   END SUBROUTINE tra_kpp
   SUBROUTINE trc_kpp( kt )          ! Dummy routine
      WRITE(*,*) 'trc_kpp: You should not have seen this print! error?', kt
   END SUBROUTINE trc_kpp
#endif

   !!======================================================================
END MODULE zdfkpp